MODULES FOR EXPERIMENTS IN STELLAR ASTROPHYSICS (MESA): BINARIES, PULSATIONS, AND EXPLOSIONS

A binary system is initialized by specifying the components and either the orbital period ${P}_{\mathrm{orb}}$ or separation a. Each component can be a point mass or a stellar model. The initial model(s) are provided by a saved MESA model or a zero-age main-sequence (ZAMS) specification. For stellar models including rotation, the initial rotational velocities of each component can be explicitly defined, or set such that the star is synchronized to the orbit at the beginning of the evolution. The orbital angular momentum of the system is

$$\begin{eqnarray}&&{J}_{\mathrm{orb}}={M}_{1}{M}_{2}\sqrt{\displaystyle \frac{{Ga}}{{M}_{1}+{M}_{2}}},\end{eqnarray} \tag{ 1 }$$

where M₁ and M₂ are the stellar masses. Evolution of M₁, M₂, and ${J}_{\mathrm{orb}}$ is used to update a using Equation (1). Masses can be modified both by Roche lobe overflow (RLOF) and winds. The total time derivatives of the component masses are given by

$$\begin{eqnarray}&&{\dot{M}}_{1}={\dot{M}}_{1,{\rm{w}}}+{\dot{M}}_{\mathrm{RLOF}},\quad {\dot{M}}_{2}={\dot{M}}_{2,{\rm{w}}}-{f}_{\mathrm{mt}}{\dot{M}}_{\mathrm{RLOF}},\end{eqnarray} \tag{ 2 }$$

where M₁ is the donor mass and M₂ the accretor mass. The stellar wind mass loss rates are ${\dot{M}}_{1,{\rm{w}}}$ and ${\dot{M}}_{2,{\rm{w}}}$ (see Paper I and Paper II) and ${\dot{M}}_{\mathrm{RLOF}}$ is the mass transfer rate from RLOF, all defined as negative. The factor ${f}_{\mathrm{mt}}$ represents the efficiency of accretion and can be used to limit accretion to the Eddington rate ${\dot{M}}_{\mathrm{Edd}}$.

To compute the rate of change of orbital angular momentum, we consider the contribution of gravitational waves, mass loss, magnetic braking, and spin–orbit (LS) coupling

$$\begin{eqnarray}&&{\dot{J}}_{\mathrm{orb}}={\dot{J}}_{\mathrm{gr}}+{\dot{J}}_{\mathrm{ml}}+{\dot{J}}_{\mathrm{mb}}+{\dot{J}}_{\mathrm{ls}},\end{eqnarray} \tag{ 3 }$$

from which the change in orbital angular momentum in one step is calculated as ${\rm{\Delta }}{J}_{\mathrm{orb}}={\dot{J}}_{\mathrm{orb}}\delta t$, where $\delta t$ is the timestep. Unless models with stellar rotation are being used, the ${\dot{J}}_{\mathrm{ls}}$ term is equal to zero, and the contribution of the individual spins of each star is not directly considered. On the other hand, the ${\dot{J}}_{\mathrm{mb}}$ term implicitly assumes a strong tide that keeps the orbit synchronized. The simultaneous usage of ${\dot{J}}_{\mathrm{mb}}$ with stellar rotation is not consistent (see Section 2.2.4). We now describe how these terms are computed.

2.2.1. Gravitational Wave Radiation

Very compact binaries can experience significant orbital decay due to the emission of gravitational waves. Observations of the Hulse–Taylor binary pulsar over three decades (Weisberg & Taylor 2005) and of the double pulsar (Kramer et al. 2006) have tested the predicted effect from general relativity to a high precision. The angular momentum loss from gravitational waves is

$$\begin{eqnarray}&&{\dot{J}}_{\mathrm{gr}}=-\displaystyle \frac{32}{5{c}^{5}}{\left(\displaystyle \frac{2\pi G}{{P}_{\mathrm{orb}}}\right)}^{7/3}\displaystyle \frac{{({M}_{1}{M}_{2})}^{2}}{{({M}_{1}+{M}_{2})}^{2/3}}.\end{eqnarray} \tag{ 4 }$$

2.2.2. Mass Loss

We assume the mass lost in a stellar wind has the specific orbital angular momentum of its star. For the case of inefficient mass transfer, angular momentum loss follows Soberman et al. (1997), where fixed fractions of the transferred mass are lost either as a fast isotropic wind from each star or a circumbinary toroid with a given radius:

$$\begin{eqnarray}{\dot{J}}_{\mathrm{ml}} & = & \left[({\dot{M}}_{1,{\rm{w}}}+{\alpha }_{\mathrm{mt}}{\dot{M}}_{\mathrm{RLOF}}){M}_{1}^{2}\right.\\ & & +\left.({\dot{M}}_{2,{\rm{w}}}+{\beta }_{\mathrm{mt}}{\dot{M}}_{\mathrm{RLOF}}){M}_{2}^{2}\right]\displaystyle \frac{{a}^{2}}{{({M}_{1}+{M}_{2})}^{2}}\displaystyle \frac{2\pi }{{P}_{\mathrm{orb}}}\\ & & +\ {\gamma }_{\mathrm{mt}}{\delta }_{\mathrm{mt}}{\dot{M}}_{\mathrm{RLOF}}\sqrt{G({M}_{1}+{M}_{2})a},\end{eqnarray} \tag{ 5 }$$

where ${\alpha }_{\mathrm{mt}}$, ${\beta }_{\mathrm{mt}}$, and ${\delta }_{\mathrm{mt}}$ are respectively the fractions of mass transferred that is lost from the vicinity of the donor, accretor and circumbinary toroid, and ${\gamma }_{\mathrm{mt}}^{2}a$ is the radius of the toroid. Ignoring winds, the efficiency of mass transfer is then given by ${f}_{\mathrm{mt}}=1-{\alpha }_{\mathrm{mt}}-{\beta }_{\mathrm{mt}}-{\delta }_{\mathrm{mt}}$. When accretion is limited to ${\dot{M}}_{\mathrm{Edd}}$, efficiency of accretion is given by

$$\begin{eqnarray}&&{f}_{\mathrm{mt}}=\mathrm{min}(1-{\alpha }_{\mathrm{mt}}-{\beta }_{\mathrm{mt}}-{\delta }_{\mathrm{mt}},| {\dot{M}}_{\mathrm{Edd}}/{\dot{M}}_{\mathrm{RLOF}}| ),\end{eqnarray} \tag{ 6 }$$

and the additional mass being lost is added to the ${\beta }_{\mathrm{mt}}{\dot{M}}_{\mathrm{RLOF}}$ term in Equation (5), i.e., it is assumed to leave the system carrying the specific orbital angular momentum of the accretor.

2.2.3. Spin–Orbit Coupling

Tidal interaction and mass transfer can significantly modify the spin angular momentum of the stars in a binary system, acting as both sources and sinks for orbital angular momentum. The impact spin–orbit interactions have on orbital evolution depends on the orbital separation and the mass ratio, with the effect being greater for tighter orbits and uneven masses. The corresponding contribution to ${\dot{J}}_{\mathrm{orb}}$ is computed by demanding conservation of the total angular momentum, accounting for losses due to the other ${\dot{J}}_{\mathrm{orb}}$ mechanisms and loss of stellar angular momentum due to winds.

In a fully conservative system, the change in orbital angular momentum in one timestep is $\delta {J}_{\mathrm{orb}}=-\delta {S}_{1}-\delta {S}_{2}$, where $\delta {S}_{1}$ and $\delta {S}_{2}$ are the changes in spin angular momenta. This needs to be corrected if mass loss is included, as winds take away angular momentum from the system. If ${S}_{1,\mathrm{lost}}$ and ${S}_{2,\mathrm{lost}}$ are the amounts of spin angular momentum removed in a step from each star due to mass loss (including winds and RLOF),

$$\begin{eqnarray}&&{\dot{J}}_{\mathrm{ls}}=\displaystyle \frac{-1}{\delta t}\left(\delta {S}_{1}-{S}_{1,\mathrm{lost}}\displaystyle \frac{{\dot{M}}_{1,{\rm{w}}}}{{\dot{M}}_{1}}+\delta {S}_{2}-{S}_{2,\mathrm{lost}}\right),\end{eqnarray} \tag{ 7 }$$

where the additional factor for the donor accounts for mass lost from the system, ignoring mass loss due to mass transfer. In the absence of RLOF this equation becomes symmetric between both stars, as then ${\dot{M}}_{1,{\rm{w}}}/{\dot{M}}_{1}=1$.

The form of Equation (7) is independent of how tides and angular momentum accretion work, as it is merely a statement on angular momentum conservation. The details of how we model these processes and their impact on the spin of each component are described in Section 2.4.

2.2.4. Magnetic Braking

The rotational velocities of low mass stars are strongly correlated with their ages (Skumanich 1972). This spin-down arises from the coupling of the stellar wind to a magnetic field. If the star is in a binary system and tidally coupled to the orbit, magnetic braking can provide a very efficient sink for orbital angular momentum (Mestel 1968; Verbunt & Zwaan 1981). We implement this effect following Rappaport et al. (1983), who assumed the star being braked is tidally synchronized:

$$\begin{eqnarray}&&{\dot{J}}_{\mathrm{mb}}=-6.82\;\times \;{10}^{34}\left(\displaystyle \frac{{M}_{1}}{{{\rm{M}}}_{\odot }}\right){\left(\displaystyle \frac{{R}_{1}}{{{\rm{R}}}_{\odot }}\right)}^{{\gamma }_{\mathrm{mb}}}{\left(\displaystyle \frac{1\ d}{{P}_{\mathrm{orb}}}\right)}^{3}\ [\mathrm{dyn}\ \mathrm{cm}],\end{eqnarray} \tag{ 8 }$$

where in the simplest approximation ${\gamma }_{\mathrm{mb}}$ = 4 (Verbunt & Zwaan 1981). A similar contribution from the accretor can be included. As tidal synchronization is assumed, this formulation is incompatible with the use of LS coupling.

It is normally assumed that once a star becomes fully convective, the dynamo process that regenerates the field will stop working or at least behave in a qualitatively different manner. Similarly, magnetic fields in stars with radiative envelopes are of a significantly different nature than those seen in stars with convective envelopes, and there is no simple way to predict even the presence of magnetism (Donati & Landstreet 2009). By default, MESAbinary only accounts for magnetic braking as long as the star being braked has a convective envelope and a radiative core, though the process might still operate outside of these conditions.

Close binary stars are defined as systems tight enough to interact through mass transfer, with the most important mechanism being RLOF. This process is commonly modeled in 1D by considering the spherical-equivalent Roche lobe radius ${R}_{\mathrm{RL}}$ of each object (Eggleton 1983)

$$\begin{eqnarray}&&{R}_{\mathrm{RL},j}=\displaystyle \frac{0.49{q}_{j}^{2/3}}{0.6{q}_{j}^{2/3}+\mathrm{ln}\left(1+{q}_{j}^{1/3}\right)}a,\end{eqnarray} \tag{ 9 }$$

where j is the index identifying each star, ${q}_{1}={M}_{1}/{M}_{2}$ and ${q}_{2}={M}_{2}/{M}_{1}$. This fit is correct up to a few percent for the full range of mass ratios, $0\lt {q}_{j}\lt \infty$. Mass transfer occurs then when the radius of a star approaches or exceeds ${R}_{\mathrm{RL}}$. Depending on several factors, once a star begins RLOF the ensuing mass transfer phase can proceed on a nuclear, thermal, or dynamical timescale.

The stability of mass transfer is normally understood in terms of mass–radius relationships (e.g., Soberman et al. 1997; Tout et al. 1997),

$$\begin{eqnarray}&&{\zeta }_{\mathrm{eq}}={\left(\displaystyle \frac{d\mathrm{ln}{R}_{1}}{d\mathrm{ln}{M}_{1}}\right)}_{\mathrm{eq}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\zeta }_{\mathrm{ad}}={\left(\displaystyle \frac{d\mathrm{ln}{R}_{1}}{d\mathrm{ln}{M}_{1}}\right)}_{\mathrm{ad}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\zeta }_{\mathrm{RL}}=\displaystyle \frac{d\mathrm{ln}{R}_{\mathrm{RL},1}}{d\mathrm{ln}{M}_{1}}.\end{eqnarray} \tag{ 12 }$$

Here, ${\zeta }_{\mathrm{eq}}$ gives the radial response of the donor to mass loss when it happens slowly enough for the star to remain in thermal equilibrium. When mass loss proceeds on a timescale much shorter than the thermal timescale of the star, but still slow enough for the star to retain hydrostatic equilibrium then the radial response will be given by ${\zeta }_{\mathrm{ad}}$. The dependency of the Roche lobe radius on mass transfer is encoded in ${\zeta }_{\mathrm{RL}}$. In general $\zeta =d\mathrm{ln}{R}_{1}/d\mathrm{ln}{M}_{1}$ is a function of ${\dot{M}}_{\mathrm{RLOF}}$, so requiring $\zeta ={\zeta }_{\mathrm{RL}}$ will determine the value of ${\dot{M}}_{\mathrm{RLOF}}$. If an overflowing star satisfies ${\zeta }_{\mathrm{eq}}\gt {\zeta }_{\mathrm{RL}}$, then it can remain inside its Roche lobe by transferring mass while retaining thermal equilibrium. If on the contrary ${\zeta }_{\mathrm{ad}}\gt {\zeta }_{\mathrm{RL}}\gt {\zeta }_{\mathrm{eq}}$, mass transfer will proceed on a thermal timescale, while for the extreme case ${\zeta }_{\mathrm{RL}}\gt {\zeta }_{\mathrm{ad}}$ the star will depart from hydrostatic equilibrium and the process will be dynamical. MESA cannot model common envelope or contact binaries.

MESAbinary provides both explicit and implicit methods to compute mass transfer rates. An explicit computation sets the value of ${\dot{M}}_{\mathrm{RLOF}}$ at the start of a step, while an implicit one begins with a guess for ${\dot{M}}_{\mathrm{RLOF}}$ and iterates until the required tolerance is reached. The composition of accreted material is set to that of the donor surface, and the specific entropy of accreted material is the same as the surface of the accretor. In models with rotation the specific angular momentum of accreted material is described in Section 2.4.

2.3.1. Explicit Methods

MESAbinary implements two mass transfer schemes: the model of Ritter (1988) which we refer to as the Ritter scheme and Kolb & Ritter (1990) which we refer to as the Kolb scheme. We use the mass ratio q₂ consistent with the Ritter scheme.

Ritter scheme: Stars have extended atmospheres therefore RLOF can take place through the L1 point even when ${R}_{1}\lt {R}_{\mathrm{RL},1}$. Ritter (1988) estimated the mass transfer rate for this case as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{RLOF}}=-{\dot{M}}_{0}\mathrm{exp}\left(\displaystyle \frac{{R}_{1}-{R}_{\mathrm{RL},1}}{{H}_{{\rm{P}},1}/\gamma ({q}_{2})}\right),\end{eqnarray} \tag{ 13 }$$

where ${H}_{{\rm{P}},1}$ is the pressure scale height at the photosphere of the donor and

$$\begin{eqnarray}&&{\dot{M}}_{0}=\displaystyle \frac{2\pi }{\mathrm{exp}(1/2)}{F}_{1}({q}_{2})\displaystyle \frac{{R}_{\mathrm{RL},1}^{3}}{{{GM}}_{1}}{\left(\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{eff}}}{{m}_{{\rm{p}}}\ {\mu }_{\mathrm{ph}}}\right)}^{3/2}{\rho }_{\mathrm{ph}},\end{eqnarray} \tag{ 14 }$$

where ${m}_{{\rm{p}}}$ is the proton mass, ${T}_{\mathrm{eff}}$ is the effective temperature of the donor, and ${\mu }_{\mathrm{ph}}$ and ${\rho }_{\mathrm{ph}}$ are the mean molecular weight and density at its photosphere. The two fitting functions are

$$\begin{eqnarray}&&{F}_{1}({q}_{2})=1.23+0.5\mathrm{log}{q}_{2},\qquad 0.5\lesssim {q}_{2}\lesssim 10,\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&\gamma ({q}_{2})\;=\\ \left\{\begin{array}{ll}0.954+0.025\mathrm{log}{q}_{2}-0.038{(\mathrm{log}{q}_{2})}^{2}, & \\ 0.04\lesssim {q}_{2}\leqslant 10.954+0.039\mathrm{log}{q}_{2}+0.114{(\mathrm{log}{q}_{2})}^{2}, & \\ \quad 1\leqslant {q}_{2}\lesssim 20. & \end{array}\right.\;\end{eqnarray} \tag{ 16 }$$

Outside the ranges of validity, ${F}_{1}({q}_{2})$ and $\gamma ({q}_{2})$ are evaluated using the value of q₂ at the edge of their respective ranges.

Kolb scheme: Kolb & Ritter (1990) extended the Ritter scheme in order to cover the case ${R}_{1}\gt {R}_{\mathrm{RL},1}$ according to

$$\begin{eqnarray}{\dot{M}}_{\mathrm{RLOF}} & = & -{\dot{M}}_{0}-2\pi {F}_{1}({q}_{2})\displaystyle \frac{{R}_{\mathrm{RL},1}^{3}}{{{GM}}_{1}}\qquad \qquad \qquad \qquad \\ & & \;\times \;{\displaystyle \int }_{{P}_{\mathrm{ph}}}^{{P}_{\mathrm{RL}}}{{\rm{\Gamma }}}_{1}^{1/2}{\left(\displaystyle \frac{2}{{{\rm{\Gamma }}}_{1}+1}\right)}^{({{\rm{\Gamma }}}_{1}+1)/(2{{\rm{\Gamma }}}_{1}-2)}{\left(\displaystyle \frac{{k}_{{\rm{B}}}T}{{m}_{{\rm{p}}}\mu }\right)}^{1/2}{dP}\end{eqnarray} \tag{ 17 }$$

where ${{\rm{\Gamma }}}_{1}$ is the first adiabatic exponent, and P_ph and P_RL are, respectively, the pressures at the photosphere and at the radius for which ${r}_{1}={R}_{\mathrm{RL},1}$.

2.3.2. Implicit Methods

Explicit schemes exhibit large jumps in ${\dot{M}}_{\mathrm{RLOF}}$ unless the timestep is severely restricted. Therefore, if one needs accurate values of ${\dot{M}}_{\mathrm{RLOF}}$ and stellar radius, this requires use of an implicit scheme. Implicit schemes also allow the calculation these quantities when there is no general closed form formula for ${\dot{M}}_{\mathrm{RLOF}}$.

These implicit methods use a bisection-based root solve to satisfy $| f({\dot{M}}_{\mathrm{RLOF}})| \lt \xi$ at the end of the step, where ξ is a given tolerance. The implicit schemes are then defined by the choice of the function $f({\dot{M}}_{\mathrm{RLOF}})$. For the Ritter and the Kolb scheme the function is chosen as

$$\begin{eqnarray}&&f({\dot{M}}_{\mathrm{RLOF}})=\displaystyle \frac{{\dot{M}}_{\mathrm{end}}-{\dot{M}}_{\mathrm{RLOF}}}{{\dot{M}}_{\mathrm{end}}},\end{eqnarray} \tag{ 18 }$$

with ${\dot{M}}_{\mathrm{end}}$ being the mass transfer rate computed at the end of each iteration.

A different implicit method is also provided. In this case, whenever the donor star overflows its Roche lobe the implicit solver will adjust the mass transfer rate until ${R}_{1}={R}_{\mathrm{RL},1}$ within some tolerance (see, e.g., Whyte & Eggleton 1980; Rappaport et al. 1982, 1983). In this case

$$\begin{eqnarray}&&f({\dot{M}}_{\mathrm{RLOF}})=\displaystyle \frac{2({R}_{1}-{R}_{\mathrm{RL},1})}{{R}_{\mathrm{RL},1}}+\xi ,\end{eqnarray} \tag{ 19 }$$

and if ${\dot{M}}_{\mathrm{RLOF}}$ is below a certain threshold and $f({\dot{M}}_{\mathrm{RLOF}})\lt -\xi$ then the system is assumed to detach and ${\dot{M}}_{\mathrm{RLOF}}$ is set to zero.

To model tidal interaction we adjusted the model of Hut (1981) to include the case of differentially rotating stars. The time evolution of the angular frequency for each component is

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Omega }}}_{i,j}}{{dt}}=\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{orb}}-{{\rm{\Omega }}}_{i,j}}{{\tau }_{\mathrm{sync},j}},\quad \displaystyle \frac{1}{{\tau }_{\mathrm{sync},j}}=\displaystyle \frac{3}{{({q}_{j}{r}_{{\rm{g}},j})}^{2}}{\left(\displaystyle \frac{k}{T}\right)}_{{\rm{c}},j}{\left(\displaystyle \frac{{R}_{j}}{a}\right)}^{6},\end{eqnarray} \tag{ 20 }$$

where $j=1,2$ is the index of each star, ${{\rm{\Omega }}}_{i,j}$ is the angular frequency at the face of cell i toward the surface, ${r}_{{\rm{g}},j}^{2}={I}_{j}/({M}_{j}{R}_{j}^{2})$ is the radius of gyration (with I_j being the moment of inertia of each star), and the ratio of the apsidal motion constant to the viscous dissipation timescale, ${(k/T)}_{{\rm{c}},j}$, is computed as in Hurley et al. (2002). Similarly to Detmers et al. (2008), we assume constant ${\tau }_{\mathrm{sync},j}$ and ${{\rm{\Omega }}}_{\mathrm{orb}}$ through a step and therefore ${\rm{\Delta }}{{\rm{\Omega }}}_{i,j}=[1-\mathrm{exp}(-\delta t/{\tau }_{\mathrm{sync},j})]({{\rm{\Omega }}}_{\mathrm{orb}}-{{\rm{\Omega }}}_{i,j})$. This extension of Hut's work to differentially rotating stars is not formally derived but merely applies his result for solid body rotators independently to each shell. The formulation of Hut (1981) can be recovered from Equation (20), by forcing solid body rotation with a large diffusion coefficient for angular momentum throughout the star. In reality tides would act mostly on the outer layers, and whether the core synchronizes or not depends on the coupling between the core and the envelope.

To compute the specific angular momentum carried by accreted material, we consider the possibility of both ballistic and Keplerian disk mass transfer (e.g., Marsh et al. 2004; de Mink et al. 2013). To distinguish which occurs, we compare the minimum distance of approach of the accretion stream (Lubow & Shu 1975; Ulrich & Burger 1976)¹¹

$$\begin{eqnarray}&&{R}_{\mathrm{min}}=0.0425a{\left({q}_{2}+{q}_{2}^{2}\right)}^{1/4},\qquad 0.0667\leqslant {q}_{2}\leqslant 15\end{eqnarray} \tag{ 21 }$$

to the radius of the accreting star. When outside the range of validity, ${R}_{\mathrm{min}}$ is computed using the value of q₂ at the respective edge. Accretion is assumed to be ballistic whenever ${R}_{2}\gt {R}_{\mathrm{min}}$ and the specific angular momentum is ${(1.7{{GM}}_{2}{R}_{\mathrm{min}})}^{1/2}$. When ${R}_{2}\lt {R}_{\mathrm{min}}$ the specific angular momentum is taken as that of a Keplerian orbit at the surface ${({{GM}}_{2}{R}_{2})}^{1/2}$.

In stars with radiative envelopes accreted material with a high mean molecular weight is expected to mix inwards due to thermohaline mixing, a process that is very sensitive to the μ-gradient (see e.g., Kippenhahn et al. 1980; Cantiello & Langer 2010). Thermohaline mixing is included in MESA (see Paper I). However, as mass with homogeneous composition is added during the accretion process, a jump is produced at the boundary between new and old material. MESAstar computes mixing coefficients explicitly at the start of each step, so this results in thermohaline mixing only operating near this boundary, leading to unphysical compositional staircases. To avoid this issue, we artificially soften the composition gradient in the outer ${({\rm{\Delta }}q)}_{\mathrm{large}}$ fraction of the star by mass. We do this starting at the surface and homogeneously mixing inwards a region of size ${({\rm{\Delta }}q)}_{\mathrm{small}}$. Then, moving toward the center, the process is repeated at each cell while linearly (with respect to mass) reducing the size of the small mixed region such that it is zero after going ${({\rm{\Delta }}q)}_{\mathrm{large}}$ inwards. All the binary models where the accretor is not a point mass are calculated using ${({\rm{\Delta }}q)}_{\mathrm{large}}=0.05$ and ${({\rm{\Delta }}q)}_{\mathrm{small}}=0.03$.

Here we describe tests designed to validate the implementation of the physics described in Section 2.2. We check orbital evolution in the presence of gravitational waves and mass loss by comparing to analytical solutions. We also verify total angular momentum conservation in calculations that include the physics of tides and spinup by accretion. To test for the thermal response of stellar models undergoing mass transfer, we compare MESAbinary results to those from the STARS code (Eggleton 1971; Pols et al. 1995; Stancliffe & Eldridge 2009).

2.6.1. Gravitational Wave Radiation

2.6.2. Inefficient Mass Transfer

An analytical expression for the evolution of orbital separation can be derived if inefficient mass transfer is the only contribution to the angular momentum evolution (Tauris & van den Heuvel 2006, p. 623). We model a $2.5\;{{\rm{M}}}_{\odot }$ main sequence (MS) star together with a $1.4\;{{\rm{M}}}_{\odot }$ point mass with an initial orbital separation of $10\;{{\rm{R}}}_{\odot }$. We choose ${\alpha }_{\mathrm{mt}}=0.03$, ${\beta }_{\mathrm{mt}}=0.95$, ${\delta }_{\mathrm{mt}}=0.01$ and ${\gamma }_{\mathrm{mt}}^{2}=2$, which give a low mass transfer efficiency of ${f}_{\mathrm{mt}}=0.01$. Such a system is representative of the evolution of an intermediate mass X-ray binary (IMXB). The model initiates mass transfer just after the end of the MS, interrupting the evolution of the star through the Hertzsprung gap and producing a low mass white dwarf (WD; ${M}_{\mathrm{He}}=0.289\;{{\rm{M}}}_{\odot }$) with a small amount of hydrogen on its surface. As the WD evolves to the cooling track, it experiences several hydrogen flashes, one of them strong enough to produce an additional phase of RLOF (see Figure 1).

Zoom In Zoom Out Reset image size

Figure 2 shows that MESAbinary computes the orbital evolution to a precision of a few parts in 10⁴. We run this system using both the Ritter and the Kolb implicit schemes to display that under some circumstances the precise choice of mass transfer scheme does not play a big role in the evolution.

Zoom In Zoom Out Reset image size

2.6.3. Spin–Orbit Coupling

We now test angular momentum conservation by ignoring all the mechanisms that remove angular momentum from the binary system. For this purpose we model an $8\;{{\rm{M}}}_{\odot }+6\;{{\rm{M}}}_{\odot }$ binary with rotating components and an initial orbital period of 1.5 days. Due to the short orbital separation we assume the initial spin periods of the two stars are equal to the orbital period. The primary undergoes RLOF during the MS, initiating a phase of mass transfer on a thermal timescale. After transferring just $0.3\;{{\rm{M}}}_{\odot }$ the accretor also fills its Roche lobe, producing a contact system. At this point we terminate the evolution.

Figure 3 shows that spin angular momentum in both components increases during the pre-interaction phase, which is due to both stars expanding on the MS while remaining tidally locked. During RLOF, the secondary is rapidly spun-up, reaching nearly $80\%$ of critical rotation before contact. The calculation of total angular momentum requires the summation of different contributions (orbital angular momentum and spin of both components). Therefore the maximum accuracy to which we can conserve angular momentum is limited by rounding errors. Figure 3 shows that conservation of angular momentum in the run is very close to machine precision.

Zoom In Zoom Out Reset image size

2.6.4. Thermal Response to Mass Loss

The fate of binary systems depends largely on the precise value of $\dot{M}$ during an interaction phase, which depends on the thermal response of the donor star to mass loss. For WDs there is a limited range of accretion rates for stable hydrogen burning (Nomoto et al. 2007; Shen & Bildsten 2007). In MS binaries the evolution of the accretor radius depends on the mass transfer rate, and expansion during the interaction phase can lead to contact or even a merger (Wellstein et al. 2001).

We calculated an $8\;{{\rm{M}}}_{\odot }+6.5\;{{\rm{M}}}_{\odot }$ binary system with an initial orbital period of 1.5 days using both MESAbinary and STARS. To minimize the modeling differences and focus on the thermal response of both components, we use an extremely simplified model that ignores internal mixing (including convective mixing). Under these conditions, the more massive star quickly depletes its central hydrogen and begins shell hydrogen burning, reaching RLOF and undergoing a phase of mass transfer on the thermal timescale. The resulting mass transfer rates are shown in Figure 4. The agreement is very good, despite mass transfer rates being computed in slightly different ways. Masses at detachment show a small difference, with the MESAbinary model terminating mass transfer when ${M}_{1}=0.952\;{{\rm{M}}}_{\odot }$ while the STARS calculation when ${M}_{1}=0.935\;{{\rm{M}}}_{\odot }$. The figure also shows the change in radius of the accreting star, with two prominent peaks at ${R}_{2}/{{\rm{R}}}_{\odot }=4.84,5.34$ for MESAbinary and $4.82,5.28$ for STARS. The larger radius of the MESAbinary model is likely associated to the slightly higher mass transfer rates.

Zoom In Zoom Out Reset image size

Although CVs span a wide range of periods, observations show a lack of systems in the range $2\;\mathrm{hr}\lt {P}_{\mathrm{orb}}\lt 3\;\mathrm{hr}$ (see, for instance, Gänsicke et al. 2009). Such a feature is commonly explained by having an angular momentum loss mechanism "turn off" or become inefficient at some point. The most popular model for such a mechanism is magnetic braking (Rappaport et al. 1983), as the magnetic field of the donor is assumed to change quickly when the star loses enough mass to become fully convective.

We now compare to the results of Howell et al. (2001), who performed a population synthesis study to explore in detail the standard scenario involving magnetic braking. In Figure 5 we show the evolution of mass transfer rates and orbital periods for a set of CV models with different component masses and orbital periods. We run all models using ${\beta }_{\mathrm{mt}}=1$ and ${\gamma }_{\mathrm{mb}}=3$ and magnetic braking is turned off when the donor star becomes fully convective. As an example the system with a $0.9\;{{\rm{M}}}_{\odot }$ donor (left panel in Figure 5) experiences a first phase of mass transfer induced by magnetic braking between ${10}^{7.1}$ and ${10}^{8.3}\;\mathrm{years}$, a non-interacting phase (the gap) between ${10}^{8.3}$ and ${10}^{8.8}\;\mathrm{years}$, and a subsequent phase of mass transfer dominated by gravitational wave radiation, reaching a minimum orbital period of about 1 hr at ${10}^{9.6}\;\mathrm{years}$. As a comparison, for the same model Howell et al. (2001) obtain a first phase of mass transfer between ${10}^{7.3}$ and ${10}^{8.4}\;\mathrm{years}$, the gap occurs between ${10}^{8.4}$ and ${10}^{8.8}\;\mathrm{years}$ and a period minimum is reached at ${10}^{9.4}\;\mathrm{years}$. Figure 5 shows that our CV models spend most time away from the observed period gap.

Zoom In Zoom Out Reset image size

In massive stars, binary interactions have dramatic effects on the evolution of both components. Kippenhahn & Weigert (1967) introduced the term "case A" to refer to a mass transfer phase occurring in systems tight enough such that RLOF starts during the MS. This results in a large amount of mass being transferred on a thermal timescale, followed by a phase of mass transfer that proceeds on the nuclear timescale until the end of core H-burning. An additional phase of thermal-timescale mass transfer then follows (the so-called "case AB"), which strips the donor and produces an almost-naked helium star.

Here we show that MESAbinary can calculate the evolution of massive interacting binaries. We reproduce one of the models from Wellstein et al. (2001), a $16\;{{\rm{M}}}_{\odot }+14\;{{\rm{M}}}_{\odot }$ system with an initial period of 3 days, using the same semiconvection efficiency of ${\alpha }_{\mathrm{sc}}=0.01$. As shown in Figures 6 and 7 this system experiences case A and AB mass transfer, and the accretor becomes a blue supergiant after core hydrogen depletion. The accretor depletes carbon before its donor.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 illustrates the prevalence of both thermohaline mixing and semiconvection in the accreting star. Newly accreted material is efficiently mixed inwards by thermohaline mixing. On the other hand the μ-gradient formed before interaction prevents the convective core from growing, with the efficiency of semiconvection controlling whether or not the star rejuvenates. Due to the choice of inefficient semiconvection, the core remains small, preventing the star from becoming a red supergiant. The star accretes a large amount of CNO-processed and helium-rich material. After being mixed through the envelope this material results in the surface being nitrogen rich and carbon depleted, with a slight enhancement in helium.

The efficiency of mass transfer plays a key role in close binary systems, but the processes by which material is lost from the system are not well-understood. In particular, whenever an accreting star approaches ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit}}=1$, it is uncertain whether accretion can continue, one option being the development of a strong wind that prevents accretion (e.g., Petrovic et al. 2005; Cantiello et al. 2007). Whenever ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit}}$ approaches one, we use an implicit method to iteratively reduce ${\dot{M}}_{2}$ until this ratio falls below a threshold.

Tides counteract the effect of spin-up from accretion. Whether or not an accreting object reaches critical rotation depends on the efficiency of tidal coupling. Here we model a $16\;{{\rm{M}}}_{\odot }+15\;{{\rm{M}}}_{\odot }$ binary system including differential stellar rotation, with an initial orbital period of 3 days and assuming initial orbital synchronization. Langer et al. (2003) argue that turbulent processes in the radiative envelope can significantly enhance tidal strength. They model the same system using the simple estimate for the synchronization timescale for a star with a convective envelope given by Zahn (1977), ${\tau }_{\mathrm{sync},j}=1\;\mathrm{year}\;\times \;{q}_{j}^{2}{(a/R)}^{6}$. For our implicit modeling of stellar winds we use a threshold of $({\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit}})=0.99$.

Figure 8 shows that MESAbinary models using both the Zahn (1977) and Hurley et al. (2002) timescales for tidal coupling. These models experience highly non-conservative phases of mass transfer, corresponding to the accreting star evolving very close to critical rotation. In particular during case AB mass transfer the accretor needs to switch from mass accretion to mass loss in order to remain sub-critical. As expected, the system with the tidal timescale from Zahn (1977) has a significantly higher mass transfer efficiency, and during the first phase of RLOF it only experiences a brief period in which the accretor reaches critical rotation. This is in broad agreement with the model by Langer et al. (2003).

Zoom In Zoom Out Reset image size

MESAbinary performs each evolution step by independently solving the structure of each component and the orbital parameters, using the same timestep $\delta t$ for each. This approach differs from STARS, which simultaneously solves for the structure of both stars and the orbit in a single Newton–Raphson solver. Our choice to solve for each star separately gives a significant amount of flexibility and simplicity, as the examples in this paper demonstrate.

The top-level algorithm for evolving a star is described in Appendix B1 of Paper II. We modified this algorithm to support the new implementation of binary interactions, which is described in detail in the MESA documentation. Additional timestep limits are imposed in MESAbinary that consider relative changes between the radius and Roche lobe radius of both components, the total orbital angular momentum, the orbital separation, and the envelope mass in the donor.

As an illustration of a large-scale calculation using MESAstar and GYRE loosely coupled, Figure 9 plots the instability strips for massive stars on and near the upper MS, for oscillation modes with harmonic degrees ${\ell }=0-3$. These strips are based on a set of 182 evolutionary tracks, each extending from the ZAMS across to a red limit at $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=3.75$, with 101 tracks spanning the initial mass range $2.5\;{{\rm{M}}}_{\odot }\leqslant M\leqslant 25\;{{\rm{M}}}_{\odot }$ in uniform logarithmic increments, and the remaining 81 tracks spanning the mass range $6\;{{\rm{M}}}_{\odot }\leqslant M\leqslant 10\;{{\rm{M}}}_{\odot }$ in uniform linear increments (the latter set is designed to adequately resolve the "fingers" discussed below). OPAL opacity tables are used with the proto-solar abundances from Asplund et al. (2009), and for simplicity we neglect any rotation or mass loss. Convection is modeled with a mixing-length parameter ${\alpha }_{\mathrm{MLT}}$ = 1.5 and an exponential overshoot parameter ${f}_{\mathrm{ov}}=0.024$, and the Schwarzschild stability criterion is assumed.

Zoom In Zoom Out Reset image size

We select points $i={i}_{1},{i}_{2},{i}_{3},\ldots$ along each of the 182 tracks (where i is the timestep index; see Section 6.4 of Paper I), chosen so that ${i}_{1}$ corresponds to the ZAMS,

$$\begin{eqnarray}&&\displaystyle \sum _{i={i}_{1}}^{{i}_{2}-1}{\left\{{\left[\displaystyle \frac{\mathrm{log}({T}_{\mathrm{eff},i+1}/{T}_{\mathrm{eff},i})}{{{\rm{\Delta }}}_{{\rm{T}}}}\right]}^{2}+{\left[\displaystyle \frac{\mathrm{log}({L}_{i+1}/{L}_{i})}{{{\rm{\Delta }}}_{{\rm{L}}}}\right]}^{2}\right\}}^{1/2}\approx 1\end{eqnarray} \tag{ 22 }$$

across the $({i}_{1},{i}_{2})$ pair, and similarly for subsequent pairs. Here, ${{\rm{\Delta }}}_{{\rm{T}}}$ and ${{\rm{\Delta }}}_{{\rm{L}}}$ are dimensionless weights which control the spacing of points in effective temperature and luminosity; we adopt the values 0.004 and 0.011, respectively, for these weights. At the selected points, GYRE searches for unstable oscillation modes with the harmonic degrees considered. First, GYRE solves the adiabatic oscillation equations to find eigenfrequencies ${\sigma }_{\mathrm{ad}}$ falling in the range extending from the asymptotic frequency of the gravity (g) mode with radial order n = 400, up to the asymptotic frequency of the pressure (p) mode with radial order n = 10. Each ${\sigma }_{\mathrm{ad}}$ is then used as an initial guess in finding a corresponding eigenfrequency $\sigma$ of the full non-adiabatic oscillation equations. The real and imaginary parts of $\sigma$ give the linear frequency ${\nu }_{\mathrm{osc}}$ and the growth e-folding time ${\tau }_{\mathrm{osc}}$ of a mode:

$$\begin{eqnarray}&&{\nu }_{\mathrm{osc}}=\displaystyle \frac{\mathrm{Re}(\sigma )}{2\pi },\qquad {\tau }_{\mathrm{osc}}=\displaystyle \frac{1}{\mathrm{Im}(\sigma )}.\end{eqnarray} \tag{ 23 }$$

If ${\tau }_{\mathrm{osc}}$ is negative, the mode is damped.

Separate strips are shown in Figure 9 for regions exhibiting unstable modes with $\mathrm{Re}(\omega )\gt 1$ and $\mathrm{Re}(\omega )\lt 1$, where

$$\begin{eqnarray}&&\omega =\sigma \sqrt{\displaystyle \frac{{R}^{3}}{{GM}}}\end{eqnarray} \tag{ 24 }$$

is the dimensionless eigenfrequency; these correspond, respectively, to the β Cephei and slowly pulsating B-type (SPB) classes of pulsating stars. In β Cephei stars during the MS phase, p and g modes with periods of a few hours and radial orders $n\approx 1-3$ are excited by a κ mechanism operating on the iron opacity bump situated in the outer envelope at $\mathrm{log}(T/{\rm{K}})\approx 5.3$ (Cox et al. 1992; Dziembowski & Pamiatnykh 1993). In SPB stars during the MS phase, g modes with periods of a few days and radial orders $n\approx 20-50$ are excited by the same mechanism (Dziembowski et al. 1993). For masses $M\gtrsim 9\;{{\rm{M}}}_{\odot }$ the strips for both classes of stars extend into the post-MS domain. During this phase, unstable modes couple with g modes trapped near the boundary of the inert helium core. In the case of the SPB stars this leads to very high overall radial orders, $n\gtrsim 100$, and ultimately limits our ability to follow the instability strips all the way to the red edge (our calculations are restricted to $n\lesssim 400$ for computational efficiency reasons). Hence, in Figure 9 we plot the red edges of the post-MS SPB strips with dotted lines, to highlight that these are not the true red edges.

Allowing for differences in adopted abundances and other modeling parameters, the instability strips plotted in Figure 9 are in general agreement with those published in the literature (e.g., Pamyatnykh 1999; Zdravkov & Pamyatnykh 2008; Saio 2011). The notable difference is the presence of fingers in the lower boundaries of our β Cephei strips for ${\ell }\geqslant 1$. Their appearance here is due to the unprecedented resolution in HR-diagram space of our stability calculations. To elucidate their origin, Figure 10 plots part of the ${\ell }=1$ frequency spectrum of an $8.5\;{{\rm{M}}}_{\odot }$ stellar model as it evolves from the ZAMS to the red edge of the main sequence (REMS), showing which modes are stable and which are unstable. The p₁ mode is unstable over the effective temperature range $4.358\geqslant \mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\geqslant 4.317$, and the g₁ mode over the cooler but overlapping range $4.341\geqslant \mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\geqslant 4.301$. The star then passes through a phase with no unstable modes, before the instability reappears in the range $4.288\geqslant \mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\geqslant 4.278$ for the g₂ mode.

Zoom In Zoom Out Reset image size

This alternation between instability and stability, seen as fingers in Figure 9, stems from the fact that the κ mechanism only excites modes whose eigenfrequencies fall in a narrow range $[{\sigma }_{\mathrm{lo}},{\sigma }_{\mathrm{hi}}]$. At frequencies $\mathrm{Re}(\sigma )\gt {\sigma }_{\mathrm{hi}}$, the pulsation period becomes comparable to the local thermal timescale in the envelope region above the iron opacity peak, and this region behaves as a damping zone, stabilizing the modes. Conversely, at frequencies $\mathrm{Re}(\sigma )\lt {\sigma }_{\mathrm{lo}}$, modes couple with gravity waves trapped in the μ-gradient zone developing at the core boundary, and are likewise damped. The intermediate stable phase in Figure 10, between $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=4.301$ and $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=4.288$ occurs when there are no modes in the $[{\sigma }_{\mathrm{lo}},{\sigma }_{\mathrm{hi}}]$ range. As the star evolves, the unstable range narrows: ${\sigma }_{\mathrm{hi}}$ decreases due to lower ${T}_{\mathrm{eff}}$, while ${\sigma }_{\mathrm{lo}}$ increases due to the growth of the μ-gradient zone.

Figure 11 shows a version of the ${\ell }=1$ panel calculated using OP opacity tables rather than OPAL tables. There is an overall shift of the instability strips toward higher luminosities, an effect already noted by Pamyatnykh (1999). The fingers persist with much the same structure, supporting the fact that they are physical effects rather than numerical artifacts.

Zoom In Zoom Out Reset image size

Returning now to Figure 9, the post-MS extension of the SPB strips has been attributed in the literature to features in the Brunt–Väisälä frequency which reflect gravity waves at the boundary of the helium core, preventing them from penetrating into the core and being dissipated by strong radiative damping. Saio et al. (2006) and Godart et al. (2009) argue that the necessary feature is an intermediate convection zone (ICZ) associated with the hydrogen-burning shell, but more recently Daszyńska-Daszkiewicz et al. (2013) have shown that even a local minimum in the Brunt–Väisälä frequency is sufficient to reflect modes. In the present case, the empirical mass threshold $M\gtrsim 9\;{{\rm{M}}}_{\odot }$ required for formation of an ICZ coincides with the lower boundaries of the SPB strip extensions. In the lowest-mass models above this threshold, the ICZ vanishes shortly after its appearance, but it leaves behind a narrow region with a steep molecular weight gradient. This gradient causes a spike in the Brunt–Väisälä frequency, which serves in a similar manner to prevent gravity waves from entering into the core and being dissipated.

The corresponding post-MS extension of the β Cephei strips was first noted by Dziembowski & Pamiatnykh (1993), but has not received much attention in the literature. Figure 9 shows that this extension has a well defined lower boundary, much like the SPB stars although situated at slightly higher masses, $M\gtrsim 10.5\;{{\rm{M}}}_{\odot }$. We have determined that the extension is also a consequence of ICZ formation; the shift to higher masses arises because it appears that multiple convection zones, rather than a single one, are necessary to reflect waves at the core boundary in the case of β Cephei pulsators.

To illustrate the updated asteroseismic capabilities of MESA, Figure 12 plots the echelle diagram for the subgiant star HD 49385, showing both the frequencies of ${\ell }=0-2$ modes measured by Deheuvels et al. (2010), and the corresponding frequencies of the best-fit model determined using the astero extension. The calculations follow the same procedure detailed in Section 3.2 of Paper II; the only significant differences are that the initial mass, helium abundance, metal abundance and mixing length parameter are refined using the downhill simplex algorithm rather than the Hooke–Jeeves algorithm; oscillation frequencies are calculated using GYRE rather than ADIPLS; and the surface corrections to frequencies are evaluated using Equation (4) of Ball & Gizon (2014) rather than with the Kjeldsen et al. (2008) scheme.

Zoom In Zoom Out Reset image size

Comparing Figure 12 against Figure 8 of Paper II reveals only small differences between the two. The ${\chi }^{2}$ of the best-fit models reported by astero is 2.3 in the former case, compared to 2.4 in the latter (cf. Table 2 of Paper II).

The instability strips presented above involved the examination of ∼11 million modes of ∼40,000 stellar models. To partially automate the process, we can leverage tight coupling between GYRE and MESAstar. This is achieved by making small modifications to the extras_check_model callback routine in MESAstar (see Appendix B.1 of Paper II), so that GYRE is run after each timestep to determine the set of eigenfrequencies $\{\sigma \}$ of a user-specified group of modes. When $\mathrm{Im}(\sigma )$ changes sign from one timestep to the next for any of these modes, indicating that an instability-strip boundary has been crossed, a search is performed to find $\mathrm{Im}$($\sigma$) $\approx$ 0.

Figure 13 presents an application of the tight coupling to the fundamental and first-overtone radial modes of the 8.5 ${{\rm{M}}}_{\odot }$ model considered in Section 3.1, showing how the growth timescales ${\tau }_{\mathrm{osc}}$ and oscillation periods ${P}_{\mathrm{osc}}$ = 1/${\nu }_{\mathrm{osc}}$ of the modes change as the star evolves from the ZAMS into the post-MS. The second-overtone radial mode remains stable, ${\tau }_{\mathrm{osc}}\lt 0$, over the range plotted. The vertical lines show where ${\tau }_{\mathrm{osc}}$ changes sign. The blue and red edges of the (${\ell }=0$) β Cephei instability strip can be seen in the upper panel of Figure 13 at $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=4.36$ and $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=4.30$, respectively. The blue edge¹² of the classical instability strip can likewise be seen in both panels at $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=3.75$. The corresponding red edge is not found because GYRE does not include a treatment of the pulsation–convection interaction—a necessary ingredient for modeling the classical red edge (see, e.g., Section 3.7.3 of Aerts et al. 2010, and references therein).

Zoom In Zoom Out Reset image size

As a further demonstration of automated instability strip location, Figure 14 plots the blue edges of the classical instability strip in the HR diagram, for fundamental and first-overtone radial modes. The edges are calculated for 51 evolutionary tracks spanning the initial mass range $1.25\;{{\rm{M}}}_{\odot }\leqslant M\leqslant 12.5\;{{\rm{M}}}_{\odot }$ in uniform logarithmic increments. At luminosities $\mathrm{log}(L/{{\rm{L}}}_{\odot })\gtrsim 2.5$ corresponding to classical Cepheid pulsators, these edges show good agreement with the set B results published by Smolec & Moskalik (2008, their Figure 1). At luminosities $\mathrm{log}(L/{{\rm{L}}}_{\odot })\lesssim 1.6$ corresponding to δ Scuti stars, the edges are somewhat cooler than results published in the literature; however, this is because we consider only fundamental and first-overtone modes, whereas the blue edge is typically set by higher overtones which are displaced toward hotter ${T}_{\mathrm{eff}}$ (e.g., Dupret et al. 2004, their Figure 1).

Zoom In Zoom Out Reset image size

Ideally, the same automated approach could be used to locate the boundaries of the non-radial (${\ell }\gt 0$) instability strips plotted in Figure 9. In practice it is very challenging to devise a robust algorithm that can unambiguously interpret the eigenfrequencies produced by GYRE. Sometimes, acoustic glitches in a model can trap modes in surface layers, where they are strongly excited; however, these modes are very sensitive to model parameters, and it is unclear whether they are physically meaningful or not.

The specific volume of cell k is

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\rho }_{k}}=\displaystyle \frac{4\pi }{3}\left(\displaystyle \frac{{r}_{k}^{3}-{r}_{k+1}^{3}}{{{dm}}_{k}}\right),\end{eqnarray} \tag{ 25 }$$

where r_k is the outer face radius, ${r}_{k+1}$ is the inner face radius, ${{dm}}_{k}$ is the cell mass, and ${\rho }_{k}$ is the cell average density (see Figure 15 for the layout of cells in MESAstar).We create an initial algebraic form of the continuity differential equation by dividing the change in the specific volume in a step by the length of time $\delta t$, using step start and end values for r_k, ${r}_{k+1}$, and ${\rho }_{k}$ and the value ${{dm}}_{k}$ which is constant during the step:

$$\begin{eqnarray}&&\displaystyle \frac{1/{\rho }_{k}-1/{\rho }_{\mathrm{start},k}}{\delta t}\\ &&=\displaystyle \frac{4\pi }{3{{dm}}_{k}}\displaystyle \frac{\left[\left({r}_{k}^{3}-{r}_{\mathrm{start},k}^{3}\right)-\left({r}_{k+1}^{3}-{r}_{\mathrm{start},k+1}^{3}\right)\right]}{\delta t}.\end{eqnarray} \tag{ 26 }$$

Next, we rewrite the right-hand side introducing variables for the time centered velocity $\hat{v}$ and the effective area A to get the final form of the mass continuity equation as used in ${\mathtt{MESAstar}}$:

$$\begin{eqnarray}&&\displaystyle \frac{1/{\rho }_{k}-1/{\rho }_{\mathrm{start},k}}{\delta t}=\displaystyle \frac{1}{{{dm}}_{k}}\left({A}_{k}{\hat{v}}_{k}-{A}_{k+1}{\hat{v}}_{k+1}\right),\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&{\hat{v}}_{k}=({v}_{k}+{v}_{\mathrm{start},k})/2\end{eqnarray} \tag{ 28 }$$

and r_k is evaluated as

$$\begin{eqnarray}&&{r}_{k}={r}_{\mathrm{start},k}+{\hat{v}}_{k}\delta t.\end{eqnarray} \tag{ 29 }$$

Algebraic simplification then shows that

$$\begin{eqnarray}&&{A}_{k}=\displaystyle \frac{4\pi }{3}\left({r}_{k}^{2}+{r}_{k}{r}_{\mathrm{start},k}+{r}_{\mathrm{start},k}^{2}\right).\end{eqnarray} \tag{ 30 }$$

To be consistent with the mass continuity equation, we use these expressions for effective area and time centered velocity in the following momentum and energy equations. It will be shown below that to get intrinsic energy conservation, we must time center the velocity and use special combinations of starting and ending radius in a couple of places, but all other terms in the equations can remain fully implicit to avoid degrading the numerical stability as would happen in a uniformly time-centered scheme.

Zoom In Zoom Out Reset image size

In ${\mathtt{MESAstar}}$, the artificial dynamic viscosity coefficient ${\eta }_{\mathrm{visc}}$ (which has the dimensions ${\rm{g}}\;{\mathrm{cm}}^{-1}\;{{\rm{s}}}^{-1}$)

$$\begin{eqnarray}&&{\eta }_{\mathrm{visc},k}={\eta }_{\mathrm{visc},\mathrm{linear},k}+{\eta }_{\mathrm{visc},\mathrm{quad},k},\end{eqnarray} \tag{ 31 }$$

where the linear term is

$$\begin{eqnarray}&&{\eta }_{\mathrm{visc},\mathrm{linear},k}=\displaystyle \frac{3}{4}{l}_{1}{\rho }_{k}{r}_{\mathrm{mid},k}{c}_{{\rm{s}},k}\end{eqnarray} \tag{ 32 }$$

and the quadratic term is

$$\begin{eqnarray}&&{\eta }_{\mathrm{visc},\mathrm{quad},k}=\displaystyle \frac{3}{4}\;{l}_{2}^{2}\;\displaystyle \frac{{\rho }_{k}^{2}4\pi {r}_{\mathrm{mid},k}^{2}}{{{dm}}_{k}}\ \mathrm{max}\left(0,{r}_{k+1}^{2}{\hat{v}}_{k+1}-{r}_{k}^{2}{\hat{v}}_{k}\right),\end{eqnarray} \tag{ 33 }$$

with ${r}_{\mathrm{mid},k}=({r}_{k+1}+{r}_{k})/2$, ${c}_{{\rm{s}},k}$ the sound speed in cell k, and l₁ (l₂) is a dimensionless coefficient for the linear (quadratic) term. The linear term is rarely used; it provides for a general damping of pulsations. The quadratic term is only nonzero in regions of compression and is the primary control for the strength of artificial viscosity. Assuming the usual case of l₁ = 0, the shock front is spread over a distance $\sim {l}_{2}{r}_{k}$. We follow Dorfi (1998) in opting for a shock spread proportional to the local radius r rather than the local cell width. This choice is dictated by the fact that step-by-step adjustments to the mesh resolution lead to dynamically changing, non-monotonic variations in cell widths of up to a factor of 2 or more between neighboring cells. Making the shock spread directly dependent on the local cell widths would produce numerically intolerable dynamically changing, non-monotonic variations in cell-to-cell values for the shock spread. Use of a local running average cell width is also ruled out by the need to keep algebraic equations dependent only on nearest neighbors to allow a block tridiagonal matrix solution. The use of a small fraction of the local radius gives a smoothly varying shock spread that avoids the numerical problems associated with using the cell width.

We define the quantity ${Q}_{\mathrm{visc},k}$ (having dimensions of energy), in cell k as

$$\begin{eqnarray}&&{Q}_{\mathrm{visc},k}={\eta }_{\mathrm{visc},k}\displaystyle \frac{{\rho }_{k}4\pi {r}_{\mathrm{mid},k}^{6}}{{{dm}}_{k}}\left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}}-\displaystyle \frac{{\hat{v}}_{k+1}}{{r}_{k+1}}\right).\end{eqnarray} \tag{ 34 }$$

The momentum equation uses ${Q}_{\mathrm{visc}}$ in an expression that defines an artificial acceleration analogous to the pressure gradient term, and the energy equation uses it to define an artificial viscous heating analogous to the mechanical work term.

The local linear momentum conservation equation at face k between inner cell k and outer cell $k-1$ is

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{k}-{v}_{\mathrm{start},k}}{\delta t}=-\displaystyle \frac{{{Gm}}_{k}}{{r}_{k}{r}_{\mathrm{start},k}}-{A}_{k}\left(\displaystyle \frac{{P}_{k-1}-{P}_{k}}{{\bar{{dm}}}_{k}}\right)+{g}_{\mathrm{visc},k},\end{eqnarray} \tag{ 35 }$$

where ${\bar{{dm}}}_{k}=({{dm}}_{k}+{{dm}}_{k-1})/2$ is the mass associated with face k, and the viscous acceleration term at face k is

$$\begin{eqnarray}&&{g}_{\mathrm{visc},k}=\displaystyle \frac{4\pi }{{r}_{k}}\left(\displaystyle \frac{{Q}_{\mathrm{visc},k-1}-{Q}_{\mathrm{visc},k}}{{\bar{{dm}}}_{k}}\right).\end{eqnarray} \tag{ 36 }$$

The use of the product ${r}_{k}{r}_{\mathrm{start},k}$ in the denominator of the gravitation term is necessary for intrinsic energy conservation as will be shown below.

The local energy conservation equation for cell k between outer face k and inner face $k+1$ is

$$\begin{eqnarray}\displaystyle \frac{{e}_{k}-{e}_{\mathrm{start},k}}{\delta t} & = & -\displaystyle \frac{{L}_{k}-{L}_{k+1}}{{{dm}}_{k}}-{P}_{k}\left(\displaystyle \frac{{A}_{k}{\hat{v}}_{k}-{A}_{k+1}{\hat{v}}_{k+1}}{{{dm}}_{k}}\right)\\ & & +{\epsilon }_{\mathrm{visc},k}+{\epsilon }_{\mathrm{nuc},k}-{\epsilon }_{\nu ,k}+{\epsilon }_{\mathrm{extra},k},\end{eqnarray} \tag{ 37 }$$

where e_k is the specific thermal energy for cell k. The viscous heating rate for cell k is

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{visc},k}=\displaystyle \frac{4\pi {Q}_{\mathrm{visc},k}}{{{dm}}_{k}}\left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}}-\displaystyle \frac{{\hat{v}}_{k+1}}{{r}_{k+1}}\right).\end{eqnarray} \tag{ 38 }$$

Energy loss from weak reaction neutrinos is already subtracted from the nuclear burning term, ${\epsilon }_{\mathrm{nuc},k}$, so only the neutrino energy loss rate from thermal processes, ${\epsilon }_{\nu ,k}$, is explicitly accounted for in Equation (37). An example of ${\epsilon }_{\mathrm{extra},k}$ is the artificial injection of energy to trigger a shock.

An alternative form of the energy equation equates the model ${dL}/{dm}$ to the expected value

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{k}-{L}_{k+1}}{{{dm}}_{k}}={\left(\displaystyle \frac{{dL}}{{dm}}\right)}_{\mathrm{expected},k},\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dL}}{{dm}}\right)}_{\mathrm{expected},k}={\epsilon }_{\mathrm{grav},k}+{\epsilon }_{\mathrm{visc},k}+{\epsilon }_{\mathrm{nuc},k}-{\epsilon }_{\nu ,k}+{\epsilon }_{\mathrm{extra},k},\end{eqnarray} \tag{ 40 }$$

and

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{grav},k}=-\displaystyle \frac{{e}_{k}-{e}_{\mathrm{start},k}}{\delta t}-{P}_{k}\left(\displaystyle \frac{{A}_{k}{\hat{v}}_{k}-{A}_{k+1}{\hat{v}}_{k+1}}{{{dm}}_{k}}\right).\end{eqnarray} \tag{ 41 }$$

Using Equation (27), the expression for ${\epsilon }_{\mathrm{grav}}$ can be rewritten

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{grav},k}=-\displaystyle \frac{{e}_{k}-{e}_{\mathrm{start},k}}{\delta t}-{P}_{k}\displaystyle \frac{(1/{\rho }_{k}-1/{\rho }_{\mathrm{start},k})}{\delta t}\end{eqnarray} \tag{ 42 }$$

thereby avoiding the use of velocities and thus be appropriate for hydrostatic cases.

The summed kinetic, potential, internal energies are

$$\begin{eqnarray}&&\mathrm{KE}=\displaystyle \sum _{k}\displaystyle \frac{1}{2}\ {\bar{{dm}}}_{k}\ {v}_{k}^{2},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\mathrm{PE}=\displaystyle \sum _{k}-\displaystyle \frac{{{Gm}}_{k}{\bar{{dm}}}_{k}}{{r}_{k}},\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&\mathrm{IE}=\displaystyle \sum _{k}{e}_{k}\ {{dm}}_{k},\end{eqnarray} \tag{ 45 }$$

and thus the total energy of the star is $E=\mathrm{KE}+\mathrm{PE}+\mathrm{IE}$. We now explicitly demonstrate that the equations we solve are formulated in such a way that the rate of change of total energy exactly equals the combined energy sources and sinks.

Multiplying Equation (35) by ${\hat{v}}_{k}{\bar{{dm}}}_{k}$ gives an equation with units of luminosity:

$$\begin{eqnarray}\displaystyle \frac{1}{2}\ {\bar{{dm}}}_{k}\ \left(\displaystyle \frac{{v}_{k}^{2}-{v}_{\mathrm{start},k}^{2}}{\delta t}\right) & = & -{{Gm}}_{k}{\bar{{dm}}}_{k}\left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}{r}_{\mathrm{start},k}}\right)\\ & & -{A}_{k}{\hat{v}}_{k}({P}_{k-1}-{P}_{k})\\ & & +4\pi \left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}}\right)({Q}_{\mathrm{visc},k-1}-{Q}_{\mathrm{visc},k}).\end{eqnarray} \tag{ 46 }$$

Using Equation (29) to eliminate ${\hat{v}}_{k}$ in the first term on the right,

$$\begin{eqnarray}&&-{{Gm}}_{k}{\bar{{dm}}}_{k}\left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}{r}_{\mathrm{start},k}}\right)=\displaystyle \frac{{{Gm}}_{k}{\bar{{dm}}}_{k}}{\delta t}\left(\displaystyle \frac{1}{{r}_{k}}-\displaystyle \frac{1}{{r}_{\mathrm{start},k}}\right),\end{eqnarray} \tag{ 47 }$$

shows that this term is the negative of the rate of change of potential energy, a result that is made possible by the use of the ${{Gm}}_{k}/{r}_{k}{r}_{\mathrm{start},k}$ in Equation (35) instead of an alternative such as ${{Gm}}_{k}/{r}_{k}^{2}$. Thus, Equation (46) can be written as

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\ {\bar{{dm}}}_{k}\ \displaystyle \frac{{v}_{k}^{2}-{v}_{\mathrm{start},k}^{2}}{\delta t}-\displaystyle \frac{{{Gm}}_{k}{\bar{{dm}}}_{k}}{\delta t}\left(\displaystyle \frac{1}{{r}_{k}}-\displaystyle \frac{1}{{r}_{\mathrm{start},k}}\right)\\ &&\quad ={-A}_{k}{\hat{v}}_{k}({P}_{k-1}-{P}_{k})+4\pi \left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}}\right)({Q}_{\mathrm{visc},k-1}-{Q}_{\mathrm{visc},k}).\end{eqnarray} \tag{ 48 }$$

Similarly, multiplying Equation (37) by ${{dm}}_{k}$ also yields an equation with units of luminosity:

$$\begin{eqnarray}\displaystyle \frac{({e}_{k}-{e}_{\mathrm{start},k})}{\delta t}{{dm}}_{k} & = & -({L}_{k}-{L}_{k+1})\\ & & -{P}_{k}\left({A}_{k}{\hat{v}}_{k}-{A}_{k+1}{\hat{v}}_{k+1}\right)\\ & & +4\pi {Q}_{\mathrm{visc},k}\left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}}-\displaystyle \frac{{\hat{v}}_{k+1}}{{r}_{k+1}}\right)\\ & & +({\epsilon }_{\mathrm{nuc},k}-{\epsilon }_{\nu ,k}+{\epsilon }_{\mathrm{extra},k})\ {{dm}}_{k}.\end{eqnarray} \tag{ 49 }$$

Adding Equations (48) and (49) and summing over the grid index k gives

$$\begin{eqnarray}&&{\displaystyle \displaystyle \sum }_{k}-\displaystyle \frac{{{Gm}}_{k}{\bar{{dm}}}_{k}}{\delta t}\left(\displaystyle \frac{1}{{r}_{k}}-\displaystyle \frac{1}{{r}_{\mathrm{start},k}}\right)\\ &&\quad +\displaystyle \frac{1}{2}\ {\bar{{dm}}}_{k}\ \displaystyle \frac{{v}_{k}^{2}-{v}_{\mathrm{start},k}^{2}}{\delta t}\\ &&\quad +\displaystyle \frac{({e}_{k}-{e}_{\mathrm{start},k})}{\delta t}{{dm}}_{k}\\ &&={\displaystyle \displaystyle \sum }_{k}-({L}_{k}-{L}_{k+1})\\ &&\quad +{P}_{k}{A}_{k+1}{\hat{v}}_{k+1}-{P}_{k-1}{A}_{k}{\hat{v}}_{k}\\ &&\quad +4\pi {Q}_{\mathrm{visc},k-1}\left(\displaystyle \frac{{\hat{v}}_{k}}{{r}_{k}}\right)-4\pi {Q}_{\mathrm{visc},k}\left(\displaystyle \frac{{\hat{v}}_{k+1}}{{r}_{k+1}}\right)\\ &&\quad +({\epsilon }_{\mathrm{nuc},k}-{\epsilon }_{\nu ,k}+{\epsilon }_{\mathrm{extra},k})\ {{dm}}_{k}.\end{eqnarray} \tag{ 50 }$$

The sum over the pressure terms is

$$\begin{eqnarray}&&{\displaystyle \displaystyle \sum }_{k}\left({P}_{k}{A}_{k+1}{\hat{v}}_{k+1}-{P}_{k-1}{A}_{k}{\hat{v}}_{k}\right)\\ &&\quad =-\left[{\left({PA}\hat{v}\right)}_{\mathrm{surface}}-{\left({PA}\hat{v}\right)}_{\mathrm{center}}\right]\\ &&\quad =-({L}_{\mathrm{acoustic},\mathrm{surface}}-{L}_{\mathrm{acoustic},\mathrm{center}}),\end{eqnarray} \tag{ 51 }$$

which cancels term by term except at the boundaries. We define ${L}_{\mathrm{acoustic},\mathrm{surface}}$ as the work done by the model on the atmosphere at the surface and ${L}_{\mathrm{acoustic},\mathrm{center}}$ as the work done on the model at the center, for example, by an artificial piston. The sum over the artificial viscosity terms leads to a similar expression, but because Q_visc vanishes at the surface and the center, the sum equals zero. That is, the energy added by artificial viscous heating in the energy equation exactly balances the loss of kinetic energy by artificial viscous acceleration in the momentum equation.

The terms on the left-hand side of Equation (50) are the difference in the total energy between the start and end of a step divided by the length of the step, in other words, the average rate of change of the total energy of the model. Therefore Equation (50) can be written as

$$\begin{eqnarray}({E}_{\mathrm{final}}-{E}_{\mathrm{initial}})/\delta t & = & -({L}_{\mathrm{surface}}-{L}_{\mathrm{center}})\\ & & -({L}_{\mathrm{acoustic},\mathrm{surface}}-{L}_{\mathrm{acoustic},\mathrm{center}})\\ & & +\displaystyle \sum _{k}\left({\epsilon }_{\mathrm{nuc},k}-{\epsilon }_{\nu ,k}+{\epsilon }_{\mathrm{extra},k}\right)\ {{dm}}_{k}.\end{eqnarray} \tag{ 52 }$$

This equation embodies conservation of energy in ${\mathtt{MESAstar}}$: the rate of change of total energy equals the combined energy sources and sinks. This demonstrates that in the given form, the algebraic equations intrinsically conserve energy in the sense that failure to get energy balance can only arise from the residual numerical errors that are inherent in using approximate solutions to the equations. This in turn means that to control energy balance errors, we can focus on reducing residuals either by changes in the Newton solver or by timestep reductions.

The Newton solver considers both the sizes of incremental changes to the variables and the sizes of residual errors for the equations. For the energy equation, the residual used by the solver is defined to be the timestep $\delta t$ times the difference between the left and right sides of Equation (37) divided by ${e}_{\mathrm{start},k}$; in other words, the residual is the error as a fraction of the specific energy at the start of the step. By adjusting tolerances for the average and maximum size of residuals, we force the Newton solver to take extra iterations to reduce the residuals which will in turn reduce the total error in energy conservation.

A second and related way to control energy conservation errors is to use the average and maximum energy residuals to adjust the timesteps. For example, if the maximum magnitude for an energy residual exceeds a specified hard limit, then the proposed solution is rejected and the step is retried with a smaller timestep. If the maximum is smaller than the hard limit but exceeds another specified limit, then there is no forced retry for this step, but the next timestep is reduced by the ratio of the limit divided by the maximum magnitude. If the maximum is smaller than both limits, then other factors determine the next timestep. Later in this section, we will show that these approaches in combination with the intrinsic conservation form of the equations yield a solution for a shock in an envelope that evolves with reasonably large timesteps while conserving energy to a high degree of accuracy.

When using hydrodynamics, we often require timesteps that are so small that we need to limit the increase in convective velocities as calculated in the standard instantaneous mixing length theory (MLT) so that they do not assume unphysically large accelerations. If convection velocities are allowed to adjust instantaneously, then our methods for artificially creating shocks will fail since however rapidly we inject energy over a limited region, convection will be able to transport the energy away. To be able to simulate shocks we need to have a way to limit convection velocity acceleration.

The primary scheme we use for this is derived from Arnett (1969) and Wood (1974). The MLT implementation in MESA has been extended to take as additional arguments the timestep and the previous convection velocity at the same mass location (${v}_{{\rm{c}},\mathrm{prev}}$). It calculates a provisional convection velocity (${v}_{{\rm{c}}0}$) using the standard instantaneous MLT, then defines a convective timescale (${\tau }_{\mathrm{MLT}}$) as the local pressure scale height (H) divided by the sum of the provisional plus previous velocities. If $\delta t$ is less than ${\tau }_{\mathrm{MLT}}$, then the next convection velocity (v_c) is only incremented from the previous one by the difference of the provisional minus the previous velocities times the ratio of the timestep divided by the timescale

$$\begin{eqnarray}&&{v}_{{\rm{c}}}={v}_{{\rm{c}},\mathrm{prev}}+\mathrm{min}\left(1,\displaystyle \frac{\delta t}{{\tau }_{\mathrm{MLT}}}\right)({v}_{{\rm{c}}0}-{v}_{{\rm{c}},\mathrm{prev}}),\end{eqnarray} \tag{ 53 }$$

where

$$\begin{eqnarray}&&{\tau }_{\mathrm{MLT}}=\displaystyle \frac{H}{({v}_{{\rm{c}}0}+{v}_{{\rm{c}},\mathrm{prev}})}.\end{eqnarray} \tag{ 54 }$$

As an alternative scheme for limiting convection acceleration, we also allow the maximum rate of change of convection velocity to be set as a fraction, ${g}_{\theta }$, of the local gravitational acceleration. If ${v}_{{\rm{c}}0}\gt {v}_{{\rm{c}},\mathrm{prev}}$, then

$$\begin{eqnarray}&&{v}_{{\rm{c}}}=\mathrm{min}({v}_{{\rm{c}}0},{v}_{{\rm{c}},\mathrm{prev}}+\delta t\ {g}_{\theta }\ g).\end{eqnarray} \tag{ 55 }$$

The final v_c is used to recalculate the convection efficiency, which is used to calculate the MLT temperature gradient.

These methods for limiting the acceleration of convective velocities reduce the energy transport rate as well as the rate of compositional mixing. Both schemes seem to give at least qualitatively reasonable results and avoid the problems of unphysically large accelerations that are possible with standard instantaneous MLT. Hopefully this ad hoc solution will soon be replaced by a quantitatively correct formulation.

MESA provides a variety of options for surface boundary conditions (see, e.g., Section 5.3 of Paper I), and several more have been added for use with hydrodynamics. The simplest allow specification of a particular value for the surface pressure, the surface temperature, or the ${T}_{\mathrm{eff}}$ if the surface is not at the photosphere. In the case of a given fixed surface pressure, the corresponding surface temperature is set using the surface luminosity and radius based on the usual blackbody relation. For the second case, where the surface temperature is fixed, the surface pressure is set to the corresponding radiation pressure. For both of these, if the surface is not at the photosphere, ${T}_{\mathrm{eff}}$ is set using the Eddington T–τ relation. Finally, for specified ${T}_{\mathrm{eff}}$ when the surface is not at the photosphere, the corresponding surface temperature is also derived using the Eddington T–τ relation, and the surface pressure is set to the radiation pressure for that temperature.

For computations involving shocks at the surface, there is an option to use boundary conditions that specify a vanishing gradient for compression at the surface and a temperature corresponding to blackbody radiation. The outermost cells ($k=1,2$) satisfy the equation

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\rho }_{1}}-\displaystyle \frac{1}{{\rho }_{\mathrm{start},1}}=\displaystyle \frac{1}{{\rho }_{2}}-\displaystyle \frac{1}{{\rho }_{\mathrm{start},2}},\end{eqnarray} \tag{ 56 }$$

which represents the vanishing of the surface compression gradient.

Finally, for computations involving interior shocks but low velocities at the surface, there is an option to use the surface pressure from the selected atmosphere prescription with the momentum equation relating the surface velocity to the surface pressure gradient. This form for the surface boundary conditions is used in the shocked massive star example in Section 6 and in the following envelope shock test.

To test the implementation, we shock the extended envelope of a 6.93 ${{\rm{M}}}_{\odot }$ asymptotic giant branch (AGB) star evolved from a 7 ${{\rm{M}}}_{\odot }$ MS star without rotation and an initial metallicity of 0.001. This case is chosen because of the uniform properties of the extended envelope (i.e., small density range, smooth density, and uniform composition). Our interest is to study the propagation of the shock, the properties of the shocked material, and the magnitude of energy conservation errors. In Section 6, we present results that mimic core-collapse supernovae.

Explosion simulations with MESA start from a converged model. The core is excised by removing inner shells of the model and setting new inner boundary conditions for mass, radius, velocity, and luminosity. For the current test, we remove the center just above the helium core at a mass of 2.40 ${{\rm{M}}}_{\odot }$ which corresponds to an inner radius of 27.2 ${{\rm{R}}}_{\odot }$. The stellar surface lies at a radius of 282.7 ${{\rm{R}}}_{\odot }$. During the following evolution, the excised region is treated as a point mass and is linked to the above layers by the inner boundary conditions which can be changed at each timestep to simulate various core events. The model grid was adjusted at each step to give higher resolution in the vicinity of the shock. The total number of cells stayed at about 1000, with cell masses dropping to about ${10}^{-4}$ of the total in approximately 100 cells around the shock.

In MESA, the artificial explosion that creates a shock can be produced in three ways: a piston, a luminosity flash, or a thermal bomb.

• 1.  
  The first option changes the inner boundary conditions for velocity and radius to mimic a piston. A core-collapse supernova can be simulated by moving the inner radius inwards (collapse) at a free-fall speed and then violently outwards (bounce and explosion). The parameterization for the piston-driven explosion is the same as that described in Woosley & Weaver (1995), and includes the infall piston time, the final inward piston radius, the initial outward piston speed, and the final piston radius.
• 2.  
  The second option increases the inner boundary luminosity over a specified time in order to deliver the desired total energy. In this approach, the inner boundary radius is fixed at all times and becomes a zero-flux inner boundary once the flash is over. In this approach, the inner boundary radius is fixed (zero velocity) and we inject the energy within the first zone of the domain.
• 3.  
  The third option deposits energy at a constant rate during a specified time and in a region bounded by two specified Lagrangian-mass coordinates. As in the second option, the inner boundary radius is fixed at all times.

The differences among these three options can alter the properties of the shocked envelope.

To benchmark MESA for these shock tests, we have used the explicit radiation-hydrodynamics code V1D (Livne 1993; Dessart et al. 2010a, 2010b). Options 1 and 3 are implemented in V1D. For the present envelope shock test, and subsequently for the explosion tests, we use option 3 in both codes. We initiate the explosion by depositing a total of 10⁴⁹ erg at a constant rate over 10 s between the Lagrangian mass coordinates of 2.40 and 2.45 ${{\rm{M}}}_{\odot }$. This energy deposited is well in excess of the initial binding energy, which is approximately $-2\;\times \;$10⁴⁷ erg. Once the energy injection is over, we save a model which is then used as initial conditions for a shock evolution simulation.

Once the stellar core has been excised, the remaining envelope has a smooth density profile, resembling a power law whose exponent is −1 at depth and decreases outwards to become about −10 at the photosphere (top row panels of Figure 16). Because convective accelerations are limited, the energy deposited increases the internal energy within the innermost 0.05 ${{\rm{M}}}_{\odot }$ of the grid. The pressure build-up leads to the sudden expansion of the innermost layers and the formation of a mildly supersonic shock (Mach number $\approx 2$). The shock propagates at a velocity in excess of 1000 $\mathrm{km}\;{{\rm{s}}}^{-1}$ initially, but slows to a few 100 $\mathrm{km}\;{{\rm{s}}}^{-1}$ by the time it reaches the stellar surface after 3 × 10⁵ s. The density contrast across this somewhat weak shock is $\approx 6$ . For a strong shock, one expects a density jump of 4 for an ideal gas with an adiabatic index of 5/3 and a value of 7 for a radiation-dominated gas ($\gamma =$ 4/3).

Zoom In Zoom Out Reset image size

This simulation is analogous to a shock-tube test. However, in the stellar context (realistic stellar envelope, realistic equation of state, spherical expansion), there is no analytical solution for comparison. We thus run the same simulation with the code V1D and include the results in Figure 16. The results agree at multiple times spanning the progression of the shock toward the stellar surface (the times used for comparison are the same to within 1% and the grid resolution is comparable). The sharpness of the shock in the two simulations differs with time and location. In V1D, the artificial viscosity has a physical spread of two grid zones, irrespective of radius, while in this MESA run, the spread is set to 0.1% of the local radius (see Section 4.2).

Since the explosion is started as a thermal bomb, the bulk of the energy is initially internal (see Figure 17). As the material expands and accelerates, the kinetic energy increases, mirroring the decrease in internal energy (essentially no energy is used to unbind the envelope). At the time of shock emergence, the internal and kinetic energies are comparable.

Zoom In Zoom Out Reset image size

In the present case, we can preserve good accuracy while still allowing timesteps an order of magnitude greater than the Courant time.¹³ The error in energy conservation at timestep i is

$$\begin{eqnarray}&&\delta {E}_{\mathrm{error},i}={E}_{i}-{E}_{i-1}-{E}_{\mathrm{sources},i},\end{eqnarray} \tag{ 57 }$$

where ${E}_{\mathrm{sources},i}$ is the right-hand side of Equation (52) multiplied by $\delta t$. The cumulative relative error in energy at a timestep n is

$$\begin{eqnarray}&&\mathrm{CRE}({t}_{n})=\displaystyle \frac{1}{{E}_{n}}\displaystyle \sum _{i=1}^{n}\delta {E}_{\mathrm{error},i}.\end{eqnarray} \tag{ 58 }$$

In the test case, after about 15,100 timesteps when the shock reaches 6.6 ${{\rm{M}}}_{\odot }$, the cumulative relative error has grown to about $-1.4\;\times \;{10}^{-6}$, corresponding to a roughly linear growth rate of about $-1\;\times \;{10}^{-10}$ per step (bottom panel in Figure 17). Note that at this stage of evolution the shock is nearing the outer edge of the envelope but has not actually broken out through the surface. Issues of shock break out are beyond the scope of the current implementation. Using the parameters selected for the test, the energy conservation with V1D is not as good as with MESA (the jumps in cumulative error correspond to times when the limit on the timestep are loosened); comparable accuracy can be obtained with V1D by reducing the explicit timestep well below the Courant limit.

Finally, to illustrate the effects of artificial viscosity, we vary the quadratic term l₂ that controls the spread of the shock in response to compression (see Equation (33)), with the explosion energy increased to 10⁵⁰ erg in order to produce a stronger shock, and otherwise the same parameters and initial conditions. The top panel of Figure 18 shows the artificial acceleration (${g}_{\mathrm{visc}}$) and energy (${\epsilon }_{\mathrm{visc}}$) terms that enter the momentum and the energy equations for ${l}_{2}\;=$ 0.001. The acceleration term is positive ahead of the shock, causing a pre-acceleration of the unshocked material, and negative behind the shock causing a deceleration of the post-shock material. The energy corresponding to those changes in momentum is balanced by the extra term for artificial viscous heating in the energy equation (${\epsilon }_{\mathrm{visc}}$). The lower panel of Figure 18 shows the expected increase in the width of the shock as we raise the parameter l₂. For the model with ${l}_{2}\;=$ 0.004, dots locate grid cells. Note that with the smallest value (${l}_{2}\;=$ 0.001), the velocity is showing small oscillations ("ringing") behind the shock indicating that we have reached a practical lower bound for the shock spread given the other parameter choices and the nature of the specific problem.

Zoom In Zoom Out Reset image size

We now test different summation methods for Equation (59) and demonstrate that improved accuracy of the summations reduces the number of timesteps with a commensurate reduction in the execution time—while producing the same answers to within the specified integration accuracy.

When the summations in Equation (59) are accumulated in IEEE 64-bit arithmetic (16 significant figures, real*8 precision in Fortran on most architectures; more specifically binary64 with round to nearest and round ties to even) the integration in Figure 19 takes 3062 timesteps using a variable-order Bader–Deuflhard integrator with a specified accuracy of ${\tau }_{\mathrm{int}}$ = 10⁻⁴ and a scaling value ${y}_{\mathrm{scale}}$ = 10⁻³. The specified accuracy ${\tau }_{\mathrm{int}}$ limits the maximum error over one timestep for any isotope. Other potential, but less demanding, choices for the meaning of ${\tau }_{\mathrm{int}}$ include limiting the average or root-mean-square error over one timestep for all isotopes. When an abundance is greater than ${y}_{\mathrm{scale}}$ a relative error is calculated, while for abundances smaller than ${y}_{\mathrm{scale}}$, the absolute error is calculated (e.g., Press et al. 1992). In essence, only abundances greater than ${y}_{\mathrm{scale}}$ can exert control on the size of the timestep.

When the summations are accumulated in IEEE 128-bit arithmetic (32 significant figures, real*16 precision in Fortran on most architectures), the same integration takes only 55 timesteps, a factor $\approx 50$ improvement in the number of timesteps, and a factor of $\approx 30$ less execution time. Both calculations returned the same answers to within the specified integration error tolerances. For tighter integration tolerances of ${\tau }_{\mathrm{int}}$ = 10⁻⁶ and ${y}_{\mathrm{scale}}$ = 10⁻⁵, the evolution with summations in IEEE 64-bit arithmetic takes 10,081 timesteps while the evolution with summations in IEEE 128-bit arithmetic takes 88 timesteps. This is a factor of $\approx 100$ improvement in the number of timesteps, a factor of $\approx 150$ in execution time, with both calculations again producing the same abundances to within the specified integration error tolerances. Both sets of integration tolerances are practical, everyday usage tolerances; they are not extreme cases of hypothetical interest only. Using low-order Rosenbrock and first-order Euler integrators also showed similar improvements in the number of timesteps when the summations were performed in IEEE 128-bit arithmetic instead of IEEE 64-bit arithmetic. We achieve a reduction in the number of timesteps and execution times regardless of the number of isotopes, choice of integrator, integration tolerances, or linear algebra solver. This improvement in efficiency is fundamentally driven by a reduction in the numerical noise of the function being integrated.

At temperatures larger than ≈5 × 10⁹ K, integrating Equation (59) can be challenging as terms in the summation usually become large and opposite in sign. As shown above, the classic symptom during an integration under these thermodynamic conditions is the integrator taking an excessive number of very small timesteps in order to satisfy the specified integration accuracy criteria. The traditional workaround to this numerical problem is abandoning a network integration at elevated temperatures and deploying equilibrium solution methods. This switching of methods raises its own numerical issues when used within the larger context of multi-dimensional simulations or stellar evolution models (see Section 5.2).

Unless precautions are taken the summation of large sets of numbers can be very inaccurate due to the accumulation of rounding errors. Methods for accurate summation within the bounds of a given arithmetic remain an active field of research (e.g., Demmel & Hida 2003; McNamee 2004; Ogita et al. 2005; Rump et al. 2008; Graillat & Ménissier-Morain 2012; Collange et al. 2014). These summation discrepancies also worsen on heterogeneous architectures—such as clusters with NVIDIA GPUs or Xeon Phi accelerators—which combine programming environments that may obey various floating-point models and offer different precision results.

The summations in Equation (59) for the neutron, proton, and α-particle abundances are especially prone to inaccuracies because every isotope in a network reacts with these three particles. We report on the summation errors for these three isotopes. Each term in the summations of Equation (59) is calculated using IEEE 64-bit arithmetic and then copied into a IEEE 128-bit variable using the Fortran promotion rules. Each IEEE 128-bit term is then imported into the MP (Brent 1978) and MPf90 (Bailey 1995) multiple precision packages. All other aspects of the integration were executed in IEEE 64-bit arithmetic. The summations are then accumulated with

• 1.  
  IEEE 64-bit: terms in the order as given;
• 2.  
  IEEE 64-bit: terms sorted by their absolute value in ascending order;
• 3.  
  IEEE 64-bit: terms sorted by their absolute value in ascending order and the Kahan (1965) algorithm, which reduces the numerical error in summation by retaining a separate variable to accumulate the errors;
• 4.  
  IEEE 128-bit: terms in the order as given;
• 5.  
  IEEE 128-bit: terms sorted by their absolute value in ascending order;
• 6.  
  MP and MPf90 100 digits: terms sorted by their absolute value in ascending order.

There are many summation methods and alternative multiple precision packages we did not deploy in these studies (e.g., Knuth 1997; Higham 2002; Li et al. 2002; Muller et al. 2010; Collange et al. 2014).

Table 4 summarizes these summation experiments. We confirm that 100 digits are sufficient to prevent errors in our multiple precision sums. Column 4 gives the minimum number of correct digits in a summation. There are three time periods in the evolution of Figure 19 where summations performed in IEEE 64-bit arithmetic greatly increase the number of timesteps taken by the integration. One is during the first rearrangement into the NSE state ending around 10⁻⁸ s, another is during the second rearrangement around 10⁻¹ s and the third time period is when the abundances do not change much ($\dot{Y}\;\approx$ 0) and reaction rates reach equilibrium. It is during these equilibrium periods where a summation in IEEE 64-bit arithmetic with the terms summed in the order given may yield only six accurate digits. As a result, 3062 timesteps are needed to complete the integration (e.g., row 2 of Table 4). It is important to note that this strategy and choice of arithmetic is commonly used by nuclear reaction networks (e.g., Timmes 1999)—and is the most inaccurate choice. Row 4 of Table 4 is an important case, sorted plus Kahan summation, because it demonstrates that a marginal improvement in the accuracy of the summation (8 minimum correct digits) has a major reduction on the number of timesteps (1174 timesteps) and execution time (a factor of $\approx 2.5$ smaller). This establishes the general trend that improved accuracy of the summations reduces the number of timesteps with a commensurate reduction in the execution time—while producing the same answers to within the specified integration accuracy.

Table 4.  Results of Summation Experiments

IEEE	Maximum	Strategy	Minimum	Number of	Ratio of
Arithmetic	Digits Compared		Correct Digitsa	Timestepsb	CPU Timesc
${\tau }_{\mathrm{int}}$ = 10−4 ${y}_{\mathrm{scale}}$ = 10−3
64-bit	16	in order given	6	3062	31.7
64-bit	16	sorted, ascending	7	2614	24.4
64-bit	16	sorted, Kahan sum	8	1141	13.1
128-bit	32	in order given	21	55	1.0
128-bit	32	sorted, ascending	22	55	1.0
${\tau }_{\mathrm{int}}$ = 10−6 ${y}_{\mathrm{scale}}$ = 10−5
64-bit	16	in order given	6	10081	156
64-bit	16	sorted, ascending	7	7972	123
64-bit	16	sorted, Kahan sum	8	7674	112
128-bit	32	in order given	21	88	1.0
128-bit	32	sorted, ascending	22	88	1.0

Notes.

^aRelative to the 100 digit sum by the MP and MPf90 multiple precision packages. ^bFor a Bader–Deuflard integrator in IEEE 64-bit arithmetic. ^cFor a single thread on one 2.7 GHz Intel Xeon E5 core with the Intel 15.0.1 Fortran compiler, and relative to the execution time for the integration with 128-bit summations with terms in the order given.

Download table as:  ASCIITypeset image

The left-hand side of Equation (59) for ${\dot{Y}}_{i}$ is a IEEE 64-bit array to be filled with one of the summations. Setting ${\dot{Y}}_{i}$ equal to one of the IEEE 128-bit summations (at least 22 digits of accuracy) or the multiple precision package summations gives the most efficient integration (55 timesteps, rows 5 and 6 in Table 4) because the Fortran precision demotion rules assure ${\dot{Y}}_{i}$ is accurate to the limit of IEEE 64-bit arithmetic. The next best strategy, but a distant second, is setting ${\dot{Y}}_{i}$ equal to the sorted, Kahan summation. The worst case is setting ${\dot{Y}}_{i}$ equal to the 64-bit arithmetic sum with the terms in the order they appear—which is a common approach (e.g., Timmes 1999).

Figure 20 shows the number of correct digits in 64-bit and 128-bit summations for $\dot{Y}$(⁴He) with the terms accumulated in the order they are given. The number of correct digits is measured against the 100 digit sum calculated by the multiple precision packages MP and MPf90. The choices for the integrator and integration tolerances are the same as in Figure 19. Figure 20 shows that the minimum number of accurate digits is usually within a few digits of the limit of IEEE 64-bit arithmetic, but degrades to 6 digits (see row 1 of Table 4) during a time period of intense isotope rearrangement. These relatively large inaccurate summations cause the right-hand side of Equation (59) to be poorly defined in IEEE 64-bit arithmetic. As a direct result, the integration of Equation (59) with IEEE 64-bit summations takes 3062 timesteps to complete. Sorting the terms in the sums in ascending order and using the Kahan summation algorithm results in an accurate digit pattern that is very similar except the number of accurate digits is improved by one or two (see row 3 of Table 4). As a direct result of the improved accuracy of the summations, the number of timesteps is reduced from 3062 to 1141.

Zoom In Zoom Out Reset image size

For the IEEE 128-bit sum relative to the 100 digit sum in Figure 20, the minimum number of accurate digits is usually near the limit of IEEE 128-bit arithmetic, but degrades to 21 digits (see row 5 of Table 4) during the second period of intense isotope rearrangement. Relative to the IEEE 64-bit summations the number of correct digits is improved by at least 15, consistent with our conversion of the IEEE 64-bit terms in the sum to IEEE 128-bit using the Fortran promotion rules. As a direct result of the improved accuracy of the summations the integration that used the IEEE 128-bit summation took only 55 timesteps to complete the evolution.

We found no notable improvements by increasing the accuracy of the summations for the Jacobian matrix used by the stiff ordinary differential equation integrators. For this problem, it is evidently more important to better define the function—the right-hand side of Equation (59)—than the Jacobian matrix holding the derivatives of the function.

Based on these experiments we currently choose to improve the accuracy of the summations by converting the terms in the order they appear from IEEE 64-bit to IEEE 128-bit using the Fortran promotion rules and adding the terms in IEEE 128-bit arithmetic. IEEE 128-bit precision is presently almost always implemented in software by a variety of techniques (e.g., double–double methods), since direct hardware support for IEEE 128-bit precision is presently rare. However, Table 4 shows the reduction in the number of timesteps from accumulating the sums in IEEE 128-bit arithmetic far exceeds the extra computational cost per addition.

At high temperatures, the traditional workaround for the numerical problem of inaccurate summations in IEEE 64-bit arithmetic is to forgo using a reaction network integration to evolve the abundances and nuclear energy generation rate and to replace it with equilibrium solution techniques. An example of such an equilibrium calculation is NSE, where a root-find for the neutron and proton chemical potentials is performed. Once these two chemical potentials are known, all the abundances can be determined from nuclear Saha equations (e.g., Clifford & Tayler 1965; Hartmann et al. 1985; Meyer et al. 1998; Nadyozhin & Yudin 2004; Seitenzahl et al. 2008; Odrzywolek 2012). Equilibrium solution methods by themselves are efficient, robust, and inexpensive.

However, combining reaction networks and equilibrium solution methods creates its own numerical issues, especially when the temperature and density are spatial and time dependent. For example, the temperature of a cell may start relatively low, move into quasi-static equilibrium (QSE) range above 3 × 10⁹ K, and then move into NSE range above 5 × 10⁹ K. Ad-hoc decision trees must be created for switching between a network integration, QSE solutions, and NSE solutions. These switches can introduce unphysical discontinuities in the abundances either from one timestep to the next or in the abundance spatial profiles from one cell to the next.

Furthermore, cells near the transition between a network integration and an equilibrium method can be unstable in the sense that the equilibrium solution can evolve a cell to lower temperatures pushing the cell into using a network integration, while the solution from the network integration can evolve the cell toward higher temperatures evolving the cell back toward using the equilibrium solution. Moreover, the reaction network used for the time integration is different (usually smaller) than the isotope listing used for the equilibrium methods. This necessitates crafting a delicate mapping between two abundance vectors, which may also introduce unphysical discontinuities. In addition, care must be taken to assure the reaction rate screening corrections used in the time integration are properly taken into account in the equilibrium solution method, otherwise a fundamental incompatibility exists between the abundance vectors.

Finally, equilibrium methods determine the composition at a fixed electron fraction ${Y}_{{\rm{e}}}$. It then becomes necessary to solve an ordinary differential equation for $\dot{{Y}_{{\rm{e}}}}$ based on weak reaction rates in order to advance the abundance solution with a time varying ${Y}_{{\rm{e}}}$ (McLaughlin et al. 1996; Townsley et al. 2009; Arcones et al. 2010, also see Section 8). Switching between integration and equilibrium methods mid-stream is a liability, not a positive asset.

The need for traditional workarounds forced by limited accuracy of the summation is now avoided. The summation experiments in Section 5.1 demonstrate that network integration can be robust and efficient, even at very high temperatures, when the accuracy of the summations is improved. We stress this is not just a solution to issues of limited accuracy. It also offers an improvement in MESA by providing a single solution methodology, network integration, that avoids the challenges of stitching together different solution methods.

The new capabilities described above allow MESAstar to use large in situ reaction networks (i.e., fully coupled to the stellar evolution rather than uncoupled co-processing). A demonstration is Type 1 X-ray bursts, a class of objects with unstable nuclear burning on the surface of a NS. These bursts are sensitive functions of accretion rate (Chen et al. 1997), accretion composition (Galloway et al. 2006), the spatial distribution of burning on the surface of the NS (Bildsten 1995), the type of burning that occurs between bursts (Galloway et al. 2008) as well as possibly other conditions, for instance "superbursts" where carbon, rather than H/He, burns (Cumming & Bildsten 2001). Here we focus on a simplified model of constant accretion rate, where the burning occurs over the whole surface of the NS. GS 1826-24 (Tanaka 1989), also known as the "clocked burster" (Ubertini et al. 1999), provides an example of such a system due to its regular Type 1 X-ray bursts.

As material is accreted at the surface of a NS it is compressed and heats the underlying material. The accreted hydrogen (from a low mass MS star (Chen et al. 1997)) burns via the hot CNO cycle. However, with high enough accretion rates the hydrogen will be accreted faster than the hot CNO cycle, which is limited by the β-decay timescale (of order minutes), can process the material. The accreted helium ignites unstably in a hydrogen rich environment, allowing rapid proton (rp) captures onto seed nuclei (Wallace & Woosley 1981). This process forms nuclei along the proton drip line up to and beyond the iron group (Schatz et al. 1999), peaking at ¹⁰⁷Te, when α-decays prevent heavier elements from being formed (Schatz et al. 2001; Fisker et al. 2008). Once the burst begins, convection will commence, mixing the freshly burnt material with the ashes of previous burning episodes (Weinberg et al. 2006).

GS 1826-24 has been studied by the Rossi X-Ray Timing Explorer (RXTE) over several years (Galloway et al. 2004, 2008). The bursts showed a decrease in the recurrence time between bursts, from 4.1 $\mathrm{hr}$ in 2000 to 3.56 $\mathrm{hr}$ in 2002, though during each observational epoch the bursts were consistent with each other. Based on the ratio of the burst energy to the persistent flux, it is assumed that the bursts are powered by hydrogen burning of solar metallicity material.

We model the NS envelope using inner boundary conditions for mass and radius of ${M}_{{\rm{c}}}=1.4\;$ ${{\rm{M}}}_{\odot }$ and ${R}_{{\rm{c}}}=11.2\;$ ${\rm{k}}{\rm{m}}$ (Heger et al. 2007), implying a gravitational redshift of $1+z=1.26$. The base of the envelope is composed of an inert layer that does not undergo reactions. The luminosity at the base of the envelope is set to $L=1.6\;\times \;{10}^{34}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$ (Woosley et al. 2004). We base our nuclear networks on the 304 species rp.net network of Fisker et al. (2008), which includes proton rich isotopes up to ¹⁰⁷Te. Isotopes above ${}^{66}\mathrm{Zn}$, which is the peak isotope in the mesa_204.net, are included due to the proton captures possible on high-Z isotopes during the peak of the burst (Fisker et al. 2006). We also include the effects of rotational mixing by setting a minimum amount of mixing in the NS envelope. This mixing, while having a physical motivation (Piro & Bildsten 2007; Keek et al. 2009), is there primarily to improve the convergence of MESA models by smoothing out the compositional gradients that form in the ashes of previous bursts. We include the post-Newtonian correction to correct the local gravity in each cell for GR effects. During the burst we allow the accretion to continue.

Our results are compared to the RXTE observations of GS 1826-24 over bursts 9-20. Time resolved spectra were binned during the bursts' rise time and decay (Galloway et al. 2008; Zamfir et al. 2012). Data output by MESA is not GR time corrected, thus we set the burst times to be $t^{\prime} =t(1+z)$, and average multiple bursts to produce a scaled light-curve.

Figure 21 shows the temperature profile during two X-ray bursts, for the rp_305 net, accreting solar metallicity material at a rate of $3.0\;\times \;{10}^{-9}\;{{\rm{M}}}_{\odot }\;{\mathrm{yr}}^{-1}$. At $t^{\prime} \approx -10\;{\rm{s}}$ the envelope ignites material and drives the formation of the first convection zone. This zone expands outwards in the envelope mixing the ashes from the burning at the base of the envelope outwards to lower pressures (Weinberg et al. 2006). As the burst decays the convection zone recedes outwards and by $t^{\prime} \approx 150\;{\rm{s}}$ the envelope returns to its pre-burst temperature profile.

Zoom In Zoom Out Reset image size

We test three reaction networks, rp_53, rp_153 and rp_305, each a modified form of that in Fisker et al. (2008). Table 5 shows that increasing the number of isotopes in the reaction network increases the recurrence time and that all (for $\dot{M}=3.0\;\times \;{10}^{-9}\;{{\rm{M}}}_{\odot }\;{\mathrm{yr}}^{-1}$) have recurrence times $\approx 1\;\mathrm{hr}$ less than that of GS 1826-24.

Table 5.  Recurrence Times of X-Ray Bursts

Model	Accretion Rate (${10}^{-9}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}$)	Composition	Recurence Time (hr)
GS 1826-24			4.0750 ± 0.0003
rp_53	3.00	2% metals	1.5 ± 0.10
rp_153	3.00	2% metals	3.3 ± 1.80
rp_305	3.00	2% metals	3.2 ± 0.07
rp_305	3.00	2% ${}^{14}{\rm{N}}$	3.0 ± 0.07
rp_305	2.40	2% metals	4.1 ± 0.30
Heger et al. (2007)	1.17	2% ${}^{14}{\rm{N}}$	5.4 ± 0.10

Download table as:  ASCIITypeset image

Figure 22 and its insert show the folded light curves for each of the three rp reaction networks plus the GS 1826-24 observations. The rise time is sensitive to the net, with the largest net matching the observed slow rise. The observed decay profile is also best matched by rp_305. Burst to burst variations of the models decrease with increasing net size and can be further reduced by increasing the temporal resolution of the models. However, increasing the size of the net reduces the variation without having to increase the temporal resolution and also highlights the impact of MESA's capability to include large nuclear networks.

Zoom In Zoom Out Reset image size

To achieve a better match to the GS 1826-24 recurrence time (see Table 5), we reduce the accretion rate to $\dot{M}=2.4\;\times \;{10}^{-9}\;{{\rm{M}}}_{\odot }\;{\mathrm{yr}}^{-1}$. However, Figure 23 shows that the light curve comparisons are not as good as for the higher $\dot{M}$.

Zoom In Zoom Out Reset image size

GS 1826-24 was also modeled by Heger et al. (2007) with accretion of hydrogen, helium and 2% ${}^{14}{\rm{N}}$. For comparison, we run a model with this same composition with the rp_305 net. Table 5 shows that the recurrence time decreases slightly when accreting 2% ${}^{14}{\rm{N}}$ rather than 2% metals. The model with metal accretion is in better agreement with both the light curve rise and decay.

We now explore adaptive nets (Woosley et al. 2004), where we allow MESA to determine which isotopes (and reactions) are necessary by assessing the available reaction pathways for the most abundant isotopes. The network is constructed by first finding those isotopes with an abundance above a threshold, ${X}_{\mathrm{keep}}$, and then introducing those isotopes which are connected by adding or removing protons, neutrons, or α particles. That determination is made via the additional parameters ${X}_{{\rm{n}}}$ (i.e., neutron reactions) and ${X}_{{\rm{p}}}$ (i.e., proton and α reactions) potentially re-adding isotopes removed with the initiating ${X}_{\mathrm{keep}}$ threshold.

Accreting solar composition material at $\dot{M}=3.0\;\times {10}^{-9}\;{{\rm{M}}}_{\odot }\;{\mathrm{yr}}^{-1}$ we follow the model to the second burst, finding a recurrence time of $3.1\;\mathrm{hr}$, comparable to that from the rp_305 net (Table 5). The adaptive net has a better rise time profile than the rp nets, while the rp_305 net has a better fit to the decay. This gives us confidence that the rp_305 net includes all relevant isotopes which drive the X-ray burst and thus is a useful approximation. For suitable values for the sensitivity of the adaptive net, the net limits itself to $\approx 400$ isotopes between bursts, which increases to $\approx 600$ isotopes during the burst. Variations of a factor 100 in the threshold parameters only change the isotope count by at most 50 isotopes and do not affect the final results.

Figure 24 shows the ${T}_{{\rm{c}}}-{\rho }_{{\rm{c}}}$ evolution of M_i = 15 and 30 ${{\rm{M}}}_{\odot }$ models from the onset of carbon burning until iron-core collapse. These non-rotating, solar metallicity models used the 204 isotope reaction network described in Section 5 and MESAstar's "Dutch" mass loss prescription with η = 0.8. These models have ≈2200 cells on the MS, ≈3500 cells as the star becomes a red supergiant, and ≈2300 cells at the onset of core collapse. At core collapse the final masses are M_f = 13.0 and 15.2 ${{\rm{M}}}_{\odot }$. The curves fall below the ${T}_{{\rm{c}}}^{3}\propto {\rho }_{{\rm{c}}}$ scaling relation because the core becomes partially electron degenerate, as indicated by tracks crossing the Fermi energy ${E}_{{\rm{F}}}/({k}_{{\rm{B}}}T)\;\approx$ 4 curve. Evolution toward lower density at nearly constant temperature signals ignition of a nuclear fuel.

Zoom In Zoom Out Reset image size

Figure 25 shows the radial velocity and ${Y}_{{\rm{e}}}$ profiles at the onset of core collapse for the M_i = 15 ${{\rm{M}}}_{\odot }$ model. Dashed curves show the results using a 22 isotope network and solid curves show the results using a 204 isotope network. Both models are evolved from the pre MS to the onset of core collapse with their respective reaction network. The vertical gray lines mark the mass of the iron core as defined by the ${Y}_{{\rm{e}}}$ jump, which is $m\;\approx$ 1.43 ${{\rm{M}}}_{\odot }$ for the 204 isotope model and $m\;\approx \;$ 1.59 ${{\rm{M}}}_{\odot }$ for the 22 isotope model. The infall speed has reached ≈1000 km s⁻¹ just inside these iron core locations.

Zoom In Zoom Out Reset image size

Figure 26 shows the thermodynamic profiles at the onset of core collapse for the M_i = 15 ${{\rm{M}}}_{\odot }$ model. Dashed curves again show the results using a 22 isotope network and solid curves show the results using a 204 isotope network. The vertical gray lines again mark the mass of the iron core as defined by the ${Y}_{{\rm{e}}}$ jump in Figure 25. The impact of these differences remains to be explored.

Zoom In Zoom Out Reset image size

Figure 27 shows the mass fraction profiles of the ten most abundant isotopes within the iron core at the onset of core collapse for the M_i = 15 ${{\rm{M}}}_{\odot }$ model evolved with the 204 isotope network. Each isotope shown dominates the NSE composition at some location within the iron core, although we stress that no NSE or QSE approximation was used; the same 204 isotope reaction network was used throughout the entire model from the pre-MS to the onset of core collapse.

Zoom In Zoom Out Reset image size

The most abundant isotopes in an NSE distribution generally have an individual ${Y}_{{\rm{e}}}$ that is within a small range of the local ${Y}_{{\rm{e}}}$. A small spread usually exists due to nuclear structure effects. For example, the dominant isotopes at the center in Figure 27 are ⁴⁹Sc and ⁴⁸Ca. These isotopes have individual ${Y}_{{\rm{e}}}$ of 0.429 and 0.417, respectively; commensurate with the central ${Y}_{{\rm{e}}}$ $\approx$ 0.428 shown in Figure 26. The dominant isotope changes as the NSE distribution adapts to the rapidly decreasing density profile and increasing ${Y}_{{\rm{e}}}$ profile. All the isotopes in the iron core eventually become part of the compact remnant after the explosion. However, the thermodynamic and composition profiles near the mass cut depend on the profiles interior to the mass cut.

The envelope shock tests described in Section 4.9 show that the hydrodynamic solver in MESA meets the basic requirements for shock propagation in a star. The AGB star model was selected because of the well behaved conditions of its envelope—a density structure that is smooth and monotonically declining toward the stellar surface, and a uniform composition.

Here we explore the more challenging conditions associated with a strong shock born deep in the stellar interior of a massive star. We study the dynamics of such a supernova shock and the explosive nucleosynthesis that takes place in the wake of the shock during the first second. The yields from explosive nucleosynthesis depend on both the energy and the power (characteristic energy deposition timescale), while the dynamics of the shock are primarily dependent on the total energy deposited.

The starting conditions for the explosion simulations are the two 15 ${{\rm{M}}}_{\odot }$ pre-supernova models discussed above; one for the approximate 22 isotope network and one for the 204 isotope network.

6.2.1. Explosion Dynamics

Since the focus here is on the dynamics rather than nucleosynthesis, we expedite the MESA simulation by using the 22 isotope network for both the pre-supernova and the explosion phases. Before triggering the explosion, we remove the iron core of the red-supergiant star by placing the inner boundary of the grid at a Lagrangian mass of 1.75 ${{\rm{M}}}_{\odot }$. We trigger the explosion by depositing 1.52 × 10⁵¹ erg at a constant rate during 1 s. The artificial viscosity is raised during the energy deposition phase (i.e., ${l}_{2}$ = 0.01), when the shock is at small radii, and lowered in the subsequent evolution until shock breakout (i.e., ${l}_{2}$ = 0.003). Since the binding energy of the envelope to be shocked is −3.2 × 10⁵⁰ erg at the time we trigger the explosion, this choice of energy deposition yields a total energy at the end of the deposition phase of 1.2 × 10⁵¹ erg. This is generally considered a standard value for a core-collapse supernova.

Figure 28 shows that the development of the explosion is analogous to the tests using the (low-density) envelope of an AGB star in Section 4.9, but with significant quantitative differences. Here, the shock born at the edge of the iron core first travels through the dense CO-rich core. At the outer edge of the He-rich shell, the shock traverses a steep density gradient to enter the low-density H-rich envelope. Hence, the shock crosses regions with densities ranging from ∼10⁶ g cm⁻³ at the edge of the iron core down to 10⁻¹⁰ g cm⁻³ at the stellar surface.

Zoom In Zoom Out Reset image size

The radii of the innermost shells are initially very small since they lie at the outer edge of the iron core. Consequently, they suffer considerable cooling from expansion. Figure 28 shows a drop in temperature from a few 10⁹ K down to ∼10⁴ K at ∼1 day. In addition, the supernova shock splits into a reverse/forward shock structure when it encounters the density drop at the transition between the He-rich core and the H-rich envelope. The reverse shock is the new feature, absent in the envelope shock test, that causes a significant deceleration of He-core material. The conversion of kinetic energy into internal energy causes this inner material to heat up, erasing the cooling effect from expansion. The innermost layers, which travel the slowest, will be shocked last. These innermost zones can evolve to temperatures ∼10⁴ K. It is in these innermost regions at late times that the differences between MESA and V1D are the largest. The offset occurs in a region of relatively high density (∼10⁻⁷ ${\rm{g}}\;{\mathrm{cm}}^{-3}$) and low temperature (∼10⁴ K). The offset in temperature between the MESA and V1D simulations at late times stems from a difference in the equation of state for metal-rich regions. MESA accounts for ionization through the OPAL equation of state table for metallicities z < 0.04. For higher metal abundances where OPAL tables are unavailable, MESA currently assumes full ionization while V1D solves for the ionization state of the gas. Note that such density/temperature regimes are normally not encountered in stellar interiors. For other quantities and/or locations/times, the agreement between MESA and V1D is excellent.

We also note that in the MESA simulation, two small spikes appear in the temperature and density profiles at <2.5 ${{\rm{M}}}_{\odot }$ at ≈1000 s after the energy deposition phase. This feature is absent in V1D because V1D uses a much larger viscous spread when the shock is in the helium-rich core ($R\lt {{\rm{R}}}_{\odot }$). One can reduce or eliminate such spikes by increasing the viscous spread, although this may visibly smear the shock when it crosses the H-rich envelope—the current choice seems a suitable compromise.

In contrast to the envelope shock test, this supernova explosion configuration raises the temperature by a factor of about ten. Consequently, because ${P}_{\mathrm{rad}}/{P}_{\mathrm{gas}}\propto {T}^{3}/\rho$, the post-shock material becomes completely radiation dominated (${P}_{\mathrm{rad}}\gg {P}_{\mathrm{gas}}$). If we neglect the binding energy and the kinetic energy of the post-shock material, the post-shock energy is of the order of the explosion energy. We indeed find a good correspondence between the post-shock temperature computed by MESA and the temperature obtained from ${({E}_{0}/{aV})}^{1/4}$ (where a is the radiation constant, E₀ is a fitting parameter, typically of the order of the explosion energy, and $V=4\pi {R}_{\mathrm{sh}}^{3}/3$ is the volume within the shock radius ${R}_{\mathrm{sh}}$). As expected, we also find that the shock accelerates (decelerates) in regions where ${\rho }_{\mathrm{sh}}{R}_{\mathrm{sh}}^{3}$ decreases (increases) outward.

6.2.2. Explosive Nucleosynthesis

Here we compare the shock nucleosynthesis results from the two independent codes, MESA and V1D. The same initial 204 isotope pre-supernova model was the starting point. Our first test case is a strong explosion triggered by injecting 1.57 × 10⁵¹ erg for 0.05 s and within 0.02 ${{\rm{M}}}_{\odot }$ of the mass cut, which is positioned at the outer edge of the iron core at 1.5 ${{\rm{M}}}_{\odot }$. The exact choice of explosion energy, deposition timescale, and mass cut is not strictly relevant.

Figures 29 and 30 compare the mass fraction profiles of MESA with a 22 isotope network, MESA with a 204 isotope network, and V1D with a 54 isotope network. The first comparison at 0.0 s shows the impact of mapping from the pre-supernova 204 isotope network to the networks used in the shock nucleosynthesis test. The next comparison at 0.05 s is at the end of the energy deposition phase. The final comparison at 42.7 s is after explosive nucleosynthesis has completed. In all cases, the silicon-rich and oxygen-rich shells are strongly influenced by the explosion; the former primarily for the production of ${}^{56}\mathrm{Ni}$ and the latter primarily for the production of ${}^{28}\mathrm{Si}$ and ${}^{32}{\rm{S}}$. The ⁵⁶Ni yields at 42.7 s are 0.092 ${{\rm{M}}}_{\odot }$ for V1D, 0.087 ${{\rm{M}}}_{\odot }$ for MESA with 22 isotopes, and 0.096 ${{\rm{M}}}_{\odot }$ for MESA with 204 isotopes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall the agreement between MESA and V1D on this strong explosion is very good. The small differences between MESA and V1D in Figures 29 and 30 are due to the difference in mapping procedures.

MESA follows rules for mapping isotopes from one network to another network: if an isotope present in the old network is also present in the new network, then the abundance from the old network is copied to the abundance for the new network. Isotopes in the new network that are not present in the old network are initially given a mass fraction of zero. MESA then separately renormalizes classes of isotopes to have the same total mass fraction in the new network as in the old network. The classes are neutrons, hydrogen, helium, carbon, nitrogen, oxygen, and other metals. This procedure guarantees that the sum of the mass fractions in a given class will be the same in the new network as in the old network. V1D's mapping procedure for isotopes is the following. For any isotope present in the V1D network but absent in the MESA input the mass fraction is set to the solar metallicity value. When an isotope included in the MESA input is absent in the V1D network, this isotope is left out in the V1D simulation. After completing the mappings, the resulting composition is renormalized so that the sum of the mass fractions is unity.

Our second test case is a lower power explosion. We inject 1.326 × 10⁵¹ erg in 1.0 s over 0.05 ${{\rm{M}}}_{\odot }$ of the mass cut, which is also positioned at the outer edge of the iron core at 1.5 ${{\rm{M}}}_{\odot }$. The total energy after the deposition phase is 10⁵¹ erg.

In Figure 31, we show the composition profiles for the eight most abundant isotopes in the inner ≈0.4 ${{\rm{M}}}_{\odot }$ at the end of the energy-deposition phase (i.e., at 1 s). The correspondence between MESA and V1D is again very good. The ${}^{56}\mathrm{Ni}$ mass fraction approaches unity—it would reach unity if we appreciably increased the power (see Figure 29 for example). Some ${}^{58}\mathrm{Ni}$ is produced in the same region, while ${}^{54}\mathrm{Fe}$ is synthesized in the layers immediately above. The ⁵⁶Ni yields at 1.0 s are 0.0041 ${{\rm{M}}}_{\odot }$ for V1D, and 0.011 ${{\rm{M}}}_{\odot }$ for MESA with 204 isotopes.

Zoom In Zoom Out Reset image size

This work shows that the power of the explosion has a significant impact on the abundance profiles. In the high power explosion, the yield of ${}^{56}\mathrm{Ni}$ is ≈10 times larger and the ${}^{4}\mathrm{He}$ is several orders of magnitude more abundant. The nucleosynthesis of the low power explosion is completed at end of deposition phase at 1.0 s, while nucleosynthesis in the high power case continues for ≈30 s. This sensitivity suggests potentially observable signatures between low and high power explosions. In addition, the explosive nucleosynthesis that takes place in core-collapse supernovae is sensitive to the way the explosion is triggered. With the approach we use (fixed power during the energy deposition phase), we find that increasing the explosion energy (at a given power), the power (at a given explosion energy), or both alters the amount of mass burnt. Moving the mass cut deeper into denser layers considerably enhances the amount of burnt material but this material may fall back rather than be ejected. Moving the mass cut further out into lower-density regions may completely quench the production of ${}^{56}\mathrm{Ni}$, in favor, for example, of ${}^{28}\mathrm{Si}$. It is thus important to keep in mind that the piston or thermal explosion trigger is artificial and that the yields from explosive nucleosynthesis bear significant uncertainties.

The independent coordinates used for writing the time-dependent structure of the star are m and t, and for fluid elements deep within the star at both timesteps, the conventional form of Equation (60) is adequate. One numeric subtlety of accretion is that the derivative at constant m cannot be evaluated for material that is not present in the star at the beginning of the timestep. However, the simplification available when ${\tau }_{\mathrm{th}}\ll {\tau }_{\mathrm{acc}}$, manifest in Equation (61), enables ${\epsilon }_{\mathrm{grav}}$ to be evaluated in the outer regions.

When T and ρ are used as independent variables, we write

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{grav}}={-C}_{P}T\left[(1-{{\rm{\nabla }}}_{\mathrm{ad}}{\chi }_{T}){\left(\displaystyle \frac{\partial \mathrm{ln}T}{\partial t}\right)}_{m}-{{\rm{\nabla }}}_{\mathrm{ad}}{\chi }_{\rho }{\left(\displaystyle \frac{\partial \mathrm{ln}\rho }{\partial t}\right)}_{m}\right]\end{eqnarray} \tag{ 62 }$$

in the interior of the star, as in Paper I, and near the surface we choose to write, using Equation (61),

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{grav}}={\epsilon }_{\mathrm{grav},\mathrm{nh}}+{\epsilon }_{\mathrm{grav},{\rm{h}}},\end{eqnarray} \tag{ 63 }$$

where

$$\begin{eqnarray}{\epsilon }_{\mathrm{grav},\mathrm{nh}} & = & {-C}_{P}T\left[(1-{{\rm{\nabla }}}_{\mathrm{ad}}{\chi }_{T}){\left(\displaystyle \frac{\partial \mathrm{ln}T}{\partial t}\right)}_{q}\right.\\ & & -\left.{{\rm{\nabla }}}_{\mathrm{ad}}{\chi }_{\rho }{\left(\displaystyle \frac{\partial \mathrm{ln}\rho }{\partial t}\right)}_{q}\right]\end{eqnarray} \tag{ 64 }$$

and

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{grav},{\rm{h}}}=\displaystyle \frac{{C}_{P}{TGm}\dot{M}}{4\pi {r}^{4}P}\left({{\rm{\nabla }}}_{\mathrm{ad}}-{{\rm{\nabla }}}_{T}\right).\end{eqnarray} \tag{ 65 }$$

Here ${{\rm{\nabla }}}_{T}=d\mathrm{ln}T/d\mathrm{ln}P$ is the T–P profile in the star, and we have used the thermodynamic derivatives ${{\rm{\nabla }}}_{\mathrm{ad}}={(\partial \mathrm{ln}T/\partial \mathrm{ln}P)}_{s}$, ${\chi }_{T}={(\partial P/\partial T)}_{\rho }$, and ${\chi }_{\rho }={(\partial P/\partial \rho )}_{T}$. There is also a transition region where a weighted combination of these forms is used, with weights varying linearly in m.

Placement of the transition is related to the mesh. The MESAstar mesh structure is unchanged from that discussed in Paper I and Paper II. An illustration of the mesh regions and an indication of the behavior of cell boundaries during accretion is shown in Figure 32. In the case of mass loss, analogous operations are performed; here, we focus only on mass gain.

Zoom In Zoom Out Reset image size

A mass $\delta M=\dot{M}({t}_{2}-{t}_{1})$ is added to the star from time t₁ to t₂. Sizes are exaggerated for clarity; there are generally many zones in each region, possibly hundreds in the newly added material. The diagram is shown in both mass coordinate, m, and homology coordinate, q. Before each timestep, MESAstar adjusts the initial mesh resolution by splitting or merging cells based on local gradient conditions, producing an adjusted resolution mesh, which we show at t₁. No simulation time elapses during that process. When constructing the mesh that will represent the star at t₂, MESAstar divides the new mesh into three regions: an inner Lagrangian region, in which the m boundaries of each cell are preserved during the timestep, a transition region, and an outer, homologous region, in which the q boundaries of each cell are the same across the timestep. The result of this operation is an expanded mass domain shown at t₂ in Figure 32.

Time derivatives appearing in the equations for physical evolution are then estimated using first order differences. In the Lagrangian mesh region, the finite difference form of ${(\partial /\partial t)}_{m}$ involves a simple same-cell difference. Similarly, in the homologous mesh region, the finite difference form of ${(\partial /\partial t)}_{q}$ involves a same-cell difference. In most cases, by design, these same-cell differences are for values whose changes, e.g., $\delta \mathrm{ln}T$, are directly available from the iterative solution of the new structure, allowing us to avoid the numerical problems inherent in subtracting two almost identical numbers. In the transition region both m and q coordinates of cells have been modified, so we cannot do a same-cell difference for either ${(\partial /\partial t)}_{m}$ or ${(\partial /\partial t)}_{q}$. Instead, we interpolate values from the model at the start of the step to corresponding locations in m or q at the end of the step.

A smooth and accurate value for ${\epsilon }_{\mathrm{grav}}$ in the transition region is important. To ensure this, the location of the transition region is selected to reduce the differences between the constant m and constant q forms of the time derivative and maintain accurate finite differencing. As a simple mechanism to control these, we limit, in units of cell size, the offset in the interpolation used to translate locations from the beginning to the end of the timestep (B. J. Miles et al. 2015, in preparation). Using the cell size implicitly takes advantage of the limits imposed by mesh controls on the maximum possible magnitude of cell-to-cell changes in key variables, including the variables of interest for ${\epsilon }_{\mathrm{grav}}$ time derivatives.

In order to demonstrate that the interface between the outer homologous region and the inner Lagrangian region provides a smooth profile that is independent of timestep size, we have repeated the test shown in Paper II Section 5.3. This test involves accretion of solar composition material onto a WD at ${10}^{-10}{{\rm{M}}}_{\odot }$ yr⁻¹. We use the same starting model as in Paper II, which was produced by accreting hydrogen-rich material through several hydrogen shell flashes on a $0.6{{\rm{M}}}_{\odot }$ WD with an initial core temperature of about 10⁷ K. As accretion proceeds, the total accumulated accreted mass, ${M}_{\mathrm{acc}}$, increases up to a maximum which causes the hydrogen flash and nova runaway, ${M}_{\mathrm{ign}}$.

Profiles near the transitions region are shown in Figure 33 at the time when ${M}_{\mathrm{acc}}/{M}_{\mathrm{ign}}=0.2$, which had the most severe discontinuity in Paper II. The first panel shows ${{\rm{\nabla }}}_{T}$ and the second panel shows $4\pi {r}^{2}\rho {H}_{{\rm{P}}}\;\times \;{\epsilon }_{\mathrm{grav}}$, where H_P is the pressure scale height. This is the amount of energy being released due to the $T\;{Ds}/{Dt}$ term in the energy equation within a scale height, and has units of luminosity. We have chosen a timestep of 5000 years, which places the homology-Lagrangian transition region in a similar place to the location of the discontinuity in Paper II.

Zoom In Zoom Out Reset image size

The orange curve in Figure 33 was computed using MESA version r4664, as used in Paper II. This displays the discontinuity due to using ${\epsilon }_{\mathrm{grav}}$ from Equations (62) and (63) with only the homologous term and a transition point a factor of five deeper in pressure than $\delta M$. The black curve shows the same simulation with the same timestep for the treatment discussed here, in which ${\epsilon }_{\mathrm{grav},\mathrm{nh}}$ is included and the transition region, indicated with open circles at either end, is placed as described in Section 7.1. We see that, away from the transitions, ${\epsilon }_{\mathrm{grav}}$ is unchanged from the values found in Paper II, in which only the homologous term was used in the exterior. We also show the result for a timestep of 1000 years is indistinguishable on this plot; the ${\epsilon }_{\mathrm{grav}}$ profile differs by a fraction of a percent at the edges of the transition region, and less elsewhere. In the current treatment the profiles are independent of timestep size.

Consider two nuclei $A\equiv (Z,N)$ and $B\equiv (Z-1,N+1)$ that have two states connected by an electron-capture transition

$$\begin{eqnarray}&&A+{e}^{-}\to B+{\nu }_{{\rm{e}}}\end{eqnarray} \tag{ 66 }$$

and beta-decay transition

$$\begin{eqnarray}&&B\to A+{e}^{-}+{\bar{\nu }}_{{\rm{e}}}.\end{eqnarray} \tag{ 67 }$$

The energy difference between the ground states can be written as

$$\begin{eqnarray}{Q}_{{\rm{g}}}=\left\{\begin{array}{ll}({M}_{A}-{M}_{B}){c}^{2} & {\rm{}}\;{\rm{for}}\;{\rm{electron}}\;{\rm{capture}},\\ ({M}_{B}-{M}_{A}){c}^{2} & {\rm{}}\;{\rm{for}}\;{\rm{beta}}\;{\rm{decay}},\end{array}\right.\end{eqnarray} \tag{ 68 }$$

where M_A and M_B are the nuclear rest masses of the ground states. The total energy difference between any two states can be written as

$$\begin{eqnarray}&&{Q}_{{ij}}={Q}_{{\rm{g}}}+{E}_{i}-{E}_{f},\end{eqnarray} \tag{ 69 }$$

where E_i and E_f are the energies of the initial and final states measured relative to the ground state. For the transitions that we consider here, ${Q}_{{\rm{g}}}\lt 0$ and ${Q}_{{ij}}\lt 0$ for electron capture and ${Q}_{{\rm{g}}}\gt 0$ and ${Q}_{{ij}}\gt 0$ for beta decay.

In this section, we use J to represent the nuclear spin. We work in the allowed approximation, which neglects all total lepton angular momentum (L = 0). This restricts us to Fermi transitions, where the total lepton spin is S = 0, and therefore the initial and final nuclear spins are equal (${J}_{i}={J}_{f}$), and Gamow–Teller transitions, where S = 1, and therefore ${J}_{i}={J}_{f},{J}_{f}\pm 1$ (excluding ${J}_{i}={J}_{f}=0$). In both cases, there is no parity change: ${\pi }_{i}{\pi }_{f}=+1$ (e.g., Commins 1973).

The total rate of the process (electron capture or beta decay) is the sum of the individual transition rates from the ith parent state to the jth daughter state, ${\lambda }_{{ij}}$, weighted by the occupation probability of the ith parent state, p_i.

$$\begin{eqnarray}&&\lambda =\displaystyle \sum _{i}{p}_{i}\displaystyle \sum _{j}\displaystyle \frac{\mathrm{ln}2}{{({ft})}_{{ij}}}{\rm{\Phi }}({\mu }_{{\rm{e}}},T,{Q}_{{ij}}),\end{eqnarray} \tag{ 70 }$$

where $({ft})$ is the comparative half-life and can be either measured experimentally or calculated from theoretical weak-interaction nuclear matrix elements. The i-sum is over all parent states and the j-sum is over all daughter states. The occupation probability is

$$\begin{eqnarray}&&{p}_{i}=\displaystyle \frac{(2{J}_{i}+1)\mathrm{exp}\left(-\beta {E}_{i}\right)}{\displaystyle \sum _{k}(2{J}_{k}+1)\mathrm{exp}(-\beta {E}_{k})},\end{eqnarray} \tag{ 71 }$$

where we define $\beta ={({k}_{{\rm{B}}}T)}^{-1}$. The quantity Φ is a phase space factor which depends on the electron chemical potential ${\mu }_{{\rm{e}}}$ (including the electron rest mass), on the temperature T, and the energy difference Q_ij. The value of Φ for electron capture is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{ec}}=\displaystyle \frac{\mathrm{exp}(-\pi \alpha Z)}{{({m}_{{\rm{e}}}{c}^{2})}^{5}}{\displaystyle \int }_{-{Q}_{{ij}}}^{\infty }\displaystyle \frac{{E}_{{\rm{e}}}^{2}{({E}_{{\rm{e}}}+{Q}_{{ij}})}^{2}}{1+\mathrm{exp}[\beta ({E}_{{\rm{e}}}-{\mu }_{{\rm{e}}})]}{{dE}}_{{\rm{e}}},\end{eqnarray} \tag{ 72 }$$

where α is the fine structure constant. For beta decay it is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\beta }=\displaystyle \frac{\mathrm{exp}(-\pi \alpha Z)}{{({m}_{{\rm{e}}}{c}^{2})}^{5}}{\displaystyle \int }_{{m}_{{\rm{e}}}{c}^{2}}^{{Q}_{{ij}}}\displaystyle \frac{{E}_{{\rm{e}}}^{2}{({E}_{{\rm{e}}}-{Q}_{{ij}})}^{2}}{1+\mathrm{exp}[-\beta ({E}_{{\rm{e}}}-{\mu }_{{\rm{e}}})]}{{dE}}_{{\rm{e}}}.\end{eqnarray} \tag{ 73 }$$

Similarly, the total rate of energy loss via neutrinos is

$$\begin{eqnarray}&&{\varepsilon }_{\nu }=\displaystyle \sum _{i}{p}_{i}\displaystyle \sum _{j}\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}\mathrm{ln}2}{{({ft})}_{{ij}}}{\rm{\Psi }}({\mu }_{{\rm{e}}},T,{Q}_{{ij}}).\end{eqnarray} \tag{ 74 }$$

The value of Ψ for electron capture is

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{ec}}=\displaystyle \frac{\mathrm{exp}(-\pi \alpha Z)}{{({m}_{{\rm{e}}}{c}^{2})}^{6}}{\displaystyle \int }_{-{Q}_{{ij}}}^{\infty }\displaystyle \frac{{E}_{{\rm{e}}}^{2}{({E}_{{\rm{e}}}+{Q}_{{ij}})}^{3}}{1+\mathrm{exp}[\beta ({E}_{{\rm{e}}}-{\mu }_{{\rm{e}}})]}{{dE}}_{{\rm{e}}},\end{eqnarray} \tag{ 75 }$$

and for beta decay it is

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\beta }=\displaystyle \frac{\mathrm{exp}(-\pi \alpha Z)}{{({m}_{{\rm{e}}}{c}^{2})}^{6}}{\displaystyle \int }_{{m}_{{\rm{e}}}{c}^{2}}^{{Q}_{{ij}}}\displaystyle \frac{{E}_{{\rm{e}}}^{2}{({E}_{{\rm{e}}}-{Q}_{{ij}})}^{3}}{1+\mathrm{exp}[-\beta ({E}_{{\rm{e}}}-{\mu }_{{\rm{e}}})]}{{dE}}_{{\rm{e}}}.\end{eqnarray} \tag{ 76 }$$

In order to implement these equations in MESA, we rewrite the integrals in terms of Fermi–Dirac integrals, following Appendix A of Schwab et al. (2015). MESA implements fast quadrature routines to evaluate integrals of this form. Each time a weak rate is needed, it is calculated on-the-fly. We discuss the computational cost of this procedure in Section 8.3.

Assuming thermal equilibrium, the energy released by weak reactions depends only on total reaction rate, total neutrino loss rate, energy difference between the nuclei, and the electron chemical potential. Therefore, the total specific heating rate from a reaction is

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ec}}=\displaystyle \frac{{n}_{A}}{\rho }\left[({Q}_{{\rm{g}}}+{\mu }_{{\rm{e}}}){\lambda }_{\mathrm{ec}}-{\varepsilon }_{\nu ,\mathrm{ec}}\right],\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{\epsilon }_{\beta }=\displaystyle \frac{{n}_{B}}{\rho }\left[({Q}_{{\rm{g}}}-{\mu }_{{\rm{e}}}){\lambda }_{\beta }-{\varepsilon }_{\nu ,\beta }\right],\end{eqnarray} \tag{ 78 }$$

where n_A and n_B are the number densities of the species undergoing electron capture and beta decay, respectively, and ρ is the total mass density.

Therefore, given a list of nuclear levels and the $({ft})$-values for the transitions between them, MESA can calculate the rates of electron capture and beta decay and the corresponding energy generation rates. Typically only a few low-lying states and the transitions between them are needed. As an example, Figure 35 shows the rates for the ${}^{23}\mathrm{Na}$–${}^{23}\mathrm{Ne}$ Urca pair.

Zoom In Zoom Out Reset image size

8.1.1. Coulomb Corrections

When ${\mu }_{{\rm{e}}}\lesssim | Q|$, only the few electrons in the tail of the Fermi–Dirac distribution have sufficient energy to overcome the energy gap and capture on A to form B. Thus, the rate of electron capture is small compared to the rate of beta decay, and so isotope B is favored in the equilibrium. When ${\mu }_{{\rm{e}}}\gtrsim | Q|$, there are only a few unoccupied states available to accept the energetic electron from the beta decay. This final state blocking means the rate of beta decay is small compared to the rate of electron capture, and so isotope A will be favored in the equilibrium.

The shift in this equilibrium can have profound consequences when it occurs in stellar interiors. It modifies the composition, reduces the electron fraction, and alters the thermal state of the plasma. We now discuss two applications of our on-the-fly treatment of the weak rates: the Urca process and accretion-induced collapse.

8.2.1. Urca Process

When the ground state to ground state transition is allowed (odd nuclei), the rates of electron capture and beta decay are both significant when ${\mu }_{{\rm{e}}}\approx | {Q}_{{\rm{g}}}|$. Since each reaction produces a neutrino which free-streams out of the star, this can lead to significant cooling. With a total number density of an Urca species ${n}_{U}={n}_{A}+{n}_{B}$, assuming the abundances are given by the detailed balance condition ${n}_{A}{\lambda }_{\mathrm{ec}}+{n}_{B}{\lambda }_{\beta }=0$, the volumetric neutrino cooling rate will be ${n}_{U}C$, where

$$\begin{eqnarray}&&C=\displaystyle \frac{{\varepsilon }_{\nu ,\mathrm{ec}}{\lambda }_{\beta }+{\varepsilon }_{\nu ,\beta }{\lambda }_{\mathrm{ec}}}{{\lambda }_{\beta }+{\lambda }_{\mathrm{ec}}}.\end{eqnarray} \tag{ 79 }$$

In the limit ${k}_{{\rm{B}}}T\ll | Q|$, the maximum value of the Urca cooling rate at a given temperature has a simple form (e.g., Tsuruta & Cameron 1970)

$$\begin{eqnarray}&&{C}_{\mathrm{max}}=\displaystyle \frac{7{\pi }^{4}\mathrm{ln}2}{120}\left(\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{({ft})}\right){\left(\displaystyle \frac{{k}_{{\rm{B}}}T}{{m}_{{\rm{e}}}{c}^{2}}\right)}^{4}{\left(\displaystyle \frac{Q}{{m}_{{\rm{e}}}{c}^{2}}\right)}^{2}\mathrm{exp}(-\pi \alpha Z).\end{eqnarray} \tag{ 80 }$$

Well-sampled rates such as those shown in Figure 35 are necessary to reproduce the correct Urca cooling rates. We illustrate this in Figure 36, which shows C_max for the ${}^{23}\mathrm{Na}$–${}^{23}\mathrm{Ne}$ Urca process for temperatures 10⁸–10⁹ K. The circles show the results using the on-the-fly treatment described in this paper; the squares show the results using the coarse tables of Oda et al. (1994). The dashed line shows the cooling rate expected from Equation (80) which is in excellent agreement with the results of the on-the-fly method. The Urca cooling rates calculated from interpolating in coarse tables severely underestimate the true cooling rate when ${k}_{{\rm{B}}}T\ll | Q|$.

Zoom In Zoom Out Reset image size

Thus, when the Urca process is important, well-resolved weak rates are necessary to correctly capture the temperature evolution of the core (Jones et al. 2013; Toki et al. 2013). Jones et al. (2013) used MESA r3709 along with a denser table described in Toki et al. (2013) to do their work. The Toki et al. (2013) table is not publicly available, so to reproduce the results of Jones et al. (2014) we save a model of an 8.8 ${{\rm{M}}}_{\odot }$ star at $\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\;{\mathrm{cm}}^{-3})=8.95$ from our run with MESA (version r3709) using the Jones et al. (2014) inlists. We then load this model into a newer MESA version (r7503) that has access to the on-the-fly weak rates and evolve this model using a network with only the Urca process reactions. During this phase other nuclear reactions are not important to the central evolution.

Figure 37 shows the central temperature and density of the core. The solid lines show the evolution using the on-the-fly rates, the dashed lines show the results when interpolating in coarse tables. The drops in temperature at $\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\;{\mathrm{cm}}^{-3})\approx 9.1$ and $\approx 9.25$ correspond to cooling from the ${}^{25}\mathrm{Mg}$–${}^{25}\mathrm{Na}$ and ${}^{23}\mathrm{Na}$–${}^{23}\mathrm{Ne}$ Urca pairs, respectively. The corresponding shifts in composition can be clearly seen in the lower panel. These results demonstrate the importance of densely sampled weak rates to the evolution of the core.

Zoom In Zoom Out Reset image size

8.2.2. Accretion-induced Collapse

When the ground state to ground state transition is forbidden (even nuclei), the first transition to become significant is typically an allowed transition into an excited state. In these cases, the beta decays from the daughter ground state are blocked and decays from daughter excited states are strongly suppressed by the Boltzmann factor. Therefore significant cooling via the Urca process does not occur. Instead, since the captures are preferentially to an excited state, significant heating occurs via gamma-ray emission as level populations relax to a thermal distribution.

Two important capture chains occur in oxygen-neon-magnesium (ONeMg) cores: ${}^{24}\mathrm{Mg}{\to }^{\;24}\mathrm{Na}{\to }^{\;24}\mathrm{Ne}$ and ${}^{20}\mathrm{Ne}{\to }^{\;20}{\rm{F}}{\to }^{\;20}{\rm{O}}$. For these sequences of captures, the excess electron energy is thermalized. These are the key reactions in electron-capture supernovae and the accretion-induced collapse (AIC) of ONeMg WDs (e.g., Miyaji et al. 1980). As the degenerate core approaches the Chandrasekhar mass, the electron captures remove the pressure support and heat the plasma. Figure 38 shows the evolution of a cold ONeMg WD (${X}_{{\rm{O}}}=0.5$, ${X}_{\mathrm{Ne}}=0.45$, ${X}_{\mathrm{Mg}}=0.05$) accreting at $\dot{M}={10}^{-6}\;{{\rm{M}}}_{\odot }\;{\mathrm{yr}}^{-1}$. The solid lines show the evolution using the on-the-fly rates described in this section; the dashed lines show the results when interpolating in coarse tables. When using the coarse tables, the electron captures on ${}^{24}\mathrm{Mg}$ do not occur until approximately a factor of two larger density. At this greater density, the energy deposition from each capture is higher and this leads to a large temperature change due to the A = 24 captures alone. In contrast, the on-the-fly rates show the behavior demonstrated in previous studies of this evolution that did not use sparse tables (e.g., Miyaji & Nomoto 1987): the A = 24 captures heat the plasma and accelerate the contraction; the A = 20 captures, due to the higher ${}^{20}\mathrm{Ne}$ abundance and a higher energy release per capture, cause a thermal runaway and the formation of an oxygen deflagration (Schwab et al. 2015).

Zoom In Zoom Out Reset image size

MESA provides the nuclear data used in the calculation of the reactions specifically discussed in this section (Tilley et al. 1998; Firestone 2007a, 2007b, 2009; Shamsuzzoha Basunia 2011; Martínez-Pinedo et al. 2014). To consider additional reactions, a list of nuclear levels and $({ft})$-values must be specified.

The expressions in Section 8.1 assume degenerate, relativistic electrons. As ${\mu }_{{\rm{e}}}$ increases, additional transitions to higher energy states of the daughter nuclei and must be included. At higher temperatures, excited states of the parent nucleus will begin to be thermally populated and captures or decays from those states and must be included. At temperatures and densities where the composition approaches NSE, these methods are particularly inappropriate, as it is necessary to consider large pools of isotopes (Juodagalvis et al. 2010).

We now describe the formalism and assumptions underlying the approach to diffusion currently present in MESA. This is followed by a discussion of the framework for numerical implementation of this formalism provided by Thoul et al. (1994) and key modifications present in the current version of the MESA diffusion routine.

9.1.1. Burgers Equations and the Low Density Limit

The Burgers equations for diffusion in an ionized plasma are derived using the Boltzmann equation for the distribution function ${F}_{s}({\boldsymbol{x}},{\boldsymbol{\xi }},t)$ for particles of type s

$$\begin{eqnarray}&&\displaystyle \frac{\partial {F}_{s}}{\partial t}+\displaystyle \sum _{i}{\xi }_{i}\displaystyle \frac{\partial {F}_{s}}{\partial {x}_{i}}+\displaystyle \sum _{i}\displaystyle \frac{{f}_{{si}}}{{m}_{s}}\displaystyle \frac{\partial {F}_{s}}{\partial {\xi }_{i}}={\left(\displaystyle \frac{{{dF}}_{s}}{{dt}}\right)}_{\mathrm{collision}},\end{eqnarray} \tag{ 81 }$$

where x_i are the components of the position vector, ${\xi }_{i}$ are the components of the velocity vector, f_si are components of the forces on particles of type s, and m_s is the mass for those particles. Throughout this section, the indices s and t refer to particle species, while i and j are used to index other quantities such as spatial components of vectors.

Burgers adopts the 13-moment approximation due to Grad (1949) as a closure scheme for taking moments of the Boltzmann equation. Burgers also assumes an approximately Maxwellian distribution function

$$\begin{eqnarray}&&{F}_{s}=\displaystyle \frac{{n}_{s}}{{\pi }^{3/2}{a}_{s}^{3}}\mathrm{exp}\left(\displaystyle \frac{-{c}_{s}^{2}}{{a}_{s}^{2}}\right)(1+{\phi }_{s}),\end{eqnarray} \tag{ 82 }$$

where ${a}_{s}={(2{k}_{{\rm{B}}}T/{m}_{s})}^{1/2}$, ${c}_{{si}}={\xi }_{i}-{u}_{{si}}$ represents the components of the deviation of the velocity from the mean flow velocity ${{\boldsymbol{u}}}_{s}$ of the species, and

$$\begin{eqnarray}&&{\phi }_{s}=\displaystyle \sum _{i,j}{B}_{{sij}}{c}_{{si}}{c}_{{sj}}+\displaystyle \sum _{i}{C}_{{si}}\left({c}_{s}^{2}-\displaystyle \frac{5}{2}{a}_{s}^{2}\right){c}_{{si}}\end{eqnarray} \tag{ 83 }$$

is the small deviation (${\phi }_{s}\ll 1$) from the Maxwellian distribution. The coefficients B_sij and C_si are defined such that the distribution function has a total of 13 free parameters corresponding to the 13 moments of the closure scheme (see Burgers 1969).

Burgers derives the collision integrals (${S}_{{st}}^{(l)}$) and cross-sections (${{\rm{\Sigma }}}_{{st}}^{({lj})}$) that result from taking moments of the right-hand side of the Boltzmann equation

$$\begin{eqnarray}&&{S}_{{st}}^{(l)}=2\pi {\displaystyle \int }_{0}^{\infty }(1-{\mathrm{cos}}^{l}{\chi }_{{st}})b\;{db},\end{eqnarray} \tag{ 84 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{st}}^{({lj})}=\displaystyle \frac{4\pi }{{\pi }^{3/2}}{\displaystyle \int }_{0}^{\infty }{dv}\;\mathrm{exp}\left(\displaystyle \frac{-{v}^{2}}{{\alpha }_{{st}}^{2}}\right)\displaystyle \frac{{v}^{2j+3}}{{\alpha }_{{st}}^{2j+4}}{S}_{{st}}^{(l)},\end{eqnarray} \tag{ 85 }$$

where ${\alpha }_{{st}}^{2}=2{k}_{{\rm{B}}}T/{\mu }_{{st}}$, ${\mu }_{{st}}={m}_{s}{m}_{t}/({m}_{s}+{m}_{t})$, v represents the relative velocity of colliding particles, and the angle of deviation ${\chi }_{{st}}$ is a function of both v and the impact parameter b that depends on the physics of the two-particle interaction between colliding particles in the gas. Burgers then defines the dimensionless coefficients z_st, ${z}_{{st}}^{\prime }$, ${z}_{{st}}^{\prime\prime }$, and ${z}_{{st}}^{\prime \prime \prime }$, along with resistance coefficients ( K_st) in terms of the collision integrals:

$$\begin{eqnarray}{K}_{{st}}=\;{K}_{{ts}} & = & \displaystyle \frac{2}{3}{n}_{s}{n}_{t}{\mu }_{{st}}{\alpha }_{{st}}{{\rm{\Sigma }}}_{{st}}^{(11)},\\ {{\rm{\Sigma }}}_{{st}}^{(12)}/{{\rm{\Sigma }}}_{{st}}^{(11)} & = & \displaystyle \frac{5}{2}(1-{z}_{{st}}),\\ {{\rm{\Sigma }}}_{{st}}^{(13)}/{{\rm{\Sigma }}}_{{st}}^{(11)} & = & \displaystyle \frac{25}{4}-\displaystyle \frac{25}{2}{z}_{{st}}+\displaystyle \frac{5}{2}{z}_{{st}}^{\prime },\\ {{\rm{\Sigma }}}_{{st}}^{(22)}/{{\rm{\Sigma }}}_{{st}}^{(11)} & = & {z}_{{st}}^{\prime\prime },\\ {{\rm{\Sigma }}}_{{st}}^{(23)}/{{\rm{\Sigma }}}_{{st}}^{(11)} & = & {z}_{{st}}^{\prime \prime \prime }.\end{eqnarray} \tag{ 86 }$$

In the "single-fluid picture" the diffusion velocities are defined with reference to the mean velocity of the gas as a whole (${\boldsymbol{u}}$), rather than with respect to the mean species velocity (${{\boldsymbol{u}}}_{s}$):

$$\begin{eqnarray}&&{u}_{{si}}=\displaystyle \frac{1}{{n}_{s}}\int d{\boldsymbol{\xi }}\;{\xi }_{i}{F}_{s},\quad {\boldsymbol{u}}=\displaystyle \frac{1}{\rho }\displaystyle \sum _{s}{\rho }_{s}{{\boldsymbol{u}}}_{s},\quad {{\boldsymbol{w}}}_{s}={{\boldsymbol{u}}}_{s}-{\boldsymbol{u}}.\end{eqnarray} \tag{ 87 }$$

Burgers defines residual heat flow vectors

$$\begin{eqnarray}&&{r}_{{si}}=\left(\displaystyle \frac{{m}_{s}}{2{n}_{s}{k}_{{\rm{B}}}T}\int {d}^{3}{\boldsymbol{\xi }}\;({\xi }_{i}-{u}_{i}){| {\boldsymbol{\xi }}-{\boldsymbol{u}}| }^{2}{F}_{s}\right)-\displaystyle \frac{5}{2}{w}_{{si}}.\end{eqnarray} \tag{ 88 }$$

As shown in Section 18 of Burgers (1969) if we assume $| {{\boldsymbol{w}}}_{s}| \ll {a}_{s}$ and the absence of magnetic fields, the basic equations of diffusion are

$$\begin{eqnarray}{\rm{\nabla }}{p}_{s}-{\rho }_{s}{\boldsymbol{g}}-{\rho }_{{es}}{\boldsymbol{E}} & = & \displaystyle \sum _{t\ne s}\;{K}_{{st}}({{\boldsymbol{w}}}_{t}-{{\boldsymbol{w}}}_{s})\\ & & +\displaystyle \sum _{t\ne s}\;{K}_{{st}}{z}_{{st}}\displaystyle \frac{{m}_{t}{{\boldsymbol{r}}}_{s}-{m}_{s}{{\boldsymbol{r}}}_{t}}{{m}_{s}+{m}_{t}},\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}\displaystyle \frac{5}{2}{n}_{s}{k}_{{\rm{B}}}{\rm{\nabla }}T & = & -\displaystyle \frac{2}{5}\;{K}_{{ss}}{z}_{{ss}}^{\prime\prime }{{\boldsymbol{r}}}_{s}-\displaystyle \frac{5}{2}\displaystyle \sum _{t\ne s}\;{K}_{{st}}{z}_{{st}}\displaystyle \frac{{m}_{t}}{{m}_{s}+{m}_{t}}({{\boldsymbol{w}}}_{t}-{{\boldsymbol{w}}}_{s})\\ & & -\displaystyle \sum _{t\ne s}\;{K}_{{st}}\left[\displaystyle \frac{3{m}_{s}^{2}+{m}_{t}^{2}{z}_{{st}}^{\prime }}{{({m}_{s}+{m}_{t})}^{2}}+\displaystyle \frac{4}{5}\displaystyle \frac{{m}_{s}{m}_{t}}{{({m}_{s}+{m}_{t})}^{2}}{z}_{{st}}^{\prime\prime }\right]{{\boldsymbol{r}}}_{s}\\ & & +\displaystyle \sum _{t\ne s}\;{K}_{{st}}\displaystyle \frac{{m}_{s}{m}_{t}}{{({m}_{s}+{m}_{t})}^{2}}\left(3+{z}_{{st}}^{\prime }-\displaystyle \frac{4}{5}{z}_{{st}}^{\prime\prime }\right){{\boldsymbol{r}}}_{t},\end{eqnarray} \tag{ 90 }$$

where ${\boldsymbol{E}}$ is the quasi-static electric field and ${\rho }_{{es}}$ is the average charge density of species s. These equations are still general, with the form of the resistance coefficients not yet fully specified. The physics of the particular types of interactions within ideal gases is fully contained in the coefficients K_st, z_st, ${z}_{{st}}^{\prime }$, ${z}_{{st}}^{\prime\prime }$, and ${z}_{{st}}^{\prime \prime \prime }$.

For ionized gases, the resistance coefficients require evaluation of collision integrals that diverge for a pure Coulomb potential. However, since the two-particle interaction potential is only truly applicable on short length scales, an integration cutoff or screened potential is commonly adopted. Burgers chooses to calculate resistance coefficients using a pure Coulomb potential truncated at the Debye radius

$$\begin{eqnarray}&&{R}_{{\rm{D}}}={\left(4\pi \displaystyle \sum _{s}\displaystyle \frac{{n}_{s}{Z}_{s}^{2}{e}^{2}}{{k}_{{\rm{B}}}T}\right)}^{-1/2},\end{eqnarray} \tag{ 91 }$$

which is assumed to be much larger than the inter-ion spacing. Indeed, for a plasma of one species, ${R}_{{\rm{D}}}/{\lambda }_{\mathrm{ion}}={(3{\rm{\Gamma }})}^{-1/2}$. Applying this form of interaction to the collision integrals, the l = 1 integrals defined in Equation (84) can be evaluated (Baalrud & Daligault 2014)

$$\begin{eqnarray}&&{S}_{{st}}^{(1)}=\displaystyle \frac{2\pi {R}_{{\rm{D}}}^{2}{\alpha }_{{st}}^{4}}{{{\rm{\Lambda }}}_{{st}}^{2}{v}^{4}}\mathrm{ln}\left[1+{{\rm{\Lambda }}}_{{st}}^{2}{\left(\displaystyle \frac{v}{{\alpha }_{{st}}}\right)}^{4}\right],\end{eqnarray} \tag{ 92 }$$

where ${{\rm{\Lambda }}}_{{st}}={\mu }_{{st}}{\alpha }_{{st}}^{2}{R}_{{\rm{D}}}/({Z}_{s}{Z}_{t}{e}^{2})$. In order to perform the integral in Equation (85), Burgers notes that the dependence of ${S}_{{st}}^{(l)}$ on v inside the logarithmic term is weak, so that we can replace v² there with its average value $\langle {v}^{2}\rangle =3{k}_{{\rm{B}}}T/{\mu }_{{st}}$. Assuming a very dilute plasma, so that ${{\rm{\Lambda }}}_{{st}}^{2}{\langle {v}^{2}\rangle }^{2}/{\alpha }_{{st}}^{4}\gg 1$, Burgers then writes

$$\begin{eqnarray}&&{S}_{{st}}^{(1)}\approx \displaystyle \frac{4\pi {R}_{{\rm{D}}}^{2}{\alpha }_{{st}}^{4}}{{{\rm{\Lambda }}}_{{st}}^{2}{v}^{4}}\mathrm{ln}\left(\displaystyle \frac{3{k}_{{\rm{B}}}{{TR}}_{{\rm{D}}}}{{Z}_{s}{Z}_{t}{e}^{2}}\right),\end{eqnarray} \tag{ 93 }$$

and the final result for the resistance coefficients follows as

$$\begin{eqnarray}&&{K}_{{st}}\approx \displaystyle \frac{16\sqrt{\pi }}{3}\displaystyle \frac{{n}_{s}{n}_{t}{Z}_{s}^{2}{Z}_{t}^{2}{e}^{4}}{{\mu }_{{st}}{\alpha }_{{st}}^{3}}\mathrm{ln}\left(\displaystyle \frac{3{k}_{{\rm{B}}}{{TR}}_{{\rm{D}}}}{{Z}_{s}{Z}_{t}{e}^{2}}\right),\end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray}&&{z}_{{st}}=\displaystyle \frac{3}{5},\qquad {z}_{{st}}^{\prime }=\displaystyle \frac{13}{10},\qquad {z}_{{st}}^{\prime\prime }=2,\qquad {z}_{{st}}^{\prime \prime \prime }=4.\end{eqnarray} \tag{ 95 }$$

With these coefficients now fully specified, Burgers diffusion equations along with constraints such as charge neutrality and current neutrality form a closed set of equations, which can be solved for ${{\boldsymbol{w}}}_{s}$, ${{\boldsymbol{r}}}_{s}$, ${\boldsymbol{E}}$, and ${\boldsymbol{g}}$ from the input of a stellar profile.

9.1.2.  MESA's Implementation of Thoul et al.'s Approach

The diffusion routine originally implemented in MESA was based on the work of Thoul et al. (1994). They start with the Burgers equations, written in a compact notation following Noerdlinger (1977, 1978) that is equivalent to Equations (89) and (90) in one dimension. However, the approach of Thoul et al. (1994) differs from Burgers' original treatment in one important respect: the resistance coefficients are based on a modified result for the collision integrals. They follow Equation (95) for the various z_st coefficients, which uses a pure Coulomb potential with a cutoff at the Debye length, but the K_st coefficients were derived from an alternative fitting of the Coulomb logarithms introduced by Iben & MacDonald (1985). For these coefficients, they define $\lambda =\mathrm{max}({R}_{{\rm{D}}},{\lambda }_{\mathrm{ion}})$, and use

$$\begin{eqnarray}{K}_{{st}} & = & \displaystyle \frac{16\sqrt{\pi }}{3}\displaystyle \frac{{n}_{s}{n}_{t}{Z}_{s}^{2}{Z}_{t}^{2}{e}^{2}}{{\mu }_{{st}}{\alpha }_{{st}}^{3}}\\ & & \;\times \;\displaystyle \frac{1.6249}{2}\mathrm{ln}\left[1+0.18769{\left(\displaystyle \frac{4{k}_{{\rm{B}}}T\lambda }{{Z}_{s}{Z}_{t}{e}^{2}}\right)}^{1.2}\right].\end{eqnarray} \tag{ 96 }$$

This expression is a fit to the numerical results of Fontaine & Michaud (1979), motivated by WD conditions where Burgers' approximations for dealing with Equation (92) are not valid (${\rm{\Gamma }}\gt 1$). Since this fit focuses on the strong coupling regime, and differs from Equation (94), these results can be incorrect in the limit of a dilute plasma as we discuss later. Nevertheless, Thoul et al. (1994) elected to use Equation (96) under all conditions, since it provides an approximately correct solution in a convenient closed form.

Using Equations (89) and (90) along with the constraints of current neutrality (${\displaystyle \sum }_{s}{\rho }_{{es}}{w}_{s}=0$) and local mass conservation (${\displaystyle \sum }_{s}{\rho }_{s}{w}_{s}=0$), Thoul et al. (1994) express an entire closed system of equations in a dimensionless matrix form suitable for numerical evaluation:

$$\begin{eqnarray}&&\displaystyle \frac{p}{{K}_{0}}\left({\alpha }_{i}\displaystyle \frac{d\mathrm{ln}p}{{dr}}+{\nu }_{i}\displaystyle \frac{d\mathrm{ln}T}{{dr}}+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j=1}{{\rm{j}}\ne {\rm{e}}}}^{S}{\gamma }_{{ij}}\displaystyle \frac{d\mathrm{ln}{C}_{j}}{{dr}}\right)=\displaystyle \sum _{j=1}^{2S+2}{{\rm{\Delta }}}_{{ij}}{W}_{j},\end{eqnarray} \tag{ 97 }$$

where S is the total number of species in the gas (including electrons) and ${C}_{j}={n}_{j}/{n}_{{\rm{e}}}$ is the concentration of the jth species. Consult Thoul et al. (1994) for definitions of K₀, ${\alpha }_{i}$, ${\nu }_{i}$, ${\gamma }_{{ij}}$, and ${{\rm{\Delta }}}_{{ij}}$. The definition of W_j is

$$\begin{eqnarray}{W}_{j}=\left\{\begin{array}{ll}{w}_{j} & {\rm{for}}\;{\rm{}}j=1\ldots S,\\ {r}_{j} & {\rm{for}}\;{\rm{}}j=S+1\ldots 2S,\\ {K}_{0}^{-1}{n}_{{\rm{e}}}{eE} & {\rm{for}}\;{\rm{}}j=2S+1,\\ {K}_{0}^{-1}{n}_{{\rm{e}}}{m}_{{\rm{p}}}g & {\rm{for}}\;{\rm{}}j=2S+2.\end{array}\right.\end{eqnarray} \tag{ 98 }$$

This is the vector containing the unknown quantities solved for after specifying K₀, ${\alpha }_{i}$, ${\nu }_{i}$, ${\gamma }_{{ij}}$, and ${{\rm{\Delta }}}_{{ij}}$. The routine provided by Thoul et al. (1994) inverts Equation (97) for one term in the left-hand side at a time so as to find the "generalized diffusion coefficients," which can be used to construct diffusion velocities or contributions from pressure, temperature, or concentrations individually.

9.1.3. Modified Coefficients and Radiative Levitation as Implemented by Hu et al.

Hu et al. (2011) extend the methods of Thoul et al. (1994) by introducing some key modifications. First, they include an extra force term due to radiative levitation, so that Equation (89) becomes

$$\begin{eqnarray}&&\displaystyle \frac{{{dp}}_{s}}{{dr}}+{\rho }_{s}(g-{g}_{\mathrm{rad},s})-{n}_{s}{\bar{Z}}_{s}{eE}\\ &&\qquad =\displaystyle \sum _{t\ne s}\;{K}_{{st}}({w}_{t}-{w}_{s})+\displaystyle \sum _{t\ne s}\;{K}_{{st}}{z}_{{st}}\displaystyle \frac{{m}_{t}{r}_{s}-{m}_{s}{r}_{t}}{{m}_{s}+{m}_{t}},\end{eqnarray} \tag{ 99 }$$

where ${g}_{\mathrm{rad},s}$ refers to the radiative acceleration on species s. ${\bar{Z}}_{s}$ is the average charge of species s, allowing an account of partial ionization so that ${n}_{s}{\bar{Z}}_{s}e={\rho }_{{es}}$. They do not modify Equation (90).¹⁵

In contrast to Thoul's original routine, Hu et al. (2011) use the resistance coefficients from Paquette et al. (1986), which were generated based on substantial improvements to Fontaine & Michaud (1979). In evaluating the collision integrals, Paquette et al. (1986) use a screened Coulomb potential of the form

$$\begin{eqnarray}&&{V}_{{st}}(r)={\bar{Z}}_{s}{\bar{Z}}_{t}{e}^{2}\displaystyle \frac{\mathrm{exp}(-r/\lambda )}{r},\end{eqnarray} \tag{ 100 }$$

where, once again, $\lambda =\mathrm{max}({R}_{{\rm{D}}},{\lambda }_{\mathrm{ion}})$. As we note below, this choice of λ makes a substantial difference in strongly coupled plasmas, where the Debye radius no longer corresponds to a distance at which other nearby charged particles can significantly screen the Coulomb field. After setting up the algebra for a matrix solution very similar to that of Thoul et al. (1994), Hu et al. (2011) solve for the vector W_j (as defined in Equation (98)) appearing in the equation

$$\begin{eqnarray}&&\displaystyle \frac{p}{{K}_{0}}\left(-\displaystyle \frac{{\alpha }_{i}{m}_{i}{g}_{\mathrm{rad},i}}{{k}_{{\rm{B}}}T}+{\alpha }_{i}\displaystyle \frac{d\mathrm{ln}p}{{dr}}\right.\\ &&\quad +\left.{\nu }_{i}\displaystyle \frac{d\mathrm{ln}T}{{dr}}+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j=1}{{\rm{j}}\ne {\rm{e}}}}^{S}{\gamma }_{{ij}}\displaystyle \frac{d\mathrm{ln}{C}_{j}}{{dr}}\right)=\displaystyle \sum _{j=1}^{2S+2}{{\rm{\Delta }}}_{{ij}}{W}_{j}.\end{eqnarray} \tag{ 101 }$$

Many of the quantities appearing in this equation are defined differently than in Thoul et al. (1994); see Hu et al. (2011) for details. We can also solve this equation directly for the vector W_j to obtain

$$\begin{eqnarray}&&\displaystyle \frac{{W}_{2S+1}}{{W}_{2S+2}}=\displaystyle \frac{{K}_{0}^{-1}{n}_{{\rm{e}}}{eE}}{{K}_{0}^{-1}{n}_{{\rm{e}}}{m}_{{\rm{p}}}g}=\displaystyle \frac{{eE}}{{m}_{{\rm{p}}}g},\end{eqnarray} \tag{ 102 }$$

the strength of the electric field relative to gravity.

In some simple cases, Burgers equations can be solved to yield an analytic expression for the electric field, providing a useful test for MESA. Starting directly with his diffusion equations, Burgers (1969) arrives at the following expressions for a pure plasma of electrons along with one species of ions (charge Ze):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\nabla }}{p}_{{\rm{e}}}+{n}_{{\rm{e}}}e{\boldsymbol{E}}}{Z\;{K}_{0}}={\boldsymbol{w}}+\displaystyle \frac{3}{5}{{\boldsymbol{r}}}_{{\rm{e}}},\end{eqnarray} \tag{ 103 }$$

$$\begin{eqnarray}&&\displaystyle \frac{5}{2}\displaystyle \frac{{n}_{{\rm{e}}}{k}_{{\rm{B}}}{\rm{\nabla }}T}{Z\;{K}_{0}}=-\displaystyle \frac{3}{2}{\boldsymbol{w}}-\left(\displaystyle \frac{2}{5}\displaystyle \frac{{K}_{\mathrm{ee}}{z}_{\mathrm{ee}}^{\prime\prime }}{{K}_{\mathrm{ie}}}+{z}_{\mathrm{ie}}^{\prime }\right){{\boldsymbol{r}}}_{{\rm{e}}},\end{eqnarray} \tag{ 104 }$$

where ${\boldsymbol{w}}={{\boldsymbol{w}}}_{{\rm{i}}}-{{\boldsymbol{w}}}_{{\rm{e}}}$. For a plasma with only one ion species in diffusion equilibrium, the constraints of current neutrality and local mass conservation give ${\boldsymbol{w}}=0$. In the case of a pure hydrogen plasma, $p=2{p}_{{\rm{e}}}$, and in hydrostatic equilibrium ${\rm{\nabla }}{p}_{{\rm{e}}}={\rm{\nabla }}p/2=\rho {\boldsymbol{g}}/2$. Hence, we can solve the above set of equations to find

$$\begin{eqnarray}&&e{\boldsymbol{E}}=-\displaystyle \frac{1}{2}{m}_{{\rm{p}}}{\boldsymbol{g}}-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{2}{5}\displaystyle \frac{{K}_{\mathrm{ee}}{z}_{\mathrm{ee}}^{\prime\prime }}{{K}_{\mathrm{ie}}}+{z}_{\mathrm{ie}}^{\prime }\right)}^{-1}{k}_{{\rm{B}}}{\rm{\nabla }}T.\end{eqnarray} \tag{ 105 }$$

The coefficient for the temperature gradient term depends directly on the nature of the resistance coefficients in the Burgers formalism, so different models of Coulomb collisions in ionized plasma will lead to different results for the electric field.

As a slight generalization of Equation (105) in one dimension, we write

$$\begin{eqnarray}&&\displaystyle \frac{{eE}}{{m}_{{\rm{p}}}g}=\displaystyle \frac{1}{2}-{\alpha }_{e}\displaystyle \frac{{k}_{{\rm{B}}}}{{m}_{{\rm{p}}}g}\displaystyle \frac{{dT}}{{dr}}.\end{eqnarray} \tag{ 106 }$$

If we calculate the coefficient ${\alpha }_{e}$ using the Burgers' formalism with Equations (95) and (96), we find

$$\begin{eqnarray}&&{\alpha }_{e}=\displaystyle \frac{3}{2}{\left(\displaystyle \frac{2}{5}\displaystyle \frac{{K}_{\mathrm{ee}}{z}_{\mathrm{ee}}^{\prime\prime }}{{K}_{\mathrm{ie}}}+{z}_{\mathrm{ie}}^{\prime }\right)}^{-1}=0.804.\end{eqnarray} \tag{ 107 }$$

A comparable analytic expression for the electric field is provided by Roussel-Dupré (1981), who applies a Boltzmann–Fokker–Planck approach to finding diffusion coefficients for trace elements in hydrogen plasma. His treatment of diffusion is more precise than the Burgers' formalism, but has the limitation of only being applicable in the case of nearly pure hydrogen with a diffusing trace element. His result for the electric field matches the form of Equation (106) with the coefficient ${\alpha }_{e}=0.703$. This provides another useful point of comparison in the specific case of nearly pure hydrogen plasmas. Below we use this analytic expression as a test of the updated resistance coefficients employed by Hu et al. (2011).

We have constructed several simple MESA test cases in order to illustrate the effects of radiative levitation and different resistance coefficients. Where possible, we compare MESA output to corresponding analytic expressions.

9.3.1. Electric Fields

By default, MESA uses the resistance coefficients provided by Paquette et al. (1986), but it can also use the resistance coefficients defined by Iben & MacDonald (1985), given here in Equation (96). In the case of a pure hydrogen star, the coefficients given in Equation (96) lead directly to Equation (107), so these coefficients are especially useful in performing simple comparisons of MESA output to a corresponding analytic expression. Due to the complicated numerical methods used to obtain the resistance coefficients of Paquette et al. (1986), it is not possible to write down a directly corresponding closed form analytic expression for the electric field, but results based on these more precise calculations compare favorably to those of Roussel-Dupré (1981) in the case of a pure hydrogen plasma. Starting with the MESA test suite, we constructed a solar mass pure hydrogen star, and we ran just long enough to turn on the diffusion routine and gather output for electric and gravitational fields. For such a star, we can compare MESA results for the electric field directly to the analytic expression given in Equation (106), with ${\alpha }_{e}=0.804$ in the solution of Burgers (1969) and ${\alpha }_{e}=0.703$ for Roussel-Dupré (1981).

Figure 39 plots the result of Equation (106) for both values of ${\alpha }_{e}$, along with the results from the diffusion routine (Equation (102)) for each type of resistance coefficients available in MESA. As expected, the curve calculated from the MESA diffusion routine output using the resistance coefficients of Iben & MacDonald (1985) closely matches the analytic expression with ${\alpha }_{e}=0.804$ as calculated by Burgers (1969) using his similar coefficients. When using the more detailed numerical calculations for the resistance coefficients provided by Paquette et al. (1986), the diffusion routine output closely resembles the more precise analytic calculation given by Roussel-Dupré (1981).

Zoom In Zoom Out Reset image size

The Sun provides another interesting test case for comparing the effects of using different resistance coefficients. An example solar model from the MESA test suite was run with different choices of the resistance coefficients. Figure 40 shows a slight difference between the electric field strengths relative to gravity given by the Paquette et al. (1986) coefficients and those by Iben & MacDonald (1985).

Zoom In Zoom Out Reset image size

9.3.2. Gravitational Fields

The MESA diffusion routine treats both the electric field and local gravitational acceleration as unknown quantities. MESA records the quantity ${W}_{2S+2}$ (Equation (98)), used to calculate the gravitational acceleration from the diffusion routine:

$$\begin{eqnarray}&&{g}_{\mathrm{diff}}=\displaystyle \frac{{K}_{0}{W}_{2S+2}}{{n}_{{\rm{e}}}{m}_{{\rm{p}}}}.\end{eqnarray} \tag{ 108 }$$

This expression for ${g}_{\mathrm{diff}}$ is independent of the simpler expression for local gravitational acceleration ${g}_{\mathrm{Gauss}}={Gm}/{r}^{2}$. Figure 41 compares g_Gauss and g_diff for a typical profile found using the example solar model from the MESA test suite. In Figure 41, a profile from a star of larger mass ($M=1.5\;{{\rm{M}}}_{\odot }$) shows disagreement between the gravity outputs in the convective core. The Burgers formalism assumes heat transfer that is correlated with temperature gradients through the residual heat flow vectors defined in Equation (88). This assumption breaks down when most of energy is transported by convection; however, the effects of diffusion in this region are completely overwhelmed by convective mixing and are therefore inconsequential.

Zoom In Zoom Out Reset image size

9.3.3. Radiative Levitation

MESA's implementation of radiative levitation is based on Hu et al. (2011). Figure 42 shows an abundance profile of a subdwarf B star model produced by MESA, where radiative levitation is responsible for the presence of ${}^{56}\mathrm{Fe}$, ${}^{58}\mathrm{Ni}$, and other metals near the surface (as also seen in Figure 3 of Hu et al. 2011).

Zoom In Zoom Out Reset image size

9.3.4. WD Sedimentation

In a cooling WD, diffusion governs sedimentation over long timescales. The assumptions behind the formalism of the Burgers equations do not hold under WD conditions.

• 1.  
  The Burgers equations assume all particle species satisfy an ideal gas equation of state. In the context of a degenerate WD both electrons and ions violate this assumption.
• 2.  
  The very dense, strongly coupled (${\rm{\Gamma }}\gt 1$) conditions of a WD call into question the validity of the two-particle scattering picture used to calculate the ion resistance coefficients.

Nevertheless, for lack of a better option, previous studies have relied on the Burgers equations with the coefficients of Paquette et al. (1986). For example, see Córsico et al. (2002).

Figure 43 shows an abundance profile produced by MESA for a CO WD after 4 Gyr of evolution, where diffusion governs sedimentation in the outer layers. The vertical lines in Figure 43 mark the outer boundaries of regions where the two concerns listed above become significant. Nearly all of the WD resides inside at least one of these regimes, and much of the interesting diffusion sedimentation occurs inside regions that are both significantly coupled and highly degenerate. Thus, improvements to the treatment of diffusion are clearly necessary before we are able to describe diffusion in WDs adequately. This MESA run turns off diffusion for ${\rm{\Gamma }}\geqslant 50$, where we expect strong coupling to substantially modify the underlying equations.

Zoom In Zoom Out Reset image size

The validity of the Boltzmann approach becomes questionable as ${\rm{\Gamma }}\gt 1$ and the ions become a liquid. Bildsten & Hall (2001) estimated the diffusion coefficient in this liquid regime by using the Stokes–Einstein relation. However, for a broad-based code such as MESA, we need to implement diffusion into the ${\rm{\Gamma }}\gt 1$ regime in a manner that allows for a smooth transition between coupling regimes.

Paquette et al. (1986) successfully described diffusion in a regime of intermediate coupling through the use of screened potentials, which are a way to account for the collective nature of interactions in a dense plasma. Though there is no rigorous reason to expect that a formalism based on the two-particle scattering picture should work well as ${\rm{\Gamma }}\to 1$, their comparison to simulations verified that this description of diffusion is very accurate for ${\rm{\Gamma }}\lesssim 1$.

Can these approximations be extrapolated to the strongly coupled regime of ${\rm{\Gamma }}\gt 1$? Baalrud & Daligault (2013) provide a method for numerically calculating resistance coefficients using a hypernetted chain (HNC) approximation from effective potentials. Figure 44 compares their HNC results (diamonds) to their Molecular Dynamics (MD) simulations of a one-component plasma (OCP, circles) for the self-diffusion coefficient ${D}^{*}$, defined by

$$\begin{eqnarray}&&{D}^{*}=\displaystyle \frac{D}{{\lambda }_{\mathrm{ion}}^{2}{\omega }_{{\rm{p}}}},\end{eqnarray} \tag{ 109 }$$

where ${\omega }_{{\rm{p}}}$ is the plasma frequency and $D=2{D}_{{ss}}^{(2)}$ (the factor of 2 in this definition ensures that if we redefine species s in terms of two subspecies s₁ and s₂, then $D={D}_{{s}_{1}{s}_{2}}^{(2)}$). The general expression for the interdiffusion coefficient is

$$\begin{eqnarray}&&{D}_{{st}}^{(2)}=\displaystyle \frac{{n}_{s}{n}_{t}}{{n}_{s}+{n}_{t}}\displaystyle \frac{{k}_{{\rm{B}}}T}{{K}_{{st}}(1-{\rm{\Delta }})},\end{eqnarray} \tag{ 110 }$$

where the $1-{\rm{\Delta }}$ term in the denominator accounts for a second order correction that can be defined using

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{{(2{{\rm{\Sigma }}}_{{st}}^{(12)}-5{{\rm{\Sigma }}}_{{st}}^{(11)})}^{2}/{{\rm{\Sigma }}}_{{st}}^{(11)}}{55{{\rm{\Sigma }}}_{{st}}^{(11)}-20{{\rm{\Sigma }}}_{{st}}^{(12)}+4{{\rm{\Sigma }}}_{{st}}^{(13)}+8{{\rm{\Sigma }}}_{{st}}^{(22)}}.\end{eqnarray} \tag{ 111 }$$

For reference, we also include a direct fit of Daligault & Murillo (2005) to the MD data of Ranganathan et al. (2003), given by

$$\begin{eqnarray}&&{D}^{*}=0.0028+0.00525{\left(\displaystyle \frac{173}{{\rm{\Gamma }}}-1\right)}^{1.154}.\end{eqnarray} \tag{ 112 }$$

The agreement between the HNC and MD simulations shows that the HNC does a better job of accounting for correlation physics in strongly coupled plasmas than a simple screened Coulomb potential and allows for a surprising (and still physically unexplained) extension of the Burgers formalism into the strongly coupled regime. This recent work allows us to go into the large Γ limit with the Burgers formalism, but the question remains as to how we obtain diffusion coefficients in a reliable manner.

Zoom In Zoom Out Reset image size

The self-diffusion coefficients from the two options in MESA are shown in Figure 44 and correlate with the MD data better than expected for the high Γ regime. In particular, the agreement is much better than that shown in Figure 2 of Baalrud & Daligault (2013) for either "cutoff" or "screened" Coulomb methods. The reason for this agreement is that both MESA implementations use the inter-ion spacing rather than the Debye length once ${\rm{\Gamma }}\gt 1/3$, which yields favorable scalings in the high Γ limit. Iben & MacDonald (1985) constructed their fitting formula based on a few numerical results for ${\rm{\Gamma }}\;\gt$ 1. Paquette et al. (1986) also showed that their formalism can be extended to ${\rm{\Gamma }}\gt 1$ as long as the inter-ion spacing is used rather than the Debye radius for the screening length.

Though MESA does not yet provide the capability of implementing resistance coefficients based on the HNC method, we hope to accomplish this in the near future by means of a table similar to that provided for the coefficients of Paquette et al. (1986). For a more thorough discussion of these methods and the likely path of application to mixtures, consult Beznogov & Yakovlev (2014). We will also need to correctly account for the electron degeneracy and the non-ideal equation of state for the ions, both of which modify the electrostatic field needed to correctly determine the forces that drive diffusion.