In-Class Assignment 3

In-Class Assignment 3#

Exploring MESA models#

Let’s look at the output from the MESA stellar evolution code for 3 different mass progenitors: 1, 8, and 15 \(M_\odot\). These models were all created using MESA Web a web-based interface to MESA for quick calculations.

To read in these files, we use py_mesa_reader

Note

py_mesa_reader is not available on PyPI, so it can be installed using pip as:

pip install git+https://github.com/wmwolf/py_mesa_reader

Credit: These exercises follow from work by Mike Zingale’s computational lectures.

import numpy as np
import matplotlib.pyplot as plt
#import mesa_reader as mr

MESA provides 2 types of output: profiles and a history. The profiles represent a snapshot of the start at one instance of time and give the stellar data as a function of radius or enclosed mass. The history provides some global quantities as a function of time throughout the entire evolution of the star. We’ll use both.

The model files are:

To make the management easier, we’ll create a container for each model that holds the history and profiles, processed by py_mesa_reader

class Model:
    def __init__(self, mass, profiles=None, history=None):
        self.mass = mass
        if history:
            self.history = mr.MesaData(history)
        self.profiles = []
        if profiles:
            for p in profiles:
                self.profiles.append(mr.MesaData(p))

Now read in all the data. For each mass we have 2 profiles and 1 history. The profiles were picked to roughly correspond to the midpoint of core H burning and core He burning.

models = []


#models.append(Model(1, profiles=["M1_default_profile8.data",
#                                 "M1_default_profile218.data"],
#                    history="M1_default_trimmed_history.data"))

#models.append(Model(8, profiles=["M8_basic_co_profile8.data",
#                                 "M8_basic_co_profile39.data"],
#                    history="M8_basic_co_trimmed_history.data"))

#models.append(Model(15, profiles=["M15_aprox21_profile8.data",
#                                  "M15_aprox21_profile19.data"],
#                    history="M15_aprox21_trimmed_history.data"))

Next, we will use these data to make some plots from the readings.

a#

HR diagram#

Lets start with an HR diagram. Below are some possible steps to do that. For reference, here are some common commands working with MESA-reader and here.

hint an efficient way to see which variables are available to plot is the Linux less command and specifically using the -S option:

less -S M1_default_trimmed_history.data 

# first, lets make sure we can access our data for a particular model history file
#one_solar_mass_model = models[0] ## for model array [0,1,2]
#one_solar_mass_model.history.Teff

# What other value do we need for the 1msun model history file to plot an HR diagram?
# find that array here

The matplotlib reference page is linked here.

# Okay, now we have two arrays for the 1msun model
# plot those data here. 

#fig, ax = plt.subplots()
# ax.loglog( fill in here)

# Plot all three models on the same HR diagram here and add a legend.
# Make sure your x-axis is going the correct direction!

Can you identify the main sequence on this plot?

We also see that only the 1 solar mass star “finished’ stellar evolution, winding up as a cooling white dwarf at the end – the path it is following is essentially a line of constant radius, since the white dwarf does not contract as it cools (it is degenerate).

b#

Central evolution#

We want to plot the history of the evolution of the central conditions in the \(\log \rho\)-\(\log T\) plane.

Now we can plot the data

## in a similar way as above, plot the three stellar model history data, this time for the central density and 
# central temperature

Note: The 1 solar mass star makes a transition from following \(T \approx \rho^{1/3}\) to \(\rho = \mbox{constant}\) when degeneracy kicks in - no density dependence!

c#

Main sequence lifetime#

We can estimate the main sequence lifetime just by looking for when the core H is all consumed.

fig, ax = plt.subplots()
#ax.semilogx(# fill in here
../_images/77d023c5f3f358490c24b973a7db92376240ba6ad42d4fa22426f7d4f4491444.png

From this plot, we see that the main sequence lifetime of the 15 \(M_\odot\) star is \(\sim 6\times 10^6\) yr, for the 8 \(M_\odot\) star \(\sim 2\times 10^7\) yr, and for the 1 \(M_\odot\) star \(\sim 9\times 10^9\) yr. Do these values match the rough estimates from HKT 1.88?

Contents