Asteroseismic Inversions for Internal Sound Speed Profiles of Main-sequence Stars with Radiative Cores

The theoretical oscillation frequencies of even the best asteroseismic models of solar-like oscillators show significant differences from observed oscillation frequencies. Structure inversions seek to use these frequency differences to infer the underlying differences in stellar structure. While used extensively to study the Sun, structure inversion results for other stars have so far been limited. Applying sound speed inversions to more stars allows us to probe stellar theory over a larger range of conditions, as well as look for overall patterns that may hint at deficits in our current understanding. To that end, we present structure inversion results for 12 main-sequence solar-type stars with masses between 1 and 1.15 M ⊙. Our inversions are able to infer differences in the isothermal sound speed in the innermost 30% by radius of our target stars. In half of our target stars, the structure of our best-fit model fully agrees with the observations. In the remainder, the inversions reveal significant differences between the sound speed profile of the star and that of the model. We find five stars where the sound speed in the core of our stellar models is too low and one star showing the opposite behavior. For the two stars in which our inversions reveal the most significant differences, we examine whether changing the microphysics of our models improves them and find that changes to nuclear reaction rates or core opacities can reduce, but do not fully resolve, the differences.


INTRODUCTION
The combination of high-precision photometric time series data from Kepler (Borucki et al. 2010), astrometric parallax data from Gaia (Gaia Collaboration et al. 2016), and high-resolution spectroscopic measurements of effective temperatures and metallicities (for example from the Kepler Follow-up Program, Furlan et al. 2018) provides an opportunity to test stellar evolution theory at unprecedented precision.
However, for stars with the highest quality asteroseismic data, there are still discrepancies between models and observations.This tension between theoretical and observed frequencies suggests that our models need to be improved, although it does not directly suggest what those improvements should be.
We aim to gain insight into the potential underlying structural differences between stellar models and observations using the technique of asteroseismic structure inversions.This technique uses the differences between the frequencies of an observed star and its model to infer localized information about the structure differences (see e.g.Gough & Thompson 1991;Gough 1993;Pijpers 2006;Bellinger et al. 2020a;Buldgen et al. 2022a).
In the case of the Sun, structure inversions have been used to study the equation of state, diffusion of heavier elements, and nuclear reaction rates in connection to the solar neutrino problem, (for a review see, for example, Basu 2016;Christensen-Dalsgaard 2021).The high precision and large number of modes observed for the Sun allow structure inversions to probe a large extent of the solar interior, from 0.06 R ⊙ to 0.96 R ⊙ .This, however, is not the case for other stars.Current asteroseismic observations are typically limited to modes of spherical degree l = 0, 1, 2, with a few l = 3 modes being observed in the best target stars.This limits the range that can be probed with local structure inversions to the nearcore region, fractional radii between ∼ 0.05 and 0.35 (Bellinger et al. 2020a).
Nevertheless, there are several examples of structure inversions performed on stars other than the Sun, including the solar analogues 16 Cyg A and 16 Cyg B (Bellinger et al. 2017;Buldgen et al. 2022b), a mainsequence star with a convective core (Bellinger et al. 2019a), and a few subgiants with mixed modes (Kosovichev & Kitiashvili 2020;Bellinger et al. 2021).Inversion techniques are also being developed for more massive stars (Vanlaer et al. 2023) and more evolved stars (Giammichele et al. 2018).By looking at a larger number of stars, we can test the theory of stellar structure and evolution under a broader range of conditions, such as different masses, metallicities, ages, and evolutionary stages.Examining several stars at once also provides the opportunity to look for overall trends that may hint at deficits in our current understanding of stars.In this work, we focus on studying the most solar-like starsmain-sequence stars with radiative cores-using structure inversions.

FORWARD MODELING
The goal of a structure inversion is to infer the differences between the actual stellar structure and that of a reference model.As the structure inversion equation is based on a linear perturbation approach, the reference model must be suitably close to the actual star.Hence, we typically use the best-fit model obtained with some modeling procedure called forward modeling.Here, we describe the forward modeling procedure used to obtain our reference model for each target star.We created a grid of 16384 tracks using the r22.05.01 version of the stellar evolution code MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2022).We vary the initial mass, initial helium mass fraction, metallicity, and mixing-length parameter using a Sobol sequence (see Appendix B of Bellinger et al. 2016;Sobol' 1967).Table 1 gives the range that was covered in each parameter.All models in this grid use metal abundances scaled to the GS98 solar composition (Grevesse & Sauval 1998), and the corresponding opacity tables from OPAL (Iglesias & Rogers 1993, 1996) in the high-temperature range, and Ferguson et al. (2005) in the low-temperature range.The equation of state data are calculated with the MESA default blend of OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS (Irwin 2004), and Skye (Jermyn et al. 2021).For details of how this blending is handled, see Jermyn et al. (2022).We use the pp_cno_extras_o18_ne22.net reaction network and take our reaction rates from JINA REACLIB (Cyburt et al. 2010) and NACRE (Angulo et al. 1999), with additional tabulated weak reaction rates (Fuller et al. 1985;Oda et al. 1994;Langanke & Martínez-Pinedo 2000).Electron screening is included via the prescription of Chugunov et al. (2007).Thermal neutrino loss rates are from Itoh et al. (1996).Convection in the models is described using the time-dependent local convection formalism of Kuhfuss (1986), which in the limit of long time steps reduces to standard mixing length theory as described in Cox & Giuli (1968).The implementation details are given in Jermyn et al. (2022).We account for atomic diffusion through gravitational settling, as described in Paxton et al. (2011).We use an Eddington-gray atmosphere and include the structure of the atmosphere out to an optical depth of τ = 10 −3 when calculating both our oscillation frequencies and structure kernels.The adiabatic frequencies of the models were computed using GYRE (Townsend & Teitler 2013;Townsend et al. 2018).
For each target star, we find reference models by fitting the observed frequencies, effective temperature, and metallicity of each star.We take our frequency data from the Kepler LEGACY sample (Lund et al. 2017) and the Kepler ages (KAGES) sample (Davies et al. 2016).In the case of 16 Cyg A and 16 Cyg B we use the frequencies given in Roxburgh ( 2017) labeled as Roxburgh (Davies).Spectroscopic measurements of the effective temperature and metallicity are from the combined stellar parameters reported by the Kepler Follow-Up program (Furlan et al. 2018, their Table 9).These values are computed by combining the results of four different spectroscopic analysis pipelines.We also adopt their suggested uncertainties of 100 K, 0.1 dex for T eff and [Fe/H] respectively.The observational parameters we consider for each star are listed in Appendix A. For each target star, we search our grid to find the model that minimizes where and Here N is the number of frequencies and the subscripts 'obs' and 'mod' refer to the observations and the model respectively.The model frequencies used to calculate χ 2 ν are first corrected for surface effects using the two-term correction from Ball & Gizon (2014).We scan our grid to find the parameters (M, Y initial , Z initial , α mlt , X c ) that minimize χ 2 fit .In the process, we interpolate in central hydrogen abundance along each track, but not we do not interpolate between the tracks.In order to reduce the computational time necessary to find a best-fit model, we consider for subsequent analysis only models that are within 6σ of the effective temperature and metallicity values, as well as within 10σ of the FLAME luminosity value from Gaia DR3 (Gaia Collaboration et al. 2016, 2022;Creevey et al. 2023).We then use these parameters, given in Appendix A, to compute the reference model that will be used for structure inversions of each star.We have made the FGONG structure files of our reference models as well as the inlists used to generate them publicly available.1

INVERSIONS
With a suitable reference model, we aim to use the differences between the frequencies of an observed star and the frequencies of the reference model to infer the underlying structure differences.We do this through the use of stellar structure kernels, which express the sensitivity of an oscillation mode frequency to a small perturbation to the structure.Mathematically, this is expressed in the kernel equation (Dziembowski et al. 1990): (5) Here i is the index of the mode which corresponds to a specific pair of radial order (n) and spherical degree (l), δν i /ν i = (ν i,obs −ν i,mod )/ν i,mod is the relative difference in frequency between the observed mode (ν i,obs ) and the corresponding mode of the reference model (ν i,mod ), f 1 and f 2 are the stellar structure variables being considered, and K i are the kernel functions of each mode.The mode kernel functions (K i ) are obtained though a linear perturbation of the oscillation equation (for more details, see Gough & Thompson 1991, Kosovichev 1999, or Thompson & Christensen-Dalsgaard 2002).Initially, mode kernels were derived in terms of the squared sound speed c 2 and density ρ (Dziembowski et al. 1990).From this expression, mode kernels for other pairs of variables have been derived, including for density and the first adiabatic exponent Γ 1 (e.g., Gough & Thompson 1991;Gough 1993), isothermal sound speed u = c 2 /Γ 1 and helium mass fraction Y (e.g., Basu & Christensen-Dalsgaard 1997;Kosovichev 1999;Buldgen et al. 2015Buldgen et al. , 2017)), and convective stability parameter A and Γ 1 (e.g., Elliott 1996;Kosovichev 1999;Buldgen et al. 2017).For more details on changing the structure variable pair, see Kosovichev (2011); Buldgen et al. (2017).Equation 5can also include a term that corrects for the surface effect; however, we instead correct for this in the calculation of the frequency differences.For the remainder of this paper, when we discuss model frequencies, we refer to the frequencies that have been corrected for surface effects.
Each oscillation mode is sensitive to many points within the star, so to obtain localized information we implement the method of optimally localized averages (OLA, Backus & Gilbert 1968, 1970) which uses a linear combination of the frequency differences.Neglecting second-order effects, Equation 5 becomes Here K is the averaging kernel and C is the cross-term kernel.These are constructed using a set of inversion coefficients c i : If K is normalized to 1 and C is small, then Equation 8 reduces to and K can be interpreted as the weight function of a mean over the structure difference δf 1 /f 1 .This is why K is called the averaging kernel.In other words, if the coefficients c i are chosen such that K has most of its amplitude around a single target radius, r 0 , then the same linear combination of frequency differences provides a localized average difference of the structure variable f 1 at that target radius.

Localized averaging kernels
To construct a localized averaging kernel, we use the method of multiplicative optimally localized averages (MOLA).For a MOLA inversion, we define a weight function, J = (r − r 0 ) 2 , that penalizes any amplitude of the averaging kernel away from the target radius.In addition to the target radius, there are two trade-off parameters that must be chosen: the error suppression parameter, µ, and the cross-term suppression parameter, β.The inversion coefficients are found by minimizing where E ij are the elements of the error-covariance matrix.Strictly speaking, there is another trade-off parameter in this formulation, the normalization of J, which we have set to 1, while other implementations often use 12.However, this only changes the relative weight of the first term in Equation 9 and can be counteracted by changing µ or β.Thus, while the optimal values of µ and β vary with the normalization of J, the inversion results do not.
The other standard method of constructing an averaging kernel is a variant of MOLA known as the method of subtractive optimally localized averages (SOLA, Pijpers & Thompson 1992, 1994).Previous works studying 16 Cyg A and 16 Cyg B have used SOLA (Bellinger et al. 2017;Buldgen et al. 2022b).We choose to use MOLA because it requires setting only two free parameters, as opposed to the three required for a SOLA inversion.Additionally, we find that MOLA is better able to suppress the amplitude of the averaging kernel at the surface.For details on the differences between MOLA and SOLA see Basu & Chaplin (2017, Chapter 10).
The next important consideration for a structure inversion is which pair of structure variables to use.We use the (u, Y ) pair because the Y kernels have low amplitude everywhere except in helium ionization zones (Basu 2003), which naturally suppresses the cross-term kernel at the radii we are targeting.The trade-off to this approach is that using Y as a structure variable requires the assumption of an equation of state.In the solar case, the error introduced by this assumption is significant in comparison to the other sources of uncertainty (Basu & Christensen-Dalsgaard 1997); however, due to the larger uncertainties on asteroseismic frequencies this is not the case for stars other than the Sun.Using Y as a structure variable requires calculating several partial derivatives of Γ 1 .To be consistent with the blend of equation of state tables used in MESA, we obtained these directly from MESA's equation of state module.

Trade-off parameters
As Equation 9shows, there are two trade-off parameters that must be chosen in the course of a structure inversion.The parameter µ controls the balance between a well-localized averaging kernel and the amplification of uncertainties.To choose an appropriate value of µ for each inversion, we utilize a set of calibration proxy models.These models are found using the process described in Section 2; however, they have slightly higher χ 2 fit values than the reference model.Since the structure of these models is known exactly, they can be used to determine how well the inversion recovers the underlying differences.We provide the details of this process in Appendix B.1.Before accepting our inversion results, we visually inspect the averaging and cross-term kernels of all target radii for each star to ensure that the averag-ing kernels are well localized and the overall amplitude of the cross-term kernels are low.

Stellar mass & radius
Another complication for inversions of stars other than the Sun is the lack of precise measurements of the stellar mass and radius.Since the frequencies of a star scale with its mean density, a mismatch in the mean density will lead to an offset in the inversion results (Basu 2003).
To minimize this, we invert for the relative difference in û = uR/GM , where R and M are the star's radius and mass, respectively, and G is the gravitational constant.This is done by using mode kernels computed in a dimensionless form and using the dimensionless frequency differences.Previously, Bellinger et al. (2021) used dimensionless frequency differences calculated by subtracting off the weighted mean of the frequency differences (for details, see Basu 2003).
We have found that this method results in correct dimensionless frequencies only when the frequency differences caused by an incorrect mean density are larger than the differences introduced by the structural variation.Whether this is true cannot be determined purely by comparing the observed and modeled frequencies.Instead, we use a new method of calculating the dimensionless frequency differences using the dependence of the large frequency separation (∆ν) on the mean density.The large frequency separation is the mean frequency difference between successive radial modes and is a proxy for the root mean density of the star ( Vandakurov 1967).The details of this method can be found in Appendix B.2.
We calculate our value of ∆ν by taking the slope of a linear fit to the relationship between the l = 0 modes and their respective radial orders and use this to calculate the dimensionless frequency differences.These corrections mean that the uncertainty of our frequency differences are no longer independent, and hence, the error co-variance matrix E, used in Equation 9 is not diagonal.We calculate it using a Monte Carlo simulation where each frequency is perturbed, 10 000 times with Gaussian noise according to their measured uncertainties.These perturbations are applied before the frequencies are corrected for the surface effect, and so this procedure also accounts for the error correlation introduced by the surface term correction.This same set of perturbed frequencies is then used to calculate the final inversion results.We take the average of this distribution as our final inversion result and report the standard deviation as the uncertainty.
This method of uncertainty estimation occurs after both the reference model and inversion parameters have been selected, and only propagates uncertainties due to the underlying frequencies.It is the same as the traditional definition of inversion uncertainties (e.g., Equation 4 of Bellinger et al. 2019a) except that it accounts for correlation introduced during the pre-processing of our frequency differences.
To validate both our method of finding reference models and our inversion results, we also obtain a reference model and inversion results using solar data that have been degraded to the level that was expected of results from Kepler.These results are given in Appendix B.3.

Overall inversion significance
For each star, we attempt structure inversions at six target radii, r 0 /R = 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, although in some cases we are only able to find suitable averaging kernels at five target radii.To quantify the disagreement between each target star and its model, across all target radii, we calculate a χ 2 inv as follows.For each star, there is a set of inversion results and their associated uncertainties v j ± u j .Since all the target radii use the same underlying data, their errors are correlated.The correlation between two target radii, r j , r k , is (Basu & Chaplin 2017).
where c i (r j ) is the inversion coefficient of the i-th mode for the j-th target radius and σ i is the relative uncertainty of the ith mode frequency.The error correlation matrix, E, is the matrix with components E rj ,r k between all different target radii.The covariance matrix then is, where U is a diagonal matrix with the uncertainty of the inversion result for each target radius.Then where V is the vector of inversion results at each target radius.This χ 2 inv summarizes the overall significance of the inversion results for each star across all target radii, with larger values indicating larger disagreement.
In summary, after finding a reference model, we calculate the surface-term-corrected dimensionless frequency differences between the target star and the model.We then use our set of calibration models to choose µ at each target radius and obtain our set of averaging kernels.With this, we use the inversion coefficients and the frequency differences to obtain our inferred localized differences in û between the observed star and our reference model.

RESULTS AND DISCUSSION
Together, the Kepler LEGACY and KAGES samples provide oscillation data for 95 stars.Since we are specifically searching for close matches to stars with radiative cores, we apply two criteria to our reference models: that they have a radiative core throughout their mainsequence evolution and that they have a χ 2 fit < 20.We obtain suitable reference models for 34 stars.Of these, 12 have enough frequencies observed (approximately 35) to form well-localized averaging kernels.
Figure 1 shows the inversion results for each of these 12 stars as a function of the target radius.We define our relative differences such that a positive inversion result indicates a sound speed that is higher in the star than in the model.We provide more information about the reference model and averaging kernels of each star in Appendix C.
The inversion results of the twelve stars in our sample can be broken down into three groups: (a) those where the û of the best-fit model is in agreement with the observations, (b) those where the û of the model is too high (resulting in an inversion result below zero), and (c) those where the û is too low (resulting in an inversion result above zero).Taking into account the uncertainties of the inversion results, we identify six stars in group a, five stars in group b, and one star in group c.Thus, half of our sample show significant differences, which suggests that there are limitations in the physics of our reference models and that these limitations most often result in internal values of û that are too low.Now we seek to understand why some of our stellar models show good agreement in û while others show significant disagreement.We search for correlations between χ 2 inv and properties of the reference model, as well as the surface rotation rate (P rot ) and magnetic activity indicator (S ph ) values for each star, as given by Santos et al. (2018).We calculate Spearman's rank correlation coefficient, ρ s , which describes the strength of the monotonic, but not necessarily linear, correlation between two variables.These results are shown in Figure 2. We use bootstrapping to obtain estimates of the uncertainty of these coefficients.
The strongest correlation is with χ 2 ratios (ρ s = 0.81).This χ 2 is a measure of how well our reference model matches the observed frequency separation ratios (r 10 , r 02 ) of the observed star.These ratios are known to be insensitive to the surface effect, (Roxburgh & Vorontsov 2003) and thus χ 2 ratios serves as a different metric for how well the internal structure of a star is reproduced in the model.The strong correlation between χ 2 ratios and χ 2 inv reaffirms that the differences found from structure inversions are internal structure differences rather than problems with the near-surface layers.
We find significant positive correlations of the discrepancies between the star and stellar model with the central abundance of 12 C (ρ s = 0.72) and 14 N (ρ s = 0.66) of the model, as well as the amount of energy generated by the CNO cycle (ρ s = 0.61).A similar-strength correlation in the opposite direction is found with the central hydrogen abundance (ρ s = −0.62).That all of these properties have a similar strength of correlation is unsurprising, as they are mutually correlated.When the central hydrogen value is lower, reactions other than the pp-chain can happen more easily.Primarily, this is an increase in energy generated by the CNO cycle.At the same time, a very small amount of energy is generated by the triple alpha process.This is not significant compared to the total energy generation of the star, but it does increase the equilibrium abundance of 12 C. Additionally, the CNO-II pathway converts 16 O into 14 N which increases the equilibrium abundance of 14 N.In general, we see a moderate correlation between the significance of the û differences inferred by inversions and more evolved main sequence stars.

Individual Stars
We now discuss the results of both our forward modeling and inversion procedures for a few individual stars.

16 Cyg A and 16 Cyg B
First, we focus on the solar analogs 16 Cyg A and 16 Cyg B. In Figure 3 we compare the frequencies (before and after surface-term corrections) and frequency separation ratios of our reference models with the observations.As these are two of the most well-studied main sequence stars in the Kepler field, they have already been studied using structure inversions by Bellinger et al. (2017); Buldgen et al. (2022b).In our results, as well as the two previous studies, there is excellent agreement between the models and the observations.In the case of 16 Cyg B, the models used in all three works are within 1σ agreement with observations.For the case of 16 Cyg A, our inversions show differences that are less than 1.5σ, which is similar to the values obtained by Bellinger et al. (2017) and Buldgen et al. (2022b).Bellinger et al. (2017) report their inferred u values, as well as their inferred values of the stellar mass and radius, which allows us to compare û values directly, as shown in Figure 4.All the points for 16 Cyg B are in good agreement.For 16 Cyg A, there is slight disagreement at a target radius of 0.25, but it is not significant.Despite the use of different reference models, a different implementation of OLA, and different inversion param- We now turn to the two stars in our sample that show the largest differences with respect to our models: KIC 6603624 and KIC 6116048.We show the frequencies and frequency separation ratios in Figure 5.Both of these stars have points where our inversions infer internal sound speed differences greater than 10 percent, and in contrast to other stars in the sample, these large differences are significant compared to their uncertainties, so first we verify that our inversions are able to recover differences of this magnitude.The δ û/û between our reference model for KIC 6116048 and our reference model for KIC 6603624 reach ∼15% in the region probed by structure inversions, and so we test our averaging kernels by attempting to recover the difference between the two models.We do this twice, once with KIC 6603624 as the reference model and then again using KIC 6116048 as the reference model.The result of these inversions are shown in Figure 6.Both sets of averaging kernels infer the correct shape of the true δ û/û curve.The averaging kernels of KIC 6603624 infer the correct value of δ û/û within the uncertainties at every target radius.This is not the case for the averaging kernels for KIC 6116048, where two points differ from the correct value by ∼ 2σ.Nevertheless, we conclude that our inversion procedure is able to recover differences around 15%.

Exploring the effects of microphysics
We now explore several changes to the microphysics in our models in an attempt to reduce the sound speed dif-ferences inferred by our inversions.A full investigation of the microphysics across all twelve of the stars studied here is beyond the scope of this work, and hence we focus on KIC 6603624 and KIC 6116048, the two stars discussed in section 4.1.2.For each star, we create three new models using the same mass, initial composition, and mixing-length parameter as our original reference model, although we allow these new models to have a different central hydrogen abundance.For the first model, motivated by the correlation to CNO energy production, we multiply the rate of the 14 N + p → 15 O + γ reaction by a factor of 0.1.For the second model, we multiply the ppII/ppIII rate 3 He + 4 He → 7 Be + γ by a factor of 0.25.For the last model, we modify the opacity by a factor of 0.85 in the parts of the model with log T > 6.7.For each of these three new tracks, we select a new reference model using the fitting procedure discussed in Section 2. Figures 7 and 8 show the results of these changes for KIC 6603624 and KIC 6116048, respectively.For each change in the microphysics we show the û profiles of each new reference model as well as the result of structure inversions.The changes to the core opacity and the 14 N + p → 15 O + γ reaction rate both increase û from the original reference model, with the opacity change resulting in a larger difference both to the central û and the frequency differences computed with respect to the observations.The change caused by modifying the 3 He + 4 He → 7 Be + γ reaction rate results in a smaller change to the û that is only apparent inside r/R < 0.07.As expected, when we apply changes that increase the internal û, the û difference inferred by the structure inversion decreases.We find better agreement in the û profile even when the fit of the model is worse than our original reference model.These changes improve our models at the deepest target radii, but have little effect at the larger radii probed by inversions.

CONCLUSIONS
Here, we have used asteroseismology to infer the detailed core structure of the best solar-type stars observed by the Kepler mission.We focused on main sequencestars with radiative cores and expanded the number of such stars studied with structure inversions from 2 to 12.After obtaining our reference models from a grid created using MESA, we use a set of calibration models to obtain our inversion parameters.We then use these inversion parameters to infer the relative difference in dimensionless squared isothermal sound speed between our reference model and the target star.In our sample, we identify three groups: those where the û of our reference model agrees with the observed star (group a, 6 stars), those where the û of our model is higher than that of the star (group b, 1 star), and those where the û of our model is lower than that of the star (group c, 5 stars).We also find significant correlations in our results, suggesting that our models of older main-sequence stars with more energy being generated by the CNO cycle have larger differences between model and star.To explore how changing the microphysics affects our inversion results, we tested the effects of changing nuclear reaction rates and core opacities, for the two stars with the most significant differences.These changes to the microphysics reduced the discrepancy between model and star at the innermost target radii.
In future work, we aim to extend our analysis to an even broader set of stars, including main-sequence stars with convective cores and more evolved stars with mixed-mode oscillations.Main-sequence stars with convective cores are particularly interesting since their dominant source of energy is the CNO cycle.Thus, they are a natural next step to explore the correlations found in this work between the inferred sound speed differences and CNO energy production.
We thank the anonymous referee for their constructive feedback that improved the manuscript considerably.The research leading to the presented results has received funding from the ERC Consolidator Grant Dipo-larSound (grant agreement #101000296).SB acknowledges NSF grant AST-2205026.SB would also like to thank the Heidelberg Institute of Theoretical Studies for their hospitality during the early days of this project.This paper includes data collected by the Kepler mission.Funding for the Kepler mission is provided by the NASA Science Mission Directorate.In addition, this work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.We have also used the gaia-kepler.funcrossmatch database created by Megan Bedell.

REFERENCE MODEL PARAMETERS
Table 2 provides the non-seismic constraints that were used to find the reference models of the 12 main sequence stars with radiative cores discussed in this work.Table 3 provides the model parameters for each reference model.

B. INVERSION DETAILS
Here we provide details on how we chose our inversion parameters and how we calculate the non-dimensional frequency differences used in our inversions.Additionally, we present the results of applying our modeling and inversion methods to degraded solar data.

B.1. Inversion Parameter Selection
For each target radius we find the value of µ that minimizes: where the angle brackets denote a mean across the set of calibration models, (δ û/û) inv is the sound speed difference inferred by the inversion, (δ û/û) True is the true sound speed difference, and the uncertainty of the inversion result is σ 2 inv = i c 2 i σ 2 i , where σ i is the relative uncertainty of the ith mode.This is not the uncertainty reported in our final results, as it does not account for uncertainty correlation introduced by our surface term and mean density corrections (see Section 3.3).In general, the term in the square brackets dominates, as with increasing values of µ the quality of the averaging kernel degrades faster than the uncertainty of the final result is reduced.We minimize M separately for each target radius, and so the value of µ can vary between different target radii of the same star.This optimization is more stable when only one variable is minimized, and so we set the cross-term trade-off parameter β = 0.
As the effect of this choice is similar for all target radii across all stars in our sample, we use the innermost target radius of KIC 6603624 as an example.Figure 9 shows the averaging kernels, cross-term kernels, and model-model inversion results for several values of β using the same value of µ.Increasing β has the expected effect of damping the cross-term kernel, however it also reduces the quality of the averaging kernel, making it less localized.This is particularly noticeable when β = 10 000.As the model-model inversions show, the results are much more sensitive to the quality of the averaging kernel than to the amplitude of the cross-term kernel, and so we conclude that setting β = 0 is justified.

B.2. Mean Density Scaling
To mitigate the effect of a difference in mean density between a star and its model, we calculate the dimensionless frequency differences before applying our structure inversions.One method of obtaining this difference was proposed in Basu (2003) and used by Bellinger et al. (2021).This approach notes that the proportional scaling with mean density shows up as a constant offset in the frequency differences.This constant offset can be approximated by taking a weighted mean of the frequency differences.This term can then be subtracted from the raw frequency differences to remove any differences due to mean density.We have found that this approximation is valid only when the frequency differences due to different mean densities are larger than the differences resulting from structure differences.Thus, in this work, we take a different approach and use the large frequency separation of the star and its reference model to calculate a dimensionless frequency difference.The dimensionless frequency is where R is the stellar radius, M is the stellar mass, and G is the gravitational constant.For two stars with stellar radii R 1 , R 2 and stellar masses M 1 , M 2 the dimensionless relative frequency difference of any given mode is with this Equation B3 reduces to .Inversion results for KIC 6116048 found using two different methods of calculating the dimensionless frequency differences.The blue points show the results when frequencies are scaled using the large frequency separations, described in Section 3, and the orange points use the differences calculated from a weighted mean, described in Basu (2003).
While this method results in different values for the dimensionless frequency differences, the effect on the inversion result is small compared to the uncertainty of the inversion result, as seen for KIC 6116048 in Figure 10.This small difference can be understood by carrying forward the effect of a small difference in mean density though the inversion procedure.
A difference of mean density shows up as a constant offset when calculating the dimensionless frequency differences.Mathematically, this is expressed as where δq/q is the offset introduced by a difference in mean density.As δq/q is independent of the frequencies, its contribution to the final inversion result will be Thus, the error introduced by a mismatch in mean density is proportional to the sum of the inversion coefficients for each target radii.We do not explicitly try to minimize this sum; however, the uncertainty of each result, which we do attempt to minimize, depends on the magnitude of the inversion coefficients.Thus, in the process of a standard inversion, we reduce the effect of a difference in mean density.Equation B7 also suggests a check to determine if a difference in mean density is the dominant difference present in our inversion results.While the sum of the coefficients will be different for each target radius, δq/q will be the same.Thus, if the mean density differences are dominating the inversion results, a plot of the inversion results at each target radii divided by the sum of the coefficients for that target radii should be a straight line.We checked this for all the stars in our sample and did not find such a constant, and so we conclude that the error introduced by a difference in the mean density is not the dominant source of difference in our inversion results.

B.3. Sun as a star
In addition to the twelve target stars, we also obtain a reference model and structure inversion results using solar data that have been degraded to the level that was expected of results from Kepler (for details, see Lund et al. 2017).Table 4 lists the parameters of the reference model obtained with these data.The parameters of our model are comparable to those found across all the pipelines used in Lund et al. (2017)  the quality of the fit in û, we compare the û profile of our reference model to that of the calibrated standard solar model S of Christensen-Dalsgaard et al. (1996), shown in Figure 11.Although we are not able to reproduce the full structure of a model calibrated with all the solar data, our reference model is a close match in the area probed by structure inversions using only lowdegree modes.These differences are of the same order of magnitude as the differences inferred between model S and the Sun (Basu et al. 2009).Thus, despite the limitations of the degraded data, we find a reference model sufficiently close for structure inversions.Using this reference model, we obtain suitable averaging kernels at four target radii and infer the difference in û using the degraded solar frequencies as shown by .At all four target radii, our structure inversions show agreement within 1σ in û.Helioseismic inversions that use non-degraded solar data do show differences between the structure of the Sun and the structure of calibrated solar models (e.g., Basu 2016); however, this results from using many more modes, with higher precision and at higher angular degrees, than are available for stars observed by Kepler.

C. FULL INVERSION RESULTS
For each target star, we attempt a structure inversion at target radii of r 0 /R = 0.05, 0.10, 0.15, 0.20, 0.25, 0.30.The target radius, however, is not necessarily the fractional radius where the averaging kernel is at its maximum value.We report the location of the maximum and FWHM of each averaging kernel and the û values inferred at each target radius of each star in Table 5. Figures 13 and 14 show the averaging and cross-term kernels for these stars, respectively.We show in Fig- ure 15 the results of model-model inversions between our reference model and one of the calibration models, as a test of the averaging kernel's ability to recover a known difference.

Figure 1 .
Figure1.Comparisons of the internal structure of stars as revealed by asteroseismology and the structures of best-fitting stellar evolution models.Relative differences are given in terms of the dimensionless squared isothermal sound speed û and span the near-core region of 0.05-0.3away from the stellar center point.The points indicate the inferred value of δ û/û between the star and the reference model at the target radius.The vertical error bars indicate the uncertainty of each inversion result from the propagation of the uncertainty of the observed frequencies.The horizontal error bars represent the full width at half maximum of the averaging kernel.The dashed horizontal line indicates complete agreement between the model and observations; points above this line imply that û of the star is larger than that of the model.The color bar indicates the statistical significance of the inferred difference, with lighter colors showing more significant results.The letter after the star's identifier indicates which group the star is in, as described in the text.The values given in the lower left of each plot indicate the mass (M ), initial helium mass fraction (Yinit), initial metallicity (Zinit), and central hydrogen mass fraction (Xc) of each reference model.We also report the overall significance of the inversion results, χ 2 inv .

Figure 3 .
Figure3.Modeling results for 16 Cyg A (top) and 16 Cyg B (bottom).In each case, the left plot shows the Frequency Échelle diagram comparing the frequencies of the reference model to the observations before applying any correction to account for surface effects.The center panel compares the reference model frequencies after applying the two-term surface correction fromBall & Gizon (2014).The color and shape indicate the spherical degree l: 0 (blue squares), 1 (black triangles), 2 (orange diamonds), and 3 (red circles).The right plot shows the frequency separation ratios r10 (pink) and r02 (light blue).In all plots, the open points represent the values from the reference model and the filled points represent the observed values.

Figure 4 .
Figure 4. Inversion results for 16 Cyg A (left) and 16 Cyg B (right).The blue points show the inversion results from this work.The orange points are the results from Bellinger et al. (2017).Since they report u we use their reported values of M and R to calculate û.

Figure 6 .Figure 7 .
Figure 6.Model-model inversions to recover the û difference between the model for KIC 6603624 and KIC 6116048.The left figure shows the result of using KIC 6603624 as the reference model, and the right figure shows the result of using KIC 6116048 as the reference model.In both plots, the black line represents the true value of δ û/û and the colored points show the result of the inversion.Different target radii are shown in different colors and correspond to the color of the averaging and cross term kernels shown in Figures 13 and 14.

Figure 9 .
Figure 9. Results of varying β for the innermost target radius of KIC 6603624.The left (center) panel shows the averaging kernels (cross-term kernels) that result from the indicated value of β.The right panel shows the result of a representative model-model inversion for each value of β.The true value of δ û/û is indicated by the vertical line.The error bars show the uncertainty on the inversion result.

Figure 12 .
Figure 12.Inversion results of the degraded solar data.All symbols and colors have the same meaning as in Figure 1.

Figure 15 .
Figure 15.Results of using the averaging and cross term kernels shown in Figures 13 and 14 to recover the difference between the reference model and one of the additional models used to calibrate our inversion parameters.The black line indicates the true δ û/û between the two models.The color of each point matches its corresponding averaging and cross term kernels in Figures 13 and 14.

Table 1 .
Grid Parameters Spearman rank correlation between the maximum significance inversion result of each star and various properties of the reference model.The color correlates with the absolute value of the correlation coefficient, which is a measure of the strength of the correlation.The value of each correlation coefficient are provided on the right side of the figure.We estimate the uncertainty of each correlation coefficient using bootstrapping.
Figure 8. Results of modifying the physics used to evolve each model of KIC 6116048.All colors and symbols have the same meaning as in Figure 7.

Table 4 .
Reference Model using Degraded Solar Data . To assessFigure 11.Relative difference in û between our reference model, calculated from degraded solar data, and the calibrated solar model S Christensen-Dalsgaard et al. (1996).The shaded region shows the area that can be probed by structure inversions using only the reduced mode set of the degraded frequency data.