In-Class Assignment 18

In-Class Assignment 18#

In-Class Only, Not Submitted for Credit, Credit: Michael Zingale.

Learning Objectives#

  • Develop a qualitative understanding of the M-R relationship for neutron stars

  • Identify the impact of the compressibility modulus \(K\) on the M-R relation and the maximum mass

import numpy as np
import matplotlib.pyplot as plt

Here’s a general integration class that implemented 4th order Runge-Kutta. It takes a rhs function of the form:

rhs(r, Y, rho_func)

that returns the derivatives \(d{\bf Y}/dr\) for each of the variables in the array \({\bf Y}\). Here, rho_func is a function rho_func(pres) that takes the pressure and returns the density according to our EOS.

class Integrator:
    """ this is a simple RK4 integrator with a uniform step size
    that can accept a list of stopping criteria to halt the integration.

    There is one list element for each variable.  An entry of "None"
    means that we have no stopping criteria for that variable.  A tuple
    of the form ("LT", val) or ("GT", val), means that we stop if
    y[n] < val (or y[n] > val)
    """

    def __init__(self, rhs, *, stop_cond=None, rhs_args=None):
        self.rhs = rhs
        self.stop_cond = stop_cond
        if rhs_args is None:
            rhs_args = []
        self.args = rhs_args

    def rk4_step(self, x, y, dx):

        dydx1 = self.rhs(x, y, *self.args)
        y1 = y + 0.5*dx*dydx1
        if self.check_stop_cond(y1) < 0: return None

        dydx2 = self.rhs(x + 0.5*dx, y1, *self.args)
        y2 = y + 0.5*dx*dydx2
        if self.check_stop_cond(y2) < 0: return None

        dydx3 = self.rhs(x + 0.5*dx, y2, *self.args)
        y3 = y + dx*dydx3
        if self.check_stop_cond(y3) < 0: return None

        dydx4 = self.rhs(x + dx, y3, *self.args)

        y_final = y + (dx/6.0)*(dydx1 + 2.0*dydx2 + 2.0*dydx3 + dydx4)
        if self.check_stop_cond(y_final) < 0: return None

        return y_final

    def check_stop_cond(self, y):
        # check our stopping conditions
        if not self.stop_cond is None:
            for n, cond in enumerate(self.stop_cond):
                if cond is None: continue
                c, v = cond
                if c == "LT":
                    if y[n] < v: return -1
                elif c == "GT":
                    if ynew[n] > v: return -1
                else:
                    sys.exit("invalid condition")

        return 0

    def integrate(self, initial_conditions, dx, xmax, tol=1.e-8):
        ic = np.array(initial_conditions)
        N = len(ic)
        sol = {}
        for n in range(N):
            sol[f"y{n}"] = [ic[n]]

        x = 0.0
        sol["x"] = [x]

        y = initial_conditions
        
        while (x < xmax and dx > tol * xmax):
            # advance for a step
            ynew = self.rk4_step(x, y, dx)

            if ynew is None:
                # we hit our stopping condition
                # cut the stepsize so we can try to get closer
                dx /= 2
                continue

            # store the solution
            x += dx
            sol["x"].append(x)
            for n in range(N):
                sol[f"y{n}"].append(ynew[n])

            y = ynew

        return sol

Some physical constants (CGS)

G = 6.67e-8
c = 3.e10
h = 6.63e-27
m_u = 1.66e-24

Here’s the righthand side function for the TOV equation

def rhs_tov(r, y, rho_func):
    m = y[0]
    P = y[1]

    rho = rho_func(P)
    
    dmdr = 4.0 * np.pi * r**2 * rho
        
    if r == 0.0:
        dPdr = 0.0
    else:
        rho_term = rho + P / c**2
        m_term = m + 4.0 * np.pi * r**3 * P / c**2
        metric_term = (1.0 - 2.0 * G * m / (r * c**2))
        dPdr = -(G / r**2) * rho_term * m_term / metric_term

    return np.array([dmdr, dPdr])

Here’s the righthand side function for Newtonian HSE

def rhs_newton(r, y, rho_func):
    m = y[0]
    P = y[1]

    rho = rho_func(P)
    
    dmdr = 4.0 * np.pi * r**2 * rho
    
    if r == 0.0:
        dPdr = 0.0
    else:
        dPdr = -G * m / r**2 * rho
        
    return np.array([dmdr, dPdr])

Degenerate neutrons#

For zero-T neutron degeneracy, the number density of neutrons is just the integral of the distribution function, \(n(p)\), over all momenta:

\[n = \int n(p) d^3 p = 4\pi \int_0^\infty n(p) p^2 dp\]

but at zero temperature the distribution function is just a step function:

\[\begin{split}n(p) = \frac{2}{h^3}\left \{ \begin{array}{cc} 1 & \mbox{if } p < p_F \\ 0 & \mbox{otherwise} \end{array}\right .\end{split}\]

This gives:

\[n = \frac{\rho}{m_u} = \frac{8\pi}{3h^3} p_F^3 = \frac{8\pi}{3} \left (\frac{m_u c}{h} \right )^3 x_F^3\]

where \(x_F = p_F / (m_u c)\).

We can then compute the pressure as:

\[P = \frac{1}{3} \int n(p) v p d^3 p\]

where we now need to use the general expression for velocity,

\[v = \frac{p}{m_u} \left [ 1 + \left ( \frac{p}{m_u c} \right)^2 \right ]^{-1/2}\]

This gives:

\[P = \frac{\pi}{3} \left ( \frac{m_u c}{h} \right )^3 m_u c^2 f(x)\]

with

\[f(x) = x (2x^2 - 3) (1 + x^2)^{1/2} + 3 \sinh^{-1}(x)\]
from scipy.optimize import brentq

def pres_deg(rho):
    # compute the Fermi momentum x = p / mc
    B = 8 * np.pi * m_u / 3.0 * (m_u * c / h)**3
    x = np.cbrt(rho / B)
    
    f = x * (2 * x**2 - 3) * np.sqrt(1 + x**2) + 3 * np.arcsinh(x)
    A = (np.pi / 3) * (m_u * c / h)**3 * m_u * c**2
    
    return A * f
    
def rho_deg(pres):
    
    A = (np.pi / 3) * (m_u * c / h)**3 * m_u * c**2
    f = pres / A
    
    # find the x that gives us this f
    x = brentq(lambda x: f - (x * (2 * x**2 - 3) * np.sqrt(1 + x**2) + 3 * np.arcsinh(x)), 0, 10)
    
    B = 8 * np.pi * m_u / 3.0 * (m_u * c / h)**3
    return B * x**3

A polytropic nuclear EOS#

A simple polytopic approximation to a nuclear equation of state (the SLy4 EOS) is:

\[P = K \left ( \frac{\rho}{\rho_0} \right )^\gamma\]

If we choose \(K = 3.5 \mathrm{MeV~fm^{-3}}\) and \(\rho_0 = 150 \mathrm{MeV~fm^{-3}}\), we approximate the behavior of the SLy4 nuclear EOS (thanks to Tianqi Zhao for these approximations) with \(\gamma = 2.5\).

In CGS units, using \(\gamma = 2.5\), we can take \(K = 4.80\times 10^{-3}~\mathrm{erg~cm^{-3}}\) and write this as:

\[P = K \rho^\gamma\]
def pres_nuc(rho):
    gamma = 2.5
    K = 4.8e-3
    
    return K * rho**gamma

def rho_nuc(pres):
    gamma = 2.5
    K = 4.8e-3
    
    return (pres / K)**(1.0 / gamma)

Notebook begins here#

a. - The “Chandrasekhar mass” of a Neutron Star#

The Chandrasekhar mass, the maximum mass that electron degenerate objects (of fixed \(\mu_{e}\)) can have before becoming gravitationally unstable and collapsing, is given by HKT 3.67:

\[M_\mathrm{Ch} = 1.44 \left ( \frac{2}{\mu} \right )^2 M_\odot\]

For electrons, the \(\mu\) was \(\mu_e\), as in white dwarfs. But, let’s use this to compute the maximum mass of a neutron star, where now \(\mu = 1\) for neutrons.

  1. Compute the maximum mass for a NS using the equation above.

# a result here

b. - Case I: Newtonian gravity, neutron degeneracy#

Now, let’s look at the M-R relation for a Neutron star where we will use our usual HSE, but with pressure derived from \(n\) degeneracy pressure.

  1. Evaluate the cell below for a M-R relation using our usual non-GR HSE, but with \(n\) degeneracy pressure.

  2. Plot the resulting M-R and label the maximum allowed NS mass from this result.

Compare with the result of a. Do they match?

We’ll loop over the central density and integrate the neutron star structure and store the mass and radius corresponding to each central density.

ns1 = Integrator(rhs_newton, rhs_args=[rho_deg],
                 stop_cond=[("LT", 0.0), ("LT", 0.0)])

R_max = 5.e6
N = 1000

mass_ns1 = []
radius_ns1 = []

for rhoc in np.logspace(np.log10(5.e14), np.log10(5.e18), 100):
    
    Pc = pres_deg(rhoc)
    
    sol = ns1.integrate([0.0, Pc], R_max/N, R_max)
    
    mass_ns1.append(sol["y0"][-1]/2.e33)
    radius_ns1.append(sol["x"][-1]/1.e5)

# a result here

We see that we are approaching that maximum mass as \(R \rightarrow 0\), as expected.

c. - Case II: GR gravity (TOV) + neutron degeneracy#

  1. Evalue the cell below, plot the resulting R vs. M plot.

  2. Compute the maximum allowed mass under these assumptions.

Is this mass less or more than the previous? Qualitatively, what lead to the change in the maximum mass from looking at the R-M profile?

ns2 = Integrator(rhs_tov, rhs_args=[rho_deg],
                 stop_cond=[("LT", 0.0), ("LT", 0.0)])

R_max = 5.e6
N = 1000

# my vector is [m, p]
mass_ns2 = []
radius_ns2 = []
for rhoc in np.logspace(np.log10(5.e14), np.log10(5.e16), 100):

    Pc = pres_deg(rhoc)

    sol = ns2.integrate([0.0, Pc], R_max/N, R_max)

    mass_ns2.append(sol["y0"][-1]/2.e33)
    radius_ns2.append(sol["x"][-1]/1.e5)

# c results here

Now, adding GR, we see that there is a maximum mass for the neutron star, but it is quite low.

d - Case III: GR gravity (TOV) + nuclear EOS#

Next, let’s use a more accuracte equation of state using an analtical form for the nuclear EoS defined above.

  1. Evaluate the cell below, plot the result.

  2. Compute the maximum NS mass allowed under these assumptions.

Qualitatively describe the new M-R profile when changing the EoS compared to the previous profile in c.

ns3 = Integrator(rhs_tov, rhs_args=[rho_nuc],
                 stop_cond=[("LT", 0.0), ("LT", 0.0)])

R_max = 5.e6
N = 1000

# my vector is [m, p]
mass_ns3 = []
radius_ns3 = []
for rhoc in np.logspace(np.log10(5.e14), np.log10(5.e16), 100):

    Pc = pres_nuc(rhoc)

    sol = ns3.integrate([0.0, Pc], R_max/N, R_max)

    mass_ns3.append(sol["y0"][-1]/2.e33)
    radius_ns3.append(sol["x"][-1]/1.e5)

# d results here

In this final case, we get a much more reasonable maximum neutron star mass. We also see that for some values of radius there are multiple masses (the M-R relation is double-valued) – this is okay. However, at the very highest central densities, we get a smaller mass than the maximum at a smaller radius – that is not physical, and is dynamically unstable.

e - The impact of the compressibility modulus#

  1. Plot the results of c and d on the same plot.

  2. Looking only at the respective maximum mass for each profile, label which relation suggests a “softer” EoS and which suggests a “stiffer” EoS.

Qualitatively describe the difference between the “softer” and “stiffer” EoS, what that means for the incompressibility of the matter and the maximum mass.

Contents