© The Authors 2025

1. Introduction

Understanding the extent of mixing within the helium core is crucial to refining estimates of core-helium-burning (CHeB) phase lifetimes and improving predictions of stellar internal structures in subsequent evolutionary stages. This knowledge plays a key role in developing accurate models of stellar populations and elucidating the chemical evolution ofgalaxies.

Notably, observations of mixed dipole modes in CHeB stars provide a valuable tool for probing their innermost regions (e.g. Montalbán et al. 2010; Bedding et al. 2011; Mosser et al. 2012a,b; Montalbán et al. 2013; Deheuvels et al. 2016; Vrard et al. 2016; Mosser et al. 2017; Deheuvels & Belkacem 2018; Mosser et al. 2018, 2024). These oscillation modes may reveal structural variations that occur at the boundaries of mixed convective regions, in areas of element ionisation, or at interfaces between layers with different chemical compositions resulting from nuclear burning. In particular, we expect that CHeB stars exhibit chemical composition gradients at the boundary between the convective and radiative core, within the hydrogen-burning shell, and at the base of the convective envelope (e.g. Ledoux 1947; Schwarzschild 1958; Castellani et al. 1971b; Kippenhahn et al. 2012). In low-mass stars with a degenerate helium core, the transition from the red giant branch (RGB) to the CHeB phase is expected to involve a succession of off-centre helium flashes (e.g. Kippenhahn et al. 2012), which induce additional chemical composition gradients within the radiative core.

In particular, as shown in Figure 1, the temperature (∇ := dlnT/dlnP), adiabatic [∇_ad := (∂lnT/∂lnP)_S], and mean molecular weight (∇_μ := dlnμ/dlnP) gradients manifest as signatures in the squared Brunt-Väisälä frequency,

Fig. 1. Brunt-Väisälä frequencies as a function of internal mass for five different star models at the beginning of the CHeB stage (top panel). The orange, green, red, and violet lines represent models calculated using the CLES, MESA, and BaSTI-IAC evolutionary codes (see Section 3), while the blue line corresponds to the fiducial barotropic model discussed in Section 3. All models are for a 1 M⊙ CHeB star with solar metallicity. The red line represents a MESA model computed with the boundary mixing prescriptions outlined in Appendix F. These models reveal common features in their Brunt-Väisälä frequencies, including a convective core, radiative core, hydrogen-burning shell, radiative envelope, and convective envelope (the latter not shown for clarity). The bottom panel shows the helium mass fraction as a function of internal mass for four of the five models.

$$\begin{array}{cc}\hfill N^2(r)& :=\frac{Gm(r)}{r^2}[\frac{1}{{\mathrm{\Gamma }}_1}\frac{dlnP(r)}{\mathit{dr}}-\frac{dln\rho (r)}{\mathit{dr}}]\hfill \\ \hfill & \approx \frac{G^2{m}^2\rho }{r^{4}P}({\mathrm{\nabla }}_{\mathrm{ad}}-\mathrm{\nabla }+{\mathrm{\nabla }}_{\mu }),\hfill \end{array}$$(1)

which determines the behaviour of oscillation modes. This frequency can also produce deviations in period spacing (ΔP, defined in e.g. Aerts et al. 2010) from the asymptotic value when the characteristic scales of structural variations within the core are comparable to or smaller than the local wavelength of the waves under investigation (see Cunha et al. 2020, for a review). These distinct signatures on the eigenfrequencies are commonly referred to as glitch signatures.

Previous studies have examined near-core mixing conditions through the asymptotic period spacing of dipole modes ΔP_a (e.g. Bossini et al. 2015; Constantino et al. 2015; Bossini et al. 2017; Noll et al. 2024), but further insights can be gained by analysing individual eigenfrequencies. Despite extensive observational and theoretical studies on these glitches (e.g. Miglio et al. 2008; Bossini et al. 2015; Cunha et al. 2015; Mosser et al. 2015; Cunha et al. 2019; Vrard et al. 2022; Hatta 2023; Cunha et al. 2024), the effects of sharp variations in the density profile on eigenfrequencies are still not fully explored. A significant challenge in the detection of glitch signatures arises from the inherent complexity of the spectra, which exhibit a densely packed pattern of mixed modes. This complicates the isolation of individual effects, such as from magnetic fields, rotation, and strong gradients in chemical composition, and it increases the difficulty in interpreting the spectra themselves. Therefore, models are needed to make the interpretation of the observed spectra more reliable.

The seismic signature of density discontinuities associated with abrupt composition changes has already been addressed by McDermott (1990), who noted that such jumps can significantly influence the pulsation mode spectrum and stability of non-radially oscillating stars. A density discontinuity can locally contribute to or even become the dominant source of buoyancy. Such a discontinuity is defined as any region where the density changes so rapidly that the density scale height becomes much smaller than the effective wavelength of any oscillation mode of interest. The role of these density discontinuities in the spectrum of gravity modes was also explored by Gabriel & Scuflaire (1979), who found that abrupt density changes can lead to mode trapping of g-modes with the proper effective wavelength. If the mode frequency associated with a density discontinuity is comparable to or larger than the ordinary g-mode frequency, these discontinuities cannot be ignored in computing the g-modespectrum.

Structural glitches give rise to periodic components in the period spacing, allowing for the recovery of information about the location and sharpness of these glitches from their periodicity and amplitude. Therefore, deviations of the period spacing from ΔP_a contain information about sharp features in the N frequency, which can be visualised through the normalised buoyancy radius

$$\begin{array}{c}\hfill \mathrm{\Phi }(r):=\frac{{\displaystyle {\int }_{r_1}^r}N/ydy}{{\displaystyle {\int }_{r_1}^{r_2}}N/ydy},N>0,r_1\le r\le r_2.\end{array}$$(2)

This quantity is a monotonically increasing function between the inner turning point, r₁, and the outer turning point, r₂, effectively serving as a radial coordinate. The coordinate Φ(r) is particularly valuable for studying the core-envelope symmetry of high-overtone modes, though this symmetry is approximate and can be disrupted by various mechanisms (Montgomery et al. 2003). Within the JWKB approximation (developed by Jeffreys, Wentzel, Kramers and Brillouin, see e.g. Unno et al. 1989; Gough 2007), the kinetic energy density per unit of normalised buoyancy radius is generally constant, except at sharp features such as composition transition zones. This property makes variations in kinetic energy density a clear indicator of mode trapping. Moreover, within the JWKB approximation, the period of the component in the period spacing (P_glitch) associated with a structural glitch located at a radius r_glitch is related to its normalised buoyancy radius, [Φ(r_glitch)], and its periodicity in radial order, (Δn_glitch), via the relation

$$\begin{array}{c}\hfill \mathrm{\Delta }{n}_{\mathrm{glitch}}\approx 1/\mathrm{\Phi }(r_{\mathrm{glitch}})\approx P_{\mathrm{glitch}}/\mathrm{\Delta }{P}_{\mathrm{a}}\end{array}$$(3)

(see e.g. Miglio et al. 2008; Cunha et al. 2015, 2019; Vrard et al. 2022; Hatta 2023). This equation is particularly useful for interpreting the periodic signals observed in the period spacing of the oscillation modes, enabling the extraction of information about the location, amplitude, and type of structural glitches.

Consequently, glitch signatures play a significant role in advancing understanding of the complex case of mixing in CHeB stars, especially in the absence of a robust mixing theory. In this context, we face challenges related to the management of an expanding convective core and the development of semi-convective layers (e.g. Castellani et al. 1971a; Straniero et al. 2003). Furthermore, evolutionary models may generate artificial glitch signatures due to numerical inaccuracies in their predictions, complicating the interpretation of the eigenfrequencies. To address these challenges, flexible and self-consistent models are essential, particularly for studying low-mass CHeB stars with high coupling factors (e.g. Matteuzzi et al. 2023, 2024), as they offer deeper insights into their internal structures compared to other red clump (RC) stars.

This paper is organised as follows. In Section 2 we describe the theoretical framework for solving the differential equations of stellar structure. Section 3 provides details on the semi-analytical model of a 1 M_⊙ RC star with solar composition and Y_c ≈ 0.9, which serves as a reference model for subsequent modifications. Section 4 presents an explanation of how we simulated different levels of mixing between adjacent zones in aself-consistent manner. The results are presented in Section 5, and Section 6 concludes the paper.

2. Parametric modelling of CHeB stars based on realistic stellar structures

The structure of a single, spherical, non-rotating, non-magnetic star can be described by the hydrostatic and mass continuity equations (see e.g. Kippenhahn et al. 2012). To fully solve these differential equations, it is essential to establish a relation between pressure and density. Furthermore, this relation must include temperature, mass fractions of chemical species, and their corresponding chemical potentials to formulate a complete system of differential equations. In particular, we must incorporate an equation of state (EOS) that connects pressure with these other variables, thereby resulting in a closed system of equations. In this paper, we employ a differential representation of an EOS. We started by defining the variable

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{local}}(r):=\frac{dlnP(r)}{dln\rho },\end{array}$$(4)

which allowed us to derive the differential equation

$$\begin{array}{c}\hfill \frac{dP(r)}{\mathit{dr}}={\gamma }_{\mathrm{local}}(r)\frac{P(r)}{\rho (r)}\frac{d\rho (r)}{\mathit{dr}}\end{array}$$(5)

that has to be solved together with the hydrostatic and mass continuity equations. All the information provided by the EOS is encapsulated in equation (4), which can be seen through a heuristic approach. For simplicity, and without loss of generality, we may consider γ_local(r)≡γ_local(T, ρ, μ). By applying the chain rule to equation (4), we are able to demonstrate that

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{local}}(T,\rho ,\mu )=\frac{{\chi }_{\rho }}{1-{\chi }_T\mathrm{\nabla }-{\chi }_{\mu }{\mathrm{\nabla }}_{\mu }},\\ \hfill \end{array}$$(6)

where ${\chi }_{\rho }:={\frac{\partial lnP}{\partial ln\rho }|}_{T,\mu }$, ${\chi }_{\mu }:={\frac{\partial lnP}{\partial ln\mu }|}_{\rho ,T}$, ${\chi }_T:={\frac{\partial lnP}{\partial lnT}|}_{\rho ,\mu }$. At this point it becomes more evident that γ_local is predominantly influenced by the thermal and chemical gradients present within the star and that in adiabatically stratified layers, γ_local ≡ Γ₁. Moreover, according to equation (5), a polytropic EOS can be derived (i.e. P ∝ ρ^γ_local) when γ_local is treated as a constant. For the purposes of the following discussion, we denote γ_local(r)≡γ(r).

2.1. Barotropic CHeB stars

As demonstrated in Section 2, a closed set of equations is obtained when we provide the temperature profile and the chemical species. However, stars can be accurately modelled as systems in barotropic equilibrium (i.e. the thermodynamic variables are stratified as the density profile), because $\mathbf{\nabla }\rho \times \mathbf{\nabla }P=\mathbf{0}$ everywhere in self-gravitating spherical stars in hydrostatic equilibrium. Furthermore, the density is a monotonic function of the radius with the exception of a small region near the surface where a density inversion due to hydrogen recombination occurs (e.g. Chitre & Shaviv 1967; Érgma 1971; Harpaz 1984). To provide a more quantitative estimate of the deviation from monotonicity, we analyse the density inversion region in a representative stellar model using the code CLES (see Section 3). Our findings indicate that a non-injective region in the function ρ(P) appears only very close to the surface (at ≈99.8% of the total stellar radius) with an extent that constitutes only 0.14% of the total stellar radius and about 10⁻⁵% of the total stellar mass. Therefore, despite the relatively high values of dP/P ≈ 44% and dρ/ρ ≈ 5% in this region, its very limited size combined with a very low mass content means that it has a negligible effect on the global structure of the star. This observation leads to the conclusion that, with a high degree of accuracy, P(r)≡P(ρ), T(r)≡T(ρ), X_i(r)≡X_i(ρ) and γ(r)≡γ(ρ). These relations can be incorporated into equation 5 to derive the differential equation

$$\begin{array}{c}\hfill \frac{dP(r)}{\mathit{dr}}=\frac{\gamma (\rho )P(\rho )}{\rho (r)}\frac{d\rho (r)}{\mathit{dr}}.\end{array}$$(7)

A general solution of equation (7) is

$$\begin{array}{c}\hfill P(\rho )=P_{c}exp\left({\displaystyle {\int }_{ln{\rho }_c}^{ln\rho }}\gamma (t)dt\right),\end{array}$$(8)

where P_c and ρ_c are the central pressure and density, respectively. Therefore, an a priori knowledge of γ(ρ), from calibrations on evolutionary models for example, allows us to find the internal structure of the star using the hydrostatic equation, the mass continuity equation, and equation (7) or (8).

2.1.1. Differential equations for barotropic stars

We normalised all the quantities with respect to reference density (ρ₁) and pressure (P₁) values. As a result, once the dimensionless quantities

$$\begin{array}{c}\hfill \{\begin{array}{c}\theta :=\frac{\rho }{{\rho }_1}\hfill \\ \xi :=r\sqrt{\frac{4\pi G{\rho }_1^2}{P_1}}\hfill \\ \beta :=\frac{P}{P_1}\hfill \\ \psi :=m\sqrt{\frac{4\pi G^3{\rho }_1^4}{P_1^3}}\hfill \end{array}\end{array}$$(9)

were defined, the hydrostatic equation, the mass continuity equation, equation (7), and equation (8) became

$$\begin{array}{c}\hfill \{\begin{array}{c}\frac{d\theta (\xi )}{d\xi }=-\frac{\psi (\xi )\theta {(\xi )}^2}{\gamma (\theta )\beta (\theta ){\xi }^2}\hfill \\ \frac{d\psi (\xi )}{d\xi }={\xi }^2\theta (\xi )\hfill \\ \beta (\theta )={\beta }_{c}exp\left({\displaystyle {\int }_{ln{\theta }_c}^{ln\theta }}\gamma (t)dt\right),\hfill \end{array}\end{array}$$(10)

with initial conditions θ(ξ = 0) = θ_c, β(ξ = 0) = β_c and ψ(ξ = 0) = 0. We numerically solved equation (10) with a Dormand and Prince Runge-Kutta fifth-order method from the centre to the surface and obtained the θ(ξ) and ψ(ξ) profiles. Finally, with such profiles, we a posteriori calculated other variables of interest, such as the Brunt-Väisälä frequency

$$\begin{array}{c}\hfill N^2(\xi )=4\pi G{\rho }_1\frac{\psi {(\xi )}^2\theta (\xi )}{\beta (\theta ){\xi }^4}[\frac{1}{\gamma (\theta )}-\frac{1}{{\mathrm{\Gamma }}_1}].\end{array}$$(11)

To mitigate numerical issues when solving equation (10) from the centre, we employed Taylor expansions to approximate the actual solutions near the centre. More details regarding these expansions are provided in Appendix A. We also note that in certain situations it may be more beneficial to solve alternative differential equations rather than equation (10). In Appendix B, we examine some of these alternative equations, with a specific focus on equations near the stellar surface. In Appendix C, we present analytical solutions to equation 10 for constant values of γ(r), which we utilised to validate our numerical solver.

2.2. Discontinuities in the internal profiles

In CHeB stars we expect that certain internal profiles (e.g. ρ) exhibit discontinuities or non-differentiable points. For example, we expect a sharp decrease in density at the boundary between convective and radiative core due to an abrupt change in the chemical profile (see Section 1). Therefore, we can establish the jump conditions at a given radial position r₀ by ensuring the continuity of pressure and mass. In fact, in CHeB stars at equilibrium, there are no internal shock waves, and the locations where the internal profiles show discontinuities form a countable set with zero Lebesgue’s measure. This ensures that the integrals of these internal profiles are continuous. For simplicity let us use the notation lim_{r → r0⁻}f(r):=f⁻ and lim_{r → r0⁺}f(r):=f⁺ for an internal profile f(r). Such a profile is continuous when f⁻ = f⁺ = f(r₀). We defined

$$\begin{array}{c}\hfill \mathrm{\Lambda }:=\frac{{\rho }^{+}}{{\rho }^{-}}\in {\mathbb{R}}^{+}\end{array}$$(12)

because we are mainly interested in jump discontinuities in the density profile. The jump conditions are then

$$\begin{array}{c}\hfill \{\begin{array}{c}{\xi }^{+}={\xi }^{-}\hfill \\ {\beta }^{+}={\beta }^{-}\hfill \\ {\psi }^{+}={\psi }^{-}\hfill \\ {\theta }^{+}={\theta }^{-}\mathrm{\Lambda }\hfill \\ {\frac{d\theta }{d\xi }|}_{r_0^{+}}={\frac{d\theta }{d\xi }|}_{r_0^{-}}{\mathrm{\Lambda }}^2\left(\frac{{\gamma }^{-}}{{\gamma }^{+}}\right)\hfill \\ {\frac{d\psi }{d\xi }|}_{r_0^{+}}={\frac{d\psi }{d\xi }|}_{r_0^{-}}\mathrm{\Lambda }\hfill \\ {(N^2)}^{+}={(N^2)}^{-}\mathrm{\Lambda }\left(\frac{{\mathrm{\Gamma }}_1^{+}-{\gamma }^{+}}{{\mathrm{\Gamma }}_1^{-}-{\gamma }^{-}}\right)\left(\frac{{\gamma }^{-}}{{\gamma }^{+}}\right)\left(\frac{{\mathrm{\Gamma }}_1^{-}}{{\mathrm{\Gamma }}_1^{+}}\right).\hfill \end{array}\end{array}$$(13)

Next, we introduce r₀ as the boundary between two zones with different normalisation constants, P₁ and ρ₁. This choice may be useful when equation 10 creates numerical issues (for example with small θ and β values; see Appendix B). We then defined as normalisation constants for the first zone (i.e. for r < r₀ ) the density and the pressure in the centre of the star (i.e. P₁ ≡ P_c and ρ₁ ≡ ρ_c). This simplifies all the equations because now θ_c ≡ 1 and β_c ≡ 1. In the second zone (i.e. where r > r₀), it is useful to define as normalisation constants the pressure and the density at the boundary between the two zones. Therefore, we have P₁ ≡ P(r₀) = P_cexp(∫₀^lnθ−γ(t) dt) and ρ₁ ≡ Λρ⁻. The new jump conditions are then

$$\begin{array}{c}\hfill \{\begin{array}{c}{\xi }^{+}={\xi }^{-}\frac{\mathrm{\Lambda }{\theta }^{-}}{\sqrt{{\beta }^{-}}}\hfill \\ {\beta }^{+}=1\hfill \\ {\psi }^{+}={\psi }^{-}\frac{{(\mathrm{\Lambda }{\theta }^{-})}^2}{{({\beta }^{-})}^{\frac{3}{2}}}\hfill \\ {\theta }^{+}=1\hfill \\ {\frac{d\theta }{d\xi }|}_{r_0^{+}}={\frac{d\theta }{d\xi }|}_{r_0^{-}}\frac{{\gamma }^{-}}{{\gamma }^{+}}\frac{\sqrt{{\beta }^{-}}}{{({\theta }^{-})}^2}\hfill \\ {\varphi }^{+}=1\hfill \\ {\frac{d\varphi }{d\xi }|}_{r_0^{+}}=({\gamma }^{+}-1){\frac{d\theta }{d\xi }|}_{r_0^{-}}\frac{{\gamma }^{-}}{{\gamma }^{+}}\frac{\sqrt{{\beta }^{-}}}{{({\theta }^{-})}^2},\hfill \end{array}\end{array}$$(14)

where we include the ϕ and $\frac{d\varphi }{d\xi }$ variables defined in Appendix B because they are useful for Section 3.3.

3. Fiducial barotropic model

In this section, we derive a fiducial barotropic model for a 1 M_⊙ star with solar composition at the onset of the CHeB stage, ensuring it incorporates all the fundamental characteristics of this evolutionary phase. As discussed in Section 2.1, knowing γ(ρ) allows us to find the internal structure of the star using the hydrostatic equation, the mass continuity equation, and equation (7) or (8). However, this information is insufficient on its own; we need detailed evolutionary models for the calibration of γ(ρ). For this purpose we employ the evolutionary codes CLES v21.0 (Code Liégeois d’Évolution Stellaire, Scuflaire et al. 2008b) and MESA v11701 (Modules for Experiments in Stellar Astrophysics; Paxton et al. 2011, 2019). In both codes we focused on models that have a mass fraction of helium in the centre of about 0.9 (i.e. Y_c ≈ 0.9). In particular, both in CLES v21.0 and in MESA v11701 we adopt as the reference solar mixture that from Asplund et al. (2009), and high- and low-temperature radiative opacity tables are computed for the solar specific metal mixture. The envelope convection is described by the mixing length theory of Cox & Giuli (1968); the corresponding α_MLT parameter, the same for all the models, is derived from the solar calibration with the same physics. Below the convective envelope, we add a diffusive undershooting (Herwig 2000) with a size parameter f = 0.02 (see Khan et al. 2018) in the MESA v11701 models, while the CLES v21.0 models incorporated a step undershooting in the form of penetrative convection. Furthermore, additional mixing over the convective core limit during the CHeB phase is treated following the formalism by Bossini et al. (2017) in the MESA v11701 models. In contrast, in the CLES v21.0 models we adopt a step overshooting of 0.50H_P in the form of penetrative convection, along with prescriptions designed to establish a semi-convective region where ∇_rad = ∇_ad (Castellani et al. 1971a). Unlike MESA v11701, CLES v21.0 begins its simulations from a zero-age horizontal branch, employing a fixed initial helium core mass of 0.45 M_⊙. As a result, these models do not undergo helium flashes, but approximately 5% by mass of helium in the core is converted into carbon (more details in Panier et al. 2025). Lastly, we also simulate a diffusive convective boundary mixing process with MESA, as detailed in Appendix F.

In Figure 2 we show γ(r) of five different models as a function of density. Three of them are computed with CLES v21.0 and MESA v11701, while another one is a model with Z = 0.012580 taken from the BaSTI-IAC library¹ explained in Hidalgo et al. (2018). It is clear that they share the same fundamental characteristics of the CHeB phase. We expect a similar situation for models of different masses and/or metallicity, provided that they retain equivalent internal structures. Specifically, these models consist of a convective core surrounded by helium-rich material, an hydrogen-burning shell at the end of the helium core and a convective envelope. The main difference being the location (in density) of such characteristics. Finally, the last γ(r) in the figure is a reference barotropic model we derive in Sections 3.1, 3.2 and 3.3, which we call the fiducial barotropic model. This model is calibrated on the evolutionary model computed with the code CLES v21.0, which we simply call the reference CLES model. We decided to divide γ(r) of the barotropic model in three different zones: the convective core (Section 3.1), the inner radiative region of the star (i.e. the radiative core, the hydrogen-burning shell and the radiative envelope, see Section 3.2), and the convective envelope (Section 3.3).

It is essential to emphasise that the primary aim of this semi-analytical modelling is to explore the details of the gravity modes spectrum. Therefore, our analysis concentrates on the characteristics of the stellar radiative regions, and we simplified the modelling of the convective envelope, as it is less relevant to our investigation.

3.1. Convective core

The first zone is the convective core of the CHeB star. We model this part by using the same ρ_c and P_c as the reference CLES model, but a different Γ₁. We assume a constant Γ₁ equal to γ(θ), i.e. we assume a polytropic EOS in this region. We calibrate the barotropic model to have the same radius and mass at the end of the convective core as the reference CLES model by keeping Γ₁ as a free parameter to fit. The results of the calibration are in Table 1, which clearly indicates that maintaining a fixed Γ₁ prevents the accurate representation of pressure and density at the end of the convective core.

3.2. Radiative region

In the second zone, we simulated the change in chemical composition at the boundary between convective and radiative core with a jump in density identical to the reference CLES model (i.e. we adopted the same Λ < 1). This approach is justified by the fact that in the perfect gas approximation, Λ ≈ μ⁺/μ⁻ because the pressure and the temperature are continuous functions at the boundary between the convective and radiative core. In this second zone, we modelled the star from the base of the radiative core to the base of the convective envelope. In particular, we took

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathrm{\Gamma }}_1=\frac{5}{3}\hfill \\ 0\le \gamma (\rho )<{\mathrm{\Gamma }}_1\hfill \\ \theta :=\frac{\rho }{\mathrm{\Lambda }{\rho }^{-}}\hfill \\ \beta :=\frac{P}{P_0}=exp\left({\displaystyle {\int }_0^{ln\theta }}\gamma (t)dt\right)\hfill \end{array}\end{array}$$(15)

with

$$\begin{array}{c}\hfill \gamma (u)=\{\begin{array}{c}Aexp[-\frac{{(u-u_P)}^2}{2{\sigma }^2}]+aexp\left(Ku\right),u_P\le u\le 0\hfill \\ \hfill \\ [A+aexp\left(K{u}_P\right)-(a^{\prime })exp\left(K^{\prime }{u}_P\right)]exp[-\frac{{(u-u_P)}^2}{2{({\sigma }^{\prime })}^2}]+\hfill \\ +(a^{\prime })exp\left(K^{\prime }u\right),u<u_P,\hfill \end{array}\end{array}$$(16)

where we define u := lnθ and u_P := lnθ_P as useful variables to simplify subsequent calculations. Therefore, we modelled γ(θ) in this zone with eight parameters: [A, θ_P, σ, a, K, σ′,K′,a′]. The normalised pressure was found when applying equation (16) to β(θ) of equation 15, that is

$$\begin{array}{c}\hfill ln\beta (\theta )=\{\begin{array}{c}A\sqrt{\frac{\pi }{2}}\sigma [\mathrm{erf}\left(\frac{ln\theta -ln{\theta }_P}{\sqrt{2}\sigma }\right)+\mathrm{erf}\left(\frac{ln{\theta }_P}{\sqrt{2}\sigma }\right)]+\hfill \\ +\frac{a}{K}({\theta }^K-1),{\theta }_P\le \theta \le 1\hfill \\ \hfill \\ A\sqrt{\frac{\pi }{2}}\sigma \left[\mathrm{erf}\left(\frac{ln{\theta }_P}{\sqrt{2}\sigma }\right)\right]+\frac{a}{K}({\theta }_P^K-1)+\hfill \\ +[A+a{\theta }_P^K-(a^{\prime }){\theta }_P^{K^{\prime }}]\sqrt{\frac{\pi }{2}}({\sigma }^{\prime })\mathrm{erf}\left(\frac{ln\theta -ln{\theta }_P}{\sqrt{2}{\sigma }^{\prime }}\right)+\hfill \\ +\frac{a^{\prime }}{K^{\prime }}({\theta }^{K^{\prime }}-{\theta }_P^{K^{\prime }}),\theta <{\theta }_P.\hfill \end{array}\end{array}$$(17)

We found initial guesses of the eight parameters by performing a fit between equation (16) and the γ(θ) of the reference CLES model with a maximum likelihood estimator. Subsequently, we modify the values of K′ and a′ to obtain the same radius and mass at the base of the convective envelope (BCE) as in the reference CLES model (see Table 1). This table further illustrates that the modelling is sufficiently accurate to yield an accurate density for the BCE, despite variations in pressure. However, in the final fiducial model we do not keep these values of K′ and a′, because they would lead to unphysical values of the total stellar mass and radius (Section 3.3).

Finally, we note that in the reference CLES model there is another drop in density located at the end of the helium core as a signature of the helium-flashes (i.e. the drop in γ of Figure 2 near ρ = 100 g/cm³). However, in our fiducial barotropic model we excluded such a signature because we explore better additional jump discontinuities in Section 5.

Fig. 2. Comparison of γ(r) as a function of density for the five different models shown in Figure 1. The dashed brown line represents the Γ1 used by CLES v21.0 for reference. All five models share the fundamental characteristics of the CHeB phase.

3.3. Convective envelope

In the reference CLES model, the density is non-differentiable at the BCE due to the implemented adiabatic step undershooting explained in Section 3. This means that γ has (with abuse of notation) a ‘jump discontinuity’ at the BCE (see the ‘jump’ in γ as illustrated in Figure 2 at ρ = 0.01 g/cm³). This characteristic is retained in our barotropic model. Moreover, in contrast to the convective core, the convective envelope (designated as the third zone in our barotropic model) exhibits low efficiency in convective transport in the near-surface layers, primarily due to their low density values. As a result, these near-surface layers are characterised by a temperature gradient that exceeds the adiabatic temperature gradient, which can be modelled with the mixing length theory (e.g. Cox & Giuli 1968). However, in our barotropic model we adopt a γ(ρ)≡Γ₁ (i.e. an adiabatic temperature gradient) everywhere in the convective envelope to simplify subsequent calculations. Additionally, we do not account for helium and hydrogen recombination and we keep a constant Γ₁ = 5/3 (similarly to the second zone, see Section 3.2). This simplification implies that the semi-analytical models we compute cannot be used to measure acoustic glitch signatures such as the helium and hydrogen partial ionisation regions. Nevertheless, our models show the expected pattern of p-mode frequencies, characterised by modes that are nearly equally spaced in frequency (see Section 5.1). Finally, we determine the total stellar mass and radius by solving the Lane-Emden equation (B.8), with the appropriate jump conditions described in equation (14), until we reach a density ρ = 0.

Due to the simplifications implemented, particularly within the convective envelope, our barotropic model does not yield the same total mass and radius as the reference CLES model; instead, it results in a star that is larger and more massive. A correct total radius or a correct total mass can be achieved when Γ₁ > 5/3 within the convective envelope. However, we decide to retain a more physically plausible value of Γ₁ = 5/3 and adjust the parameters [K′,a′]² to ensure that the barotropic model attains the correct total mass and radius. The change in radius, mass, density and pressure at the BCE resulting from the introduction of the new [K′,a′] parameters are presented in Table 1.

4. Smooth transitions in discontinuous density profiles

In real stars, both macroscopic (see Appendix F) and microscopic (e.g. Michaud et al. 1984, 2010) diffusion processes are expected to smooth gradients in pressure, temperature, concentration and density. They tend to smooth out such gradients when present, as in the case of jump discontinuities in density. Therefore, to make our barotropic model even more realistic, we need a function that smoothly joins two adjacent zones having otherwise a jump discontinuity in density. This must be done in a self-consistent manner and by simulating different levels of mixing between the two zones.

In the fiducial barotropic model of Section 3, we have a jump discontinuity in the density at the boundary between convective and radiative core. We can smoothly join these two zones using a sigmoid function S(P) in the lnP − lnρ plane. This means that ρ ≡ ρ(ρ_I, ρ_R, S), where ρ_I(r) is the density profile in the convective core and ρ_R(r) the density profile in the radiative part of the star. A convenient choice for S(P) is that it vanishes in the radiative part and becomes one in the convective core. Consequently, one can define ρ(ρ_I, ρ_R, S)≡ρ(P) as

$$\begin{array}{c}\hfill ln\rho (P):=ln{\rho }_R+(ln{\rho }_I-ln{\rho }_R)S(P),\end{array}$$(18)

the derivative of which is

$$\begin{array}{c}\hfill \frac{dln\rho }{dlnP}=\frac{1-S(P)}{{\gamma }_R[{\rho }_R(P)]}+\frac{S(P)}{{\mathrm{\Gamma }}_{1,c}}+(ln{\rho }_I-ln{\rho }_R)\frac{dS(P)}{dlnP}.\end{array}$$(19)

Finally we obtained that

$$\begin{array}{cc}\hfill \frac{1}{\gamma (\rho )}& =\frac{1-S[P(\rho )]}{{\gamma }_R\{{\rho }_R[P(\rho )]\}}+\frac{S[P(\rho )]}{{\mathrm{\Gamma }}_{1,c}}+\hfill \\ \hfill & +\{ln{\rho }_I[P(\rho )]-ln{\rho }_R[P(\rho )]\}\frac{dS[P(\rho )]}{dlnP},\hfill \end{array}$$(20)

with P(ρ) obtained by inverting equation (18) numerically. Therefore, the new stellar structure is obtained by solving equation (10) with the new P(ρ) and γ(ρ) functions obtained from equations (18) and (20), respectively. As expected, lim_S → 0γ = γ_R, lim_S → 1γ = Γ_1, c, and lim_{S → 1/2}γ has a local minimum whose value depends on the functional form of S(P).

From this point on, we decided to use the error function as a functional form of S(P):

$$\begin{array}{c}\hfill \{\begin{array}{c}S(P):=\frac{1+\mathrm{erf}[\alpha (lnP-ln{P}_0)]}{2}\hfill \\ \frac{dS(P)}{dlnP}=\frac{\alpha }{\sqrt{\pi }}exp[-{\alpha }^2{(lnP-ln{P}_0)}^2],\hfill \end{array}\end{array}$$(21)

where α ≥ 0 is a constant whose increasing value makes the slope of S(P) steeper, and P₀ is the pressure at the boundary between convective and radiative core in the fiducial model of Section 3. With such choices, the higher is α, the lower is the amount of mixing between the two consecutive zones. Indeed, S(P) in equation (21) tends to a Heaviside function for α → ∞; thus, we tend to the fiducial model of Section 3, while for α = 0 we have a fixed S(P) = 1/2 that corresponds to the maximum mixing. Furthermore, this choice for P₀ resulted in ρ_I(P₀)≡ρ⁻, ρ_R(P₀)≡Λρ⁻ and

$$\begin{array}{c}\hfill {\displaystyle \underset{P\to P_0}{lim}}\frac{1}{\gamma }=\frac{1}{2{\gamma }_R(\mathrm{\Lambda }{\rho }^{-})}+\frac{1}{2{\mathrm{\Gamma }}_{1,c}}+\frac{\alpha }{\sqrt{\pi }}ln\left(\frac{1}{\mathrm{\Lambda }}\right).\end{array}$$(22)

In summary, in the limit for α → ∞, equation (20) tends to

$$\begin{array}{c}\hfill {\displaystyle \underset{\alpha \to \infty }{lim}}\gamma (\rho )=\{\begin{array}{c}{\mathrm{\Gamma }}_{1,c},P>P_0\hfill \\ 0,P=P_0\hfill \\ {\gamma }_R(\rho ),P<P_0,\hfill \end{array}\end{array}$$(23)

and equation (18) tends to

$$\begin{array}{c}\hfill {\displaystyle \underset{\alpha \to \infty }{lim}}\rho (P)=\{\begin{array}{c}{\rho }_I(P),P>P_0\hfill \\ {\rho }_R(P),P<P_0,\hfill \end{array}\end{array}$$(24)

which means that we tend to the fiducial model of Section 3, as expected.

To obtain the inverse of ρ(P), we needed to evaluate lnρ_I(P) and lnρ_R(P). From Section 3.1 we obtained

$$\begin{array}{c}\hfill ln{\rho }_I(P)=ln{\rho }_c+\frac{lnP-ln{P}_c}{{\mathrm{\Gamma }}_{1,c}},\end{array}$$(25)

but lnρ_R(P) does not have an analytical solution. However,

$$\begin{array}{c}\hfill ln{P}_R(\rho )\approx ln{P}_0+\frac{a}{K}\frac{{\rho }^K-{\rho }_{+}^K}{{\rho }_{+}^K}\end{array}$$(26)

is a very good approximation of equation (17) when ρ ≳ 0.24Λρ⁻, because in that range their relative difference is |1 − P_true/P_approx|< 10⁻¹⁰%. Therefore, the inverse is well approximated by

$$\begin{array}{c}\hfill ln{\rho }_R(P)\approx ln(\mathrm{\Lambda }{\rho }^{-})+\frac{ln[1+\frac{K}{a}ln\left(\frac{P}{P_0}\right)]}{K}\end{array}$$(27)

when ρ_R ≳ 0.24Λρ⁻. Finally we obtained that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{1}{\gamma (\rho )}& \approx \frac{1-S[P(\rho )]}{a+K[lnP(\rho )-ln{P}_0]}+\frac{S[P(\rho )]}{{\mathrm{\Gamma }}_{1,c}}\hfill \\ \hfill & +ln\left\{\frac{{\rho }_c}{\mathrm{\Lambda }{\rho }^{-}}{\left[\frac{P(\rho )}{P_c}\right]}^{\frac{1}{{\mathrm{\Gamma }}_{1,c}}}{\left[\frac{a}{a+K(lnP-ln{P}_0)}\right]}^{\frac{1}{K}}\right\}\frac{dS[P(\rho )]}{dlnP}\hfill \end{array}\end{array}$$(28)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill ln\rho (P)& \approx ln(\mathrm{\Lambda }{\rho }^{-})+\frac{ln[1+\frac{K}{a}ln\left(\frac{P}{P_0}\right)]}{K}\hfill \\ \hfill & +ln\left\{\frac{{\rho }_c}{\mathrm{\Lambda }{\rho }^{-}}{\left(\frac{P}{P_c}\right)}^{\frac{1}{{\mathrm{\Gamma }}_{1,c}}}{\left[\frac{a}{a+K(lnP-ln{P}_0)}\right]}^{\frac{1}{K}}\right\}S(P)\hfill \\ \hfill \end{array}\end{array}$$(29)

when ρ_R ≳ 0.24Λρ⁻, and from a numerical evaluation of P(ρ) we obtained γ(ρ).

For |α(lnP − lnP₀)| > > 1, we have γ ≈ Γ_1, c for ρ ≈ ρ⁻ and γ ≈ γ_R for ρ ≈ ρ⁺. We note that to a good approximation

$$\begin{array}{c}\hfill \begin{array}{cc}& P(\rho =0.24\mathrm{\Lambda }{\rho }^{-})\approx P_R({\rho }_R=0.24\mathrm{\Lambda }{\rho }^{-})\mathrm{for}\alpha \ge 1.885\hfill \\ \hfill & \gamma (\rho =0.24\mathrm{\Lambda }{\rho }^{-})\approx {\gamma }_R({\rho }_R=0.24\mathrm{\Lambda }{\rho }^{-})\mathrm{for}\alpha \ge 2.074\hfill \end{array}\end{array}$$(30)

since the ratio of the two sides of the equation deviates from unity by less than 10⁻⁶%. This means that the approximation in equation (27), which is valid when ρ_R ≳ 0.24Λρ⁻, is also valid when ρ ≳ 0.24Λρ⁻ if α ≥ 2.074. Moreover, a consequence of equation (30) is that γ(ρ ≤ 0.24Λρ⁻)≈γ_R(ρ_R ≤ 0.24Λρ⁻) and P(ρ ≤ 0.24Λρ⁻)≈P_R(ρ_R ≤ 0.24Λρ⁻) if α ≥ 2.074. Therefore, we decided to solve equations (28) and (29) when ρ ≥ 0.24Λρ⁻, and when ρ < 0.24Λρ⁻ (i.e. until the end of the second zone), we solved γ ≈ γ_R using equation (16) and

$$\begin{array}{c}\hfill \begin{array}{cc}& lnP(\theta )\approx ln{P}_R=ln{P}_1+{\displaystyle {\int }_{ln{\theta }_1}^{ln\theta }}{\gamma }_R(t)dt\hfill \\ \hfill & =ln{P}_1+\{\begin{array}{c}A\sqrt{\frac{\pi }{2}}\sigma [\mathrm{erf}\left(\frac{ln\theta -ln{\theta }_P}{\sqrt{2}\sigma }\right)+\mathrm{erf}\left(\frac{ln{\theta }_P-ln{\theta }_1}{\sqrt{2}\sigma }\right)]\hfill \\ +\frac{a}{K}({\theta }^K-{\theta }_1^K),{\theta }_P\le \theta <{\theta }_1\hfill \\ \hfill \\ A\sqrt{\frac{\pi }{2}}\sigma \left[\mathrm{erf}\left(\frac{ln{\theta }_P-ln{\theta }_1}{\sqrt{2}\sigma }\right)\right]+\frac{a}{K}({\theta }_P^K-{\theta }_1^K)\hfill \\ +[A+a{\theta }_P^K-(a^{\prime }){\theta }_P^{K^{\prime }}]\sqrt{\frac{\pi }{2}}({\sigma }^{\prime })\mathrm{erf}\left(\frac{ln\theta -ln{\theta }_P}{\sqrt{2}{\sigma }^{\prime }}\right)\hfill \\ +\frac{a^{\prime }}{K^{\prime }}({\theta }^{K^{\prime }}-{\theta }_P^{K^{\prime }}),\theta <{\theta }_P.\hfill \end{array}\hfill \end{array}\end{array}$$(31)

Here, θ := ρ/(Λρ⁻), θ₁ := 0.24 and lnP₁ := lnP(θ₁), which is found using equations (28) and (29). Finally, we chose Γ₁ such that

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1=\{\begin{array}{c}{\mathrm{\Gamma }}_{1,c}\mathrm{for}\rho \ge {\rho }^{-}\hfill \\ (\frac{5}{3}-{\mathrm{\Gamma }}_{1,c})\frac{ln\left(\rho \right)-ln(\mathrm{\Lambda }{\rho }^{-})}{ln(\mathrm{\Lambda })}+\frac{5}{3}\mathrm{for}\mathrm{\Lambda }{\rho }^{-}<\rho <{\rho }^{-}\hfill \\ \frac{5}{3}\mathrm{for}\rho \le \mathrm{\Lambda }{\rho }^{-}.\hfill \end{array}\end{array}$$(32)

With this formalism we can explore different behaviours at the boundary between convective and radiative core by varying the parameters Λ and α. In the top panel of Figure 3, we illustrate three scenarios: one with a jump discontinuity in density (blue lines), and two models that share the same Λ as the blue case but differ in α values (orange and green lines). The bottom panel displays the corresponding Brunt-Väisälä frequencies as functions of internal mass, highlighting that a higher α leads to a sharper bell-shaped structure. Furthermore, the same formalism can be applied to smoothly connect two adjacent zones that would otherwise exhibit a jump discontinuity in density within the radiative core. In this context, the primary distinction lies in the fact that γ_I(ρ)≠Γ_1, c, and the validity of the approximation presented in the equations must be verified.

Fig. 3. Comparison between three different density profiles (top panel) at the boundary between the convective and the radiative core and their corresponding Brunt-Väisälä frequencies (bottom panel) as functions of internal mass. In particular, we compare the fiducial barotropic model (blue) with two models that smoothly join the two zones (orange and green). The orange line incorporates a lower α value than the green line. The blue arrow corresponds to a δ-distribution. The FWHM of the peaks in N presented here are significantly lower than 0.1Hp.

5. Results

In this section we explain the main results concerning our fiducial barotropic model (Section 5.1) in comparison with a different boundary of the convective core (Section 5.2), different glitches in the radiative core (Section 5.3), and different degree of jump discontinuities near the boundary of the convective core (Section 5.4). It is important to note that some of these glitches are primarily used as test cases, as within the observable region of the oscillation spectrum they would appear too smooth to be identified as glitch signatures.

We compute the adiabatic eigenfrequencies, normalised inertia (E_norm, see the definition in e.g. Aerts et al. 2010) and eigenfunctions of radial (ℓ = 0) and non-radial (ℓ = 1 − 3) modes using the code GYRE (version 6.0.1, Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020), but we verify independently the calculations with the tool LOSC (Scuflaire et al. 2008a) in the case of the fiducial model. We also compute the uncoupled g-modes, called γ-modes, through the prescriptions provided by Ong & Basu (2020) and included in GYRE. As becomes clear in the subsequent sections, this simplifies the detection of buoyancy glitches. Finally, we simulate four-year-long Kepler observations (lightcurves and power spectral densities) with the code AADG3 (AsteroFLAG Artificial Dataset Generator, version 3.0.2; Ball et al. 2018, and references therein), with information on mode lifetimes taken from the work of Vrard et al. (2018) on CHeB stars in the Kepler field.

5.1. Fiducial model

In this section, we highlight the results of the fiducial barotropic model of Section 3.

In Figure 4 we show the Brunt-Väisälä frequency as a function of the internal mass, with a specification for the position of the convective core, the hydrogen-burning shell and the convective envelope (top panel). In the boundary between the convective and radiative core there is a change in Γ₁, a change in γ(ρ) and a jump discontinuity in density. This implies that N is a δ-distribution (blue arrow in the figure) at the position corresponding to the boundary. In the BCE, the density profile is non-differentiable due to a ‘jump’ in γ(ρ), and thus N has a jump discontinuity here.

Fig. 4. Brunt-Väisälä frequency as a function of internal mass (top left), period spacing as a function of eigenfrequencies in the observable region of the oscillation spectrum (top right), and the corresponding simulated PSD (bottom) for the fiducial model discussed in Section 5.1. The cyan dashed line indicates the νmax of the reference CLES model. Along with the simulated PSD, we also show the normalised inertia of radial and non-radial modes (coloured dots connected by coloured lines).

As discussed at the beginning of Section 5, we assess our model by computing adiabatic eigenfrequencies that serve to simulate four-year-long Kepler observations. The corresponding power spectral density (PSD) as a function of frequency is shown in the bottom panel of Figure 4. This panel includes a cyan dashed line representing the ν_max, derived from the reference CLES model, alongside the normalised inertiae as a function of eigenfrequencies for modes with ℓ = 0, 1, 2,  and 3. As expected (e.g. Montalbán et al. 2010; Mosser et al. 2011, 2012b), the radial modes, along with the minima of inertia for the non-radial modes, are evenly spaced in frequency with a ⟨Δν⟩≈5μHz. Additionally, a ‘forest’ of dipole mixed modes appears between two consecutive radial modes, a characteristic observed in actual low-mass CHeB stars. We further expect that, in these stars, the period spacing (ΔP) of the observable dipole mixed modes approaches the asymptotic value (ΔP_a), which is notably higher than that of RGB stars (e.g. Montalbán et al. 2010; Bedding et al. 2011). Our fiducial model confirms this expectation, as illustrated in the upper right panel of Figure 4. It is apparent that as the frequency decreases, ΔP approaches the asymptotic value of ΔP_a = 248.34 s due to the increasing radial order of the g-modes. Moreover, each minimum in ΔP corresponds to a minimum of inertia, linked to a pressure-dominated mode, and the higher the frequency the lower the ΔP. We also calculate the coupling parameter (q) at ν_max using the prescriptions of Takata (2016), yielding q ≈ 0.24, which aligns with the observational findings (e.g. Vrard et al. 2016; Mosser et al. 2017).

In the hydrogen-burning shell, we observe that the Brunt-Väisälä frequency is non-differentiable at the point where θ = θ_P, a condition resulting from the specific choice of γ(ρ) employed in equation (16). This non-differentiability has the potential to introduce an artificial glitch in the period spacing of the mixed modes, as noted in previous studies (e.g. Cunha et al. 2015). Therefore, it warrants careful examination. In Appendix D, we argue that the small deviations in ΔP_a noted in our model are likely attributable to numerical inaccuracies in the computation of eigenfrequencies rather than the non-differentiable point itself. Notably, the maximum observed peak-to-peak relative difference in period spacing of the γ-modes is at most 0.2 %, which is significantly lower than the typical observational uncertainties (e.g. ≈3 s in the asymptotic period spacing of CHeB stars in the Kepler database, Vrard et al. 2016, 2022). This indicates that the glitches identified in our model do not compromise the interpretation of observed power spectral densities. Additionally, this suggests that the interpretation of other glitches, which is discussed in subsequent sections, are unlikely to be affected by this specific instance of non-differentiability or by numerical errors. Furthermore, it is important to point out that the hydrogen-burning shell, which is characterised by a half width at half maximum (HWHM) of approximately 0.3H_P in our model, is incapable of producing a glitch signature in the eigenfrequencies. Any potential evidence of such a signature may only be detected in the very early stages of the CHeB phase, particularly when the hydrogen-burning shell is relativelythin.

5.2. Different boundary of the convective core

In this section, we explore the impact on frequencies of a model similar to the fiducial barotropic model, but with a smooth transition in the density profile between the convective and radiative core (Section 4). In Figure 5 we show in the top panel its N profile (red line) as a function of the internal mass in comparison with the fiducial model (blue line). This new profile is very similar to the fiducial model, but it contains a bell-shaped structure centred at the boundary instead of a δ-distribution. We expect a value of the asymptotic period spacing of the new model ΔP_a, Smooth lower than the one of the fiducial model ΔP_{a, Fiducial}, because the new model has a slighlty greater g-cavity than the fiducial model and it has higher values of N near the boundary. Moreover, we expect that such a bell-shaped structure in the N profile becomes a glitch, which creates a sinusoidal behaviour of ΔP around ΔP_a, Smooth with a decreasing amplitude at increasing period, and dips in the period spacing that are evenly spaced in period, with a periodicity compatible with the buoyancy radius of the glitch (Cunha et al. 2019, 2024). To verify that the model behaves as predicted by theory, we calculate eigenmodes also outside the observable region of thespectrum.

Fig. 5. Comparison of N profiles as functions of internal mass (top panel) and period spacings of the γ-modes as functions of the eigenperiods (bottom panel). The blue line represents the fiducial model, while the red model features a bell-shaped structure in the N profile at the boundary between the convective and radiative core (see Section 4). The periodicity of the glitch corresponds to the location of the bell-shaped structure. The half width at half maximum of the bell-shaped structure presented here is 0.0012HP.

In the bottom panel we show the period spacings of the γ-modes as functions of the eigenperiods for this new model (in red) and the fiducial model (in blue). As expected, there is a very small difference between the two asymptotic period spacings, that is, ΔP_a, Smooth ≈ ΔP_{a, Fiducial} − 0.65 s. Moreover, we observe in the new model a periodicity of approximately 23 hours with a decreasing amplitude at increasing eigenperiods. When expressed in terms of radial orders, this same periodicity corresponds to Δn ≈ 334, which aligns with the expected value derived from applying equation (3) to the new model. According to the same equation, we also infer that this glitch signature correlates to a structural variation situated at Φ(r_glitch)≈0.0030. This finding is consistent with the normalised buoyancy radius of the bell-shaped structure in N, which peaks at Φ ≈ 0.0015 and possesses a FWHM of 0.0023. To further support the relation between the bell-shaped structure in N and the observed glitch signature, we compare the Brunt-Väisälä length scale (H_N) with the local wavelength of g-modes (λ_g, as discussed in Section 1) near ν_max.

We find that H_N of the bell-shaped structure is significantly smaller than λ_g at Φ(r_glitch)≈0.0030 (see Figure 6). It is important to highlight that within the observable region of the oscillation spectrum, this periodic component would appear too smooth to be identified as a glitch signature. Nonetheless, it serves as a useful case study for testing, and we cannot rule out the possibility that signatures of convective boundary mixing may be observed in more evolved CHeB stars, where glitches at the boundary between the convective and radiative core could potentially be more pronounced (see Section 5.4). Finally, the small glitch explained in Section 5.1 adds to the above glitch, generating relative fluctuations much lower than 1 % around the main behaviour of the period spacing.

Fig. 6. Ratio of the local wavelength of g-modes (λg) near νmax to the Brunt-Väisälä length scale (HN), shown in black, as a function of the normalised buoyancy radius (Φ) for a model characterised by a bell-shaped structure in the N profile at the boundary between the convective and radiative core (refer to Section 4). The dashed black line indicates the threshold where λg/HN = 1. Additionally, we present the N profile (red line) as a function of Φ.

5.3. Glitches in the radiative core

In this section, we explore the impact on frequencies of two different types of glitches in the radiative core. In Section 5.3.1 we study the impact of a jump discontinuity in density, whereas in Section 5.3.2 we smooth the same discontinuity using the method of Section 4.

5.3.1. δ-distribution structural glitch

In Figure 7, the top panel displays the N profile as a function of internal mass for a model similar to that presented in Appendix E. This model features a continuous, albeit non-differentiable, density function at the boundary between the convective and radiative core (see Appendix E for more details), and it introduces a new jump discontinuity in density that is compatible with a subflash event occurring within the radiative core or with the chemical gradient found at the end of an overshoot region that possesses radiative thermal stratification. This conclusion is drawn from the observation that the jump present in the reference CLES model at the outer edge of the helium core closely resembles that at the boundary between the convective and radiative core. We have selected a position of r = 0.03 R_⊙ for the jump, as it correlates well with a main flash event identified in other evolutionary codes (Kippenhahn et al. 2012). The total stellar radius of our barotropic models is 9.9457 R_⊙, which provides a reference scale for interpreting the position discussed above. The resulting N profile (red line) closely resembles the fiducial model (blue line); however, the red model is characterised by a jump discontinuity at the boundary between the convective and the radiative core rather than a δ-distribution. Moreover, a δ-distribution is observed at the location of the jump discontinuity in density within the radiative core, with the bell-shaped peak associated with the hydrogen-burning shell occurring at a lower radius compared to the fiducial model. Overall, the integral ∫_r1^xN/r dr at a fixed position x ≤ r₂ in this new model resembles that of the Appendix E mentioned earlier. However, the jump discontinuity in density within the radiative core leads to a reduction in ∫_r1^xN/r dr for positions above the discontinuity when compared to the other case. This cumulative effect drives the value of the integral ∫_r1^r₂N/r dr below that of the fiducial model. Consequently, we expect a value of the asymptotic period spacing ΔP_a, Discont higher than ΔP_{a, Fiducial}. Furthermore, we expect that the δ-distribution in the N profile will manifest as a glitch, resulting in periodic deviations of ΔP relative to ΔP_a, Discont (e.g. Cunha et al. 2015). These deviations will exhibit a decreasing amplitude with increasing frequency³, creating evenly spaced dips in period that align with the periodicity associated with the normalised buoyancy radius of the glitch. Additionally, some modes are expected to become trapped in regions close to the density discontinuity as a result of the abrupt change in N. To verify that the model behaves as predicted by theory, we calculate eigenmodes also outside the observable region of the spectrum.

Fig. 7. Comparison of the N profiles as functions of internal mass (top panel) and the period spacings of the γ-modes as functions of the eigenperiods (bottom panel). The blue line represents the fiducial model, while the red model incorporates a δ-distribution glitch in the N profile located within the radiative core. The fiducial model is not displayed in the bottom panel, as it would be indistinguishable from the black line, which represents the ΔPa of the red model. Notably, the periodicity in the bottom panel corresponds to the location of the δ-distribution, and the maximum period spacings measured closely match ΔPa, until Discont. This value reflects the asymptotic ΔP we would obtain if the g-cavity extends from the position of the glitch to the outer boundary.

In the bottom panel, the period spacing of the γ-modes as a function of the eigenperiods for this new model. The fiducial model is not displayed in the bottom panel, as it would be indistinguishable from ΔP_a, Discont, which represents the ΔP_a of the red model. As expected, there is a small difference between the two asymptotic period spacings, specifically ΔP_a, Discont ≈ ΔP_{a, Fiducial} + 0.46 s. Moreover, the period spacing exhibits a periodicity of about 51 minutes; when expressed in terms of radial orders, this periodicity translates to Δn ≈ 12.4, in accordance with the expected value derived from the application of equation (3) to the model. This analysis also suggests that the glitch signature corresponds to a structural variation situated at Φ(r_glitch)≈0.080, which is compatible with the normalised buoyancy radius of the δ-distribution. Notably, the maximum period spacings measured closely match ΔP_{a, until Discont}. This value reflects the asymptotic ΔP we would obtain if the g-cavity extends from the position of the glitch to the outer boundary. It is evident that the inferred ΔP_a from observations would resemble ΔP_{a, until Discont} rather than the true ΔP_a, as expected (e.g. Cunha et al. 2015). Finally, the small glitch explained in Section 5.1 adds to the above glitch, generating relative fluctuations much lower than 1 % around the main behaviour of the period spacing.

After testing our model, in the top panel of Figure 8 we compare the period spacing of the γ-modes (in blue) with the period spacing of the mixed dipole modes (in red) in the observable region of the spectrum. Notably, some minima in the period spacing result from mode trapping rather than the coupling between p-modes and g-modes. This suggests that certain trapped modes may be detectable in actual data. However, not all of these minima show observable amplitudes, as modes with higher inertia tend to have lower amplitudes (e.g. Dupret et al. 2009; Grosjean et al. 2014), leading to decreased detectability even when the period spacing is low. Consequently, the measured ΔP is likely to differ from the true value if some modes are absent. This situation becomes even clearer in the bottom panel of Figure 8, where we present the PSD (along with the normalised inertiae) for the model featuring a density discontinuity within the radiative core.

Fig. 8. Period spacing (top panel) as a function of eigenfrequencies for the model discussed in Section 5.3.1. We compare the γ-modes (shown in blue) with mixed dipole modes (in red) in proximity to the observable region of the oscillation spectrum, whose νmax is the cyan line. The black and lime lines are the same as Figure 7. Some minima in the period spacing result from mode trapping rather than the coupling between p-modes and g-modes, and not all of these minima exhibit observable amplitudes. This is further illustrated in the bottom panel, which presents the corresponding simulated PSD along with the normalised inertia of radial and non-radial modes (coloured dots connected by coloured lines).

5.3.2. Bell-shaped structural glitch

In the top panel of Figure 9, we show the N profile as a function of the internal mass for a model similar to the model of Section 5.3.1, but with a smooth transition in density instead of a jump discontinuity (i.e. with α = 10², see Section 4). This new model has a N profile (red line) very similar to the fiducial model (blue line), but it contains a bell-shaped structure where there is the smooth change in density. This structure increases the integral ∫_r1^r₂N/r dr compared to the fiducial model. Thus, we expect a value of the asymptotic period spacing ΔP_a, Discont lower than ΔP_{a, Fiducial}. Moreover, we expect that such bell-shaped structure in the N profile becomes a glitch, because its width is much lower than the local wavelength of the waves (Cunha et al. 2019, 2024). Therefore, we expect a sinusoidal behaviour of ΔP around ΔP_a, Smooth with a decreasing amplitude at increasing period, and dips in the period spacing that are evenly spaced in period, with a periodicity compatible with the normalised buoyancy radius of the glitch. This structure is aligned identically with the glitch discussed in Section 5.3.1, as we have placed the peak of the bell function at a location analogous to the δ-distribution. To verify that the model behaves as predicted by theory, we calculate eigenmodes also outside the observable region of the spectrum.

Fig. 9. Comparison of N profiles as functions of internal mass (top panel) and period spacings of the γ-modes as functions of the eigenperiods (bottom panel). The blue line represents the fiducial model, while the red model features a bell-shaped structure in the N profile within the radiative core. The half width at half maximum of the bell-shaped structure in N presented here is 0.02HP.

In the bottom panel, the period spacing of the γ-modes as a function of the eigenperiods for this new model. The fiducial model is not displayed in the bottom panel, as it would be indistinguishable from ΔP_a, Smooth, which represents the asymptotic period spacing of the model with a smooth transition. As expected, there is a small difference between the two asymptotic period spacings, that is ΔP_a, Smooth ≈ ΔP_{a, Fiducial} − 0.86 s. The period spacing shows a periodicity of approximately 51 minutes with a decreasing amplitude at increasing eigenperiods. This periodicity is compatible with the observed Δn ≈ 12.4 and corresponds to a normalised buoyancy radius of 0.081, as described in equation (3). This indicates consistency with the normalised buoyancy radius of the bell-shaped structure in N, which peaks at Φ ≈ 0.082 and possesses a FWHM ≈ 0.011. We further reinforce the alignment of these findings with the characteristics of the bell-shaped structure by checking the relation between the bell-shaped structure and the observed glitch signature similarly to Section 5.2.

We find the presence of two structure variations where H_N is significantly smaller than λ_g. These two glitches are located at Φ ≈ 0.078 and Φ ≈ 0.085 (see Figure 10). Therefore, contrary to the case with the δ-distribution, we expect a glitch signature in the period spacing with visible beats attributable to these two adjacent sharp variations, which is clearly visible in the bottom panel of the figure. Furthermore, applying a Fouriertransform to the period spacing data demonstrates that the two primary frequencies responsible for the beats correspond with the aforementioned peaks, collectively accounting for the observed periodicity. Finally, the small glitch explained in Section 5.1 adds to the above glitch, generating relative fluctuations much lower than 1 % around the main behaviour of the period spacing.

Fig. 10. Ratio of the local wavelength of g-modes (λg) near νmax to the Brunt-Väisälä length scale (HN), shown in black, as a function of the normalised buoyancy radius (Φ) for a model characterised by a bell-shaped structure in the N profile within the radiative core (refer to Section 4). The dashed black line indicates the threshold where λg/HN = 1. Additionally, we present the N profile (red line) as a function of Φ.

A more detailed characterisation of the detectability of the modes and the ability to differentiate between varying degrees of smoothing of the bell-shape structure is presented in Figure 11. This figure compares two period spacings of the mixed dipole modes within the observable region of the spectrum derived from different levels of smoothing (i.e. α = 10² and α = 10⁴). Notably, we observe clear differences in their period spacings, even in areas where dipole modes are detectable. This finding highlights promising directions for interpreting glitch signatures in high-quality asteroseismic data, including thepossibility to discern the sharpness of the density jump. However, one has to remember that not all of these minima show observable amplitudes and the measured ΔP is likely to differ from the true value if some modes are absent (as explained inSection 5.3.1).

Fig. 11. Comparison of two period spacings of mixed dipole modes within the observable region of the spectrum derived from different levels of smoothing applied in the radiative core. The cyan dashed line indicates νmax, while the black dashed line denotes the asymptotic ΔP. Additionally, the grey lines correspond to the radial modes. The shown ΔPa and radial modes are representative of both models, as they are indistinguishable in this scale. Notably, we observe clear differences in their period spacings, even in areas where dipole modes are detectable.

5.4. Evolved CHeB stars

All the barotropic models examined are derived from evolutionary models at the beginning of the CHeB phase (i.e. with Y_c ≈ 0.9). However, as this evolutionary phase progresses, a growing density gradient is expected at the boundary between the convective and radiative core (see e.g. Kippenhahn et al. 2012). Moreover, if the overshoot region exhibits radiative stratification, this density discontinuity is expected to occur within the radiative core, potentially giving rise to glitches within the observable region of the oscillation spectrum. In this section, we investigate the impact on the eigenfrequencies of convective overshooting beyond the boundary between the convective and radiative core in the context of more evolved stellar models. The overshoot region we model possesses radiative thermal stratification. To mitigate complications arising from the appearance of semi-convective layers, we evolve stellar models using CLES and MESA until reaching Y_c ≈ 0.6. We then establish this new density gradient at the boundary as the new Λ, while maintaining the same total mass and radius⁴ as the fiducial model of Section 3.

The approach for modelling the additional structural glitches follows the methodologies outlined in Sections 5.2 and 5.3.

The main results concerning the observability of structural glitches across the seven models are summarised in Table 2. This table details the periodicity of the detected glitches, if present, along with their inferred locations. In all these cases the observed glitches align with the predictions made by equation (3), but a strict interpretation of this equation is not always useful in determining the observability of structural variations. For instance, a structural glitch situated precisely at the boundary between the convective and radiative core would still yield a glitch signature in the eigenfrequencies, albeit with an infinite periodicity. Therefore, even if a structural variation exists, it would remain undetectable. A similar condition arises for the case where α = 10; however, in this case, no glitch signature is observable because the structural variation lacks sufficient sharpness. This scenario also pertains to the hydrogen-burning shell discussed in Section 5.1. Finally, when the parameter α assumes a finite value, beats emerge in the period spacing similarly to Section 5.3.2.

The amplitude of the glitch signatures presented in this section is higher than that observed in Sections 5.2 and 5.3 due to the higher Λ density jump. Such behaviour is consistent with the theoretical expectations (e.g. Miglio et al. 2008; Cunha et al. 2015). Consequently, it may be possible to observe these glitch signatures in the eigenfrequencies depending on the helium content in the centre of the CHeB star, as this content affects the density jump at the core boundary. Additionally, the observability of such glitches is influenced by the total radius of the CHeB star and/or its effective temperature, as both parameters alter ν_max for a given total stellar mass. Ultimately, one has to consider that the observable region of the oscillation spectrum is limited, and not all modes within this range are detectable. Indeed, there is an inverse relationship between the inertia of a mode and its amplitude (e.g. Dupret et al. 2009). Consequently, the inability to detect certain modes could hinder accurate identification of the location of the structural glitches and of the period spacing. Therefore, even when distinct differences among models are evident in the period spacing of the mixed dipole modes within the observable region of the spectrum (see Figures 12 and 13), interference from glitch signatures with p-like dipole modes may further reduce the observable amplitude of the dipole modes, potentially resulting in their absence in the observed frequencies. These collective insights underline the complexities encountered in detecting glitch signatures in actual data (e.g. Vrard et al. 2022).

Fig. 12. Comparison of two period spacings of mixed dipole modes within the observable region of the spectrum derived from different conditions at the core boundary. In particular, we show a model with a density jump (blue solid line), and a model with a smooth transition (orange solid line). The dashed lines denote the asymptotic period spacings for the model with a density jump (in black), and for the other model (in red). Furthermore, the grey lines correspond to the radial modes associated with both models, as they are indistinguishable in this scale. Notably, we observe clear differences in their period spacings, even in areas where dipole modes are detectable.

Fig. 13. Comparison of two period spacings of mixed dipole modes within the observable region of the spectrum. In particular, we show a model with a density jump located at the boundary between the convective and radiative core (blue solid line), and a model with a density jump located within the radiative core (orange solid line). The dashed lines denote the asymptotic period spacings for the model with a density jump (in black), and for the other model (in red). Furthermore, the grey lines correspond to the radial modes associated with the model with a density jump in the radiative core. Notably, we observe clear differences in their period spacings, even in areas where dipole modes are detectable.

6. Summary and conclusions

In this paper, we have conducted a theoretical analysis of how structural variations adjacent to the convective core and chemical composition gradients within the radiative core influence the period spacing of mixed dipole modes and γ-modes in low-mass CHeB stars. These variations are also expected to occur within semi-convective layers, within the hydrogen-burning shell, and at the base of the convective envelope. Additionally, in low-mass stars with a degenerate helium core, the transition from RGB to CHeB is characterised by a succession of off-centre helium flashes, which induce chemical composition gradients in the radiative core. To investigate the impact of density discontinuities and associated structural glitches on the period spacings of these oscillation modes, we developed semi-analytical models of low-mass CHeB stars. These models were calibrated using the advanced evolutionary codes BaSTI-IAC, CLES, and MESA. The distinct physical prescriptions of these codes allowed us to identify common relevant features for calibration while maintaining control over the type of structural glitch introduced and keeping a realistic representation of stellar interiors. We first established a fiducial model based on a solar mass CHeB star at solar metallicity with Y_c ≈ 0.9 featuring a realistic g-cavity. Subsequently, we explored the effects of various discontinuities in density, non-differentiable points within the density profile, and bell-shaped glitches in the Brunt–Väisälä frequency on the adiabatic eigenfrequencies. The key findings from our analysis are summarised as follows:

• Jump discontinuities in the density profile result in distinct glitch signatures in the period spacings of mixed modes. These glitches exhibit periodic behaviours that are closely tied to the normalised buoyancy radius associated with the discontinuities. While certain trapped modes may be observable, their detectability is influenced by their inertia, with higher inertia modes exhibiting lower amplitudes. This may lead to discrepancies between measured and true values of ΔP due to possible missing modes in observations. Notably, our analysis indicates that the inferred period spacing from observational data is more likely to represent the maximum measured value rather than the actual asymptotic one. However, they are typically close to each other.

• The comparison between models featuring smooth transitions and those with discontinuities highlighted differences in the periodic behaviours. Specifically, this finding suggests that smooth transitions can have an impact comparable to that of jump discontinuities (i.e. be characterised as structural glitches) if the scale height of the transition is much shorter than the local wavelength of the waves under investigation, thus presenting opportunities to detect not only the position and amplitude of the glitches but also theirsharpness.

• Our simulations of four-year-long Kepler observations reinforced the necessity of incorporating realistic stellar interior models to predict oscillation frequencies accurately. The resulting PSDs and period spacing patterns closely resemble observed data, providing a promising avenue for verifying our theoretical models against real-world observations.

• Additionally, we investigated more evolved CHeB stars with Y_c ≈ 0.6. Our results indicate that the observability of glitch signatures is influenced by the helium content at the centre of the star, its total radius, and/or its effective temperature. We also find that not all modes within the observable region of the spectrum are detectable, and the presence of glitch signatures near p-like mixed modes can decrease their observable amplitudes, complicating the correct identification of glitch signatures in actual data.

This work establishes a solid foundation for future asteroseismic studies aimed at probing the internal structures of stars. Our models enable realistic predictions of how each sharp structural variation impacts the observed power spectral density. This alignment not only validates our theoretical approach but also suggests promising directions for interpreting glitch signatures in high-quality asteroseismic data, such as that obtained from the Kepler mission, as well as from future space missions, for example, the upcoming PLATO mission (Rauer et al. 2025) and the HAYDN project (Miglio et al. 2021).

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author.

Acknowledgments

We are grateful to Margarida Cunha, Arlette Noels and Richard Scuflaire for useful discussions. MM, AM, and WEvR acknowledge support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, https://www.asterochronometry.eu, G.A. n. 772293). GB acknowledges fundings from the Fonds National de la Recherche Scientifique (FNRS) as a postdoctoral researcher. LP acknowledges fundings from the FNRS as a PhD student. We thank the anonymous referee for the helpful comments that improved the quality of the manuscript.

References

Appendix A: Taylor series solutions near the centre of barotropic stars

In Section 2.1.1 we discuss the set of differential equations we use to solve the internal structure of a barotropic star knowing the initial conditions. However, we cannot solve the differential equations numerically from the exact centre, because there the term ψ/ξ² would lead the Runge-Kutta method to errors. To avoid this problem, we can begin the numerical evaluation of the equations near the centre using Taylor expansions. For simplicity, let us assume for now that γ(θ) = γ_c(θ/θ_c)^B near the centre of the star. We then obtained

$$\begin{array}{c}\hfill \{\begin{array}{c}\theta (\xi )={\theta }_c-\frac{{\theta }_c^3{\xi }^2}{6{\beta }_c{\gamma }_c}+\frac{{\theta }_c^5{\xi }^4}{360{\beta }_c^2{\gamma }_c^2}[13-5(B+{\gamma }_c)]+o({\xi }^6)\hfill \\ \frac{d\theta (\xi )}{d\xi }=-\frac{{\theta }_c^3\xi }{3{\beta }_c{\gamma }_c}+\frac{{\theta }_c^5{\xi }^3}{90{\beta }_c^2{\gamma }_c^2}[13-5(B+{\gamma }_c)]+o({\xi }^5)\hfill \\ \psi (\xi )=\frac{{\theta }_c{\xi }^3}{3}-\frac{{\theta }_c^3{\xi }^5}{30{\beta }_c{\gamma }_c}+\frac{{\theta }_c^5{\xi }^7}{2520{\beta }_c^2{\gamma }_c^2}[13-5(B+{\gamma }_c)]+o({\xi }^9),\hfill \end{array}\end{array}$$(A.1)

with which we completely solve the numerical issue. Moreover, it can be inferred from equation A.1 that, in close proximity to the centre of a barotropic star, γ(r) can be regarded as approximately constant. Therefore, very near the centre, the star can be effectively modelled as a polytrope.

Appendix B: Differential equations near the surface of barotropic stars

The set of differential equations we use in Section 2.1.1 creates numerical issues at both the centre of the star (as discussed in Appendix A) and at its surface. The issue observed at the surface arises from the ratio θ²/β, where both θ and β approach zero. Therefore, it is essential to formulate an alternative set of differential equations that avoids the issue associated with θ²/β in this region. To achieve this, we introduced a Lane-Emden-like variable (ω) defined as

$$\begin{array}{c}\hfill \frac{\partial \omega (\theta )}{\partial \theta }:=\frac{\beta (\theta )}{{\theta }^2}\end{array}$$(B.1)

or, equivalently,

$$\begin{array}{c}\hfill \omega (\theta ):={\omega }_c+{\displaystyle {\int }_{{\theta }_c}^{\theta }}\frac{\beta (t)}{t^2}dt.\end{array}$$(B.2)

By substituting equation B.1 (or equation B.2) along with equation 10 into the expression for dω/dξ, we obtained

$$\begin{array}{c}\hfill \frac{d\omega }{d\xi }=-\frac{\psi (\xi )}{\gamma (\omega ){\xi }^2}.\end{array}$$(B.3)

The resulting general differential equation, which combines equation 10 with equation B.3, takes the form

$$\begin{array}{c}\hfill \frac{d^2\omega }{d{\xi }^2}+{\left(\frac{d\omega }{d\xi }\right)}^2\frac{1}{\gamma (\omega )}\frac{\partial \gamma (\omega )}{\partial \omega }+\frac{2}{\xi }\frac{d\omega }{d\xi }+\frac{\theta (\omega )}{\gamma (\omega )}=0,\end{array}$$(B.4)

where the new initial conditions are ω(ξ = 0) = ω_c and ${\frac{d\omega }{d\xi }|}_{\xi =0}=0$. This new equation overcomes the mentioned numerical issue, provided that γ ≠ 0, because as θ approaches zero we do not have to deal with an indeterminate form.

Equation B.4 simplifies when γ(ω) is a constant equal to γ_c because we obtained the Lane-Emden-like equation

$$\begin{array}{c}\hfill \frac{d^2\omega }{d{\xi }^2}+\frac{2}{\xi }\frac{d\omega }{d\xi }+\frac{\theta (\omega )}{{\gamma }_c}=0,\end{array}$$(B.5)

with

$$\begin{array}{c}\hfill \beta (\theta )={\beta }_c{\left(\frac{\theta }{{\theta }_c}\right)}^{{\gamma }_c}\end{array}$$(B.6)

and

$$\begin{array}{c}\hfill \theta (\omega )=\{\begin{array}{c}{\left[\frac{({\gamma }_c-1){\theta }_c^{{\gamma }_c}}{{\beta }_c}\omega \right]}^{\frac{1}{{\gamma }_c-1}}:={\varphi }^{\frac{1}{{\gamma }_c-1}},\mathrm{for}{\gamma }_c\ne 1\hfill \\ {\theta }_{c}exp\left(\frac{{\theta }_c}{{\beta }_c}\omega \right),\mathrm{for}{\gamma }_c=1.\hfill \end{array}\end{array}$$(B.7)

From equation B.7 one obtains that ϕ_c = θ_c^{γ_c − 1} for γ_c ≠ 1, and that ω_c = 0 for γ_c = 1. Therefore, equation B.5 can be divided into the two differential equations

$$\begin{array}{c}\hfill \{\begin{array}{c}\frac{d^2\varphi }{d{\xi }^2}+\frac{2}{\xi }\frac{d\varphi }{d\xi }+\frac{({\gamma }_c-1){\theta }_c^{{\gamma }_c}}{{\beta }_c{\gamma }_c}{\varphi }^{\frac{1}{{\gamma }_c-1}}=0,\varphi \in [{\theta }_c^{{\gamma }_c-1},0],{\gamma }_c\ne 1\hfill \\ \frac{d^2\omega }{d{\xi }^2}+\frac{2}{\xi }\frac{d\omega }{d\xi }+{\theta }_{c}exp\left(\frac{{\theta }_c}{{\beta }_c}\omega \right)=0,\omega \in [0,-\infty ),{\gamma }_c=1,\hfill \end{array}\end{array}$$(B.8)

and from $(\varphi ,\frac{d\varphi }{d\xi })$ or $(\omega ,\frac{d\omega }{d\xi })$ we could a posteriori evaluate

$$\begin{array}{c}\hfill \psi (\xi )=\{\begin{array}{c}\frac{{\beta }_c{\gamma }_c}{({\gamma }_c-1){\theta }_c^{{\gamma }_c}}(-{\xi }^2\frac{d\varphi }{d\xi }),{\gamma }_c\ne 1\hfill \\ -{\xi }^2\frac{d\omega }{d\xi },{\gamma }_c=1,\hfill \end{array}\end{array}$$(B.9)

and

$$\begin{array}{c}\hfill N^2(\xi )=\{\begin{array}{c}4\pi G{\rho }_1\frac{{\beta }_c}{{\theta }_c^{{\gamma }_c}\varphi (\xi )}{\left(\frac{{\gamma }_c}{{\gamma }_c-1}\frac{d\varphi }{d\xi }\right)}^2(\frac{1}{{\gamma }_c}-\frac{1}{{\mathrm{\Gamma }}_1}),{\gamma }_c\ne 1\hfill \\ 4\pi G{\rho }_1\frac{{\theta }_c}{{\beta }_c}{\left(\frac{d\omega }{d\xi }\right)}^2\left(\frac{{\mathrm{\Gamma }}_1-1}{{\mathrm{\Gamma }}_1}\right),{\gamma }_c=1.\hfill \end{array}\end{array}$$(B.10)

Appendix C: Numerical solver verification

We validate the numerical solver thanks to known analytical solutions to equation 10 when γ(ρ) = γ_c is a constant. One solution is the Plummer sphere (i.e. γ_c = 6/5), which is

$$\begin{array}{c}\hfill \{\begin{array}{c}\theta (\xi )={\theta }_c{(1+\frac{{\theta }_c^2{\xi }^2}{18{\beta }_c})}^{-\frac{5}{2}}\hfill \\ \beta (\theta )={\beta }_c{\left(\frac{\theta }{{\theta }_c}\right)}^{\frac{6}{5}}={\beta }_c{(1+\frac{{\theta }_c^2{\xi }^2}{18{\beta }_c})}^{-3}\hfill \\ \psi (\xi )=\frac{{\theta }_c{\xi }^3}{3}{(1+\frac{{\theta }_c^2{\xi }^2}{18{\beta }_c})}^{-\frac{3}{2}}.\hfill \end{array}\end{array}$$(C.1)

Another solution is for γ_c = 2, which is

$$\begin{array}{c}\hfill \{\begin{array}{c}\theta (\xi )=\frac{\sqrt{2{\beta }_c}}{\xi }sin\left(\frac{{\theta }_c\xi }{\sqrt{2{\beta }_c}}\right)\hfill \\ \beta (\theta )={\beta }_c{\left(\frac{\theta }{{\theta }_c}\right)}^2=\frac{2{\beta }_c^2}{{\theta }_c^2{\xi }^2}{sin}^2\left(\frac{{\theta }_c\xi }{\sqrt{2{\beta }_c}}\right)\hfill \\ \psi (\xi )=\frac{{(2{\beta }_c)}^{\frac{3}{2}}}{{\theta }_c^2}[sin\left(\frac{{\theta }_c\xi }{\sqrt{2{\beta }_c}}\right)-\left(\frac{{\theta }_c\xi }{\sqrt{2{\beta }_c}}\right)cos\left(\frac{{\theta }_c\xi }{\sqrt{2{\beta }_c}}\right)].\hfill \end{array}\end{array}$$(C.2)

Appendix D: Tests on artificial glitches in the fiducial model

Upon analysing ΔP at frequencies lower than ν_max, we identify a behaviour resembling a small amplitude oscillation centred in ΔP_a. The periodicity of this oscillation is approximately 11.6 minutes, while the corresponding periodicity in radial order is Δn ≈ 2.80. These periodicities are consistent with expectations based on equation 3. However, it is important to note that the Brunt-Väisälä frequency at the normalised buoyancy radius associated with Δn ≈ 2.80 is smooth, which implies that it should not affect the eigenfrequencies with a glitch signature (e.g. Cunha et al. 2015). Notably, this periodicity is accompanied by an alias observed at Δn ≈ 1.555. When this alias is taken into consideration, we find that the associated normalised buoyancy radius is in good agreement with the characteristics of the hydrogen-burning shell, which peaks at Φ ≈ 0.641 and has a full width at half maximum (FWHM) of 0.126. Nonetheless, it is crucial to remain aware that this periodicity may also arise from numerical inaccuracies in the computation of the adiabatic frequencies and/or from the presence of the aforementioned non-differentiable point.

A way to solve this problem is by a comparison between the local wavelength of the eigenfunctions near the shell and the width of the glitch (or the scale of variation of N, e.g. Cunha et al. 2019), because it shows whether there is a period trapping in the region and what the main source of the trapping is. A simpler way to see buoyancy glitches, and conversely to solve the issue, is to measure the ΔP of the γ-modes. As can be seen from the top panel of the Figure 4 and from Figure D.1, the small amplitude oscillations centred in the asymptotic period spacing ΔP_a are more visible when the modulation caused by the coupling with the p-modes is absent. While the aforementioned periodicity in radial order of these oscillations persists in the γ-modes, we must also account for a newly prominent glitch signature with periodicity Δn ≈ 2, as evidenced by a Fourier transform applied to the period spacing of the oscillation modes. According to equation 3, we would expect a glitch structure at Φ(r_glitch)≈0.5. However, the N function is smooth at this location, which means that it cannot generate a glitch signature in the eigenfrequencies. It is now apparent that the glitch at Δn ≈ 1.555 aligns with the non-differentiable point, as at this specific location the scale of variation of N is much lower than the local wavelength of the g-modes near ν_max. At this point, we can devise a simplified model that eliminates such non-differentiable point, albeit at the expense of a less realistic g-cavity. In particular, we employed the equations

Fig. D.1. Comparison between normalised period spacings for the γ-modes as functions of frequency. The grey area represents the observable region of the spectrum and the cyan line the νmax. It is clear that the fiducial model (in black) and the simplified model (in red) have similar properties of the g-cavity.

$$\begin{array}{c}\hfill \{\begin{array}{c}\gamma (\theta )=Aexp[-\frac{{(ln\theta -ln{\theta }_P)}^2}{2{\sigma }^2}]+c,\theta \le 1\hfill \\ \hfill \\ ln\beta (\theta )=A\sqrt{\frac{\pi }{2}}\sigma [\mathrm{erf}\left(\frac{ln\theta -ln{\theta }_P}{\sqrt{2}\sigma }\right)+\mathrm{erf}\left(\frac{ln{\theta }_P}{\sqrt{2}\sigma }\right)]+\hfill \\ +cln\theta ,\theta \le 1\hfill \end{array}\end{array}$$(D.1)

instead of the equations 16 and 17, with c > 0 a constant chosen in order to have the correct total mass of the star. Figure D.1 illustrates the normalised period spacings for the γ-modes from both the fiducial model (in black) and the simplified model. The results demonstrate that the observed glitch signature with periodicity Δn ≈ 2 remains intact, and that this oscillation appears to be at least compatible, if not predominant, with the glitch signature originating from the non-differentiable point. To further investigate this, we employed a finer grid and utilised different numerical solvers. Despite these modifications, the same glitches in the period spacing were consistently observed. In particular, increasing the number of grid points by an order of magnitude results in a relative change in period spacing of the order of 10⁻²%, while switching to a lower order numerical solver results in a relative change of the order of 10⁻⁵%. Once the model has converged, these deviations from the mean period spacing do not decrease as the number of grid points increases. Therefore, the main deviations from the asymptotic value could arise from numerical inaccuracies in the computation of eigenfrequencies, rather than from the presence of the non-differentiable point. Notably, the peak-to-peak relative difference in the period spacing of the modes is at most 0.2 % (in the high-frequency regime), which is significantly lower than typical observational uncertainties (e.g. ≈3 s in the asymptotic period spacing of CHeB stars in the Kepler database, Vrard et al. 2016, 2022). Consequently, this glitch does not compromise the interpretation of the observed power spectral distributions. Furthermore, it confirms that the other glitches (which are addressed in Sections 5.2, 5.3, and 5.4) can be examined without concern for interference from this particular glitch, as those will likely be much more significant.

Appendix E: Non-differentiable continuous density function at the boundary between convective and radiative core

In Figure E.1 we show in the top panel the N profile as a function of the internal mass for a model similar to the fiducial barotropic model of Section 3, but with a continuous density function that is non-differentiable at the boundary between the convective and the radiative core. This has been done by choosing Λ = 1 at the boundary. This new model has a N profile (red line) very similar to the fiducial model (blue line), but it contains a jump discontinuity at the boundary instead of a δ-distribution, and the bell-shaped peak related to the H-burning shell is located at a lower radius than in the fiducial model. Therefore, ∫_r1^r₂N/r dr is higher in the new model and we expect a value of the asymptotic period spacing ΔP_a, NO lower than ΔP_{a, Fiducial}.

Fig. E.1. Comparison of N profiles as functions of internal mass (top panel) and period spacings of the γ-modes as functions of the eigenfrequencies (bottom panel). The blue line represents the fiducial model of Section 3, while the red model features a step-like structure in the N profile at the boundary between the convective and radiative core instead of a δ-distribution.

In the bottom panel, the period spacings are presented as functions of the eigenfrequencies for both the fiducial model (in blue) and the new model (in red). Although the two asymptotic period spacings differ, the difference is minimal (≈0.76 s). Notably, the small glitch discussed in Section 5.1 now shows an increased peak-to-peak relative difference in the period spacing of the modes. This could be related to the change of the N profile at the boundary. However, it is important to note that the JWKB approximation used by Cunha et al. (2019, 2024) cannot be applied in this context to infer the properties of the glitches. Despite this increased difference, the relative difference remains at most 1 % in the observable region of the spectrum. As a result, it is very difficult to detect such a glitch even with the highest quality data available, especially when analysing mixed dipole modes.

Appendix F: Modelling of the convective boundary mixing

To address the convective boundary mixing problem, we employ phenomenological modelling with the same MESA tool applied in Section 3. For these simulations, we focus on a 1 M_⊙ CHeB star with solar metallicity, employing identical undershooting prescriptions as described in Section 3. Specifically, we implement an exponential overshooting prescription (e.g. Maeder 1975; Freytag et al. 1996; Herwig 2000), in which the overshoot region is defined by a temperature gradient ∇ = ∇_rad. The chosen prescription modifies the diffusive mixing coefficient D_OV near the boundary of the convective core such that

$$\begin{array}{c}\hfill D_{\mathrm{OV}}=D_{0}exp\left(\frac{-2(r-r_0)}{0.035H_{p,cc}}\right),\end{array}$$(F.1)

where H_p, cc represents the pressure scale height at the convective core boundary, r is the radial coordinate, D₀ is the mixing coefficient at r₀ = r_cc − 0.005H_p, cc, and r_cc denotes the radius of the convective core boundary. Additionally, we implement a minimum mixing coefficient D_min = 1 cm²/s in the internal structure of the star.

All Tables

All Figures

Fig. 1. Brunt-Väisälä frequencies as a function of internal mass for five different star models at the beginning of the CHeB stage (top panel). The orange, green, red, and violet lines represent models calculated using the CLES, MESA, and BaSTI-IAC evolutionary codes (see Section 3), while the blue line corresponds to the fiducial barotropic model discussed in Section 3. All models are for a 1 M⊙ CHeB star with solar metallicity. The red line represents a MESA model computed with the boundary mixing prescriptions outlined in Appendix F. These models reveal common features in their Brunt-Väisälä frequencies, including a convective core, radiative core, hydrogen-burning shell, radiative envelope, and convective envelope (the latter not shown for clarity). The bottom panel shows the helium mass fraction as a function of internal mass for four of the five models.

In the text

	Fig. 2. Comparison of γ(r) as a function of density for the five different models shown in Figure 1. The dashed brown line represents the Γ1 used by CLES v21.0 for reference. All five models share the fundamental characteristics of the CHeB phase.
In the text

Fig. 3. Comparison between three different density profiles (top panel) at the boundary between the convective and the radiative core and their corresponding Brunt-Väisälä frequencies (bottom panel) as functions of internal mass. In particular, we compare the fiducial barotropic model (blue) with two models that smoothly join the two zones (orange and green). The orange line incorporates a lower α value than the green line. The blue arrow corresponds to a δ-distribution. The FWHM of the peaks in N presented here are significantly lower than 0.1Hp.

In the text

Fig. 4. Brunt-Väisälä frequency as a function of internal mass (top left), period spacing as a function of eigenfrequencies in the observable region of the oscillation spectrum (top right), and the corresponding simulated PSD (bottom) for the fiducial model discussed in Section 5.1. The cyan dashed line indicates the νmax of the reference CLES model. Along with the simulated PSD, we also show the normalised inertia of radial and non-radial modes (coloured dots connected by coloured lines).

In the text

Fig. 5. Comparison of N profiles as functions of internal mass (top panel) and period spacings of the γ-modes as functions of the eigenperiods (bottom panel). The blue line represents the fiducial model, while the red model features a bell-shaped structure in the N profile at the boundary between the convective and radiative core (see Section 4). The periodicity of the glitch corresponds to the location of the bell-shaped structure. The half width at half maximum of the bell-shaped structure presented here is 0.0012HP.

In the text

Fig. 6. Ratio of the local wavelength of g-modes (λg) near νmax to the Brunt-Väisälä length scale (HN), shown in black, as a function of the normalised buoyancy radius (Φ) for a model characterised by a bell-shaped structure in the N profile at the boundary between the convective and radiative core (refer to Section 4). The dashed black line indicates the threshold where λg/HN = 1. Additionally, we present the N profile (red line) as a function of Φ.

In the text

Fig. 7. Comparison of the N profiles as functions of internal mass (top panel) and the period spacings of the γ-modes as functions of the eigenperiods (bottom panel). The blue line represents the fiducial model, while the red model incorporates a δ-distribution glitch in the N profile located within the radiative core. The fiducial model is not displayed in the bottom panel, as it would be indistinguishable from the black line, which represents the ΔPa of the red model. Notably, the periodicity in the bottom panel corresponds to the location of the δ-distribution, and the maximum period spacings measured closely match ΔPa, until Discont. This value reflects the asymptotic ΔP we would obtain if the g-cavity extends from the position of the glitch to the outer boundary.

In the text

Fig. 8. Period spacing (top panel) as a function of eigenfrequencies for the model discussed in Section 5.3.1. We compare the γ-modes (shown in blue) with mixed dipole modes (in red) in proximity to the observable region of the oscillation spectrum, whose νmax is the cyan line. The black and lime lines are the same as Figure 7. Some minima in the period spacing result from mode trapping rather than the coupling between p-modes and g-modes, and not all of these minima exhibit observable amplitudes. This is further illustrated in the bottom panel, which presents the corresponding simulated PSD along with the normalised inertia of radial and non-radial modes (coloured dots connected by coloured lines).

In the text

	Fig. 9. Comparison of N profiles as functions of internal mass (top panel) and period spacings of the γ-modes as functions of the eigenperiods (bottom panel). The blue line represents the fiducial model, while the red model features a bell-shaped structure in the N profile within the radiative core. The half width at half maximum of the bell-shaped structure in N presented here is 0.02HP.
In the text

Fig. 10. Ratio of the local wavelength of g-modes (λg) near νmax to the Brunt-Väisälä length scale (HN), shown in black, as a function of the normalised buoyancy radius (Φ) for a model characterised by a bell-shaped structure in the N profile within the radiative core (refer to Section 4). The dashed black line indicates the threshold where λg/HN = 1. Additionally, we present the N profile (red line) as a function of Φ.

In the text

Fig. 11. Comparison of two period spacings of mixed dipole modes within the observable region of the spectrum derived from different levels of smoothing applied in the radiative core. The cyan dashed line indicates νmax, while the black dashed line denotes the asymptotic ΔP. Additionally, the grey lines correspond to the radial modes. The shown ΔPa and radial modes are representative of both models, as they are indistinguishable in this scale. Notably, we observe clear differences in their period spacings, even in areas where dipole modes are detectable.

In the text

Fig. 12. Comparison of two period spacings of mixed dipole modes within the observable region of the spectrum derived from different conditions at the core boundary. In particular, we show a model with a density jump (blue solid line), and a model with a smooth transition (orange solid line). The dashed lines denote the asymptotic period spacings for the model with a density jump (in black), and for the other model (in red). Furthermore, the grey lines correspond to the radial modes associated with both models, as they are indistinguishable in this scale. Notably, we observe clear differences in their period spacings, even in areas where dipole modes are detectable.

In the text

Fig. 13. Comparison of two period spacings of mixed dipole modes within the observable region of the spectrum. In particular, we show a model with a density jump located at the boundary between the convective and radiative core (blue solid line), and a model with a density jump located within the radiative core (orange solid line). The dashed lines denote the asymptotic period spacings for the model with a density jump (in black), and for the other model (in red). Furthermore, the grey lines correspond to the radial modes associated with the model with a density jump in the radiative core. Notably, we observe clear differences in their period spacings, even in areas where dipole modes are detectable.

In the text

	Fig. D.1. Comparison between normalised period spacings for the γ-modes as functions of frequency. The grey area represents the observable region of the spectrum and the cyan line the νmax. It is clear that the fiducial model (in black) and the simplified model (in red) have similar properties of the g-cavity.
In the text

	Fig. E.1. Comparison of N profiles as functions of internal mass (top panel) and period spacings of the γ-modes as functions of the eigenfrequencies (bottom panel). The blue line represents the fiducial model of Section 3, while the red model features a step-like structure in the N profile at the boundary between the convective and radiative core instead of a δ-distribution.
In the text