In-Class Assignment 9

In-Class Assignment 9#

Stellar Energy Sources#

Learning Objectives#

explore thermonuclear burning in various stellar models
identify the dominant burning sources in stars during H-burning
compare crossover times in dominant burning categories
visualizing main-sequence stellar models using tulips

To produce plots in the optional exercises we will use mesaplot and Tulips - A Tool for Understanding the Lives, Interiors, and Physics of Stars and the Paper here. Install each by performing the following commands to install them using pip:

mesaPlot

Note

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

pip install git+https://github.com/rjfarmer/mesaplot

Tulips Installation

Note

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

pip install astro-tulips

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# only uncomment if you were able to successfully install
#import tulips

# note the capital `P`
#import mesaPlot as mp

### only try below if imports above to not work
##import sys 
# you may have to run the below, try importing without first
# this will be different for you, find its location using `pip show mesaplot`
## sys.path.append('/opt/anaconda3/lib/python3.11/site-packages/mesaPlot') 

a. Thermonuclear Burning Accross Mass Ranges#

Download the following model files locally:

\(1 M_\odot\): 1m_h_dep_history.data;
\(3 M_\odot\): 3m_h_dep_history.data;
\(20 M_\odot\): 20m_he_dep_history.data;

Then, using the \(1 M_\odot\) and \(3M_\odot\) model MESA history data:

Plot the log luminosity (\(L_{\rm{nuc}}\)) for the pp and cno burning categories as a function of the model_number (x-axis) to determine which dominates in each model.

Note: The variables pp and cno are \(\rm{log}_{10}(L_{\rm{nuc,pp}})\) etc.

To help guide follow the temporal evolution, also plot the center_h1 as a function of model number for both models on the same plot and label.

Set your ylim to (-2,) - when you leave the upper limit unset it will auto set the ymax.

Identify the approximate model_number where the dominant burning crossover occurs for the \(3M_{\odot}\) model and plot a vertical line (plt.axvline) denoting this on the same plot.
On a new plot, plot \(\textrm{log}~(L_{\rm{nuc}}/\rho X^2)\) as a function of the central \(T_{6,c}\) (\(T_{c}/T^{6}\rm{K}\)).

Hint: Set xlim to (5,30) and ylim to (-6,6)

On this plot, plot the approximate central \(T_{6}\) where the CNO burning category begins to dominate and label using plt.axvline.
Lastly, print the central \(T_{6}\) at the model number for the \(3M_{\odot}\) model found for the crossover in part 3. Do they match?
In a few sentences, describe the dominant burning category for the 1M model and 3M models based off these two plots.

Hint: Compare with HKT Figure 6.10!

## load MESA data here
#one_m_ms_history = pd.read_csv('###',sep=r'\s+',header=4)
#three_m_ms_history = pd.read_csv('####
#list(one_m_ms_history)

# load one m burning info here
#one_m_ms_model_number = one_m_ms_history['###  # model number
#one_m_ms_pp = ### L_nuc due to pp chains
#one_m_ms_cno = # L_nuc due to cno chains
#one_m_ms_center_h1 = one_m_ms_history['##

# load three m burning info here

## 1-3 result here
#plt.title('1$M_\odot$ and 3$M_\odot$ MESA Models During H-burning')

#plt.plot(##,
         ##,label=r'$1M_\odot - pp-chains$')

#plt.plot(##
         ##,label=r'$1M_\odot - CNO-chains$')

#plt.plot(##
         ##,label=r'$1M_\odot - ^{1}\rm{H}$')


#plt.axvline(###,lw=1,color='goldenrod',ls='-.') ## crossover model number

#plt.ylim(


#plt.legend()
#plt.xlabel(r'$\rm{Model \ Number}$')
#plt.ylabel(r'$\rm{log}_{10} (L_{\rm{nuc}}) \ (\rm{erg \ g}^{-1} \ \rm{s}^{-1})$')

# load central rho and T here
#one_m_ms_center_T = 10**(### central T
#three_m_ms_center_T = 10**(### central T

#one_m_ms_center_rho = 10**(### central rho
#three_m_ms_center_rho = 10**(### central rho

## 4 result here
#plt.title('1$M_\odot$ and 3$M_\odot$ MESA Models During H-burning')

#plt.plot(##/1e6,
#         np.log10(##),label=r'$1M_\odot - pp$')

#plt.plot(##/1e6,
#         np.log10(##),label=r'$1M_\odot - CNO$')

# approximate present location of the Sun
#plt.text(15, -5.5, "⊙")

#plt.axvline(##,color='dodgerblue',lw=1,ls='-.')

#plt.ylim(
#plt.xlim(

#plt.legend()
#plt.ylabel(r'$\rm{log}~L_{\rm{nuc}}/\rho X^2$')
#plt.xlabel(r'$T_{6,\rm{central}}$')

# compute approximate crossover temperature here
#crossover_index = np.where(three_m_ms_history['model_number']==XXX)[0]
#print(crossover_index,three_m_ms_center_T[crossover_index]/1e6)

In a few sentences, describe the dominant burning category for the 1M model and 3M models based off these two plots. Does the approximate crossover point match your estimates from Plot 1 and 2?

Your thoughtful answer here.

b. Thermonuclear Burning In a Massive Star#

Using the \(20M_\odot\) model MESA history data (20m_he_dep_history.data):

Plot the log luminosity (\(L_{\rm{nuc}}\)) for the pp and cno burning categories as a function of the model_number (x-axis) to determine which dominates in each model.

Set your xlim to (370,600) to focus on the main-sequency

To help guide follow the temporal evolution, also plot the center_h1 as a function of model number for both models on the same plot and label.

Are there any crossover points observed in this model between dominant burning categories?

Plot \(\textrm{log}~(L_{\rm{nuc}}/\rho X^2)\) as a function of the central \(T_{6,c}\) (\(T_{c}/T^{6}\rm{K}\)) for both categories.

Hint: set xlim to (30,100)

Conclude the dominant burning source on the main-sequence and state if there are any crossover points in a few words.

# load 20m data
#twenty_m_ms_history = pd.read_csv('###',sep=r'\s+',header=4)
#twenty_m_ms_model_number = twenty_m_ms_history['### model number
#twenty_m_ms_pp = ### L_nuc due to pp chains

## 1-3 result here
#plt.title('20$M_\odot$ MESA Models During H-burning')

#plt.plot(,label=r'$20M_\odot - pp-chains$')

#plt.plot(,label=r'$20M_\odot - CNO-chains$')

#plt.plot(,label=r'$20M_\odot - ^{1}\rm{H}$')

#plt.ylim(
#plt.xlim(

#plt.legend()
#plt.xlabel(r'$\rm{Model \ Number}$')
#plt.ylabel(r'$L_{\rm{nuc}} \ (\rm{erg \ g}^{-1} \ \rm{s}^{-1})$')

#twenth_m_ms_center_T = 10**(twenty_m_ms_history['log_center_T']) # central T
#twenty_m_ms_center_rho = 10**(twenty_m_ms_history['log_center_Rho']) # central rho

## 4 result here
#plt.title('20$M_\odot$ MESA Models During H-burning')

#plt.plot(#
#         label=r'$20M_\odot - pp$')

#plt.plot(,
#         label=r'$20M_\odot - CNO$')


#plt.xlim(30,

#plt.legend()
#plt.ylabel(r'$\rm{log}~L_{\rm{nuc}}/\rho X^2$')
#plt.xlabel(r'$T_{6,\rm{central}}$')

c. (Optional!) Visualizing the Main-Sequence using Tulips#

Import tulips and mesaPlot here if you did not import above and then:

Download this \(1M_{\odot}\) model evolved from MS to this H-shell burning phase:

\(1 M_\odot\): 1m_thermohaline_history.data;

Create a folder in your data directory or anywhere really, named tulips and copy your 1m_thermohaline_history.data into the directory and rename it as history.data.
Then set your directory path for tulips. It will look in this directory for a file called history.data.
Create your data object for analysis using mesaPlot:

m1 = mp.MESA()

Load the history data using this object:

m1.loadHistory(f=`YOUR_PATH_NAME`)

Now, using tulips to produce an energy_and_mixing plot using your history data. It will create a series of these plots based on your defined step size N.
Lastly, open and view the plot, can you identify where the model is convective or radiative on the main-sequence?
(Bonus), Can you see when are where the model experiences thermohaline/double diffusive mixing?

## 1- 4 here
#SINGLE_M1_DIR = "./data/tulips"  ## the location of your renamed `history.data` file

#m1 = mp.MESA()
#m1.loadHistory(f=SINGLE_M1_DIR) # load the history data file for plotting
#N = 1 # series index

# compute and plot the viz
#tulips.energy_and_mixing(m1, time_ind=(0,-1, N), fps=10, fig_size=(8,6),
#                         show_total_mass=True,show_hrd_ticks_and_labels=True, show_mix=True, show_mix_legend=True, 
#                         output_fname="energy_and_mixing")
#plt.show()

# open and watch the visualization here
#from IPython.display import Video

#Video("energy_and_mixing.mp4", embed=True, width=700, height=600)

Your thoughtful response to 6 (and 7 possibly) here.