Asteroseismic Investigation on KIC 10526294 to Probe Convective Core Overshoot Mixing

The convective overshoot mixing is an important mixing process in the stellar interior. The overshoot mixing outside a convective core significantly affects stellar evolution by refueling the burning core, prolonging the stellar lifetime, and resulting in a larger helium core. Observations of the extension of the main sequence in the HR diagram confirm the existence of such a mixing process in the stellar interior (see, e.g., Maeder & Mermilliod 1981, and the references therein). However, the precise mechanisms through which convective overshooting mixes stellar material remain unclear.

The classical overshoot models are nonlocal mixing-length theories (MLTs; e.g., Shaviv & Salpeter 1973; Maeder 1975; Bressan et al. 1981; Zahn 1991) that were developed to determine the penetrating distance, which refers to how far a fluid element can travel across the convective boundary. The penetrating distance represents the extended convective mixing range of the convective core, and it is expected that the convective core should be fully mixed with the extended (overshoot) region. However, the nonlocal mixing-length models have been found to be self-inconsistent, and the extent of overshooting depends on the model assumption, specifically the adiabatic or radiative stratification in the overshoot region (Renzini 1987). A nearly radiative stratification results in negligible overshooting, while a nearly adiabatic stratification leads to strong overshooting. In models with a nearly adiabatic stratified overshoot region (e.g., Zahn 1991), the temperature gradient undergoes an abrupt change from adiabatic to radiative. This abrupt variation in the temperature gradient can amplify the amplitude of oscillatory signals in frequencies originating from the base of the solar convection zone (Gough 1990). Analysis of solar p-mode frequencies has indicated a small overshoot length below the base of the solar convection zone (Basu & Antia 1994; Christensen-Dalsgaard et al. 1995). Furthermore, Christensen-Dalsgaard et al. (2011) found that the required temperature gradient profile is even smoother than the profile predicted by the solar model without overshoot, which cannot be produced by the nonlocal mixing-length overshoot models.

At present, another widely adopted treatment of the overshoot mixing is to handle it as a diffusion process with an exponential decreasing diffusion coefficient in the overshoot region (e.g., Freytag et al. 1982; Herwig 2000; Zhang & Li 2012a; Zhang 2013; Li 2017; Zhang et al. 2019). This exponentially decreasing behavior of the diffusion coefficient is based on numerical simulations (e.g., Freytag et al. 1982; Kupka et al. 2018) that show exponentially decreasing mean vertical convective velocity in the overshoot region, and stellar turbulent convection models (e.g., Xiong 1989; Xiong & Deng 2002; Zhang & Li 2012b; Li 2017; Xiong 2021) that show exponentially decreasing turbulent kinetic energy and turbulent dissipation. The stellar turbulent convection models are based on hydrodynamic equations, focusing on the turbulent variables (e.g., turbulent kinetic energy, heat flux, temperature variance, turbulent dissipation, and chemical abundance flux) in statistical equilibrium with some closure models of diffusion and dissipation (Xiong 1985, 1989; Canuto 1993, 1997; Canuto & Dubovikov 1998; Kupka 1999; Xiong & Deng 2002; Li 2007; Li 2012; Zhang & Li 2012b; Zhang 2013; Li 2017; Xiong 2021; Ahlborn et al. 2022; Kupka et al. 2022). The turbulent convection models are more reasonable and in better agreement than the classical overshoot models when compared with numerical simulations (e.g., Singh et al. 1995; Kupka 1999; Brummell et al. 2002; Kupka & Muthsam 2017). Recently, the properties and structure of the overshoot region predicted by the stellar turbulent convection model (Zhang & Li 2012b) have been found to be in good agreement with numerical simulations (Cai 2020a, 2020b, 2020c). Those models (e.g., Xiong & Deng 2002; Zhang & Li 2012a, 2012b; Xiong 2021) predict a smooth transition of the temperature gradient profile near the convective boundary, which is favored by the investigation of the solar p-mode frequencies (Christensen-Dalsgaard et al. 2011).

In the treatment of the overshoot mixing as a diffusion process, the diffusion coefficient is determined by an initial value multiplied by an exponentially decreasing factor. In the widely used model of Herwig (2000), the initial value of the diffusion coefficient is set to be a typical value obtained from the MLT in the convection zone near the convective boundary. On the contrary, Zhang (2013) predicted that the convective characteristic length in the overshoot region differs from that in the convection zone, resulting in a significant decrease of the diffusion coefficient near the convective boundary, followed by an exponential decrease in the bulk of the overshoot region (see, e.g., Meng & Zhang 2014). Recently, considering the potential rapid decrease of the diffusion coefficient near the convective boundary, Zhang et al. (2022) developed a model for overshoot mixing. This model sets the initial value of the diffusion coefficient as the typical value in the convection zone multiplied by a parameter C. This model recovers the model of Herwig (2000) when C = 1, while C ≪ 1 corresponds to a quick decrease of the diffusion coefficient near the convective boundary. Helioseismic and asteroseismic investigations (Zhang et al. 2022) on solar models and 13 Kepler solar-like stars have not placed any constraints on the parameter C. In this work, we plan to conduct asteroseismic investigations on stars exhibiting g-mode oscillations to explore the possibility of a rapid decrease.

Asteroseismology, which studies wave oscillations throughout the stellar interior, is a powerful tool for probing stellar structure and investigating physical processes within stars. Two types of asteroseismic investigations are commonly conducted to explore convective core overshoot. The first type investigates the ratios of small- to large-frequency separations in stellar p-mode oscillations, while the second type examines stellar g-mode and mixed-mode oscillations, as both are sensitive to the structure of the stellar core. Deheuvels & Michel (2011) provided detailed modeling of avoided crossings to constrain the stellar properties of HD 49385. Yang et al. (2015) analyzed the frequency separation ratios to identify a significant overshoot distance in KIC 2837475. Deheuvels et al. (2016) studied the ratio of small- to large-frequency separations in 24 CoRoT and Kepler low-mass stars and observed a trend of increasing overshoot region size with stellar mass. Moravveji et al. (2015) modeled the oscillations of the slowly pulsating B-type (SPB) star KIC 10526294 and found that the diffusion overshoot mixing model of Herwig (2000) provided a better fit to the observed oscillation periods than the classical nonlocal mixing-length overshoot model. Similar results were obtained by Moravveji et al. (2016) in their asteroseismic investigation of another intermediate-rotating SPB star, KIC 7760680, and they suggested that KIC 10526294 and KIC 7760680 are ideal targets for testing theories of convection. Wu & Li (2019) examined the oscillations of the SPB star HD 50230 and constrained the stellar parameters and overshoot mixing. Wu et al. (2020) also investigated the SPB star KIC 8324482 and identified a weak convective core overshoot mixing. Michielsen et al. (2019) demonstrated the capability of asteroseismology in probing thermal stratification in the overshoot region. Recently, Michielsen et al. (2021) analyzed the SPB star KIC 7760680 and found that a radiative stratification in the overshoot region is more favorable than an adiabatic-dominated stratification. Noll et al. (2021) investigated KIC 10273246 and traced the evolution of mixed-mode frequencies to constrain core overshoot. They found that the data could not distinguish between exponential decreasing diffusion and constant diffusion in the overshoot region. Pedersen et al. (2014) modeled 26 SPB stars with different internal mixing process patterns and investigated them using asteroseismology. They found that internal mixing profiles in SPB stars exhibit radial stratification. No fixed pattern emerged as the best fit for all target stars, and exponential decreasing diffusion of overshoot did not yield better results than constant diffusion of overshoot.

In this paper, we conduct an asteroseismic investigation of the SPB star KIC 10526294 to explore core overshoot, particularly focusing on the rapid decrease of the diffusion coefficient near the convective boundary. Section 2 provides a brief introduction to our target, KIC 10526294. In Section 3, we describe the overshoot mixing model adopted in detail. The input physics used in the calculations of the stellar models are outlined in Section 4. The numerical results are presented in Sections 5 and 6. We discuss the obtained results in Section 7. Finally, we summarize our conclusions in Section 8.

KIC 10526294 has been classified as an SPB star, and a detailed analysis of its 4 yr light curve revealed an oscillatory period range of approximately 1–2 days, including 19 individual eigenfrequencies, with a mean period spacing of about 0.06 day (Pápics et al. 2014). The analysis of its spectra provided the atmospheric parameters of the star: an effective temperature of T_eff = 11550 ± 500 K, a surface gravitational acceleration of $\mathrm{log}g=4.1\pm 0.2$ (in cgs units), and a metallicity of $Z={0.016}_{-0.007}^{+0.013}$ (Pápics et al. 2014). The rotation period of KIC 10526294 is estimated to be approximately 188 days based on the average splitting of stellar oscillation frequencies (Pápics et al. 2014). The internal rotation profile of KIC 10526294 has been analyzed by Triana et al. (2015) using frequency splitting (Aerts et al. 2010). They found that the average rotation rate near the core overshoot region is approximately 163 ± 89 nHz, and a mild counter-rotating region in the envelope toward the surface is preferred when assuming a smooth profile of rotation rate in the stellar interior.

Asteroseismic investigations based on standard stellar models (Pápics et al. 2014) have indicated that KIC 10526294 is a young star with a stellar mass of about 3.2 solar masses, and it has an upper limit for core overshoot of α_ov ≤ 0.15 for a classical full mixing model or f_ov ≤ 0.015 for the diffusion overshoot model of Herwig (2000). Further comprehensive asteroseismic investigations on stellar models of KIC 10526294 have been conducted by Moravveji et al. (2015), and subsequent investigations based on a statistical approach have been performed by Aerts et al. (2018) and Pedersen et al. (2014). Moravveji et al. (2015) examined stellar models on a denser grid, considering an additional global mixing, testing different compositions, and comparing classical fully overshoot mixing with the diffusion overshoot model of Herwig (2000). They found that the effects of the additional global mixing are significant in fitting an oscillation-like pattern of the period spacing. The χ² value of the best stellar model with the diffusion overshoot model is less than half of the χ² value of the best stellar model with classical overshoot fully mixing, indicating that the diffusion overshoot model is more reasonable than the classical overshoot models. All of these analyses of KIC 10526294 highlight it as an excellent target for studying overshoot and testing the theory of diffusion overshoot mixing. Therefore, in this work, we continue to model KIC 10526294 using the diffusion overshoot mixing approach and aim to investigate the possible rapid variation of the diffusion coefficient near the convective boundary. We hope that such an analysis can provide further insights into the current developments of stellar convection theory.

The adopted model for convective core overshoot mixing is based on an exponentially decreasing diffusion coefficient, as developed by Zhang et al. (2022). The model is described by the following equation:

$$\begin{eqnarray}&&{D}_{\mathrm{ov}}={{CD}}_{0}{\left(\displaystyle \frac{P}{{P}_{\mathrm{CB}}}\right)}^{\theta },\end{eqnarray} \tag{ 1 }$$

where D_ov represents the diffusion coefficient for overshoot mixing, P denotes pressure, P_CB represents the pressure at the boundary of the convective core, D₀ corresponds to the typical diffusion coefficient near the boundary of the convective core, and C and θ are parameters that characterize the properties of the model. The specific details of these parameters will be introduced later. To determine D₀, the typical value of the diffusion coefficient near the boundary of the convective core, the MLT is employed. In this approach, D₀ is taken as D(r_*), where r_* is a location inside the convective core chosen to avoid the local properties of MLT. In this case, r_* = r_CB − d, with r_CB representing the convective boundary and d denoting a nonzero distance, which is set to 0.1 H_P by default. A graphical representation of the diffusion coefficient is presented in Figure 1.

Zoom In Zoom Out Reset image size

The meanings of the model have been discussed in detail in Zhang et al. (2022). Here, we summarize the main properties of the model. The parameter θ represents the factor responsible for the exponential decrease in D, primarily determining the effective length of the overshoot mixing. The parameter C describes the potential quick variation of the diffusion coefficient near the convective boundary. When C ≪ 1, it corresponds to a significant rapid decrease in D. This rapid decrease is assumed to occur within a very thin layer, and we do not trace the precise details of this rapid decline. Consequently, combining Equation (1) with D in the convection zone leads to a discontinuity at the boundary. C also secondarily determines the effective length of the overshoot mixing. The model reverts to the widely used formulation of Herwig (2000) when C = 1 and θ = 2/f_ov. Here, C = 1 signifies the absence of a quick variation in D near the boundary, which is an implicit assumption in the formula proposed by Herwig (2000).

Stellar models and oscillations are calculated using the YNEV code (Zhang 2015). The AGSS09 (Asplund et al. 2009) solar composition, revised for Ne (Young 2018), is adopted for elements heavier than He in the stellar models. This composition is essentially the same as the recently revised solar composition proposed by Asplund et al. (2021). The opacity in the stellar interior is interpolated from the OPAL opacity tables (Iglesias & Rogers 1996) based on the adopted composition. The low-temperature opacity tables (Ferguson et al. 2005) are smoothly connected to the OPAL tables. The OPAL equation-of-state tables (Rogers & Nayfonov 2002) are used to interpolate pressure P, its derivatives with respect to density and temperature ${(\partial \mathrm{ln}P/\partial \mathrm{ln}\rho )}_{T}$, ${(\partial \mathrm{ln}P/\partial \mathrm{ln}T)}_{\rho }$, and the adiabatic temperature gradient ${{\rm{\nabla }}}_{\mathrm{ad}}={(\partial \mathrm{ln}T/\partial \mathrm{ln}P)}_{S}$. Other quantities such as the adiabatic compression index Γ₁ and specific capacities c_V , c_P are obtained using thermodynamic relations. Nuclear reaction rates are based on NACRE (Angulo et al. 1999) and are enhanced by weak screening (Salpeter 1954). The temperature gradient in the convective core/zones is calculated using the MLT with α = 1.8, which is calibrated in the solar model. Since the envelope of a main-sequence 3M ⊙ star is radiative-dominated, the effects of varying α are not significant. The convective heat flux in the overshoot region is ignored. Apart from the convective overshoot mixing and a global diffusion with $\mathrm{log}{D}_{\mathrm{ext}}=1.75$, no other nonstandard physical processes such as microscopic diffusion, mass loss, or accretion are taken into account.

In order to reduce the number of stellar models, we adopt a universal helium enhancement law given by Y = Y_p + (ΔY/ΔZ)Z for the initial stellar abundances. Here, Y_p represents the primordial helium abundance, and we adopt the value Y_p = 0.2485 based on the Big Bang nucleosynthesis result (Cyburt et al. 2008). The helium enhancement ratio ΔY/ΔZ is set to 2.0 by comparing the enhancement law with the standard chemical mixture observed in OB stars in the solar neighborhood (Y = 0.276, Z = 0.014; Nieva & Przybilla 2012). However, there is evidence suggesting that the linear helium enhancement law might be too simplistic, and the helium enhancement ratio could vary over a wide range (see, e.g., Verma et al. 2019, and references therein). In asteroseismic investigations of SPB stars, an additional global diffusion process has been found to be widely present in order to explain the observed period spacing patterns (see, e.g., Pedersen et al. 2014; Moravveji et al. 2015, 2016; Wu & Li 2019; Wu et al. 2020). For KIC 10526294, we adopt a global extra diffusion with $\mathrm{log}{D}_{\mathrm{ext}}=1.75$ in our calculations, based on the results obtained by Moravveji et al. (2015). We will discuss the effects of varying these two constraints, namely the reduced freedom in the initial stellar abundances and the fixed value of $\mathrm{log}{D}_{\mathrm{ext}}$.

All stellar models in this work evolve from the pre-main-sequence (PMS) stage with a central temperature of T_C = 10⁵ K to the end of the main sequence, defined as the point when the central hydrogen abundance decreases to 0.01. The stellar models typically consist of around 2000 mesh points. The time step is controlled to ensure that the variation of hydrogen abundance in the stellar core is less than 0.005 within a single time step. As a result, there are approximately 200 models along the evolutionary track.

We have calculated 11,988 stellar evolutionary tracks using a coarse grid with varied stellar parameters covering various ranges, as listed in Table 1. This coarse grid, with relatively low resolutions of the stellar parameters, provides an initial exploration of the parameter space where the computed stellar models exhibit good agreement with observations. Subsequently, a second step involves calculating additional models with a finer grid based on the preliminary parameter space to search for the best-fitting models.

Table 1. Stellar Parameters and Their Ranges for the Grids of Stellar Models

	From	To	Step	Points
coarse grid				11988
stellar mass M (M⊙)	2.70	3.60	0.05	37
metallicity Z	0.008	0.040	0.008	9
$\mathrm{log}C$	−6	0	2	4
$\mathrm{log}\theta$	1.4	3.0	0.2	9
fine grid				233393
S = M + 16Z	3.40	3.70	0.01	31
Z	0.008	0.040	0.002	17
$\mathrm{log}C$	−6	0	0.5	13
$\mathrm{log}\theta$	1.30	2.30	0.02	∼20

Download table as:  ASCIITypeset image

A stellar evolutionary track is determined by the stellar parameters at each mesh point on the grid. Along each track, the oscillation frequencies depend on the stellar age t. We evaluate the stellar models using the χ² statistic defined as follows:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{{N}_{\mathrm{obs}}}\displaystyle \sum _{k=1}^{{N}_{\mathrm{obs}}}{\left(\displaystyle \frac{{P}_{k}-{P}_{k,\mathrm{obs}}}{\sigma {P}_{k,\mathrm{obs}}}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where N_obs = 19 represents the number of observed stellar oscillatory periods. The observed periods and their uncertainties are obtained from Moravveji et al. (2015). As a stellar evolutionary track is determined by the stellar parameters, and the interior structure of a star on the track is determined by its age, the value of χ² becomes a function of ($\mathrm{log}C$, $\mathrm{log}\theta$, M, Z, t). Along a given evolutionary track, we can identify the best-fitting model at a particular age, denoted as t_m , which reproduces the observations with the minimum value of χ² (${\chi }_{m}^{2}$). Both ${\chi }_{m}^{2}$ and t_m depend on the stellar parameters ($\mathrm{log}C$, $\mathrm{log}\theta$, M, and Z).

No significant dependence of ${\chi }_{{\rm{m}}}^{2}$ on individual parameters C, M, or Z has been observed in the coarse grid. However, as depicted in the top panel of Figure 2, χ² exhibits a dependence on a linear combination of M and Z. The optimal range for M and Z is indicated by the two dashed lines, defined as 3.4 ≤ S = M/M_⊙ + 16Z ≤ 3.7. The dependence of χ² on $\mathrm{log}C$ and $\mathrm{log}\theta$ is shown in the bottom panel of Figure 2. Mesh points that lie below or near the dashed line, represented by $\mathrm{log}\theta =0.05\mathrm{log}C+2.3$, generally yield lower χ² values, suggesting an optimal range for these two parameters. However, two mesh points, specifically ($\mathrm{log}\theta =1.6$, $\mathrm{log}C=-1$) and ($\mathrm{log}\theta =1.8$, $\mathrm{log}C=0$), show local minima of χ². Considering the sparsity of the parameter space in the coarse grid, we decided to include all points below the dashed line as the optimal region. This approach allows us to further constrain the models using a fine grid generated based on the two criteria outlined in Figure 2. To recap, the optimal parameter space is summarized as 3.40 ≤ S ≤ 3.70 and $1.3\leqslant \mathrm{log}\theta \leqslant 0.05\mathrm{log}C+2.3$. We computed a total of 233,393 stellar evolutionary tracks in the fine grid, covering the specified parameter ranges with small steps for each stellar parameter, as listed in Table 1. In the fine grid, S, Z, $\mathrm{log}C$, and $\mathrm{log}\theta$ are evenly spaced. Since S is a function of stellar mass, the step size for M is not fixed, and the range of M in the fine grid is from 2.76–3.572. Subsequently, we constrained the range of each stellar parameter (M, Z, $\mathrm{log}\theta$, $\mathrm{log}C$) and some global variables (age, effective temperature, radius, central hydrogen abundance, etc.) of the stellar models by analyzing the distribution of ${\chi }_{{\rm{m}}}^{2}$ for each of these parameters.

Zoom In Zoom Out Reset image size

6.1. Suggested Ranges of Stellar Parameters and Global Variables

The minimum value of ${\chi }_{{\rm{m}}}^{2}$ among all of the evolutionary tracks in the fine grid is ${\chi }_{\mathrm{best}}^{2}=495$. Using an interpolation of oscillation frequencies (to be described later) with respect to the stellar parameters $\mathrm{log}\theta$, M, Z, and stellar age, we identified the best model with ${\chi }_{\mathrm{best}}^{2}=442$. Following the approach of Zhang et al. (2022), where it is assumed that χ² values outside the suggested range of stellar parameters should be significantly larger than those within the suggested range, we defined the suggested range as ${\chi }_{{\rm{m}}}^{2}\lt 2{\chi }_{\mathrm{best}}^{2}$. Therefore, the range of stellar parameters satisfying $\mathrm{log}{\chi }_{{\rm{m}}}^{2}\lt 2.95$ is considered the suggested range.

In Figure 3(a), it is shown that the suggested mass range for KIC 10526294 is 2.9–3.4 M_⊙, with the best value being approximately 3.0 M_⊙ based on the interpolation of oscillation frequencies. Within the considered range of 0.008 ≤ Z ≤ 0.040, Figure 3(b) demonstrates that the minimum value of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ increases as Z increases, with the minimum almost reaching the upper limit of the suggested range, 2.95, at Z = 0.010. Therefore, the suggested range for metallicity is Z ≥ 0.01. However, a higher Z-value larger than 0.02 is more likely since the minimum value of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ is lower and exhibits minimal changes for Z ≥ 0.02. This preference for a higher Z-value is also supported by the nonadiabatic investigation on mode excitation (Moravveji et al. 2015).

Zoom In Zoom Out Reset image size

In Figure 3(b), it is evident that the model with the lowest ${\chi }_{m}^{2}$ is obtained when Z = 0.04. The reason for not using a higher upper limit for Z is that the OPAL EOS tables implemented in the YNEV code are only valid for 0 ≤ Z ≤ 0.04. Hence, models with Z > 0.04 are not included in this analysis to avoid extrapolations on EOS variables. Although the best model has Z = 0.04, it can be observed from Figure 3(b) that the minimum value of ${\chi }_{{\rm{m}}}^{2}$ changes little as Z increases beyond 0.03. Consequently, setting a higher value for Z is unlikely to significantly reduce ${\chi }_{{\rm{m}}}^{2}$. Moreover, the spectral results by Pápics et al. (2014) indicate that a Z-value exceeding 0.04 is unlikely for KIC 10526294.

The two-variable investigation on M and Z is depicted in the top panel of Figure 4. The color of each data point represents the minimum value of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ for fixed M and Z and different $\mathrm{log}C$ and $\mathrm{log}\theta$. By applying the criterion ${\chi }_{{\rm{m}}}^{2}\lt 2{\chi }_{\mathrm{best}}^{2}$, a tight constraint on M and Z is revealed as follows:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{M}{{M}_{\odot }}+16Z & = & 3.50\pm 0.03\quad (\mathrm{for}\ Z\geqslant 0.016)\\ \displaystyle \frac{M}{{M}_{\odot }}+26Z & = & 3.66\pm 0.03\quad (\mathrm{for}\ Z\lt 0.016).\end{array}\end{eqnarray} \tag{ 3 }$$

This constraint aligns with the findings of Moravveji et al. (2015). Negative correlations between mass and metallicity have been observed in asteroseismic modeling of other intermediate-mass stars as well, such as in Ausseloos et al. (2004) and Wu et al. (2020). A plausible explanation, as suggested by Wu et al. (2020), is that the combination of parameters M and Z should result in an appropriate convective core to match the observed period spacing pattern. In comparison to Moravveji et al. (2015), the suggested mass in this study is slightly lower. This discrepancy is likely due to the wider range of Z-values employed in our fine grid. The constraint on M and Z indicates a lower stellar mass for higher Z. In Moravveji et al. (2015), the maximum Z-value in their fine grid is 0.02, which is lower than ours. Their best model exhibits a stellar mass of approximately 3.25 M_⊙ and Z = 0.014, consistent with the constraint in Equation (3) discovered in this work.

Zoom In Zoom Out Reset image size

The overshoot parameter C is of particular interest, and thus we conducted a detailed investigation of the best stellar model for a given $\mathrm{log}C$ by interpolating the oscillation frequencies. Interpolation of the frequencies from the calculated stellar models is necessary for the given values of $\mathrm{log}C$, $\mathrm{log}\theta$, M, Z, and the reduced stellar age τ = t/t_ms (where t_ms represents the stellar lifetime on the main sequence). The interpolation process involves two steps. In the first step, we interpolate the frequencies at an arbitrary τ_* on a stellar evolutionary track determined by the stellar parameters at a mesh point on the grid. Hermite interpolation with piecewise cubic polynomials is employed for this purpose. The derivatives of the frequencies with respect to τ are calculated using a five-point Taylor expansion. The second step involves interpolating the frequencies using arbitrary $\mathrm{log}{\theta }_{* }$, S_*, Z_*, and τ_*. For any given point $(\mathrm{log}{\theta }_{* },{S}_{* },{Z}_{* })$ within the considered range of the fine grid, we identify the elemental cell that contains it. We then interpolate the frequencies at τ_* on the eight stellar evolutionary tracks corresponding to the eight apexes of the cell, following the first step. Subsequently, trilinear interpolation is performed on $\mathrm{log}\theta$, S, and Z to obtain the frequencies with the parameters $\mathrm{log}{\theta }_{* }$, S_*, Z_*, and the reduced age τ_*. This interpolation scheme ensures the continuity of frequencies throughout the entire space of stellar parameters and ages. By iterating based on the continuously interpolated frequencies, we can determine the minimum value of $\mathrm{log}{\chi }^{2}$ for a given $\mathrm{log}C$. To assess the accuracy of the interpolation scheme, we compare the differences in oscillation periods δ P and period spacings δΔP between the results of interpolation and those obtained from actual calculations of stellar evolution and oscillation for the best model. The differences, displayed in Figure 5, are at a level of approximately 1 s, confirming the validity of the interpolation scheme.

Zoom In Zoom Out Reset image size

The distribution of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ for the mesh points on the grid with respect to $\mathrm{log}C$ is presented in Figure 3(d). The empty circles indicate the minimum values of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ for given $\mathrm{log}C$ obtained using the interpolation method described earlier. It can be observed that the minima of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ for all considered $\mathrm{log}C$ values are below 2.95, and the minimum value of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ at $\mathrm{log}C=-6$ is close to 2.95. Thus, the suggested range for $\mathrm{log}C$ is $\mathrm{log}C\geqslant -6$. An important result from Figure 3(d) is that the overall minimum of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ does not occur at $\mathrm{log}C=0$, but rather around $\mathrm{log}C=-4$, with minimal variation in the range $-5\leqslant \mathrm{log}C\leqslant -3$. This finding suggests the existence of a rapid decrease in the diffusion coefficient in a thin layer near the convective boundary, which is a crucial aspect of the adopted overshoot mixing model. The diffusion mixing model proposed by Herwig (2000), corresponding to $\mathrm{log}C=0$, does not yield the best result. Another noteworthy observation is that overshoot mixing must be taken into account because the minima of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ rapidly increase as $\mathrm{log}C$ decreases when $\mathrm{log}C\lt -5$. As the overshoot mixing is proportional to C, very small values of C indicate negligible overshoot.

However, it should be noted that $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ is not highly sensitive to changes in $\mathrm{log}C$. The reduction in ${\chi }_{{\rm{m}}}^{2}$ is only approximately 30% when $\mathrm{log}C$ varies from 0 to −4. According to the adopted criterion, which sets ${\chi }_{{\rm{m}}}^{2}\lt 2{\chi }_{\mathrm{best}}^{2}$ as the suggested range, C = 1 still falls within the suggested range. Therefore, the results provide weak evidence for the rapid decrease of the diffusion coefficient near the convective boundary. Further investigations involving a larger sample size are required to confirm this observation.

The results of the two-variable investigation on $\mathrm{log}C$ and $\mathrm{log}\theta$ are presented in the bottom panel of Figure 4. Based on Figure 3(c) and the criterion ${\chi }_{{\rm{m}}}^{2}\lt 2{\chi }_{\mathrm{best}}^{2}$, the suggested range for $\mathrm{log}\theta$ is $1.7\geqslant \mathrm{log}\theta \geqslant 2.1$. In the classical diffusion mixing model, θ is linked to f_ov through the relation f_ov θ = 2. Consequently, the suggested range for f_ov is from 0.016–0.040. It is important to note that the broad range of f_ov is due to the relatively loose constraint on C. However, if we focus on the case of C = 1, which corresponds to the diffusion mixing model proposed by Herwig (2000), the result depicted in Figure 4 indicates a tight constraint on $\mathrm{log}\theta$ as 2.12 ± 0.01 for $\mathrm{log}{\chi }_{{\rm{m}}}^{2}\lt 2.95$, corresponding to f_ov = 0.016. This finding is consistent with the results reported by Moravveji et al. (2015) and also aligns with the recommendation made by Herwig (2000).

The distributions of $\mathrm{log}{\chi }_{{\rm{m}}}^{2}$ are shown in Figure 3(e)–(i) for various stellar global variables, including the stellar age t, effective temperature T_eff, radius R, center hydrogen abundance X_C , period spacing ΔP, mass fraction of the convective core mass M_C /M, center temperature T_C , and center density ρ_C . The suggested range and the best value for these variables are listed in Table 2. Additionally, the suggested range for the surface gravitation acceleration is determined to be $4.15\leqslant \mathrm{log}g\leqslant 4.25$ in cgs units. The results obtained for T_eff and $\mathrm{log}g$ are consistent with the spectral analysis conducted by Pápics et al. (2014). It is worth noting that the period spacing is the global variable that exhibits the tightest constraint. This result is expected, as high-order g-mode oscillations tend to exhibit nearly uniform spacing in periods.

Table 2. Suggested Range and the Best Values of the Stellar Parameters and Global Variables, in cgs Units

Global Variable	Suggested Range	Best
M/M⊙	2.9 ∼3.4	3.0
Z	0.01 ∼0.04	0.04
$\mathrm{log}\theta$	1.7 ∼2.1	1.9
$\mathrm{log}C$	≥-6	−4
t/Gyr	0.05 ∼0.09	0.08
$\mathrm{log}{T}_{\mathrm{eff}}$	4.02 ∼4.16	4.04
R/R⊙	2.12 ∼2.38	2.36
XC	0.56 ∼0.65	0.57
ΔP/s	5500 ∼5600	5570
MC /M	0.18 ∼0.22	0.19
$\mathrm{log}{T}_{C}$	7.37 ∼7.41	7.37
$\mathrm{log}{\rho }_{C}$	1.52 ∼1.56	1.52

Download table as:  ASCIITypeset image

6.2. The Differences between the Models and the Observations

The minimum value of χ² is still significantly larger than unity, indicating a substantial residual in the oscillation periods. In Figure 6, the period spacings of the lowest-χ² model are represented by black dots, while the observed period spacings are denoted by gray empty circles. A comparison between the two reveals that the model's period spacing exhibits a generally smooth pattern, whereas the observations display irregular oscillatory behavior. The relative differences between the observed oscillation periods and those of the lowest-χ² model are depicted by the gray symbols in Figure 7. This plot demonstrates that the irregular oscillatory behavior of the period spacings arises from similar variations in the oscillation periods themselves. To gain insight into the observed pattern of oscillation periods, it is useful to investigate potential adjustments to the stellar models by analyzing the residuals between a reference stellar model and the observed periods.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Cowling approximation simplifies the stellar adiabatic oscillation equations to second-order equations, allowing for analytical investigations using the JWKB approximation for high-order modes (see, e.g., Gough 1993; Montgomery et al. 2003). By combining this approximation with the asymptotic expression of the g-mode period (Tassoul 1980), given by P_k ≈ 2π²Π₀(k + k₀)/L (where k represents the mode node and $L=\sqrt{l(l+1)}$), and employing the variational method, we can derive the variations in the oscillation period resulting from a perturbation of the Brunt-Väisälä frequency N as

$$\begin{eqnarray}&&\displaystyle \frac{\delta {P}_{k}}{{P}_{k}}\approx -{\displaystyle \int }_{0}^{1}z(u)\{1+\sin [2\pi (k+{k}_{0})u]\}{du},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&z=\displaystyle \frac{\delta N}{N},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&x=\displaystyle \frac{r}{R},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{{\rm{\Pi }}}_{x}}^{-1}={\int }_{0}^{x}\displaystyle \frac{\left|N\right|}{x}{dx},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{0}={{\rm{\Pi }}}_{x}(x=1),\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&u=\displaystyle \frac{{{{\rm{\Pi }}}_{x}}^{-1}}{{{{\rm{\Pi }}}_{0}}^{-1}}.\end{eqnarray} \tag{ 9 }$$

Equation (4) differs slightly from the result presented in Miglio et al. (2009). The detailed derivation of Equation (4) and its numerical validation are provided in the Appendix. It has been shown that this equation accurately reproduces both the qualitative and quantitative aspects of the numerical results concerning the variation of the oscillation period.

Considering a perturbation of N in a single region with u₁ ≤ u ≤ u₂, let z(u₁) = z(u₂) = 0 for continuity and expand z to Fourier sine series:

$$\begin{eqnarray}&&z(u)=\left\{\begin{array}{c}0,u\lt {u}_{1}\\ \displaystyle \sum _{n=1}^{\infty }{A}_{n}\sin \left(\displaystyle \frac{u-{u}_{1}}{{u}_{2}-{u}_{1}}n\pi \right),{u}_{1}\leqslant u\leqslant {u}_{2}\\ 0,u\gt {u}_{2}\end{array}\right.,\end{eqnarray} \tag{ 10 }$$

where A is the coefficients of the sine series. The variation of oscillation period based on Equation (4) is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\delta {P}_{k}}{{P}_{k}}\approx -\displaystyle \frac{({u}_{2}-{u}_{1})}{\pi }\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{A}_{n}}{n}\left\{[1-{\left(-1\right)}^{n}]\right.\\ \left.+\,\displaystyle \frac{\sin [2\pi (k+{k}_{0}){u}_{1}]-{\left(-1\right)}^{n}\sin [2\pi (k+{k}_{0}){u}_{2}]}{1-{\left[\tfrac{2(k+{k}_{0})({u}_{2}-{u}_{1})}{n}\right]}^{2}}\right\}.\end{array}\end{eqnarray} \tag{ 11 }$$

When a perturbation of N exists in many detached regions, the total δ P_k /P_k is the sum of δ P_k /P_k in each region.

By utilizing the above formula, it becomes possible to investigate the potential modifications to the N profile that may correspond to the characteristics of the oscillatory component in the observed period spacing. Using a reference model to calculate P_k , the objective is to find the values of u₁, u₂, and A_n for all regions that minimize the following function:

$$\begin{eqnarray}&&Q=\displaystyle \sum _{k}{\left(\displaystyle \frac{\delta {P}_{k}}{{P}_{k}}-\displaystyle \frac{{P}_{k,\mathrm{obs}}-{P}_{k}}{{P}_{k}}\right)}^{2}.\end{eqnarray} \tag{ 12 }$$

In this equation, δ P_obs,k is not used as the denominator since the uncertainty arising from Equation (4) based on the asymptotic relation is not comparable to the very small observational uncertainties. Considering δ P_obs,k as the denominator would introduce significant random variations to Q, rendering it unhelpful for analysis. Despite not incorporating δ P_obs,k , the above equation is sufficient to capture the observational characteristics of the period spacing through variations in the reference model. In this problem, all A_n can be determined using the linear least-squares method; thus, only u₁ and u₂ for all regions need to be determined.

In our case, we have selected the stellar model with the following parameters as the reference model: $\mathrm{log}C=-4$, $\mathrm{log}\theta =1.932$, M = 2.983M_⊙, Z = 0.03171, and τ = 0.1942. For this model, k₀ is fixed at 1.0. Alternatively, (k + k₀) in Equations (4) and (11) can be replaced with P_k /ΔP, where ΔP = 2π²Π₀/L represents the asymptotic period spacing. The differences between the two approaches are negligible. Perturbations on N within a single region are assumed, and only the lowest-order term in the Fourier sine series is taken into account. By scanning through different values of u₁ and u₂, it was found that the profile of δ N/N that minimizes Q is determined by A₁ = 0.0755 within the range 0.2273 ≤ u ≤ 0.2388. The revised profile of N is illustrated in Figure 8. It should be noted that the three observational modes with the highest periods are not taken into account, as their inclusion would result in Fourier sine series coefficients that are too large and violate the assumption of perturbation.

Zoom In Zoom Out Reset image size

The results of the relative variations in the oscillation periods with the δ N/N profile are depicted in Figure 7. In this context, P represents the oscillation period of the reference model, P_rev signifies the oscillation period of the reference model with the revised N profile, and δ P = P_rev − P. The revised N profile leads to a reduction in the overall residual of the oscillation periods, i.e., ${\sum }_{k}\left|{P}_{\mathrm{rev},k}-{P}_{\mathrm{obs},k}\right|=3208\,{\rm{s}}$, which is lower than ${\sum }_{k}\left|{P}_{k}-{P}_{\mathrm{obs},k}\right|=3707\,{\rm{s}}$. However, the revised N profile yields a higher χ² = 736 compared to χ² = 692 for the reference model. The analytical approximation Equation (4) is not accurate enough to be comparable with the very small observational uncertainties (ΔP_obs/P_obs is less than ~10⁻⁴), so that χ² is out of control. The differences in oscillation periods between observations and the reference model, P_obs/P − 1, are illustrated by the gray symbols, while the error bars are multiplied by three to enhance their visibility. The analytical δ P/P based on Equation (11) is represented by the blue symbols, and the numerical δ P/P calculated using the oscillation code is shown by the red symbols. The numerical δ P/P generally reproduces the analytical results, validating the reasonability of the approximation given by Equation (4). However, the numerical δ P/P does not perfectly replicate P_obs/P − 1 since the former exhibits regular oscillatory behavior while the latter is irregular. This indicates that the irregular oscillatory behavior observed in the period spacings cannot be attributed to a simple sine function of δ N/N in a single region. Nevertheless, the numerical δ P/P qualitatively reveals the oscillatory nature of the observational period spacings, as the positions of the peaks and valleys in the numerical δ P/P align with those in P_obs/P − 1. Similar results can be observed in the period spacings of the model with the revised N profile, illustrated in Figure 6 as black triangles connected by a dashed line. It can be found that the model with the revised N profile exhibits an oscillatory component, and its period spacings are closer to the observations than those of the reference model within the range of 1.3 < P/day < 1.9. It is evident from Equation (4) that the period (measured by the node number k as shown in Figure 7) of the oscillatory behavior is sensitive to the location of the revised N profile. Therefore, the oscillatory behavior observed in the period spacings suggests the existence of an additional abundance gradient near u = 0.23 or r = 0.21R.

The above investigation suggests that the stellar models lack additional physical processes that could cause variations in N near r = 0.21R. Some potential mechanisms directly associated with this matter include mixing and inhomogeneous accretion during the PMS phase. Exploring these possible mechanisms would require extensive calculations of numerous stellar models. However, investigating these aspects is beyond the scope of this paper, and we intend to address them in future research.

7.1. Implications on the Convective Overshoot

7.2. Effects of Varied Initial Abundance and Extra Global Mixing

In this study, we examined convective core overshoot mixing using asteroseismology of the SPB star KIC 10526294. The selected overshoot mixing model is based on the work of Zhang et al. (2022). The power-law diffusion coefficient is derived from numerical simulations and turbulent convection models, which provide a power-law description of the turbulent velocity (as discussed in Section 1). The parameter C is determined based on the theoretical model of convective overshoot mixing presented by Zhang (2013), which predicts a different characteristic length in the overshoot region compared to the convection zone. The chosen overshoot mixing model aligns with the widely used model proposed by Herwig (2000) when the parameter C is set to 1.

The oscillation periods of the stellar models for KIC 10526294, generated by varying stellar parameters, are compared with the observed periods. Initially, a coarse grid with large parameter steps is employed to explore the preliminary parameter space. Subsequently, a finer grid with smaller parameter steps within the preliminary parameter space is used for a detailed analysis to constrain the overshoot parameters and stellar parameters. To address the issue of resolution in stellar parameters and age (e.g., Wu & Li 2016), an interpolation technique is applied to the oscillation frequencies, considering all parameters and stellar age. Hermite cubic polynomial interpolation is utilized for stellar age along a stellar evolutionary track, while linear interpolation is employed for stellar parameters and overshoot parameters. The estimated error in the oscillation period resulting from the interpolation scheme is within a few seconds.

Based on the comparison between the model oscillation periods and the observations, the stellar parameters and overshoot parameters are constrained. The suggested parameter ranges are listed in Table 2. While there are slight differences in the stellar parameters compared to the results of Moravveji et al. (2015), our models tend to prefer a lower stellar mass. This discrepancy may be attributed to the broader range considered for the stellar parameters Z. The minimum value of χ² in our results is approximately 442, which is lower than that obtained by Moravveji et al. (2015). The improvement can be attributed to two main factors: the incorporation of the parameter C variation and the interpolation of the oscillation frequencies using the parameters and stellar age.

Although the recommended range for the parameter C still includes unity, the most intriguing finding is that the best value of C is 10⁻⁴, significantly smaller than unity. This result could be viewed as weak evidence supporting the prediction of the turbulent convective mixing model proposed by Zhang (2013), which suggests a rapid decrease in the diffusion coefficient of convective mixing near the convective boundary.

Furthermore, an analysis of the residuals in the oscillation periods between a reference model and the observations is conducted. A formula describing the variation of the oscillation periods due to perturbations in the Brunt-Väisälä frequency is derived and quantitatively validated. The location of the perturbation in the N profile can be determined by comparing the model and observational periods, utilizing Fourier expansion for the perturbation and the linear least-squares method for the expansion coefficients. It is found that the oscillatory behavior in the residuals of the period spacing may be related to anomalies in the N profile. This method is expected to be applicable for analyzing the residuals of oscillation periods in general, comparing models with observations.

The authors offer many thanks to the anonymous referee for providing valuable comments, which improved this work. Many thanks to OpenAI (2020) for checking and correcting some syntactical errors in the original manuscript. Fruitful discussions with Prof. Jørgen Christensen-Dalsgaard are highly appreciated. Funding for Yunnan Observatories is co-sponsored by the National Key R&D Program of China (grant No. 2021YFA1600400 / 2021YFA1600402), the Strategic Priority Research Program of the Chinese Academy of Sciences (grant No. XDB 41000000), the National Natural Science Foundation of China (grant Nos. 11773064, 11873084, 12133011, 12273104, and 12288102), the Foundation of the Chinese Academy of Sciences (Light of West China Program and Youth Innovation Promotion Association), and the Yunnan Ten Thousand Talents Plan Young & Elite Talents Project, and International Centre of Supernovae, Yunnan Key Laboratory (No. 202302AN360001).

Equation (4) is obtained using the variational method for a perturbation on the linear operator of the stellar oscillation. With the Cowling approximation and the JWKB approximation, the variational method shows that the variation of the g-mode oscillation frequency is given by (see, e.g., Gough 1993; Montgomery et al. 2003)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\delta {\omega }_{k}}{{\omega }_{k}} & \approx & \displaystyle \frac{2}{{\displaystyle \int }_{0}^{R}\left|N\right|{dr}^{\prime} /r^{\prime} }{\displaystyle \int }_{0}^{R}\displaystyle \frac{\delta N}{N}\displaystyle \frac{N}{r}{dr}{\sin }^{2}\\ & & \times \,\left(\displaystyle \frac{{P}_{k}}{{\rm{\Delta }}P}\pi \displaystyle \frac{{\displaystyle \int }_{0}^{r}\left|N\right|{dr}^{\prime} /r^{\prime} }{{\displaystyle \int }_{0}^{R}\left|N\right|{dr}^{\prime} /r^{\prime} }+\displaystyle \frac{\pi }{4}\right).\end{array}\end{eqnarray} \tag{ A1 }$$

Here we use P_k /ΔP to replace k − 1/2 in the original formula in Montgomery et al. (2003) in order to better fit the model oscillation periods to the asymptotic relation, i.e., P_k ≈ 2π²Π₀(k + k₀)/L (Tassoul 1980), where k₀ is generally not −1/2. Therefore the variation of the oscillation period is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\delta {P}_{k}}{{P}_{k}} & \approx & -\displaystyle \frac{\delta {\omega }_{k}}{{\omega }_{k}}\approx -2{{\rm{\Pi }}}_{0}{\displaystyle \int }_{0}^{1}\displaystyle \frac{\delta N}{N}\displaystyle \frac{N}{x}{dx}{\sin }^{2}\left(\pi \displaystyle \frac{{P}_{k}}{{\rm{\Delta }}P}\displaystyle \frac{{{{\rm{\Pi }}}_{x}}^{-1}}{{{{\rm{\Pi }}}_{0}}^{-1}}+\displaystyle \frac{\pi }{4}\right)\\ & = & -2{{\rm{\Pi }}}_{0}{\displaystyle \int }_{0}^{1}\displaystyle \frac{\delta N}{N}\displaystyle \frac{N}{x}{dx}\displaystyle \frac{1}{2}\left[1-\cos \left(2\pi \displaystyle \frac{{P}_{k}}{{\rm{\Delta }}P}\displaystyle \frac{{{{\rm{\Pi }}}_{x}}^{-1}}{{{{\rm{\Pi }}}_{0}}^{-1}}+\displaystyle \frac{\pi }{2}\right)\right]\\ & = & -{{\rm{\Pi }}}_{0}{\displaystyle \int }_{0}^{{{{\rm{\Pi }}}_{0}}^{-1}}\displaystyle \frac{\delta N}{N}d{{{\rm{\Pi }}}_{r}}^{-1}\left[1+\sin \left(2\pi \displaystyle \frac{{P}_{k}}{{\rm{\Delta }}P}\displaystyle \frac{{{{\rm{\Pi }}}_{x}}^{-1}}{{{{\rm{\Pi }}}_{0}}^{-1}}\right)\right]\\ & = & -{\displaystyle \int }_{0}^{1}z(u)\left[1+\sin \left(2\pi \displaystyle \frac{{P}_{k}}{{\rm{\Delta }}P}u\right)\right]{du}\\ & \approx & -{\displaystyle \int }_{0}^{1}z(u)\{1+\sin [2\pi (k+{k}_{0})u]\}{du}.\end{array}\end{eqnarray} \tag{ A2 }$$

We validate Equation (A2) in three cases: (i) a constant z = δ N/N within a region, (ii) a linear function of z within a region, and (iii) a sine function of z within a region. The reference model used for validation has a stellar mass of 2.983 solar masses, Z = 0.0317, $\mathrm{log}\theta =1.932$, $\mathrm{log}C=-4$, and τ = 0.1942. However, it should be noted that Equation (A2) is independent of specific stellar models. To quantitatively validate Equation (A2), we calculate the l = 1 oscillation periods P_k for models with different perturbations of z applied to the reference model. We then compare the quantity P_k /P_k,0 − 1 with Equation (A2), where P_k,0 represents the oscillation period of the reference model without any perturbations.

In case 1, the perturbation of N is

$$\begin{eqnarray}&&z(u)=\left\{\begin{array}{c}0,u\lt {u}_{1}\\ A,{u}_{1}\leqslant u\leqslant {u}_{2}\\ 0,u\gt {u}_{2}\end{array}\right..\end{eqnarray} \tag{ A3 }$$

Based on Equation (A2), the resulting variation of oscillation period is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\delta {P}_{k}}{{P}_{k}} & \approx & \displaystyle \frac{A}{2\pi (k+{k}_{0})}\{\cos [2\pi (k+{k}_{0}){u}_{2}]\\ & & -\,\cos [2\pi (k+{k}_{0}){u}_{1}]\}-A({u}_{2}-{u}_{1}).\end{array}\end{eqnarray} \tag{ A4 }$$

We have set u₁ = 0.25, u₂ = 0.3, and A = 0.02. The results are shown in Figure 9(a).

Zoom In Zoom Out Reset image size

In case 2, the perturbation of N is

$$\begin{eqnarray}&&z(u)=\left\{\begin{array}{c}0,u\lt {u}_{1}\\ A\displaystyle \frac{u-{u}_{1}}{{u}_{2}-{u}_{1}},{u}_{1}\leqslant u\leqslant {u}_{2}\\ 0,u\gt {u}_{2}\\ \end{array}\right..\end{eqnarray} \tag{ A5 }$$

Based on Equation (A2), the resulting variation of the oscillation period is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\delta {P}_{k}}{{P}_{k}}\approx \displaystyle \frac{A}{2\pi (k+{k}_{0})}\{\cos [2\pi (k+{k}_{0}){u}_{2}]\\ \quad -\,\displaystyle \frac{\sin [2\pi (k+{k}_{0}){u}_{2}]-\sin [2\pi (k+{k}_{0}){u}_{1}]}{2}\pi (k+{k}_{0})({u}_{2}\\ \quad -\,{u}_{1})\}-\displaystyle \frac{1}{2}A({u}_{2}-{u}_{1}).\end{array}\end{eqnarray} \tag{ A6 }$$

We have set u₁ = 0.4, u₂ = 0.45, and A = 0.02. The results are shown in Figure 9(b).

In case 3, the perturbation of N is

$$\begin{eqnarray}&&z(u)=\left\{\begin{array}{c}0,u\lt {u}_{1}\\ A\sin \left(\displaystyle \frac{u-{u}_{1}}{{u}_{2}-{u}_{1}}\pi \right),{u}_{1}\leqslant u\leqslant {u}_{2}\\ 0,u\gt {u}_{2}\\ \end{array}\right..\end{eqnarray} \tag{ A7 }$$

Based on Equation (A2), the resulting variation of oscillation period is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\delta {P}_{k}}{{P}_{k}}\approx -\displaystyle \frac{{A}_{n}({u}_{2}-{u}_{1})}{\pi }\\ \times \,\left\{2+\displaystyle \frac{\sin [2\pi (k+{k}_{0}){u}_{1}]+\sin [2\pi (k+{k}_{0}){u}_{2}]}{1-{\left[2(k+{k}_{0})({u}_{2}-{u}_{1})\right]}^{2}}\right\}.\end{array}\end{eqnarray} \tag{ A8 }$$

We have set u₁ = 0.5, u₂ = 0.55, and A = 0.02. The results are shown in Figure 9(c).

We find that Equation (A2) well reproduces the variations of the oscillation periods in all cases of perturbations of N. The differences between the numerical results and the analytical results are relatively significant in the low wavenumber region, while they are small in the high wavenumber region. This discrepancy is likely due to the limitations of the JWKB approximations and the asymptotic expression of g-mode periods, which may not hold well in the low wavenumber region.