Group Exercise 3.ipynb

Kernel: Python 3 (ipykernel)

Neutron star

Reference: Silbar & Reddy, Neutron Stars for Undergraduates (https://arxiv.org/abs/nucl-th/0309041)

Background

In this exercise, we solve for the mass-radius relation for neutron stars. We assume that the neutron star material is in hydrostatic equilibrium, i.e., assuming that the outward pressure and inward gravitational forces balance. In the Newtonian limit of weak gravitational fields and nonrelativistic velocities, the equation for hydrostatic equilibrium is

$$\frac{dp}{dr}= - \rho(r) \frac{G m(r)}{r^2} \qquad (1)$$

where $p(r)$ is the pressure, $\rho(r)$ is the density, $m(r)$ is the enclosed mass within a radius $r$, and $G$ is Newton's constant. Here and henceforth we work in cgs units, where $G = 6.673 \times 10^{-8} \; {\rm dyne} \; {\rm cm}^2/{\rm g}^2$ and the speed of light is $c = 3\times 10^{10} \; {\rm cm}/{\rm s}$.

We have also assumed spherical symmetry, such that all quantities depend on $r$ only.

The enclosed mass $m(r)$ can be calculated by integrating the density $\rho(r)$ over the volume of a sphere of radius $r$

$$m(r) = 4\pi \int_0^r dr^\prime r^{\prime 2} \rho(r^\prime) \; .$$

However, it is more helpful to take the derivative of this relation and relate $m(r)$ and $\rho(r)$ through a differential equation

$$\frac{dm}{dr} = 4\pi r^2 \rho(r) \; . \qquad (2)$$

In General Relativity, Eq. (1) is modified to have the following form

$$\frac{dp}{dr}= - \rho(r) \frac{G m(r)}{r^2} \left( 1 + \frac{p(r)}{\rho(r) c^2} \right) \left( 1 + \frac{4\pi r^3 p(r)}{m(r) c^2} \right) \left( 1 - \frac{2G m(r)}{r c^2} \right)^{-1}\qquad (3)$$

known as the Tolman-Oppenheimer-Volkov equation. This equation looks very similar to Eq. (1) except for the extra correction factors being multiplied on the right-hand side.

To solve for the structure of a neutron star, we need to solve for

$$p(r), \; m(r) , \; \rho(r)$$

by solving either Eqs. (1) and (2) in the Newtonian limit, or Eqs. (1) and (3) in full General Relativity. Since we have three unknown functions, we must also impose a third relation, known as the equation of state, which is a relationship between pressure and density.

Here, we adopt a quadratic equation of state following the reference given above (see Eq. (96) therein and converting to cgs units)

$$p(r) = K_0 \rho^2 \qquad (4)$$

where the constant $K_0 = 2.03 \times 10^5$ for $p$, $\rho$ given in cgs units.

In solving these equations, we need to impose boundary conditions. The strategy is to fix the boundary conditions at the origin, $r=0$, and solve for $p(r)$ and $m(r)$ for $r > 0$. The boundary conditions are

$$p(0) = p_0 , \quad m(0) = 0 \, .$$

Since $dp/dr$ is negative, the pressure will decrease with $r$. Where the pressure goes to zero, that defines the radius of the neutron star $R$, i.e.,

$$p(R) = 0$$

The total mass of the neutron star $M$ is the mass enclosed within $R$, i.e.,

$$m(R) = M$$

The values of $M$-$R$ that one obtains, as a function of the initial conditions, is the mass-radius relation.

Part A

Solve the system of coupled ODEs (1) and (2) corresponding to neutron stars in the Newtonian limit, taking a central density $\rho_0 = 10^{15} \; {\rm g}/{\rm cm}^3$.

Compute the total mass $M$ and radius $R$. Note: Values of $M$ will be in units of grams. Convert $M$ to units of solar masses. Also, $R$ will be in units of ${\rm cm}$; convert $R$ to units of ${\rm km}$.

Next, repeat the calculation above in the case of full General relativity, by solving the coupled ODEs (2) and (3), with the same initial conditions.

Make a log-log plot that shows $\rho(r)$ as a function of $r$ for the two cases (Newtonian and General Relativity).

Hints:

1. You can write your own ODE solver, or you can use solve_ivp (see documentation here). For solve_ivp, note that here the independent variable t is $r$, and the dependent variable y is $$y = \left(\begin{array}{c} p(r) \\ m(r) \end{array}\right)\, .$$

2. Since Eq. (1) is ill-defined exactly at $r=0$ (dividing by 0), impose your boundary condition at a small nonzero radius $r_0$. Take $$p(r_0) = p_0 = K_0 \rho_0^2 \, , \quad m(r_0) = m_0 = \frac{4\pi}{3} \rho_0 r_0^3$$ for your initial conditions, and solve for $p(r)$ and $m(r)$ for $r > r_0$.

3. You want to stop your ODE solver when the pressure drops to zero, i.e., you go from $r=r_0$ to $r=R$, where $R$ is defined by $p(R)=0$. This can be accomplished in solve_ivp using the events keyword, as follows. First, define a function event_function(t,y) that returns the pressure. Second, include the line of code

event_function.terminal = True

which tells solve_ivp to stop the integration if the pressure changes sign. Third, include the keyword events=event_func in solve_ivp.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Define constants
G = 6.673e-8 # dyne-cm^2/g^2
c = 3e10 # cm/s
M_sun = 1.989e33 # solar mass in grams
K0 = 2.03e5

# Your code here

Part B

Consider a range of central densty values defined by the following array:

rho0_list = np.logspace(14,16)

i.e., values of $\rho_0$ between $10^{14}$. For all values of $\rho_0$ given above, compute the values of $M$ (in solar masses) and $R$ (in ${\rm km}$) for both cases of Newtonian and General Relativistic neutron stars. Make a plot of $M$ (on the y-axis) and $R$ (on the x-axis), labeling the two cases. Use a log-scale for the y-axis and a linear scale for the x-axis.

Determine the maximum neutron star mass $M$ according to your calculations in the case of full General Relativity.

Hint: You should find that the Newton radius $R$ is independent of $M$, and that there is no maximum mass.

# Your code here