Fossil Signatures of Main-sequence Convective Core Overshoot Estimated through Asteroseismic Analyses

Some physical processes that occur during a star's main-sequence evolution also affect its post-main-sequence evolution. It is well known that stars with masses above approximately 1.1 M ⊙ have well-mixed convective cores on the main sequence; however, the structure of the star in the neighborhood of the convective core regions is currently underconstrained. We use asteroseismology to study the properties of the stellar core, in particular convective boundary mixing through convective overshoot, in such intermediate-mass stars. These core regions are poorly constrained by the acoustic (p) mode oscillations observed for cool main-sequence stars. Consequently, we seek fossil signatures of main-sequence core properties during the subgiant and early first-ascent red giant phases of evolution. During these stages of stellar evolution, modes of mixed character that sample the deep interior can be observed. These modes sample the parts of the stars that are affected by the main-sequence structure of these regions. We model the global and near-core properties of 62 subgiant and early first-ascent red giant branch stars observed by the Kepler, K2, and TESS space missions. We find that the effective overshoot parameter, α ov,eff, increases from M = 1.0 M ⊙ to M = 1.2 M ⊙ before flattening out, although we note that the relationship between α ov,eff and mass will depend on the incorporated modeling choices of internal physics and nuclear reaction network. We also situate these results within existing studies of main-sequence convective core boundaries.


INTRODUCTION
Physical processes occurring within a star during its mainsequence dictate its future evolutionary history.Stars with masses above about 1.1 M ⊙ host well-mixed convective cores, but the structure of the regions near the outside of the convective core is currently not well understood.Physical processes such as convective overshooting, the process by which parcels of convective fluid pass the classical convective boundary because of their momentum, are thought to extend the well-mixed convective cores past the classical boundary, usually defined by the Schwartzchild or Ledoux criterion.Determining how best to model overshooting in 1-D stellar evolutionary codes is important, since incorporating core overshoot into stellar models increases the amount of hydrogen available to a main-sequence star, thereby increasing its main-sequence lifetime and altering the star's evolution, see Figure 1.A better understanding of this convective boundary mixing will better anchor our calibration of absolute ages in stellar modelling, in turn clarifying the ages of Corresponding author: Christopher Lindsay christopher.lindsay@yale.eduother astrophysical systems of interest, such as Milky Way progenitors (Chaplin et al. 2020) and exoplanet hosts (Huber et al. 2019).
In addition to changing a star's internal structure and evolutionary history, convective overshoot from the stellar core also changes the observable global properties of stars, such as their Radius and Effective Temperature.Previous studies used observed global properties of stars and have shown that convective overshooting is necessary in order for stellar models to reproduce observations of eclipsing binaries (e.g.Schroder et al. 1997;Pols et al. 1997;Ribas et al. 2000;Claret 2007;Claret & Torres 2018, 2019;Claret et al. 2021;Constantino & Baraffe 2018;Costa et al. 2019) and the color-magnitude diagrams of clusters (e.g.Maeder & Mermilliod 1981;Aparicio et al. 1990;Bertelli et al. 1992;Demarque et al. 1994;VandenBerg et al. 2006;Rosenfield et al. 2017).Extra mixing beyond the convective core is also known to emerge from numerical hydrodynamics, although internal gravity waves also contribute to the convective boundary mixing (e.g.Higl et al. 2021).See Anders & Pedersen (2023) for a recent review of observational and theoretical constraints on mixing processes at the convective boundaries in main-sequence stars.Different investigations find that the amounts of core overshoot needed to match a star's observable quantities depend on global stellar properties.For example, Claret & Torres (2016) studied a sample of eclipsing binary stars with observed masses, radii, temperatures, and elemental abundances and found that the size of the overshoot region has a positive dependence on the stellar mass.This result remains debated due to uncertainties in calibrating convective overshoot using eclipsing binaries (Constantino & Baraffe 2018;Claret & Torres 2019) and the exact relationship between overshoot and stellar properties such as mass, metallicity, and evolutionary state remains uncertain.Even so, published grids of stellar evolution models also frequently make use of core overshooting prescriptions, typically scaling the amount of overshooting with stellar mass (Demarque et al. 2004;Pietrinferni et al. 2004;Bressan et al. 2012).
Asteroseismology, or the study of stellar oscillations, has provided the means to directly study the interiors of stars.The global oscillation properties of solar-like oscillators, the frequency of maximum oscillation power(ν max ), and the large frequency separation (∆ν ), have been used extensively to determine global stellar parameters (Yu et al. 2018).However, they cannot be used to probe the deep stellar interiors we are interested in.Additionally, even modelling individual p-mode (pressure-mode) oscillation frequencies cannot give us detailed insights into the near-core structure of mainsequence stars with convective cores.In a previous paper (Lindsay et al. 2023), we have shown that pure p-modes are not able to sample the near-core regions directly, making in-ferences about the amount of core overshooting occurring in main-sequence stars difficult for some targets.Instead of analyzing the global asteroseismic properties (ν max or ∆ν) or the p-mode oscillations of main-sequence stars, we study dipolar oscillation modes which exhibit mixed character, which only arise after a star leave the main sequence.
After a star depletes its reserves of hydrogen in the core, its core begins to contract while its envelope begins to expand.At this stage of evolution, hydrogen is burned in a thin shell surrounding the now inert helium core.Since the convective motions in the core have also ceased by the time the star becomes a subgiant, the chemical composition gradient in these overshooting regions, having been frozen in at the main-sequence turnoff, serve as a fossil signature of the main sequence structure around the convective core.At the same time, the contraction of the stellar core is accompanied by an expansion of the envelope outside the hydrogen burning shell, leading to a large density contrast between the stellar core and envelope, allowing the study of these inner layers through mixed-mode asteroseismology (see Hekker & Christensen-Dalsgaard 2017, for a review of evolved star asteroseismology).Thus, studying the structure of subgiant and early first-ascent red giant branch stars can answer questions about main-sequence processes occurring above convective cores, since the main-sequence structural details of the star will be frozen into the star's inert core and the star's structure render their inspection observationally feasible through the asteroseismic analysis of mixed modes.This idea of using mixed modes observed in evolved stars to study interior processes occurring during the main sequence was first proposed in Deheuvels & Michel (2011).
The effects of main-sequence core overshoot on the structural properties of subgiants can be seen in Figure 2. Panel a of the Figure 2 shows the propagation diagrams of two main-sequence stellar models with the same mass and chemical compositions, but with different amounts of convective core overshoot.While significant difference can be seen in the deep, near-core layers, these are inaccessible to inspection using non-radial p-modes.After a main-sequence star evolves to a subgiant (panel b of Figure 2), a large density contrast develops between the stellar core and envelope.Some of the stars we study in this work are early red giant branch stars (see Figure 3) but, as shown in panel c of Figure 2, structural differences between stellar models with and without core overshooting remain for stars with log(g) values larger than 3.0, which is the case for all the stars in our sample.These differences don't remain forever, as after the stars evolve past the red giant branch luminosity bump, the structures of the stellar models with and without core overshooting are the same (panel d of Figure 2).
During the subgiant and red giant stages of evolution, ℓ = 1 oscillation modes with frequencies near ν max may propagate in two different regions: the envelope, which supports pressure modes (p-modes where the restoring force is pressure, with frequencies ν ≥ N and S ℓ=1 ) and the core, which supports gravity modes (g-modes where the restoring force is buoyancy, with frequencies ν ≤ N and S ℓ=1 ).Panels b, c, and d in Figure 2 show stellar model structures which may support ℓ = 1 mixed modes with frequencies near ν max , with the p-mode region above the N and S ℓ=1 curves, and the g-mode region below the N and S ℓ=1 curves.The observed mixed modes couple these two regions, with g-like character in the core and p-like character in the envelope (Scuflaire 1974;Aizenman et al. 1977).The observed stellar oscillations will depend strongly on the details of how the boundry between the core g-mode region and envelope p-mode region is configured, so we use individual mixed-mode oscillations to investigate the properties and structures around stellar cores in this work.
Space based photometry missions such as CoRoT (Baglin et al. 2006), Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2015) which observe many stars over long temporal baselines have made possible the observational detection of mixed modes in many stars possible.Since these mixed modes sample the interior layers of evolved stars near the core-envelope boundary, they can, and have been, used to constrain the amplitude of overshooting above convective cores.In particular, Deheuvels & Michel (2011) found from the CoRoT data of the solar-like oscillator HD49385 that the oscillation spectrum of the star can only be properly explained by an avoided crossing (a characteristic feature of onresonance mixed modes), whose shape constrains the amount of core overshooting above the stellar core during the main sequence to either very small or moderate values.Deheuvels et al. (2016) also made seismic estimates of the extent of convective cores in 8 low mass main sequence stars using data from Kepler.Viani & Basu (2020) examined 9 intermediate mass main-sequence stars from the Kepler LEGACY sample (Lund et al. 2017;Silva Aguirre et al. 2017) and found through asteroseismic modelling that the amplitude of convective overshoot from the main-sequence core increases with stellar mass.Noll et al. (2021) then studied the Kepler subgiant, KIC10273246, and found again that accounting for core overshooting improved their models' agreement with the observed oscillation mode frequencies.Building on these studies, we now make similar measurements from an analysis using a grid based modelling approach to determine the convective core boundary properties of a larger sample of subgiants and early first-ascent red giant branch stars observed by Kepler and TESS.
In this paper, we analyze a sample of 62 subgiants and early first-ascent red giant branch stars to determined how much convective core overshoot was present in these stars during their main sequence.Our objective in doing so is to determine relationships between global stellar properties and near core-mixing processes.Incorporating convective overshooting also changes the global properties of stars, such as their effective temperature, so we combine the asteroseismic data for our sample of stars with spectroscopic observables avaliable from the literature.The rest of this paper is organized as follows.We discuss our sample of subgiant stars in section 2 and describe our grid of models in section 3.In section 4 we explain how we compare the observed spectroscopic and asteroseismic properties of the subgiant target stars to our model grid.In section 5 we show how the mainsequence core properties of the stars in our sample depend on the global stellar parameters and place our results into the context of other studies focused on determining the relationship between core overshoot properties and stellar parameters.
2. THE SAMPLE We study a total of 62 stars in this work, 36 observed during the Kepler mission, 8 observed during the K2 mission, and 18 observed with TESS, of which 12 were observed in the TESS Southern Continuous Viewing Zone (CVZ) during its first year of operations.We selected our sample of stars to include subgiants and early first-ascent red giant branch stars whose frequency échelle diagrams display avoided crossings, indicating the presence of mixed-character modes whose frequencies sample the core/envelope boundary.Figure 3 shows a Kiel diagram (log(g) vs. Effective Temperature) of our sample.The different missions observed the stars for varying lengths of time, consequently, the quality of the power spectra used to determine the oscillation frequencies varies tremendously.

Kepler Stars
The nominal Kepler mission ran for just over 4 years, giving the best possible scenario for determining the asteroseismic mode frequencies we require for our fitting procedure.For this work, like in Ong et al. (2021a), we study a sample of stars observed with short cadence which was already examined using a grid-based modeling approach (Li et al. 2020a,b).The mode frequencies used in our analysis were measured in Li et al. (2020b).The global asteroseismic parameters, ∆ν and ν max , for these targets were derived in Serenelli et al. (2017), while the spectroscopic observables T eff and [Fe/H] were taken from Table 1 of Li et al. (2020a).When available, we also used stellar luminosity (L) measurements derived from the Gaia mission Gaia Collaboration et al. (2018).

K2 Stars
The time series photometry available from the K2 mission are only about 75 days long in each campaign (Howell et al. 2014), with 60 second sampling at short cadence.This results in significantly degraded frequency resolution, and, since the photometric noise was also higher due to decreased pointing stability, only 8 K2 subgiants and early first-ascent red giant branch stars studied in Ong et al. (2021a) showed significant oscillation power excess.The oscillation frequencies and spectroscopic properties for these stars are the same as those used in Ong et al. (2021a) (see section 3.3), although the mode frequencies for some of these stars were reanalyzed in González-Cuesta et al. (2023).

TESS Stars
During its nominal mission, stars observed by TESS can, in the worst case, be observed for only a single sector ( 2710 −1 10 0 m/M 10 0 N, S , (in units of ν max ) a. Main Sequence, both models at X H = 0.5   days) which is not much longer than the average oscillation mode lifetime.Therefore, in the frequency domain, the limited spectral resolution of the power spectrum is comparable to the mode line widths.In addition, the TESS pixels are much larger than Kepler's meaning the oscillation mode frequencies derived from TESS photometry are more sensitive to noise and suffer from more contamination when compared with Kepler-derived mode frequencies.Since the detection and extraction of oscillation mode frequencies is not the focus of this work, our TESS sample of stars is relatively small and made from studies already communicated through the asteroseismology community.The spectroscopic and asteroseismic observables of our sample of 6 stars observed by the nominal TESS mission (β Hyi, δ Eri, η Cep, ν Ind, TOI 197, and HD 38529) is further described in section 3.2 of Ong et al. (2021a).
In addition to these 6 stars, our sample includes 12 subgiants and early first-ascent red giant branch stars observed by TESS in its Southern Continuous Viewing Zone.The oscillation mode frequencies for these stars were fitted against the power spectra using 7 different peakbagging pipelines (J.M. J. Ong, in preparation).The spectroscopic properties (T eff and [Fe/H]) were taken from the 17th data release of the Sloan Digital Sky Survey (Abdurro'uf et al. 2022), while luminosity measurements were found with SED measurements in conjunction with GAIA parallaxes.

OUR GRID OF MODELS
There are two main classes of techniques to model the interiors of subgiants and early first-ascent red giant branch stars using asteroseismic and spectroscopic data.One that uses large scale grid searches to study many targets (as in Mc-Keever et al. 2019;Jørgensen et al. 2020;Li et al. 2020a;Nsamba et al. 2021;Ong et al. 2021b,a), and another, where the parameter space of the grid is adapted to each target individually (as in Ball et al. 2018Ball et al. , 2020;;Huber et al. 2019;Chaplin et al. 2020;Noll et al. 2021).These two approaches are often combined, with the results of coarser grid searches used to restrict the parameter space for further individual study.While boutique modelling involving optimization based parameter searches or individual dense grids created separately for each target is a good method for determining the properties of individual stars, it can be very slow and computationally expensive.Instead, we used a large scale grid search method for this study.
The grid we use in our fitting procedure is largely based on that used in Ong et al. (2021a).It consists of stellar models created with MESA version r12778 (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019)).These models use the solar chemical abundance mixture detailed in Grevesse & Sauval (1998), an Eddington gray atmospheric boundary condition, and the mixing-length prescription of Cox & Giuli (1968).Element diffusion was handled using the formulation of Thoul et al. (1994), and we included a mass-dependent scaling prefactor (see Viani et al. 2018).The parameters of the models in the grid, Sobol-sampled in almost an identical fashion to that used in Ong et al. (2021a), were the initial stellar mass (M ∈ [0.8M ⊙ , 2.0M ⊙ ), the mixing length parameter (α mlt ∈ [1.3, 2.2]), initial metallicity ([Fe/H] 0 ∈ [−1, 0.5]), and initial helium abundance (Y 0 ∈ [0.176, 0.32]).We construct models with sub-primordial (Y primordial ∼ 0.248) to avoid edgeeffects from our grid's construction affecting the fitting results, but we down-weight the likelihoods of stellar models with, Y 0 < Y primordial , as we explain in section 4.
Convective overshoot from the stellar core was incorporated according to MESA's implementation of overmixing with a step profile (cf.§2 of Lindsay et al. 2022).Unlike the grid used in Ong et al. (2021a), we Sobol-sample the step overmixing parameter, α ov , between 0 and 0.6, as shown in Figure 4.The result of MESA's implementation of step overmixing is that the well-mixed stellar core extends beyond the Schwarzschild boundary by a distance r ov , given by where H p is the pressure scale height at the convective boundary and r cz is the radius of the core convection zone.This ensures that when the convective core is very small, the overshooting region does not become unphysically large, as would be the case for overshooting by α ov H p .Another feature of the MESA implementation of step overmixing is that the overshoot region does not start exactly at the Schwarzschild boundary (estimated by where ∇ radiative = ∇ adiabatic ) but rather begins at a location f 0 H p into the convective core from the Schwarzschild boundary.This is because the mixing coefficient approaches 0 at the Schwarzschild boundary.In our grid, f 0 is set to 0.005.In order to avoid confusion when referencing the step overshooting parameter, and to make it easier to compare our results with studies that used other stellar evolutionary codes, we save the two different overshoot parameters in our grid.
1. Input Overshoot, α ov : This is the value entered into the MESA controls inlist when creating the stellar model tracks.We varied this parameter in our grid from 0 to 0.6.
2. Effective Overshoot, α ov, eff : The effective overshoot parameter is calculated directly from the MESAgenerated profile file for the time step at which the MESA-defined convective core mass was at its maximum.We define the effective overshoot parameter as α ov, eff = (r well-mixed − r cz )/H p where r well-mixed is the well-mixed core boundary, defined by looking for the central-most grid point where the gradient of the mean molecular weight gradient, For the subsequent analysis of the effective overshoot parameter, α ov, eff is set to 0 in our grid if the corresponding model tracks maintain a convective core for less than % of their main-sequence lifetime.This cutoff ensures that the overshoot parameters for models which do not maintain a convective core for a significant amount of time is 0. The 30% limit was chosen to ensure that models with varying initial compositions, but with masses ≲ 1.0M ⊙ , all maintain α ov, eff = 0 (see Figure 6).
The models we use in this work use MESA's 'basic.net'nuclear network, which considers lithium and beryllium to be at equilibrium and therefore is known to overestimate the modelled size of the convective during the main sequence, when compared with a stellar model calculated using a full nuclear network (see Noll & Deheuvels 2023).We tested how using a full nuclear network (MESA's 'h burn.net') would affect our calculation of the effective overshoot parameter (α ov, eff ) for a given model track and found that since α ov, eff is calculated using the difference between the well-mixed core boundary and the convection zone boundary, and both r well-mixed and r cz are decreased by the usage of a full nuclear network, our determination of α ov, eff for a given model track is only modestly altered by the choice of nuclear networks.The choice of nuclear network changes the overall size of the convective core, since with 'basic.net' the core abundance of lithium and beryllium is assumed to be at equilibrium, which in turn means the overall core energy production from the protonproton chain could be incorrect.Thus, the choice of network will change the oscillation frequencies of the stellar models, and may therefore alter the inferred stellar parameters obtained through asteroseismic modelling.
The grid input parameters which we vary, M, α mlt , [Fe/H] 0 , Y 0 , and α ov , are distributed uniformly using joint Sobol sequences of length 16382 over the ranges described previously.Since each of the 16382 tracks we calculate has a slightly different combination of input parameter values, we are able to obtain more precise estimates of our output parameters, when compared to if we used an equi-sampled grid.Each of these tracks were evolved using MESA from the premain sequence until the point where ∆ν = 9µHz (following Ong et al. 2021a).For each model in our grid, the oscillation mode frequencies are calculated using GYRE version 6.0 (Townsend & Teitler 2013).The radial and quadrupole (ℓ = 0 and ℓ = 2) p-mode frequencies are calculated within ±6∆ν of ν max .For the dipole (ℓ = 1 modes, we calculated both π-mode and γ-mode frequencies and mixed-mode coupling matrices according to the mode isolation construction of Ong & Basu (2020).π-modes refer to the pure p-modes whose frequencies only depend on the p-mode cavity of the stellar model, while γ-modes refer to the pure g-modes whose frequencies only depend on the g-mode cavity; mixed modes are those linear combinations of π and γ modes that are also eigenfunctions of the wave operator.Following Ong et al. (2021a), we compute γ-mode frequencies and matrix elements for γ modes from a lower-bound frequency of ν max − 7∆ν up to the n g = 1 γ-mode.The ℓ = 0 and ℓ = 2 modes, as well as the ℓ = 1 π− modes were corrected for inadequate modelling of the near surface layers (surface effects or surface term), before comparison with the observed mode frequencies, as discussed further in section 4. these include the spectroscopic parameters effective temperature (T eff ), metallicity ([Fe/H]), and when available, Luminosity (L).These variables are used to define the spectroscopic likelihood.The global asteroseismic observables, ∆ν and ν max , are also used just in the down-selection of the large grid.

Down-selecting the Complete Grid
The first step in our modelling procedure is to search for all models in the complete grid with properties close to the target star's observed values of T eff , [Fe/H], ∆ν, and ν max .We calculate a penalty function based on individual χ 2 values for each observed parameter, P, given by, and the total χ 2 is given by If luminosity is available for the target star, we include χ 2 L in the calculation of χ downselect .We down-select the complete grid to models satisfying χ 2 downselect < 10 2 .This greatly reduces computational cost by restricting the number of models requiring an expensive seismic likelihood evaluation from more than 14 million (the whole grid) to around ten thousand models per target star.

Spectroscopic Likelihoods
To quantify how well the models in the down-selected grid matches the spectroscopic observables of the target star, we calculate the following spectroscopic cost function, following Eq.( 2).For the 6 target stars (KIC6442183, KIC11137075, KIC11414712, δ Eri, EPIC212478598, and EPIC246305274) for which we could not find reliable luminosity measurements, we omitted χ 2 L from the calculation of χ 2 spec , obtaining We divide χ 2 spec by the number of spectroscopic parameters we use in order to ensure that the combination of spectroscopic constraints are weighted the same as the asteroseismic constraints, described in the next subsection.
The spectroscopic likelihood for every model in the downselected grid is finally calculated as (6)

Seismic Likelihoods
In order to calculate a seismic cost function, we must quantify how well a stellar model's oscillation mode frequencies match the observed frequencies of a given target.To do this, the model's set of oscillation modes must be compared to the set of observed modes, which entails matching each of the observed oscillation mode frequencies to a corresponding model mode frequency.The model modes also must be corrected for the surface effect, a frequency-dependent error in stellar model oscillation mode frequencies caused by our inability to model the near-surface layers of a star in one dimension.
We match the model oscillation mode frequencies to the observed frequencies in two different ways, depending on the angular degree (ℓ) of the mode.The radial and quadrupole modes (ℓ = 0 and 2) are matched based on their inferred values of radial order n p .We infer the n p values of the observed modes based on the observed mode frequency, the observed ∆ν as n p, obs = (ν obs /∆ ν ) − (ℓ/2).The n p values for a models' oscillation frequencies are returned from GYRE (Townsend & Teitler 2013) along with the mode frequencies.The matched ℓ = 0 and 2 modes are then used to determine the coefficients of the two-term surface term from Ball & Gizon ( 2014) by minimizing the quantity, where n denotes the radial order of the mode, ℓ denotes the angular degree of the mode, N is the total number of modes, ν obs, nℓ is the observed mode frequency, σ ν,obs, nℓ is that mode's associated frequency error, ν model, nℓ is the uncorrected model mode frequency, while δν surf, nℓ is the two term parametric correction of Ball & Gizon (2014).This correction is given by, where ν nℓ is the model mode frequency, I nℓ is the model mode inertia, ν 0 is set to ν max , and the coefficients (a −1 and a 3 ) are chosen to minimize the quantity in Eq. ( 7).
The above procedure gives us surface corrected model mode frequencies for ℓ = 0 and ℓ = 2 modes, and we can now define ν model, corr = ν model, uncorr + δν surf .For the ℓ = 1 modes though, following Ong et al. (2021c), we apply the surface correction to the dipole π modes, and recover surfacecorrected model mixed modes by coupling these surface corrected π modes to the γ modes, which remain unaffected by surface effects (Ong & Basu 2020).Surface-corrected dipole (ℓ = 1) mode frequencies are matched to the observed mode frequencies using an iterative nearest-neighbor search.
Next, for each combination of model in the down-selected grid and target, we calculate a seismic cost function with the form, seis is weighted by 1/N ν , the total number of matched oscillation modes in the set.We also weight χ 2 seis by W ν , a surface-term-weight, due to limitations involved with correcting stellar model frequencies for near-surface effects using only a power-law based frequency shift (Ball & Gizon 2014).The surface term makes the model mode frequencies larger than the observed frequencies, and this difference increases with frequency.The model of the surface term shown in Eq. ( 8) does not take this into account, and it is possible that after the surface term correction, a model that has frequencies lower than the observed frequencies would end up with a smaller χ 2 seis than a model with the expected behavior of the surface term.To account for this, we add the surfaceterm-weight, W ν , in Eq. ( 9) to give a larger weight to models whose uncorrected frequencies were larger than the observed frequencies.To do this, we assign W ν a value of 0.5 if all model mode frequencies are higher than their corresponding observed frequencies.If this condition is not satisfied, we set W ν to 1.
σ ν, eff in Eq. ( 9) accounts for the systematic error in the modeling due to grid undersampling (following Li et al. 2020a;Ong et al. 2021a).We first calculate the seismic cost function (Eq.( 9)) with σ ν, eff set to 0 to identify the best seismic fit model for each target.Using this best seismic fit model, we set σ ν, eff in Eq. ( 9) to the root-mean-squared dif-ference between the observed mode frequencies and the best fitting model's surface corrected mode frequencies.The values of σ ν, eff we find for each target in this work are reported in Table 3 and Table 4.This gives an indication about the asteroseismic goodness-of-fit that is reached for each star.
As an illustrative example for how we calculate σ ν, eff , the best mode match for the target KIC4346201 is shown in Figure 5.Note that the parameters of this best seismic fit model are not the best-fit parameters we report in our results, but are rather incorporated (along with all the parameters of all the other models in the grid) into the results through taking the likelihood weighted means of our parameters of interest, as described in the next subsection.The fact that the best seismic fit model's frequencies represent the best seismic match to the observed mode frequencies simply means that the seismic likelihood of this model is higher than the others.Additionally, we can see in Figure 5 that the most g-mode-dominated dipole frequency (left-most ℓ = 1 mode) does not agree with its corresponding model mode frequency, an issue which could be eliminated by using an optimization approach to fitting the modes (Noll et al. 2021;Noll & Deheuvels 2023), rather than the grid-based approach we take in this work.
For models with interiors that match observed stars, it is a known property of the surface term that the frequency differences between the observed and model p-modes should be largest at high frequencies, and lowest at low frequencies.We follow Basu & Kinnane (2018), Ong et al. (2021b), and the appendix B1 and B2 procedures of Cunha et al. (2021) in accounting for this by adding another penalty function to χ 2 seis with the form, This term is calculated using the 4 lowest frequency radial (ℓ = 0) modes and the 4 lowest frequency quadrupole (ℓ = 2) modes.The factor 1/4 comes from the number of modes that we are summing over, and the factor 1/10 gives a lower weight to χ 2 low n compared to χ 2 seis , calculated in Eq. ( 9) (following the procedure of Ong et al. (2021a)).We do this because we do not want to double count the 4 modes with the lowest frequencies, χ 2 low n is incorporated for the purpose of regularization, and is not meant to influence the overall posterior distribution.
For each of the 62 targets in our sample, we calculate the two terms (χ 2 seis and χ 2 low n ) of the seismic cost function.We then calculate a seismic likelihood (L seis ) for each model in the target star's associated down-selected grid as (11)

Estimating Stellar Parameters
For each of the stars in our sample, we determine total likelihoods, L tot , for each model in the down-selected grid by combining the spectroscopic and seismic normalized likelihoods as with Y p = 0.248 (Steigman 2010).The total likelihood of a given model incorporates this weight term in order to penalize models with initial helium abundance lower than the primordial value of helium, see Silva Aguirre et al. (2017).The total likelihoods are normalized by dividing each likelihood by the sum of all the models' total likelihoods, leaving us with an associated likelihood value, L tot for every model in the target star's down-selected grid.These likelihoods are used to infer our results, stellar parameters, especially the amount of convective core overshooting, for the 62 subgiant and early first-ascent red giant branch stars in our sample.

Best Fit Parameters for our Targets
To obtain the best fit parameters and parameter errors for each of the targets in our sample (section 2), we calculate the likelihood weighted means of each parameter saved in the grid discussed in section 3 10,000 times for different realizations of the spectroscopic parameters.The likelihood weighted mean for a given parameter, P, is calculated as, where N is the number of models in a target's downselected grid and P model is the model parameter.We draw 10,000 sample values of each spectroscopic parameter from a normal distribution with the mean equal to the literature-reported value of the parameter and standard deviation equal to the reported error, as was done in Ong et al. (2021a).Each sampling of the spectroscopic parameters results in a different likelihood weighted mean for our model output parameters, and hence, performing a Monte-Carlo over the spectroscopic parameters allows us to determine posterior distributions of our output parameters.The 16th, 50th, and 84th percentiles of the posterior distributions for stellar mass (M), radius (R), luminosity (L), temperature (T eff ), age, initial metallicity ([FeH] 0 ), initial helium abundance (Y 0 ), mixing length (α mlt ), and effective overshoot pa- rameter (α ov, eff ) are listed in tables Table 1 and Table 2 for the stars observed by Kepler and TESS respectively.Previous modelling studies of main sequence stars and subgiants (Lebreton & Goupil 2014;Li et al. 2020a;Noll et al. 2021) have shown that a high level of degeneracy exists between the initial helium abundnace and mass, due to an apparent anti-correlation between mass and Y 0 .Analysis of the posterior distributions for our inferred stellar parameters, however, does not show this strong degeneracy between mass and Y 0 .This is likely due to our inclusion of luminosity in our likelihood function, which is not included in Lebreton & Goupil (2014) and Li et al. (2020a).

The relationship between α ov, eff and Mass
Figure 6 shows the results for the effective overshoot parameter (α ov, eff ) as a function of the stellar mass.The errors bars for both parameters show the 16th and 84th percentiles of the output parameter posterior distributions.We can see that for sufficiently low mass stars (M ≤ 1.0M ⊙ ) the modelled α ov, eff is 0. This is because the stellar models at this mass do not maintain a convective core for a sufficiently long amount of time (> 30% of their main sequence lifetimes.)At slightly higher masses, from M ≈ 1.0M ⊙ to M ≈ 1.15M ⊙ , the value of the effective overshoot parameter increases with stellar mass from α ov, eff ≲ 0.05 to α ov, eff ≈ 0.1 -0.15.At masses higher than 1.2M ⊙ , our results for α ov, eff are seemingly independent of stellar mass.The points in Figure 6 are colored by our initial metallicity ([Fe/H] 0 ) results.We do not see a strong relationship between our results for α ov, eff and [Fe/H] 0 .
To place our results shown in Figure 6 into context, we compare them with other studies of convective core overshooting in Figure 7 and Figure 8.We note that in Figure 7 and Figure 8, we plot our values of α ov, eff along side the input overshoot values the modelling work in Deheuvels et al. and Noll & Deheuvels (2023).Our construction of α ov, eff (see section 3) is meant to correspond to the input overshooting values in the other stellar evolution codes.We have defined α ov, eff (which in general is time-dependent) such that the thickness of the overshoot region is equal to α ov, eff times H p no matter what α ov is.We have however restricted our attention to its value at when the convective core is at its maximum size.Since the overshoot region thickness is the physical quantity ultimately being modified in the stellar model by overshooting, it is α ov, eff rather than input α ov which is constrained by the mode frequencies.
We employ this construction of α ov, eff since some other studies (Claret & Torres 2018, 2019) calculate the size of the overshoot region as α ov times the pressure scale height at the convective boundary (H p ), instead of calculating the overshoot region size as α ov multiplied by the minimum of H P and the radius of the convective core (as is done in MESA, see section 3 and Deheuvels et al. (2016)).Additionally MESA's overshooting regions begin slightly interior to the convective boundary, unlike in the other stellar modelling codes used in the other works, such as YREC (Demarque et al. 2008).
Overall, the results from this work agrees with the previous overshooting versus mass results from Deheuvels et al. (2016), Noll et al. (2021), andNoll &Deheuvels (2023), as well as with the results at the low mass end of the sample from Viani & Basu (2020).The spread in our reported values of α ov, eff for stars in the mass range of 1.1M ⊙ ≤ M ≤ 1.4 is also similar to the spread in the overshoot amplitude results of Deheuvels et al. (2016), for the same range of masses.
On the higher mass end of their sample (M ≥ 1.3M ⊙ ), Viani & Basu (2020) reported significantly larger values of effective overshoot parameter, when compared with our results for stars in our sample with comparable masses.We identify two possible reasons for this discrepancy.First, the models calculated in Viani & Basu (2020) had set the temperature gradient, ∇ to ∇ ad in the convective overshooting regions above the models' convective cores.This is analogous to the "Step Penetrative Overshoot" we describe in section 2 of Lindsay et al. (2022) or the "convective penetration" described in Anders & Pedersen (2023).As shown in figure 8 of Viani & Basu (2020), the effective overshoot parameter results are larger, for a given stellar mass, if the stellar models used in the analysis maintained an adiabatic temperature gradient in the overshooting region.This temperature gradient difference alone cannot fully explain the large discrepancy between our results and the results of Viani & Basu (2020).The second major difference between the two modelling methods is that our modelling procedure sampled α ov, input far more densely (see Figure 4) when compared to the procedure of Viani & Basu (2020), which only allowed about 8 different α ov, input values per target star.The denser sampling of both α ov, input and initial mass in our work results in more precise results of stellar mass and overshoot parameter, whereas the sparse sampling of α ov, input in Viani & Basu (2020).means that one model with a high overshooting parameter could dominate the likelihood function and significantly increase the output α ov, eff result.It should also be noted that the Viani & Basu (2020) models did not include the gravitational settling of heavy elements, and this too could increase their estimate of overshoot, however, that is unlikely to explain the large difference.2019), and we show how our results compare to theirs in Figure 8.We note that the masses of the stars in their sample are higher than most of the stars in our sample.The reported exponential overshoot parameters, f ov from Mombarg et al. (2021) are multiplied by 10 before plotting in Figure 8, as Mombarg et al. (2021) uses the conversion factor of 10 in their analysis.We also include the stellar mass and overshooting amplitude results of Claret & Torres (2018) and Claret & Torres (2019) in Figure 8.The stars studied in these works were higher-mass main-sequence stars in double line eclipsing binary systems, with M ≥ 1.2M ⊙ .Claret & Torres (2018) and Claret & Torres (2019) modelled the stars' masses and core overshoot parameters by comparing the observed masses, radii, and effective temperatures of the stars to stellar models made with varying amounts of convective core overshoot.In order to convert the exponential overshoot parameters, f ov , reported in Claret & Torres (2019)  Comparing our results (red points in Figure 8) to the results of Claret & Torres (2018), Claret & Torres (2019), and Mombarg et al. (2021) shows significant differences in stellar mass/overshooting amplitude relationship.The combined results of Claret & Torres (2018) and Claret & Torres (2019) show the overshoot amplitude increasing with stellar mass from α ov, eff ≈ 0 at M = 1.2M ⊙ to α ov, eff ≳ 0.15 at M = 2.0M ⊙ .On the other hand, for the same mass range, Mombarg et al. (2021) finds that the step overshooting parameters for their sample of γ Doradus stars are larger in magnitude (α ov, eff > 0.1) and less dependent on stellar mass.
Although our sample of stars do not go past masses of M ≈ 1.5M ⊙ , we report higher values of α ov, eff for stars in the mass range of 1.2 to 1.5 M ⊙ range when compared to the eclipsing binary work of Claret & Torres (2018), Claret & Torres (2019).The majority of the γ Doradus stars studied by Mombarg et al. (2021) appear to have higher values of α ov, eff compared to the results from our sample at similar masses.This is not an apples to apples comparison though, since we are using detailed asteroseismic data of individual oscillation modes while Claret & Torres (2018) and Claret & Torres (2019) produced their overshooting results based on the masses, radii, and effective temperatures of binary stars, without asteroseismology, while Mombarg et al. (2021) trained a machine learning model to return stellar masses and overshoot parameters based on observed period spacings between adjacent modes.In addition, the grid of models used in this work sampled the input parameters more densely, and covered a wider range of parameters when compared to the grids of models used in the eclipsing binary and γ Doradus work.Finally, the models calculated by Claret & Torres (2018), Claret &Torres (2019), andMombarg et al. (2021) implemented penetrative overshoot which, as discussed previously, would cause higher output values of α ov, eff , when compared with an analysis done with models with a radiative temperature gradient in the overshooting region.
We note that our reported relationship between overshooting amplitude and stellar mass also depends on stellar modelling choices beyond the input parameters described in sec-tion 3, such as the choice of nuclear network.For example, Noll & Deheuvels (2023) find that changing the choice of nuclear network could change the size of the convective core on the main sequence.Since, like overshooting amplitude, the convective core size is mass dependent, the choice of nuclear network used in the model grid could alter the resultant relationship determined by fitting asteroseismic observations to the model mode frequencies calculated from the models of that grid.Quantifying the magnitude of this dependence on nuclear network will be the aim of subsequent work.

SUMMARY AND CONCLUSION
We have conducted a study analyzing 62 subgiant and early first-ascent red giant branch stars with high quality asteroseismic data from Kepler and TESS.The goal of this work was to better understand convective overshooting above main sequence star convective cores.Previous studies have analyzed main sequence convective core properties using asteroseismology (Deheuvels et al. 2016;Noll et al. 2021;Noll & Deheuvels 2023), but the total number of targets that are bright enough and in the mass range of interest is quite small.Additionally, Lindsay et al. (2023) showed that core regions of stars are poorly constrained by non-radial p-mode oscillations on the main sequence.Therefore, we analyze subgiants and early first-ascent red giant branch stars, which are much brighter, and whose dipolar mixed-mode oscillation frequencies sample the interior structures of subgiants.The main sequence structure of stars with convective cores alter the stellar structure, and these effects persist after the main sequence, when the star evolves into a subgiant (see Figure 2).
Our sample consisted of 44 stars observed during the Kepler and K2 missions, as well as 18 stars observed with TESS.The spectroscopic parameters, as well the individual oscillation mode frequencies, were taken from various different literature sources (see section 2).To determine the amount of convective core overshooting, as well as the other properties of the stars in our sample, we implement a grid-based modelling technique, as opposed to performing a boutique modelling of every star in our sample.The grid we use densely Sobol-samples stellar mass, mixing length, initial helium abundance, metallicity, and the input overshoot parameter, α ov, input .We match the observed mode frequencies of each star in our sample to the surface corrected model mode frequencies for the models in our grid.Using the matched mode frequencies, we calculate a seis-mic likelihood (Eq.( 11)) for each combination of star and model.These seismic likelihoods are combined with a spectroscopic likelihood (Eq.( 6)) to produce a combined 'total' likelihood (Eq.( 12)), which we use to determine the likelihood weighted mean parameter values for all the different quantities in our model grid.The posterior distributions of our output quantities were determined by taking 10,000 samples of the spectroscopic observables and recalculating the likelihood weighted mean parameter values for each random sample.
These modelling results of the stars observed by Kepler are detailed in Table 1 while Table 2 contains our results for the stars observed by TESS.Our results for stellar mass and effective overshoot parameter, α ov, eff , are visualized in Figure 6 which shows that α ov, eff increases with stellar mass from M = 1.0M ⊙ to M = 1.2M ⊙ .For stars with masses greater than 1.2 M ⊙ our results indicate a weak, or no correlation between stellar mass and overshooting amplitude.We compare these findings to previous work in section 5 while discussing how differences between our results and previous works can be explained by differences in modelling procedure.

Figure 1 .
Figure 1.Two evolutionary tracks showing the effects of incorporating core overshooting into a stellar model.Both tracks show M = 1.2M ⊙ , Z = Z ⊙ models evolved using MESA version r12778 from the pre-main sequence to the red giant branch.The track with overshooting parameter α ov = 0.4, shown in red, shows that incorporating overshoot changes the evolution leading up the subgiant branch.

Figure 2 .
Figure2.Propagation diagrams showing the interior structural evolution of two 1.3M ⊙ stellar models with and without overshooting.The main-sequence structure of the models (at a central hydrogen fraction of 0.5) are shown in panel a.Both models' future subgiant structures (at T eff = 5500 K) are show in panel b.The models' structure on during their first-ascent red giant branch stage (log(g) = 3.0) are shown in panel c, where there is still a significant difference between the stellar models with and without overshooting.The models' structure after the red giant branch bump (log(g) = 2.25) is shown in panel d.In all panels, the Brunt-Väisälä (buoyancy) frequency (N) is indicated with the solid lines, while the ℓ = 1 Lamb frequency (S ℓ=1 ) is indicated with the dotted lines.The horizontal dot-dashed lines show the frequency of maximum oscillation power, ν max .The models incorporating convective core overshoot with overshoot parameter α ov = 0.3 are shown in red, while the models without overshoot are shown in black.Both models have the same initial metallicity ([Fe/H] 0 = −0.025),initial helium abundance (Y 0 = 0.284), and mixing length parameter α MLT = 1.75.

Figure 3 .
Figure 3.Our sample of subgiants and early first-ascent red giant branch stars studied in this work shown on a Kiel (log(g) vs. T eff ) diagram.The different symbols refer to which mission observed the star.The background gray curves show solar-abundance evolutionary tracks with masses between 1.0 and 1.8 M ⊙ in steps of 0.2M ⊙ .The tracks are for modes without core overshoot.

Figure 4 .
Figure 4. Comparison of sampling of the mass-overshoot plane in the model grid used in this work (gray points -Sobol sampling in a finite range) against that of the grid used in Ong et al. (2021a) (orange points -random sampling over the mass-overshoot relation of Viani & Basu 2020, shown with the blue points).

Figure 5 .
Figure 5. Échelle diagram showing the best mode match for the subgiant target KIC4346201 used to determine σ ν, eff .The best seismic model is determined by minimizing Eq. (9) for all models in the down-selected grid (with σ ν, eff set to 0).The observed modes are shown with open circles, and the surface-corrected modes of the model are shown with dots.Each observed mode is connected to its corresponding model mode with a line.The colors of the points and lines indicate the angular degree of the modes (ℓ = 0 in blue, ℓ = 1 in orange, and ℓ = 2 in gray.)

Figure 6 .
Figure 6.Effective overshoot parameter α ov, eff versus stellar mass for our sample of 62 subgiants and early first-ascent red giant branch stars.The points are colored by the modelling results for initial metallicity.Each point's position and color represents the 50th percentiles of the given parameter for each target.The error bars are given by the 16th and 84th percentiles of each parameter's posterior distribution.

Figure 7 .
Figure 7. α ov, eff versus stellar mass results (red points) in comparison with other asteroseismic studies of convective core overshoot.

Figure 8 .
Figure 8.Our α ov, eff versus stellar mass results (red points) in comparison with other studies of convective core overshoot, including those shown in Figure 7 as well as the eclipsing binary work of Claret & Torres (2018) and Claret & Torres (2019), shown in black triangles.The machine learning results for γ Doradus stars from Mombarg et al. (2021) are shown in maroon triangles.
Mombarg et al. (2019) andMombarg et al. (2021) studied convective core overshoot properties in intermediate mass stars through asteroseismology of γ Doradus stars.Mombarg et al. (2021) used machine learning techniques to study core overshooting in the same sample of stars from Mombarg et al. ( to step overshooting parameters, α ov , we follow Claret & Torres (2019) and multiply the reported values of f ov by 11.36 (11.36 f ov = α ov ) before plotting those points in Figure 8.When a star is reported in both Claret & Torres (2018) and Claret & Torres (2019), we plot the value from Claret & Torres (2019).

Table 1 .
Modelling results for the stars in our sample observed by the Kepler and K2 missions.

Table 3 .
σ ν, eff values for all of the Kepler targets in our sample.A higher value of σ ν, eff indicates that our model grid has a higher degree of effective undersampling for that target.

Table 4 .
σ ν, eff values for all of the TESS targets in our sample.A higher value of σ ν, eff indicates that our model grid has a higher degree of effective undersampling for that target.