Low-to intermediate-mass stars (up to ∼ 8 M ⊙ on the main sequence), which account for the vast majority of stellar objects found in the Universe ( ∼ 98 % of them), ultimately end their lives as fading white dwarfs (see, e.g., the review of Fontaine and Wesemael 2000). White dwarfs are the remnants of former stellar cores, keeping engraved in their interior the imprint of past events and processes that occurred during stellar evolution (Salaris et al., 2010; De Gerónimo et al., 2017). Such “stellar fossils” therefore hold many keys to reach a deeper understanding of stars and their constitutive physics. In addition, they are important tools for, e.g., galactic archaeology, in particular for dating stellar populations through the use of “white-dwarf cosmochronology” (Fontaine et al., 2001), and their mergers are responsible for producing type Ia supernovae that are used as standard candles in Cosmology (Domínguez et al., 2001). While being of importance for several fields in Astrophysics, some relying on accurate white-dwarf models, the internal structure of these stars, in particular their chemical stratification, remains poorly known due to significant uncertainties in the physics involved in stellar modelling (Salaris 2009; Fields et al., 2016). Ways to derive constraints on white-dwarf internal structure, which in principle is within the reach of asteroseismology, are therefore highly desirable.

Inferring the properties of white dwarf stars from their pulsations was one of the very first application of asteroseismology to stars other than the Sun, driven at that time by the discoveries of a growing number of pulsators of the kind (see Fontaine and Brassard 2008 for a complete overview of the historical developments in this field). Pulsators in this evolutionary stage are found in three main flavors as white dwarf stars cross the GW Vir (DOV), V777 Her (DBV), or ZZ Cet (DAV) instability strips when they cool down. Other groups of white-dwarf variables have also recently been identified, suspected, or predicted, adding further potential to an already promising situation. The oscillation modes detected in these stars are low degree (typically ℓ = 1 and 2), low-to-mid radial order g-modes with periods of a few minutes to about half an hour. Their periods are sensitive to the global stellar structure (and stellar fundamental parameters), internal chemical stratification, and rotation, as demonstrated by thorough explorations of the pulsation properties of DA (Brassard et al., 1991, Brassard et al., 1992a, Brassard et al., 1992b, Brassard et al., 1992c), DB (Bradley et al., 1993), and GW Vir (Kawaler and Bradley 1994) stars. This led to the first quantitative asteroseismic ¹ inferences of white-dwarf properties, starting with the analysis of the hot GW Vir star PG 1159–035 (Winget et al., 1991), followed by numerous attempts. However, it is mostly from the year 2000 and beyond that the technique improved, triggered in part by the increasing power of computers allowing for more systematic forward modeling approaches (e.g., Metcalfe et al., 2000; Metcalfe 2001). Fontaine and Brassard (2008) and Giammichele et al. (2017a) provide very complete reviews and critical discussions–which we do not repeat here–of the various quantitative asteroseismic determinations that have been claimed for pulsating white dwarfs in the past 20 years. Most analyses were (and still very often are) focused on deriving global parameters and envelope thickness of the white dwarf, but one must acknowledge early efforts to probe the core, in particular to measure the ¹² C(α,γ)¹⁶ O nuclear reaction rate from asteroseismology (Metcalfe et al., 2001), a result that we however considered as premature due to issues with the models used for the analysis (Fontaine and Brassard 2002).

In a nutshell, not all seismic analyses of white-dwarf pulsators rest on the same solid grounds, due to significant differences in the quality of the equilibrium models used for the task and due to how the optimal model is ultimately found. Readers new to the field should always keep a critical mind to this respect. Various methods have been used, ranging from basic trial-and-errors modeling to systematic searches in more-or-less constrained parameter spaces involving up to hundreds of thousands of model computations through grids or optimization algorithms (e.g., Pech et al., 2006; Castanheira and Kepler 2008; Giammichele et al., 2016). In the past 15 years, the field has been strongly influenced by the rising belief that using full evolutionary models (or internal C/O profiles predicted by such models) for white dwarf asteroseismology was the way to go (e.g., Córsico et al., 2013; Romero et al., 2013). This strategy has however significant drawbacks because it implicitly assumes that evolution structures are accurate enough for that purpose, a questionable assumption considering the already mentioned uncertainties in the physics involved. Recently, De Gerónimo et al. (2017) explored the effects on core chemical stratification of various known uncertainties that plague central and shell helium-burning phases preceding the white-dwarf stage, and indeed the impact on pulsation periods can be quite large (up to 20–30 s period shifts according to the authors). The magnitude of such variations naturally raises concerns about the reliability of any approach that considers fixed (i.e., assumed perfectly known) composition profiles in the white-dwarf core. More generally, evolution-based approaches for seismic modeling do not allow, by construction, possibilities to provide independent measurements of several aspects of the internal structure, contrary to more flexible models. An additional caveat of this method is that it imposes the use of relatively coarse grids to explore rather limited (low dimension) model parameter spaces, due to large computation times. Grid-based approaches also have inherent drawbacks due to grid-resolution limitations as discussed by Charpinet et al. (2015).

In the meantime, efforts carried out by our group to exploit white-dwarf asteroseismology have taken a different route, based mainly on static white-dwarf models (as opposed to structures derived from evolution). Model grids and, later on, grid-free optimization tools were developed to achieve optimal seismic fits of the observed period spectra. Early applications of these methods include a derivation of the envelope structure of the DBV star GD358 (Fontaine and Brassard 2002), probes of the internal rotation of white dwarfs (Charpinet et al., 2009; Fontaine et al., 2013), and estimations of their bulk core composition (Giammichele et al., 2013a; Giammichele et al., 2013b, Giammichele et al., 2014, Giammichele et al., 2015, Giammichele et al., 2016). More recently, after establishing the period sensitivity of some deeply confined g-modes to the inner core chemical stratification that provide a chance to infer independent constraints on it from asteroseismology, we engaged into a significant extension of our seismic analyses in order to incorporate the detailed core structure in the fitting process. These developments were reported in detail in Giammichele et al. (2017a, 2017b, and 2018), Charpinet et al. (2019a), Charpinet et al. (2021). In the present paper, we review these latest developments (Section 2), discuss available and new results of the method applied to a handful of white dwarf pulsators (Section 3), and discuss some implications for our understanding of stellar modeling in the late stages of evolution (Section 4).

The foundations of our method are by themselves not new, as they rely on the widely-used forward-modeling approach, which in asteroseismology consists of searching for the best possible match between a set of pulsations periods detected in the star under consideration and periods computed from stellar models. Our implementation of this strategy originated from the pioneering work of Brassard et al. (2001), carried out in the context of pulsating hot subdwarf stars, and has been considerably improved since then, as detailed in Charpinet et al. (2005); Van Grootel et al. (2013a). Its adaptation to white dwarf asteroseismology, following the initial work by Fontaine et al. (2001), led to notable results (e.g., Charpinet et al., 2009), but mostly boomed over the last decade with developments carried out by Giammichele et al. (2014, 2015, 2016, 2017a, 2017b, 2018) that we summarize in subsequent subsections.

The technique relies on a double-optimization scheme that incorporates, along with the search in model parameter space, a nested minimization step that simultaneously matches observed and modeled periods (see Figure 1). This first, combinatorial optimization is performed each time a model of given parameters is evaluated, ensuring that the best possible fit for that model–which may not be a good fit in absolute terms, if the model is not a good representation of the star–is found. This step is necessary because in most cases observed periods are not univocally identified to a single mode of given degree ℓ and radial order k. With this first optimization, a mode identification is therefore achieved based on its ability to globally best-match the period spectrum, and an often-lacking pre-identification of the modes is no longer a requirement. Note that in cases where a mode identification is available a priori (from the observation of rotational splittings, through the analysis of period spacings, or from other independent methods such as multicolour photometry), these constraints can easily be incorporated in the combinatorial optimization process. The quality of the period match is evaluated quantitatively through a so-called merit function. This function could take many different forms (as long as it can be minimized or maximized), but using a formulation that is conveniently related to χ ²-statistics (for error estimates) has been our preferred choice, to date. The merit function, S ², that we use in the white-dwarf context is defined in its most general form by the relation S 2 = ∑ i = 1 N P obs , i − P th , i σ i 2 , (1) where N is the number of detected modes, P _th,i is the model period matched to the observed period P _obs,i , and σ _i is the error associated to the measured period. In practice, model uncertainties dominate by orders of magnitude over the usually tiny observational errors associated to period measurements (in particular with current long-baseline data from space) and in most cases we simply disregard σ _i (i.e., we set σ _i = 1) in the definition of the merit function. Again, we stress that the first optimisation step is a combinatorial minimization of S ² amongst all possible (P _obs,i , P _th,i )-pairs given a set of observed periods and a seismic model. The obtained minimized value, S ̄ 2 , is then considered as a function of the {a _n } model parameters that need to be minimized as part of the second optimization step. Once converged to the best-fit solution, i.e. the global minimum of S ̄ 2 ( a 1 , … , a n ) , we normally obtain the optimal seismic solution that best-reproduces the observed periods, which consists of the structure and fundamental parameters of the best-fit stellar model and the corresponding (ℓ, k)-identification of the modes.

An essential component of our approach to quantitative asteroseismology is the estimation of errors associated to the inferred properties. This requires that not only the optimal model shall be found, but also that its surroundings in parameter space be thoroughly explored in order to assess its statistical significance. Brassard et al. (2001) originally made these assessments from approximating the region localized around the S ̄ 2 -minimum by quadratic forms and applying criteria of χ ²-statistics to set confidence levels. Starting from Van Grootel et al. (2013a), we adopted a more sophisticated approach based on evaluating density of probability distributions from the likelihood function L a 1 , … , a n ∝ exp − 1 2 S ̄ 2 . (2)

The density of probability for each parameter is then obtained through marginalization P a i d a i ∝ d a i ∫ L a 1 , … , a n d a 1 … d a i − 1 d a i + 1 … d a n with ∫ P a i d a i = 1 . (3)

Both the search for an optimal model and the sampling of the Likelihood function to estimate errors are very computationally intensive tasks that require the evaluation of many seismic models. Efficient strategies to solve this generally difficult optimization problem are therefore needed. There are two factors on which we can act to mitigate this difficulty. The first one is by using stellar models that are fast to calculate and which permits a wide exploration of the model parameter space. The second is by adopting clever optimization techniques to reduce the total number of model evaluations. Disregarding these technical constraints usually leads to oversimplifying artificially the problem, e.g., by limiting the number of free parameters only to keep computation times acceptable and therefore implicitly disregarding entire sets of potentially valid solutions.

An important aspect of our strategy is to rely on full (but static, as opposed to evolutionary) white dwarf models. This approach has two main advantages: First, various parametrization of the internal structure can more easily be introduced, thus permitting more systematic explorations of a wide variety of stellar configurations. Second, computation of such models is fast (of the order of seconds on modern hardware, typically). These are key assets when pursuing the forward modeling optimization route described previously. Beyond these technical justifications, it is also worth recalling that timescales for oscillation waves to travel across the star (roughly the pulsation periods, which are of the order of minutes in white dwarfs) are much shorter than any timescale related to evolution (cooling) or mixing processes (e.g., microscopic diffusion; but with the notable exception of convective mixing that can have turnover timescales similar to the oscillation periods; Montgomery 2005). This fact implies that for the propagating waves, the background in which they travel is effectively static to a very high level of accuracy. The pulsations mode properties reflect this static configuration, but are not directly sensitive to the mechanisms that produce the structures during evolution. For instance, they can “see” a chemical transition of some shape at, e.g., the envelope or core boundary, since they are sensitive to μ-gradients, whatever the process that caused this chemical stratification. Pulsation modes are “agnostic” instruments allowing us to measure the structure of a star at a given time of its evolution, considering that by themselves they are fully independent of the evolution processes taking place on much longer timescales. The use of static models for asteroseismology (of white-dwarfs in particular, but this applies to many other types of stars), is therefore well justified in this context.

Realistic static models of white dwarfs can be computed from the fact that in these cooling objects the approximate relationship between the local luminosity and integrated mass, L(r) ∝ M(r), is quite accurate and, consequently, may be used to obtain an excellent estimate of the luminosity profile. This is true in particular for the lowest temperature white dwarfs including the pulsating ZZ Ceti stars. Note that Timmes et al. (2018) pointed out that this relation is no longer satisfied for the hottest white dwarfs (in particular the hot DBVs) which still have significant cooling occurring from neutrino emission. However, ways exist to overcome such complications by adapting the models and optimization strategy (see Charpinet et al., 2019a for further details concerning this specific case).

In order to define a full static white dwarf model, one must specify two fundamental parameters, the surface gravity log  g (or equivalently the mass, via the mass-radius relation) and the effective temperature T _eff, as well as a set of structural parameters that controls the envelope layering and the core stratification. Thus far, we have restricted ourselves to the main chemical constituents forming typical white dwarfs, namely oxygen, carbon, helium, and hydrogen. Massive oxygen-neon white dwarfs are not currently covered by our different flavors of static models, but there is no hard point, in principle, preventing us to do so in the future. Additional parameters associated to components of the constitutive physics could also be considered too, such as the convective efficiency via one or another mixing-length theory prescription. However, Fontaine and Brassard (2008) showed that the periods of low-order g-modes in pulsating white-dwarf models practically do not depend on the assumed convective efficiency. This parameter can therefore be fixed (to the so-called ML2/α = 1 version calibrated by Van Grootel et al., 2013b, in our case).

The envelope layering is generally dealt with a set of parameters–D(H), D(He), Pf _H, and Pf _He–that indicates the position (D-parameters) and shape (Pf-parameters) of the H and He composition profiles at the transitions. These parameters fully define a sigmoïd function inspired by analytical expressions derived when diffusive equilibrium is reached. A more detailed description of this method is provided by Giammichele et al. (2017a). For the specific case of DB pulsators, which are thought to have experienced C/He mixing during a late thermal pulse before the pre-white dwarf stage, which gravitational settling does not have time to erase, a double-layered envelope structure that can account for this possibility is more appropriate. This case is described in the Supplementary Informations of Giammichele et al. (2018).

Our treatment of the core is the part that has seen the most significant improvements recently. The most recent developments are described in the next subsection. We recall here that various assumptions can been considered when specifying the core structure in our static models. The simplest approach is to assume an homogeneous C/O core requiring only the mass fraction of C (or O) to be specified. These simple models were used in all of our asteroseismic studies of white dwarf pulsators up to (and including) the analyses of Giammichele et al. (2016) that helped determine the bulk composition in some white dwarf stars. This approach is however crude and does not fully exploit the asteroseismic potential of white dwarf pulsators. Another implemented option is to incorporate a pre-calculated C/O profile that mimic the features of an evolutionary-based chemical profile. This more sophisticated treatment originally proposed by Metcalfe (2005) only employs two free parameters, the constant central abundance of oxygen and the fractional mass at which the constant oxygen mass fraction starts dropping. Everything else is predetermined by scaling the shape of the drop from an evolutionary track at the specific white dwarf mass of 0.61M _⊙ computed by Salaris et al. (1997). Over the years, these “physically motivated” C/O profiles were used for the study of numerous DAV and DBV by Bischoff-Kim and Metcalfe (2011) and Bischoff-Kim et al. (2014, 2019). While this may sound as an attractive alternative (in line with other studies based on evolution models), the lack of flexibility in defining the core configuration becomes a major limitation. Quite notably, this parametrization fails at reproducing the core composition profiles from evolutionary models at other masses and is unable to reproduce the triple chemical transition at the edge of the core, where oxygen, carbon, and helium are expected to coexist, according to current evolutionary sequences (see, e.g., Althaus et al., 2010).

Finally, once all components (parameters) defining the white dwarf structure have been specified, stellar structure equations are solved to determine the corresponding equilibrium state. This is done with the Montreal stellar structure and evolution code STELUM (Bédard et al., 2021; see also Section 4.2), which is the most recent evolution of our white-dwarf and hot-subdwarf modeling tool. The code is based on a Galerkin finite element solver and implements current state-of-the-art physics (EOS, opacities) to describe the star interior. Importantly, special care is taken by the code to produce models with smooth profiles for all physical quantities. Profile smoothness, which may not be as important in a stellar evolution context, is absolutely key for calculating pulsation periods accurately (see Brassard et al., 1991, Brassard et al., 1992a, Brassard et al., 1992b, Brassard et al., 1992c).

The novelty introduced by Giammichele et al. (2017a) in our static modeling capabilities has been to consider more realistic depth dependent core compositions, as obtained from detailed evolution calculations, without sacrificing flexibility. Determining the best (optimal) composition stratification given a set of observed pulsations periods sensitive to the core structure can be defined as an optimal design (or shape optimization) problem. These are encountered in many domains, such as aerodynamics, hydrodynamics, acoustics, electromagnetism, just to name a few. The typical problem is to find the shape that is optimal in that it minimizes some merit function, while satisfying given constraints. For example, in airfoil shape design, the problem is entirely defined by a set of parameters controlling the shape that are adjusted to minimize a merit function, typically the inverse of the lift-to-drag ratio. Our situation, the “design” of the chemical profile optimally matching the pulsation period spectrum through minimization of the S ̄ 2 -value defined previously, is analogous. It is important to realize that the number of parameters used to define the shape (either be it of an airfoil or of a white-dwarf core composition profile) has nothing to do with the number of quantities involved in the merit function. Every single shape of the chemical (or airfoil) profile leads to a single merit function value, independent of the adopted parametrization of the shape.

If the chosen parametrization is not flexible enough, meaning it cannot represent all possible shapes, then the solution will generally be sub-optimal, i.e., the true optimal shape can never be reached. Conversely, if the number of design parameters is too large, the underlying optimization problem becomes impractical due to prohibitive computation times. Various strategies, with their advantages and drawbacks, exists for the shape parametrization, involving either analytical, discrete, polynomial or spline representations (Gallart 2004). The most commonly used methods are based on spline representations which provide excellent compromises between flexibility and sparsity of the shape parameters. With this approach, a set of well-positioned weighting parameters, called control points, is used as the design variables and spline interpolation is applied to produce the profiles in between points. The great advantages of using splines are: it allows for a reduction of the number of control points (fewer parameters to optimize) while retaining a high level of flexibility (large variety of attainable shapes); the perturbation of one control point only has local effects on the overall shape (increased independence of the shape parameters), and it produces curves with C ²-continuity (guarantying a smooth shape without high frequency oscillations). The latter property, in particular, is key for our problem of asteroseismology because smoothness in the composition profiles is essential to avoid artificial perturbations of the pulsation period spectrum, which is very sensitive to steep gradients.

In our context, amongst the many possibilities, we adopt spline representations to reproduce the oxygen chemical profile in a white-dwarf core based either on Akima (Akima 1970) or Fritsch-Carlson (Fritsch and Carlson 1980) splines. Both are continuously differentiable sub-spline interpolations using piece-wise third-order polynomials. Only the data from the next neighbor points are used to determine the coefficients of the interpolation polynomial, improving locality. Fritsch-Carlson assumes in addition a monotonic (either always decreasing or increasing) behavior of the function, which is generally relevant in our case. Both schemes are roughly equivalent, but the Fritsch-Carlson splines tend to be smoother near sharp chemical transitions. Compared to other more classical splines, such as cubic splines, both Akima and Fritsch-Carlson splines are immune to oscillations (also called “Gibbs noise”) in the vicinity of outliers or when gradients change abruptly. Their key advantage is to provide a stable and robust way to construct a wide variety of chemical stratification in the core.

When adopting this prescription, eight parameters are necessary to fully define a chemical stratification in the core that can resemble that of a white dwarf, as evolution models predict. It is expected from such calculations that the oxygen profile decreases outward following a double-ramp structure (see, e.g., Salaris et al., 2010). From our perspective, this double-ramp profile is obtained (see Figure 1 of Giammichele et al., 2017a) by specifying the position t ₁ and t ₂ (in fractional mass depth log  q) of the two transitions, their respective half-width Δt ₁ and Δt ₂ (in log  q), and the amount of oxygen (in mass fraction) at key locations, that is core(O) at the center, t1(O) and t2(O) at the bottom of the t ₁ transition and top of the t ₂ transition, respectively, and envl(O) at the bottom of the t ₂ transition. This last parameter accounts for specific cases involving DB white dwarfs for which it is believed that the mass fraction of oxygen may not drop to zero at the edge of the C/O core because of some dredge-up that occurred in the pre-white dwarf stage. In the context of DA white dwarfs, the envl(O) parameter is simply set to zero and the problem reduces to a seven-parameter shape optimization.

Figure 2 illustrates the flexibility of this approach by closely reproducing a set of configurations obtained from detailed evolution calculations. The C, O, He profiles are those derived by De Gerónimo et al. (2017) from a series of evolutionary models assuming different input physics for the helium-burning core (i.e., changing nuclear rates, overshoot at the boundary of the convective core) and the subsequent helium-shell burning phase (modifying the number of thermal pulses). As mentioned previously, these illustrate the significant changes occurring in the white dwarf interior just by taking into account typical uncertainties associated with these processes. Such considerations apart, the goal of the present exercise was to mimic each of these configurations with our static models incorporating the new parametrized core structure. This is achieved with no difficulty in all cases. In a context of detailed asteroseismic studies, it means that all of these configurations would belong to the searched parameter space and one of them could therefore be isolated if it happens to be the optimal representation of the real star. However, the power of the method extends well beyond choosing amongst predicted evolution profiles. A wide variety of configurations, including many that are unreachable with current evolution models, can be generated through this technique and tested for their ability to match the observed pulsation properties. Hence, this flexibility provides the opportunity to let the pulsation modes themselves determine the internal composition profiles in a way that is mostly independent of stellar evolution.

The approach described previously, in particular the inclusion of parametrized shapes to construct the chemical profiles in the core, defines a multi-dimensional global optimization problem where the quantity S ̄ 2 has to be minimized in the model parameter space. This is quite difficult to perform with brute-force methods such as non-adaptive, uniform grid computations commonly used for white-dwarf asteroseismology. With such grids, the number of model evaluations would grow as N ^d , where N is the number of grid nodes and d the dimension (number of parameters) of the problem. Even a coarse grid made of only N = 10 grid-points per dimension for a problem limited to optimize only the core stratification (7 parameters) and the star fundamental parameters (2 parameters), thus assuming a fixed envelope structure, would require 10⁹ model computations. This is not feasible with contemporary computational resources (it would require more than ∼ 1300 years of CPU time for our example) and, as a matter of fact, grid-based approaches have been limited so-far to low-dimension (few parameter) problems ² . We recall that this limitation often imposes to over constrain the search (by adding assumptions that fix parts of the model structure) in studies that are based upon grids. Grids also have additional caveats that make them dangerous to use for optimization problems. Resolution is one of them. If the grid is too coarse for the problem investigated, the true global solution might be missed in favor of a different sub-optimal secondary minimum closer to a grid node. The results could therefore be misleading. This issue can be mitigated to some extent by using interpolation schemes between grid-points, as long as variations between adjacent models are smooth enough (which is not always the case, especially in evolved stars due to mode trapping effects). A critical discussion of the limitations and dangers of model grids for white-dwarf asteroseismology has been provided by Charpinet et al. (2015). This should be kept in mind as grids are still widely used in the field.

In order to address this complex global optimization problem, we developed a specific package that, amongst other tools, includes an hybrid real-coded genetic algorithm (RCGA; e.g., Eshelman and Schaffer 1993) multi-modal optimizer (a code named LUCY) whose role is to locate efficiently, in a single pass, the global optimum of the function to optimize, as well as eventual secondary optima. The basic principle of genetic algorithms (GA; Holland 1975; Goldberg 1989), which otherwise exist in many different flavors (see the recent review of Lambora et al., 2019), is to evolve a population of solutions, with each individual solution defined by a set of genes (coded from the parameters) and a merit value (evaluated from the function to optimize), through a set of selection, mating (crossover), and mutation (random sampling) operators. A selection pressure in favor of more fitted solutions, along with the mixing of their properties (mating) ensures the overall migration of the population toward best-fit regions of the parameter space, while randomness (through mutation) preserves diversity and the ability to explore all regions of the search domain. GA have proven to be highly effective to address hard optimization problems, while keeping the number of function evaluations relatively low. They are less effective in converging rapidly to the solution with precision, however, which is why hybridation (in our case by using the simplex method; Nelder and Mead 1965) is often considered in the final stage of the optimization run. The implementation in LUCY has several critical properties for the purposes of asteroseismology. First, the code uses real-coded genes (as opposed to more traditional binary-coded GAs) that allow to search parameter space continuously, i.e., there is no such thing as a resolution since the code generates solutions (models) with any value for their parameters in the specified range. Second, the code is multi-modal, by allowing subsets of its generated population to exist and evolve independently, and can therefore locate local optima at the same time as it robustly converges toward the global optimum. This capacity is fundamental to assess unicity of the solution, an aspect that is often overlooked in asteroseismic studies. Finally, the code is massively parallel and can perform hundreds-to-thousands of model evaluations simultaneously (when run on high-performance computing facilities), which is also a key property for the success of our approach. Just to give an idea of the computing power needed when white-dwarf asteroseismology is the objective, a typical optimization run with LUCY requires from several hundred thousand to nearly a million of seismic model evaluations (i.e., sequential operations as illustrated in Figure 1). This is many orders of magnitude less than an equivalent (and unfeasible) grid computation, but still a heavy calculation task that typically requires several hours using a few hundred CPU cores working in parallel on modern hardware. An application of the LUCY optimizer to a “noisy” 2D test function that would be challenging for traditional methods is presented in Figure 3, demonstrating how the code can locate multiple optima (including of course the global optimum) and perform some sampling in relevant regions of parameter space (near most significant optima) to estimate the local shape of the merit function. We see that the 3 lowest minima of interest are easily recovered, while other secondary deeps (up to the maximum number of subgroups of independent solutions that the code can sustain given the initial population size) are also spotted. If this function was a χ ²-type distribution, such as in a real asteroseismic optimization exercise, parameter marginalization would lead to the probability distribution of each parameter (and therefore the error estimate), as illustrated in the lower left panel (and see Giammichele et al., 2016).

Our approach to tackle quantitative asteroseismology (of white-dwarf stars, in particular) has therefore a number of key methodological advantages that we summarize once more here. Among these, objectivity is an essential requirement that we always seek to achieve. To be objective, the search for the optimal seismic model has to derive from a wide and thorough exploration of model parameter space (criterion 1), using a robust method for finding the best-fit solution(s) (criterion 2). Uniqueness of the optimal solution should be stated in light of the quality of fit achieved for the best solution and its significance when compared to secondary optima (criterion 3). Finally, consistency with other independent constraints (e.g., from spectroscopy or parallax measurements) must be checked as well (criterion 4). Several of these requirements, in particular criteria 1 and 3, are often overlooked in the literature. Romero et al. (2012) and Castanheira and Kepler (2008), for instance, due to the use of grids of evolutionary models, generally impose a highly constrained search often limited to 3 parameters only. Uniqueness of the solution, often related to the achieved quality of fit (in an absolute sense) are rarely evaluated in this context. Most methods also use rather coarse grids (with the inherent problems associated to grids mentioned previously) using a relatively small number of model evaluations. If a high number of model calculations is not by itself a guarantee to find the global minimum and assess its unicity, it is necessary to ensure this last criterion is met. It should be obvious that a global search for best-fit models in a wide parameter space, owing to also assess errors and uniqueness, cannot be achieved with only a few hundred of seismic model calculations. Hundred of thousands of such models are actually needed to effectively sample the function that is optimized. The numerical tools that we developed are designed to give us the means to reach all of these criteria of objectivity. They do not rely on model grids and are computationally efficient to allow for a detailed exploration of parameter space.

This approach still has some limitations which one should also keep in mind. In particular it is not fully model independent. Whilst dependency on stellar evolution has been highly mitigated through the introduction of flexible static models, the latter still depend on constitutive ingredients, such as radiative and conductive opacities, equation of state, convection, just to name a few, which all have their underlying uncertainties. It should be clear, however, that the same limitations apply to all methods found in the literature, in particular those based on full evolutionary models which make use of the same ingredients for the input physics. Because of this, a seismic solution can still be biased, in the sense that whatever the precision achieved for the period fit and parameter determination, the solution may not be accurate, i.e., prone to systematics. Quality of fit (in an absolute sense) and consistency with independent constraints (criterion 4) become particularly important in this context, as they offer the only way to evaluate solution accuracy. There exist other sources of systematics linked to the numerical schemes themselves, such as how composition profiles are approximated with splines for instance, or to various assumptions entering the definition of the static models. Neglecting or not neutrino cooling or the presence of minor species is one of them, but other sources of (not necessarily known) uncertainties in the constitutive ingredients of the models could also influence the solution. These can, to some extent, be treated when evaluating realistic errors (as opposed to formal errors of the fit that can be very small) to the derived stellar parameters. This issue will be discussed further in Section 3 (see also the various discussions provided in Giammichele et al., 2018; Charpinet et al., 2019a).

Using the new seismic DA white dwarf models, we applied our optimization technique to 4 ZZ Ceti stars that we estimated were promising candidates: KIC 11911480 (Greiss et al., 2014), L19-2 (Sullivan and Sullivan 2000), SDSS J1136 + 0409 (Pyrzas et al., 2015), and EPIC 220347759. These stars were selected on the basis of the simplicity of their pulsation properties, a general characteristic of stars located near the blue edge of the ZZ Ceti instability strip, and exhibit as such pulsation periods in the range 100–500 s. These stars have high quality photometric lightcurves obtained from space with the Kepler spacecraft, except for L19-2 whose analysis relies on data from a former ground based multi-site campaign (Sullivan and Sullivan 2000). We present in this work, a selection of the results obtained from the seismic fits performed for these 4 DAV stars, while the details (fit analysis, errors assessments, comparison with other works, comparison with independent observational constraints, and analysis of rotation) will be disclosed elsewhere (Giammichele et al., 2022; in preparation).

Five independent periods were isolated for the seismic inferences of KIC 11911480 and L19-2, while SDSS J1136 + 0409, and EPIC 220347759 were analyzed based on eight and seven independent modes, respectively. These stars being close to the blue edge of the DAV instability strip, show few, but sufficient, low-degree and low-order g-modes, to obtain non-ambiguous seismic solutions in all cases.

We find that the seismic fits uncovered reproduce very closely the measured frequencies for all 4 stars, although less so for SDSS J1136 + 0409, as illustrated in Figure 4. While still reaching the precision of observations, the quality of fit for SDSS J1136 + 0409 is significantly lower than those obtained for the other stars. This case shows that our approach is not by nature an under-constrained problem that could match perfectly any period spectra. Another important point to highlight is that matching the periods at this precision with a model just reflects the fact that internal errors of the fit are very small, and that we are not limited, for example, by the flexibility of the parametrization or by insufficient grid resolution. Of course it should be clear that an almost perfect fit does not necessarily imply that the parameters inferred from it are perfectly constrained, for the reason that additional sources of uncertainty or systematics must also be taken into account (see below).

Our seismic estimates of the most relevant parameters derived for the 4 stars considered are provided in Table 1. In all cases but L19-2, we find seismic solutions that show a very good consistency with independent measurements of T _eff and log  g from spectroscopy. The latter were obtained from the most recent grids of DA model atmospheres (Bergeron et al., 1995; Gianninas et al., 2011) and by applying the 1D/3D correction from Tremblay et al. (2013). The apparent discrepancy for L19-2 is due to the fact that the optimal seismic model rather belongs to a family of solutions defining a long valley mainly oriented along the effective temperature axis, meaning that the observed pulsation modes in L19-2 do not constrain very tightly that parameter. It is important to note that, as was mentioned earlier, an almost perfect fit as obtained for L19-2 does not prevent possible degeneracies in some parameters of the solution.

The seismic solutions uncovered lead to a complete mapping of the chemical stratification inside these stars. Some of the values defining these composition profiles are given in Table 1, including the estimated masses for the hydrogen-rich envelope (log  q(H) = log[1 − M(H)/M _*]), the helium mantle, log  q(He), and the central homogeneous C/O core, log  q(core). We also indicate the oxygen mass fraction in that central region, O(core). In addition, we show in Figures 5–8 the derived composition profiles (with estimated uncertainties) for these objects. As a general rule of thumb, we typically find that the best determined profile parameters are the positions of the main chemical transitions–which is not surprising as they have a direct and usually strong impact on the period spectrum through mode trapping effects –, followed by central oxygen mass fraction in the core, which is also rather tightly estimated. In contrast, the inter-shell mixed C-O-He regions are usually the least constrained ones (larger uncertainties), with some variations from one target to the other (see also Charpinet et al., 2019a).

Overall, our results point towards rather thick hydrogen envelops and helium mantles, by mass. Remarkably, the derived central core extent is found to be quite similar for the four stars, ranging from log  q(core) ∼ − 0.63 to ∼ − 0.82. We emphasize that a similar core size has also been obtained from the seismic analysis of the DB white dwarf pulsator KIC 08626021 for which log  q(core) = −0.72 ± 0.03 (Giammichele et al., 2018). These values lead to a mass for the central homogeneous C/O core–mainly produced by the former helium burning core–that is ∼ 40 % larger than predicted by standard stellar evolution calculations. This striking consistency obtained from our analyses of 5 different white dwarfs, so far, strongly strengthen the conclusion that the helium burning phase generally produces larger cores than currently expected. Similarly, we find a convergence from all the white dwarfs we have analyzed toward a larger (by ∼ 15 % ) mass fraction of oxygen produced in the core compared to current model predictions. This will be further discussed in Section 4, below.

Quantifying proper uncertainties is a mandatory exercise in evaluating what information can be extracted from the seismic solutions. Giammichele et al. (2018), and later on Charpinet et al. (2019a), attempted to evaluate additional sources of uncertainty (systematics or external errors), that could overcome the (often very small) internal errors estimated from the period fit propagated to the derived parameters. Internal chemical profiles and global parameters inferred do, to some extent, depend specifically on the particular components of the constitutive physics that went into building any model. That is, the results must be sensitive to, for example, our choice of radiative opacity, conductive opacity, or equation of state for the fully ionized interior as well as the partially ionized non-ideal envelope.

Beside those basic physics constituents, many other choices come down in the building of one single model. In the evolutionary approach, each WD chemical structure is the result of stacking numerous physical processes through the evolution of progenitor stars from the ZAMS to the thermally pulsing and mass-loss phases on the AGB, and from time-dependent element diffusion in the WD stage. Accepted standard prescriptions are routinely used without proper evaluation of systematics associated to them, and how those uncertainties can propagate to the final stages of evolution. Uncertain convective boundary mixing during the He-core burning and the thermally pulsating AGB phase (e.g., Constantino et al., 2015), coupled with the occurrence of possible extra-mixing episodes generate major variations on the extent of the convective boundary, and therefore on the resulting WD chemical transitions. Efforts to secure an Initial-Final Mass Relationship (IFMR) for stellar evolution models are weakened by the lack of proper prescription for convective boundary mixing, uncertainties on mass loss rates and effects of metallicity or He content dependence (Salaris et al., 2009).

Efforts have been made in the past few years to quantifying some uncertainties inherent to WD progenitor evolution and assess their impact on parameters of the resulting model (De Gerónimo et al., 2017, De Gerónimo et al., 2019; Pepper et al., 2022). These authors have tested the impact of the occurrence of thermal pulses (TPs) in WD progenitors, the uncertainty in the ¹² C(α,γ)¹⁶ O cross section, and the occurrence of extra mixing during core He burning. For example, they found that the ¹² C(α,γ)¹⁶ O reaction rate, as one of the main sources of uncertainties, is alone, responsible for varying the central oxygen abundance from 26% up to 45%. They also estimated that comparing models derived from having experienced 0 or 10 TPs leads to average variations in the period spectrum of the order of 10–20 s, with absolute difference going up to 30 s. This is by no mean a small effect on the pulsation periods of white dwarf stars.

Such large period variations are obviously important uncertainties propagating to parameters obtained from optimal fits based on evolution models. In addition, we recall that other sources of systematic errors are still left unaccounted for. The He content of the final WD model, for instance, can vary by a factor up to 3–4 (Romero et al., 2012), but is usually assumed fixed in model grids used for asteroseismology. The thickness of the H layer (the principal inference that comes out from seismic analyses built upon evolution models) is based on an artificial procedure, which does not reconcile with a consistent standard WD formation. The suppression of breathing pulses, occurring at the end of the core He burning, which we will go into more details in the next section, is also an overlooked model prescription that could be source of uncertainties. Numerical smoothness of the models is yet another critical point for accurate pulsation period calculations. Worth noting in this context, Salaris et al. (2013) evaluated that the sole effect of different numerical implementations of the stellar evolution equations, while having rigorously the same physics modules input, can lead to a difference of 2–8% in WD cooling times, hence without even taking into account different physics prescriptions, nor the probable accumulation of numerical noise over the complete stellar evolution.

Figure 9 illustrates the impact of various global (non-adiabatic effects, presence of ²²Ne traces, impact of neutrino cooling) and localized (in that they affect the core stratification; e.g., a varying number of TPs, uncertain rate for the ¹² C(α,γ)¹⁶ O reaction, extra mixing) sources of uncertainties on low radial-order g-mode pulsation spectra representative of a pulsating white dwarf near the blue edge (i.e., showing low radial-order g-modes). We emphasize that in this comparison, static models were used to provide a strictly identical structure (and parameters), except for the part of the model that is modified to reproduce the effect being investigated. For instance, the model corresponding to an evolution that experienced no TP is exactly the same as the model corresponding to 10 TPs, except for the internal chemical profiles which were reconstructed from the evolution models shown in Figure 2. This ensures that only the signature coming from the effect under consideration is measured, a disentanglement that is usually difficult to achieve when comparing evolution models. From Figure 9, we clearly see how the global effects considered tend to be minor in magnitude compared to other sources of uncertainties coming from the core structure generated by past evolution. It is also important to note that the periods are almost uniformly shifted in those cases, meaning that trapping/confinement properties of the modes are essentially preserved (we point out that this was explicitly verified by Timmes et al. (2018), in the context of the impact of neutrino cooling). A consequence is that these global sources of uncertainties are unlikely to significantly bias model quantities that can be derived from asteroseismology thanks to their sensitivity to mode trapping, namely the positions of steep chemical gradients inside the white-dwarf star (Charpinet et al., 2019a showed that these are indeed robust determinations). In contrast, processes affecting locally the core chemical profiles can greatly modify periods on an individual basis. It is clear from just a glimpse at Figure 9 that the period distribution itself is strongly affected, meaning that the trapping/confinement status of the modes are changed. Therefore, we warn that adopting a fixed, potentially inaccurate chemical stratification for the core structure could be disruptive regarding the reliability of the seismic solution uncovered.

Table 2 is an effort to gather and sum-up all possible contributions that could contribute to the external error budget when using our static seismic models with the parameterized cores or when using approaches based on evolutionary models. The two approaches are not equivalent in that respect, as the use of static models can in many occasions mitigate the uncertainties associated to a given process by an increased flexibility on how the model is designed. An example of this comes from the parametrization of the core chemical stratification which, in our static models, can handle any profile generated by the processes shaping the core of a white dwarf, including those deriving from the uncertainties associated to such processes. Deviations from the expected profiles (i.e., the uncertainties in the processes that produce them) are already part of the search space for a best-fit seismic solution in that case. From the table, we clearly see that the more physical processes we add to the making of a pulsating WD model, the more uncertainties pile up. What is especially important to realize is that many of the processes involved in the making of a WD model through an entire evolutionary loop are not thoroughly asserted yet. And for those already evaluated, the impact on the error budget is already an order of magnitude greater than the remaining external uncertainties estimated for the static seismic models.

Beside a proper assessment of uncertainties, another essential aspect when performing a seismic analysis, is the criterion used to quantify the goodness of the fit. Recently, we have seen number of studies that rely on the Bayes Information Criterion (BIC) to evaluate and compare quality of different fits, which has the advantage of taking into account different numbers of data points and free parameters. As we mentioned before, finding an optimal fit does not guarantee that it is a good or even an acceptable fit. The same reasoning goes into comparing different fits with a BIC value. There is absolutely no guarantee that a better BIC value will yield a better fit (i.e., if the comparison model is a bad fit to the periods to start with, having a better BIC value does not guarantee that the new model indeed provides a good fit, merely that the fit is better). We therefore advocate here for the use of more robust criteria and more thorough analyses to ensure the conclusiveness of the optimal solution and quality of the parameters extracted.

As an example of why this is a critical point, we illustrate in Figure 10 the different orders of magnitude achieved in terms of quality of fit (here quantitatively evaluated as the mean and maximal differences achieved between the matched observed and computed periods) for typical DAV pulsators at the blue-edge using both methods. We show in the process some relevant precision limits that those fits should be compared with. First of all, we argue that caution should be used when the quality of fit reaches close to half the mean period spacing of the ℓ = 1 and 2 g-mode series, as it is the case for many DAV analyses performed with evolutionary models. This is an indication that the fit might be inconclusive to properly identify modes and derive reliable structural parameters, because the “best” seismic solution then becomes dangerously close to the worst cases that might be achieved when matching the periods. In addition to that, the impact on pulsation periods of the uncertainties associated with the core stratification (see again Figure 9), which would normally add up to the quality of fit, is also approaching the limit of half the value of the mean period spacing. To this account, we find Figure 10 to be particularly suggestive that what currently limits evolutionary models to achieve seismic solutions with a better quality of fit is the adoption of a core chemical structure that derives from evolution without considering its associated uncertainties. On average, the quality of fit obtained from evolutionary models remain close to two or three orders of magnitude larger than the actual precision of the measured periods, while static models provide a significant improvement in terms of quality of fit to that respect.

Our most constant and robust finding from the analyses of KIC08626021 by Giammichele et al. (2018) and of the 4 DA pulsators reported in Section 3, is that the central homogeneous core in white-dwarf stars is significantly more massive (by ∼ 40 % , with − log[1 − M _core/M _*] ∼ 0.6–0.8) and more oxygen-rich (by ∼ 15 % , with X(O) ∼ 0.80–0.90) than predicted by canonical evolution models (Salaris et al., 2010). This finding is highlighted in Figure 11, where we gather the core oxygen profiles derived for these 5 white-dwarf stars from our seismic modeling.

This has implications at various levels, starting with our basic understanding of the core helium burning phase that shape these layers. We discuss specifically this aspect in the next Section 4.2. It also has impact on the estimation of white-dwarf cooling ages which depend on the internal C/O composition ratio. To this respect, our seismic estimates provide a powerful calibration of the internal composition profiles that is of direct benefit for white-dwarf cosmochronology (Fontaine et al., 2001). The total C/O ratio–that we find lower than expected–also has an impact for the study of type Ia supernovae (Domínguez et al., 2001). Less ⁵⁶Ni is produced during the supernova explosion when this ratio is lower, with a direct influence on the light curve of type Ia supernovae which is governed by the decay of ⁵⁶Ni. These supernovae are often used as standard candles in Cosmology, notably for testing the cosmological equation of state which might then benefit from a more accurate seismic calibration of white-dwarf core composition. Finally, we point out that our findings should also impact the field of white-dwarf asteroseismology itself, when methods based on full evolutionary models or models incorporating profiles from evolution calculations are used to estimate fundamental and structural parameters of white dwarf pulsators. The core composition and structure has a significant impact–in fact much larger than other effects of lesser importance, such as neutrino cooling (Charpinet et al., 2019a) or the presence of minor species like ²²Ne–on the g-mode period spectrum. Therefore, assuming a canonical core structure derived by evolution calculations is likely to introduce a significant bias in the obtained seismic solutions.

The internal chemical stratification inferred from our white dwarf seismic models naturally raise the question of how such profiles can be produced by stellar evolution processes. Significant differences exist between these profiles and those predicted by standard evolutionary models and reconciling the two may indeed be challenging, possibly requiring deep introspection into our current understanding and modelling of stellar interiors and evolution. De Gerónimo et al. (2019) has been among the first to explore this issue by attempting to reproduce the seismic model uncovered by Giammichele et al. (2018) for the DBV pulsator KIC08626021. Exploring with some of the constitutive physics making up their evolution models (mostly varying parameters and values of already well identified and implemented processes), they could not easily reach a structure similar to the one found by Giammichele et al. (2018) unless–they find–unrealistic adjustments are imposed to the constitutive physics. Their underlying conclusion was therefore that it is not possible to reproduce the seismic model of KIC08626021, thus casting doubts about its validity. We have, however, to point out that De Gerónimo et al. (2019) did not consider all possible reasons, even amongst already known phenomena, that might have led to the differences suggested by asteroseismology.

As an illustration of this, we hereafter focus specifically on the problem of producing more oxygen at the white dwarf centre and overall larger cores than predicted by evolution. On that matter, De Gerónimo et al. (2019) reported that the rate of the ¹² C(α,γ)¹⁶ O reaction would have to be increased by at least a factor 100 to reach the mass fraction of oxygen found by Giammichele et al. (2018), which by far sounds unrealistic. However, raising the rate of nuclear reactions is not the only lever that can affect the amount of oxygen ultimately produced in the helium-burning core. In order to illustrate this point, we present a full evolutionary sequence computed with our most recent version of the STELUM package ³ , our specialized stellar structure code used for asteroseismology (among other applications) and evolution of evolved compact stars. STELUM cannot, in its current state, evolve stars through the Asymptotic Giant Branch (AGB; i.e., the usual path followed by low-to-intermediate mass stars prior the white dwarf stage), but it can produce Extreme Horizontal Branch (EHB) models and evolve them through the helium-core and helium-shell burning phases, and then let them contract and cool on the white-dwarf sequence. Therefore, if STELUM cannot currently account for all the processes that ultimately shape the white-dwarf internal chemical stratification (those occurring during the AGB), it is useful to explore the evolution of helium burning cores that determine white dwarf most central layers, considering that these cores are expected to be similar, in structure and evolution, for all horizontal branch stars.

A full description of STELUM is beyond the scope of the present paper. We refer the interested reader to Bédard et al. (2021; see also Brassard and Fontaine 2015) for a detailed presentation of this numerical instrument and the historical background behind it. Yet, it is worth mentioning that one of the most important and defining particularity of STELUM is that it uses the so-called Galerkin finite-element scheme to solve the stellar structure equations, unlike most other stellar evolution codes that rely on finite-difference methods. Brassard et al. (1992b) showed, in the context of solving stellar pulsation equations, that Galerkin finite elements outperform finite-difference schemes in terms of stability, accuracy, and effectiveness. Another particularity of STELUM, which is most relevant in our present discussion, is that it treats element transport, including ordinary convection, convective overshoot, semi-convection, chemical and thermal diffusion, gravitational settling, thermohaline convection, stellar winds and external accretion as time-dependent diffusive processes that are all considered simultaneously. Special care is taken to resolve these processes over their highly varying timescales during evolution, meaning that substantial computational effort (and time) is invested by the code into computing many (often thousands of) intermediate models between each evolution time step. This implementation turns out to be important for following the evolution of convective helium burning cores. In this specific context, the mixing processes that we find to be essential are convection, convective overshoot, and gravitational settling. Convection is treated in the framework of the mixing length theory using our so-called ML1/α = 1 flavor (Tassoul et al., 1990) for the sequence described below. Convective overshoot at the edge of the convective core is modeled as suggested by Freytag et al. (1996), using an exponentially decaying diffusion coefficient. In the calculated evolution sequence, the overshoot parameter was set to its “standard”, 3D simulation-calibrated value f _ov = 0.02 (Herwig et al., 1997). Gravitational settling and diffusion in general (including thermal and chemical diffusion) are treated with the formalism of Burgers (1969) using diffusion coefficients evaluated from the results of Fontaine et al. (2015) based on an improved version of the method described in Paquette et al. (1986). Finally, the nuclear reaction rates adopted for these calculations are those of Angulo et al. (1999).

Figure 12 shows, in the log  g − T _eff plane, our reference STELUM evolutionary track computed from an initial structure of 0.48 M _⊙ on the Zero Age Extreme Horizontal Branch (ZAEHB) to a T _eff ∼ 8500 K white dwarf. The model evolves through the typical core helium-burning phase (EHB)–the main focus of the present discussion –, followed by the shell helium-burning (post-EHB), pre-white dwarf and white dwarf stages. Figure 13 shows the evolution of the central composition and core size from ZAEHB until full exhaustion of helium at the center of the star. Figure 14 illustrates how the core helium, carbon, and oxygen mass fraction profiles evolve during core helium burning (left and middle panels), and during the post-EHB and white dwarf stages (right panel).

During central helium burning, the core experience three main phases. In a first early stage that last for ∼ 25 − 30 Myr, the core (defined as the region of mixed C/O/He composition underneath the nearly pure helium mantle) is fully convective and grows with time due to the production of carbon (more opaque than helium) transported above the purely convective layers by overshoot. This well-known outward migration proceeds smoothly until injecting carbon above the Schwarzschild limit becomes strongly destabilizing. This point is reached when a minimum in the radiative temperature gradient coincide with the core boundary, which happens when the central helium mass fraction Y _c ∼ 0.7 in our model. Evolution models often diverge at that point, depending on the treatments implemented to handle convection and the core boundary (see, e.g., Paxton et al., 2013, Paxton et al., 2018, Paxton et al., 2019; Constantino et al., 2015; Li et al., 2018 and references therein). Some prescriptions allow for a second convection zone to develop above the central one while others introduce prescriptions to produce a stabilizing partial mixing (or “semi-convective”) zone above the fully convective core. In calculations carried out with STELUM, no explicit prescription for semi-convection or partial mixing is used. However, we find that a partial mixing zone does indeed naturally develop from our implementation, resulting primarily from the competition between convective mixing, overshooting, and gravitational settling above the fully convective core.

The second stage of the core evolution is therefore characterized in our model by the development of a partial mixing region that further extend the overall core size over time. During that period, occurring from ∼ 30 Myr to ∼ 173 Myr, a mild composition gradient builds up between the fully mixed convective core and the core-to-mantle sharp boundary. A closer look at the corresponding model structures indicate that this gradient tends to achieve convective neutrality in this region ⁴ . During that period, the ¹² C(α,γ)¹⁶ O reaction becomes more prominent and oxygen is produced at an increasing rate in the core. At roughly 160 Myr, oxygen becomes the dominant element (in mass fraction) in that region.

When Y _c reaches ∼ 0.05 (at an age of ∼ 172 Myr), the core enters a third phase characterized by another instability of its outer boundary. The convective core undergoes a series of short-lived (a few tens of thousand year long) episodes, referred to as “breathing pulses” during which its size grows significantly, transporting fresh helium from the mantle into the core. The existence of these breathing pulses in real stars has been a long debated issue which is not settled yet (Sweigart and Demarque 1973; Castellani et al., 1985; Caloi 1989; Caloi and Mazzitelli 1993; Salaris and Cassisi 2017; Li et al., 2018; Ostrowski et al., 2021; Li and Li 2021). In most evolution codes, in particular those that currently produce canonical white-dwarf models, actions are usually taken to prevent this instability to occur and the pulses are suppressed (Caputo et al., 1989; Dorman and Rood 1993; Cassisi et al., 2001; Constantino et al., 2016, Constantino et al., 2017). In our STELUM models, the breathing pulses do not appear to have the characteristics of a numerical instability. They build up progressively and are well-resolved events. We therefore let them develop as a natural evolution of the stellar core, in particular since they have interesting implications. These episodes, that occur three times in our sequence, have two main effects on the star. First, each pulse pushes further out the boundary of the mixed C/O/He region, therefore increasing its overall size. Second, injection of helium into the core has a rejuvenating effect on the star which moves back to a (T _eff, log  g) position that it had earlier in its evolution (the visible loops in Figure 12). Consequently, the sequence reaches full helium exhaustion at the center after ∼ 220 Myr. This age for the whole core helium-burning phase is roughly 40 Myr (20%) larger than the age that would have been obtained without the breathing pulses. Another consequence of these pulses is that more oxygen is ultimately produced in the stellar core. Without the pulses, oxygen would top at a mass fraction of ∼ 0.7 at the end of core helium burning, the value that is typically obtained from standard evolution codes. With the breathing pulses, oxygen reaches a mass fraction of ∼ 0.86 at the time helium is exhausted. The latter value is remarkably consistent with the amount of oxygen estimated from our white-dwarf seismic models, a result that we interpret as evidence that breathing pulses are not model artefacts and that helium-burning cores do experience such events (we point out that Li and Li 2021 has reached similar conclusions). As a corollary, we suggest that arguments developed in the past to support the suppression of these pulses in evolution models–e.g., from the study of globular clusters (Caputo et al., 1989)–should be re-investigated.

In summary, it is not necessary to modify the ¹² C(α,γ)¹⁶ O reaction rate to obtain more oxygen in white-dwarf cores at levels comparable to the seismically determined ones. Breathing pulses during the late stage of core helium-burning evolution lead to such oxygen enrichment and appear to be the key process that is generally missing (or dismissed) in current stellar evolution models used to predict white-dwarf interior structures. More specifically, the breathing pulses were suppressed in the models discussed by De Gerónimo et al. (2019) and they were not considered as an alternative to solve some of the discrepancies these authors find between their evolution calculations and the seismic model obtained by Giammichele et al. (2018). We also find that the breathing pulses, along with the careful numerical treatment of mixing in the STELUM evolution code, contribute to produce cores that are larger than those obtained from canonical evolution calculations, thus closing or narrowing the gap between seismic and evolutionary models on that respect too. Overall, this demonstrates that reconciling our white-dwarf seismic models with evolution is not an impossible task, as claimed by De Gerónimo et al. (2019), but requires to review more critically how star interiors and evolution are modeled during the ultimate stages, keeping open the possibility that key events or processes may have been overlooked.

Stellar structures derived solely from observed oscillation frequencies using direct modeling or inversion techniques (that we commonly refer to as seismic models) have also been obtained for other types of stars. In particular, helioseismology has led to the construction of a seismic model of the Sun which shows puzzling differences with evolution calculations (see, e.g., Gough et al., 1996). An example of this is the sound speed profile at the base of the Sun’s convective envelope that appears to be shallower than expected from standard evolutionary models. While the seismic Sun still cannot be fully reproduced by evolution models, it is not an argument that is used to invalidate the seismic model itself, but rather to question the treatment of physics in the evolution calculations. White-dwarf seismic models, as the ones we presented in this paper and in Giammichele et al. (2018), have a vocation in all accounts similar to the establishment of the referential seismic Sun, with the added observation that uncertainties in the evolutionary path leading to white dwarfs are certainly far greater than the uncertainties in evolution that led to the Sun structure. We are therefore convinced that building seismic models of white dwarfs, and more generally of evolved compact stars (including hot-subdwarf pulsators; see Charpinet et al., 2019b and references therein), will be amongst the most prolific ways to shed light on the grey zones that still plague our understanding of the late stages of stellar evolution. For that purpose, the context is particularly favorable on the observation side, as large amounts of high-quality seismic data collected from space with instruments such as Kepler and the still operating TESS are available. We are committed to fully exploit these data in the upcoming years and refine our tools to do so.