MESA

Introduction to MESA

MESA

Modules for Experiments in Stellar Astrophysics
modules can be used separately in other codes
modules: rates (nuclear reaction rates), eos (equation of state), mlt (mixing-length theory), kap (opacities), …
star module is a full stellar evolution code

Equations of stellar structure

\begin{align} \PD{r}{m} &= \frac{1}{4 \pi r^2 \rho}\\ \PD{p}{m} &= - \frac{G m}{4 \pi r^2}\\ \PD{l}{m} &= \epsilon_\text{n} - \epsilon_\nu - c_p \PD{T}{t} + \frac{\delta}{\rho} \PD{p}{t}\\ \PD{T}{m} &= - \frac{G m T}{4 \pi r^4 p} \nabla\\ \PD{X_i}{t} &= \frac{m_i}{\rho} \left(\sum_j r_{ji} - \sum_k r_{ik} \right) \end{align}

Running the code

./mk: compile (only needed once per directory)
./rn: run code from initial condition
Ctrl+c: interrupt code at any time
./re x100: restart from a checkpoint
(look in photos directory for number)

Control using inlists

plain text files control most features of MESA
using the namelist feature of Fortran 90
unset options have a (reasonable) default value
everything after ! is a comment (no effect on MESA)
sections between &section_name and /
each entry has a fixed data type:
- logical: .true. or .false.
- string: 'a string'
- integer: 42
- floating-point number: 2.99792458d8

Some important inlist entries

pgstar_flag: show plot windows during simulation
initial_mass: zero age main sequence (ZAMS) mass (in $M_\odot$)
initial_z: initial metallicity
profile_interval: frequency of output to profile files
(similar for history, terminal, photo)
profile_columns_file: name of file with list of columns to include in profile
TRho_Profile_win_flag: show window with $T$–$\rho$ profile
(also Hr_win_flag, Abundance_win_flag, …)
complete list in $MESA_DIR/star/defaults/*.defaults

Some stopping criteria

star_H_mass_min_limit: stop when total hydrogen mass (in $M_\odot$) drops below this value

stop when central $^4\mathrm{He}$ mass fraction drops below this value

xa_central_lower_limit_species(1) = 'he4'
xa_central_lower_limit(1) = 0.25

stop when $L_\mathrm{nuc}/L_\mathrm{total}>0.99$

Lnuc_div_L_zams_limit = 0.99d0
stop_near_zams = .true.

Useful diagrams

HR diagram

Abundance profile

$T-\rho\,$ diagram

MESA output files

plain text files in the LOGS directory
frequency of output controlled in inlist files
units usually in CGS or solar units
(check comments in columns file)
logarithmic values use $\log_{10}$
quick viewing on console:
```
less -S LOGS/filename.data
```
(-S prevents line wrapping)
move with arrow keys, exit with q

History file

LOGS/history.data
records development of scalar quantities (e.g., total mass or luminosity) for many/all steps

Profile file

full stellar profile at a specific time
lines correspond to mass shells in the code
number of mass shells changes between steps (adaptive mesh refinement)
ordered from the surface to the core
make sure to select the right columns before a long simulation

Intensive vs. extensive quantities

meaning of a value at index $k$ depends on the type of quantity
intensive quantity: independent of size of system
(e.g., $T$, $P$, $\rho$, $\epsilon_\mathrm{nuc}$)
extensive quantity: depends on size of system
(e.g. $m$, $r$, $L$)

Cells in MESA

intensive quantities are given at the cell centre
extensive quantities are given at the cell interface
values can be translated from centre to interface using: \begin{align} \alpha_k &= \frac{dm_{k-1}}{dm_{k-1}+dm_k}\\ P^\mathrm{edge}_k &= \alpha_k P^\mathrm{centre}_k + (1-\alpha_k) P^\mathrm{centre}_{k-1} \end{align}
distance between cell centres: $$\overline{dm}_k = \frac{dm_k + dm_{k-1}}{2}$$

figure explaining the cell indices in MESA

Paxton et al. (2011)

Reading and plotting data

Python

general-purpose programming language
efficient handling of numerical arrays using numpy
plotting using matplotlib
free software
interactive shell ipython

ipython --pylab

Simple plotting

x = linspace(0,5)
plot(x, exp(x), color='green')
xlabel('x')
ylabel('exp(x)')
yscale('log')
title('exponential')

save you figure with the save button (or savefig(…))
many formats supported (pdf, png, eps, svg, …)

Arrays in Python

efficient array implementation using numpy package
indexing starts at 0
operators are always applied element-wise
index operator [] can be used to extract sub-arrays
Try these examples:
```
x = arange(10)
x[5:7]
x[1:]
x[:-1]
x[1::2]
x[::-1]
```

Reading MESA profiles

make sure mesa.py is in the current directory

read the profile into an object

import mesa
prof = mesa.profile('LOGS/profile22.data')
plot(prof.mass, prof.pressure)

plot mass coordinate against pressure
```
plot(prof.mass, prof.pressure)
```

all columns from the file are accessible as attributes

prof.columns # list of columns
prof.he4 # helium 4 mass fraction (array)
prof.star_age # stellar age in years (scalar)

Reading MESA profiles in Matlab

make sure mesa_profile.m is in the current directory
```
prof = mesa_profile('LOGS/profile22.data')
```
rest of the interface is pretty much identical to the Python function

First assignment (due 16/10/2017)

Hydrostatic equilibrium (HSE)

exercise sheet: printed/Blackboard/git
one-page written report: marked, but does not influence your grade
first 4 groups give a presentation (15 min each, incl. questions)
part of your final grade
some ideas for the content:
- introduce HSE
- discuss all tasks
- show your HSE plot (Task c)
- technical details of analysing and plotting data
- all groups have a different star: introduce it with a meaningful plot

Finite differences

How do we find the derivative of a function given at fixed values $x_i$?

forward difference: $f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}$
backward difference: $f'(x_{i+1}) = \frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}$
central difference: $f'(\frac{1}{2}(x_{i+1}+x_i)) = \frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}$