© A. Noll et al. 2021

1. Introduction

One of the most important current open questions in stellar physics is the extent of convective cores. Several physical processes are known to extend the convective core boundaries beyond the standard Schwarzschild limit. The most frequently quoted are overshooting of ascending blobs of fluids due to their inertia, rotational mixing or semi-convection. All these processes remain poorly described by theory, and the way they interact is understood even less. They are therefore generally modeled, in stellar evolution codes, as an extension of the mixed core over a distance d_ov, which is often referred to as the distance of overshoot, even though other processes can contribute as well. In practice, this can either be achieved by simply extending the fully mixed central region, or by treating overshooting as a diffusive process, following a behavior found in numerical simulations (Freytag et al. 1996). In both cases, a free parameter controlling the extent of the additional mixing is required. Observations are therefore necessary to better constrain those processes.

Initial constraints have been obtained thanks to the study of the HR diagram of clusters (see e.g., Maeder & Mermilliod 1981, VandenBerg et al. 2006), and the modeling of eclipsing binaries (Claret & Torres 2016). Most of those studies favor adding overshooting, to various extents. Typically, d_ov is around 0.2 H_p, where H_p is the pressure scale height. Claret & Torres (2016) found that α_ov, the ratio between d_ov and H_p, increases with mass for stars under 2 M_⊙ before reaching a plateau. However, this result is still debated (Constantino & Baraffe 2018, Claret & Torres 2019).

Over the last decade, asteroseismology allowed us to probe the structure of stellar cores. Thanks to the data of CoRoT (Baglin et al. 2006), Kepler (Borucki et al. 2010) and now TESS (Ricker et al. 2014) missions, we have been able to precisely measure the oscillation frequencies of numerous pulsators. The study of pressure (p) modes in low-mass main sequence (MS) stars, showed the need for core overshooting to correctly reproduce the observed frequencies (Goupil et al. 2011, Deheuvels et al. 2010, Silva Aguirre et al. 2013). Deheuvels et al. (2016), modeling several MS stars, found that α_ov increases with the mass. Moreover, gravity (g) mode pulsators, like slowly-pulsating B (SPB) stars, are interesting targets to constrain the additional mixing around convective cores. Indeed, gravity modes probe the inner chemical structure of the star and allow detailed investigation of the convective core extensions. Moravveji et al. (2015, 2016), when modeling SPB stars, found that overshoot was necessary, and they favored diffusive overshooting over a simple extension of the central mixed region.

Post-main-sequence stars are another way to put constraints on the amount of overshooting. Once the central hydrogen is exhausted, nuclear energy production stops, leaving an inert radiative helium core. This core then contracts, heating the surrounding hydrogen layers of the star until shell burning starts. At that moment, the star begins its subgiant phase, and evolves on a nuclear timescale for masses below about 1.5 solar masses (M_⊙). For stars that are close to the terminal-age main sequence (TAMS), the structure and the evolution remain highly influenced by the properties of the MS convective core. Interestingly, the star begins to exhibit mixed modes at that moment. These modes behave like gravity modes in the internal propagation cavity and pressure modes in the outer one. Thus, they allow us to finely probe the deepest layers of the star, all the while being observable. This and the proximity of the subgiant to the TAMS make the mixed modes of young subgiants valuable data in studying the extension of convective cores.

Another particularity of mixed modes is their very fast evolution, compared to the nuclear evolution timescale of the subgiant phase. Indeed, mixed mode frequencies change dramatically over the course of a few million years. This makes their seismic modeling challenging. Recently, increasing efforts have been made to model subgiants (Huber et al. 2019; Stokholm et al. 2019; Metcalfe et al. 2020; Deheuvels et al. 2020; Li et al. 2019, 2020), driven by both their great physical interest and the sudden increase of seismic data for these stars. Most of those works focused on finding the optimal stellar parameters for one or several subgiants. So far, few studies have used subgiants as tools to test stellar physics, mainly due to the challenges of their modeling, as mentioned above.

Deheuvels & Michel (2011) successfully constrained α_ov from a subgiant observed by CoRoT, HD 49385, which exhibits only one g-dominated mode and is therefore very close to the TAMS. They found that either no overshooting, or a model with α_ov = 0.19 were giving equally good results. In this work, we modeled a young subgiant, KIC10273246, which was observed by Kepler over almost 1000 days. That star exhibits two g-dominated modes, which allows us to better constrain its inner structure. We performed a thorough seismic modeling of the star, in order to precisely estimate its stellar parameters and to place constraints on the extension of its MS convective core.

In Sect. 2, we show the utility of having access to two g-dominated mixed modes in young subgiants. In Sect. 3, we present the surface observables of KIC10273246 and perform a fresh analysis of its oscillation spectrum using the full Kepler data set. We then describe, in Sect. 4, the modeling technique that we adopted, which is an improved version of the method developed by Deheuvels & Michel (2011). Section 5 presents our optimal stellar models and the constraints that were obtained from the extent of the MS convective core for KIC10273246. We discuss these results in Sect. 6, and Sect. 7 is dedicated to our conclusions.

2. Probing potential of mixed modes

Just after the main sequence, the oscillation spectra of solar-like pulsators show the presence of mixed modes, which are due to the coupling between the observed p-modes and low radial-order g-modes (n_g = 1, 2, 3, n_g being the number of nodes in the g-mode cavity). The frequency spacing between low-order g-modes is large (several times the large separation of p modes), so that only a few are in the observable frequency window during the subgiant phase. Moreover, with n_g being low, the pure g modes that couple to p modes do not follow an asymptotic behavior (as described in Shibahashi 1979, Tassoul 1980). The oscillation spectra of subgiants therefore constrast with those of more evolved stars, which typically have more g-dominated modes than p-dominated modes, and for which n_g is of the order of several tens (e.g., Mosser et al. 2012).

The frequencies of mixed modes are mostly determined by two structural features of the star. The first is the g-mode (G) cavity, which is delimited by the Brunt-Väisälä frequency N. The second is the evanescent zone between the g-mode and p-mode (P) cavities, the latter being delimited by the Lamb frequency S_l.

The G cavity affects the frequency of the g-mode that is involved in the mixed mode frequency. G-mode frequencies, in the asymptotic theory, can be approximated by

$$\begin{array}{c}\hfill {\nu }_{n,l}\approx \frac{\sqrt{l(l+1)}}{2{\pi }^2(n-1/2)}{\displaystyle {\int }_{r_1}^{r_2}\frac{N}{r}\mathrm{d}r,}\end{array}$$(1)

l being the degree of the mode, r₁ and r₂ the turning points in the G cavity, and r the local radius of the star. In our case, n_g is low for the observed modes, so the asymptotic expression given in Eq. (1) should not apply. However, it has been shown that it can provide qualitative information about the behavior of the mixed mode frequencies (Deheuvels et al. 2010). It tells us that g-dominated modes should give strong constraints on the Brunt-Väisälä frequency in the G cavity. One can write it in the following form (e.g., Kippenhahn et al. 2012):

$$\begin{array}{c}\hfill N^2=\frac{g\delta }{H_p}({\mathrm{\nabla }}_{\mathrm{ad}}-\mathrm{\nabla }+\frac{\varphi }{\delta }{\mathrm{\nabla }}_{\mu }),\end{array}$$(2)

where g is the local gravity, δ = −(∂lnρ/∂lnT)_P, μ, ϕ = (∂lnρ/∂lnμ)_P, T, ∇_ad = (∂lnT/∂lnP)_ad, ∇ = ∂lnT/∂lnP, and ∇_μ = ∂lnμ/∂lnP. The Brunt-Väisälä frequency consequently carries information about both the thermal (two first terms in parentheses) and compositional structure (last term) of the star.

The evanescent zone affects the coupling between the two cavities, whose strength is linked to the size of this region and the value of N inside it (e.g., Unno et al. 1989). Using a toy model, Deheuvels & Michel (2011) showed that a strong coupling induces a shift of the l ≥ 1 p-dominated frequencies that are close to a g-mode. The p-dominated frequencies therefore provide complementary information about the internal structure of the subgiant.

In this work, we investigated whether having several g-dominated modes in the observed oscillation spectrum could offer more information regarding the extension of the MS core. From above, we know that the frequencies of the g-dominated mixed modes are related to the N/r integral between the turning points of the G cavity. Figure 1 shows the propagation diagram of a subgiant close to the TAMS, highlighting the frequencies of the first two g-dominated l = 1 mixed modes, that is, those that arise due to the coupling of p modes with g modes of radial orders n_g = 1 and 2. We denote their turning points in the G cavity as r_i1, r_o1 for the mode with n_g = 1, and r_i2, r_o2 for the mode with n_g = 2. The difference between the two frequencies is thus mainly related to the Brunt-Väisälä frequency value between r_o1 and r_o2 (as the one in the [r_i1, r_i2] region is negligible). This region, as it can be seen in Fig. 1, is dominated by the μ-gradient contribution. This gradient is related to the characteristics of the hydrogen-burning shell, especially the nuclear energy generation, and thus its temperature and composition. It has been shown that a H-burning shell structure depends on the MS core features, and especially on the amount of core overshooting. One can see this in Fig. 5 of Deheuvels & Michel (2010), which exhibits two Brunt-Väisälä profiles of stars having the same evolutionary stage and position in the HR diagram, but computed with different α_ov. The two profiles differ mainly in the peak caused by the μ-gradient, and these structural differences are large enough to have a significant impact on the eigenfrequencies of the star.

Fig. 1. Typical propagation diagram of a low-mass subgiant star. The Brunt-Väisälä frequency is represented in blue, delimiting the g-mode cavity (light blue). The Lamb frequency, in orange, delimits the p-mode cavity (light orange). Two g-dominated mixed mode angular frequencies, with ng = 1, 2, are represented (solid lines in propagation zones, dotted lines in evanescent zones). The G cavity turning points are noted as ri1, ri2 and ro1, ro2. Finally, the thermal and chemical contributions to the Brunt-Väisälä frequency are represented in green (dashed) and red (dot-dashed), respectively.

For all those reasons, stars with two visible g-dominated modes are therefore expected to be interesting targets on which to place constraints on the efficiency of core overshooting. That criterion led us to the choice of KIC10273246, a subgiant with two g-dominated modes and a mass of 1.26 ± 0.10 M_⊙ (Creevey et al. 2012), which places it safely in the mass range where stars are expected to have a convective core during the MS.

3. Observational properties of KIC 10273246

3.1. Surface constraints

3.1.1. Constraints from spectroscopy

Surface constraints were obtained for KIC10273246 by Creevey et al. (2012). The authors used different algorithms on the same spectra obtained with the FIES spectrograph. For our target, they found effective temperatures (T_eff) ranging from 5933 to 6380 K. A weighted mean gives us a value of 6150 ± 100 K, which we have used to constrain our seismic modeling. The star was also found to have to have a sub-solar metallicity, with [Fe/H] = −0.13 ± 0.1 dex.

3.1.2. Constraints from broadband photometry

To obtain a reliable value of the luminosity of the star, we performed a spectral energy distribution (SED) fit, following the procedure of Stassun & Torres (2016). We extracted photometric values using the VizieR photometry tool. Those data come from the NUV filter from GALEX (Martin et al. 2005), the B_T and V_T filters from Tycho-2 (Høg et al. 2000), the J,H and K_s filters from 2MASS (Skrutskie et al. 2006), the gri filters from SDSS (Skrutskie et al. 2006), the W1-W4 filters from WISE (Wright et al. 2010), and the G magnitude from Gaia (Evans et al. 2018). The atmosphere model comes from the Kurucz atmosphere grid (Kurucz 2005), with the surface gravity (log g) derived from the seismic relations (see Sect. 3.2.4), and the metallicity coming from spectroscopic measurements. We then fit the photometry points to the spectrum, with the T_eff and extinction A_v as free parameters. We also used the spectroscopic data from Creevey et al. (2012) and the extinction from Green et al. (2019) as priors. With a reduced χ² of 0.7, we found T_eff = 6000 ± 33 K, and $A_{\mathit{v}}=0.{00}_{-0.000}^{+0.037}\mathrm{mag}$. The fit spectrum and the photometric data are represented in Fig. 2. Finally, we integrated the flux over all the wavelengths and used the distance from Gaia to obtain the luminosity of the star. According to Zinn et al. (2019), a parallax bias exists in the Kepler field, which depends on the G-band magnitude and the pseudo-color ν_eff (effective wavenumber of the photon flux distribution in the Gaia band) of the star. We found ϖ − ϖ_Gaia = 39.15 ± 9.46 μas, which gives L = 5.74 ± 0.17 L_⊙. This result is, as expected, lower than the Gaia archive value (5.92 ± 0.13 L_⊙) due to the parallax offset.

Fig. 2. Best fit of the SED using Kurucz atmosphere models. The orange points represent the observations with the corresponding error bars, and the blue curve represents the best fit SED model.

3.2. Seismic constraints

3.2.1. Preparation of Kepler light curve

The subgiant KIC10273246 was observed with Kepler between quarters Q0 and Q11 (total duration of 978 days) in short cadence (58.85 s). An early seismic analysis of the target was performed by Campante et al. (2011) using the first four quarters of Kepler observations (325 days of data). They estimated the frequencies of oscillation modes of degrees l = 0, 1, and 2 over eight overtones. We revisited this analysis using the complete Kepler data set. The light curve of the star was processed using the Kepler pipeline developed by Jenkins et al. (2010). Corrections from outliers, occasional drifts and jumps were performed following the method of García et al. (2011). The power density spectrum (PSD) was then obtained by applying the Lomb-Scargle periodogram (Lomb 1976, Scargle 1982).

The PSD is shown in the shape of an échelle diagram in Fig. 3. We recall that the échelle diagram is built by dividing the PSD in consecutive chunks with a length corresponding to the large separation of acoustic modes Δν, and piling them up. Here, we used the estimate of Δν = 48.2 μHz obtained by Campante et al. (2011). The main interest of échelle diagrams is that acoustic modes of the same degree align in nearly straight ridges, which eases mode identification. The neighboring l = 0 and l = 2 ridges are readily identified on the left part of the échelle diagram (modes indicated by crosses and triangles, respectively, in Fig. 3).

Fig. 3. Échelle diagram of KIC10273246, folded with Δν = 48.2 μHz. For clarity, the spectrum was smoothed over a 0.2-μHz boxcar. The white symbols indicate the frequencies that have been extracted for modes of degree l = 0 (crosses), l = 1 (diamonds), and l = 2 (triangles) in Sect. 3.2.2. The two l = 1 g-dominated mixed modes are circled in red.

The ridge of l = 1 modes (indicated by diamonds in Fig. 3) deviates from the nearly-vertical line that is expected for purely acoustic modes. This behavior is known to arise for dipolar modes in the presence of avoided crossings between low-order g modes and p modes in subgiants. Each avoided crossing is characterized by the presence of an additional mode, which lies away from the ridge of theoretical purely acoustic l = 1 modes (which would be a nearly vertical line at an abscissa of about 35 μHz in Fig. 3). This mode is most strongly trapped in the core and is thus g-dominated. The neighboring l = 1 modes are p-dominated, but their frequencies are nevertheless affected by the presence of the g-dominated mode. The modes with frequencies larger than the g-dominated mode are shifted to higher frequencies (to the right in the échelle diagram) and those with frequencies below the g-dominated mode are shifted to lower frequencies (to the left in the échelle diagram). These features are clearly visible in Fig. 3, corresponding to two l = 1 avoided crossings. The l = 1 g-dominated modes associated to these avoided crossings are circled in red in Fig. 3.

3.2.2. Extraction of oscillation mode parameters

To extract the parameters of the oscillation modes, we followed the method of Appourchaux et al. (2008). Here, we briefly recall the main steps of the procedure and refer the reader to that paper for more details.

Prior to fitting the individual oscillation modes, we modeled the background of the PSD. The contribution from granulation was modeled by two Harvey-like profiles, following the prescription of Karoff et al. (2013), and we added a white noise component to account for photon noise. The overall contribution from the oscillations was modeled as a Gaussian function. We fit this model to the PSD using maximum-likelihood estimation (MLE). The central frequency of the Gaussian function gives an estimate of the frequency of maximum power of the oscillations ν_max. To determine the error on this quantity, we subdivided the Kepler light curve in ten chunks of equal duration, and fit the background model on the PSD calculated with these time series. The error on ν_max was taken as the standard deviation of the measurements of this quantity for each chunk. We thus obtained ν_max = 843 ± 20 μHz. The PSD was then divided by the optimal background model.

We then performed a fit of the oscillation modes, which were modeled as Lorentzian functions to account for their finite lifetimes. Each mode profile of degree l, radial order n, and azimuthal order m was characterized by its central frequency ν_n, l, m, its height H_n, l, m and its line width Γ_n, l, m. Since dipolar modes have a mixed character, it cannot be assumed that they share similar heights and line widths with the neighboring radial modes, as is often done for main sequence solar-like pulsators. Most quadrupolar modes are expected to be p-dominated, owing to the weak coupling between the P and G cavities for these modes. We therefore assumed that the l = 2 modes have the same heights and widths as their closest l = 0 modes, with the exception of one g-dominated l = 2 mode, which is discussed below and in Sect. 3.2.3. Non-radial modes are split into multiplets by rotation. Owing to the slow rotation of subgiants, the effects of rotation on the mode frequencies can be found by applying a first-order perturbation analysis. The components of a rotational multiplet are thus expected to be equally spaced by the rotational splittings δν_n, l. We also assumed that they share similar line widths, and that their height ratios depend only on the inclination angle i of the star following the expressions given by Gizon & Solanki (2003). In principle, mixed modes can have different rotational splittings, because they probe the rotation at different depths in the star. This has been used to probe the internal rotation of subgiants (e.g., Deheuvels et al. 2014).

To test whether individual rotational splittings can be measured in KIC10273246, we first performed local fits of the non-radial modes. Following the method described by Deheuvels et al. (2015), we fit each mode using two different models: one under the H₀ hypothesis (no rotation, so that each mode is modeled as a single Lorentzian), and one under the H₁ hypothesis (rotation is considered and each mode is modeled as a set of 2l + 1 Lorentzians separated by the rotational splitting). It is clear that hypothesis H₁ necessarily provides better fits to the data than hypothesis H₀ since it involves two additional free parameters (inclination angle and rotational splitting). The significance of hypothesis H₁ can be tested using the likelihoods ℓ₀ and ℓ₁ of the best fits obtained under the H₀ and H₁ hypotheses, respectively. As shown by Wilks (1938), the quantity ΔΛ ≡ 2(lnℓ₁ − lnℓ₀) follows the distribution of a χ² with Δn degrees of freedom, where Δn is the difference between the number of free parameters involved in hypotheses H₁ and H₀ (here, Δn = 2)¹. For each multiplet, we thus obtained a value of ΔΛ. The false-alarm probability was then given by the p-value p = P(χ²(2 d.o.f.)≥ΔΛ), which corresponds to the probability that a mode under the null hypothesis can produce such a high value of ΔΛ.

For dipolar modes, the lowest p-value that we found is 0.08, which is too high to consider the measurement as significant. This means that we cannot reliably extract individual rotational splittings for dipolar modes in this star. The most likely explanation is that the modes have large line widths compared to the rotational splitting. For quadrupolar modes, only one mode (the one with a frequency around 779.4 μHz) was found to have a low p-value, of about 4 × 10⁻⁵, which shows a very high significance level. A rotational splitting of 0.53 ± 0.03 μHz was obtained for this mode (see Fig. 4). This mode is in fact a g-dominated mixed mode, as we show in Sect. 3.2.3.

Fig. 4. Oscillation spectrum of KIC102773246 in the vicinity of a quadrupolar mode that was found to be significantly split by rotation (see Sect. 3.2.2). The thick red curve corresponds to our best-fit model of the spectrum. Two quadrupolar mixed modes are visible (around 779.4 μHz and 783.9 μHz) and one radial mode (around 785.6 μHz).

We then performed a global fit of the modes (all the modes are fit simultaneously). Since individual splittings cannot be measured, we assumed a common rotational splitting for all l = 1 and l = 2 modes (except for the aforementioned l = 2 mode around 779.4 μHz). Since most non-radial modes are p-dominated, we expect the common rotational splitting to essentially measure the rotation in the envelope. The best fit corresponds to a rotational splitting of δν = 0.45 ± 0.02 μHz for non-radial modes and an inclination angle of i = 55 ± 6°. As was done for local fits, we also performed an additional fit of the modes without including the effects of rotation (null hypothesis). We could therefore estimate the p-value corresponding to the measurement of a mean rotational splitting. We found p ∼ 10⁻⁴, which indicates a high level of significance. Our results are compatible with the estimates of Campante et al. (2011), who had found i ≳ 20° for this star, and optimal values of the rotational splitting slightly below 0.5 μHz.

The best-fit parameters for the oscillation modes (frequencies, heights, and line widths) are given in Table A.1. The uncertainties of the fit dipolar mode frequencies range from 0.08 to 0.50 μHz. The measured mode frequencies are in quite good agreement with the ones found by Campante et al. (2011). Discrepancies at the level of 3σ were found for only two modes (the dipole mode around 1055 μHz and the quadrupole mode around 880 μHz). Using the complete Kepler data set enabled us to detect l = 0 and l = 2 modes over three additional radial overtones compared to Campante et al. (2011). Our results are also in very good agreement with the recent measurements of mode frequencies for KIC10273246 by Li et al. (2020) using the complete Kepler data set (agreement at the level of 2σ or better for all oscillation modes).

3.2.3. Detection of an l = 2 mixed mode

We mentioned above that the l = 2 mode with a frequency of about 779.4 μHz is the only mode for which an individual rotational splitting could be measured. This mode also has other distinctive features. It is separated from the closest radial mode by 6.1 ± 0.2 μHz. By comparison, for the other radial orders, the average separation between the l = 2 mode and the neighboring l = 0 mode is 4.4 μHz, with a standard deviation of 0.4 μHz. This suggests that this mode might be an l = 2 mixed mode, the frequency of which is modified by the coupling between the p- and g-mode cavities. This hypothesis is strengthened by the fact that it has a short line width (0.26 ± 0.08 μHz) compared to the width of the neighboring l = 2 modes (between 1.7 and 2.4 μHz). Indeed, if the mode under study is a g-dominated mixed mode, it should have a higher inertia than p-dominated l = 2 modes, and therefore a shorter line width. Figure 4 shows the profile of the radial mode that is closest to the l = 2 mode under study. There appears to be an additional mode in the left wing of the radial mode, at a frequency of about 783.9 μHz. To determine the significance of this mode, we performed local fits assuming either its presence (H₁ hypothesis) or absence (H₀ hypothesis). We found a p-value of 0.01, indicating a reasonable significance level. This also supports the identification of the l = 2 mode at 779.4 μHz as a mixed mode. In this case, the additional mode at 783.9 μHz would also be an l = 2 mixed mode undergoing an avoided crossing with its close neighbor. As is shown in Sect. 5, the best-fit models for KIC10273246 do show a pair of mixed modes in the vicinity of these two modes, which confirms our identification.

3.2.4. First estimates of stellar parameters using seismic scaling relations

To obtain first estimates of the global stellar parameters of the star, we used seismic scaling relations, which relate the global seismic parameters Δν and ν_max to stellar properties such as the mass, radius and surface gravity (Brown et al. 1991). These relations could be derived because ν_max scales to an equally good approximation as the acoustic cut-off frequency (Brown et al. 1991; Stello et al. 2008; Belkacem et al. 2011).

To estimate the asymptotic large separation of acoustic modes, we followed the prescription of Mosser et al. (2013). We fit an expression of the type

$$\begin{array}{c}\hfill {\nu }_{n,0}=[n+\frac{\alpha }{2}{(n-n_{\max })}^2+{\epsilon }_p]\mathrm{\Delta }{\nu }_{\mathrm{obs}}\end{array}$$(3)

to the observed radial modes, where Δν_obs is the observed large separation around ν_max, α measures the curvature corresponding the to the second-order term in the asymptotic development, ε_p is an offset, and n_max = ν_max/Δν_obs. We thus obtained Δν_obs = 48.47 ± 0.02 μHz, which translates into an asymptotic large separation of Δν_as = 50.63 ± 0.02 μHz, following Mosser et al. (2013).

Using our estimates of Δν_as, ν_max from Sect. 3.2.2, and T_eff from Sect. 3.1, we could apply seismic scaling relations to derive preliminary estimates of the star’s mass, radius, and surface gravity. We obtained M = 1.24 ± 0.12 M_⊙, R = 2.10 ± 0.07 R_⊙, and log g = 3.88 ± 0.03.

4. Seismic modeling method

4.1. Physics of the models

We used MESA v10108 (Paxton et al. 2015) evolution models, with OPAL equation of states and opacity tables (Rogers et al. 1996, Iglesias & Rogers 1996), with the solar mixture from Asplund et al. (2009). The models were computed with an Eddington-gray atmosphere. The convection regions were treated using the standard mixing length theory (MLT) as prescribed in Cox & Giuli (1968), with a free parameter α_conv corresponding to the ratio between the mixing length and the pressure scale height. Microscopic diffusion was taken into account, unless otherwise specified, using the Burgers formalism (Burgers 1969) and diffusion coefficients from Stanton & Murillo (2016). However, radiative accelerations have not been included in the computed models, as the increase in computational time could not be afforded in this study. The impact of those processes are discussed in Sect. 6.2.

As Gabriel et al. (2014) stated, for stars that have a growing convective core, it is necessary to use the Ledoux criterion to determine the radius of the convective core R_cc. This way, we avoid the creation of unphysical convective zones outside the core in strongly chemically stratified regions, which may have an impact on the composition profile of the star and thus on its evolution. Moreover, we used the predictive mixing scheme (Paxton et al. 2018).

Core overshooting was modeled as a step extension of the convective core, over a distance

$$\begin{array}{c}\hfill d_{\mathrm{ov}}={\alpha }_{\mathrm{ov}}min(H_p,R_{\mathrm{cc}}/{\alpha }_{\mathrm{conv}}),\end{array}$$(4)

where d_ov is the distance of instant mixing overshooting, H_p the pressure scale height, and α_ov a free parameter quantifying the phenomenon. Equation (4) replaces the traditional expression d_ov = α_ovH_p in order to prevent d_ov from becoming unphysically large when the core is small (H_p → ∞ when r → 0). It is important to note that this definition varies from one evolution code to another (see, e.g., Eq. (1) of Deheuvels & Michel 2011 for CESAM2K). Low-mass stars have small convective cores, therefore those differences must be kept in mind when comparing models coming from different codes. Additionally, the impact on our results of using a diffusive overshooting, as proposed by Freytag et al. (1996), is discussed in Sect. 6.1.

The adiabatic oscillations of the models were computed using ADIPLS (Christensen-Dalsgaard 2008), and the surface effects were corrected for using the cubic term of the prescription of Ball & Gizon (2014).

4.2. Why modeling subgiants is challenging

The frequencies of g-dominated mixed modes evolve over a very short time, with a non-linear change of several μHz per million years, which is much larger than the usual uncertainties coming from Kepler data. As this timescale is several orders of magnitude shorter than the typical nuclear evolution time of low-mass subgiants, reproducing the mixed modes with a traditional grid technique requires extremely small steps in mass and age. This makes this method prohibitive when the number of free parameters is large, as is required to test the model physics. Interpolation in age is possible (Li et al. 2020), but somewhat difficult for l = 2 g-dominated modes, which we observed in KIC10273246. Interpolation across tracks (as used e.g., in AIMS, Rendle et al. 2019) could mitigate the need for extremely fine steps in mass, but needs to be tested for subgiants, especially regarding the extreme sensitivity of the g-dominated frequencies to the masses of the models. Additionally, an “on-the-fly” optimization technique may perform badly due to the highly non-linear behavior of the mixed modes, especially during the computation of the derivatives in the case of a gradient-descent kind of algorithm.

To overcome those difficulties, a new approach is necessary. We thus developed a nested optimization, where we optimize the physical parameters of models (e.g., metallicity, initial helium abundance etc.) that have been optimized in mass and age beforehand. This way, we can handle those two sensitive parameters using a dedicated and specific procedure, separately from the other ones for which a more traditional technique is possible. This modeling method originates from Deheuvels & Michel (2011) and has been adapted to make it more robust. In the following, we recall the basic principles of this method and highlight the differences with the one used in the present study.

4.3. Optimization in mass and age

In that part of the optimization process, we compute models with only two free parameters, the mass and the age of the star, the rest being fixed. The optimization of those two parameters can be made easier thanks to the fact that, if all the other physical parameters (such as metallicity, mixing-length parameter...) are fixed, reproducing only Δν and the frequency ν_g of a g-mode is enough to precisely constrain the mass and the age.

A physical justification of that approach can be found in Deheuvels & Michel (2011). We remind the reader of it here using a HR-diagram represented in Fig. 5. It shows the iso-Δν line, as $L\propto T_{\text{eff}}^5$ for models with the same large separation, and the iso-ν_g line, computed from stellar evolution models. The two lines meet at a unique point, that can be reached by tuning only the mass (i.e., choosing the “right” evolution path) and age (i.e., stopping at the right moment on that path). In concrete terms, our first step is, at a given mass, to find the age that correctly reproduces the ν_g frequency.

Fig. 5. HR-diagram representing the evolution tracks of stellar models with masses varying from 1.2 M⊙ (lightest gray) to 1.3 M⊙ (darkest gray) and otherwise identical physics. Each evolution is stopped when νg is correctly reproduced.

As we only see mixed modes and not pure g-modes, we cannot directly measure ν_g. A possible solution would be to choose a g-dominated mode (i.e., a non-radial mode far from its ridge) frequency. Unfortunately, such a frequency does not evolve monotonously with age, as can be seen in the upper panel of Fig. 6. We thus preferred to look at the distance between that g-dominated mode and its closest radial mode, which we call δν.

Fig. 6. Evolution of δν (top panel) and ν1, 11 (bottom panel) with age for a 1.3 M⊙ star, after the main sequence. Here, like in our modeling of KIC10273246, δν is defined as ν1, 11 − ν0, 13, and the observed value is represented by the dotted line. The plot has been strongly magnified in order to see the 1-σ uncertainties from the data.

As we can see in the top panel of Fig. 6, this value always decreases with age, but it also keeps the interesting properties of the mixed modes as it evolves very quickly during an avoided crossing, allowing us to tightly constrain the age. This step would be equivalent to following a unique evolution path in Fig. 5 and stopping it when it crosses the iso-ν_g line.

We then optimize on those “good-age” models in order to correctly reproduce the large separation. In practice, to do this we minimize the χ² of only the radial modes, which we define as

$$\begin{array}{c}\hfill {\chi }_{\mathrm{rad}}^2={\displaystyle \sum _n\frac{{({\nu }_{0,n}^{\mathrm{mod}}-{\nu }_{0,n}^{\mathrm{obs}})}^2}{{\sigma }_{0,n}^2}.}\end{array}$$(5)

We do not take into account the non-radial modes at this stage to eliminate the complexity of behavior of the mixed modes. This approach differs from the one followed by Deheuvels & Michel (2011), who at this stage searched for models that minimized the difference between the observed average large separation and the one coming from the models. By using all the radial modes instead here, we found that the optimization process is more robust regarding the correction of near-surface effects. It may be observed that the behavior of Δν (and, in turn, of the radial frequencies) is close to linear when varying the mass. Then, a simple Newton-type algorithm (such as a Levenberg-Marquard algorithm) is enough to quickly find the optimal mass. This step would then be equivalent to the right evolution path that leads to the meeting points of the two iso-lines on Fig. 5.

Figure 7 shows the échelle diagram of a model that we can obtain after that first step, with arbitrary physical parameters: metallicity [Fe/H] = −0.2 dex, mixing-length parameter α_conv = 1.5, initial helium abundance Y₀ = 0.28. We can see that the radial modes and δν = ν_1, 11 − ν_0, 13 (the proxy for ν_g) are, by construction, correctly reproduced. However, the other frequencies are far from the observed ones. Especially, the g-dominated mode ν_1, 18 is about 10 μHz away from the observations. Thus, to find a better matching model, we adjust the other parameters.

Fig. 7. Échelle diagram of a model optimized in mass and age (open symbols) and of the observed frequencies (full, with their 3-σ error bars). The radial and dipolar modes are represented by crosses and diamonds, respectively, with their radial order indicated.

4.4. Optimizing the other parameters

Now that we have a method to correctly find the mass and the age of a star at a given physics, we must find the other parameters of the stars, this time taking into account all the observational constraints. Thus, we define a new χ² as

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=1}^{N_{\mathrm{obs}}}\frac{{(x_i^{\mathrm{obs}}-x_i^{\mathrm{mod}})}^2}{{\sigma }_i^2}=\sum _{i=1}^{N_{\mathrm{obs}}}{\mathrm{\Delta }}_i,}\end{array}$$(6)

where N_obs is the total number of observational constraints, both seismic and non-seismic, $x_i^{\text{obs}},x_i^{\text{mod}}$ the values of those observed constraints or their computed equivalent, and σ_i their respective uncertainties. We also introduced the quantities ${\mathrm{\Delta }}_i\equiv {(x_i^{\text{obs}}-x_i^{\text{mod}})}^2/{\sigma }_i^2$, which indicate the contributions of every observable to the χ², to be used later.

As those parameters have a lower impact on the frequencies than the mass and age, it is now possible to use more traditional approaches. One possibility is to compute grids of models, where each model of the grid is optimized in mass and age. Another option is to perform an optimization using an iterative method, where again each iteration consists of an optimization of the mass and age. To model KIC10273246, we opted for a hybrid method, which is described in the following section.

4.5. Fitting procedure adopted for KIC10273246

For the modeling of KIC10273246, we left the initial metallicity [Z/X]₀, the mixing-length parameter α_conv, the initial helium abundance Y₀, and, of course, the overshoot parameter α_ov as free parameters. At first, to explore the global behavior of the χ², we computed a very loose grid ([Fe/H] between −0.2 and 0.2, step 0.1; Y₀ between 0.24 and 0.28, step 0.02; α_conv between 1.5 and 2.1, step 0.2 and α_ov between 0.0 and 0.2, step 0.05). We recall that each model of this grid is optimized in mass and age as explained in Sect. 4.3. The purpose of this loose grid was to investigate whether double solutions or local minima exist. No such features have been found. Moreover, those grids allowed us to identify the region of the best parameters.

We thereafter refined those parameters. As mentioned in Sect. 4.4, the optimization of [Z/X], α_conv, Y₀, α_ov can be performed either through a grid approach or an iterative procedure. We therefore conducted several tests, using stellar models as mock observations, to determine which method is preferable. We found the best robustness when following a hybrid approach: for given values of Y₀ (0.26 through 0.31, step 0.01) and α_ov (0 through 0.25, step 0.05), we conducted iterative optimizations with the Levenberg-Marquardt algorithm to find optimal values of [Fe/H] and α_conv. This method differs from the one used in Deheuvels & Michel (2011) where a single grid was used for all the free parameters.

Among those models, we considered only those that were compatible with the observational estimates of the chemical enrichment law ΔY₀/ΔZ₀. Typical values quoted for ΔY₀/ΔZ₀ range from 1.4 to 4 (e.g., Chiappini et al. 2002, Balser 2006, Casagrande et al. 2006). We consequently had a conservative approach and took into account all models with ΔY₀/ΔZ₀ < 5.

5. Results

In this section, we describe the general characteristics of the best models, before commenting the constraints on α_ov. We finally investigate at the internal structures of the best models and the constraints brought by the mixed modes.

5.1. General characteristics of the best models

Following the method described in Sect. 4.5, we obtained optimal models for each value of α_ov, whose characteristics are in Table 1. The best model, with α_ov = 0.15, has a reduced χ² of 3.2. The échelle diagram of this model is represented in Fig. 8. Also, surface observables are consistent with the ones found in the literature, or with the SED fitting previously described: we find a less than 1-σ difference for the effective temperature T_eff, the metallicity [Fe/H], and the luminosity L. We found a radius and a mass that are consistent with the seismic scaling relations as well. We can note that, as expected from the mass-based prediction, all good models had a convective core during the MS. This supports our choice of using this star to constrain α_ov.

Fig. 8. Échelle diagram of the best model, with αov = 0.15. Symbols and colors are identical to those in Fig. 7.

We note that our best-fit models are significantly less massive and are older than those of Li et al. (2020), who also performed a seismic modeling of KIC20173246 and found M = 1.49 ± 0.08 M_⊙ and an age of 2.84 ± 0.60 Gyr. These discrepancies could be partially explained by the different assumptions made on the input physics. For instance, Li et al. (2020) considered a solar-calibrated mixing length parameter, while we left this parameter free in our study. Also, contrary to us, Li et al. (2020) neglected element diffusion and adopted the mixture of Grevesse & Sauval (1998). Finally, we stress that the agreement with the observed dipole mode frequencies, in particular for the g-dominated mode ν_1, 18, is significantly better in the present study than it is for the best-fit models of Li et al. (2020) (compare Fig. 10 of Li et al. 2020 to Fig. 8 of the present paper). These mismatches between models and observations for dipole modes are acknowledged by Li et al. (2020), and the authors attribute them to an imprecise modeling of the core structure.

For each combination of (α_ov, Y₀), our minimization using the LM algorithm can be used to derive uncertainties in the stellar parameters. The error bars in the free parameters of the fit ([Fe/H] and α_conv) are obtained as the diagonal coefficients of the inverse of the Hessian matrix. The uncertainties on the other parameters can then be obtained using Eq. (10) of Deheuvels et al. (2016). We thus obtain very small error bars, of the order of 0.007 for [Fe/H] and 0.002 for α_conv, which translates into uncertainties of approximately 0.004 M_⊙ for the stellar mass and 0.04 Gyr for the age. This means that for a given combination of (α_ov, Y₀), the available observables provide very strong constraints on the stellar parameters. By comparison, we find that optimal models with different Y₀ can yield similar agreement with the observations (statistically equivalent χ²) but have quite different stellar parameters. This degeneracy of stellar models with respect to Y₀ is addressed in more detail in Sect. 6.3. It thus seems that the uncertainties in the stellar parameters are dominated by the model degeneracy in Y₀. We thus used the optimal models with different Y₀ to estimate uncertainties in the stellar parameters. For a given α_ov, we fit a second order polynomial to the χ² curve and retrieved the interval of values corresponding to ${\chi }_{\min }^2+1$. This gave us the 1-σ uncertainties, which are reported in Table 1.

5.2. Constraints on core overshooting

Figure 9 shows the variation in the χ² of the optimal models with α_ov. We can see that adding overshooting allows us to reproduce the observed frequencies significantly better, with the χ² of the models without overshoot and with α_ov = 0.15 being 315 and 127, respectively.

Fig. 9. χ2 of the best model for every value of overshooting. The colored regions indicate the contributions to the χ2 (i.e., the sum of Δi) of the surface observables (Teff, L and [Fe/H]) and the frequencies according to their degrees.

To better understand which frequencies allow us to favor models with overshoot, we investigated the contributions of the different observables to the χ² (Δ_i in Eq. (6)). We denote the χ² contributions of the observables coming from the optimal models without overshoot as ${\mathrm{\Delta }}_i^{\text{nov}}$ and the equivalent from models with α_ov = 0.15 as ${\mathrm{\Delta }}_i^{\text{ov}}$. Figure 10 represents the differences ${\mathrm{\Delta }}_i^{\text{nov}}-{\mathrm{\Delta }}_i^{\text{ov}}$, positive values meaning that the observable is better reproduced by the α_ov = 0.15 model. As expected, we can see that the main χ² difference is due to the dipolar modes, which have a mixed behavior. However, we observe that the g-dominated modes (indicated by the dotted vertical lines) hardly contribute to distinguishing the models with and without overshooting. Both types of model fit the g-dominated frequencies well. The main contributors to the χ² differences are in fact the l = 1 p-dominated modes in the neighborhood of the g-dominated modes. As explained in Sect. 2, the frequencies of these modes are mainly influenced by the coupling between the P and the G cavities. The intensity of that coupling thus accounts for the main part of the differences between the models. We note that all those models correctly reproduce ν_1, 18, as the high sensitivity of this g-dominated mode strongly constrains the region of parameters of the models with the smallest χ². The major role played by the dipolar modes in the constraints on α_ov is also illustrated in Fig. 9, where the colored regions indicate the contributions to the χ² of the surface observables and modes depending on their degree.

Fig. 10. Difference of the χ2 contributions of the different observables, between αov = 0.0 and 0.15 models. The two dipolar g-dominated modes are represented by dotted vertical lines. (a), (b), and (c) are Teff, L, and [Fe/H], respectively.

Moreover, Fig. 9 indicates that the contribution to the χ² of the l = 2 modes hardly changes with α_ov. This was partly expected, because their evanescent zone is larger than that of the dipole modes, making the coupling between the G and P cavities weaker. Most of the detectable modes are therefore very strongly p-dominated and do not constrain the deep structure of the star, hence α_ov. Yet, one g-dominated l = 2 mode was detected (ν_2, 10, see Sect. 3.2.3). It is interesting to see that, in a similar way to the previous paragraph, its frequency is equally well reproduced by models with and without overshooting. One can see this in Fig. 10, where ${\mathrm{\Delta }}_i^{\text{nov}}-{\mathrm{\Delta }}_i^{\text{ov}}$ of that mode is less than 1-σ. On the other hand, the (2, 11) mode, whose frequency is close enough to ν_2, 10 to be influenced by the coupling, varies substantially with α_ov. Figure 10 shows a 3-σ difference, despite the high 0.65 μHz observational uncertainty, confirming the key role of the coupling in the constraint on α_ov. Interestingly, however, while the α_ov = 0.15 model better reproduces the dipolar modes, the (2, 11) mode is better fit in the model without overshooting. Nevertheless, its large observational uncertainty prevents it from being too constraining.

Finally, we notice that adding a larger amount of overshooting (α_ov > 0.15) strongly worsens the quality of the fit, placing a strong maximum limit on the value of α_ov. To better understand this behavior, we investigate the seismic constraints on the internal structure of the models in the next section.

5.3. Constraints on the internal structure from mixed modes

5.3.1. Constraints on central density

A tight link is expected between the g-dominated frequencies and the central density ρ_c. This comes from the relation between the g-mode frequencies and the Brunt-Väisälä frequency, which approximately scales as ρ_c (see Eq. (15) from Deheuvels & Michel 2011). To verify this, we investigated the constraints placed on ρ_c by the frequency of a g-dominated dipolar mode. For this purpose, we considered the values of ρ_c in the models computed in the loose grids defined in Sect. 4.5, in which models are optimized in mass and age in order to reproduce Δν and the frequency of a g-dominated mode (here ν_1, 11, as described in Sect. 4.3). We observed that, despite the wide range of parameters considered, the standard deviation of ρ_c among those models is as low as 32.4 g cm⁻³, which represents around 1% of the mean central density $\stackrel{\sim }{{\rho }_c}=2100{\mathrm{g}\mathrm{cm}}^{-3}$. This standard deviation even decreases to 10 g cm⁻³ if we only keep the 200 best models, illustrating the tight relation between ρ_c and the frequency of a g-dominated mixed mode.

This plays a role in the increase of χ² for α_ov > 0.15. To illustrate this point, we computed models with the same values of [Z/X], Y₀, and α_conv, but with different values of α_ov. Each of these models was optimized in mass and age, as described above. Figure 11 shows the evolution of ρ_c with age for those models. One can see that they all reach approximately the same value of central density, $\stackrel{\sim }{{\rho }_c}$, in accordance with the previous paragraph. Moreover, the intensity of the ρ_c jump that is due to the post-MS core contraction increases with α_ov. We can explain this by the fact that for bigger cores, the layers where hydrogen remains are more distant from the center. They are colder and require a stronger core contraction, hence a stronger jump in ρ_c, to reach the fusion temperature. Therefore, when α_ov increases, models that reach $\stackrel{\sim }{{\rho }_c}$ are closer to the TAMS. When the model gets too close to the TAMS, this impacts the whole stellar structure. Internally, the μ-gradient in the core is very much affected, because the nuclear reactions in the H-burning shell did not have enough time to smooth its shape. In the outer layers, the convective envelope, which expands during the subgiant phase, has a different size. Those processes alter the frequencies (the g- and p-dominated ones, respectively), which are thereby not compatible with the observations.

Fig. 11. Evolution of ρc with age, for models with different αov that have been optimized in mass and age in order to reproduce Δν and νg.

5.3.2. Constraints on the Brunt-Väisälä profile.

Based on Eq. (1), we expect the frequency of mixed modes to place constraints on the integral of the Brunt-Väisälä frequency in the G cavity. To investigate this, in Fig. 12 we plot the N² profiles for the models of the loose grid defined in Sect. 4.5, which all reproduce the low-frequency g-dominated mode (ν_1, 11) and Δν. We observe that reproducing both already strongly constrains the part of the Brunt-Väisälä frequency dominated by the thermal term (∇_ad − ∇), which corresponds to the most central layers (r < 0.05 R_⊙). This was expected because of the 1/r factor in the Eq. (1) integral. On the contrary, the part of N² that is dominated by the μ-gradient changes significantly within the grid.

Fig. 12. N2 profiles of the loose grid models. The left and right panel profiles are colored in relation to the difference between the models and the observations of ν1, 18 and δ, respectively. Those differences are normalized by the observational uncertainties. The blue horizontal lines represent the observed g-dominated frequencies.

We expect that part to be strongly determined by the dipolar modes. We therefore investigated the constraints brought by the two most determining seismic features (see Sect. 2): the coupling intensity and the frequency of pure g-modes. As those two are not directly measurable, we used observational proxies. The intensity of the coupling can be quantified by δ ≡ ν_1, 12 − ν_1, 11, which is from Deheuvels & Michel (2011). That value is the difference between the low-frequency g-dominated mode (ν_1, 11) and the following dipolar p-dominated mode (ν_1, 12). Thus, δ increases with the coupling. The frequency of a pure g-mode is measured through ν_1, 18, which is the high-frequency g-dominated mode. For those two values, we color-coded, in Fig. 12, the profiles of the Brunt-Väisälä and Lamb frequencies based on their agreement with the observations (left panel for δ and right panel for ν_1, 18).

One can see on the right panel that models correctly reproducing the coupling (i.e., dark profiles) have very similar H-burning shell positions (N² peak around r = 0.05 R_⊙). However, the Brunt-Väisälä profiles become more degenerate for higher r: several different profiles can reproduce δ within 1-σ. This degeneracy is lifted thanks to the high-frequency g-dominated mode: on the left panel, models closely reproducing ν_1, 18 all have similar Brunt-Väisälä profiles. This corroborates our theoretical approach in Sect. 2: the high-frequency g-dominated mode adding tight constraints on the shape of the μ-gradient. Important gains in structural constraints are therefore obtained from having a second detectable g-dominated mode.

6. Discussion

6.1. Diffusive overshooting

During this study, overshooting was modeled as a step extension of the convective core. However, Freytag et al. (1996) proposed another prescription, based on results coming from 2D simulations of A-stars and white dwarfs. Overshoot is then modeled as a diffusive process, with a coefficient D_ov exponentially decaying from the boundary of core, following the expression

$$\begin{array}{c}\hfill D_{\mathrm{ov}}=D_{\mathrm{conv}}exp[-\frac{2(r-R_{\mathrm{cc}})}{f_{\mathrm{ov}}{H}_p}],\end{array}$$(7)

with D_conv being the MLT derived coefficient taken just below R_cc, and f_ov a free parameter that tunes the length scale of the overshooting region.

In order to compare the results coming from the two types of overshooting, we first found the f_ov that is equivalent to a given value of α_ov. In order to do this, we searched for the value of f_ov that gives models reaching the TAMS at the same age as models computed with a step overshooting and α_ov = 0.15. We found f_ov = 0.01. After that, we modeled it using a method similar to the one used in Sect. 5 and compared the best model with the one computed with a step overshoot.

As we can see in Fig. 13, the differences between the frequencies and observables of the best models with step and diffusive overshoot are mainly less than 1-σ. We note the exception of the g-dominated ν_2, 10 mode, which is better reproduced by the diffusive overshoot model. However, its impact on the global χ² is counter-balanced by the generally better reproduced dipolar frequencies of the step overshoot model. Moreover, the difference between the characteristics of the two models are within the uncertainties of Table 1. Therefore, we cannot discriminate between the two kinds of modelings with the current set of data.

Fig. 13. Differences in the χ2 contributions of the observables coming from the best models with step and diffusive overshoot.

6.2. Effect of microscopic diffusion

The models presented in Sect. 5 of our study include both gravitational settling and chemical diffusion. Such processes, which happen to be necessary in order to correctly reproduce the helioseismic observations (Christensen-Dalsgaard et al. 1993), are expected to have an impact on our results for two main reasons. The first is the sinking during the main sequence of heavier elements because of the gravitational settling. This reduces the hydrogen abundance in the core and shortens the main sequence, which eventually decreases the age of models with same mean density. The second is the smoothing of the structure of the subgiant. High μ-gradient peaks, like the one produced by the withdrawal of the convective core, are strongly smoothed out by the chemical diffusion, which impacts the mixed mode frequencies. Thus, it is interesting to see how gravitational settling and chemical diffusion change the characteristics and the quality of the fit. We therefore modeled the star following a similar methodology, but without including those two processes. We found that the best model also has α_ov = 0.15, but provides a significantly worse agreement with the observations than the best diffusion model, with ${\chi }_{\text{nodiff}}^2-{\chi }_{\text{diff}}^2=71$. It is more massive (1.29 M_⊙) and older (4.85 Gyr), as was expected from the fact that heavy elements do not sink during the MS. We note that the surface observables are less well reproduced, with a best model that is too cold (5874 K, which represents 1.84 times the observational uncertainty) and has too low a luminosity (5.0 L_⊙, which represents 4.3 σ).

Moreover, similarly to what we found in Sect. 5, the quality of the fit improves as α_ov increases for α_ov ≤ 0.15. However, this is less significant (${\chi }_{{\alpha }_{\text{ov}}=0}^2-{\chi }_{{\alpha }_{\text{ov}}=0.15}^2=24$). For higher values, the quality of the fit strongly worsens, in a comparable way to what has been found with gravitational settling and chemical diffusion.

Figure 14 illustrates the differences in the Brunt-Väisälä profiles between the two best models, with and without both diffusive processes. One can see the high ∇_μ peak (at r = 0.06 R_⊙) is, as expected, smoothed out by the chemical diffusion. Otherwise, the two profiles are remarkably similar, despite the different physics of the models, which highlights the robustness of the constraints coming from the mixed modes.

Fig. 14. Brunt-Väisälä profiles of the best models with gravitational settling and chemical diffusion (blue) and without those two processes (orange).

Finally, the effects of gravitational settling are expected to be somewhat counter-balanced in the envelope by radiative accelerations, which can have a significant impact on stars with masses greater than 1.1 M_⊙ (Deal et al. 2018). However, including the effects of radiative accelerations in the models increases the computational time of stellar evolution calculation by orders of magnitude, and it could not be afforded in the present study. To test the impact of this process on our modeling, we computed a model that takes into account radiative accelerations, with the same parameters as the best model for α_ov = 0.15. We obtained slightly different large separations (Δν_rad − Δν_norad = 0.12 μHz) but very similar frequencies, once normalized by the large separation. Radiative accelerations are therefore not expected to change the conclusions of this study.

6.3. Helium degeneracy

To model KIC10273246, we initially performed optimizations with fixed α_ov and considering Y₀, α_conv and [Fe/H] as free parameters. In this case, we observed an unwanted sensitivity of the Y₀ parameter to the guess value of our optimization process. This led us to the hybrid approach described in Sect. 4.5, using optimizations with fixed values of Y₀ and varying [Fe/H], α_conv. We found that optimal models with different values of Y₀ indeed have surprisingly close frequencies, despite their wide range of mass. This is illustrated in Fig. 15, which shows the échelle diagrams of the best model with Y₀ = 0.28 (blue) and the best model with Y₀ = 0.24 (orange), both of which have α_ov = 0.15. Those models have quite different characteristics, as reported in Table 2. However, their frequencies are almost indistinguishable, despite the very small uncertainties on the mode frequencies from Kepler data. Only the g-dominated l = 2 mode allows us to slightly favor the Y₀ = 0.28 model. Such degeneracy is related to the anti-correlation between mass and Y₀, that has been observed in MS stars (see e.g., Lebreton & Goupil 2014) as well as subgiant stars (Li et al. 2020). Additionally, we note that no monotonic behavior has been found between the age and Y₀. We therefore conclude that the seismic modeling of subgiants, despite bringing strong constraints on the deep structure of the star, does not lift the degeneracy between Y₀ and the mass.

Fig. 15. Échelle diagrams of the best model (blue, Y0 = 0.28) and the best model for Y0 = 0.24 (orange). 3-σ uncertainties from the observations are represented by black bars.

6.4. Internal rotation

We mentioned in Sect. 3.2.2 that a rotational splitting of 0.53 ± 0.03 μHz could be measured with a high level of significance for the l = 2 mode at 779.4 μHz. This is obviously not enough to probe the internal rotation in detail. However, since this mode is g-dominated, it can be used to place approximate constraints on the rotation in the core of the star.

Using our best-fit model from Sect. 5, we were able to compute the rotational kernel K(r), which relates the splitting δν_s of this mode to the rotation profile Ω(r):

$$\begin{array}{c}\hfill \delta {\nu }_s={\displaystyle {\int }_0^{R}K(r)\mathrm{\Omega }(r)/(2\pi )\mathrm{d}r.}\end{array}$$(8)

This can be re-written as δν_s = K_g⟨Ω_g⟩+K_p⟨Ω_p⟩, where ⟨Ω_g⟩ and ⟨Ω_p⟩ are average rotation rates in the g- and p-mode cavities, respectively, and K_g (resp. K_p) corresponds to the integral of K(r) in the g-mode (resp. p-mode) cavities. For the l = 2 mode under study, using our best-fit stellar model we found that 84% of the kernel energy is enclosed in the g-mode cavity, which confirms that the mode is indeed g-dominated.

Campante et al. (2011) found a clear rotational modulation in the Kepler light curve of KIC10273246. They thus estimated the surface rotation rate of the star to about 0.5 μHz (rotation period of about 23 days). This value is comparable to the average rotational splitting of 0.45 ± 0.02 μHz that we obtained in this study. As mentioned in Sect. 3.2.2, this average splitting is dominated by the contribution of p-dominated modes, to the extent that it essentially measures the envelope rotation rate. Taken together, these two measurements suggest a low rotation contrast within the p-mode cavity, which is in line with the conclusions of Benomar et al. (2015) and Nielsen et al. (2015) for main-sequence solar-like pulsators.

The splitting measured for the l = 2 g-dominated mode is close to the rotation rate inferred for the envelope, which suggests a low amount of radial differential rotation in the star. If we take ⟨Ω_p⟩/(2π)≈0.45 μHz, we obtain a core rotation rate of about 0.65 μHz (rotation period of about 18 days). Clearly, more rotational splittings would be required to precisely measure the core rotation rate. However, our results indicate that KIC10273246 could be rotating nearly as a solid-body, like the two Kepler subgiants whose internal rotation profiles were recently measured (Deheuvels et al. 2020).

7. Conclusion

In this study, we performed a seismic modeling of KIC10273246, a subgiant observed by Kepler, and obtained strong constraints on the size of its MS convective core. We chose this target because it exhibits two dipolar g-dominated modes, which we expected to bring stronger constraints on the internal structure. We extracted the mode parameters from the oscillation spectrum of the star using the full Kepler data set and thus updated the mode frequencies that were previously obtained by Campante et al. (2011).

The seismic modeling of subgiants is notoriously complex. We here elaborated on the algorithm proposed by Deheuvels & Michel (2011). This method consists of a two-step approach. The purpose of the first step is to find the mass and age that match the large separation of p modes and the frequency of one g-dominated mixed mode. The second step optimizes the other free parameters ([Fe/H], Y₀, α_conv and α_ov) to reproduce the other frequencies as closely as possible. In this work, we improved this algorithm to make it more robust. This enabled us to perform a detailed seismic modeling of KIC10273246, with a particular emphasis on the investigation of the size of the MS convective core.

We found models in good agreement with the observations, with a reduced χ² of 3.2 for the best model, and with surface observables that are reproduced to within less than 1σ. One key result of this study is that models with core overshooting during the MS reproduce the observations significantly better, with an optimal value of α_ov = 0.15. For higher values of α_ov, the quality of the fit significantly worsens. We found that such models are very close to the TAMS. Their internal structure thus differs from that of the lower-α_ov solutions, and their seismic properties show strong mismatch with the observations. We tested the robustness of our conclusions by considering other choices for the input physics. No significant difference was found when modeling core overshooting as a diffusive process. Models computed without microscopic diffusion also favor models with α_ov = 0.15, albeit less significantly, and show a strong mismatch compared with the observations for higher values of α_ov. However, they yield poorer agreement with the seismic and surface observables compared to the models computed with microscopic diffusion. This study thus confirms the high potential of young subgiants with mixed modes to measure the extent of the MS convective cores.

We also investigated the information conveyed by the mixed modes about the core structure. We showed that the combined knowledge of the large separation Δν and the frequency of one g-dominated mixed mode is enough to estimate the central density ρ_c to a precision of about 1%. This helps us understand why models with a greater amount of core overshooting (α_ov > 0.15) are not compatible with the observations. Because of their larger MS convective core, they have a higher ρ_c just after the end of the MS, and they thus reach the optimal central density closer to the TAMS. We then studied the roles of the different mixed mode frequencies in determining the profile of the Brunt-Väisälä frequency inside the star. While the first g-dominated mixed mode strongly constrains the thermal part, the second one helps constrain the part of the Brunt-Väisälä frequency that is dominated by the μ-gradient. We therefore confirm that having access to two g-dominated mixed modes helps better characterize the Brunt-Väisälä profile.

Also, despite the strong constraints that were obtained on the internal structure, we noted the existence of a degeneracy between the stellar mass and the initial helium abundance Y₀. This degeneracy, which is already well known for MS stars (e.g., Lebreton & Goupil 2014), is not lifted by the mixed modes. We find that it is in fact the main source of uncertainties in the determination of the stellar parameters. This should be kept in mind when modeling subgiants. Current modeling techniques, such as traditional grid-based methods, tend to miss a significant fraction of the best-fit models because of the size of the mesh. In such conditions, the degeneracy between Y₀ and mass could be explored only partially, thus causing the uncertainties on the stellar parameters to be underestimated.

As a byproduct of this work, we obtained partial constraints on the internal rotation of KIC10273246. We were not able not measure individual rotational splittings for the dipolar mixed modes, but we obtained a splitting of 0.53 ± 0.03 μHz for the only g-dominated l = 2 mixed mode in the spectrum of the star. Interestingly, this value is close to the surface rotation rate of 0.5 μHz that was found for this star by Campante et al. (2011) using photometric data from Kepler. This suggests that this star might be rotating nearly as a solid-body, similarly to the two young subgiants recently studied by Deheuvels et al. (2020).

This work highlights the large potential of the seismic modeling of young subgiants to indirectly obtain constraints on the core structure of the star during the MS. The next step will be to use this method on a larger sample of stars drawn from the targets observed with Kepler and TESS, and therefore place quantitative constraints on the overshooting process in low-mass stars. The data from the upcoming PLATO mission (Rauer et al. 2014) will add a large amount of potential targets for this type of analysis.

Moreover, we show in this study that detecting several g-dominated dipole modes places stronger constraints on the shape of the Brunt-Väisälä profile, and therefore on the μ-gradient in the stellar core. It could thus be anticipated that more evolved subgiants, which show a larger number of g-dominated mixed modes, would be more favorable targets for our purpose. However, these stars are also further from the end of the MS, and a worry is that the chemical composition in the core might be less dependent on the properties of the MS core. We plan, therefore, to study this effect in a subsequent study.

Acknowledgments

We thank Morgan Deal for enlightening discussions about microscopic diffusion. We also thank the anonymous referee for comments that improved the clarity of this paper. We acknowledge support from from the project BEAMING ANR-18-CE31- 0001 of the French National Research Agency (ANR) and from the Centre National d’Etudes Spatiales (CNES).

References

Appendix A: Observed frequencies

All Tables

All Figures

Fig. 1. Typical propagation diagram of a low-mass subgiant star. The Brunt-Väisälä frequency is represented in blue, delimiting the g-mode cavity (light blue). The Lamb frequency, in orange, delimits the p-mode cavity (light orange). Two g-dominated mixed mode angular frequencies, with ng = 1, 2, are represented (solid lines in propagation zones, dotted lines in evanescent zones). The G cavity turning points are noted as ri1, ri2 and ro1, ro2. Finally, the thermal and chemical contributions to the Brunt-Väisälä frequency are represented in green (dashed) and red (dot-dashed), respectively.

In the text

	Fig. 2. Best fit of the SED using Kurucz atmosphere models. The orange points represent the observations with the corresponding error bars, and the blue curve represents the best fit SED model.
In the text

	Fig. 3. Échelle diagram of KIC10273246, folded with Δν = 48.2 μHz. For clarity, the spectrum was smoothed over a 0.2-μHz boxcar. The white symbols indicate the frequencies that have been extracted for modes of degree l = 0 (crosses), l = 1 (diamonds), and l = 2 (triangles) in Sect. 3.2.2. The two l = 1 g-dominated mixed modes are circled in red.
In the text

	Fig. 4. Oscillation spectrum of KIC102773246 in the vicinity of a quadrupolar mode that was found to be significantly split by rotation (see Sect. 3.2.2). The thick red curve corresponds to our best-fit model of the spectrum. Two quadrupolar mixed modes are visible (around 779.4 μHz and 783.9 μHz) and one radial mode (around 785.6 μHz).
In the text

	Fig. 5. HR-diagram representing the evolution tracks of stellar models with masses varying from 1.2 M⊙ (lightest gray) to 1.3 M⊙ (darkest gray) and otherwise identical physics. Each evolution is stopped when νg is correctly reproduced.
In the text

	Fig. 6. Evolution of δν (top panel) and ν1, 11 (bottom panel) with age for a 1.3 M⊙ star, after the main sequence. Here, like in our modeling of KIC10273246, δν is defined as ν1, 11 − ν0, 13, and the observed value is represented by the dotted line. The plot has been strongly magnified in order to see the 1-σ uncertainties from the data.
In the text

	Fig. 7. Échelle diagram of a model optimized in mass and age (open symbols) and of the observed frequencies (full, with their 3-σ error bars). The radial and dipolar modes are represented by crosses and diamonds, respectively, with their radial order indicated.
In the text

	Fig. 8. Échelle diagram of the best model, with αov = 0.15. Symbols and colors are identical to those in Fig. 7.
In the text

	Fig. 9. χ2 of the best model for every value of overshooting. The colored regions indicate the contributions to the χ2 (i.e., the sum of Δi) of the surface observables (Teff, L and [Fe/H]) and the frequencies according to their degrees.
In the text

	Fig. 10. Difference of the χ2 contributions of the different observables, between αov = 0.0 and 0.15 models. The two dipolar g-dominated modes are represented by dotted vertical lines. (a), (b), and (c) are Teff, L, and [Fe/H], respectively.
In the text

	Fig. 11. Evolution of ρc with age, for models with different αov that have been optimized in mass and age in order to reproduce Δν and νg.
In the text

	Fig. 12. N2 profiles of the loose grid models. The left and right panel profiles are colored in relation to the difference between the models and the observations of ν1, 18 and δ, respectively. Those differences are normalized by the observational uncertainties. The blue horizontal lines represent the observed g-dominated frequencies.
In the text

	Fig. 13. Differences in the χ2 contributions of the observables coming from the best models with step and diffusive overshoot.
In the text

	Fig. 14. Brunt-Väisälä profiles of the best models with gravitational settling and chemical diffusion (blue) and without those two processes (orange).
In the text

	Fig. 15. Échelle diagrams of the best model (blue, Y0 = 0.28) and the best model for Y0 = 0.24 (orange). 3-σ uncertainties from the observations are represented by black bars.
In the text