In-Class Assignment 5

In-Class Assignment 5#

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

EOS regimes#

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

import numpy as np
import matplotlib.pyplot as plt

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

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}\]

The model files are again:

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

# compute the constant B as defined above here

# create a function that takes in a mu_e, and rho array
# then, computes x_F, then returns T for E_F = kT

# create an array for rho here to evaluate this line
# rho = np.logspace(-5, 10,###

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

# create a function here that takes in rho and mu and returns 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

Plotting the boundaries#

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.

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

#ax.plot(rho, ###

Plotting the boundaries + MESA#

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

To plot the results all on the same plot, you will need to match the size of rho to the number of points for MESA.

One strategy would be to

  1. load the MESA data using pandas,

  2. find the number of cells,

  3. then create a new rho array for use with the above functions that matches the size of the history central \(\rho\) array.

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

#ax.plot(rho, ###
Contents