MODULES FOR EXPERIMENTS IN STELLAR ASTROPHYSICS (MESA)

The MESA module const contains mathematical, physical, and astronomical constants relevant to stellar astrophysics in cgs units. The primary source for physical constants is the CODATA Recommended Values of the Fundamental Physical Constants (Mohr et al. 2008). Values for the solar age, mass, radius, and luminosity are taken from Bahcall et al. (2005).

The MESA module chem is a collection of data, functions, and subroutines to manage the chemical elements and their isotopes. It contains basic information about the chemical elements and their isotopes from hydrogen through uranium. It includes routines for translating between atomic weights and numbers and isotope names. It contains full listings of solar abundances on several scales (Anders & Grevesse 1989; Grevesse & Noels 1993; Grevesse & Sauval 1998; Lodders 2003; Asplund et al. 2005). Module chem contains a framework for the user to provide an arbitrary set of species in a text file.

The EOS is delivered by the eos module. It works with density, ρ, and temperature, T, as independent variables. These are the natural variables in a Helmholtz free energy formulation of the thermodynamics. However, as some calculations are more naturally performed using pressure, P, and T (as in a Gibbs free energy formulation), a simple root find can provide ρ given the desired P_gas = P − aT⁴/3 and T. While conceptually simple, this can impose a substantial computational overhead if done for each eos call. To alleviate this computational burden, the root finds are pre-processed, creating a set of tables indexed by P_gas and T. As a result, the runtime cost of evaluating eos using P_gas and T is the same as for using ρ and T, as long as the P_gas − T requests are within the pre-computed ranges. When outside those ranges, the root find is performed during runtime, slowing the computations.

The MESA ρ–T tables are based on the 2005 update of the OPAL EOS tables (Rogers & Nayfonov 2002). To extend to lower temperatures and densities, we use the SCVH tables (Saumon et al. 1995) and construct a smooth transition between these tables in the overlapping region that we define (shown by the blue dotted lines in Figure 1). The limited thermodynamic information available from these EOSs restricts our blending to the output quantities listed in Table 3. The resulting MESA tables are more finely gridded than the original tables (so that no information is lost) and are provided at six X and three Z values: X = (0.0, 0.2, 0.4, 0.6, 0.8, 1.0) and Z = (0.0, 0.02, 0.04) in keeping with the OPAL tables, allowing for Helium rich compositions. In order to save space, the MESA tables are not rectangular in the independent variables. Instead, the region occupied by usual stellar models is roughly rectangular in the stellar modeling motivated variables, log T and log Q = log ρ − 2log T + 12. The range in log T is from 2.1 to 8.2 in steps of 0.02 and the range in log Q is from −10.0 to 5.69 in steps of 0.03. Partials with log T and log Q are derived from the interpolating polynomials, while partials with respect to log ρ then follow. The resulting region of these MESA tables is that inside of the dashed black lines of Figure 1. The MESA P_gas–T tables are rectangular in log T and log W = log P_gas − 4log T over a range −17.2 ⩽ log W ⩽ −2.9 and 2.1 ⩽ log T ⩽ 8.2.

Zoom In Zoom Out Reset image size

Outside the region covered by the MESA tables, the HELM (Timmes & Swesty 2000), and PC (Potekhin & Chabrier 2010) EOSs are employed. Both HELM and PC assume complete ionization and were explicitly constructed from a free energy approach, guaranteeing thermodynamic consistency. In nearly all cases, the full ionization assumption is appropriate since the OPAL and SCVH tables are used at those cooler temperatures where partial ionization is significant.¹⁴ Since the MESA tables are only constructed for Z ⩽ 0.04, eos uses HELM and PC for Z>0.04 in the whole ρ–T plane.

HELM was constructed for high temperatures (up to log T = 13) and densities (up to log ρ = 15), and accounts for the onset of electron–positron pair production at high temperatures. The dotted red line in Figure 1 shows where the number of electrons per baryon has doubled due to pair production. The domination of pairs in the plasma creates a region where Γ₁ < 4/3 (to the left of the dashed red line). The blending region to HELM (from the MESA tables) is shown by the black dashed lines in Figure 1 and can be modified by the user. In this transition region, the blend of the two EOSs is performed in a way that preserves thermodynamic consistency. Therefore, if each separate EOS satisfies Maxwell's relations, the blend will also satisfy them. To accomplish this, we linearly sum the EOS quantities Q_i (i.e., P, E, S and their partial derivatives with respect to ρ and T) needed to satisfy Maxwell's relations (Timmes & Swesty 2000).¹⁵ The blend is calculated by defining the boundary limits, inside of which we define a fractional "distance," F, from the boundary. As F varies from 0 to 1, we use the smoothing function S = (1 − cos(Fπ))/2 and for each of the nine quantities we construct Q_i = SQ^A_i + (1 − S)Q^B_i, where Q^A_i and Q^B_i are the outputs from the two EOSs. We then use these to rederive the thermodynamic quantities (χ_ρ, χ_T, C_P, ∇_ad, Γ₃, Γ₁) delivered by the eos routine.

In late stages of the cooling of white dwarfs, the ions in the core will crystallize. For pure oxygen, the crystallization limit corresponds to a value of the Coulomb coupling parameter, Γ ≈ 175, shown by the black dashed line in Figure 1. In this region, we use the PC EOS, which accounts for the modified thermodynamics of a crystal, as well as carefully handling mixtures (e.g., carbon and oxygen). The blend between PC and HELM (as shown by the dot-dashed red lines in Figure 1) is performed in the same manner as described above. In the dense liquid realm, the blending region is defined by the Coulomb coupling parameter, $\Gamma =\overline{\rm Z}^2e^2/a_i k_BT$, where a_i is the mean ion spacing, and $\overline{\rm Z}$ is the average ion charge. The default choice is PC for Γ>80 and HELM for Γ < 40. The PC EOS is not constructed for arbitrarily low densities, forcing a transition to HELM at log ρ < 2.8, with the blend beginning at log ρ = 3.7. These boundaries may be re-defined by the user if needed.

In addition to the two independent variables, the eos module requires as input X, Z, $\overline{{\rm A}}$ (the mass-averaged atomic weight of metals) and $\overline{\rm Z}$ (the mass-averaged atomic charge of metals). When operating in the regime where the PC EOS is implemented, the mass fractions for all isotopes with mass fractions above a specified minimum are needed (default is 0.01), allowing PC to correctly handle isotope mixtures. It returns a total of 16 quantities (listed in Table 3) as well as the partial derivatives of each quantity with respect to the independent variables. The tables are interpolated in the independent variables using bicubic splines from interp\₂d with partial derivatives determined from the splines. Separate quadratic interpolations are performed in X and Z.

The construction of eos tables as outlined above is the default option for MESA but the eos module has the flexibility to accept tables from any source so long as the tables conform to the MESA standard format. For example, the comparison with the Stellar Code Calibration project (Weiss et al. 2007) described in Section 7.1.2 utilizes tables constructed using FreeEOS.¹⁶ The FreeEOS code does not cover the same range of ρ and T as SCVH+OPAL+HELM but the eos module is designed with this flexibility in mind: the table dimensions are specified in the table headers and the module dynamically allocates arrays of the appropriate size to hold them when the tables are read in.

Since not all EOS sources may be in the tabular form desired by eos, we have created a module, other_eos, that provides the user an opportunity to incorporate their own EOS and use it with the stellar evolution module MESAstar.

The pre-processor make_kap resides within the kap module and constructs the MESA opacity tables by combining radiative opacities with the electron conduction opacities from Cassisi et al. (2007). In the rare circumstances where ρ or T are outside the region covered by Cassisi et al. (2007) (−6 ⩽ log ρ ⩽ 9.75 and 3 ⩽ log T ⩽ 9), the Iben (1975) fit to the Hubbard & Lampe (1969) electron conduction opacity is used for non-degenerate cases, while the Yakovlev & Urpin (1980) fits are used for degenerate cases. Radiative opacities are taken from Ferguson et al. (2005) for 2.7 ⩽ log T ⩽ 4.5 and OPAL (Iglesias & Rogers 1993, 1996) for 3.75 ⩽ log T ⩽ 8.7. The low T opacities of Ferguson et al. (2005) include the effects of molecules and grains on the radiative opacity. Tables from OP (Seaton et al. 2005) can be used in place of OPAL as the table format is identical. The radiative opacity is dominated by Compton scattering for log T>8.7 and is calculated using the equations of Buchler & Yueh (1976) up to a density of 10⁶ g cm^-3. We use the HELM EOS to calculate the increasing number of electrons and positrons per baryon when pair production becomes prevalent, an important opacity enhancement.

The OPAL tables with fixed metal distributions are called Type 1 (Iglesias & Rogers 1993, 1996) and cover the region 0.0 ⩽ X ⩽ 1 − Z and 0.0 ⩽ Z ⩽ 0.1. Additionally, there is support for the OPAL Type 2 (Iglesias & Rogers 1996) tables that allow for varying amounts of C and O beyond that accounted for by Z; these are needed during helium burning and beyond. These have a range 0.0 ⩽ X ⩽ 0.7 and 0.0 ⩽ Z ⩽ 0.1.

The resulting kap tables cover the large range 2.7 ⩽ log T ⩽ 10.3 and −8 ⩽ log R ⩽ 8 (R = ρ/T³₆, so log R = log ρ − 3log T + 18), as shown by the heavy orange and black lines in Figure 2. The MESA release includes MESA opacity tables derived from Type 1 and 2 OPAL tables, tables from OP, and Ferguson et al. (2005). The heavy orange lines delineate the boundaries where we use existing tables to make the MESA opacity table. The blended regions in Figure 2 are where two distinct sources of radiative opacities exist for the same parameters, requiring a smoothing function that blends them in a manner adequate for derivatives. The blend is calculated at a fixed log R by defining the upper (log T_U) and lower (log T_L) boundaries of the blending region in log T space, where κ_U (κ_L) is the opacity source above (below) the blend. We perform the interpolation by defining F = (log T − log T_L)/(log T_U − log T_L), and using a smooth function S = (1 − cos(Fπ))/2 for

$$\begin{equation} \log \kappa =S \log \kappa _U(R,T)+(1-S)\log \kappa _L(R,T). \end{equation} \tag{ 1 }$$

At high temperatures, the blend from Compton to OPAL (or OP) has log T_U = 8.7 and log T_L = 8.2. At low temperatures, the blend between Ferguson et al. (2005) and OPAL has log T_U = 4.5 and log T_L = 3.75.

Zoom In Zoom Out Reset image size

The absence of tabulated radiative opacities for log R>1 and log T < 8.2 (the region below the heavy dashed line in Figure 2) leads us to use the radiative opacity at log R = 1 (for a specific log T) when combining with the electron conduction opacities. This introduces errors in the MESA opacity table between log R = 1 and the region to the right of the dashed blue line in Figure 2 where conductive opacities become dominant. However, as we show in Figure 3, main-sequence stars are always efficiently convective in this region of parameter space, alleviating the issue.

Zoom In Zoom Out Reset image size

The module kap gives the user the resulting opacities by interpolating in log T and log R with bicubic splines from interp\₂d. The user has the option of either linear or cubic interpolation in X and Z and can specify whether to use the fixed metal (Type 1) tables or the varying C and O (Type 2) tables. In the latter case, the user must specify the reference C and O mass fractions, usually corresponding to the C and O in the initial composition.

For requests outside the log T and log R boundaries, the following is done. The region to the left of log R = −8 and below log T = 8.7 is electron scattering dominated, so the cross-section per electron is density independent. However, the increasing importance of the Compton effect as the temperature increases (which is incorporated in the OPAL/OP tabulated opacities) must be included, so we use the opacity from the table at log R = −8 at the appropriate value of log T. For higher temperatures (log T>8.7) electron–positron pairs become prevalent, as exhibited by the red dashed line that shows where the number of positrons and electrons from pair production exceeds the number of electrons from ionization. MESA incorporates the enhancement to the opacity from these increasing numbers of leptons per baryon.

At the end of a star's life, low enough entropies can be reached that an opacity for log R>8 is needed. When kap is called in this region, we simply use the value at log R = 8 for the same log T. For regions where Z>0.1, the table at Z = 0.1 is used.

The resulting opacities for Z = 0.019 and Y = 0.275 are shown in Figure 3, both as a color code, and as contours relative to the electron scattering opacity, κ₀ = 0.2(1 + X) cm² g^-1. The orange lines show (top to bottom) where log R = −8, log R = 1 and log R = 8. We show a few stellar profiles for main-sequence stars as marked. The green parts of the line are where heat transfer is dominated by heat transport, requiring an opacity, whereas the light blue parts of the line are where the model is convective. As is evident, nearly all of the stellar cases of interest (shown by the green-blue lines) are safely within the boundaries or the MESA tables. The lack of radiative opacities in the higher density region to the right of log R = 1 implies opacity uncertainties until the dark blue line is reached (where the conductive opacity takes over). However, the stellar models are convectively efficient in this region, so that the poor value for κ does not impact the result as long as the convective zone's existence is independent of the opacity (the typical case for these stars, where the ionization zone causes the convection).

It is also possible to generate a new set of kap readable opacity tables using the make_kap pre-processor. The requirements are high-temperature radiative opacities in the standard OPAL format and low-temperature radiative opacities in the number and format provided by Ferguson et al. (2005).¹⁷ Specific high-temperature radiative opacities can be made by using the OPAL site¹⁸ or the Opacity Project site.¹⁹

Since not all opacity sources can be placed in the tabular form desired by kap, we have created a module, other_kap, that provides the user an opportunity to incorporate their own opacity source. A simple flag tells MESAstar to call other_kap rather than kap, allowing for experiments with new opacity schemes and physics updates. The first example of such an implementation that has now become a MESA module is karo. It was developed to study the stellar evolution effects of dust-driven winds in carbon-rich stars, using the Rosseland opacities of Lederer & Aringer (2009) and the hydro-dynamical wind models of Mattsson et al. (2010).

The rates module contains thermonuclear reaction rates from Caughlan & Fowler (1988, CF88) and Angulo et al. (1999, NACRE), with preference given to the NACRE rate when available. The reaction rate library includes more than 300 rates for elements up to nickel, and includes the weak reactions needed for hydrogen burning (e.g., positron emissions, electron captures), as well as neutron–proton conversions and a few other electron and neutron capture reactions. Significant updates to the NACRE rates have been included for ¹⁴N(p,γ)¹⁵O (Imbriani et al. 2004), triple-α (Fynbo et al. 2005), ¹⁴N(α, γ)¹⁸F (Görres et al. 2000), and ¹²C(α, γ)¹⁶O (Kunz et al. 2002). In these special cases, the rate can be selected from CF88, NACRE, or the newer reference by the user at run time.

The weaklib module calculates lepton captures and β-decay rates for the high densities and temperatures encountered in late stages of stellar evolution. The rates are based on the tabulations of Fuller et al. (1985), Oda et al. (1994), and Langanke & Martínez-Pinedo (2000) for isotopes with 45 < A < 65. The most recent tabulations of Langanke & Martínez-Pinedo (2000) take precedence, followed by Oda et al. (1994), then Fuller et al. (1985). The user can override this to create tables using any combination of these or other sources.

The screen module calculates electron screening factors for thermonuclear reactions in both the weak and strong regimes. The treatment has two options. One is based on Dewitt et al. (1973) and Graboske et al. (1973). The other²⁰ combines Graboske et al. (1973) in the weak regime and Alastuey & Jancovici (1978) with plasma parameters from Itoh et al. (1979) in the strong regime.

The neu module calculates energy-loss rates and their derivatives from neutrinos generated by a range of processes including plasmon decay, pair annihilation, bremsstrahlung, recombination, and photo-neutrinos (i.e., neutrino pair production in Compton scattering). It is based on the publicly available routine (see footnote ²⁰, derived from the fitting formulae of Itoh et al. (1996).

The net module implements nuclear reaction networks and is derived from publicly available code²⁰. It includes a "basic" network of eight isotopes: ¹H, ³He, ⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, and ²⁴Mg, and extended networks for more detailed calculations including coverage of hot CNO reactions, α-capture chains, (α,p)+(p,γ) reactions, and heavy-ion reactions (Timmes 1999). In addition to using existing networks, the user can create a new network by listing the desired isotopes and reactions in a data file that is read at run time. The amount of heat deposited in the plasma by reactions is derived from the nuclear masses in chem, taken from the JINA Reaclib database (Rauscher & Thielemann 2000; Sakharuk et al. 2006; Cyburt et al. 2010), and accounts for positron annihilations and energy lost to weak neutrinos, using Bahcall (1997, 2002) for the hydrogen burning reactions. The list of approximately 350 reactions is stored in a data file that catalogs the reaction name, the input and output species, and their heat release.

The jina module is an alternative nuclear network module that specializes in large networks. It is based on the "netjina" package by Ed Brown and uses the JINA Reaclib database for thermonuclear reaction rates (Rauscher & Thielemann 2000; Sakharuk et al. 2006),²¹ the rates and weaklib modules for weak interactions, and screen for electron screening. Most importantly, it allows the user to create large nuclear networks by specifying the list of isotopes to consider. All nuclear reactions (both strong and weak) linking the isotopes in the set are automatically included in the network. In all, jina covers more than 76,000 nuclear reactions involving more than 4500 isotopes. The jina module is slower than net for small networks but the flexibility and capacity to handle large networks make it advantageous in some cases.

Both net and jina include one-zone burn routines that operate on a user-defined initial composition, nuclear network, and a trajectory comprising density and temperature as a function of time. The one-zone burn routines interface with mtx and num, enabling the use of the sparse matrix solver, which substantially improves performance compared to the dense matrix solver for networks of more than a few hundred isotopes. Figures 4 and 5 demonstrate these one-zone routines operating on conditions appropriate for the Sun on the main sequence and for a core He-burning star, respectively. Both examples were evolved at fixed density and temperature for 10 Gyr. Each figure compares three networks of varying size in terms of the mass fractions and net energy generation per unit mass, e_nuc (erg g⁻¹), produced. Tables 4 and 5 complement Figures 4 and 5, respectively, by listing the final values from each of the one-zone burn simulations. These comparisons indicate that the eight isotope network produces results that agree with larger networks to 4–5 significant figures in net energy generation per unit mass and generally 2–3 significant figures in the mass fractions of various isotopes. Hence, the eight isotope network is sufficiently accurate to describe the energy generation for hydrogen and helium burning.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mlt module implements the standard mixing length theory (MLT) of convection as presented by Cox & Giuli (1968, chap. 14). There are options both for computing the actual temperature gradient, ∇_T, when the total luminosity, L, is specified and for computing the convective luminosity, L_conv, when ∇_T is specified. The mlt module calculates diffusion coefficients for those codes, such as MESAstar, that treat convective mixing of elements as a diffusive process. The quantities listed in Table 6 and their partial derivatives with respect to several physical variables are returned by the mlt module.

In addition to the standard MLT of Cox and Giuli, the mlt module includes the option to use the modified MLT of Henyey et al. (1965). Whereas the standard MLT assumes high optical depths and no radiative losses, the Henyey et al. (1965) variation allows the convective efficiency to vary with the opaqueness of the convective element, an important effect for convective zones near the outer layers of stars. If the Henyey et al. (1965) option is used, the parameters ν (a mixing length velocity multiplier) and y (a parameter that sets the temperature gradient in a rising bubble) may be set by the user. They default to the recommended values of y = 1/3 and ν = 8.

Toward the center of a star, the commonly used definition of the pressure scale height, λ_P = P/gρ, diverges as g → 0. Therefore, we provide the option of using the alternate definition of Eggleton (1971), λ'_P = (P/Gρ²)^1/2, when λ'_P < λ_P. At the center of the star, λ'_P ∼ R.

The dependence of mixing instabilities on composition gradients (e.g., semiconvection and thermohaline mixing) are presently not taken into account by the mlt module. To include this capability is a high priority for the MESA project.

As described in Section 6.2, MESAstar treats convective mixing as a time-dependent, diffusive process with a diffusion coefficient, D, determined by the mlt module. In the absence of a three-dimensional hydrodynamical treatment of convection it is necessary to account for the hydrodynamical mixing instabilities at convective boundaries, termed overshoot mixing, via a parametric model. After the MLT calculations have been performed, MESAstar sets the overshoot mixing diffusion coefficient

$$\begin{equation} D_{\rm OV} = D_{\rm conv,0} \exp \left(-{2z\over f\lambda _{P,0}}\right), \end{equation} \tag{ 2 }$$

where D_conv,0 is the MLT derived diffusion coefficient at a user-defined location near the Schwarzschild boundary, λ_P,0 is the pressure scale height at that location, z is the distance in the radiative layer away from that location, and f is an adjustable parameter (Herwig 2000). In MESAstar the adjustable parameter, f, may have different values at the upper and lower convective boundaries for non-burning, H-burning, He-burning, and metal-burning convection zones.

Parameters are provided to allow the user to set a lower limit on D_OV below which overshoot mixing is neglected and to limit the region of the star over which overshoot mixing will be considered. So as to model the ¹³C pocket needed for s-process nucleosynthesis, MESAstar also allows an increase in the overshooting parameter at the bottom of the convective envelope during the third dredge-up compared to the inter-pulse value (Lugaro et al. 2003). There is also an option to change the value of overshoot mixing at the bottom of the asymptotic giant branch (AGB) thermal pulse-driven convection zone compared to the standard value chosen for the bottom of the He-burning convection zone.

As described in Section 6.2, the pressure, P_s, and temperature, T_s, at the top of the outermost cell in MESAstar must be set by an atmospheric model. This is done by the atm module, which uses M, R, and L to provide P_s and T_s. It also gives partial derivatives of T_s and P_s with respect to the input variables. The atm module assumes the plane parallel limit, so that the relevant variables are g = GM/R² and T_eff⁴ = L/4πσ_SBR². With some options, the user must specify the optical depth τ_s to the base of the atmosphere, whereas in other cases, the atm module has an implicit value. Three methods are supplied by atm: direct integrations, interpolations in model atmosphere tables, and a "simple" recipe.

The integrations of the hydrostatic balance equation, dP_gas/dτ = g/κ − (a/3)dT⁴/dτ, with dτ = −κρdr are performed using either the relation T⁴(τ) = 3T⁴_eff(τ + 2/3)/4 (Eddington 1926), or the specific T–τ relation of Krishna Swamy (1966). These integrations start at τ = 10⁻⁵ and end at a user specified stopping point, τ_s, which defaults to τ_s = 2/3 (0.312) for Eddington (Krishna Swamy).²² The routine integrates the gas pressure and then adds the radiation pressure at the stopping point to get P_s.

The MESA model atmosphere tables come in two forms. The MESA photospheric tables (which return T_s ≡ T_eff and assume that τ_s ≈ 1) cover log Z/Z_☉ = −4 to +0.5 assuming the Grevesse & Noels (1993) Solar abundance mixture. They span log(g) = −0.5 to 5.5 at 0.5 dex intervals and T_eff= 2000–50,000 K at 250 K intervals. They are constructed, in precedence order, with, first, the PHOENIX (Hauschildt et al. 1999a, 1999b) model atmospheres (which span log(g) = −0.5 to 5.5 and T_eff = 2000 to 10,000 K); and second, the Castelli & Kurucz (2003) model atmospheres (which span log(g) = 0–5 and T_eff = 3500–50,000 K). In regions where neither table is available, we generate the MESA table entry using the integrations described above with the Eddington T–τ relation. The second MESA table is for solar metallicity and gives P_s and T_s at τ_s = 100. It is primarily for the evolution of low-mass stars, brown dwarfs, and giant planets. It is constructed from Castelli & Kurucz (2003) and, for T_eff < 3000 K, the COND model atmospheres (Allard et al. 2001) which assume gravitational settling of those elements that form dust, depleting those elements from the photosphere. Figure 6 shows the regions where the different sources are used, and in those regions where there are no published results, we use the integration of the Eddington T–τ relation.

Zoom In Zoom Out Reset image size

Finally, there is a simple option where the user specifies τ_s and we use the constant opacity, κ_s, solution of radiative diffusion,

$$\begin{equation} P_s = \frac{\tau _s g}{\kappa _s} \left[1 + 1.6\times 10^{-4} \kappa _s\left(\frac{L/{L_{\odot }}}{M/\,{M_{\odot }}}\right)\right], \end{equation} \tag{ 3 }$$

where the factor in square brackets accounts for the nonzero radiation pressure (see, e.g., Cox & Giuli 1968, Section 20.1). The temperature is simply given by the Eddington relation. The user can either specify κ_s or it will be calculated in an iterative manner using the initial value of P_s from an initial guess at κ_s (usually given by MESAstar as the value in the outermost cell; see Section 6.2). In addition, the atm module has the option to revert to Equation (3) if a model wanders outside the range of the currently used model atmosphere tables or if the atmosphere integration fails for any reason.

MESA diffusion calculates particle diffusion and gravitational settling by solving Burger's equations using the method and diffusion coefficients of Thoul et al. (1994). The transport of material is computed using the semi-implicit, finite difference scheme described by Iben & MacDonald (1985). Radiative levitation is not presently included. The diffusion module treats the elements present in the stellar model as belonging to "classes" defined by the user in terms of ranges of atomic masses. For each class, the user specifies a representative isotope, and all members of that class are treated identically with their diffusion velocities determined by the representative isotope, and the diffusion equation solved with the mass fraction in that class. The caller can either specify the ionic charge for each class at each cell in the model or has the charge calculated by the ionization module, which estimates the typical ionic charge as a function of T, ρ, and free electrons per nucleon from Paquette et al. (1986).

The lines in Figure 7 plot four classes (H, He, O, and Fe) with a solar model from MESAstar and compares where possible to the results from Figure 9 of Thoul et al. (1994), shown by the filled green circles. The agreement is excellent for H, O, and Fe (when we fix Fe to have the Z = 21 ionization state chosen by Thoul et al. 1994). Thoul et al. (1994) did not exhibit the He velocity, so we have no comparison. For Fe, we also show the diffusion velocity when ionization finds a changing ionization state in the Z = 16, 17, 18 region (shown by the upper dot-dashed blue line), highlighting the need to better determine the Fe ionization state (Gorshkov & Baturin 2008). We also compared the diffusion output to the recent calculations of Gorshkov & Baturin (2008), finding agreement at better than 5% for the Fe case at Z = 26 and for O.

Zoom In Zoom Out Reset image size

The diffusion calculation can be restricted to areas where the temperature is above some minimum value, or where the mass fraction of a diffusing element is above some minimum value, aiding the convergence of solutions in a variety of environments. The physics implementation is presently limited to regions where the Coulomb coupling parameter, Γ, is less than unity. At present, this inhibits an accurate calculation for segregation and settling of the remaining envelope H and He envelope on a cooling white dwarf.

The complex, nonlinear behavior of stellar structure and evolution models makes it difficult to disentangle the effects of model components (e.g., EOS, opacities, boundary conditions, etc.) when comparing results of separate codes. By design, the modularity of MESA allows individual physics modules to be incorporated into an existing stellar evolution code, tested, and then compared against the prior implementation of comparable physics in the same code.

During the development of MESA, several MESA modules were integrated into the Dartmouth Stellar Evolution Program (DSEP; Dotter et al. 2007). This section reports the results of using four MESA modules, eos, kap, atm, and mlt, in DSEP to compute the evolution of a 1.0 M_☉ star with initial values of X = 0.70 and Z = 0.02. The star was evolved from the fully convective pre-main sequence (PMS) to the onset of the core He flash. This was done six times: once, as the control case, using only DSEP routines and no MESA modules; next, using each of four MESA modules individually; and, finally, using the four MESA modules at the same time in DSEP.

DSEP employs a ρ(P, T) EOS and so the MESA P_gas–T tables were used during the eos test. Though DSEP and kap use the same sources for radiative opacities, they differ in interpolation methods and the treatment of electron conduction opacities (see Bjork & Chaboyer 2006, for a thorough list of the physics in DSEP). When atm was tested, we used the Eddington gray atmosphere model integrated to τ = 2/3. DSEP uses the Henyey et al. (1965) modification of the MLT, which is available in mlt, and assumes that convective regions are instantaneously mixed to a uniform composition.

DSEP tracks employing either the atm or the mlt modules produce results that agree with the DSEP-only track to about one part in 10⁴. DSEP tracks employing the kap and eos modules exhibit some difference when compared to the DSEP-only track but, even in these cases, the main-sequence lifetime differs by less than 0.3% and T_eff differs by less than 10 K along the main sequence. As shown in Figure 8, the largest discrepancy between the DSEP-only track and the one that employs all four MESA modules appears in the T_c–ρ_c diagram when ρ_c>3 × 10⁴ g cm^-3, corresponding to the growing helium core in the center of the red giant. Above log ρ_c = 4, the track employing MESA modules is hotter than the DSEP-only track by ∼0.02 in log T_c at constant log ρ_c. The center of the model has entered the region of electron degeneracy and electron conduction has become an important source of opacity. The majority of the difference is due to the EOS whereas the opacity difference amounts to about −0.005 in log T_c, in the opposite direction to the EOS. The hotter conditions produced by the eos module are likely the cause for the slightly shorter red giant branch (RGB) lifetime that can be seen in Figure 8.

Zoom In Zoom Out Reset image size

MESAstar receives basic input from two Fortran namelist files. One file specifies the type of evolutionary calculation to be performed, the type of input model to use, the source of EOS and opacity data, the chemical composition and nuclear network, and other properties of the input model. The second file specifies the controls and options to be applied during the evolution.

There are two ways to start a new evolutionary sequence with MESAstar. The first is to use a saved model from a previous run. A variety of saved models are distributed with MESA as a convenience. These saved models fall into three general categories: (1) zero age main sequence (ZAMS) models for Z = 0.02 with 32 masses between 0.08 and 100 M_☉ (MESAstar will automatically interpolate any mass within this range); (2) very low mass, PMS models for Z = 0.02 and masses from 0.001 to 0.025 M_☉; and (3) white dwarf models for Z = 0.02 with He cores of 0.15–0.45 M_☉, C/O cores of 0.496–1.025 M_☉, and O/Ne cores of 1.259–1.376 M_☉. The user can also create saved models for essentially any purpose through available controls.

The second way to start a new evolution is to create a PMS model by specifying the mass, M, a uniform composition, a luminosity, and a central temperature, T_c low enough that nuclear burning is inconsequential (T_c = 9 × 10⁵ K by default). For a fixed T_c and composition, the total mass depends only on the central density, ρ_c. An initial guess for ρ_c is made by using the n = 1.5 polytrope, which is appropriate for a fully convective star, but we do not assume that the star is fully convective during the subsequent search for a converged PMS model. Instead, MESAstar uses the mlt, eos, and Newton solver from num to search for a ρ_c that gives a model of the desired mass. The PMS routine presently creates starting models for 0.02 ⩽ M/ M_☉ ⩽ 50. Beyond these limits we find challenges converging the generated PMS model within the MESAstar evolutionary loop. For such cases, it is currently better to generate a starting model within the acceptable mass range, save it, relax it to a new mass with a specified mass gain or loss (see Section 6.6), and save that model.

MESAstar has the ability to create a binary file of its complete current state, called a photo, at user-specified timestep intervals. Restarting from a photo ensures no differences in the ensuing evolution. When restarting from a photo, many controls and options can be changed. A photo is different than a saved model in that a saved model is a text file containing a minimal description of the structure and composition but does not have enough information to allow a perfect restart. However, saved models are not tied to a particular version of the code and therefore are suitable for long-term use or sharing with other users.

There are two additional types of output files, logs and profiles. A log records evolutionary properties over time such as stellar age, current mass, and a wide array of other quantities. A profile records model properties at a specified timestep at each zone from surface to center. MESAstar can also output models in the FGONG format²³ for use with stellar pulsation codes and se output for nucleosynthesis post-processing with NuGrid codes.²⁴ Finally, a few simple lines of user-supplied code allows for saving variables or combinations of variables that are not in the list of predefined options.

MESAstar builds one-dimensional, spherically symmetric models by dividing the structure into cells, anywhere from hundreds to thousands depending on the complexity of nuclear burning, gradients of state variables, composition, and various tolerances. Cells are numbered starting with one at the surface and increasing inward. MESAstar does not require the structure equations to be solved separately from the composition equations (operator splitting). Instead, it simultaneously solves the full set of coupled equations for all cells from the surface to the center. The solution of the equations is done by the Newton solver from num using either banded or sparse matrix routines from mtx. The partial derivatives for use by the solver are calculated analytically using the partials returned by modules such as eos, kap, and net.

Each cell has some variables that are mass-averaged and others that are defined at the outer face, as shown in Figure 9. This way of defining the variables is a consequence of the finite volume, flux conservation formulation of the equations, and improves stability and efficiency (Sugimoto et al. 1981). The inner boundary of the innermost cell is usually the center of the star and, therefore, has radius, luminosity, and velocity equal to zero. Nonzero center values can be used for applications that remove the underlying star (e.g., the envelope of a neutron star), in which case the user must define the values of M_c and L_c at the inner radius R_c. The cell mass-averaged variables are density ρ_k, temperature T_k, and mass fraction vector X_i,k. The boundary variables are mass interior to the face m_k, radius r_k, luminosity L_k, and velocity v_k. In addition to these basic variables, composite variables are calculated for every cell and face, such as _nuc, κ, σ_k, and F_k (see Table 1 for variable definitions). All variables are evaluated at time t + δt unless otherwise specified.

Zoom In Zoom Out Reset image size

The density evolution of cell k is determined by a finite volume form of the mass conservation equation

$$\begin{equation} \rho _k = \frac{dm_k}{(4/3) \pi \big(r_k^3 - r_{k+1}^3\big)} \enskip. \end{equation} \tag{ 4 }$$

For the innermost cell, r_k+1 is replaced by the inner boundary condition which is typically zero but can be nonzero for some applications. We reformulate many of our equations to improve numerical stability of the linear algebra and minimize round-off errors. We thus rewrite Equation (4) as

$$\begin{equation} \ln r_k = \frac{1}{3} \ln \left[ r_{k+1}^3 + \frac{3}{4\pi }\frac{dm_k}{\rho _k} \right] \enskip. \end{equation} \tag{ 5 }$$

The velocity of face k is zero unless the hydrodynamics option is activated, in which case

$$\begin{equation} v_k = r_k \frac{d(\ln r_k)}{dt} \end{equation} \tag{ 6 }$$

is the Lagrangian time derivative of the radius at face k. For enhanced numerical stability, we rescale this equation by dividing by the local sound speed.

The pressure P_k is set by momentum conservation at interior cell boundaries,

$$\begin{eqnarray} P_{k-1} - P_k &=& {\overline{dm}_k} \left[ \left(\frac{{\it dP}}{dm}\right)_{\rm hydrostatic} + \left(\frac{{\it dP}}{dm}\right)_{\rm hydrodynamic} \right] \nonumber \\ &=& {\overline{dm}_k} \left[ -\frac{Gm_k}{4 \pi r_k^4} - \frac{a_k}{4 \pi r_k^2} \right] \enskip, \end{eqnarray} \tag{ 7 }$$

where $\overline{dm}_k = 0.5 (dm_{k-1} + dm_k)$, and a_k is the Lagrangian acceleration at face k, evaluated by the change in v_k over the timestep δt. The acceleration is set to zero if the hydrodynamic option is not used. Similarly, the temperature of interior cells T_k is set by energy transport across interior cell boundaries,

$$\begin{equation} T_{k-1} - T_k = {\overline{dm}_k} \left[ \nabla _{T,k} \left(\frac{{\it dP}}{dm}\right)_{\rm hydrostatic} \frac{\overline{T}_k}{\overline{P}_k} \right] \enskip, \end{equation} \tag{ 8 }$$

where ∇_T,k = dln T/dln P at face k from the MESA module mlt (see Section 5.1), $\overline{T}_k = (T_{k-1} dm_k + T_k dm_{k-1})/(dm_k + dm_{k-1})$ is the temperature interpolated by mass at face k, and $\overline{P}_k = (P_{k-1} dm_k + P_k dm_{k-1})/(dm_k + dm_{k-1})$ is the pressure interpolated by mass at face k. For enhanced numerical stability, we rescale Equation (7) by dividing by $\overline{P}_k$ and Equation (8) by dividing by $\overline{T}_k$.

The pressure and temperature boundary conditions are constructed by using P_s and T_s from the MESA module atm (see Section 5.3). The difference in pressure and temperature from the surface to the center of the first cell is found from hydrostatic equilibrium and ∇_T by

$$\begin{eqnarray} {\it dP}_s &=& \frac{G m_1 dm_1/2}{4 \pi r_1^4}, \nonumber \\ dT_s &=& {\it dP}_s \nabla _{T,1} \frac{T_1}{P_1} \enskip. \end{eqnarray} \tag{ 9 }$$

The boundary conditions are then

$$\begin{eqnarray} \ln T_1 &=& \ln (T_s+ {\it dT}_s), \nonumber \\ \ln P_1 &=& \ln (P_s+ {\it dP}_s) \enskip. \end{eqnarray} \tag{ 10 }$$

These implicit equations for P₁ and T₁ are solved together with the regular structure and composition equations.

Our finite volume form of energy conservation for cell k is

$$\begin{equation} L_k - L_{k+1} = dm_k(\epsilon _{\rm nuc} - \epsilon _{\nu,{\rm thermal}} + \epsilon _{\rm grav}) \enskip, \end{equation} \tag{ 11 }$$

where _nuc (from module net or jina) is the total nuclear reaction specific energy generation rate minus the nuclear reaction neutrino-loss rate, and _ν,thermal (from module neu) is the specific thermal neutrino-loss rate. The _grav term is the specific rate of change of gravitational energy due to contraction or expansion,

$$\begin{equation} \epsilon _{\rm grav} = -T \frac{ds}{dt} = - T C_{\rm P} \left[(1 - \nabla _{\rm ad} \chi _T) \frac{d\ln T}{dt} - \nabla _{\rm ad} \chi _{\rho } \frac{d\ln \rho }{dt} \right] \!, \end{equation} \tag{ 12 }$$

where dln T/dt and dln ρ/dt are Lagrangian time derivatives at cell center by mass, and the other symbols are defined in Tables 3 and 6. For the innermost cell, L_k+1 is replaced by the inner boundary condition which is typically zero but can be nonzero, L_c, in specific applications. For additional numerical stability, we rescale Equation (11) by dividing by a scale factor that is typically the surface luminosity of the previous model.

The equation for mass fraction X_i,k of species i in cell k is

$$\begin{eqnarray} X_{i,k}(t + \delta t) - X_{i,k}(t) &=& dX_{\rm burn} + dX_{\rm mix} \nonumber \\ &=& \frac{dX_{i,k}}{dt} \delta t + (F_{i,k+1} - F_{i,k}) \frac{\delta t}{dm_k},\quad\ \end{eqnarray} \tag{ 13 }$$

where dX_i,k/dt is the rate of change from nuclear reactions reported by net or jina, F_i,k is the mass of species i flowing across face k

$$\begin{equation} F_{i,k} = (X_{i,k} - X_{i,k-1}) \frac{\sigma _k}{\overline{dm}_k} \enskip, \end{equation} \tag{ 14 }$$

where σ_k is the Lagrangian diffusion coefficient from the combined effects of convection (Section 5.1) and overshoot mixing (Section 5.2). For numerical stability, σ_k is calculated at the beginning of the timestep and held constant during the implicit solver iterations. This assumption accommodates the non-local overshooting algorithm and significantly improves the numerical convergence. It leads to a small inconsistency between the mixing boundary and the convection boundary as calculated at the end of the timestep.

Equations (5), (7), (8), (11), (14), and, optionally Equation (6), are by default solved fully coupled and simultaneously with a first-order backward differencing time integration.

The generalized Newton–Raphson scheme is represented by

$$\begin{equation} 0 = \vec{F}(\vec{y}) = \vec{F}(\vec{y}_i + \delta \vec{y}_i) = \vec{F}(\vec{y}_i) + \left[ \frac{{\it d\vec{F}}}{d\vec{y}} \right]_i \delta \vec{y}_i + O(\delta \vec{y}_i^{\,2}), \end{equation} \tag{ 15 }$$

where y_i is a trial solution, $\vec{F}(\vec{y}_i)$ is the residual, $\delta \vec{y}_i$ is the correction, and $[d\vec{F}/d\vec{y}]_i$ is the Jacobian matrix.

MESAstar uses the previous model, modified by remeshing, mass change, and element diffusion, as the initial trial solution for the Newton–Raphson solver. This is generally successful because we use analytic Jacobians and have sophisticated timestep controls (see Section 6.4). The use of analytic Jacobians in MESAstar requires that each of the MESA modules provides not just the required output quantities but also quality, preferentially analytic, partial derivatives with respect to the input quantities. At each timestep, MESAstar converges on a final solution by iteratively improving upon the trial solution. We calculate the residuals, construct a Jacobian matrix, and solve the resulting system of linear equations with the solvers in mtx to find the corrections to the variables.

The trial solution is accepted when the corrections and residuals meet a specifiable set of comprehensive convergence criteria. In most cases, the solver is able to satisfy these limits in two or three iterations. However, under difficult circumstances like the He core flash or advanced nuclear burning in massive stars, MESAstar can automatically adjust the convergence criteria. The corrections to the variables will, generally, not produce zero residuals because the system of equations is nonlinear. In some cases, the corrections might make the residuals larger rather than smaller. In such cases, the length of the correction vector is reduced by a line search scheme²⁵ until they improve the residuals. In principle, the residuals can be made arbitrarily small, but this may take a prohibitively large number of iterations. In practice, the use of the line search scheme helps the convergence rate in many cases, but cannot ensure convergence in all cases.

If convergence cannot be achieved with the current timestep, then MESAstar will first try again with a reduced timestep (a "retry") anticipating that a smaller timestep will reduce the nonlinearity. If the retry fails, MESAstar will return to the previous model and with a smaller timestep than it used to get to the current model (a "backup"). If the backup fails, MESAstar will continue to reduce the timestep until either the model converges or the timestep reaches some pre-defined minimum, in which case the evolutionary sequence is terminated.

Timestep selection is a crucial part of stellar evolution. The timesteps should be small enough to allow convergence in relatively few iterations, but large enough to allow efficient evolutions. Changes to the timestep should also provide for rapid responses to varying structure or composition conditions, but need to be carefully controlled to avoid overcorrections that can reduce the convergence rate.

MESAstar does timestep selection as a two stage process. The first stage proposes a new timestep using a scheme based on digital control theory (Soderlind & Wang 2006). The second stage implements a wide range of tests that can reduce the proposed timestep if certain selected properties of the model are changing faster than specified. For the first stage, we use a low-pass filter. The control variable v_c is the unweighted average over all cells of the relative changes in ln ρ, ln T, and ln R. The target value v_t is 10⁻⁴ by default. For improved stability and response, the low-pass filter method uses the previous two results. Let δt_i−1, δt_i, and δt_i+1 be the previous, current, and next timestep, respectively, while v_c,i−1 and v_c,i are the previous and current values of v_c. The timestep for model i + 1 is then determined by

$$\begin{equation} \delta t_{i+1} = \delta t_i f \left[ \frac{f(v_{t}/v_{c,i}) f(v_{t}/v_{c,i-1})}{f(dt_i/dt_{i-1})} \right]^{1/4}, \end{equation} \tag{ 16 }$$

where f(x) = 1 + 2tan⁻¹[0.5(x − 1)]. The control scheme implemented by Equation (16) allows rapid changes in timestep without undesirable fluctuations.

The timestep proposed by this low-pass filtering scheme can be reduced according to a variety of special tests that have hard and soft limits. If a change exceeds its specified hard limit, the current solution is rejected, and the code is forced to do a retry or a backup. If a change exceeds its specified soft limit, the next timestep is reduced proportionally. Examples of special tests include limits on the maximum absolute or relative changes in mesh structure, composition variables, nuclear burning rate, T_eff, L, M, T_c, ρ_c, and integrated luminosity from various types of nuclear burning.

MESAstar checks the structure and composition profiles of the model at the beginning of each timestep and, if necessary, adjusts the mesh. Cells may be split into two or more pieces, or they may be made larger by merging two or more adjacent cells. The overall remeshing algorithm is designed such that most cells are not changed during a typical remesh. This minimizes numerical diffusion and tends to help convergence. Remeshing is divided into a planning stage and an adjustment stage.

The planning stage determines which cells to split or merge based on allowed changes between adjacent cells. Mesh revisions minimize the number of splits and maximize the number of merges while ensuring that the magnitudes, Δ, of differences between any two adjacent cells are below specific thresholds: Δlog P < θ_P, Δlog T < θ_T, and Δlog [X(⁴He) + X(⁴He₀)] < θ_He, where X(⁴He) is the helium mass fraction and X(⁴He₀) sets an effective lower limit on the sensitivity to the helium abundance. The default thresholds are θ_P = 1/30, θ_T = 1/80, θ_He = 1/20, and X(⁴He₀) = 0.01. Options are available for specifying allowed changes between cells for other mass fractions, Δ∇_ad and Δlog(T/(T + T₀)) for arbitrary T₀.

Local reductions in the magnitude of allowed changes will place higher resolution in desired regions of the star. For example, the default is to increase resolution in regions of nuclear burning having Δlog _nuc large compared to Δlog P. This increase takes effect at a minimum log _nuc = −2 and increases to a maximum factor of four in resolution for log _nuc ⩾ 4. The size and range of enhancement can also be set for various specific types of burning. Similarly, it is possible to increase resolution near the boundaries of convection zones over a distance measured in units of the pressure scale height. Different enhancements and distances can be specified for above and below the upper and lower boundaries of zones with or without burning. There are also options to increase spatial resolution in regions having Δlog X_i large compared to Δlog P, or near locations where there are spatial gradients in the most abundant species. Finally, further splitting is done as necessary to limit the relative sizes of adjacent cells.

The adjustment stage executes the remesh plan. Cells to be split are constructed by first performing a monotonicity-preserving cubic interpolation (Steffen 1990) in mass to obtain the luminosities and enclosed volumes at the new cell boundaries. The new densities are then calculated from the new cell masses and volumes, as shown in Equation (4). Next, new composition mass fraction vectors are calculated. For cells being merged, this is straightforward. For cells being split, neighboring cells are used to form a linear approximation of mass fraction for each species as a function of mass coordinate within the cell. The slopes are adjusted so that the mass fractions sum to one everywhere, and the functions are integrated over the new cell mass to determine the abundances.

Finally, the method for calculating the new temperature varies according to electron degeneracy. As the electrons become degenerate (i.e., η>0), split cells simply inherit their temperature while merged cells take on the mass average of their constituent temperatures. If the electrons are not degenerate (i.e., η < 0), then a reconstruction parabola is created for the specific internal energy profile of the parent and its neighbor cells (Stiriba 2003). The parabola is integrated over the new cell to find its total internal energy. The new cell temperature is determined by repeatedly calling the eos module using the new composition and density with trial temperatures until the desired internal energy is found.

Mass adjustment for mass loss or accretion is done at each timestep before solving the equations for stellar structure and composition. MESAstar offers a variety of ways to set the rate of mass change $\dot{M}$. A constant mass accretion or mass-loss rate may be specified in the input files (see Section 6.1). Implementations of Reimers (1975) for red giants, Blöcker (1995) for AGB stars, de Jager et al. (1988) for a range of stars in the H-R diagram, mass loss for massive stars by Glebbeek et al. (2009), Vink et al. (2001), Nugis & Lamers (2000), and Nieuwenhuijzen & de Jager (1990), supersonic mass loss inspired by Prialnik & Kovetz (1995), and super-Eddington mass loss (Paczynski & Proszynski 1986) are available options. An arbitrary mass accretion or mass-loss scheme may be implemented by writing a new module. An example of such a routine provided with MESAstar is Mattsson et al. (2010) mass loss for carbon stars. Finally, one may write a wrapper program that calculates $\dot{M}$ for each timestep and then calls the MESAstar module.

Since MESAstar allows for simulations with a fixed (and unmodeled) inner mass, M_c, the total mass is M = M_c + M_m, where M_m is the modeled mass. For cell k, MESAstar stores the relative cell mass dq_k = dm_k/M_m and the relative mass interior to a cell face q_k = m_k/M_m = 1 − ∑^i=k−1_i=1dq_i (see Figure 9). Rather than evaluate dq_k as q_k − q_k+1, it is essential to define q in terms of dq to maintain accuracy (Lesaffre et al. 2006). For example, in the outer envelope of a star where the q_k approach 1, the dq_k can be 10⁻¹² or smaller. Subtraction of two adjacent q_k to find a dq_k leads to a intolerable loss of precision.

After a change in mass, δM, has been determined, the mass structure of the stellar model is modified. This procedure changes the mass location of some cells and revises the composition of those cells to match their new location. It does not add or remove cells, nor does it change the initial trial solution for the structure variables such as ρ, T, r, or L. The mass structure is divided into an inner (usually the central regions of the star), an intermediate, and an outer region (usually the stellar envelope). The boundaries of the inner and outer regions are initially set according to temperature, with defaults of log T = 6 for the inner boundary and log T = 5 for the outer boundary. This range is automatically expanded, for enhanced numerical stability, if the mass in the intermediate region is not significantly larger than δM. The range is first enlarged by moving the outer boundary to the surface. If the enclosed mass in the intermediate region is still too small, then the inner boundary can be moved inward subject to certain limits. One limit is that the inner boundary does not cross a region of the model where the composition changes rapidly. Another limit is that the fractional mass of the intermediate region cannot change by more than a factor of two from its previous value nor exceed 10% of the total mass.

Once the regions have been defined, the dq_k are updated. In the inner region the dq_k are rescaled by M/(M + δM). Thus, dm_k, m_k, and X_k have the same values before and after a change in mass to eliminate the possibility of unwanted numerical mixing in the center. In the outer region, cells retain the same value of dq_k to improve convergence in the high entropy regions of the star (Sugimoto et al. 1981). The dq_k in the intermediate region are scaled so that ∑dq_k = 1. The composition of cells in the intermediate and outer regions is then updated. In the case of mass accretion, the composition of the outermost cells whose enclosed mass totals δM is set to match the specified accretion abundances. Cells that were part of the old structure have their compositions set to match the previous composition.

We examined the resolution convergence properties of a 1 M_☉ model by varying the parameters for mesh refinement and timestepping. The mesh refinement parameter multiplies the limits for variable changes across mesh cells and is closely correlated with the cell size. The timestepping parameter controls the tolerance of the cell average of the relative changes between timesteps in log ρ, log T, and log R (see Section 6.4) and is closely correlated with the timestep. We varied the mesh refinement and timestepping controls in tandem through a parameter C, which is a multiplicative factor on their default values of 1 and 10⁻⁴, respectively. Therefore, C is anti-correlated with the time and space resolution.

Table 7 and Figure 10 detail the convergence properties with C of a solar metallicity, 1.0 M_☉ model with an η_R = 0.5 Reimers mass-loss model (Reimers 1975; see Section 6.6). These calculations begin at the ZAMS and are terminated at 11.0 Gyr, when the model stars are turning off the main sequence. As a measure of convergence, we use the difference, ξ, between a quantity at a given resolution and the quantity at the highest resolution considered (C = 1/16). In order to determine how convergence depends on resolution (|ξ| ∝ C^α), we determine the order of convergence, α, for increasingly resolved pairs in Table 7:

$$\begin{equation} \alpha = \log \left(\frac{\xi _{\rm fine}}{\xi _{\rm coarse}} \right) \Big / \log \left(\frac{\rm C_{fine}}{\rm C_{coarse}}\right) \enskip. \end{equation} \tag{ 17 }$$

The convergence orders show that all values converge linearly at large values of C and display super-linear convergence (α∼ 1.6) at smaller values of C. These convergence orders are plausible given that we use a first-order time integration scheme and a finite volume differencing scheme that is second-order accurate on uniform grids.

Zoom In Zoom Out Reset image size

Table 8 and Figure 11 detail the same stellar models as a function of C except the calculations are stopped at L = 100 L_☉, when the stars are on the RGB. Table 8 suggests the age of the star converges linearly at larger values of C and super-linearly (α∼ 1.5) at smaller values of C. However, T_eff, T_c, ρ_c, M, and M_He all display oscillatory behavior about the C = 1/16 solution, suggesting factors other than spacetime resolution are dominating the error at this stage of the evolution. Such factors could be limits in the precision attained by interpolation in the various tables (e.g., 4 significant figures for opacities, ∼6 significant figures for the OPAL and SCVH EOS), or small changes in boundary conditions.

Zoom In Zoom Out Reset image size

The lower panel of Figure 12 shows the sound speed profile for the 100 L_☉ model with C = 1/16, our highest resolution case. The helium core and convective zone boundaries are labeled. The mass interior to the convective zone boundary is 0.95 M_☉. The upper panel of Figure 12 shows the convergence properties of the sound speed profile with resolution in the hydrogen layer. Both the C = 1 and C = 1/2 profiles have sound speeds that are smaller than the C = 1/16 profile, while the C = 1/4 and C = 1/8 profiles have larger sound speeds. Note that the difference between the various profiles becomes less as the parameter C is made smaller, indicating that the convergence rate is becoming smaller. This suggests that factors other than spacetime resolution are dominating the convergence rates at this stage of the evolution. Again, this could be due to the precision attained by table interpolations, small changes in boundary conditions, or differencing errors.

Zoom In Zoom Out Reset image size

MESA modules are designed to be thread-safe (see Section 2), thereby enabling parallel execution. Table 9 lists the execution times in seconds of several specific tasks from four identical evolutionary calculations, each with a different number of threads (one thread per core). The essential difference between modules that scale as the inverse of the number of threads (eos and net) and those that do not (e.g., kap and neu) is the ratio of the overhead associated with parallel execution to the actual time required for each module to perform its calculation. The time required to complete serial tasks sets a lower limit on the total execution time of MESAstar.

By default, MESAstar provides alphanumeric output at regular intervals. In addition, it provides the option for concurrent graphical output. PGstar uses the PGPLOT²⁶ library to create on-screen plots or images in PNG format that can be post-processed into animations of an evolutionary sequence. A wide variety of visualization options are provided and these are all configurable through the PGstar inlist. For example, the PGstar window can simultaneously hold: an H-R diagram, a T_c–ρ_c diagram, and interior profiles of physical variables, such as nuclear energy generation, and composition. Animation is very useful for visualizing complex, time-dependent processes. For example, view the short selection from the MESA Web site that shows the He core flash in a 1 M_☉ star.²⁷

MESAstar has sufficiently broad input physics to compute the evolution of low-mass stars and substellar objects down to Jupiter's mass (≈10⁻³ M_☉), as well as complete evolutionary sequences of low-mass stars (M ≲ 2 M_☉) from the PMS to the white dwarf cooling curve without any intervention. Figure 13 shows evolutionary tracks in the H-R diagram for 1 and 1.25 M_☉ models with Z = 0.01.²⁸ The 1.25 M_☉ model exhibits a late He-shell flash during the pre-white dwarf phase.

Zoom In Zoom Out Reset image size

Figure 14 provides further examples, spanning 0.9–2 M_☉ at Z = 0.02; for clarity, the PMS portion of the tracks was removed and the runs were terminated after the models left the thermally pulsating asymptotic giant branch (TP-AGB). The bottom panel shows the evolution in the T_c–ρ_c plane, exhibiting the convergence of the 0.9, 1.2, and 1.5 M_☉ models to nearly identical, degenerate, helium cores when on the RGB. The 2 M_☉ model ignites He at a lower level of degeneracy.

Zoom In Zoom Out Reset image size

Though the required MESAstar timesteps get short (≈ hours) during the off-center He core flash in the 0.9, 1.2 and 1.5 M_☉ models, the stellar model does not become dynamic; the entropy change timescales are always longer than the local dynamical time (Thomas 1967; Serenelli & Weiss 2005; Shen & Bildsten 2009; Mocak et al. 2009). The reduction of hydrostatic pressure in the core at the onset of flashes leads to the adiabatic expansion of the core, visible as the drop in T_c at constant degeneracy. Successive He flashes (Thomas 1967; Serenelli & Weiss 2005) work their way into the core over a 2 × 10⁶ year timescale, eventually heating it (at nearly constant pressure; see the dashed red line in the bottom panel of Figure 14) to ignition and arrival onto the horizontal branch. The further evolution during He burning and the thermal pulses is seen in the bottom panel, where the small changes in the C/O core during thermal pulses on the AGB are resolved in MESAstar.

We start the more detailed calculations and comparisons to previous work in Section 7.1.1 by displaying the MESAstar PMS evolution of low-mass stars, brown dwarfs and giant planets, and comparing to prior results of Baraffe et al. (2003).

7.1.1. Low-mass Pre-main Sequence Stars, Contracting Brown Dwarfs and Giant Planets

The PMS evolution of low-mass stars (Burrows et al. 1989; D'Antona & Mazzitelli 1994; Baraffe et al. 1998) and gravitationally contracting brown dwarfs and giant planets (Burrows et al. 1997; Chabrier et al. 2000; Chabrier & Baraffe 2000; Burrows et al. 2001) has been studied extensively. The lasting importance of these problems motivates us to ensure that MESAstar can successfully perform these evolutions.

Figure 15 shows the evolution of PMS stars with masses of 0.08 M_☉ < M < 1 M_☉. Each solid line starts on the PMS Hayashi track for a fixed mass, M, and ends at an age of 1 Gyr. All stars with M ⩾ 0.2 M_☉ have reached the ZAMS (L = L_nuc) by this time. The red circles show the location where the ⁷Li is depleted by a factor of 100, but only for M < 0.5 M_☉. Stars with M>0.5 M_☉ deplete their ⁷Li after the core becomes radiative (D'Antona & Mazzitelli 1994; Chabrier et al. 1996; Bildsten et al. 1997), adding an uncertain dependence on convective overshoot that we do not investigate here.

Zoom In Zoom Out Reset image size

Lower-mass stars (M < 0.3 M_☉) remain fully convective throughout their PMS and arrival on the main sequence. In Figure 16 we show the evolution in the T_c − ρ_c plane for these objects. The light solid line in the upper left denotes the T_c ∝ ρ^1/3_c relation expected during the Kelvin–Helmholtz contraction for a fixed mass non-degenerate star. Deviations from this relation occur when electron degeneracy occurs, which is shown by the gray line at η ≈ _F/k_BT = 4, roughly where the electron degeneracy has increased the electron pressure to twice that of an ideal electron gas. We extend the mass range down to M = 0.01 M_☉ to reveal the distinction between main sequence stars and brown dwarfs. That distinction becomes clearer in Figure 17 which shows the L evolution for a range of stars with M < 0.3 M_☉. Only the M>0.08 M_☉ stars asymptote at late times to a constant L, whereas the others continue to fade. We also exhibit the expected scaling for a contracting fully convective star with a constant T_eff, L ∝ t^−2/3. At the D mass fraction of this calculation, X_D = 3 × 10⁻⁵, the onset of D burning provides some luminosity for a finite time, causing the evident kink at early times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

MESAstar models evolved from the PMS Hayashi line to an age of 10 Gyr with masses ranging from 0.09 to 0.001 M_☉ are compared with the models of Baraffe et al. (2003, BCBAH) in Figure 18. Ages in increasing powers of 10 are marked by filled circles along each track from 1 Myr to 10 Gyr. For comparison, separate points from the BCBAH evolutionary models are plotted as plus symbols ("+"). So as to match the choice of BCBAH, we set the D mass fraction at X_D = 2 × 10⁻⁵ (Chabrier et al. 2000). Evolution at the youngest ages is uncertain due to different assumptions regarding D burning but beyond 10 Myr the MESAstar and the BCBAH models overlap at almost every point. Note that the BCBAH models were only evolved to 5 Gyr for the two lowest masses shown in Figure 18.

Zoom In Zoom Out Reset image size

7.1.2. Code Comparisons of 0.8 M_☉ and 1 M_☉ Models

The Stellar Code Calibration Project (Weiss et al. 2007) was created to provide insight into the consistency of results obtained from different state-of-the-art stellar evolution codes. The contributors performed a series of stellar evolution calculations with the physics choices held constant to the greatest extent possible. This section compares MESAstar models with published results from that project for two specific cases. The comparison codes are BaSTI/FRANEC (Pietrinferni et al. 2004), DSEP (Dotter et al. 2007), and GARSTEC (Weiss & Schlattl 2008). MESAstar models lie within the range exhibited by BaSTI/FRANEC, DSEP, and GARSTEC in these comparisons.

Two examples are shown here, 0.8 M_☉, Z = 10⁻⁴ evolutionary tracks (Figure 19) and 1 M_☉, Z = 0.02 evolutionary tracks (Figure 20). In both cases the models were evolved from the pre-MS to the onset of the He core flash. The models assume, as much as possible, the same nuclear reaction rates (NACRE), opacities (OPAL and Alexander & Ferguson 1994), equation of state (FreeEOS), and mixing length (α_MLT = 1.6). These tests do not represent the best models for the various codes. Instead, the goal of the comparisons was to see how consistent the codes would be when using simple assumptions and comparable input physics (Weiss et al. 2007). While the agreement is good in most respects, in temporal resolution there is a discrepancy.

Zoom In Zoom Out Reset image size

The H-R diagram of the 0.8 M_☉, Z = 10⁻⁴ case in Figure 19 is nearly identical for all four tracks except near the main sequence turnoff, where DSEP and BaSTI/FRANEC are hotter than MESAstar and GARSTEC. These models have essentially no convection during the main sequence and there is remarkably little scatter during this phase. In the T_c-L plane, the models are almost indistinguishable until they enter the red giant phase at L ≈ 10 L_☉, where the central temperatures differ slightly only to re-converge at maximum luminosity on the RGB. Finally, the lifetime–luminosity plane indicates that the four codes split into two pairs with one pair shorter lived by ≈5% than the other pair. It is beyond our scope here to explain the reasons for these differences; the purpose of the present comparison is to indicate that MESAstar produces results that are consistent with the range exhibited among the other three codes.

Convection plays a more prominent role in the 1 M_☉, Z = 0.02 case, as shown in Figure 20, and the scatter is greater than in the Z = 10⁻⁴ case. The BaSTI/FRANEC model is hotter than the other three models. Treatment of convection and, in particular, the resolution of the surface convection zone is primarily responsible for the spread seen in the main sequence portion of the tracks.

Zoom In Zoom Out Reset image size

In both cases, the central conditions are very similar until the models become red giants. In the Z = 0.02 case, the range of lifetimes is somewhat reduced compared to the Z = 10⁻⁴ case with BaSTI/FRANEC and MESAstar shortest (though the order is reversed with respect to the Z = 10⁻⁴ case) followed by DSEP and then GARSTEC.

7.1.3. The MESAstar Solar Model

MESAstar performs a Solar model calibration by iterating on the difference between the final model and the adopted solar parameters of L_☉ and R_☉ (from Bahcall et al. 2005) and the surface value of Z_s/X_s from Grevesse & Sauval (1998) at 4.57 Gyr. This is done by iteratively varying α_MLT and the initial Y_i and Z_i values (all for the abundance ratios of Grevesse & Sauval 1998), while including diffusion.

The properties of the converged model (which reaches the desired parameters to better than one part in 10⁵) are shown in Table 10, and match the measured depth, R_CZ, of the surface convection zone within 1σ and the surface helium abundance, Y_s, within 2σ (Bahcall et al. 2005). The difference between the model and the helioseismologically inferred solar sound speed profile is compared with similar results from Bahcall et al. (1998, BBP98) and Serenelli et al. (2009, S09) in Figure 21 demonstrating that MESAstar is capable of stellar evolution calculations at the level of one part in 10³ demonstrated by others (Basu & Antia 2008).

Zoom In Zoom Out Reset image size

Table 10. MESAstar Standard Solar Model at 4.57 Gyr

Quantity	Value
Converged input parameters
αMLT	1.9179113764
Yi	0.2744267987
Zi	0.0191292323
Properties of converged model
(Z/X)s	0.02293
Xs	0.73973
Ys	0.24331
Zs	0.01696
Xc	0.33550
Zc	0.02125
RCZ/R☉	0.71398
log ρc	2.18644
log Pc	17.3695
log Tc	7.19518
rms[(cModel − c☉)/c☉]	0.00093

Download table as:  ASCIITypeset image

MESAstar can calculate the evolution of intermediate-mass stars (2 ≲ M/ M_☉ ≲ 10) through the He-core burning phase and the advanced He-shell burning AGB phase. MESAstar produces results compatible with published results from existing stellar evolution codes.

We start by showing in Figure 22 a grid of MESAstar evolutionary tracks with masses ranging from 2 to 10 M_☉ with Z = 0.02. The top panel shows the evolution in the H-R diagram, while the bottom panel shows the evolution in the T_c − ρ_c diagram. The 8 and 10 M_☉ models start to ignite carbon burning off center, whereas the 2–7 M_☉ models do not and will presumably go on to form C/O core white dwarfs. The lack of a complete treatment in MESAstar of liquid diffusion inhibits our ability to verify the resulting white dwarf cooling sequences from MESAstar at this time.

Zoom In Zoom Out Reset image size

7.2.1. Comparison of EVOL and MESAstar

We compare M_i = 2 M_☉, Z = 0.01 stellar models from MESAstar and EVOL (Blöcker 1995; Herwig 2004; Herwig & Austin 2004) starting from the PMS to the tip of the thermal pulse AGB (TP-AGB). Both codes employed the exponentially decaying overshoot mixing treatment described by Herwig (2000; see Section 5.2) at all convective boundaries with f = 0.014, except during the third dredge-up where we adopt f = 0.126 at the bottom of the convective envelope to account for the formation of a ¹³C pocket, and at the bottom of the He-shell flash convection zone we use f = 0.008 (Herwig 2005).

In both codes we use the mass-loss formula of Blöcker (1995, see Section 6.6). Thermal pulses start at a slightly lower core mass, and hence luminosity, in the EVOL model. In order to maintain similar envelope mass evolution through the TP-AGB, the parameter η_Bl in the mass-loss formula was set to 0.05 in MESAstar and 0.1 in EVOL. Every effort has been made to tailor the MESAstar model to the EVOL model. However, the AGB evolution is very sensitive to the initial core mass, which depends on the mixing assumptions and their numerical implementation during the preceding He-core burning phase. Consequently, small differences on the TP-AGB are unavoidable when comparing tracks from two codes.

As shown in Figure 23, the EVOL and MESAstar tracks compare well in the H-R diagram. Table 11 shows that key properties differ by less than 5%. MESAstar has the ability to impose a minimum size on convection zones below which overshoot mixing is ignored. EVOL does not have such limits, leading to more mixing of He into the core and, hence, the ≈4% larger age of the EVOL sequence at the first thermal pulse.

Zoom In Zoom Out Reset image size

The TP-AGB is characterized by recurrent thermonuclear instabilities of the He-shell, leading to complex mixing and nucleosynthesis. These processes are properly represented in MESAstar calculations, as revealed in Figure 24. The ability of MESAstar to calculate the evolution of stellar parameters in a smooth and continuous manner even during the advanced thermal pulse phases and beyond is demonstrated in Figure 25. The top panel shows the evolution in the H-R diagram, whereas the bottom panel shows the evolution of the conditions in the C/O core. The adiabatic cooling in the C/O core that occurs during the He flash (due to the pressure dropping at the surface of the C/O core) is evident in the downturns that are parallel to the line of constant degeneracy (which is also the adiabatic slope). The overall trend of increasing ρ_c reflects the growing C/O core mass, which for this model is shown in the top panel of Figure 24.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 11. Comparison of MESAstar and EVOL Models with M_i = 2 M_☉, Z = 0.01

Quantity	MESAstar	EVOL
Main-sequence lifetime (Gyr)	0.939	0.962
Deepest penetration of first dredge-up ( M☉)	0.328	0.327
H-free core mass at the end of He-core burning ( M☉)	0.466	0.454
Core mass at first thermal pulse ( M☉)	0.504	0.481
Age at first thermal pulse (Gyr)	1.269	1.328
Core mass at second thermal pulse with DUP ( M☉)	0.563	0.563
Following interpulse time (1000 yr)	116	106
Following pulse-to-pulse core growth (10−3 M☉)	6.4	6.9
Dredge-up mass at following pulse (10−3 M☉)	1.1	1.3

Download table as:  ASCIITypeset image

An example of the evolution of convection zones, shell burning, and total luminosities as well as core boundaries for two subsequent thermal pulses is shown in Figure 26 as a function of model number; compare to Figure 3 in Herwig (2005). Quantitative comparison of interpulse time, core growth, and dredge-up amount (see Table 11) shows excellent agreement between the MESAstar thermal pulses and the equivalent pulses in the EVOL sequence.

Zoom In Zoom Out Reset image size

Another important property of the He-shell flashes is the intershell abundance as a result of the convective mixing and burning. Again, the comparison of results from both codes shows good agreement, which is expected since they both implement the same overshooting mixing assumptions for the He-shell flash convection zone. A consequence of the third dredge-up is the gradual increase of the envelope C/O ratio as thermal pulses repeatedly occur. Since the 2 M_☉ models do not experience hot-bottom burning, the evolution of this ratio is an effective probe of the cumulative efficiency of the third dredge-up in these simulations.²⁹ The top panel of Figure 27 shows the surface C/O ratio evolution according to EVOL (dashed red line) and MESAstar (solid black line). They are in good agreement, e.g., in terms of the time period over which the third dredge-up occurs and the amount by which C/O increases. The mass-loss history over the same time period, shown in the bottom panel of Figure 27, is similar by design.

Zoom In Zoom Out Reset image size

7.2.2. Interior Structure of Slowly Pulsating B Stars and Beta Cepheids

The advent of space-based asteroseismology for main-sequence B stars with the CoRot (Degroote et al. 2009) and Kepler (Gilliland et al. 2010) satellites is probing the slowly pulsating B stars (SPBs, M ≈ 3–8 M_☉) and the more massive (M ≈ 7–20 M_☉) β Cepheids (Degroote et al. 2009, 2010). These stars are all undergoing main-sequence H burning and are unstably pulsating due to the κ mechanism from the Fe-group opacity bump at T ≈ 2 × 10⁵ K (Dziembowski et al. 1993). The observed modes have finite amplitudes deep in the stellar core, demanding a full interior model for mode frequency (and stability) prediction (Dziembowski et al. 1993; Pamyatnykh et al. 2004).

These papers provide a few specific models that allow a direct comparison to the MESAstar prediction of the Brunt–Väisälä frequency

$$\begin{equation} N^2=g\left({1\over \Gamma _1} {d\ln P \over dr}-{d\ln \rho \over dr}\right)=g\left(-{g\over c_s^2}-{d\ln \rho \over dr}\right), \end{equation} \tag{ 18 }$$

where c²_s = Γ₁P/ρ is the adiabatic sound speed, and we used hydrostatic balance, dP/dr = −ρg. Numerically, these are obtained by interpolating the sound speed at the cell boundary, whereas dln ρ/dr is estimated by numerical differencing and then smoothed. This method naturally captures the extra restoring force from composition gradients, especially relevant in these evolving stars that leave a He-rich radiative region above the retreating convective core during the main sequence.

Our first comparison is to Dziembowski et al.'s (1993) M = 4 M_☉ main-sequence star with Z = 0.02 at a time when the hydrogen abundance in the convective core is X_c = 0.37. With no overshoot from the convective core, Dziembowski et al. (1993) found log L/L_☉ = 2.51 and log T_eff = 4.142 whereas MESAstar gives log L/L_☉ = 2.50 and log T_eff = 4.125. The top panel in Figure 28 compares the MESAstar results (solid line) to the values (green circles) from Figure 3 of Dziembowski et al. (1993). The agreement is remarkable as an integral test of MESAstar. The bottom panel of Figure 28 is a comparison to the more massive M = 9.858 M_☉ main-sequence star with Z = 0.015 from Figure 5 of Pamyatnykh et al. (2004) at an age (15.7 Myr) when X_c = 0.2414 with log L/L_☉ = 3.969 and log T_eff = 4.3553. MESAstar gave log L/L_☉ = 3.966, log T_eff = 4.358, and an age of 16.4 Myr at the same value of X_c. These comparisons highlight the readiness of MESAstar for adiabatic asteroseismological studies of main-sequence stars.

Zoom In Zoom Out Reset image size

To explore MESAstar's results in this mass range, models of 15 M_☉, 20 M_☉, and 25 M_☉ of solar metallicity and 1000 M_☉ of zero metallicity were evolved from the Hayashi track to the onset of core collapse. Nuclear reactions are treated with the 21 isotope reaction network, inspired by the 19 isotope network in Weaver et al. (1978), that is capable of efficiently generating accurate nuclear energy generation rates from hydrogen burning through silicon burning (see Section 4.5). This network includes linkages for PP-I, steady-state CNO cycles, a standard α-chain, heavy ion reactions, and aspects of photodisintegration into ⁵⁴Fe. Atmospheres are treated as a τ=2/3 Eddington gray surface as described in Section 5.3. Mass loss for the solar metallicity stars uses the combined results of Glebbeek et al. (2009), Vink et al. (2001), Nugis & Lamers (2000), and Nieuwenhuijzen & de Jager (1990), as described in Section 6.6. These massive star models are non-rotating, use no semi-convection, employ a mixing length parameter of α_MLT = 1.6, and adopt f = 0.01 for exponential diffusive overshoot (see Section 5.2) for convective regions that are either burning hydrogen or are not burning.

Most of this section consists of comparisons to results from other stellar evolution codes. However, for consistency (and completeness), we show in Figure 29 the H-R diagram and central condition evolution of 10–100 M_☉ stars from the PMS to the end of core helium burning. Though these are stars with Z = 0.02, we turned off mass loss during this calculation so that the plot would be easier to read and of some pedagogical use. The tendency of T_c to scale with ρ^1/3_c (also a constant radiation entropy) during these stages of evolution is expected from hydrostatic balance with only a mildly changing mean molecular weight. The rest of the calculations in this section included mass loss as described above.

Zoom In Zoom Out Reset image size

7.3.1. 25 M_☉ Model Comparisons

Figure 30 shows the T_c–ρ_c evolution in M_i = 25 M_☉ solar metallicity models from MESAstar, KEPLER (A. Heger 2010, private communication), Hirschi et al. (2004), and FRANEC (Limongi & Chieffi 2006) from helium burning until iron-core collapse. The curves fall below the T_c ∝ ρ^1/3_c scaling relation as the mean molecular weight increases due to the subsequent burning stages. The curves are also punctuated with non-monotonic behavior when nuclear fuels are first ignited in shells. Figure 30 shows that MESAstar produces core evolutionary tracks consistent with other pre-supernova efforts. The bump in the MESAstar curve around carbon burning is due to the development of central convection whereas the other codes do not (although see Figure 2 of Limongi et al. 2000). The development of a convective core during carbon burning depends on the carbon abundance left over from core helium burning (Limongi et al. 2000).

Zoom In Zoom Out Reset image size

The mass fraction profiles of the inner 2.5 M_☉ of this M_i = 25 M_☉ model are shown in Figure 31 at the onset of core collapse. At the time of these plots, the infall speed has reached ≈1000 km s^-1 just inside the iron core (at m = 1.5 M_☉) and the electron fraction, Y_e, has dropped below ≈0.48. The oxygen shell lies at 1.88 ⩽ m/ M_☉ ⩽ 2.5, the silicon shell between 1.61 ⩽ m/ M_☉ ⩽ 1.88, and the iron core at m ⩽ 1.61 M_☉. Figure 32 shows T, ρ, S, the radial velocity, the infall timescale, and Y_e of this inner 2.5 M_☉. Note the entropy decrements at the oxygen, silicon, and iron core boundaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 33 summarizes the history of the inner 7 M_☉ of this M_i = 25 M_☉ model as a function of interior mass (left y-axis). Evolution is measured by the logarithm of time (in years) remaining until the death of the star as a supernova (x-axis), which reveals the late burning stages. Levels of red and blue shading indicate the magnitude of the net energy generation (nuclear energy generation minus neutrino losses), with red reflecting positive values and blue indicating negative ones. The vertical lines indicate regions that are fully convective. Note the appearance of a convective envelope characteristic of a red supergiant late during helium burning. Abundance profiles of key isotopes during the major burning stages are shown (right y-axis). The hydrogen core shrinks toward the end of hydrogen burning, and the helium core grows as helium is depleted. The total mass shrinks to about M = 12 M_☉ due to mass loss.

Zoom In Zoom Out Reset image size

Table 12. Massive Star Core Burning Lifetime Comparison

Core Burning	Lifetime (years)
Element	HMM	WHW	LSC	MESA
Mi = 15 M☉
H	1.13	1.11	1.07	1.14	×107
He	1.34	1.97	1.4	1.25	×106
C	3.92	2.03	2.6	4.23	×103
Ne	3.08	0.732	2.00	3.61
O	2.43	2.58	2.43	4.10
Si	2.14	5.01	2.14	0.810	×10−2
Mi = 20 M☉
H	7.95	8.13	7.48	8.01	×106
He	8.75	11.7	9.3	8.10	×105
C	9.56	9.76	14.5	13.5	×103
Ne	0.193	0.599	1.46	0.916
O	0.476	1.25	0.72	0.751
Si	9.52	31.5	3.50	3.32	×10−3
Mi = 25 M☉
H	6.55	6.706	5.936	6.38	×106
He	6.85	8.395	6.85	6.30	×105
C	3.17	5.222	9.72	9.07	×102
Ne	0.882	0.891	0.77	0.202
O	0.318	0.402	0.33	0.402
Si	3.34	2.01	3.41	3.10	×10−3

References. HMM: Hirschi et al. (2004); WHW: Woosley et al. (2002); LSC: Limongi et al. (2000); MESA: this paper.

Download table as:  ASCIITypeset image

7.3.2. Comparison of 15, 20, and 25 M_☉ Models

7.3.3. 1000 M_☉ Metal-free Star Capabilities

We close this section with a demonstration of MESAstar's capabilities by describing the unlikely scenario of a purely metal-free stellar evolution of an M_i = 1000 M_☉ star. The T_c–ρ_c trajectory for a 1000 M_☉, zero metallicity, zero mass-loss model is shown in the left panel of Figure 34. The starting time point is in the lower left corner and the final model, at the onset of core collapse, is in the upper right at very high values of T_c and ρ_c. Fluid elements in the region to the left of the red dashed line have Γ₁ < 4/3. When the center enters this region, the central portions of the star become dynamically unstable and begin to contract. However, the entire star does not collapse because the infalling regions become denser and hotter, causing the central region to leave the Γ₁ < 4/3 region and the infall to slow. Now another part of the star moves into the Γ₁< 4/3 region and begins to infall at high velocity. The net result is that the region where Γ₁< 4/3 starts at the center and moves outward. The right panel of Figure 34 shows the material speed and Γ₁ profiles for the final model, where the infalling region is now at m ≈ 480 M_☉.

Zoom In Zoom Out Reset image size

The global history of the 1000 M_☉ model as a function of time is shown in the left panel of Figure 35. A convective envelope appears during late helium burning. Abundance profiles of key isotopes during the major burning stages are shown (right y-axis). Note the short carbon burning era. At late times the core photodisintegrates to ⁴He instead of creating ⁵⁶Ni because of the lower central densities encountered in these supermassive progenitors. This also partially causes the large endothermic central regions of the star.

Zoom In Zoom Out Reset image size

MESAstar can be used to examine how a star responds to mass loss or accretion (see Section 6.6). This opens up a large variety of possible applications, including accretion onto white dwarfs for classical novae and thermonuclear supernovae, mass transfer in tight stellar binaries, and learning the response of a star to sudden mass loss. We show two examples where MESAstar's results can be compared to previous work. The first is a mass-transfer scenario relevant to P_orb < 2 hr cataclysmic variables, and the second is the response of a neutron star to accretion of pure He.

7.4.1. Mass Transfer in a Binary

To illustrate MESAstar's ability to calculate the impact of mass loss on a star, we model the evolution of a compact binary consisting of a Roche Lobe filling low-mass ZAMS (M < 0.2 M_☉) model and an accreting white dwarf with M_WD = 0.6 M_☉. These short orbital period (P_orb < 2 hr) cataclysmic variables are the end points of these mass transferring systems (Patterson 1998; Kolb & Baraffe 1999) and are now being discovered in large numbers in the SDSS database (more than 100 studied by Gänsicke et al. 2009).

We model the parameters of the binary system and the Roche lobe overflow triggered mass transfer rate $\dot{M}$ as in Madhusudhan et al. (2006). So as to compare to the previous work of Kolb & Baraffe (1999), we presume angular momentum losses from gravitational wave emission and keep the accreting WD mass fixed at its initial value, M_WD = 0.6 M_☉. The evolution of the donor star is carried out by MESAstar, using the τ = 100 atmosphere tables from atm. The evolution shown in Figure 36 is followed for over 6 Gyr until the donor has been reduced to a brown dwarf remnant of M ≈ 0.03 M_☉ (see Table 14). During that time, the binary period drops to a minimum value of 67.4 minutes and then increases, independent of the initial donor mass. This plot is very similar to Figure 1 of Kolb & Baraffe (1999). We also show in Table 14 the evolution in time of the main properties of the donor star and mass transfer rate of the M_i = 0.21 M_☉ model. Again, this agrees with the results in Table 2 of Kolb & Baraffe (1999). The prime differences can be attributed to a slightly different R(M) relation.

Zoom In Zoom Out Reset image size

Table 14. Mass Transfer History for M_i = 0.21 M_☉ and M_WD = 0.6 M_☉

Time (Gyr)	Porb(hr)	M/M☉	Teff (K)	log(L/L☉)	R/R☉	$\log \dot{M}$
0.00	2.1319	0.2100	3278	−2.2688	0.2279	−10.24
0.25	2.0962	0.1987	3262	−2.3041	0.2209	−10.39
0.50	2.0367	0.1887	3242	−2.3467	0.2129	−10.40
0.75	1.9770	0.1787	3217	−2.3939	0.2049	−10.41
1.00	1.9181	0.1693	3191	−2.4414	0.1971	−10.42
1.25	1.8560	0.1599	3167	−2.4904	0.1891	−10.44
1.50	1.7938	0.1510	3142	−2.5404	0.1814	−10.45
1.75	1.7299	0.1421	3113	−2.5952	0.1735	−10.46
2.00	1.6684	0.1336	3074	−2.6563	0.1659	−10.47
2.25	1.6097	0.1252	3018	−2.7277	0.1585	−10.48
2.50	1.5475	0.1170	2965	−2.8004	0.1510	−10.50
2.75	1.4829	0.1093	2916	−2.8737	0.1435	−10.51
3.00	1.4151	0.1016	2854	−2.9590	0.1358	−10.52
3.25	1.3484	0.0942	2774	−3.0577	0.1283	−10.53
3.50	1.2804	0.0868	2666	−3.1796	0.1207	−10.54
3.75	1.2172	0.0796	2525	−3.3273	0.1135	−10.54
4.00	1.1663	0.0726	2355	−3.4987	0.1071	−10.55
4.25	1.1345	0.0659	2161	−3.6915	0.1019	−10.58
4.50	1.1227	0.0596	1963	−3.8930	0.0980	−10.63
4.75	1.1291	0.0541	1771	−4.0949	0.0953	−10.70
5.00	1.1487	0.0495	1595	−4.2924	0.0937	−10.78
5.25	1.1752	0.0457	1445	−4.4725	0.0927	−10.86
5.50	1.2051	0.0426	1314	−4.6436	0.0922	−10.94
5.75	1.2379	0.0399	1201	−4.8016	0.0919	−11.02
6.00	1.2706	0.0377	1110	−4.9398	0.0918	−11.09
6.25	1.3035	0.0358	1028	−5.0720	0.0918	−11.16
6.50	1.3343	0.0343	969	−5.1749	0.0919	−11.22
6.75	1.3636	0.0328	921	−5.2630	0.0920	−11.29
7.00	1.3902	0.0316	879	−5.3425	0.0921	−11.34

Download table as:  ASCIITypeset image

7.4.2. Rapid Helium Accretion onto a Neutron Star

The outer envelope of an accreting neutron star is modeled in MESAstar by using non-zero boundary conditions M_c and L_c (see discussion in Section 6.2) at a finite radius R_c. This allows for a time-dependent calculation of the thermonuclear instability that yields Type I X-ray bursts (Strohmayer & Bildsten 2006) for those accretion rates where the burning is thermally unstable ($\dot{M}\le 10^{-8}\,M_\odot \ {\rm yr^{-1}}$). Such calculations have been performed with the KEPLER code (Woosley et al. 2004; Cyburt et al. 2010) and prove very valuable in direct comparisons to observed Type I X-ray burst recurrence times and light curves, especially for the H-rich accreting "clocked burster" GS 1826-24 (Heger et al. 2007). We focus here on pure He accretion, relevant to neutron stars in ultra-compact binaries, such as 4U 1820-30 (Cumming 2003).

For these simulations we set M_c = 1.4 M_☉, R_c = 10 km, L_c = 3.6 × 10³⁴ erg s^-1, and g = 2.39 × 10¹⁴ cm s^-2 (correcting for the gravitational redshift). The initial model consisted of 3 × 10²⁵g of pure ⁵⁶Fe and accreted pure He at $\dot{M}=3\times 10^{-9}\,M_\odot \ {\rm yr^{-1}}$. We require a slightly higher value of core luminosity $L_c/\dot{M}\approx 0.19 \ {\rm keV \ nucleon^{-1}}$ to reach the same ignition column depth (5 × 10⁸ g cm^-2) as Weinberg et al. (2006). We used 31 species in the nuclear reaction network, including the ¹²C bypass reaction chain ¹²C(p,γ)¹³N(α,p)¹⁶O and elements (²³Na, ²⁷Al, ³¹P, ³⁵Cl, and ³⁹K) that can appear as intermediates in (α, p)(p, γ) reactions and serve as the proton source for the ¹²C bypass (Weinberg et al. 2006).

Figure 37 shows a snapshot of the time history of the helium burning luminosity, L_He, which is periodic at the Type I burst recurrence time of 9.56 hr. This luminosity, as well as L, very quickly exceeds L_Edd, in which case we allow for mass loss via a wind (Paczynski & Proszynski 1986). We arbitrarily set our time coordinate to zero at the time of maximum luminosity, L, in the second burst after the start of accretion. The peak for L_He is at t = −0.0269 s, and L>L_Edd for the time interval −0.0047 s < t < 1.2169 s.³⁰

Zoom In Zoom Out Reset image size

Figure 38 shows the evolving temperature profile during the convective burning runaway, where time increases upward. Though not for the same ignition depth, this plot is very similar to Figure 2 of Weinberg et al. (2006), including the evolution of the location of the top of the convective zone (open squares). Weinberg et al. (2006) discussed in detail the onset of heat transport in the outer, thin, radiative layer that allows for the retreat of the top of the convective zone. This MESAstar result is the first numerical confirmation of this transition for a pure helium accretor and demonstrates our ability to obtain excellent time and mass resolution as shown in Figure 39. By using nonzero center boundary conditions so that the dq variables (see Section 6.6) cover only the relatively small envelope mass, we reach a mass resolution of ≈1.5 × 10⁻²⁰ M_☉. The timestep adjustment algorithms (Section 6.4) provide a smooth change from timesteps of almost an hour between bursts down to millisecond steps at peak luminosity (see the middle panel in Figure 39). The secular increase in the number of cells is to track the accumulation of the pile of ashes from each burst. The evolution of abundances at the base of the convective zone is shown in Figure 40 and exhibits the presence of the isotopes ³⁵Cl and ³⁹K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have performed simulations at lower accretion rates but these become dynamical events where the temperature rises on a local dynamical timescale and are beyond the present scope of MESAstar. While multi-dimensional hydrodynamical codes (e.g., Zingale et al. 2009) may be needed to follow the details of such an event, MESAstar can be used for studying the longer timescale, hydrostatic evolution leading up to the point where hydrodynamic effects become dominant.

Stellar evolution calculations (i.e., stellar evolution tracks and detailed information about the evolution of internal and global properties) are a basic tool that enable a broad range of research in astrophysics. Areas that critically depend on high-fidelity and modern stellar evolution include asteroseismology, nuclear astrophysics, galactic chemical evolution and population synthesis, compact objects, supernovae, stellar populations, stellar hydrodynamics, and stellar activity. New observational capabilities are emerging in these fields that place a high demand on exploration of stellar dependencies on metallicity and age. So, even though one-dimensional stellar evolution is a mature discipline, we continue to ask new questions of stars. The emergence of demand requires the construction of a general, modern stellar evolution code that combines the following advantages.

• 1.  
  Openess: should be open to any researcher, both to advance the pace of scientific discovery, but also to share the load of updating physics, fine-tuning, and further development.
• 2.  
  Modularity: should provide independent, reusable modules.
• 3.  
  Wide applicability: should be capable of calculating the evolution of stars in a wide range of environments, including low and massive stars, binaries, accreting, mass-losing stars, early and advanced phases of evolution etc. This will enable multi-problem, multi-object physics validation.
• 4.  
  Modern techniques: should employ modern numerical approaches, including high-order interpolation schemes, advanced AMR, simultaneous operator solution; should support well-defined interfaces for related applications, e.g., atmospheres, wind simulations, nucleosynthesis simulations, and hydrodynamics.
• 5.  
  Microphysics: should allow for up-to-date, wide-ranging, flexible, and modular micro-physics.
• 6.  
  Performance: should parallelize on present and future shared-memory, multi-core/thread and possibly hybrid architectures so that performance continues to grow within the new computational paradigm.

A tool that combines the above features is a significant research and education resource for stellar astrophysics. We acknowledge that some important aspects of stars are truly three dimensional, such as convection, rotation, and magnetism. Those applications remain in the realm of research frontiers with evolving understanding and insights, quite often profound. However, much remains to be gained scientifically (and pedagogically) by accurate one-dimensional calculations, and this is the present focus of MESA.

The MESA code library project is open. It explicitly invites participation from anybody (researchers, students, interested amateurs). Participation in MESA can take a wide range of forms, from just using a MESA release for a science project, to testing and debugging (i.e., report bugs, find fixes and submit them for inclusion into the next release) as well as taking on responsibility for the continued stewardship of certain aspects (modules) of the code. The participation of experienced stellar evolution experts is very welcome.

Users are encouraged to add to the capabilities of MESA, which will remain a community resource. However, use of MESA requires adherence to the "MESA code of conduct".

• 1.  
  That all publications and presentations (research, educational, or outreach) deriving from the use of MESA acknowledge this paper and the MESA Web site.
• 2.  
  That user modifications and additions are given back to the community.
• 3.  
  That users alert the MESA Council (see below) about their publications, either pre-release or at the time of publication.
• 4.  
  That users make available in a timely fashion (e.g., online at the MESA Web site) all information needed for others to recreate their MESA results—"open know how" to match "open source."
• 5.  
  That users agree to help others learn MESA, giving back as the project progresses.

Users are requested to identify themselves by name, email contact, and home location.

The MESA project began as an initiative to construct a reliable computational tool for stellar structure and evolution that takes full advantage of modern processor architectures, algorithms, and community engagement. The release of MESA has forced some explicit thinking of what structure is needed so as to achieve the mission of stewarding MESA in its use for scientific research, education, and outreach, while also enabling the development of new tools and ideas. The MESA operating principles are simple: be open in your scientific discussions, give credit to all contributors, and be prepared to give back to the community of users. We hope that this creates an environment where the young are encouraged to become engaged in a career-enhancing manner.

We have established the MESA Council that consists of those engaged in working toward the shared missions outlined here:

• 1.  
  Steward MESA. There are many ways, this will be done: supporting the contributors, maintaining the Web access and Web site updates, seeking enabling funding, holding yearly working groups that allow for continued engagement, documenting MESA development in the peer-reviewed literature, and sustaining advanced development.
• 2.  
  Interface with the user community. This starts with answering questions from users, developing a way to accept new code in an integrated fashion, maintain a user registry, and identify new MESA Council members from those most active and engaged in the intelligent use of MESA.
• 3.  
  Enable scientific research and education with MESA. Promote MESA and its goals, e.g., through scientific contributions at relevant conferences. Identify science opportunities that match MESA capabilities and facilitate and encourage appropriate collaborative activities. Track the science carried out by the community with MESA.