WDEC: A Code for Modeling White Dwarf Structure and Pulsations

Modeling stars is complex, but perhaps the least difficult stars to model are the white dwarfs, at least in some mass and temperature range. They lack extended atmospheres (greatly simplifying energy transport) but are not so compact that their interior equations of state enter the ill-constrained realm. After an initial phase of contraction, white dwarfs reach a stage supported by electron degeneracy pressure, where the release of gravitational energy through contraction becomes a negligible contribution to the luminosity, the former being transferred internally to the free electrons. Nuclear fusion also becomes negligible, if at all present. The White Dwarf Evolution Code (WDEC) utilizes the fact that white dwarfs settle into a phase of stellar evolution where physical processes change on long timescales in order to compute white dwarf models using numerical methods that compute quickly. The average WDEC model runs in 10–15 s on a standard laptop computer.

The speed comes at a price, however. WDEC does not evolve the chemical profiles. They are given as an input and held fixed throughout the computation of a model. It has been argued that this is unphysical, and yet the output models are fully consistent solutions of the equations of stellar structure. WDEC will yield physical models if one feeds it physical chemical profiles. These can come from the output of stellar evolution models that do work out the time-dependent diffusion of elements.

It is worth discussing what the value of the WDEC is and how it has been used. As described in Section 2, the WDEC initially referred purely to a code that made models of white dwarfs. Early on, however, it was paired with a pulsation code and the two were used together to not only build white dwarf models but also calculate their oscillation modes. Today, the two codes are packaged together, and it is the package that we refer to as "WDEC" in this paper. As mentioned before, its main strength is its speed; relative ease of use is another advantage. This has allowed the use of WDEC in the asteroseismic fitting of pulsation spectra of white dwarfs.

Since their discovery (Landolt 1968), observations of pulsating white dwarfs underwent two revolutions that led to the wealth of pulsation data we have today and are still gathering. The first was the creation of the Whole Earth Telescope (Nather et al. 1990), which allows uninterrupted observations of pulsating white dwarfs over the course of weeks. The pulsation periods vary between 2 and 30 minutes, and multi-night observations are required in order to precisely measure the periods of oscillations. More recently, the space mission Kepler and its successor K2 have offered even more precise measurement of periods and revealed lower amplitude modes. Efforts have been ongoing to capitalize on that data and to try to infer what it tells us about the interior structure of white dwarf stars.

This is a computationally intensive process where individual models that compute quickly are an asset. As our goal is to unveil interior structure, the fact that the models do not compute the chemical profiles through the time-dependent diffusion of elements is not a setback. On the contrary, the philosophy is to assume as little prior knowledge as possible regarding interior structure (even though we do utilize what we know from stellar evolution calculations), and allow the observed periods to guide our determination of the interior structure. The efforts using WDEC started early on (e.g., Bradley & Winget 1994). Later they were pursued with the advent of faster computers (e.g., Metcalfe et al. 2000) and have continued since (e.g., Bischoff-Kim et al. 2014). In parallel, it is worth noting the recent efforts and success of a team in Canada and France who are using their own code but with the same approach (e.g., Giammichele et al. 2017).

We begin with a targeted history of the development of the code, to give some background to the current architecture of the code and also point to some key publications. In Section 3, we give an update on what physics are included in the code and on their implementation. In Section 4, we briefly discuss updates on pulsation calculations. We follow with a few validation tests of the code in Section 5, and with a summary in Section 6. We also point the reader to where to find the code and additional resources.

The history provided in this section is an important background to understanding the structure of the code and also some of its limitations. We do not discuss here the updates to the physics made throughout the years. For that, it is sufficient to present the physics currently included in the code (Section 3).

WDEC was originally put together by Don Lamb who fused together the evolutionary code developed over the years at the University of Rochester (Kutter & Savedoff 1969, and references therein) and the white dwarf envelope code written by Gilles Fontaine as part of a global effort led by their PhD supervisor Hugh Van Horn. The envelope code was "stitched" onto an evolving interior structure so as to provide a better description of the surface layers, including, in particular, the outer superficial convection zones that develop during white dwarf cooling.

A brief description of the envelope code is provided in Fontaine & van Horn (1976). An equally brief description of the resulting upgraded Rochester WDEC was provided by Lamb & van Horn (1975) including an interesting discussion of key processes occurring during white dwarf evolution such as neutrino cooling, convection, and crystallization. Lamb & van Horn (1975) presented the then most realistic evolutionary description of a white dwarf, using the example of a one-solar-mass pure C structure. The next major user of the code was Don Winget, who adapted it for his investigations of pulsating hydrogen atmosphere white dwarfs during his own PhD studies at the University of Rochester in the late seventies. The term WDEC was coined at that time.

Further details of the initial version of WDEC can be found in the PhD theses of Lamb (1974) and Fontaine (1974). Further developments and upgrades of the code migrated to the University of Texas at Austin. From that group, Wood (1992) gives a good overview of the code as it stood in the early 1990s.

WDEC was initially developed to run on the limited memory available on supercomputers in the 1970s. That meant that the models had a limited number of shells. That low resolution proved insufficient when the model output was used in the computation of p-mode non-radial oscillations expected in white dwarf stars. It also caused problems in non-adiabatic calculations through the lack of a sufficient number of shells in the driving/damping region.

The situation was remedied by the use of an intermediate code, called the "prep" code for a lack of a better word. The prep code added shells to the models through interpolation and also computed key quantities required by the pulsation code. In subsequent years, the stellar structure code computed models with increasing resolution, but the prep code retained its role. Output from the prep code would then be fed into the pulsation code to calculate periods of oscillation.

The process from providing input to obtaining pulsation periods required a series of manual steps. This became inadequate when the code started to be used in modern asteroseismic fitting involving the computation of thousands of models. As a part of his PhD dissertation, Travis Metcalfe re-packaged the different pieces of the code into a single Fortran function that could be called with a set of parameters as well as a list of periods to fit and that would return a single number, a goodness-of-fit parameter. That version of the code was folded into a Genetic Algorithm engine (GA) used to find the set of parameters that minimized the goodness-of-fit parameter, indicating that the periods of the model were a best match to the list of observed periods provided (Metcalfe & Charbonneau 2003).

When the present lead author inherited the code, she took a different approach. With observed period spectra that included an increasing number of modes, the routine that was matching the calculated periods to the observed periods turned out to be non-trivial to properly implement. Also, we wanted the ability to quickly match different sets of periods. In that respect, it appeared that saving the lists of periods produced in the process of minimizing the fitness parameter was a valuable thing to do. The grids of models could then be reused at will to match different sets of observed periods. The fitting routine, a possible source of fitting errors, became a separate process that could be better controlled. Today, WDEC does not include a matching routine. It simply makes white dwarf models and computes their pulsation periods.

As of this writing, the final round of key improvements to the code included rewriting the source code in modern Fortran and in a modular form, so that we could more easily integrate new physics into the code. Updates included interfacing WDEC with MESA (version 8118), Modules for Experiments In Stellar Astrophysics (Paxton et al. 2011). MESA is widely known in the stellar astrophysics community as a stellar evolution code. At its core, however, it remains a library of Fortran modules that can be integrated in codes to fit one's purpose. The new WDEC uses the MESA equation of states and opacities routines.

The core architecture of the evolution part of the code even today is still well described in the initial instrument paper (Lamb & van Horn 1975). In particular, the paper describes integration schemes used and boundary conditions. Much of the physics, however, has been updated over the years. Some of these changes were documented in publications that featured work using the code, but in the absence of a dedicated paper detailing these changes and in light of the most recent changes, it is worth gathering that information in the present work.

Numerically, the code treats two regions of the interior separately and then stitches them together (quite literally, the subroutine that does that is called "stitch"). The outer region is described by the envelope code referred to above. The location of the boundary between the two regions can be adapted to one's need. Generally, it is best to place it as close to the surface of the model as possible, as the envelope is not as complete as the core in terms of physics (for instance, it does not include neutrino cooling). A convenient and sufficient place to place the core-envelope boundary is .99 M_*.

We strive to make the models continuous across this boundary. MESA has refined the smooth transitioning between the different equation of state tables. The transition from core to envelope is continuous in that respect. The equation of state and opacity tables used in MESA are thoroughly described in Paxton et al. (2011). To help put things in perspective, we graph in Figure 1 the region the white dwarf models occupy in the ρ–T plane and indicate the relevant equation of state tables.

Zoom In Zoom Out Reset image size

We detail below what physics are included in each part of the models (core and envelope).

3.1. Neutrino Emission

We stress again that the term "core" refers to the region of the model below a purely numerical boundary. The core is the only region of the model where we treat neutrino emission. Neutrino emission (and any other energy generation) is set to zero in the envelope. For average-mass white dwarfs with effective temperature greater than 24000 K, neutrino emission dominates the cooling and is an important piece of the physics (Winget et al. 2004). In its current version, the code includes the prescriptions for neutrino emission through a variety of physical processes listed in Table 1. Of these processes, when neutrino emission is present, plasmon neutrino emission dominates (Figure 2).

Zoom In Zoom Out Reset image size

Table 1.  Neutrino Emission Processes

Process	References
Photoneutrinos	Itoh et al. (1989) (+ Errata Naoki et al. 1990)
Pair neutrinos	Itoh et al. (1989) (+ Errata Naoki et al. 1990)
Plasmon neutrinos	Itoh et al. (1996)
Recombination neutrinos	Kohyama et al. (1993)
Neutrino-pair Bremsstrahlung
Liquid Metal	Itoh & Kohyama (1983)
Low-temperature quantum corrections	Itoh et al. (1984a)
Crystal Lattice	Itoh et al. (1984c)
Phonon contributions	Itoh et al. (1984b)
Partially degenerate electrons	Munakata et al. (1987)

Download table as:  ASCIITypeset image

Another source of cooling included in the models is that hypothetically due to axions. The specific kind of axions considered are those that would be produced in electron bremsstrahlung events with the emission of an axion instead of a photon. That is the production process for axions one would expect in a white dwarf interior, where free electrons abound. The axion emission rates are implemented using routines from Nakagawa et al. (1988). By default, the axion mass, a free parameter, is set to zero (no axion cooling). The axion emission rates are dependent on the axion mass. Current constraints obtained by observing the effect of cooling on the pulsation periods of white dwarfs point to m_acos²θ < 30 meV (Bischoff-Kim et al. 2008; Córsico et al. 2016).

3.2. Treatment of Convection

The envelope is where convection zones can form in the white dwarf models. The treatment of convection is standard mixing length theory, following Bohm & Cassinelli (1971). The theory involves one parameter, α, tied to the convective efficiency. Convection is important below 30000 K for helium atmosphere white dwarfs and below 14000 K for hydrogen atmosphere white dwarfs (3). The former being the temperature range where helium in the atmosphere of the white dwarf is partially ionized and the latter where hydrogen is partially ionized.

For carbon and oxygen core white dwarfs, different calibration schemes can be used to assign a value to α and we implemented two different ones in the new version of WDEC. In recent years, it has become possible to model convection in three-dimensional (3D) hydrodynamical simulations. The results of the simulation can then be used to choose values of α that lead to convection zones of proper depth for one-dimensional (1D) models of any given effective temperature and surface gravity. Such calibration was carried out for hydrogen atmosphere white dwarfs (Tremblay et al. 2015).

Another method uses pulsating white dwarfs (Provencal et al. 2015). Pulsating white dwarfs that are on the hot end of the pulsating range (blue edge objects) are observed to produce light variations that are sinusoidal in nature. Cooler white dwarfs have light variations that strongly depart from sinusoidal curves. Such nonlinearities can be understood to be caused by deep convection zones. The thermal timescale at the base of their convection zones can be inferred from the shape of the pulses in the light curve. We have used these results to set the value of α in the code.

Finally, for numerical experiments, it is also possible to use α as a free parameter when running models. Larger values lead to deeper convection zones for any given atmospheric conditions in the model.

Figure 3 summarizes the values of α that get computed by the code using the two calibration schemes above. Note that we obtain different values of α depending on the calibration method used (with agreement for hydrogen atmosphere white dwarfs around 12000 K, a temperature at which a number of the pulsating white dwarfs are found). The differences can either be due to the differing methods, or discrepancies in the models, or both. It is worth discussing the methods used in the calibration further. We check for possible discrepancies in the models in Section 5.

Zoom In Zoom Out Reset image size

The α based on the 3D hydrodynamical simulations are calculated using fitting formulae provided in Tremblay et al. (2015). According to the authors, they are good to within 5%. One difficulty we had to overcome in the implementation of the formulae is the fact that they are given in terms of independent variables (T_eff, log g). The WDEC accepts the effective temperature and the stellar mass as input, the surface gravity being an output of the code. In order to obtain a log g to feed into the formulae at the beginning of the execution, we produced tables of surface gravities for models of differing temperatures and masses over a wide range of these parameters and produced fitting formulae to find the surface gravity based on these tables. This is only done for hydrogen atmosphere white dwarfs, as the Tremblay et al. simulations only apply to hydrogen atmospheres. The log g obtained are well within 1% of the values obtained for the model at the end of execution.

The α based on the modeling of the nonlinearities in the light curves is calculated by tuning α in the WDEC models until the thermal timescale a the base of the resulting convection zones matches the thermal timescales determined by fitting the shape of the light curves of pulsating white dwarfs (Provencal et al. 2015). The result of such fitting, based on 13 white dwarfs, is shown by filled circles in Figure 3. The dashed parabola is a quadratic fit to these points. The values of α are calculated using the equation of that parabola, for temperatures below 15000 K. If the value of alpha comes up negative, α is set to an arbitrarily low value. At those temperatures, the depth of the convection zone is not very sensitive to α (Figure 4), as convection becomes adiabatic (Tassoul et al. 1990). Above 15000 K, α is arbitrarily set to 1. At those higher temperatures, hydrogen is fully ionized and the convection zones disappear, for any reasonable value of α.

Zoom In Zoom Out Reset image size

Regardless of the method adopted to determine it, the chosen α is further scaled by a factor that is near unity in order to obtain convection zones that have the same depth as that determined by other codes, for a given α, log g, and T_eff. In effect, this allows WDEC to speak the same language as other codes when it comes to convection.

We apply a similar method for helium atmosphere white dwarfs. We list the helium atmosphere white dwarfs with "measured" timescales at the base of their convection zones and the corresponding α in Table 2 (Provencal et al. 2015). The data for these stars is more limited and noisier and a curve fit is not appropriate. Instead, we use the average of the different values of α listed in Table 2 (0.96).

Table 2.  Properties of Helium Atmosphere White Dwarfs Used in the Calibration of Convection

White Dwarf	Teff (K)	log g	Mass (M⊙)	τtha (s)	α
PG1115	25,000	7.91	0.561	580	1.20
GD358	24,000	7.80	0.505	580	0.96
EC04207	25,970	7.91	0.563	90	0.85
PG1351	26,000	7.91	0.563	90	0.85
EC20058	25,500	8.01	0.615	50	0.67
WDJ1929	30,000	7.89	0.563	10	1.24

Note.

^aThermal timescale at the base of the convection zone.

Download table as:  ASCIITypeset image

The two methods are based on different approaches that are of equal value and validity, and yet do not yield the same values for α. We leave it up to the user to adopt the method they feel more comfortable with, or to treat α as a free parameter. For helium core white dwarfs, we do not propose any calibration, and the user has to provide an input value for α.

3.3. Parameterization of Chemical Profiles

As mentioned in the introduction, WDEC does not calculate element abundances, nor does it carry out time-dependent diffusion calculations. The chemical profiles must be specified as an input and are held fixed throughout the computation of a model. One version of the code makes carbon/oxygen core white dwarf models, while another makes helium core white dwarfs.

3.3.1. Carbon/Oxygen Cores

In that version of the code, the user is called upon to furnish an oxygen abundance profile, the helium abundance in the region of the model transitioning from a mix of carbon and oxygen to pure helium, and if desired, where the transition from pure helium to pure hydrogen is to take place. There are also parameters associated with the shape of the transition zone from the C/O region to the pure helium region. The code determines the carbon abundance profile and the shape of the He/H transition zone. The latter is calculated assuming the helium and the hydrogen are in diffusive equilibrium (Arcoragi & Fontaine 1980), using the "exact" (non-trace element approximation) given in Equation (22) of Althaus et al. (2003). If one wishes to model a helium atmosphere white dwarf, then one simply sets the location of the base of the hydrogen layer, −log M_H (a free parameter set in an input file), to 20 or any reasonable value that puts the hydrogen beyond the surface of the model.

The chemical profiles span the core and the envelope. While the modern equation of state tables allow for any composition one may fancy, the code is set up to accept chemical profiles that are consistent with what we know from stellar evolution calculations (e.g., Althaus et al. 2005, 2010). As a result, there is an inner region where a mix of oxygen, carbon, and helium is expected and an outer region where a mix of helium and hydrogen is expected (for a helium atmosphere white dwarf, that region can be pure helium). It is currently not possible to have all four elements present at any point in the model. This is illustrated in Figure 5.

Zoom In Zoom Out Reset image size

3.3.2. Helium Cores

In essence, the helium core version of the code is the same as above, except the abundance of carbon and oxygen are set to zero in the core, and there is a single transition zone from helium to hydrogen. That transition can happen in stages (partial diffusion, see Figure 6). The definition of the composition profiles in this version of the code is entirely parameterized. The user is required to provide the locations of the transition from pure helium in the core to a homogeneous mix of helium and hydrogen, and the transition from mixed helium/hydrogen to pure hydrogen. Other parameters define the abundance of helium in the mixed helium/hydrogen region and the shape of the transitions. This is illustrated in Figure 6.

Zoom In Zoom Out Reset image size

If one so chooses, the WDEC in its current form will calculate oscillation modes for a model. Periods of oscillation are calculated using the adiabatic code "cjhanro" by Carl Hansen, best described in (Kawaler et al. 1985) and references therein. Aside from being rewritten using Fortran90 syntax, no substantive changes were made to the code.

4.1. New Method to Compute the B − V Frequency

5.1. Convection

We checked the treatment of convection in WDEC by using a separate, independent envelope code (the "Warsaw" envelope code; Paczyński 1969, 1970; Pamyatnykh 1999) and compared the depth of the convection zone we were getting for our hydrogen and helium atmospheres. To run this test, we operated the WDEC in the mode where the values of α are determined using the nonlinear light curve fitting calibration described in Section 3.2 (see also Figure 3). Most importantly, this is a mode where the value of α varies according to the surface gravity and effective temperature of the model. We chose this mode because it is likely the one adopted by most users most of the time.

The results of such a test are shown in Figure 7. We checked a wide range of masses (0.35–1.0 M_⊙) and show the results for the 0.60 and 1.0 M_⊙ models, which are representative of the type of agreement we obtained. The vertical axis is set up such that the mass coordinate on the vertical axis is high when the convection zone is thin and low when the convection zone is thick. The surfaces of the models are near $-\mathrm{log}(1-{M}_{\mathrm{conv}}/{M}_{* })=20$.

Zoom In Zoom Out Reset image size

For the helium atmosphere white dwarfs, the agreement between the two codes is good over the entire range of temperatures over which convection takes place. For hydrogen atmosphere white dwarfs, the agreement is not as good for temperatures higher than ∼11000 K. We note that for those thinner convection zones, the thickness of the convection zone is more sensitive to α, and a small change in this parameter can lead to a much thinner convection zone (or one that numerically goes down to zero).

5.2. Pulsations

Another way to check our models is to compare them to models independently produced. We used a set of models from a group at the Universidad Nacional de La Plata (Romero et al. 2012). Their codes (LPCODE and LP-PUL) are described in Althaus et al. (2005) and Córsico & Althaus (2006) and references therein. We reproduced one of their models with WDEC and compare essential properties in this section. Periods of oscillation are a sensitive probe of the interior structure of models. Asteroseismology is based on the idea that observed periods can be used to infer the interior structure of stars. If two independent codes give correct solutions of the equations of stellar structure and of non-radial oscillations, models with identical stellar parameters and chemical profiles should result in identical sets of periods. Comparing pulsational properties is a sensitive and comprehensive test.

We compared the periods obtained for the LPCODE model shown in Figure 8 to those calculated for the WDEC model. The results are shown in Figure 9. For higher radial overtones (k > 10), periods differ by less than 2%, while at lower radial overtones, they differ by as much as nearly 6%. The agreement at larger radial overtones is a sign that the models are structurally very similar (and also a sign that the period computations are yielding consistent results). The discrepancies at lower radial overtones are a testimony to how sensitive the modes are to chemical transitions. The seemingly minute differences seen in Figure 8 translate into significant differences in the pulsation periods. Asteroseismology rests on the sensitivity of low k modes in particular to help us infer the interior structure of stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

WDEC offers a fast and fairly easy way to produce models of white dwarfs. It can also be used to obtain g-mode oscillation modes for the models. In this paper, we detailed a recent overhaul made to this code, most significantly the inclusion of MESA modules for equations of states and opacities, and improvements on the treatment of convection, using simulations and data that have become available.

The code is open source and may be obtained from GitHub (Bischoff-Kim & Montgomery 2018, Codebase: https://github.com/kim554/wdec). Further documentation, including a user manual, may be found in the GitHub repository.

We wish to thank Don Winget, Steve Kawaler, and Paul Bradley for helping us find sources for the discussion on the background of the code. Special thanks to Gilles Fontaine for placing the code in its historical context and for providing pointed feedback as referee. Winget and Fontaine were also a resource in helping us understand some of the detailed aspects of the code.

Software: "Warsaw" envelope code; (Paczyński 1969, 1970; Pamyatnykh 1999), LPCODE (Althaus et al. 2005), MESA (v8118; Paxton et al. 2011, 2013, 2015), LP-PUL (Córsico & Althaus et al. 2006).