In-Class Assignment 3

In-Class Assignment 3#

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

Work in groups of 2

Learning Objectives#

  • explore EOS regimes

  • label various EOS regimes

  • plot stellar models on EOS regimes

  • compare density profiles esimtates to stellar models

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

Physical constants, all in CGS

m_e = 9.11e-28
m_u = 1.66e-24
c = 3.e10
k = 1.38e-16
h = 6.63e-27
a = 5.67e-15
e = 4.8e-10 # esu

# crystallization
Gamma_C = 170

EOS regimes#

We can look at the \(\rho\)-\(T\) plane to understand where each of the components of the EOS dominate.

Some composition info. We’ll specify

\[\mu_e = \left [ \sum_k \frac{Z_k X_k}{A_k} \right ]^{-1}\]
\[\mu_I = \left [ \sum_k \frac{X_k}{A_k} \right ]^{-1}\]

and \(\mu\) such that

\[\frac{1}{\mu} = \frac{1}{\mu_e} + \frac{1}{\mu_I}\]
mu_e = 2
mu_I = 1.0

mu = 1.0 / (1.0/mu_e + 1.0/mu_I)

Z = 1.0

Degeneracy vs. Ideal gas#

We find the (dimensionless) Fermi momentum (\(x_F = p_F / (mc)\)) via

\[n_e = \frac{\rho}{\mu_e m_u} = \int_0^{x_F} x^2 dx\]

where we assume complete degeneracy (so the distribution function is just a step-function). This gives

\[x_F = \left ( \frac{\rho}{\mu_e B} \right )^{1/3}\]

where the constant \(B\) is:

\[B = \frac{8\pi}{3} m_u \left (\frac{m_e c}{h} \right )^3\]

We then define the transition to degeracy by comparing to the thermal energy, \(k T\)

\[E_F = m_e c^2 \left ( \sqrt{1 + x_F^2} - 1 \right ) \sim k T\]

Here’s a function that computes this line

# Fermi momentum constant (from number density integral)
B = # 

def deg_ideal(rho, mu_e=2.0):
    """temperature on the line separating degeneracy and ideal gas as a function
 of rho"""

    # dimensionless Fermi momentum
    x_F = (rho / mu_e / B)**(1./3.)

    # set E_F = k T
    T = (1./k) * m_e * c**2 * (np.sqrt(1 + x_F**2) - 1)

    return T

  Cell In[4], line 2
    B = #
        ^
SyntaxError: invalid syntax

Ideal gas vs. radiation#

A simple way to determine whether an ideal case or radiation dominates is to set their energy densities equal. At the temperatures where we expect this to occur, we will assume that the electrons are not degenerate, so we will use the full \(\mu\) for the composition in the ideal gas instead of \(\mu_I\).

Our line is:

\[\frac{\rho k T}{\mu m_u} \sim \frac{1}{3}a T^4\]

Here’s a function that computes this

def rad_ideal(rho, mu=0.5):
    """temperature on the line separating radiation and ideal gas as a function 
of rho"""

    # T = complete here
    return T

Finally, at low temperatures and high densities, the ions can crystalize. The comparison where is the thermal energy to the Coulomb energy:

\[\Gamma = \frac{kT}{Z^2 e^2 / a}\]

where we take \(\Gamma = 170\) as the point where crystalization sets in.

The separation between ions, \(a\) is found via:

\[V = n_I^{-1} = \frac{\mu_I m_u}{\rho} = \frac{4\pi}{3} a^3\]

Putting this together, the expression for the boundary of crystalization is computed below

def crystallization(rho):
    """temperature where crystallization of ions sets in"""

    T = ((rho / mu_I) * (Z**2 * e**2 / (Gamma_C * k))**3 * (4.0 * np.pi / 3.0 / 
m_u))**(1./3.)
    return T

a. Plotting the boundaries#

Download the following model files (or reuse them from last time!):

  • \(1 M_\odot\) history data M1_default_trimmed_history.data

  • \(1 M_\odot\) - main sequence profile data M1_default_profile8.data

    1. Plot and label the three lines defined above. Label the radiation dominated, ideal gas, the electron degenerate regime (including a line denoting the density from non-relativistic to relativistic regimes), and lastly, label the crystallization regime.

#rho = np.logspace(-5, 10, 2000)

#fig = plt.figure()
#ax = fig.add_subplot(111)

#ax.plot(rho,##

# annotate NR vs ER degeneracy line here
#ax.vlines(*


#ax.text(1.e-2, 1.e9, 

# label the location of the Sun here
#ax.text(150,#)

#ax.set_xscale("log")
#ax.set_yscale("log")

#ax.set_xlabel(r"$\rho$ (g/cc)")
#ax.set_ylabel("$T$ (K)")
#ax.grid(ls=":")

b. Plotting the boundaries + MESA#

Here, produce the same plot but overlay the results of the 1M MESA stellar model history data from ICA1.

  1. load the MESA data using pandas,

  2. copy the above labelled plot below and plot the 1Msun MESA model on the same plot

Descirbe in a few sentences to trajectory of the 1M stellar model in the diagram.

#one_m_sun_history = pd.read_csv('data/M1_default_trimmed_history.data',sep=r'\s+',header=4)
#one_m_sun.columns

# obtain density and temperature data here

# combined plot here

c. Stellar density profile#

In ICA2, we made an a stellar density profile of the form:

\[\rho(r) = \rho_c \left [ 1 - \left (\frac{r}{R_\star} \right )^3 \right ]\]

where \(\rho_c\) is the central density.

Below we will test if this profile is an appropriate for a 1 solar mass model on the main-sequence.

  1. load the MESA profile data using pandas M1_default_profile8.data this is the stellar structure profile on the main-sequence

  2. next, produce an equation density profile above using the central density from the profile data as well as the max radius for \(\rho_c\) and \(R_\star\).

  3. Plot both on the same plot then answer the following below.

Descirbe in a few sentences if the analytical profile is appropriate for this stellar model?

Hint: Such analytical or polytropic EOS can be great approximations for stellar structure but typically not for the entire star except for in the cases of fully convective or white dwarf stars - not a star on the MS!

#load `M1_default_profile8.data` profile data here with pandas

# set the radius and density variables here, we only have log rho so youll need to convert!

# compute max / central density here and Rstar or maximum radius here

# use these above rho_central and Rstar values to define our rho(r) profile from the analytical equation

# finally, plot both on the same plot

Descirbe in a few sentences if the analytical profile is appropriate for this stellar model?

Contents