© ESO 2021

1. Introduction

Thanks to the advent of long space-based photometric surveys such as CoRoT (Auvergne et al. 2009), Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2014, 2015), asteroseismology has become the most efficient technique for testing the theory of stellar structure and evolution. In this respect it serves as the golden path to analyse the dynamical processes acting in stellar interiors. One of the key results of the space-based photometry revolution is the detection of ‘mixed modes’ in post-main-sequence stars that behave like acoustic modes in the outer layers and gravity modes in the deep layers. Their dual nature enabled the determination of the internal rotation of subgiant and red giant stars (Deheuvels et al. 2012, 2014, 2015, 2020; Beck et al. 2012; Mosser et al. 2012; Di Mauro et al. 2016, 2018; Gehan et al. 2018) that pointed towards a very efficient angular momentum (AM) transport process absent from theoretical stellar models (Eggenberger et al. 2012, 2017, 2019a; Ceillier et al. 2013; Marques et al. 2013; Cantiello et al. 2014; Spada et al. 2016).

Various transport mechanism candidates have been suggested to reproduce asteroseismic data, from internal gravity waves (Pinçon et al. 2017) and mixed modes (Belkacem et al. 2015a,b) to magnetic instabilities (Spada et al. 2016; Fuller et al. 2019; Eggenberger et al. 2019b; den Hartogh et al. 2020) and internal fossil magnetic fields (Kissin & Thompson 2015; Takahashi & Langer 2021). These processes are still largely investigated by the community with the use of stellar evolution codes including a proper treatment of AM transport during evolution. The question is quite tedious as even non-standard processes invoked to reproduce the solar rotation profile deduced from helioseismic measurements (e.g. Eggenberger et al. 2005, 2019c) are not efficient enough in later stages (Cantiello et al. 2014; den Hartogh et al. 2019). However, some of these studies have only focused on reproducing the core rotation of red giants, while tighter constraints can be drawn from the degree of differential rotation in the star (Deheuvels et al. 2014) and the localization of the transition between the slow rotating upper layers and the faster rotating core (Di Mauro et al. 2018). The signature of the location of this transition is closely linked to the pulsation properties of the star and the extent of the gravity and pressure cavities. Determining precisely the location of the transition in rotation and discriminating between various forms of the internal rotation profile of these red giants is crucial to constraining the physical nature of the missing transport mechanism.

In this study we investigate the potential of parametric rotation profiles, coupled to a Markov chain Monte Carlo (MCMC) algorithm to determine the properties of the internal rotation profile of the red giant star Kepler-56 (KIC 6448890) studied by Huber et al. (2013). In the case of Kepler-56, we see a good example of strong coupling between the two pulsation cavities. This case has been recently studied by Takata (2016a,b), who derived the asymptotic properties of mixed modes in this regime. Previous studies referring to this case can be found in Mosser et al. (2012), Goupil et al. (2013). The case of strong coupling is characterized by one p-mode coupled to a whole range of g-modes of neighbouring frequencies.

As a good depiction of the internal structure is required for a detailed analysis of the internal rotation, we start by carrying out a detailed structural modelling of the star, presented in Sect. 2. In Sect. 3 we describe the MCMC analysis and demonstrate its robustness on synthetic data, while the detailed analysis of Kepler-56 is carried out in Sect. 4 and the consequences of our results for the physical nature of the missing transport mechanism in stellar interiors are discussed in Sect. 4.2. The conclusion and the potential of our approach for other targets is discussed in Sect. 5.

2. Stellar models and properties

To be able to study the internal rotation of Kepler-56, we first need to determine a reliable model of its internal structure, so that the first-order perturbation analysis of the detected rotational splittings is meaningful. To this end, we combine global and local minimization techniques using seismic and non-seismic constraints and mean-density inversions to obtain a robust model of our target.

2.1. Observational constraints

The star was modelled based on asteroseismic and classical constraints as in Huber et al. (2013). We summarize in Table 1 the global seismic indexes as well as the classical constraints for Kepler-56, while a table with the full seismic data adopted is given in Table A.1.

The luminosity presented in Table 1 was computed from the following formula:

(1)

Here m_λ, BC_λ, and A_λ are the magnitude, bolometric correction, and extinction in a given band λ. We use the 2MASS K-band magnitude properties. The bolometric correction is estimated using the code written by Casagrande & VandenBerg (2014, 2018), and the extinction is inferred with the Green et al. (2018) dust map. A value of M_bol, ⊙ = 4.75 is adopted for the solar bolometric magnitude. Using Gaia DR3 (Gaia Collaboration 2021), we obtain a luminosity value of L = 8.602 ± 0.363 L_⊙. This value is preliminary, however, as we could not take into account the systematic correction of the Gaia DR3 catalogue, as well as for a correction of the zero-point parallax offset.

2.2. Stellar modelling

The stellar evolutionary models used in this study and the corresponding oscillation frequencies were computed with the Code Liégeois d’Evolution Stellaire (CLES, Scuflaire et al. 2008a) and the Liège OScillation Code (LOSC, Scuflaire et al. 2008b).

Our modelling is divided into three distinct phases. First, we carry out the modelling using classical constraints and the individual radial frequencies using the AIMS software (Reese 2016; Rendle et al. 2019). Second, this initial step is used to determine an inverted mean density value, which is then used in the third modelling step using a Levenberg–Marquardt minimization technique (Moré 1978) and aiming at reproducing the whole oscillation spectrum of Kepler-56, computing the evolutionary models on the fly rather than using a predetermined grid of models coupled with an interpolation method, as in AIMS.

This last modelling step is a requirement to make sure that the properties of the mixed modes are well reproduced and to avoid non-linearities in the variational formulations used for the analysis of the rotational splittings of these modes. Looking at Eq. (4), we can see that the integral relation between the rotational splittings and the rotation profile depends on the kernel function. In practice, this function also depends on the properties of the eigenfunction. In the case of mixed modes, a small shift in frequency can be associated with a large variation of the eigenfunction, and can lead to biased estimates of the internal rotation (see discussion in Sect. 3.1).

The grid properties for the AIMS modelling are summarized in Table 2. This grid has been tailored to encompass the results for the mass and radius given by the seismic scaling relations and the surface metallicity given in Table 1. We used the AGSS09 (Asplund et al. 2009) solar abundances with a corrected abundance of Ne/O denoted AGSS09Ne (Landi & Testa 2015; Young 2018) for AIMS, while the Levenberg–Marquardt modelling is using the classical AGSS09 abundances¹. For both modelling steps, the models use the FreeEOS equation of state (Irwin 2012), the OPAL opacities (Iglesias & Rogers 1996), the T(τ) relation from Model-C of Vernazza et al. (1981) for the atmosphere, the mixing length theory of convection implemented as in Cox & Giuli (1968), and the nuclear reaction rates of Adelberger et al. (2011). The constraints and free parameters of both the AIMS and Levenberg–Marquardt modelling steps are summarized in Table 2. In this table X₀ denotes the initial hydrogen mass fraction, Z₀ the initial metal mass fraction, [Fe/H] the observed metallicity, L the observed luminosity, and M the mass. The convection is controlled by two parameters, the classical mixing-length parameter α_MLT and α_over characterizing the length of the core overshooting region with the same formalism as the MLT. The temperature gradient in the overshooting region is assumed adiabatic in both modelling steps. The mixing-length parameter was kept at a solar calibrated value for all tracks of the AIMS grid. While fitting the individual radial modes, we took into account the impact of surface effects by using the two-term surface correction of Ball & Gizon (2014) in AIMS. Uniform priors on the width of the grids were used for the MCMC modelling with AIMS. We provide the grid properties in Table 2.

The inversion procedure of the mean density is then used to derive a model-independent mass interval that is used as a safeguard for the following modelling steps. In practice, as this inversion is only computed using radial modes following the approach of Buldgen (2019), we do not need to worry about potential non-linearities at this stage. The results for the inverted mean density and the model-independent mass interval derived from Gaia and spectroscopic data is given in Table 3. The model-independent mass interval was derived by combining the inverted mean density to the stellar radius obtained though the expression of the luminosity of a black body using the Gaia luminosity and spectroscopic stellar effective temperature. As we can see, the modelling results lie well within this mass interval and confirm the robustness of our procedure.

In the last modelling step we used the constraints and free parameters listed in the right column of Table 2. The modelling follows a similar approach to that introduced by Deheuvels & Michel (2011), but here the same number of free parameters are used as constraints. Before attempting to reproduce the dipolar mixed modes, a run of the Levenberg–Marquardt minimization aiming at reproducing the individual frequency ratios was carried out, serving as an intermediate step. These ratios are sensitive to the stellar mass for a given mean density (Montalbán et al. 2010), and are defined as

(2)

where ⟨Δν⟩ is the average large frequency separation deduced from a linear fit of the radial modes (using the definitions in Reese et al. 2012). Once a good agreement is obtained for the frequency ratios, we turn our attention to reproducing the g-dominated dipolar frequency of the avoided crossing of lowest frequency. Specifically, here we ensure that the g-dominated mode of frequency 192.402 μHz, denoted ν_cross in Table 2, is reproduced well by the model. The final verification step consists in the analysis of the Echelle diagram for the model, which ensures that we reproduce with a sufficiently good accuracy the global oscillation pattern.

The fit was carried out starting from various mass values to better explore the parameter space and we determined the uncertainties in our modelling from an analysis of a χ² map where we considered the constraints given in Col. 2 of Table 2 and the individual frequency ratios r₀₂ in the evaluation of the χ². The map was carried out using a grid of five different masses, three metallicities, three initial hydrogen mass fraction, and three α_MLT, all centred around the value found for our best fitting model. The models were limited in age by the Gaia luminosity given in Table 1. A proxy for the age (the mass of the helium core) was used to avoid problems linked to the extreme age differences between tracks caused by the initial hydrogen, metallicity, or mass. Finally, an interpolation in mass and helium core mass was made to produce the χ² map shown in Fig. 1.

Fig. 1. Mapping of the parameter space around the model carried out with the Levenberg–Marquardt technique. The black circle represents the position of the best fitting model found with the Levenberg–Marquardt minimization technique.

A valley of minima can clearly be seen in Fig. 1 with our model as the overall minimum. A second clear minimum can be seen in the models with masses of 1.38 M_⊙, but the χ² value for our optimal model is two orders of magnitude lower and appears as a clear minimum in the valley observed in the parameter space. This study confirms the robustness of our stellar model and of the determined stellar parameters.

2.3. Stellar properties

Ultimately, the final solution of our structural modelling is given in Table 3, while the Hertzsprung-Russell diagram of the star modelled with the Levenberg–Marquardt minimization technique is presented in Fig. 2. The agreement obtained in terms of pulsation frequencies is illustrated by the Echelle diagram shown in Fig. 3 for the whole pulsation spectrum.

Fig. 2. Hertzsprung-Russell diagram of the evolution of Kepler-56 from the pre-MS to the modelled actual state of the star. The grey area corresponds to the constraints on the effective temperature and luminosity listed in Table 1.

Fig. 3. Echelle diagram of Kepler-56 comparing observed oscillation frequencies (listed in Table A.1) represented as squared points and theoretical oscillation frequencies represented as circle points. The difference between observed and modelled data is quantified, in the legend, as a χ2 for each spherical order ℓ. Finally, the error bars on the observed frequencies are included in the data points and not clearly distinguishable here.

The classical stellar properties obtained with the Levenberg–Marquardt method are within a 1σ difference with respect to the values presented by Huber et al. (2013). A direct comparison of the models of the two studies using a χ² function is irrelevant here since different constraints were used in the fit. However, we can compare the agreement in terms of individual frequencies for both the model computed with AIMS, fitting only the radial modes, and the model computed with the Levenberg–Marquardt technique, taking into account non-radial oscillations. Unsurprisingly, we then see that the latter has a much better χ² than the former, by an order of magnitude for the ℓ = 2 modes and by a factor 3 for the ℓ = 1 modes.

We note that the Echelle diagram shown in Fig. 3 includes the empirical surface corrections of Ball & Gizon (2014) for the ℓ = 0 and the p-dominated ℓ = 1 modes. These corrections were calibrated from acoustic oscillations of main-sequence stars. By fitting only the lowest frequency modes we ensured a limited impact of the surface effect correction on the final results, while avoiding additional parameters for the modelling. A comparison with Fig. S8 of Huber et al. (2013) shows that our solution is as good as theirs, given that they included surface corrections in their direct modelling of the individual frequencies. Overall, this model shows good agreement in individual frequencies and in frequency ratios r₀₂.

We also see that our optimal mass value agrees well with the model-independent interval defined using the mean density inversion. Our global parameters agree within one sigma with the results of the detailed modelling of Huber et al. (2013), despite the use of the different stellar evolution codes (Ventura et al. 2008) and reference solar abundances (Grevesse & Sauval 1998). This lends us further confidence in the robustness of our results, and that the agreement with the individual frequencies is suitable enough to carry out the analysis of the rotational properties of the mixed modes.

3. MCMC analysis

In this section we detail our procedure regarding the MCMC analysis of the rotation profile of post-main-sequence stars. In Sect. 3.1 we introduce the theoretical brackgound of our approach and the basis equations applied in the MCMC analysis. In Sect. 3.2 we detail our methodology for our MCMC analysis. Finally, in Sect. 3.3 we prove the robustness of our method by testing the recovery of synthetic pre-computed rotation profiles from their associated rotational splittings.

3.1. Rotational splittings and rotational kernels

In the case of a slow rotator, the effects of rotation can be treated as a perturbation of the non-rotating spherically symmetric equilibrium state. To the first order the frequencies of the oscillations will be given by

(3)

where is the frequency of the modes taking into account rotation and is the frequency of the non-rotating state. The difference between and is called the rotational splitting (Ledoux 1951). It provides a direct way to determine the internal rotation of stars in asteroseismology. Rotational splittings depend on two quantities, first the rotation profile of the star Ω(r) and then the so-called rotational kernel K_n, ℓ(r). The behaviour of the kernel is directly linked to the physical properties of the oscillation mode with which it is associated. Assuming the rotational profile to be spherically symmetric, a rotational splitting can be expressed mathematically as

(4)

In practice there is a similar integral expression for each observed rotational splitting. The information they provide is however degenerate, as it is intimately bound to the nature of the rotation kernel. It is also worth noting that they must satisfy mathematical inequalities (Reese 2015), and thus each splitting is not mathematically independent of the others.

Regarding these constraints, Kepler-56 is an excellent target since ten splittings were observed and are presented in Table A.1. Each splitting carries information about the rotation profile of the star. However, due to the extreme asymmetry of the splitting at 205.437 μHz, we exclude it from our first-order analyses.

An important aspect of the modelling of RGB stars is the strong non-linear behaviour of the modes as a slight shift in frequency may induce a large change in the eigenfunctions, and thus a complete misunderstanding of the information carried by the rotation kernels. The non-linearity is controlled in our modelling with direct reproduction of the cavities coupling through the values ν_cross and the verification of the Echelle diagram. Ultimately, the non-linearity can be illustrated by looking directly at the amplitude of kernels, as shown in Fig. 4 where kernels of a given mode are compared for two different models, one built only fitting radial oscillation modes and the other fitting the whole oscillation spectrum. The clear difference of amplitude can be seen to directly impact their behaviour, illustrating perfectly the structural dependency in the splittings and the importance of taking into account constraints on dipolar mixed modes to obtain a reliable fit of the internal structure before going through a thorough rotational analysis of a red giant star.

Fig. 4. Kernels of the same gravity dominated modes from the best fitting model found with AIMS (in orange) fitting the radial modes only and the model found with the Levenberg–Marquardt minimization technique (in blue) taking into account constraints on dipolar mixed modes.

The splittings of the whole ℓ = 1 spectrum is illustrated in Appendix B for the models obtained with AIMS and the Levenberg–Marquardt minimization technique for a given rotation profile. This figure shows the impact of differences in the rotational kernels, induced by slight differences in the stellar models, directly impacting the splittings. Appendix B illustrates perfectly the need to reproduce the whole oscillation spectrum, here by the means of the Levenberg–Marquardt minimization technique to avoid biased inferences due to an improper reproduction of the oscillation cavities. Finally, this figure shows that the use of the AIMS model could induce a strong variation in the derived surface and core rotation, in this case with a slower surface rotation and a faster rotation expected in the core.

3.2. The MCMC method

Markov chain Monte Carlo methods are a class of algorithms that sample a probability distribution. The result from such an algorithm is, ideally, the posterior probability distribution of the variables considered. To perform our analysis we used the emcee package for python 3 (Foreman-Mackey et al. 2013) with the Affine Invariant MCMC Ensemble sampler and the parallel tempering approach implemented in the package ptemcee (Vousden et al. 2016; Foreman-Mackey et al. 2013).

The MCMC analysis was based on Eq. (4). Using LOSC and InversionKit (Reese & Zharkov 2016), we extracted the kernels of individual pulsation modes for the observed splittings presented in Table A.1.

Parametric rotation profiles were then injected in Eq. (4) to get an estimate of the splittings. The role of the MCMC analysis is to optimize the free parameters of the rotation profiles to reproduce the observations. To do this, each step of the MCMC is decomposed as follows:

• For each walker, splittings are computed following Eq. (4) using the previously extracted kernels and a rotation profile computed with a given set of parameters.

• The splittings obtained by integration are then compared to the observed ones by means of the logarithm of the likelihood ℒ:

  (5)

  with

  (6)

  where σ_δνobs are the uncertainties on the observed splittings.

• New parameters of the rotation profile are then adopted, based on the new value of lnℒ, and the loop starts back from its first step until the number of iterations requested is achieved.

The results obtained at the end correspond to the median of the parameter distributions; the errors correspond to the first 15.9% and 84.1% of the distributions. The viability of the MCMC approach is verified through the computation of the autocorrelation for each free parameter and an analysis of the walkers trajectories. The walkers correspond to chains of stochastic processes randomly sampling the probability distribution according to the algorithm.

Since we are probing a highly degenerate parameter space, we implemented a parallel tempering approach with ptemcee and emcee to treat these complex problems. We based the number of iterations N on the observed autocorrelation of each sample to guarantee that the autocorrelation of all the MCMC results presented below are at least around N/50, as recommended by Foreman-Mackey et al. (2013) to have acceptable sampling. The number of walkers, temperatures, and the ladder of temperature were chosen in order to reduce the autocorrelation and increase the overall quality of the probability distribution obtained, thus reducing the required number of iterations and the MCMC running time. In the following section we detail the technical characteristics of each MCMC run.

3.3. Confirmation and verification of the method

Before applying this method to the actual target, Kepler-56, its robustness had to be investigated. To do so, we carried out tests on artificial data using the following methodology. First, the viability of the integrator is checked by comparing the splittings found by our integrator to the one computed with InversionKit for a given rotation profile. Synthetic splittings are then created with a rotation profile for a given set of parameters. The MCMC is run as described in Sect. 3.2. The objective of this set-up is to prove that the MCMC is capable of recovering the parameters of the rotation profile used to build the synthetic splittings.

The results obtained on such tests depend on the parametric rotation profile used and especially the number of input free parameters. Assuming a large number of free parameters may well lead to degeneracies as the information on the rotation profile given by the splittings is limited. Thus, the degree of customization of the synthetic rotation profiles is intrinsically limited by the number of observed splittings and by the nature of the modes for which they are observed. In our analysis we tailor the tests to the case of the dataset of Kepler-56, meaning that we use the exact same modes with their observational uncertainties to prove that the method is viable for this particular target.

The parametric profiles used find some basis in more physical analyses of AM transport processes, while trying to keep the number of free parameters as limited as possible. For example, we introduce the simple definition

(7)

where r_BCE is the radius at the base of the convective envelope, and α and Ω_core are the two free parameters of this profile. In the following analysis, we refer to this rotation profile as the power-law profile.

In this case the rotation is assumed as a solid body in the radiative zone and differential in the convective envelope, following a functional dependency in r^−α. The particular case of α = 1 corresponds to the prescription of Kissin & Thompson (2015) and Takahashi & Langer (2021) regarding the expected rotation profile in stars with large-scale fossil magnetic fields.

For the MCMC analysis, this profile shows no degeneracies, thus the Affine Invariant MCMC Ensemble sampler of emcee samples the parameter space extremely efficiently. Each MCMC analysis with the power-law rotation profile uses 20 walkers for 5000 iterations and a burn-in of 200 iterations. The priors are uniform for both α and Ω_core. The free parameter Ω_core is assumed to be positive. Depending on the cases studied α is limited differently; during the MCMC check-up α is explored form 0 to 10 and for the study of Kepler-56 the limits of α are detailed in Sect. 4.1.

A prescription aimed at mimicking the profiles of Eggenberger et al. (2012) is also part of our investigation through the following rotation profile:

(8)

Here r_norm is an arbitrary normalization constant (r_norm = 0.00804), and Ω_core, Ω_surf, and σ are the free parameters. In our analysis we refer to this rotation profile as the Gaussian profile. This profile predicts a differential rotation in the radiative zone with a solid-body rotation in the envelope due to the highly efficient transport of AM by convection.

For the MCMC analysis of the Gaussian profile, we used uniform priors on each free parameter. Due to the highly degenerate parameters space and the highly multimodal expected posterior distribution we used the parallel tempering approach of ptemcee with 40 walkers, 8 temperatures, 2000 iterations and 100 steps of burn-in for each results of the Gaussian rotation profile presented below. All the positive values of the surface rotation Ω_surf and (Ω_core − Ω_surf) were explored to study all the possible scenarios, and in particular discard solid-body rotation. Finally, with σ we explored a transition in the central area of the star by setting its minimum to 0.01 to avoid divergences, and its maximum to 4 expecting a transition in the rotation profile close the chemical composition gradient following Eggenberger et al. (2012). We also do not expect to have any constraints after σ = 4 because the splittings of ℓ = 1 modes are not affected by the rotation of this region, only observations of higher degrees splittings can give hints of the rotation for the rest of the radiative zone.

The position of the transition from the rapid internal to the slow surface rotation is of particular interest for this profile. We expect an abrupt transition, due to the important chemical composition gradient located at the peak of the Brunt–Vaïsälä frequency close to the hydrogen-burning shell. This behaviour is expected from a variety of physical mechanisms, bound to a turbulent nature, such as magnetic instabilities (Eggenberger et al. 2019c) and internal gravity waves (Pinçon et al. 2017) that would be inhibited by the effects of mean molecular weights gradients.

In this context we tried to characterize the slope of the transition with several rotation profiles. Finally, we present a last simple rotation parametrization as

(9)

where Ω_core, Ω_surf, and r_TR are the free parameters. In our analysis we refer to this rotation profile as the step profile. This profile models the transition from core to surface rotation as a discontinuity to ensure the expected abrupt transition.

For the MCMC analysis of the step rotation profile, uniform priors are used for all free parameters. As in the Gaussian rotation profile, the parameter space shows a strong degeneracy between the parameters, thus we use the parallel tempering implemented by ptemcee to properly sample the prior distribution. For each MCMC run with the step rotation profile we used 40 walkers, 8 different temperatures, and 2000 iterations with a burn-in of 100 steps. All the positive values Ω_surf and Ω_core were explored. Finally, r_TR is explored from 0.001 to 0.03 R_⋆ expecting a transition in the rotation profile close the chemical composition gradient following Eggenberger et al. (2012).

The objective with the last two rotation profiles is to define whether the slope of the transition should be close to a discontinuity or much smoother. Other profiles were tested, but either included too many free parameters (which led to degeneracies in the solutions) or were affected by poor convergence. As a result we decided to limit ourselves to the profiles given by Eqs. (7)–(9). For the chosen rotation profiles the solution from the MCMC check-up and the rotation profile they should recover is shown in Fig. 5.

Fig. 5. Illustration of the different rotation profiles found with the MCMC analysis (in color), their uncertainties as a shaded area and the rotation profiles they should replicate in black.

The MCMC check-up shows no degeneracy for the parameters of the power-law function and recovers easily such a profile. The step and Gaussian function are able to recover the input rotation profile within their error bars even if a degeneracy appears between the parameters controlling the transition and the core rotation. For all the rotation profiles tested, the surface rotation in particular is always found without any degeneracy, and thus this quantity is very well constrained by our method.

With the step and Gaussian rotation profile the limit of the method presented here can already be seen for a star without observed ℓ = 2 splittings like Kepler-56. We could expect for both rotation profiles a relationship or a correlation between the position of the transition and the core rotation. However, as illustrated in the middle left panel of both Figs. 6 and C.1 representing the posterior distribution in check-up for respectively the step and Gaussian rotation profile, after σ = 1.5 or r_TR = 0.012 R_⋆ the position of the transition and the core rotation are close to being independent parameters. This independence is probably caused by a lack of sensitivity in the kernels after σ = 1.5 or r_TR = 0.012 R_⋆; the amplitude of the kernels is not sufficient to impact the splittings very much, thus this area of the parameter space is not constrained and the MCMC can freely sample this region finding the same solution, independently from the position of the transition. The solution obtained after σ = 1.5 or r_TR = 0.012 R_⋆ is at an almost constant Ω_core, meaning that a major peak in the distribution is obtained around this value, as seen in the middle panel of Figs. 6 and C.1. This peak is a feature of the uniform sampling made in the region with the lack of resolution, and thus it has no real physical justification.

Fig. 6. Parameter probability distribution obtained with the check-up of MCMC analysis on the Gaussian function. The red solid line represents the median of the distribution, the red dotted lines correspond to the 1-σ uncertainties assuming a Gaussian distribution. The solid green line represent the maximum of the likelihood solution while the solid orange correspond to the initial input parameters that the check-up should recover. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

For the Gaussian rotation profile a second peak with a lower amplitude located around the maximum likelihood and initial input parameters can be seen in the middle panel of Fig. 6. This peak and its equivalent in the distribution of σ shows that the input parameters can be recovered at the limit of the one sigma error bars if the exploration of the parameter space is dominated by a region where the kernels have no resolving power.

The posterior distribution of the step rotation profile is very similar to the posterior distribution of Gaussian rotation profile. The input parameters are also recovered within the errors bars even if the distributions shows signs of multimodality especially for Ω_core. The multimodality observed in the posterior distribution can be attributed to the behaviour of the rotation kernels and the shape of the rotation profile. Namely, the sharp, discontinuous transition in the rotation profile will be placed at the location where it has a maximum impact on the splittings, near the local maximum of the kernels.

As mentioned before, the surface rotation is well constrained by our method; the posterior distribution of the surface rotation, illustrated in the bottom right panel of Figs. 6 and C.1, shows one Gaussian-type peak centred on the input parameters.

4. Kepler-56

The MCMC analysis was run with the splittings presented in Table A.1. We detail the results obtained for each rotation profile in Sect. 4.1. We then discuss the implications of this analysis in the context of a missing AM transport process during the evolution of the star in Sect. 4.2.

4.1. Results from the MCMC analysis

As a starting point to the analysis we carried out a standard inversion procedure using the Substractive Optimally Localized Averages (SOLA, Pijpers & Thompson 1994) method to derive the core rotation of Kepler-56, to have a measurement that is independent from the MCMC analysis. We derived an internal rotation in the core of Kepler-56 of 1100 ± 100 nHz from the SOLA inversion. This result will be compared to the final result of the MCMC analysis, keeping in mind that both techniques follow very different approaches and provide complementary views of the inversion problem. Before presenting the actual results of the MCMC analysis some additional information should be given on the rotation profiles and the final probability distributions of their parameters.

The power-law profiles were separated into two cases, a first MCMC analysis (called limited power law) was done limiting the free parameter α to 1.5 following the prescription of Kissin & Thompson (2015), Takahashi & Langer (2021), which assumes that α lies between 1 and 1.5. The second MCMC analysis on this profile (called unlimited power law) leaves all the free parameters unlimited in their values. In both cases the final distribution of parameters shows no degeneracy and a well-defined solution, the parameter probability distribution of the unlimited power-law function, is illustrated in Fig. C.1. The best fitting MCMC results and the median of the distribution are presented in Table 4 for both rotation profiles.

As mentioned in Sect. 3.3, the Gaussian rotation profile shows a correlation between the parameters controlling the position of the transition and the core rotation. The posterior probability distribution obtained using Kepler-56 data (Fig. 7) shows the same general behaviour as the check-up. A clear peak in the probability of the parameter controlling the core rotation can be seen, while the position of the transition is poorly constrained. The final parameters values obtained from the MCMC analysis are presented in Table 4.

Fig. 7. Triangle diagram presenting the posterior distribution from the MCMC analysis on the Gaussian function with the data of Kepler-56. The red solid line represents the median, the red dotted lines it’s one-sigma uncertainties and the green solid line the maximum of the likelihood solution. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

The last rotation profile studied is the step rotation profile. The posterior parameter probability distributions obtained with Kepler-56 data are illustrated in Fig. 8. The distributions show a strong degeneracy with several pairs of parameters maximizing the likelihood. As expected, the surface rotation is efficiently constrained independently from the core rotation or the position of the transition.

Fig. 8. Triangle diagram presenting the result from the MCMC analysis on the step function with the data of Kepler-56. The probability distribution of the position of the transition is cut at 0.030 R⋆ due to negligible amplitude of kernels after this mark. The red solid line represents the median, the red dotted lines it’s one-sigma uncertainties and the green solid line the maximum of the likelihood solution. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

For the sake of completeness, we also compare the values of the Bayesian Information Criterion (BIC) indicator, defined as

(10)

where k is the number of free parameters, n the number of constraints, and the maximum of the likelihood.

The distribution obtained for the limited power law shows that no clear minimum is found in the parameter space allowed. This rules out this profile as a potential solution for the case of Kepler-56; we showed with synthetic data that it should be possible to retrieve such a profile from the available dataset of Kepler-56. For the unlimited case a minimum is found with a χ² = 5.98. In this case the minimization algorithm tends to compensate the advanced position of the transition by making a sharper rotation profile in the convective envelope. The resulting rotation profile is physically unlikely due to the extremely slow rotation found at the surface and the strong dependency in ∼r^−4.55 of the rotation in the convective zone, in disagreement with a local conservation of AM. All modelled splittings for each of the parametric profiles studied are presented and compared to the observed splittings in Fig. 9.

Fig. 9. Splittings observed and presented in Table A.1 in black and compared to the modelled splittings with the different rotation profiles obtained with the median of the MCMC final distribution of parameters.

On the other hand, the step and Gaussian rotation profiles reproduce with more accuracy the observed splittings with respectively χ² = 2.64 and χ² = 2.72 and a ΔBIC ≈ 40 compared to the power-law rotation profile. Despite the clear χ² and BIC advantage of these profiles, the introduction of a third free parameter controlling the position of the transition also leads to a degeneracy between the free parameters as illustrated in Figs. 7 and 8. Even with this degeneracy, major peaks in the parameter probability distribution came out of the noise.

Looking at the posterior probabilities, we observe a behaviour similar to that in our check-ups using artificial data. This leads us to believe that a strong degeneracy exists between the various parameters of the profiles we tested using our MCMC approach. This resolution limit is a direct result of the behaviour of the rotation kernels, which show low amplitudes above 0.03 R_⋆, meaning that the rotation profile of Kepler-56 cannot be constrained from the rotation splittings above this limit without additional independent constraints.

In the specific case of Kepler-56 the surface rotation determined from starspots could have played such a role. However, we find it to be very precisely constrained by the MCMC technique, at a mean value of nHz. We note that the Ω_surf found with our MCMC analysis corresponds to the rotation in the acoustic cavity. The presence of a shear layer at the surface of the star or a breaking of the spherical symmetry in the rotation profile in the envelope can impact what we deem to be surface rotation, which should be considered when comparing measurements from other techniques. Huber et al. (2013) reported a stellar surface rotation with the value of 156 ± 6 nHz. Using 156 nHz as input for our analysis we were unable to reproduce the ratio of the p to g dominated dipolar modes, as already noted by Klion & Quataert (2017). However, while they could reproduce the minimum ratio of the splittings with a surface rotation of 78 nHz, we find that we require an even slower rotation. The origin of these differences could be found in the fact that we carried out a detailed structural modelling of Kepler-56 and analysed the entire oscillation spectrum, while they restricted themselves to the minimum ratio of the g to p dominated splittings. Huber et al. (2013) discussed the value they obtained, stating that the period is very close to the duration of the observation quarters of the Kepler spacecraft, and that the value given could be a harmonic of the actual rotational frequency. As noted above, our measurement cannot be directly compared to that of Huber et al. (2013) as other effects, such as a breaking of spherical symmetry of rotation in the acoustic cavity or a surface shear layer, could reduce the disagreements with the value obtained from starspots by Huber et al. (2013). The strong discrepancy between core and surface rotation (seen in the ratios of the splittings of p-dominated to g-dominated modes, as discussed in Eggenberger et al. 2012 and visible in our case in Fig. 9) can also support this result.

Meanwhile, our MCMC technique provides a core rotation value that is in agreement with the SOLA inversion results and with the value provided by Klion & Quataert (2017). The large number of g dominated modes strongly constrains this parameter to a narrow range of values in our modelling. Had the surface rotation not been so widely different, we could have hoped to more precisely locate the position of the transition in rotation in the radiative zone. The observed splittings and those modelled with the Gaussian and step rotation profiles are illustrated in Fig. 9.

We note that all the observed splittings illustrated in Fig. 9 are reproduced within their one sigma uncertainty by both the step and the Gaussian rotation profile. The robustness of our approach with respect to issues in the observational dataset was also tested. We modified the splitting at ν = 221.464 nHz of about 5% and checked whether the results obtained with our method were impacted by this modification. We found our method to be less sensitive to such issues, unlike classical inversion techniques such as the SOLA method (Pijpers & Thompson 1994; Pijpers 1997). The final rotation profiles associated with the splittings presented in Fig. 9 are illustrated in Fig. 10.

Fig. 10. Three different profiles obtained with the median of the MCMC final distribution of parameters compared to the Brünt–Väisälä frequency in black. The dotted lines represent the best MCMC solutions for each rotation profile with the Kepler 56 data.

4.2. Impact on the missing AM transport process

Characterizing the rotation profile of evolved stars is of course of prime importance to determine the physical nature of the efficient AM transport process needed to correctly reproduce the asteroseismic measurements available for these stars. The present results on the internal rotation of Kepler-56 show that a rigid rotation in the radiative interior followed by a power-law transition in the convective envelope, as theorized by Kissin & Thompson (2015), Takahashi & Langer (2021), is incompatible with the seismic observations. Our findings seem to favour a transition possibly located in the deep radiative layers. While our results discard large-scale magnetic fields imposing a rotation as prescribed in Kissin & Thompson (2015), Takahashi & Langer (2021), they tend to favour AM transport by magnetic instabilities, which are strongly affected by the presence of steep chemical composition gradients (e.g. den Hartogh et al. 2019; Eggenberger et al. 2019b), here located at the hydrogen-burning shell. Our parametric profiles, which mimick the results of evolutionary computations using such processes, reproduce very well the observed rotational splittings of Kepler-56. However, the quality of the data does not allow us to very precisely locate the transition in the radiative zone in the rotation profile as a result of the uncertainties on the surface rotation measured from starspots.

5. Conclusion

In this study we carried out a detailed analysis of the properties of Kepler-56, a well-known early red giant branch exoplanet-host star observed by Kepler (Huber et al. 2013). Our goal was to study the internal rotation of the star using seismic data and to see if we could provide constraints on the physical nature of the missing AM transport mechanism acting inside post-main-sequence stars. To do so we needed to discriminate between rotational profiles showing a solid-body rotation in the whole radiative region from profiles showing a transition in rotation in these layers.

We started by carrying out an extensive seismic modelling procedure of the structure of Kepler-56. We combined global minimization techniques with the AIMS software (Reese 2016; Rendle et al. 2019), local minimization following the approach of Deheuvels & Michel (2011), and seismic inversions following Reese et al. (2012), Buldgen (2019). This led us to obtain an excellent agreement with the seismic and non-seismic constraints and allowed us to consider that the modelling was robust enough to carry out an extensive analysis of the rotational properties of the star. We started by testing our methodology on synthetic data with the exact same number and type of oscillation modes as for Kepler-56, and the same uncertainties on the rotational splittings. The results of this analysis proved the robustness of the method in the considered case, and gave us confidence in applying our technique to the real target. We then showed that we are able to discriminate between various types of rotation profiles, and rejected the hypothesis of solid-body rotation in the whole radiative zone of Kepler-56 followed by a slowly decreasing power-law profile in radius, as theorized by Kissin & Thompson (2015), Takahashi & Langer (2021). Our analysis shows that a transition in the rotation profile located in the radiative region, close to the hydrogen-burning shell is favoured. Parametric models including such a sharp transition close to the shell are able to reproduce the observed splittings of Kepler-56 with excellent agreement. The stellar envelope rotation is efficiently constrained to a value of nHz, which is a factor of 6 lower than the surface rotation suggested by starspot measurements reported by Huber et al. (2013). A shear layer at the surface of the star, or a breaking of the spherical symmetry of the rotation profile in the convective envelope, as in the Sun, could potentially reduce this discrepancy.

By ruling out the expected profiles in the case where large-scale magnetic fields ensure the angular momentum transport, we do not support their presence as potential solutions for the current observed discrepancies between theoretical and observed rotation properties of red giants. However, in peculiar cases large-scale magnetic fields could play a role in the transport of the angular momentum, but there appears to be at least another process present. This process could be found in magnetic instabilities, which would lead to sharp transitions in the internal rotation profile located in regions of steep chemical composition gradients inside the star.

The main limitation of the present method is that it requires very high-quality data, a large number of observed splittings, and an extensive modelling of the stellar structure. Moreover, our method is limited by the parametric description of the rotation profiles inside the star. However, it is definitely applicable to some of the best Kepler targets (see e.g. Di Mauro et al. 2016, 2018), and also potentially to TESS targets with long observation durations, some of which having even better datasets than Kepler-56.

Acknowledgments

G. B. acknowledges funding from the SNF AMBIZIONE grant No. 185805 (Seismic inversions and modelling of transport processes in stars). P. E. and S. J. A. J. S. have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 833925, project STAREX). A. M. and J. M. acknowledge support from the European Research Council (ERC grant agreement No. 772293 for the project ASTEROCHRONOMETRY). This article used an adapted version of InversionKit, a software developed within the HELAS and SPACEINN networks, funded by the European Commissions’s Sixth and Seventh Framework Programmes.

References

Appendix A: Seismic data

Appendix B: Illustration of the non-linear behaviour in the stellar modelling

Fig. B.1. Splittings obtained with the two different models presented in Sect. 2, in blue the model obtained with Levenberg–Marquardt minimization technique and in green, the model obtained with AIMS. Splittings were computed, here, with a given rotation profile, thus, the differences seen between the two models can be attributed to the rotational kernel, linked to the structure of the star.

Appendix C: Additional figures

Fig. C.1. Parameter probability distribution obtained with the check-up of MCMC analysis on the step function. The red solid line represents the median of the distribution, the red dotted lines correspond to the 1 − σ uncertainties assuming a Gaussian distribution. The solid green line represent the maximum of the likelihood solution while the solid orange correspond to the initial input parameters that the check-up should recover. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

Fig. C.2. Triangle diagram presenting the result from the MCMC analysis on the power law function with the data of Kepler-56. The probability distribution of α is cut at α = 8 due to a low amplitude noise impacting the computation of the median. The red solid line represents the median, the red dotted lines its one-sigma uncertainties and the green solid line the best MCMC solution.

All Tables

All Figures

	Fig. 1. Mapping of the parameter space around the model carried out with the Levenberg–Marquardt technique. The black circle represents the position of the best fitting model found with the Levenberg–Marquardt minimization technique.
In the text

	Fig. 2. Hertzsprung-Russell diagram of the evolution of Kepler-56 from the pre-MS to the modelled actual state of the star. The grey area corresponds to the constraints on the effective temperature and luminosity listed in Table 1.
In the text

Fig. 3. Echelle diagram of Kepler-56 comparing observed oscillation frequencies (listed in Table A.1) represented as squared points and theoretical oscillation frequencies represented as circle points. The difference between observed and modelled data is quantified, in the legend, as a χ2 for each spherical order ℓ. Finally, the error bars on the observed frequencies are included in the data points and not clearly distinguishable here.

In the text

	Fig. 4. Kernels of the same gravity dominated modes from the best fitting model found with AIMS (in orange) fitting the radial modes only and the model found with the Levenberg–Marquardt minimization technique (in blue) taking into account constraints on dipolar mixed modes.
In the text

	Fig. 5. Illustration of the different rotation profiles found with the MCMC analysis (in color), their uncertainties as a shaded area and the rotation profiles they should replicate in black.
In the text

Fig. 6. Parameter probability distribution obtained with the check-up of MCMC analysis on the Gaussian function. The red solid line represents the median of the distribution, the red dotted lines correspond to the 1-σ uncertainties assuming a Gaussian distribution. The solid green line represent the maximum of the likelihood solution while the solid orange correspond to the initial input parameters that the check-up should recover. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

In the text

Fig. 7. Triangle diagram presenting the posterior distribution from the MCMC analysis on the Gaussian function with the data of Kepler-56. The red solid line represents the median, the red dotted lines it’s one-sigma uncertainties and the green solid line the maximum of the likelihood solution. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

In the text

Fig. 8. Triangle diagram presenting the result from the MCMC analysis on the step function with the data of Kepler-56. The probability distribution of the position of the transition is cut at 0.030 R⋆ due to negligible amplitude of kernels after this mark. The red solid line represents the median, the red dotted lines it’s one-sigma uncertainties and the green solid line the maximum of the likelihood solution. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

In the text

	Fig. 9. Splittings observed and presented in Table A.1 in black and compared to the modelled splittings with the different rotation profiles obtained with the median of the MCMC final distribution of parameters.
In the text

	Fig. 10. Three different profiles obtained with the median of the MCMC final distribution of parameters compared to the Brünt–Väisälä frequency in black. The dotted lines represent the best MCMC solutions for each rotation profile with the Kepler 56 data.
In the text

	Fig. B.1. Splittings obtained with the two different models presented in Sect. 2, in blue the model obtained with Levenberg–Marquardt minimization technique and in green, the model obtained with AIMS. Splittings were computed, here, with a given rotation profile, thus, the differences seen between the two models can be attributed to the rotational kernel, linked to the structure of the star.
In the text

Fig. C.1. Parameter probability distribution obtained with the check-up of MCMC analysis on the step function. The red solid line represents the median of the distribution, the red dotted lines correspond to the 1 − σ uncertainties assuming a Gaussian distribution. The solid green line represent the maximum of the likelihood solution while the solid orange correspond to the initial input parameters that the check-up should recover. Only a small fraction of the sampler is displayed to avoid overweighted figures, however, the histograms were constructed with the full sample.

In the text

	Fig. C.2. Triangle diagram presenting the result from the MCMC analysis on the power law function with the data of Kepler-56. The probability distribution of α is cut at α = 8 due to a low amplitude noise impacting the computation of the median. The red solid line represents the median, the red dotted lines its one-sigma uncertainties and the green solid line the best MCMC solution.
In the text