Stars located below the instability strip in the Hertzsprung-Russell diagram have external convective zones and are usually called solar-like stars. Turbulent motions inside these convective layers excite gravito-acoustic modes (Goldreich and Keeley, 1977) that are used in asteroseismology to characterize the internal structure and dynamics of the stars. With the development of continuous (from months to years) high-precision photometry from space with several missions such as CoRoT (Baglin et al., 2006), Kepler (Borucki et al., 2010), K2 (Howell et al., 2014), and TESS (Ricker et al., 2014) hundreds of main-sequence stars and dozens of thousands red giants have been observed.

For many of these solar-like stars, the fundamental parameters have been estimated using the seismic measurements. Using the so-called “global asteroseismic scaling relations” (e.g., Brown et al., 1991; Kjeldsen and Bedding, 1995; Chaplin et al., 2011), masses and radii are obtained with a typical precision of 10 and 5%, respectively (e.g., Kallinger et al., 2010). These estimates are model-independent because they rely on three observables: the effective temperature of the star, T_eff, the frequency of the maximum power where the modes are located, ν_max, and the large frequency spacing, Δν. The observed data from the Sun and other stars are used as references. Better results can be obtained using stellar models to derive mass, radius, and also age. In the case of grid-based modeling, when fitting spectroscopic observables and global seismic parameters all together to a pre-calculated grid of models, we can reach a precision better than 3% in radius, 6% in mass, and 23% in age (Serenelli et al., 2017). The improvement is even better when incorporating the information from the individual mode frequencies with systematic uncertainties of 1% in radius, 3% in mass, and 15% in age (e.g., Mathur et al., 2012; Lebreton et al., 2014a,b; Metcalfe et al., 2014; Silva Aguirre et al., 2015; Creevey et al., 2017).

The precision reached by the asteroseismic inferences when the information of the individual modes is used, is among the best for field stars. Therefore, they are particularly interesting for other research fields such as galactic-archaeology studies (e.g., Miglio et al., 2013, 2015; Anders et al., 2016) as well as to properly characterize the exoplanet radius through the transit method (e.g., Huber et al., 2013; Van Eylen et al., 2018a,b).

The combined action of convection and differential rotation in main-sequence solar-like stars can produce magnetic fields of dynamo origin (e.g., Brun and Browning, 2017). Under the action of stellar winds, stars of masses below ~ 1.4M_⊙ slow down during their evolution (e.g., Skumanich, 1972; Kawaler, 1988; Matt et al., 2012; García et al., 2014; Réville et al., 2015). Sometimes, these magnetic fields change regularly, producing magnetic cycles similar to the 11-year cycle of the Sun in which the poles change polarity, requiring 22 years to complete a full magnetic cycle. The perturbations induced by the magnetic fields modify the properties of the acoustic modes measured in the Sun (e.g., Woodard and Noyes, 1985; Pallé et al., 1989; Broomhall et al., 2014) and in other stars (e.g., García et al., 2010; Salabert et al., 2016; Kiefer et al., 2017; Karoff et al., 2018; Santos et al., 2018). In particular, mode frequencies are shifted toward higher frequencies as the magnetic cycle evolves toward the maximum activity.

The shifts in the mode frequencies can introduce systematic errors in the inferred global stellar parameters, specially if the star is only observed during the maximum of an on-going magnetic activity cycle. However, if the frequency shifts induced by the magnetic field smoothly increase with the frequency, as it is in fact the dominant change in the Sun, this effect will be partially or perhaps mostly masked by the filtering of the so-called surface effects (see Howe et al., 2017 where the time variation in the surface term due to the solar activity cycle was analyzed in detail).

Such filtering is required at present in any model fitting procedure to deal with the unknown physics of the upper layers of the star (e.g., Kjeldsen et al., 2008; Ball and Gizon, 2017). This is perhaps one of the reasons why these potential systematic uncertainties have not drawn great attention to the asteroseismic community.

However, there are additional terms in the frequency shifts induced by the magnetic activity, already present in the Sun and that seem to be more important in other stars. As it is going to be described below, the magnitude of the shifts in some stars not only are much larger than in the Sun but the shifts of the different low-degree modes are clearly different. This can introduce a bias in the stellar parameter determination. Moreover, as uncovered by Salabert et al. (2018), in some stars the variation of the frequency shifts with frequency does not follow a simple smooth function but has a more sinusoidal shape, suggesting that the magnetic perturbation is occurring in deeper layers compared to the Sun. This frequency changes are not filtered with the surface correction terms.

In this paper we study the impact of both effects in the extraction of the global stellar parameters providing additional uncertainties that takes into account these biases induced by the magnetic field on the frequencies. The observed and simulated data used are described in section 2 whereas in section 3 we explain in detail the frequency shifts considered in this work. In section 4 we describe the methodology used to perform the model fitting. The results are discussed in section 5 and some conclusions are presented in section 6.

We have selected two main-sequence stars studied by Salabert et al. (2018) and that show frequency changes due to magnetic activity. The first one is KIC8006161 that has about 1M_⊙ and a radiative core. The second one is KIC9139163 which is a more massive star, of about 1.4M_⊙, and presumably with a convective core. The former star seems an ideal choice for a solar analog because it has a large range of detected modes, including some ℓ = 3, and it shows frequency changes correlated with the photometric magnetic activity that depends on angular degree ℓ (Santos et al., 2018) and frequency ν (Salabert et al., 2018). The frequency range of the acoustic modes detected in KIC9139163 is smaller, but it corresponds to the typical number of modes we have for an F8 main-sequence star without mixed modes. This star is still perfect for our analysis as it shows both angular degree and frequency dependent shifts (see Santos et al., 2018; Salabert et al., 2018, respectively). The spectroscopic parameters and frequency sets used for both stars are given in Table 1. For KIC8006161 the full set of frequencies reported by Appourchaux et al. (2012) includes ℓ = 0 modes with radial orders in the range n = 15−30, ℓ = 1 modes with n = 13−29, ℓ = 2 modes with n = 15−29 and ℓ = 3 modes with n = 19−24. For KIC9139163 the range of frequencies considered was: ℓ = 0, n = 11−29, ℓ = 1, n = 11−29, and ℓ = 2, n = 10−28. The values of the frequency of maximum power ν_max and the large separation Δν were not used in the fitting procedure but are indicated in Table 1 for guidance.

The effect or bias that magnetic activity can introduce on the determination of the fundamental stellar parameters will be influenced by the observational constraints (spectroscopic and asteroseismic inputs) and the stellar evolution and oscillation codes considered (underlying physics, free parameters considered, etc.). Since observations and theoretical models can be improved in the near future we have found convenient to partially remove such effects by considering in addition to the real stars, two models as proxy stars. They are constructed from the best-fit model for each star (obtained as indicated in section 4). We used the output of the best-fit model as input parameter for computing the frequencies of the corresponding proxy star. In the minimization procedure the same set of (n, ℓ) modes were considered. Although these proxies will be free of systematic errors other than the bias introduced by the frequency shifts from the simulated magnetic activity (to be discussed in the next section), assuming these targets free of observed errors would be unrealistic. So we have adopted the criteria described below.

Considering the theoretical frequencies for the proxies would assume that near-future models will introduce systematic errors much lower than the observed ones. Presumably this would not be the case due to the complex physics required to model the upper layers both in the stellar structure models and in the oscillation laws for sonic, turbulent, non-adiabatic motions. Although some progress has already been done for improving the physics of the uppermost layers (e.g., Trampedach et al., 2017), we can only expect that in the coming years surface effects would introduce diminished discrepancies between theoretical and observational frequencies, but remaining higher than the observational errors. Hence it is convenient that our proxy stars include some departure from the models which gives rise to similar surface effects. In particular we computed the eigenfrequencies of the proxy stars with a different surface boundary condition, namely δP = 0 rather than the standard match to the analytical solution for an isothermal atmosphere at the top of the model. Blue points in Figure 1 correspond to the differences between the proxy stars (the same models with δP = 0 as boundary condition) and the model frequencies (computed with the standard boundary conditions). Upper panel is for KIC8006161 and bottom panel for KIC9139163. As expected these differences show a typical smooth frequency function with δω → 0 as ω → 0 as any other surface effect. In the same figure we have compared such effect with the surface terms found for the actual stars after the standard model fit procedure detailed in the section 4 was carried out. The red points correspond to the differences between the observed frequencies and the ones from the best models. Note that in our simulation the surface errors are of opposite sign but we have not found any indication of this being relevant in the final results.

Regarding the frequency errors in the proxies we have considered two cases: one case with the observed errors and the other one with half of the errors. Similarly, theoretical spectroscopic parameters were used but with errors half the observed ones for each star; these inputs are indicated in Table 1. The error in the frequencies are given by a combination of the frequency resolution, the signal-to-noise ratio of the spectrum, and the width of the modes (see e.g., Libbrecht, 1992; Toutain and Appourchaux, 1994). Maintaining the width and the SNR constant the error will be reduced as T where T is the length of the observations. Thus, to reduce the error by a factor of 2, it is necessary to measure four times longer.

Model fitting is based on a grid of stellar models evolved from the pre-main sequence to the TAMS using the MESA code (Paxton et al., 2011, 2013, 2015), version 10398. The OPAL opacities (Iglesias and Rogers, 1996), the GS98 metallicity mixture (Grevesse and Sauval, 1998) and the exponential prescription for the overshooting given by Herwig (2000) were used, otherwise the standard input physics from MESA was applied. Eigenfrequencies were computed in the adiabatic approximation using the ADIPLS code (Christensen-Dalsgaard, 2008).

As specifically implemented in MESA, the overshooting formulation corresponds to that introduced by Herwig (2000) and it is a simplified version of the formulation given by Freytag et al. (1996). Here the overshooting parameter f_ov is defined by the equation

where D_ov is the diffusion coefficient, D₀ the diffusion coefficient at the transition layer between the standard convection zone and its overshooting extension, z is the distance from this transition layer, H_v the velocity scale height of the overshooting convective elements at z = 0, and H_p the pressure scale height at the same point. To establish the transition layer, first the Schwarzschild criterium is used to determine the edge of the convection zone and then the z = 0 point is shifted a radius 0.01H_p (or f_ovH_p/2 if f_ov < 0.02) into the convection zone where D₀ is computed. From z = 0 and toward the radiative zone, the overshooting formulation is used.

For convective cores lower than 0.001 stellar masses overshooting is never added. The same f_ov value was used for the core and the envelope and throughout the evolutionary sequence.

Different grids of models were used for each star. For KIC8006161 (and the related proxy) the grid is composed of evolution sequences with masses M from 0.94M_⊙ to 1.06M_⊙ with a step of ΔM = 0.005M_⊙, initial abundances [M/H] from −0.12 to 0.40 with a step of 0.05 and mixing length parameters α from 1.5 to 2.2 with a step of Δα = 0.1. Since these stars are expected to have convective core only at the beginning of the main sequence, and adding overshooting can cause longer lived convective cores, no overshooting was considered in this case.

The initial metallicity Z and helium abundance Y were derived from [M/H], constrained by taking a Galactic chemical evolution model with ΔY/ΔZ = (Y_⊙−Y₀)/Z_⊙ fixed. Assuming a primordial helium abundance of Y₀ = 0.249 and initial solar values of Y_⊙ = 0.2744 and Z_⊙ = 0.0191 (consistent with the opacities and GS98 abundances considered above) a value of ΔY/ΔZ = 1.33 is obtained. A surface solar metallicity of (Z/X)_⊙ = 0.0229 was used to derive values of Z and Y from the [M/H] interval.

For a typical evolutionary sequence in the initial grid, we save about 100 models from the ZAMS to the TAMS. Owing to the very rapid change in the dynamical time scale of the models, tdyn=(R3/GM)1/2, such grid is too coarse in the time steps. Nevertheless, the dimensionless frequencies of the p-modes change so slowly that interpolations between models introduce errors much lower than the observational ones. This procedure was discussed in more detail in Pérez Hernández et al. (2016) and was found safe and consumes relatively less time.

For KIC9139163 the procedure is similar except that the grid is composed of masses from 1.32M_⊙ to 1.60M_⊙ with a step of ΔM = 0.01M_⊙ and initial abundances [M/H] from −0.10 to 0.40 with a step of 0.05. In this case overshooting was included, the parameter fov goes from 0 to 0.04 with a step of 0.01.

A χ² minimization, including p-mode frequencies and spectroscopic data, was applied to the grid of models. The general procedure is similar to that described in Pérez Hernández et al. (2016), but for these stars all the eigenfrequencies are approximately in the p-mode asymptotic regime and hence a simplified procedure can be done. Specifically we minimize the function

Regarding the spectroscopic parameters, we have included the effective temperature T_eff the surface gravity logg, the surface metallicity Z/X, and the luminosity L. Thus,

where δT_eff, δ(Z/X), δg, and δL correspond to differences between the observations and the models whereas σ_Teff, σ_ZX, σ_g, and σ_L are their respective observational errors. Their values are given in Table 1.

The term χ²_freq in Equation (3) corresponds to the frequency differences between the models and the observations after removing a smooth function of frequency in order to filter out surface effects not considered in the modeling. This surface term is computed only using radial oscillations, and after scaling the frequency differences with the dimensionless energy I_nl, namely,

where B_i are constant coefficients, P_i a Legendre polynomial of order i, and x corresponds to ω linearly scaled to the interval [−1, 1]. A value of k = 2 was adopted for the same reason than in Pérez Hernández et al. (2016).

Whence the surface term is determined, we consider radial as well as non-randial modes for computing the corresponding minimization function

where j = 1, …N runs for all the modes and σ_ωj is the relative error in the frequency ω_j.

The term Ij-1S(ωj) subtracted to the relative frequency differences in Equation (6) includes a constant coefficient which contains information on the mean density of the star. In fact, introducing the dynamical time t_dyn = (R³/GM)^1/2 one has δω_n0/ω_n0 = δt_dyn/t_dyn+δσ_n0/σ_n0 where σ_nl = ω_nl/t_dyn are the dimensionless frequencies. As in Pérez Hernández et al. (2016) we restore it in the minimization procedure adding a new χ² term, but increasing the error compared to the formal one derived from individual frequencies since this term can also include some influence from surface effects as discussed in Pérez Hernández et al. (2016). Specifically we introduce the quantity

where B₀₀ is an offset fixed such that ν^model−ν^obs would be positive within errors for the lower radial modes. The parameter σ_B0 is the error associated to B₀ and was taken as σ_B0 = B₀₀. Typically we found values σB0~10-3, at least one order of magnitude higher than the formal error found in a typical fit to Equation (5).

To estimate the uncertainty in the output parameters we assumed normally distributed uncertainties for the observed frequencies, for σ_B0, and for the spectroscopic parameters. We then search for the model with the minimum χ² in every realization and compute mean values and standard deviations.

Quantifying the influence of the frequency shifts caused by the magnetic activity on the determination of the stellar parameters with asteroseismic techniques cannot be expressed with a simple mathematical rule because most of the techniques rely on some kind of model fitting that includes different constraints on a multiple parameter problem. For instance, quite different results could be obtained if one deals with the overshooting by introducing a free parameter or if it is fixed a priori. The same can be true when considering the helium and metallicity abundances, they can be implemented as two independent parameters or they can be coupled through a ΔY/ΔZ enrichment law. Moreover, if rather than a direct comparison between observed and theoretical frequencies, derived quantities are used to determine some specific parameters (e.g., see Verma et al., 2019 for using the acoustic glitches to determine the helium abundance) then the sensibility to the magnetic perturbations can also be different. In this work we have not tried to consider a broad range of situations but rather limit the analysis to some few examples, which we believe can be representative of some model fitting procedures. Although generalizing our quantitative results to other cases should not be correct, we hope that our results can at least serve as order-of-magnitude guidance for a large range of cases.

If observations are limited to periods of months, magnetic activity can cause mode frequencies to depart from their raw values by tenths of μHz or even more. We found that an ℓ-dependent shift due to magnetic activity similar to that reported in Santos et al. (2018) can introduce errors in the age of the order of 10% at most, and of a few percents in the mass and radius. These figures can be higher than formal uncertainties and in particular we found that the uncertainties were more relevant for the 1M_⊙ example than for the 1.4M_⊙ case. In principle one can think that this result could be due to the different quality of the input parameters in particular because KIC8006161 has a larger range of frequencies, including some ℓ = 3 modes, but we have checked that limiting the frequency modes of this star to a range similar to the observed one for KIC9139163 leads to very similar results. Hence it seems that the different positions in the HR-diagram should be the cause of the higher influence of the magnetic activity on the determination of the stellar parameters (relative to errors).

As long as surface uncertainties remain significant in the models, the frequency dependence of the frequencies shifts induced by the magnetic activity would be partially masked by the filtering procedures required in the modeling. However, as shown by Salabert et al. (2018), the frequency shifts can include an oscillatory component of a period similar to that corresponding to the acoustic depth of the He II zone, which would not be suppressed by the surface filtering. These shifts can also give rise to a misdetermination of some global stellar parameters, as for instance the helium abundance. But we have found that the uncertainties introduced seem to be below the 3% level in all the parameter analyzed.

Reducing the uncertainties on the frequencies to half of the observed one improves the determination of the stellar parameters in one of the cases (KIC8006161) and potentially would allow to flag incorrect determination of the best-fit model as they seem to have high χ² values.

FPH has made the simulations and do the theoretical work. RG, SM, AS, and CR have worked on the observational side and obtained the variations with activity of the oscillation modes. All the authors contributed to the discussion written in the article and to the conclusions.