Modules for Experiments in Stellar Astrophysics (MESA): Pulsating Variable Stars, Rotation, Convective Boundaries, and Energy Conservation

Since the energy density of radial pulsations drops rapidly going inward from a star's surface, a full stellar model reaching to the center is frequently not necessary. The use of RSP is currently restricted to cases in which pulsations are determined by the structure of the envelope and are independent of the detailed structure of the core. RSP begins by building a chemically homogeneous envelope from given stellar parameters (M, L, T_eff, X, and Z). These parameters can be freely chosen and need not originate from a MESAstar model. It is not yet possible to directly import an envelope from MESAstar into RSP primarily because of the different treatments of convection (a version of MLT in MESAstar vs. detailed time evolution of turbulence in RSP). Tighter integration of MESAstar and RSP is a future project.

Specifications for the initial model include the number of cells and the temperature at the base (see Figure 2). This inner boundary temperature is defined by a chosen temperature (RSP_T_inner ≃ 2 × 10⁶ K) that should be set hot enough so that the eigenvector amplitudes generated in the following stability analysis go to zero and cool enough to exclude regions of nuclear burning and justify the assumption of chemical homogeneity. The model is divided into inner and outer regions at a specified anchor temperature. In the outer region, cells have the same mass; in the inner region cell masses grow by a constant factor so that the innermost cells are significantly larger than the ones at the surface. The anchor temperature should be in the part of the model driving the pulsations. For example, for pulsations in the classical instability strip a value of T_anchor = 11,000 K is typical. In the case of Z-bump pulsations a higher temperature would be appropriate. Proper choice of the number of outer cells and placement of the anchor are necessary to ensure that the driving region is well resolved.

The initial model builder iteratively constructs an envelope in hydrostatic equilibrium that satisfies the RSP equations. Starting from the outer radius determined by L and T_eff, this process involves selection of a cell mass to be used in the outer part of the envelope and a scale factor that is used to progressively increase cell masses in the inner region. Those choices must match the desired number of cells, both N and N_outer, and also satisfy the surface boundary conditions and the required temperatures at the anchor location and at the inner boundary. The model builder is a complex multistage iterative procedure that works well for the range of cases presented in the following but may fail when applied outside of that range.

The LNA analysis is performed on the initial model using a full linearization of the RSP equations (for details see Smolec 2009). These include time-dependent convection, moving beyond the frozen-in convection approximation made in software instruments like GYRE. This yields the eigenmodes, periods, and growth rates. The eigenvectors are used to perturb the initial model for the time evolution.

Figure 3 shows amplitudes of radial displacements and differential work for the first three eigenmodes of the classical Cepheid model in the MESA test_suite. A common resolution for exploratory model surveys, N = 150 and N_outer = 40, is adopted. For T ≳ 5 × 10⁵ K, the displacements and differential work of all three radial eigenmodes are negligible, indicating that the extent of the computational domain is sufficient.

Zoom In Zoom Out Reset image size

The initial static model is perturbed with a linear combination of the velocity eigenvectors of the three lowest-order radial modes. More specifically, the velocity eigenvectors are scaled to have a surface value of 1. RSP_fraction_1st_overtone and RSP_fraction_2nd_overtone multiply the 1O and 2O eigenvectors, respectively. The F-mode eigenvector is then multiplied by (1-RSP_fraction_1st_overtone - RSP_fraction_2nd_overtone). The linear combination of these three scaled eigenvectors is then multiplied by the surface velocity RSP_kick_vsurf_km_per_sec.

The time integration commences with a constant time step (RSP_target_steps_per_cycle) and continues for a specified number of pulsation cycles (RSP_max_num_periods). A new cycle begins when the model passes through a maximum radius. Controls allow filtering out secondary maxima in the radial velocity curve.

Figure 4 shows Γ, the fractional growth of the kinetic energy per pulsation period near the start of a time integration, where

$$\begin{eqnarray}&&{\rm{\Gamma }}=2({E}_{{\rm{k}},\max }^{i+1}-{E}_{{\rm{k}},\max }^{i})/({E}_{{\rm{k}},\max }^{i+1}+{E}_{{\rm{k}},\max }^{i}),\end{eqnarray} \tag{ 5 }$$

and ${E}_{{\rm{k}},\max }^{i}$ is the maximum kinetic energy of the envelope during pulsation cycle i. Agreement between these three time integrations and the corresponding LNA analyses is satisfactory. Similarly, the pulsation periods match the linear values during the low-amplitude phase of development. Consistency between the time integrations and LNA analyses forms the basis for interpreting the nonlinear results.

Zoom In Zoom Out Reset image size

Which of the perturbed modes attains large-amplitude pulsations in the nonlinear regime may depend on the initial conditions (e.g., Smolec 2014). Figure 5 shows the results of longer time integrations for the RR Lyrae model shown in Figure 4. The upper triplet of panels, case (a), is for a 1O-mode initialization with a 4.5 km s⁻¹ amplitude. The middle triplet panel, case (b), is for an F-mode initialization with 4.5 km s⁻¹ amplitude. The lower triplet, case (c), is for an F-mode initialization with a 9.5 km s⁻¹ amplitude.

Zoom In Zoom Out Reset image size

For case (a), the pulsations converge toward a single, 1O-mode pulsation. After a ≃500-cycle transient phase, the pulsation period and radius amplitude barely change and Γ ≃ 0. For case (b), the model has not converged to a single-periodic mode after 4000 cycles. Despite the pure F-mode initialization, at ≃300 cycles the pulsation switches toward the 1O-mode. This does not prove that the model cannot pulsate in the F-mode, as case (c) demonstrates. After a transient phase with beating F and 1O-modes, the 1O-mode decays and the single-periodic F-mode pulsation grows to saturation.

Figure 5 is an example of two different single-mode solutions whose selection depends on the initial conditions. Two stars can have the same physical parameters but pulsate in different modes depending on their evolutionary history.

The final state in the nonlinear regime is usually a single-periodic oscillation. The shape of the light and radial velocity curves may depend on the values of the eight free parameters listed in Table 3. In Table 4 set A corresponds to the simplest convection model. Set B adds radiative cooling, set C adds turbulent pressure and turbulent flux, and set D includes these effects simultaneously. The parameter α_m has little effect on the shape of the light curve but strongly affects its amplitude. Figure 6 shows the effect of varying α_m on the T_eff = 6000 K Cepheid model. The free parameter α_m may thus be used to match the observed amplitude. For sets A–D, α_m was adjusted so that models with different sets have similar amplitudes.

Zoom In Zoom Out Reset image size

Figure 7 shows the effect of these parameter sets on the shapes of the bolometric light, photosphere radial velocity, and radius variation curves for a saturated F-mode RR Lyrae and two saturated F-mode classical Cepheid models. The pulsation periods, radial velocity curves, and radius variation curves show only small differences. For the RR Lyrae models, there are differences in the fine structure of the light curves. For example, the bump before minimum light is weaker when turbulent pressure and turbulent flux are included (sets C and D). The shape of the light curve near maximum light also differs for both the RR Lyrae and Cepheid models.

Zoom In Zoom Out Reset image size

Figure 8 shows the convective luminosity profiles for the models of Figure 7. Depending on pulsation phase, one convective region (darker hues) extends from the surface cells down to cell ≃90. This convective region is associated with partial ionization of H and He. Another convective region (lighter hues) lies deeper in the envelope, centered at cell ≃110, and is associated with the second ionization of He. In most of the models these two convective regions merge at pulsation phase ≃0.5 during maximum contraction when both convective regions are at their strongest and most extended. In the cooler models, the first convective region carries nearly all of the luminosity throughout a pulsation cycle. In the hotter models, this convective region becomes very weak at pulsation phase ≃0.8 (before maximum expansion) and is barely resolved as it progresses deeper into the envelope. This behavior is pronounced in the RR Lyrae models with radiative cooling (sets B and D), as cooling contributes to damping the turbulent energy and hence the near disappearance of the convective region.

Zoom In Zoom Out Reset image size

There are two parameters, α_cut and C_q, that control the Stellingwerf (1975) artificial viscosity, entering into the definition of p_av (see Table 1). Figure 9 shows the effect of these parameters on light curves for the T_eff = 6000 K Cepheid model. The value of C_q plays a minor role. For C_q = 4, the top panel shows that a different light curve develops only if α_cut = 0.0, corresponding to artificial viscosity acting for very small compressions, which leads to excessive dissipation that quenches the pulsation amplitude. For α_cut ≥ 0.01, the light curves are similar and roughly have the same pulsation amplitude. When α_cut = 0.1, artificial viscosity turns on only for strong shocks, seen as the wiggle on the ascending branch of the light curve. This choice (α_cut = 0.1) numerically captures shocks without excessive dissipation, barely affecting the light-curve shape and amplitude. While an artificial viscosity modifies the velocity structure in the envelope at each epoch, we find that these differences are smaller than the differences for the bolometric light curves. The bottom panel shows that the Tscharnuter & Winkler (1979) form of artificial viscosity yields light curves with the same amplitude and qualitatively the same shape. Small differences are apparent at a shock-prone phase shortly before the maximum brightness.

Zoom In Zoom Out Reset image size

Figure 10 shows the sensitivity of the bolometric light curve to the total number of cells N, number of cells above the anchor N_outer, anchor location T_anchor, and inner boundary location T_inner for an RR Lyrae model (M = 0.65 M_⊙, L = 50 L_⊙, T_eff = 7000 K, X = 0.75, Z = 0.0014) with convective parameter set A. For classical pulsators, T_inner is typically placed at 2 × 10⁶ K, and common choices for T_anchor are 11,000 K or 15,000 K. The top panel shows that the light curves are weakly sensitive to the choice of T_inner and T_anchor for this RR Lyrae model. Light and radial velocity curves are usually the most sensitive to N and N_outer. The bottom panel shows that this effect is small for this RR Lyrae model. Section 2.4.4 shows a case with a much larger sensitivity.

Zoom In Zoom Out Reset image size

The default value of 600 time steps per pulsation cycle works well for most cases, but smaller time steps are recommended for models that include radiative cooling, turbulent pressure, or turbulent flux, or develop violent pulsations (e.g., the chaotic models of Section 2.4.3). We stress that there is no unique choice of grid or time step that will work for all applications or guarantees convergence. All nonlinear modeling of variable stars should be accompanied by sensitivity and convergence tests.

We consider two sequences of RR Lyrae-type models. The first sequence has M = 0.65 M_⊙, L = 45 L_⊙, and [Fe/H] = −1.0 (X = 0.75, Z = 0.0014), with T_eff varying in 100 K steps for convective sets A–D. Figure 11 shows a gallery of I-band light curves from this sequence. The top panel shows F-mode pulsators (commonly known as RRab stars), and the bottom panel shows 1O-mode pulsators (known as RRc stars). The latter have smaller amplitudes and are less nonlinear in shape (i.e., more sinusoidal) than the F-mode light curves. Models with different convective settings differ the most near minimum and maximum light. For example, F-mode models with convective set B develop a bump preceding minimum light that is absent in the light curves with convective set D. On the other hand, F-mode models with cooler T_eff from convective set D develop broad, double-peaked light-curve maxima that are absent in models from convective set B.

Zoom In Zoom Out Reset image size

To compare the overall morphology of I-band light curves from these sequences with OGLE observations, we perform a Fourier decomposition of the synthetic light curves

$$\begin{eqnarray}&&I(t)={A}_{0}+\displaystyle \sum _{k}{A}_{k}\sin (2\pi {kft}+{\phi }_{k}),\end{eqnarray} \tag{ 6 }$$

where f is the pulsation frequency and A_k and ϕ_k are amplitudes and phases, respectively. We then construct the amplitude ratios R_k1 and epoch-independent phase differences φ_k1 (Simon & Lee 1981):

$$\begin{eqnarray}&&{R}_{k1}=\displaystyle \frac{{A}_{k}}{{A}_{1}},\quad {\varphi }_{k1}={\phi }_{k}-k{\phi }_{1}.\end{eqnarray} \tag{ 7 }$$

Observationally derived values of R_k1 and φ_k1 are taken from the OGLE catalog (Soszyński et al. 2014) and shown in Figure 12 by gray circles for RRab (F-mode) stars and blue circles for RRc (1O-mode) stars. The observations show that the Fourier parameters follow progressions with pulsation period, traced by the highest density of data points, but with significant scatter. Fourier parameters from the model I-band light curves are shown with colored symbols. Left panels are for the first sequence of models. Right panels are for the second sequence, which has M = 0.65 M_⊙, convective set B, T_eff varying in 100 K steps, and either L = (40, 45, 50) L_⊙ at [Fe/H] = −1.0 or L = 45 L_⊙ at [Fe/H] = −1.5.

Zoom In Zoom Out Reset image size

The left panels show that F-mode pulsators with convective sets A and B progress similarly. Models with convective sets C and D also progress similarly but are qualitatively different from those with convective sets A and B. These differences are more pronounced for cooler, longer-period models. The overall match of the F-mode decompositions with the OGLE RRab stars (gray circles) is reasonable but shows some systematic differences. For example, the model values of φ₃₁ are larger than the observed values. However, the model physical parameters except T_eff are fixed, while this is not the case for the OGLE RRab stars.

The match between the 1O-mode models and the OGLE RRc stars (blue circles) is worse. The amplitudes and the amplitude ratios are systematically too large. The Fourier phases are systematically too small. This may indicate that different convective parameters are needed to reproduce the observed light-curve shapes of F-mode and 1O-mode pulsators. Note that the 1O-mode instability strip with convective set C is narrower than the F-mode instability strip at the luminosity considered. Models from set C, where the 1O-mode is linearly unstable and the integration is initialized with a 1O-mode velocity perturbation, all switch to an F-mode pulsation.

The right panels in Figure 12 show that the light-curve shapes are sensitive to the physical parameters. By varying the luminosity and metallicity in a narrow range, the model sequences match the OGLE values. However, no sequence considered follows the OGLE progression, because RR Lyrae in the Galactic bulge are characterized by mass, luminosity, and metallicity distributions that cannot be reproduced with a single sequence.

Cepheids display a feature known as the Hertzsprung progression (Hertzsprung 1926). A secondary bump in their light and radial velocity curves appears near minimum light on the descending branch when P ≃ 5 days. The bump moves toward earlier phases on the descending branch as the period increases and is coincident with maximum light when P ≃ 10 days. The bump then moves onto the ascending branch for longer periods and disappears at P ≃ 20 days. The bump is driven by a 2:1 resonance between the F-mode and a damped 2O-mode (e.g., Simon & Schmidt 1976; Buchler et al. 1990). This behavior is reflected in Cepheids' Fourier parameters, which follow more complex progressions with pulsation period than the RR Lyrae stars.

We consider 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, and 10.0 M_⊙ models with X = 0.736, Z = 0.008, and convective parameter sets A–D. Further details are in the test suite example 5M_cepheid_blue_loop. Figure 13 shows the evolutionary tracks during the core He-burning, blue-loop phase for a selection of masses. Blue and red edges of the instability strip computed with RSP for convective sets B and D are shown. Edges for convective sets A and C largely overlap those for convective sets B and D, respectively. Nonlinear pulsation models are computed with a ΔT_eff = 300 K or 500 K offset from the blue edge.

Zoom In Zoom Out Reset image size

Figure 14 shows I-band light curves and radial velocity curves for the ΔT_eff = 300 K models. Due to the shift in the location of blue edge models, models from convective sets B and D have different T_eff and consequently different pulsation periods. Figure 15 compares Fourier parameters of the Cepheid models with those derived from LMC Cepheid light curves (left panel; Ulaczyk et al. 2013; Soszyński et al. 2015) and Galactic F-mode Cepheid radial velocity curves (right panel; Storm et al. 2011). The I-band light curves are more sensitive to the convective parameters than the radial velocity curves for both ΔT_eff offsets. As with the RR Lyrae models, the Fourier parameters for convective sets A and B are similar to each other, as are those for sets C and D. The radial velocity Fourier parameters follow tighter progressions than the I-band light curves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Fourier phases in the left panels of Figure 15 are systematically larger than the observationally inferred values for P ≲ 10 days, with the difference being larger for the cooler ΔT_eff = 500 K models. For P ≳ 10 days the model Fourier phases are systematically smaller. Large discrepancies for the radial velocity curves are absent in the right panels of Figure 15, except for the amplitudes and R₂₁ ratio at the shortest periods. The projection (or p) factor is the ratio of the pulsation velocity to the radial velocity deduced from spectral line-profile observations, dependent on at least rotation and gravity darkening (see Section 4), and plays a role in the amplitudes. We use p = 1.3, close to the average value of determinations based on eclipsing binary Cepheid systems (Pilecki et al. 2018) and interferometric methods (e.g., Breitfelder et al. 2016).

This brief survey is not exhaustive, as only a few masses and two model sequences are explored. The M−L relation is also important for Cepheids, as evolutionary tracks depend on overshooting, rotation, and metallicity.

Type II Cepheids are more similar to RR Lyrae stars than to classical Cepheids owing to their lower masses (M ≃ 0.5 M_⊙), Population II chemical composition, and evolutionary history. Their masses are similar to RR Lyrae, but they cross the instability strip at larger luminosities and pulsate with longer periods (P > 1 day). Type II Cepheids are F-mode pulsators except for a few double-mode stars pulsating simultaneously in the F- and 1O-modes (Smolec et al. 2018; Udalski et al. 2018) and two recently discovered 1O-mode pulsators (Soszyński et al. 2019). BL Her variables are a subclass of Type II Cepheids with P ≲ 4 days.

Nonlinear radiative models of Type II Cepheids revealed a variety of complex dynamics, including period-doubled and deterministic chaos pulsations (e.g., Buchler & Kovacs 1987; Kovacs & Buchler 1988; Buchler & Moskalik 1992). With convective pulsation models, modulated pulsations were also found (Smolec & Moskalik 2012). Buchler & Moskalik (1992) discovered period-doubled pulsations in their survey of BL Her models and predicted that they should be observed in BL Her variables. The period doubling is caused by a 3:2 resonance of the F- and 1O-modes and nonlinear phase synchronization (Moskalik & Buchler 1990); pulsations repeat only after two cycles of the F-mode.

Soszyński et al. (2011) report a period-doubled BL Her star, OGLE-BLG-T2CEP-279. We adopt M = 0.6 M_⊙, L = 184 L_⊙, T_eff = 6050 K, X = 0.76, Z = 0.01, and convective parameter set A. These are nearly the same physical parameters Smolec et al. (2012) chose for a T2CEP-279 model survey. The observed and model light curves are shown in Figure 16. The model amplitudes of the bump at minimum light of the period-doubled cycle are larger and the shape of the light maximum is more pronounced relative to the observed features. The model period, P_model = 2.6976 days, is longer than the observed period, P_obs = 2.3993 days. Still, the light curves are qualitatively similar, devoid of fine-tuning, and demonstrate that RSP can be used to model specific stars.

Zoom In Zoom Out Reset image size

Figure 17 shows the radius variation for two models with M = 0.55 M_⊙, L = 136 L_⊙, X = 0.76, Z = 0.0001, and convective set A but with a reduced eddy viscosity, α_m = 0.05 (yielding unrealistic light curves). The two models differ only in T_eff = 5932 K (top panel) and T_eff = 6410 K (bottom panel). Figure 18 shows the return map of maximum radii for these two models, plotting the maximum radii for each pulsation cycle ${R}_{\max }^{n}$ versus the preceding ${R}_{\max }^{n-1}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the cooler T_eff = 5932 K model, the top panel of Figure 17 shows a cyclic modulation in the envelope of the period-doubled pulsation. The modulation period of ≃57 days is longer than the pulsation period of ≃2.3 days. The return map in the top panel of Figure 18 is constructed from ≃8000 pulsation cycles and shows two loops, corresponding to alternating smaller and larger maximum radii. Since the modulation period is not commensurate with the pulsation period, the return maps develop a locus of points that form the closed lobes. Light-curve modulation is common in RR Lyrae stars (Blazhko effect), and periodic pulsation modulation was recently discovered in BL Her variables (Smolec et al. 2018).

For the hotter T_eff = 6410 K model, the radius variation appears irregular in the bottom panel of Figure 17. The return map in the bottom panel of Figure 18 reveals a strange attractor, an example of deterministic chaos in nonlinear models. Tracing time series from models is simple, but tracing chaotic dynamics in observations is difficult (e.g., Plachy et al. 2018). Chaotic dynamics is reported in a few type II Cepheids with longer periods, in the RV Tau variable star range, and in semiregular variable stars (e.g., Buchler et al. 1996, 2004; Kollath et al. 1998; Plachy et al. 2018). While the T_eff = 6410 K model has a shorter period, in the BL Her range, such models may provide insight into chaotic dynamics in pulsating stars (see Plachy et al. 2013; Smolec & Moskalik 2014).

Very low mass stars do not enter the classical instability strip within a Hubble time. However, mass loss from the more massive component in a close interacting binary can lead to a low-mass star that then evolves through the instability strip. The Binary Evolution Pulsator (BEP, OGLE-BLG-RRLYR-02792) is the prototype of this new class of pulsators (Pietrzynski et al. 2012; Smolec et al. 2013). The BEP's variability is similar to an F-mode RR Lyrae pulsator but with a dynamical mass of ≃0.26 M_⊙.

Following Smolec et al. (2013), we adopt M = 0.26 M_⊙, L = 33 L_⊙, T_eff = 6910 K, X = 0.7, Z = 0.01, and convective set A. Figure 19 shows I-band light curves and radial velocity curves for the BEP models. Results for a coarse grid (N/N_outer = 150/40, gold curves) qualitatively match these observations and have a period of 0.6373 days close to the observed 0.6275-day period. We also show results for a medium grid (300/120, orange curves) and a fine grid (300/120, red curves). The differences are most pronounced around maximum and minimum light. The shape of the light and radial velocity curves approach convergence only on grids with ≳600 cells. The amplitude of the model curves is sensitive to convective parameters and can be fine-tuned for a better match with observations.

Zoom In Zoom Out Reset image size

The origin of BLAPs, introduced in Section 1, is unknown. They have been modeled as ≃0.3 M_⊙ shell H-burning stars that are progenitors of low-mass WDs and ≃1.0 M_⊙ stars undergoing core He burning (Pietrukowicz et al. 2017; Byrne & Jeffery 2018; Romero et al. 2018; Wu & Li 2018). Though mass loss in a close interacting binary must be invoked for both hypotheses, none of the BLAPs are known to be in a binary. Figure 1 shows that BLAPs are located near sdBVs in the H-R diagram. The latter have noncanonical abundance profiles that are strongly affected by radiative levitation (see Section 6.2). Following the linear study of Romero et al. (2018), we adopt Z = 0.05, to account for the increased envelope metallicity caused by radiative levitation.

We explored nonlinear models with M = 1.0 M_⊙, T_eff = 30,000 K, L = 430 L_⊙, and three convection parameter sets. These envelope models used N = 280, N_outer = 140, T_anchor = 2 × 10⁵ K, and T_inner = 6 × 10⁶ K. The LNA analyses show that the F- and 1O-modes are unstable. Figure 20 compares the OGLE-BLAP-011 and 1O-mode (P ≃ 35 minutes) model I-band light curves. The qualitative agreement is reasonable. Figure 20 shows that the model radial velocity curves have amplitudes of ≃200 km s⁻¹, a pronounced temporal asymmetry, and light maxima that are narrower than light minima.

Zoom In Zoom Out Reset image size

These BLAP models lie far off the classical instability strip so that the initial model builder (see Section 2.2.1), which is optimized for classical pulsators, failed to relax the initial models to complete hydrostatic equilibrium. An option is to switch off the relaxation process and commence the time integration with a near-hydrostatic-equilibrium initial model. The price to pay is that the LNA growth rates are only indicative, not accurate.

Figure 21 shows the evolution of Γ for the (α, α_m) = (1.5, 0.1) model in Figure 20. Before ≃10 days, Γ fluctuates around a mean of ≃7.5 × 10⁻³. The fluctuations diminish with time, but Γ remains above the LNA value of 6.25 × 10⁻³ up to ≃30 days. Between 30 and 50 days Γ diminishes and approaches zero, the sign that the nonlinear pulsation saturates at its terminal amplitude, yielding the results of Figure 20. In contrast to Γ, the period of the nonlinear pulsations remains close to the LNA period. In cases where the initial model cannot be relaxed, the initial Γ should not be expected to match the LNA analysis.

Zoom In Zoom Out Reset image size

High-amplitude δ Scuti (HADS) pulsators are defined to have V-band light-curve amplitudes greater than 0.1 mag. HADS pulsators lie close to the MS (see Figure 1), where growth rates are usually much smaller than those for RR Lyrae stars or classical Cepheids. This implies that long time integrations are needed to drive nonlinear pulsations to saturation.

We consider a stellar model with M = 2 M_⊙, L = 30 L_⊙, T_eff = 6900 K, X = 0.7, Z = 0.01, and convective set A. This represents a star evolving toward the Hertzsprung gap. The LNA analysis of the initial model reveals that the F- and 1O-modes are linearly unstable, with growth rates of 1 × 10⁻⁶ and 6 × 10⁻⁵, respectively.

Section 3 emphasizes the importance of numerical energy conservation. Figure 22 shows the evolution of the relative cumulative error in the energy for a 250,000-cycle integration (≃150 million time steps, at 600 steps per cycle). The relative cumulative error grows from ≃3 × 10⁻¹¹ after about 15 cycles to ≃3 × 10⁻⁷ after 250,000 cycles. The inset figure shows that the per-step relative error in the energy scatters around −2.5 × 10⁻¹⁵ but is systematically different than zero.

Zoom In Zoom Out Reset image size

Figure 23 shows that after 70,000 cycles the asymmetric V-band light curves have an amplitude of ≃0.2 mag. The dominant pulsation period is 0.127 days, close to the LNA period of 0.1269 days for the 1O-mode. The amplitude, though, varies cyclically over about four pulsation cycles, reflecting the presence of the F-mode. Continuing the integration to 130,000 cycles leads to a more saturated light curve with an amplitude of ≃0.3 mag and an unchanged dominant pulsation period. Extending the integration to 250,000 cycles does not lead to any significant changes.

Zoom In Zoom Out Reset image size

The robustness in the light curves between 130,000 and 250,000 cycles might suggest that the long-term behavior is that of a double-mode HADS model. However, additional integrations with different initial perturbations, chosen to more adequately sample the neighboring amplitude phase space, are necessary to assess the possible mode selection. Figure 24 shows the evolution of these integrations in the amplitude–amplitude diagram using the analytical signal method (e.g., Kolláth et al. 2002). The amplitude behavior of the model integration shown in Figure 23 is traced out by the red line. After initially rapid evolution, all the amplitude trajectories develop an arc along which evolution slows markedly. The model in Figure 23 and neighboring trajectories bend to the left toward smaller F-mode amplitudes. Their likely final state is that of a single-periodic 1O pulsation. In contrast, the two rightmost trajectories bend to the right toward larger, more dominant F-mode amplitudes. Despite the limited sampling of phase space in Figure 24, we cautiously conclude that the most likely long-term outcome of this δ Sct model is a single-periodic 1O-mode or F-mode pulsator depending on the initial conditions (see Section 2.2.4).

Zoom In Zoom Out Reset image size

Because MESA works on a Lagrangian mesh, it handles accretion and mass loss in a two-stage process (Paper III, Section 7). In stage I the masses of certain cells are increased or decreased as needed to give the desired end-of-step total mass, but no attempt is made to ensure energy conservation at this stage. In stage II the model thus produced is evolved in time by an amount $\delta t$. This separation of stages means that the time step is only taken for a model of fixed mass. However, because the energy of the model changes in stage I, a correction must be added to stage II to make the overall step consistent with energy conservation. Thus, we introduce a new source term ${\epsilon }_{\dot{M}}$ that accounts for the heating associated with mass changes.

The change in the mass of cell k during stage I results from the difference between the outward²² mass flux F_m through each cell face. This flux obeys

$$\begin{eqnarray}&&{{dm}}_{\mathrm{mid},k}-{{dm}}_{\mathrm{start},k}=\delta t\left({F}_{m,k+1}-{F}_{m,k}\right)\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&{M}_{\mathrm{mid}}-{M}_{\mathrm{start}}=-\delta {{tF}}_{m,1}=\delta t\dot{M}.\end{eqnarray} \tag{ 14 }$$

During stage I the temperature, density, and velocity of each cell are held fixed, but the composition is updated to track the flow of material between cells. The subscript start is used for quantities at the start of stage I. The subscript mid is used for quantities at end of stage I, which is the start of stage II. No subscript is used for quantities evaluated at the end of the time step (after stage II). There is no mass change during stage II, so dm_k = dm_mid,k. In the following we write dm_k rather than dm_mid,k.

The change in energy for cell k during stage I is

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{mid},k}-{E}_{\mathrm{start},k} & = & {{dm}}_{k}{{ \mathcal E }}_{\mathrm{mid},k}-{{dm}}_{\mathrm{start},k}{{ \mathcal E }}_{\mathrm{start},k}\\ & = & {{dm}}_{k}\left({{ \mathcal E }}_{\mathrm{mid},k}-{{ \mathcal E }}_{\mathrm{start},k}\right)\\ & & +({{dm}}_{k}-{{dm}}_{\mathrm{start},k}){{ \mathcal E }}_{\mathrm{start},k},\end{array}\end{eqnarray} \tag{ 15 }$$

where ${{ \mathcal E }}_{k}$ is the total specific energy of cell k, given by the sum of specific potential, kinetic, and internal energies. Neglecting changes in specific energy owing to changes in composition, the difference in ${{ \mathcal E }}_{k}$ across stage I is

$$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{mid},k}-{{ \mathcal E }}_{\mathrm{start},k}={\left(\displaystyle \frac{{{Gm}}_{k,{\rm{C}}}}{{r}_{k,{\rm{C}}}}\right)}_{\mathrm{start}}-{\left(\displaystyle \frac{{{Gm}}_{k,{\rm{C}}}}{{r}_{k,{\rm{C}}}}\right)}_{\mathrm{mid}}.\end{eqnarray} \tag{ 16 }$$

where m_k,C and r_k,C are the mass coordinate and radius, respectively, at the center of mass of the cell (Paper IV).

We now introduce _k,I, the effective source term in cell k during stage I. This is defined by writing the change in energy in flux-conservative form as

$$\begin{eqnarray}&&{E}_{\mathrm{mid},k}-{E}_{\mathrm{start},k}=\delta t\left({{dm}}_{k}{\epsilon }_{k,{\rm{I}}}+{F}_{e,k+1,{\rm{I}}}-{F}_{e,k,{\rm{I}}}\right),\end{eqnarray} \tag{ 17 }$$

where F_e,k,I is the outward flux of energy across face k owing to work and material passing through face k. Inserting Equation (15) and rearranging we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\epsilon }_{k,{\rm{I}}} & \equiv & \displaystyle \frac{{{ \mathcal E }}_{\mathrm{mid},k}-{{ \mathcal E }}_{\mathrm{start},k}}{\delta t}+\displaystyle \frac{{{ \mathcal E }}_{\mathrm{start},k}}{\delta t}\left(1-\displaystyle \frac{{{dm}}_{\mathrm{start},k}}{{{dm}}_{k}}\right)\\ & & -\displaystyle \frac{{F}_{e,k+1,{\rm{I}}}-{F}_{e,k,{\rm{I}}}}{{{dm}}_{k}}.\end{array}\end{eqnarray} \tag{ 18 }$$

The energy flux is

$$\begin{eqnarray}&&{F}_{e,k,{\rm{I}}}={p}_{\mathrm{face},k}{{ \mathcal A }}_{k}{\hat{v}}_{k,{\rm{I}}}+{{ \mathcal E }}_{\mathrm{face},k}{F}_{m,k},\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}&&{\hat{v}}_{k,{\rm{I}}}\equiv \displaystyle \frac{{r}_{\mathrm{mid},k}-{r}_{\mathrm{start},k}}{\delta t}\,,\end{eqnarray} \tag{ 20 }$$

and the face values ${{ \mathcal E }}_{\mathrm{face}}$ are interpolated²³ from the cell values ${{ \mathcal E }}_{\mathrm{start}}$. The radial coordinate after state I is calculated using the updated cell masses, holding cell densities fixed. The change in total energy of the model during stage I is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{k}\left({E}_{\mathrm{mid},k}-{E}_{\mathrm{start},k}\right) & = & -\delta t\left({p}_{\mathrm{face},1}{{ \mathcal A }}_{1}{\hat{v}}_{1,{\rm{I}}}+{{ \mathcal E }}_{\mathrm{face},1}{F}_{m,1}\right)\\ & & +\delta t\displaystyle \sum _{k}{{dm}}_{k}{\epsilon }_{k,{\rm{I}}}.\end{array}\end{eqnarray} \tag{ 21 }$$

The term in parentheses on the right-hand side accounts for the energy of new material entering or leaving the model and the work done in the process. The additional _k,I term implies the need for a corrective source term ${\epsilon }_{\dot{M}}$ that must be added during stage II so that energy is properly accounted for.

To determine this new source, we consider the ratio of the thermal timescale,

$$\begin{eqnarray}&&{\tau }_{\mathrm{th},k}\equiv {\left|\displaystyle \frac{{x}_{m}{c}_{{\rm{p}}}T({{\rm{\nabla }}}_{\mathrm{rad}}-{{\rm{\nabla }}}_{\mathrm{ad}})}{L}\right|}_{\mathrm{start},k},\end{eqnarray} \tag{ 22 }$$

to the mass-change timescale

$$\begin{eqnarray}&&{\tau }_{{\rm{m}},k}\equiv {\left|\displaystyle \frac{{x}_{m}}{{F}_{m}}\right|}_{\mathrm{start},k},\end{eqnarray} \tag{ 23 }$$

where x_m,k is the mass above face k. When this ratio is small, the thermodynamic state of material is a function primarily of depth (Sugimoto 1970; Sugimoto & Nomoto 1975). This means that as material moves from cell to cell it is a good approximation to the time evolution to suppose that it adopts the state of whichever cell it is in (Paper III). In the opposing limit (τ_th ≫ τ_m) the entropy of material adjusts minimally as it moves from cell to cell.

We account for the effects of the thermal and mass-change timescales by tracking the heating of material as it moves from cell to cell in stage I and estimating what fraction of that heat is released as part of L versus carried by the material. Consider the path taken by an infinitesimal fluid element as it moves from face k to face k + 1 due to accretion. Over the course of this adjustment the material changes state and releases some heat. We take this heat to be given by

$$\begin{eqnarray}&&{{dq}}_{k}=-\delta t{\epsilon }_{k,{\rm{I}}}\displaystyle \frac{{{dm}}_{k}}{{{dm}}_{\mathrm{pass},k}},\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{{dm}}_{\mathrm{pass},k}={{dm}}_{\mathrm{start},k}+\delta t\left(\max (0,{F}_{m,k+1})-\min (0,{F}_{m,k})\right)\end{eqnarray} \tag{ 25 }$$

is the total mass of material that at any point during stage I was inside cell k. This evenly distributes the heating that occurs in a cell over all material that starts in, ends in, or passes through the cell. When the thermal time is long, however, the fluid does not have a chance to release all of this heat before it has finished crossing the cell. We estimate the fraction of this heat it releases as the leak fraction

$$\begin{eqnarray}&&{f}_{\mathrm{leak},k}\equiv \min \left(1,\displaystyle \frac{{\tau }_{{\rm{m}},k}}{{\tau }_{\mathrm{th},k}}\right).\end{eqnarray} \tag{ 26 }$$

This parameterizes the extent to which material flowing through a cell follows the implied energy gradient (f_leak,k ≈ 1) versus evolving adiabatically (f_leak,k ≪ 1).

We define δE_k,j to be the amount of energy that the material that ends up in cell j had not leaked by the time it reached face k. This is zero for k = 1 and for k > j, and for all other faces it is given by

$$\begin{eqnarray}&&\delta {E}_{k+1,j}=(1-{f}_{\mathrm{leak},k})\left(\delta {E}_{k,j}+{{dq}}_{k}{{ \mathcal M }}_{k,j}\right),\end{eqnarray} \tag{ 27 }$$

where ${{ \mathcal M }}_{k,j}$ is the amount of material that ends in cell j that passes through cell k during stage I. The heat that is actually released in cell k is then given by

$$\begin{eqnarray}&&{\epsilon }_{\dot{M},k}=\displaystyle \frac{\delta {E}_{k,k}+{\displaystyle \sum }_{j\ne k}{f}_{\mathrm{leak},k}\left(\delta {E}_{k,j}+{{dq}}_{k}{{ \mathcal M }}_{k,j}\right)}{\delta {{tdm}}_{k}}.\end{eqnarray} \tag{ 28 }$$

This is our new corrective source term. The same procedure may be used in the case of mass loss, but with δE_N+1,j = 0 instead of δE_1,j and −1 in the subscript on the left-hand side of Equation (27) rather than +1. In the limit of long thermal times most heat is retained and the resulting evolution is adiabatic. In the limit of short thermal times we recover the results of Paper III.

Along the lines of Equation (17), the energy change of cell k during stage II may now be written in flux-conservative form as

$$\begin{eqnarray}&&{E}_{k}-{E}_{\mathrm{mid},k}=\delta t\left[{{dm}}_{k}({\epsilon }_{\dot{M},k}+\epsilon )+{F}_{e,k+1,\mathrm{II}}-{F}_{e,k,\mathrm{II}}\right],\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{F}_{e,k,\mathrm{II}}={L}_{k}+{p}_{\mathrm{face},k}{{ \mathcal A }}_{k}{\hat{v}}_{k,\mathrm{II}}\end{eqnarray} \tag{ 30 }$$

and

$$\begin{eqnarray}&&{\hat{v}}_{k,\mathrm{II}}\equiv \displaystyle \frac{{r}_{k}-{r}_{\mathrm{mid},k}}{\delta t}.\end{eqnarray} \tag{ 31 }$$

MESA solves Equation (29) when using the dedt-form with mass changes. Because the time evolution is implicit, all sources are evaluated at the end of stage II. Nuclear burning is evaluated as if material spends the whole step in the cell in which it ends stage I. While this is usually a good approximation, it may break down when $\delta t\dot{M}$ is large relative to the masses of cells in burning regions.

When $\dot{M}$ ≥ 0, the above procedure for redistributing energy is conservative, so that

$$\begin{eqnarray}&&\displaystyle \sum _{k}{\epsilon }_{k,{\rm{I}}}\delta {{tdm}}_{k}+\displaystyle \sum _{k}{\epsilon }_{\dot{M},k}\delta {{tdm}}_{k}=0.\end{eqnarray} \tag{ 32 }$$

The first sum represents the effect of stage I; the second sum represents the effect of stage II. When $\dot{M}$ < 0, the equality is instead

$$\begin{eqnarray}&&\displaystyle \sum _{k}{\epsilon }_{k,{\rm{I}}}\delta {{tdm}}_{k}+\displaystyle \sum _{k}{\epsilon }_{\dot{M},k}\delta {{tdm}}_{k}=-\delta {E}_{\mathrm{1,0}},\end{eqnarray} \tag{ 33 }$$

with the additional term δE_1,0 being the energy carried out of the star. This term is explicitly accounted for in the MESA energy budget, so in both cases energy is properly accounted for.

To examine the behavior of ${\epsilon }_{\dot{M}}$, we modeled a 0.33 M_⊙ WD accreting He and N at a rate of 10⁻⁵ M_⊙ yr⁻¹. The first panel of Figure 26 shows a profile of ${\epsilon }_{\dot{M}}$ along with the analytic accretion heating calculations of Townsley & Bildsten (2004). For the most part the two agree closely. The second panel shows the mass integral of the same inward from the surface, while the third shows their ratio. Around x_m ≈ 10²⁶ g, τ_th becomes long relative to τ_m, and the two prescriptions differ because that of Townsley & Bildsten (2004) is only applicable where τ_th ≪ τ_m. The ${\epsilon }_{\dot{M}}$ term handles both limits.

Zoom In Zoom Out Reset image size

To demonstrate improved energy conservation during mass changes with the dedt-form, we model accretion onto a 0.3 M_⊙ He WD with an initial $\mathrm{log}({T}_{{\rm{c}}}/{\rm{K}})=6.7$. The accretion rate is fixed at 10⁻¹⁰ M_⊙ yr⁻¹. Nuclear reactions are disabled throughout the run. This is repeated with the dLdm-form. The relative cumulative error in energy conservation is shown in the top panel of Figure 27, and the relative error in each step is shown in the bottom panel. Near the beginning of the run there is a period where the dLdm-form performs better; however, at those early times both forms do a good job, conserving energy to 1 part in ∼10⁵. At later times the dedt-form produces less error, staying below 1 part in 10⁴, while the dLdm-form yields cumulative errors greater than 1 part in 10³.

Zoom In Zoom Out Reset image size

We prefer using the dedt-form because of its improved energy conservation and handling of long thermal times. We now explore the consequences for stellar evolution of these different prescriptions. The top panel of Figure 28 shows the luminosity for the same case as Figure 27, again for both the dedt-form and the dLdm-form. The difference is shown in the bottom panel. The largest differences are at early times as the models adjust to the accretion. After that, both yield results similar to a few percent, and the relative difference only improves as the luminosity increases. Figure 29 shows ${T}_{{\rm{c}}}$ for the same two runs as a function of time. The differences are small. This is because the core lies beneath the region where cell masses are adjusted significantly, so the precise handling of mass changes only matters for the core insofar as the core temperature is sensitive to the luminosity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Mass changes are ubiquitous in binary stellar evolution. To demonstrate the effect of the dedt-form, we model the evolution of a binary system with an 8 M_⊙ primary and a 6.5 M_⊙ secondary with an initial orbital period of 3 days. This is the same example as in Figure 4 of Paper III. The relative cumulative error in energy conservation is shown in the top panel of Figure 30. The model is also run with the dLdm-form. The errors are shown independently for the primary and the secondary. Figure 31 compares the mass transfer rate for this system computed using each energy equation. The differences are typically of order 1%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use the gravity-darkening model of Espinosa Lara & Rieutord (2011, hereafter ELR), where it is assumed that the radiative flux is directed antiparallel to the effective surface gravity. At a point on the stellar surface with polar angle θ, we find the value of the scaled photosphere radius $\tilde{r}=R/{R}_{{\rm{e}}}$ by solving

$$\begin{eqnarray}&&\displaystyle \frac{1}{\tilde{r}}+\displaystyle \frac{{\omega }^{2}}{2}{\tilde{r}}^{2}{\sin }^{2}\theta =1+\displaystyle \frac{{\omega }^{2}}{2}.\end{eqnarray} \tag{ 40 }$$

We then solve

$$\begin{eqnarray}&&\cos \vartheta +\mathrm{ln}\ \tan (\vartheta /2)=\displaystyle \frac{1}{3}{\omega }^{2}{\tilde{r}}^{3}{\cos }^{3}\theta +\cos \theta +\mathrm{ln}\ \tan (\theta /2),\end{eqnarray} \tag{ 41 }$$

for the modified angular variable ϑ. Using ELR Equation (31), we use this value to obtain the local T_eff.

Figure 35 shows the variation of ${\tilde{T}}_{\mathrm{eff}}={T}_{\mathrm{eff}}{[L/(4\pi \sigma {R}_{{\rm{e}}}^{2})]}^{-1/4}$ over a range of θ for a series of curves with different values of ω. When ω = 0, ${\tilde{T}}_{\mathrm{eff}}=1$ for all θ. When ω = 1, ${\tilde{T}}_{\mathrm{eff}}$ varies by nearly a factor of 2 between the pole and the equator.

Zoom In Zoom Out Reset image size

We are interested in the projected—the directional average over the surface along the LOS—T_eff and L. The two parameters governing the problem are ω and the inclination angle, i, of the LOS with respect to the rotation axis of the star: i = 90° when the LOS is in the plane of the equator. We denote the LOS unit vector $\hat{l}(i)$ and the projected surface area Σ_proj. Figure 36 shows a grid of Roche equipotential surfaces for different values of ω and i. The color describes the variation of ${\tilde{T}}_{\mathrm{eff}}$ over the surface.

Zoom In Zoom Out Reset image size

To calculate the luminosity projected along the LOS, L_proj, requires the surface integral

$$\begin{eqnarray}&&{L}_{\mathrm{proj}}=4\times \mathop{\iint }\limits_{d{\boldsymbol{\Sigma }}\cdot \hat{l}\gt 0}{Fd}{\boldsymbol{\Sigma }}\cdot \hat{l},\end{eqnarray} \tag{ 42 }$$

where only the flux projected toward the observer, i.e., $d{\boldsymbol{\Sigma }}\cdot \hat{l}\gt 0$, is kept. The emergent specific intensity from each surface element is assumed to be isotropic. Once L_proj and Σ_proj are known, the projected T_eff can be obtained from the Stefan–Boltzmann law

$$\begin{eqnarray}&&{T}_{\mathrm{eff},\mathrm{proj}}={\left(\displaystyle \frac{{L}_{\mathrm{proj}}}{\sigma {{\rm{\Sigma }}}_{\mathrm{proj}}}\right)}^{1/4}.\end{eqnarray} \tag{ 43 }$$

As noted by Georgy et al. (2014), the ratio of L_proj/L and, likewise, the ratio of T_eff,proj/T_eff are geometric factors that depend only on ω and i. We define gravity-darkening coefficients

$$\begin{eqnarray}&&{C}_{L}(\omega ,i)\equiv {L}_{\mathrm{proj}}(\omega ,i)/L\end{eqnarray} \tag{ 44 }$$

and

$$\begin{eqnarray}&&{C}_{T}(\omega ,i)\equiv {T}_{\mathrm{eff},\mathrm{proj}}(\omega ,i)/{T}_{\mathrm{eff}}.\end{eqnarray} \tag{ 45 }$$

C_T and C_L are tabulated in MESA for the valid domain of ω and i; values are readily obtained via bicubic interpolation. The projected L and T_eff as viewed from the pole and the equator can be output in the MESA history file; other inclination angles can be accessed via run_star_extras.

Figure 37 shows how C_T and C_L vary over the (ω,i)-plane. These can be compared with the top panels of Figure 2 in Georgy et al. (2014). The variation of C_L is greater than that of C_T owing to the one-fourth-power relationship between L and T_eff. At ω = 1, where the geometric factors are the largest, C_L varies from −20% to +50%, while C_T varies by only about ±2.5%.

Zoom In Zoom Out Reset image size

Figure 37 shows a slight, but noticeable, decrease in C_T for ω ≃ 1, whereas the comparable figures from Georgy et al. (2014) do not. When we calculate the coefficients using oblate spheroids instead of Roche equipotential surfaces, we see no such decrease in C_T for near-critical rotation.

Figure 38 demonstrates the effect of gravity darkening on three of the 3 M_⊙ tracks from Figure 34 in the H-R diagram. In each panel of Figure 38 the nonrotating track is shown for reference as the dotted line. The magnitude of the gravity-darkening effect increases with ω and is more substantial for L than for T_eff. Relative to the intrinsic track, the polar projection is brighter and hotter, and the equatorial projection is cooler and fainter.

Zoom In Zoom Out Reset image size

The control structure for blending the various EOS sources in Figure 50 considers a particular source EOS (PTEH, PC, HELM, etc.) to determine what fraction f of the final result comes from that source. A recursive calling structure is used:

• Level 1: If f_PTEH = 0, then return result of level 2. If f_PTEH = 1, then evaluate PTEH and return. Otherwise, call level 2, blend result with that of PTEH, and return.
• Level 2: If f_PC = 0, then return result of level 3. If f_PC = 1, then evaluate PC and return result to level 1. Otherwise, call level 3, blend result with that of PC, and return the blended result to level 1.
• Level 3: If f_PTEH,Z = 0 (i.e., Z ≲ 0.04), then call level 4 with DT2 and return the result to level 2. If f_PTEH,Z = 1 (i.e., Z ≳ 0.04), then call level 4 with PTEH and return the result to level 2. Otherwise (i.e., Z near 0.04), call level 4 with each of DT2 and PTEH, blend the result in Z and return the blended result to level 2.
• Level 4: If f_PTEH/DT2 = 0, then return result of level 5. If f_PTEH/DT2 = 1, then evaluate PTEH/DT2²⁷ and return result to level 3. Otherwise, call level 5, blend the result with that of PTEH/DT2, and return the blended result to level 3.
• Level 5: If f_ELM = 0, then return result of level 6. If f_ELM = 1, then evaluate ELM and return result to level 4. Otherwise, call level 6, blend the result with that of ELM, and return the blended result to level 4.
• Level 6: Evaluate HELM and return result to level 5.

The blends use a quintic polynomial with zero slope at boundaries. The partial derivatives of the blend polynomial are included in the calculation of the final result to maintain high numerical consistency in the blend region and across EOS region boundaries.

There remain challenges to providing a broad coverage EOS given the current need to combine multiple sources, some of which may not provide the necessary thermodynamic information. For example, when the $\mathrm{log}(\rho$/g cm⁻³) blend region extended from 3.0 to 3.1, a negative χ_ρ,gas resulted because $\mathrm{log}({p}_{\mathrm{gas}}$/erg cm⁻³) from DT2 were slightly greater than $\mathrm{log}({p}_{\mathrm{gas}}$/erg cm⁻³) from ELM. This could give a drop in p_gas as $\mathrm{log}(\rho$/g cm⁻³) transitions from DT2 to ELM. Extending the blend region from 2.98 to 3.12 is enough to ensure χ_ρ,gas > 0 in the DT2-to-ELM transition.

To check the numerical accuracy of partials using the new options, we have compared their results with iteratively acquired high-precision numerical derivatives (Ridders 1982; Press et al. 1992). For example, Figure 52 shows the relative difference between the eos derivative and a Richardson iterative numerical derivative across the ρ–T plane for χ_ρ,gas. The right panel shows that new options give relative errors of ≃10⁻⁷, while the left panel shows that the previous options have larger errors, particularly in the OPAL/SCVH regions. As mentioned in Section 3, those larger errors limit the ability of the MESAstar Newton–Raphson solver to reduce residuals, and that in turn can lead to increased errors in numerical energy conservation.

Zoom In Zoom Out Reset image size

The first law of thermodynamics is an exact differential, which implies

$$\begin{eqnarray}\begin{array}{rcl}p & = & {\left.{\rho }^{2}\displaystyle \frac{\partial e}{\partial \rho }\right|}_{T,{Y}_{i}}+{\left.T\displaystyle \frac{\partial p}{\partial T}\right|}_{\rho ,{Y}_{i}},\\ {\left.\displaystyle \frac{\partial e}{\partial T}\right|}_{\rho ,{Y}_{i}} & = & {\left.T\displaystyle \frac{\partial s}{\partial T}\right|}_{\rho ,{Y}_{i}},\\ {\left.-\displaystyle \frac{\partial s}{\partial \rho }\right|}_{T,{Y}_{i}} & = & \displaystyle \frac{1}{{\rho }^{2}}{\left.\displaystyle \frac{\partial p}{\partial T}\right|}_{\rho ,{Y}_{i}}.\end{array}\end{eqnarray} \tag{ 47 }$$

An EOS is thermodynamically consistent if these relations are satisfied. Thermodynamic inconsistency may manifest itself as an artificial buildup or decay of the entropy during what should be an adiabatic flow. Models that are sensitive to the entropy may suffer inaccuracies if thermodynamic consistency is systematically violated over sufficiently long timescales. Equation (47) may be recast in a form suitable for evaluating numerical inconsistencies

$$\begin{eqnarray}\begin{array}{rcl}{\mathtt{dpe}} & = & {\left.\displaystyle \frac{{\rho }^{2}}{p}\displaystyle \frac{\partial e}{\partial \rho }\right|}_{T,{Y}_{i}}+{\left.\displaystyle \frac{T}{p}\displaystyle \frac{\partial p}{\partial T}\right|}_{\rho ,{Y}_{i}}-1,\\ {\mathtt{dse}} & = & \left.T{\left(\displaystyle \frac{\partial s}{\partial T}\right|}_{\rho ,{Y}_{i}}/{\left.\displaystyle \frac{\partial e}{\partial T}\right|}_{\rho ,{Y}_{i}}\right)-1,\\ {\mathtt{dsp}} & = & \left.-{\rho }^{2}{\left(\displaystyle \frac{\partial s}{\partial \rho }\right|}_{T,{Y}_{i}}/{\left.\displaystyle \frac{\partial p}{\partial T}\right|}_{\rho ,{Y}_{i}}\right)-1.\end{array}\end{eqnarray} \tag{ 48 }$$

Ideally, dpe, dse, and dsp are zero. Figure 53 shows the first thermodynamic consistency quantity, dpe, across the ρ–T plane for an older (MESA r8845) and current version of eos. The other two thermodynamic consistency metrics, dse and dsp, show similar magnitudes. In general, the thermodynamic consistency with DT2 and ELM is reduced relative to the older options directly using OPAL/SCVH and HELM. This is because the thermodynamic consistency relations can only be approximated by bicubic splines. If possible, Hermite interpolation of the Helmholtz free energy would make use of partial derivatives from the EOS and guarantee thermodynamic consistency (e.g., Timmes & Swesty 2000). However, a bicubic Hermite interpolation produces discontinuities in second-derivative quantities (e.g., $\partial {p}_{\mathrm{gas}}/\partial \rho {| }_{T,{Y}_{i}}$) that are problematic for MESA's Newton–Raphson solver, and a biquintic Hermite interpolation requires partial derivatives that are unavailable from some constituents of the EOS patchwork. So for now, we compromise by providing the previous options that give better thermodynamic consistency and the new options that provide better numerical accuracy of partials. The impact of the thermodynamic inconsistencies of the new approach must be evaluated on a case-by-case basis.

Zoom In Zoom Out Reset image size

The Joint Institute for Nuclear Astrophysics (JINA) REACLIB library, which provides the default nuclear reaction rates for MESA, has been updated from the jina_reaclib_results_v2.2 snapshot to the default snapshot.²⁸ This update includes changes to the fitting formula for a few neutron capture rates that previously returned erroneous values at T ≳ 8 × 10⁹ K. We have modified the default snapshot to include a missing, temperature-independent ${}^{26}\mathrm{Al}{\to }^{26}\mathrm{Mg}$ weak reaction rate. Reaction rates between the ²⁶Al ground and meta-stable states now use Gupta & Meyer (2001). Nuclear partition functions use JINA winvne_v2.0.dat table²⁹ and JINA masslib_library_5.data³⁰ provides atomic masses.

A cell's temperature may exceed $\mathrm{log}(T/{\rm{K}})$ = 10.0 in shocks and explosive burning, which is beyond the range of validity for the fits to the reaction rates and the partition functions. Previously, MESA extrapolated for $\mathrm{log}(T/{\rm{K}})$ > 10.0, leading to erroneous reaction rates. Reaction rates are now set equal to their $\mathrm{log}(T/{\rm{K}})$ = 10.0 values when $\mathrm{log}(T/{\rm{K}})$ > 10.0.

MESA now includes the option to use the electron-capture and β-decay rates from Suzuki et al. (2016), which cover sd-shell nuclei with A = 17–28. The primary application for these tables is the evolution of high-density oxygen-neon cores (e.g., Miyaji et al. 1980; Miyaji & Nomoto 1987; Jones et al. 2013). Compared to the on-the-fly weak-rate approach described in Paper III, these tabulated rates are less computationally expensive and more accurate at high temperatures (T ≳ 10⁹ K) but less accurate at low temperatures (T ≲ 10⁸ K) and do not allow explicit updating of nuclear physics.

The plasma coupling parameter of two reactants Γ_i,j and the ion sphere radius a_i are

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{i,j}=\displaystyle \frac{{Z}_{i}{Z}_{j}{e}^{2}}{0.5\left({a}_{i}+{a}_{j}\right){k}_{{\rm{B}}}T},\quad {a}_{i}={\left[\displaystyle \frac{3{Z}_{i}}{4\pi \left({Z}_{i}{n}_{i}+{Z}_{j}{n}_{j}\right)}\right]}^{1/3},\end{eqnarray} \tag{ 49 }$$

where Z_i is the atomic charge of isotope i, e is the electron charge, ${k}_{{\rm{B}}}$ is the Boltzmann constant, and n_i is the ion number density of species i. MESA applies screening factors to correct nuclear reaction rates for plasma interactions (e.g., Salpeter 1954). Previous versions defaulted to using one set of expressions for the weak screening regime (Γ_i,j ≤ 0.3; Dewitt et al. 1973; Graboske et al. 1973) and another set of expressions for the strong screening regime (Γ_i,j ≥ 0.8, Alastuey & Jancovici 1978; Itoh et al. 1979). A linear blend of the weak and strong screening factors is used in the 0.3 < Γ_i,j < 0.8 intermediate regime. Silicon-burning reactions in the cores of massive stars often operate in this intermediate regime. The numerical blending provides a smooth and continuous function that equals the nonblended values and their derivatives at the edges. A new default for screening, based on Chugunov et al. (2007), includes a physical parameterization for the intermediate screening regime and reduces to the familiar weak and strong limits at small and large Γ values. We extend the Chugunov et al. (2007) one-component plasma results to a multicomponent plasma following Itoh et al. (1979), where the Z_i are replaced with the average charge $\bar{Z}$.

Figure 54 compares three MESA nuclear reaction rate screening options on the evolution of a 25 M_⊙ model. The differences are small from H burning to the onset of core collapse. However, the Chugunov et al. (2007) implementation takes ≈20% fewer time steps, retries, and backups, which indicates a numerically smoother solution.

Zoom In Zoom Out Reset image size