from astropy import constants as const
from astropy import units as u
import matplotlib.pyplot as plt
import numba
import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp

plt.rcParams.update({
    'text.usetex': True,
    'font.sans-serif': ['Helvetica']
})

Introduction¶

Sec. II A: Newtonian structure equation
Sec. II B: General relativistic corrections
Sec. III: White dwarf stars
Sec. IV: Pure neutron stars

The Tolman-Oppenheimer-Volkov equation¶

A Newtonian formulation¶

$$\begin{cases} \dfrac{dp}{dr} &= -\dfrac{G\rho(r)M(r)}{r^2} \\ \dfrac{dM(r)}{dr} &= 4\pi r^2 \rho(r) \\ \displaystyle M(r) &= \displaystyle 4\pi\int_0^r \rho(s) s^2 \ ds \end{cases} $$

Sometimes we may want to express in terms of energy density: $$ \varepsilon(r) = \rho(r) c^2 $$

General relativistic corrections¶

$$ \frac{dp}{dr} = -\frac{G\varepsilon(r)M(r)}{r^2c^2}\left[1 + \frac{p(r)}{\varepsilon(r)}\right]\left[1 + \frac{4\pi r^3p(r)}{M(r)c^2}\right]\left[1 - \frac{2GM(r)}{rc^2}\right]^{-1} $$

White dwarf stars¶

A few facts¶

Newtonian structure equation is sufficient
Rather light mass $\approx 1.4M_\odot$
Small radius $\approx 10^4$ km (Earth already has a radius of $0.64\times10^4$ km)
Formation: Burned up all nuclear fuel $\rightarrow$ Electrons fall to lowest orbits $\rightarrow$ Pauli principle and degeneracy pressure stabilize against gravitational force

Fermi gas model for electrons¶

$$ dn = \frac{4\pi k^2 dk}{(2\pi\hbar)^3} \hspace{10pt}\Longrightarrow\hspace{10pt} n = \frac{k_F^3}{3\pi^2\hbar^3} $$

The symbol $k$ is momentum, *not* wavenumber. So instead of $p = \hbar k$, here we should write $k = \hbar k_\mathrm{wave}$, which is why $n$ has a unit of "per unit volume".

Assume each electron is neutralized by a proton, and accompanied by one or a few neutrons, then the mass density: $$ \rho = nm_N \frac{A}{Z} \tag{8}$$ Hence, $$ k_F = \hbar\left(\frac{3\pi^2\rho}{m_N}\frac{Z}{A}\right)^{1/3} \tag{9} $$

A convenient unit for the Fermi momentum of electrons: $$x \equiv \frac{k_F}{m_ec}$$

A group of constants that will keep appearing: $$ \varepsilon_0 \equiv \frac{m_e^4c^5}{\pi^2\hbar^3} $$

epsilon_0 = (const.m_e**4 * const.c**5) / (np.pi**2 * const.hbar**3)
epsilon_0.to('erg/cm^3')

To work out the E.O.S., $p = p(\varepsilon_\mathrm{elec})$, we need energy density and pressure: $$\varepsilon_\mathrm{elec}(x) = \varepsilon_0\cdot\int_0^x (1+u^2)^{1/2}u^2 \ du \tag{10}$$ $$p(x) = \varepsilon_0\cdot\frac{1}{3}\int_0^x (1+u^2)^{-1/2}u^4 \ du \tag{13}$$ Notice both quantities have the same unit.

Mathematica can give us the integrals: $$\varepsilon_\mathrm{elec}(x) = \varepsilon_0\cdot\frac{1}{8}\left[(2x^3 + x)\sqrt{1+x^2} - \sinh^{-1}(x)\right] \tag{10}$$ $$p(x) = \varepsilon_0\cdot\frac{1}{24}\left[(2x^3 - 3x)\sqrt{1+x^2} + 3\sinh^{-1}(x)\right] \tag{13}$$

def energy_density(x, include_constant=True):
    quant = (1/8) * ((2 * x**3 + x) * np.sqrt(1 + x**2) - np.arcsinh(x))
    return epsilon_0 * quant if include_constant else quant

def pressure(x, include_constant=True):
    quant = (1/24) * ((2 * x**3 - 3 * x) * np.sqrt(1 + x**2) + 3 * np.arcsinh(x))
    return epsilon_0 * quant if include_constant else quant

Equation-of-state $p(\varepsilon_0)$

There is no simple way to express $p$ as a function of $\varepsilon_\mathrm{elec}$. But we can plot the relation numerically, or look at both relativistic and classical limits.

Relativistic limit: $x \gg 1$ $$\begin{aligned} p(x) &= \varepsilon_0\cdot\frac{1}{3}\int_0^x (1 + u^2)^{-1/2}u^4 \ du \\ &\rightarrow \varepsilon_0\cdot\frac{1}{3}\int_0^x u^{-1}u^4 \ du \\ &= \varepsilon_0\cdot\frac{x^4}{12} \end{aligned} $$

Then using $$ k_F = \hbar\left(\frac{3\pi^2\rho}{m_N}\frac{Z}{A}\right)^{1/3} \tag{9} $$ and $\varepsilon_\mathrm{elec} = \rho c^2$ (from $E = mc^2$), we get $$ p = \frac{\hbar c}{12\pi^2}\left(\frac{3\pi^2}{m_Nc^2}\frac{Z}{A}\right)^{4/3}\cdot\varepsilon_\mathrm{elec}^{4/3} = K_\mathrm{rel}\varepsilon_\mathrm{elec}^{4/3} $$

def relativistic_pressure(x, include_constant=True):
    quant = x**4 / 12
    return epsilon_0 * quant if include_constant else quant

def K_rel(A=1, Z=0.5):
    K = (const.hbar * const.c) / (12 * np.pi**2)
    K *= ((3 * np.pi**2 * Z) / (const.u * const.c**2 * A))**(4/3)
    return K

def relativistic_EOS(epsilon, A=1, Z=0.5):
    return K_rel(A, Z) * epsilon**(4/3)

Non-relativistic limit: $x \ll 1$ $$ p(x) \rightarrow \varepsilon_0\cdot\frac{x^5}{15}$$ And the E.O.S. is $$ p = \frac{\hbar^2}{15\pi^2m_e}\left(\frac{3\pi^2}{m_Nc^2}\frac{Z}{A}\right)^{5/3}\cdot\varepsilon_\mathrm{elec}^{5/3} = K_\mathrm{nonrel}\varepsilon_\mathrm{elec}^{5/3} $$

def non_relativistic_pressure(x, include_constant=True):
    quant = x**5 / 15
    return epsilon_0 * quant if include_constant else quant

def K_nonrel(A=1, Z=0.5):
    K = (const.hbar**2) / (15 * np.pi**2 * const.m_e)
    K *= ((3 * np.pi**2 * Z) / (const.u * const.c**2 * A))**(5/3)
    return K

def non_relativistic_EOS(epsilon, A=1, Z=0.5):
    return K_nonrel(A, Z) * epsilon**(5/3)

fig, ax = plt.subplots(dpi=160, figsize=(3.5, 3.0))
x = np.linspace(0, 10, 100) # fermi momentum
unit = 'erg/cm^3'
ax.plot(x, pressure(x).to(unit), color='black', label='exact')
ax.plot(x, relativistic_pressure(x).to(unit), c='blue', lw=0.6, ls='dashed', label=r'relativistic ($\gamma=\frac{4}{3}$)')
ax.plot(x, non_relativistic_pressure(x).to(unit), c='red', lw=0.6, ls='dashed', label=r'classical ($\gamma=\frac{5}{3}$)')
ax.legend()
ax.set_xlabel(r'Fermi momentum $x \equiv k_F/m_ec$')
ax.set_ylabel(r'Pressure $p(x)$ (erg/cm$^3$)')
ax.set_xscale('log')
ax.set_yscale('log')
plt.show()

findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.

And just to show give us a sense of the unit erg/cm$^3$, I present the conversion at $10^{27}$ erg/cm$^3$ into eV/fm$^3$:

(1e27 * u.erg / u.cm**3).to('eV/fm^3')

fig, ax = plt.subplots(ncols=2, dpi=150, figsize=(8, 3.5))
x = np.linspace(0, 2, 100)
c = 0
ax[c].plot(energy_density(x).to('erg/cm^3'), pressure(x).to('erg/cm^3'), c='red')
ax[c].set_title('Attempt to reproduce Fig. 2')
ax[c].set_xlabel(r'$\varepsilon_\mathrm{elec}$ (erg/cm$^3$)')
ax[c].set_ylabel(r'$p$ (erg/cm$^3$)')

c = 1
ax[c].plot(energy_density(x, include_constant=False), pressure(x, include_constant=False), c='red')
ax[c].set_title(r'Same, but as multiples of $\varepsilon_0$')
ax[c].set_xlabel(r'$\varepsilon_\mathrm{elec}/\varepsilon_0$')
ax[c].set_ylabel(r'$p/\varepsilon_0$')
plt.show()

Mathematica can give us the integrals: $$\varepsilon_\mathrm{elec}(x) = \varepsilon_0\cdot\frac{1}{8}\left[(2x^3 + x)\sqrt{1+x^2} - \sinh^{-1}(x)\right] \tag{10}$$ $$p(x) = \varepsilon_0\cdot\frac{1}{24}\left[(2x^3 - 3x)\sqrt{1+x^2} + 3\sinh^{-1}(x)\right] \tag{13}$$

The structure equations for a polytrope¶

Polytrope: $$ p(\varepsilon) = K\varepsilon^\gamma $$

To facilitate numerical calculations, we rewrite the structure equations for a polytrope: $$ \begin{cases} \dfrac{d\overline{p}(r)}{dr} &= -\dfrac{\alpha \overline{p}(r)^{1/\gamma}\overline{M}(r)}{r^2} \\ \dfrac{d\overline{M}(r)}{dr} &= \beta r^2\overline{p}(r)^{1/\gamma} \end{cases} $$ The only dimension ($r$) remains in this ODE is length, and we choose to use kilometer for convenient.

@numba.jit
def ode(r, y, alpha=1.0, beta=1.0, gamma=1.0):
    p, M = y[0], y[1]
    deriv_p = -alpha * p**(1/gamma) * M / r**2
    deriv_M = beta * r**2 * p*(1/gamma)
    return [deriv_p, deriv_M]

Our goal is to solve for the dimensionless $\overline{p}(r) = p(r) / \epsilon_0$ and $\overline{M}(r) = M(r) / M_\odot$.

Notice that $\epsilon_0$ is **not** the $\varepsilon_0$ in earlier section, but a *free parameter* whose value is only estimated to give reasonable solutions.

Both $\alpha$ and $\beta$ are derived from $\epsilon_0$: $$\alpha(\epsilon_0; K, \gamma) = \frac{R_0}{(K\epsilon_0^{\gamma-1})^{1/\gamma}} \tag{24}$$ $$\beta(\epsilon_0; K, \gamma) = \frac{4\pi\epsilon_0}{M_\odot c^2(K\epsilon_0^{\gamma-1})^{1/\gamma}} \tag{24}$$ Here, $R_0 = GM_\odot/c^2 \approx 1.477$ km.

R_0 = const.G * const.M_sun / const.c**2

def get_alpha(epsilon_est, K, gamma):
    alpha = R_0 / (K * epsilon_est**(gamma - 1))**(1/gamma)
    return alpha.to('km')

def get_beta(epsilon_est, K, gamma):
    beta = (4 * np.pi * epsilon_est) / (const.M_sun * const.c**2)
    beta /= (K * epsilon_est**(gamma - 1))**(1/gamma)
    return beta.to('1/km^3')

Integrating the polytrope numerically¶

We will use two initial conditions:

$\overline{M}(0) = 0$ for obvious reason;
$\overline{p}(0)$ is the pressure at the core. The choice for this value is not so trivial.

Useful reminder: $\overline{p}(r) = p(r)/\epsilon_0$ or $p(r) = \epsilon_0\overline{p}(r)$

Usually we pick $\epsilon_0$ such that $\alpha(\epsilon_0)$ and $\beta(\epsilon_0)$ have comparable orders of magnitude (for numerical stability when solving ODE).

We don't really know the core pressure $p(0)$. But we can try different $p(0)$, hence $\overline{p}(0)$, and see if the $R$ and $M$ it predicts match our observation.

The relativistic case $k_F \gg m_e$¶

gamma = 4/3

epsilon_est = (const.M_sun * const.c**2) / (1e4 * u.km)**3
epsilon_est.to('erg/cm^3')

epsilon_est = 1.0e39 * u.erg / u.cm**3
epsilon_est.to('erg/cm^3')

alpha = get_alpha(epsilon_est, K_rel(), gamma)
beta = get_beta(epsilon_est, K_rel(), gamma)
display(alpha, beta)

sol = solve_ivp(ode, t_span=[1e-5, 0.2e5], y0=[1e-16, 0.0],
                #args=[1.473, 52.46, gamma],
                args=[alpha, beta, gamma],
                rtol=1e-8, atol=1e-14,
                dense_output=True,
               )

fig, ax = plt.subplots(ncols=2, dpi=160, figsize=(6.0, 3.0))
# r = np.linspace(0, 1e5, 100)
r = np.linspace(0, 2e4, 100)
c = 0
ax[c].plot(r, sol.sol(r)[0])
ax[c].set_title('pressure')

c = 1
ax[c].plot(r, sol.sol(r)[1])
ax[c].set_title('mass')

for axe in ax:
    axe.set_xlabel('Radius (km)')
plt.show()

In Table I, the $\epsilon_0$ used is $$ \epsilon_0 = \frac{M_\odot c^2}{(10~\mathrm{km})^3} $$

((const.M_sun * const.c**2) / (10 * u.km)**3).to('erg/cm^3')

The non-relativistic case $k_F \ll m_e$¶

As core pressure becomes smaller, the electron gas is no longer relativistic. Since the electron gas can support less mass, this results in larger dwarfs.

gamma = 5/3

$$ \epsilon_0 = 2.488\times10^{37}~\mathrm{erg/cm^3} = 0.01392M_\odot c^2/\mathrm{km}^3 \tag{36}$$

Power law ansatz

Which I called "asymptotic matching": $$ \overline{\epsilon}(\overline{p}) = A_\mathrm{nonrel}\overline{p}^{3/5} + A_\mathrm{rel}\overline{p}^{3/4} $$

Pure neutron star, fermi gas equation of state¶

General relativistic corrections must be included, i.e. TOV equation $$ \frac{dp}{dr} = -\frac{G\varepsilon(r)M(r)}{r^2c^2}\left[1 + \frac{p(r)}{\varepsilon(r)}\right]\left[1 + \frac{4\pi r^3p(r)}{M(r)c^2}\right]\left[1 - \frac{2GM(r)}{rc^2}\right]^{-1} $$
Real neutron stars consists not just neutrons, but also a small portion of protons and electrons
Fermi gas model ignores the strong nucleon-nucleon interactions, which contributes to the energy density

The non-relativistic case $k_F \ll m_n$¶

gamma = 5/3

def K_nonrel(A=1, Z=0.5):
    K = (const.hbar**2) / (15 * np.pi**2 * const.m_n)
    K *= ((3 * np.pi**2 * Z) / (const.m_n * const.c**2 * A))**(5/3)
    return K

K_nonrel(A=1, Z=1).to('cm^2 / erg^(2/3)')

Typical neutron stars: $$ R = 10~\mathrm{km} \hspace{30pt} M = 1.0 M_\odot $$

The relativistic case $k_F \gg m_n$¶

No, it's not $\gamma = 4/3$. Instead, we have $\gamma = 1$.

Would find $R = 50$ km!

Basically, this is incorrect --- try the power law ansatz.

The Fermi gas equation of state for arbitrary relativity¶

Power law ansatz

$$ \overline{\epsilon}(\overline{p}) = A_\mathrm{nonrel}\overline{p}^{3/5} + A_\mathrm{rel}\overline{p} $$

Parameters can be found be standard fitting routine: $$ A_\mathrm{nonrel} = 2.4216 \hspace{30pt} A_\mathrm{rel} = 2.8663 \tag{46} $$

With a choice of $$\epsilon_0 =\frac{m_n^4c^5}{(3\pi^2\hbar)^3} = 5.346\times10^{36}~\mathrm{erg/cm^3} \tag{45} $$ the authors found for a starting value of $\overline{p}(0) = 10^{-2}$: $$ R = 15.0~\mathrm{km} \hspace{30pt} M = 1.037M_\odot \hspace{30pt}\textsf{(Newtonian equations)} \tag{47} $$ $$ R = 13.4~\mathrm{km} \hspace{30pt} M = 0.717M\odot \hspace{30pt}\textsf{(full TOV equation)} \tag{48} $$

Why is there a maximum mass?¶

Maximum mass $\rightarrow$ exists a maximum pressure to support the mass (otherwise collapse)
Thermal contribution toward pressure is negligible in cold stars
So pressure increment $dp$ mainly comes from increase in energy density $d\varepsilon$: $$ \left(\frac{c_s}{c}\right)^2 = \frac{dp}{d\varepsilon} \tag{92} $$
Special relativity: $c_s < c$
Conclusion: increase in $\varepsilon$ eventually will not return sufficient increase in $p$ to support the gravity!