Part II: Compact Objects | Chapter 2

Neutron Stars

The densest observable objects in the Universe: the TOV equation, nuclear equations of state, pulsar timing, magnetars, and glitch phenomena

Overview

Neutron stars are born in the gravitational collapse of massive stellar cores during core-collapse supernovae. With masses of $1.2\text{--}2.2\,M_\odot$ compressed into radii of just $\sim 10\text{--}13$ km, they achieve nuclear densities $\rho \sim 10^{17}\text{--}10^{18}$ kg m$^{-3}$, making them natural laboratories for the physics of ultra-dense matter. The prediction of neutron stars by Baade and Zwicky (1934) and the theoretical framework of Oppenheimer and Volkoff (1939) preceded their observational discovery as pulsars by Jocelyn Bell Burnell and Antony Hewish in 1967.

In this chapter we derive the Tolman-Oppenheimer-Volkoff equation from general relativity, explore the nuclear equation of state and its implications for the maximum neutron star mass, analyze pulsar spin-down and the P-Pdot diagram, and discuss magnetars and the glitch phenomenon.

1. The Tolman-Oppenheimer-Volkoff Equation

The internal structure of a neutron star requires general relativity because the gravitational field is strong: the compactness parameter $\mathcal{C} = GM/(Rc^2) \sim 0.2$ is far from the Newtonian limit. The TOV equation generalizes Newtonian hydrostatic equilibrium to the Schwarzschild metric.

1.1 Starting from the Schwarzschild Interior Solution

For a static, spherically symmetric perfect fluid, the metric takes the form:

$$ds^2 = -e^{2\Phi(r)}c^2\,dt^2 + \left(1 - \frac{2Gm(r)}{rc^2}\right)^{-1}dr^2 + r^2\,d\Omega^2$$

The Einstein field equations $G_{\mu\nu} = 8\pi G T_{\mu\nu}/c^4$ for a perfect fluid with energy density $\epsilon$ and pressure $P$ yield the mass equation:

$$\frac{dm}{dr} = \frac{4\pi r^2 \epsilon}{c^2}$$

and the pressure equation, which is the TOV equation:

$$\boxed{\frac{dP}{dr} = -\frac{G\epsilon\, m}{r^2 c^2}\left(1 + \frac{P}{\epsilon}\right)\left(1 + \frac{4\pi r^3 P}{mc^2}\right)\left(1 - \frac{2Gm}{rc^2}\right)^{-1}}$$

1.2 Physical Interpretation of the Correction Factors

Compared to the Newtonian hydrostatic equilibrium equation $dP/dr = -G\rho m/r^2$, the TOV equation contains three general-relativistic correction factors:

(1 + P/ε): In GR, pressure contributes to the gravitational mass-energy. This factor accounts for the "weight of pressure" — pressure itself gravitates.

(1 + 4πr$^3$P/mc$^2$): The pressure of the overlying layers contributes to the enclosed gravitating mass.

(1 - 2Gm/rc$^2$)$^{-1}$: The metric correction factor from the curvature of spacetime; it increases the effective gravitational pull compared to flat space.

All three factors are greater than unity, meaning that general relativity makes gravity stronger than in the Newtonian case. This is why the maximum mass of a neutron star is lower than a naive Newtonian estimate would suggest.

2. The Nuclear Equation of State

The relationship between pressure and density at supra-nuclear densities is one of the major unsolved problems in physics. The EOS determines the mass-radius relation, the maximum mass, and the tidal deformability of neutron stars.

2.1 Nuclear Saturation Properties

Nuclear matter saturates at density $n_0 \approx 0.16$ fm$^{-3}$(corresponding to $\rho_0 \approx 2.7 \times 10^{17}$ kg m$^{-3}$) with binding energy per nucleon $E/A \approx -16$ MeV. The energy per nucleon near saturation can be expanded as:

$$\frac{E}{A}(n, \delta) = \frac{E}{A}(n_0, 0) + \frac{K_0}{18}\left(\frac{n - n_0}{n_0}\right)^2 + S(n)\,\delta^2 + \cdots$$

where $K_0 \approx 240$ MeV is the nuclear incompressibility, $\delta = (n_n - n_p)/(n_n + n_p)$is the isospin asymmetry, and $S(n)$ is the density-dependent symmetry energy. For pure neutron matter, $\delta = 1$.

2.2 Beta Equilibrium

In neutron star matter, the composition is determined by beta equilibrium, where the reactions $n \to p + e^- + \bar{\nu}_e$ and $p + e^- \to n + \nu_e$are in balance. The equilibrium condition is:

$$\mu_n = \mu_p + \mu_e$$

Combined with charge neutrality $n_p = n_e$, this yields a proton fraction of $x_p \approx 4\text{--}10\%$ at nuclear density, rising slowly with increasing density. The low proton fraction means the pressure is dominated by neutron degeneracy and nuclear interactions.

2.3 Exotic Matter at High Densities

At densities exceeding $2\text{--}3\,n_0$, various exotic phases may appear: hyperons ($\Lambda$, $\Sigma$, $\Xi$), a Bose-Einstein condensate of kaons, or a transition to deconfined quark matter (up, down, strange quarks). Each possibility produces a distinct EOS and hence different mass-radius predictions. The detection of $\sim 2\,M_\odot$ neutron stars (PSR J1614-2230, PSR J0348+0432, PSR J0740+6620) has placed strong constraints on the EOS, ruling out many models that predict too-soft equations of state at high density.

3. Pulsar Spin-Down and the P-Pdot Diagram

Pulsars are rapidly rotating, highly magnetized neutron stars that emit beamed electromagnetic radiation. The observed pulsed emission arises from the lighthouse effect as the beam sweeps across the observer's line of sight.

3.1 Magnetic Dipole Radiation

A rotating magnetic dipole radiates energy at a rate given by the Larmor formula for a time-varying dipole moment. For a neutron star with surface magnetic field $B$, radius $R$, and angular velocity $\Omega = 2\pi/P$:

$$\dot{E}_{\text{dipole}} = -\frac{B^2 R^6 \Omega^4}{6c^3}\sin^2\alpha$$

where $\alpha$ is the angle between the rotation and magnetic axes. Setting this equal to the rotational kinetic energy loss rate $\dot{E}_{\text{rot}} = I\Omega\dot{\Omega}$:

$$I\Omega\dot{\Omega} = -\frac{B^2 R^6 \Omega^4}{6c^3}\sin^2\alpha$$

This gives the spin-down law:

$$\boxed{\dot{\Omega} = -\frac{B^2 R^6 \sin^2\alpha}{6Ic^3}\,\Omega^3}$$

3.2 Characteristic Age and Magnetic Field

From $\dot{\Omega} \propto \Omega^n$ with braking index $n = 3$ (magnetic dipole), the characteristic age is:

$$\boxed{\tau_c = \frac{P}{2\dot{P}}}$$

and the surface magnetic field strength is estimated as:

$$\boxed{B_s = 3.2 \times 10^{19} \sqrt{P\dot{P}} \;\text{G}}$$

assuming $R = 10$ km, $I = 10^{45}$ g cm$^2$, and $\sin\alpha = 1$. These relations define contours of constant $\tau_c$and constant $B_s$ on the P-Pdot diagram, the fundamental diagnostic tool for pulsar populations.

4. Magnetars

Magnetars are neutron stars with ultra-strong magnetic fields, $B \sim 10^{14}\text{--}10^{15}$ G (compared to $\sim 10^{12}$ G for typical pulsars). They manifest as soft gamma repeaters (SGRs) and anomalous X-ray pulsars (AXPs).

4.1 Magnetic Energy Budget

The total magnetic energy stored in a magnetar is:

$$E_B = \frac{B^2}{8\pi}\cdot\frac{4}{3}\pi R^3 \approx 1.6 \times 10^{46}\left(\frac{B}{10^{15}\,\text{G}}\right)^2\left(\frac{R}{10\,\text{km}}\right)^3 \;\text{erg}$$

This exceeds the rotational energy by orders of magnitude for slowly-rotating magnetars (with $P \sim 2\text{--}12$ s). The persistent X-ray luminosity of magnetars ($L_X \sim 10^{34}\text{--}10^{36}$ erg/s) is therefore powered by magnetic field decay rather than spin-down, as proposed by Duncan and Thompson (1992).

4.2 Giant Flares

Magnetars occasionally produce giant flares with peak luminosities exceeding $10^{44}$ erg/s. The 27 December 2004 giant flare from SGR 1806-20 released $\sim 2 \times 10^{46}$ erg in about 0.2 seconds. These events are attributed to catastrophic rearrangements of the magnetic field (starquakes), analogous to coronal mass ejections on the Sun but powered by fields $10^{10}$ times stronger. The oscillation spectrum of the decaying tail reveals neutron star seismology modes that constrain the EOS.

5. Pulsar Glitches and Superfluidity

Pulsar glitches are sudden increases in the rotation rate, $\Delta\Omega/\Omega \sim 10^{-9}\text{--}10^{-6}$, observed in many pulsars. The Vela pulsar glitches approximately every 3 years.

5.1 The Two-Component Model

The standard explanation involves a superfluid component in the neutron star interior. At the densities found in neutron star interiors, neutrons form Cooper pairs and become superfluid. The superfluid rotates by forming quantized vortices, each carrying angular momentum $\hbar$ per neutron. The areal density of vortices is:

$$n_v = \frac{2m_n \Omega}{\pi\hbar} \approx 6 \times 10^3 \left(\frac{P}{33\,\text{ms}}\right)^{-1}\;\text{cm}^{-2}$$

As the star spins down, the superfluid component (which is decoupled from the crust on short timescales) maintains its rotation while the crust slows. The angular velocity lag builds until the vortices catastrophically unpin and transfer angular momentum to the crust, producing the observed glitch:

$$\boxed{\frac{\Delta\Omega}{\Omega} = \frac{I_{\text{sf}}}{I_{\text{total}}}\frac{\Delta\Omega_{\text{lag}}}{\Omega}}$$

The cumulative glitch activity of the Vela pulsar requires that at least 1.6% of the star's moment of inertia resides in the superfluid component. This constraint, combined with EOS models, implies that the superfluid extends through much of the inner crust and possibly into the outer core.

5.2 Glitch Recovery

After a glitch, pulsars typically exhibit exponential recovery toward the pre-glitch spin-down rate, with timescales ranging from days to months. The recovery timescale probes the coupling between the superfluid interior and the normal crust, providing unique constraints on the viscosity and mutual friction in neutron star interiors.

Applications

Gravitational Wave Constraints on the EOS

The gravitational wave signal from the binary neutron star merger GW170817 (detected by LIGO/Virgo on 17 August 2017) provided the first measurement of the tidal deformability parameter $\Lambda$, which directly constrains the EOS. Softer equations of state produce larger, more deformable neutron stars. The measured value $\tilde{\Lambda} = 300^{+420}_{-230}$ rules out very stiff EOS models and implies neutron star radii of $R \approx 11\text{--}13$ km.

Millisecond Pulsars as Precision Clocks

Millisecond pulsars (MSPs), spun up by accretion from a companion, have rotational stabilities rivaling atomic clocks. Pulsar timing arrays (PTAs) monitor ensembles of MSPs to search for nanohertz gravitational waves from supermassive black hole binaries. The NANOGrav, EPTA, PPTA, and IPTA collaborations reported evidence for a stochastic gravitational wave background in 2023, opening a new window on the Universe.

Historical Notes

Walter Baade and Fritz Zwicky predicted neutron stars in 1934, just two years after the discovery of the neutron by James Chadwick. They proposed that neutron stars are formed in supernovae and represent the final state of stellar collapse. Oppenheimer and Volkoff (1939) computed the first general-relativistic models. However, most astronomers considered neutron stars unobservable until Jocelyn Bell Burnell discovered the first pulsar (PSR B1919+21) in 1967 during her PhD work. Tommy Gold immediately identified pulsars as rotating neutron stars. The discovery of the Crab pulsar (with $P = 33$ ms) in the remnant of the 1054 CE supernova confirmed the connection. Russell Hulse and Joseph Taylor's discovery of the binary pulsar PSR B1913+16 in 1974 provided the first indirect evidence for gravitational wave emission, earning them the 1993 Nobel Prize.

Neutron Star Internal Structure

A neutron star has a layered structure with dramatically different physics at each depth:

Outer crust ($\rho < 4 \times 10^{14}$ kg m$^{-3}$): A lattice of neutron-rich nuclei immersed in a degenerate electron gas. The composition progresses from $^{56}\text{Fe}$ at the surface through progressively more neutron-rich isotopes with depth, as electron capture ($e^- + p \to n + \nu_e$) drives neutronization.

Inner crust ($4 \times 10^{14}\text{--}2 \times 10^{17}$ kg m$^{-3}$): Above the neutron drip density, free neutrons appear alongside the nuclear lattice. These neutrons form a superfluid, responsible for the glitch phenomenon. Near the crust-core boundary, exotic "nuclear pasta" phases may form (sheets, rods, bubbles) as nuclear matter transitions from isolated nuclei to uniform matter.

Outer core ($2 \times 10^{17}\text{--}5 \times 10^{17}$ kg m$^{-3}$): Uniform nuclear matter consisting of neutrons, protons, electrons, and muons in beta equilibrium. Both neutrons and protons are superfluid/superconducting.

Inner core ($\rho > 5 \times 10^{17}$ kg m$^{-3}$): The composition is unknown. Possibilities include hyperons, kaon or pion condensates, or deconfined quark matter. The inner core may contain a color-superconducting phase of quarks. This is the regime where the EOS is least constrained and where observations of neutron star masses and radii have the greatest discriminating power.

Computational Exploration

The following simulation solves the TOV equation for several nuclear equations of state, produces the mass-radius diagram, constructs the P-Pdot diagram for the known pulsar population, and models pulsar spin-down evolution.

Neutron Star Structure, P-Pdot Diagram, and Spin-Down Evolution

Python

script.py332 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(14, 11))
fig.patch.set_facecolor('#0a0a0a')
for ax in axes.flat:
    ax.set_facecolor('#0f0f0f')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('#333')

# Constants
G = 6.674e-11
c_light = 3.0e8
M_sun = 1.989e30
hbar = 1.055e-34
m_n = 1.675e-27
MeV = 1.602e-13  # Joules
fm = 1e-15

# --- Plot 1: TOV solutions for different EOS ---
ax1 = axes[0, 0]

def tov_solve(eos_func, rho_c, dr=50.0, r_max=30000.0):
    """Solve TOV equation given EOS: P(rho), rho(P)"""
    r = dr
    eps = rho_c * c_light**2  # energy density
    P = eos_func(rho_c)
    m = (4.0/3.0) * np.pi * dr**3 * eps / c_light**2

r_arr = [r]
    m_arr = [m]
    P_arr = [P]

while r < r_max and P > 0:
        # TOV equation
        dPdr = -G * eps * m / (r**2 * c_light**2)
        dPdr *= (1 + P / eps)
        dPdr *= (1 + 4 * np.pi * r**3 * P / (m * c_light**2))
        factor = 1 - 2 * G * m / (r * c_light**2)
        if factor <= 0:
            break
        dPdr /= factor

P += dPdr * dr
        if P <= 0:
            break

# Get density from EOS (invert)
        rho = P / (eos_func(1.0) / 1.0)  # Simple polytropic inversion
        eps = rho * c_light**2

dmdr = 4 * np.pi * r**2 * eps / c_light**2
        m += dmdr * dr

r += dr
        r_arr.append(r)
        m_arr.append(m)
        P_arr.append(P)

return r / 1e3, m / M_sun  # R in km, M in solar masses

# Polytropic EOS models: P = K * rho^Gamma
def make_polytropic_eos(K, Gamma):
    def eos(rho):
        return K * rho**Gamma
    return eos

# Different stiffness models
eos_models = [
    ("Soft EOS", 5e-4, 2.0, '#44bbff'),
    ("Medium EOS", 2e-3, 2.2, '#ffa500'),
    ("Stiff EOS", 8e-3, 2.5, '#ff4444'),
]

for label, K, Gamma, color in eos_models:
    rho_c_values = np.logspace(17.5, 19.0, 80)
    M_list = []
    R_list = []

for rho_c in rho_c_values:
        eos = make_polytropic_eos(K, Gamma)
        # Simplified TOV integration
        r = 100.0
        m = (4.0/3.0) * np.pi * r**3 * rho_c
        P = K * rho_c**Gamma
        dr = 100.0
        rho = rho_c

while P > 0 and r < 50000:
            eps = rho * c_light**2
            if m <= 0:
                break
            dPdr = -G * eps * m / (r**2 * c_light**2)
            dPdr *= (1 + P / eps) if eps > 0 else 1
            dPdr *= (1 + 4*np.pi*r**3*P/(m*c_light**2)) if m > 0 else 1
            factor = 1 - 2*G*m/(r*c_light**2)
            if factor <= 0.01:
                break
            dPdr /= factor

P = P + dPdr * dr
            if P <= 0:
                break
            rho = (P / K)**(1.0/Gamma)
            m += 4*np.pi*r**2*rho*dr
            r += dr

M_list.append(m / M_sun)
        R_list.append(r / 1e3)

M_arr = np.array(M_list)
    R_arr = np.array(R_list)
    valid = (M_arr > 0.5) & (M_arr < 3.5) & (R_arr > 5) & (R_arr < 25)
    if np.any(valid):
        ax1.plot(R_arr[valid], M_arr[valid], color=color, linewidth=2.5, label=label)

# GW170817 constraint band
ax1.axhspan(1.17, 1.6, alpha=0.1, color='#ffa500', label='GW170817 range')
ax1.axhline(y=2.01, color='#ff4444', linestyle='--', linewidth=1.5, alpha=0.7, label='J0348+0432 (2.01 M_sun)')
ax1.axhline(y=2.08, color='#ff6644', linestyle=':', linewidth=1.5, alpha=0.7, label='J0740+6620 (2.08 M_sun)')

# NICER radius measurements
ax1.errorbar([13.0], [2.08], xerr=[[1.5],[1.5]], yerr=[[0.07],[0.07]], fmt='s', color='#ffdd00', markersize=8, capsize=4, label='NICER J0740')
ax1.errorbar([12.4], [1.44], xerr=[[1.0],[1.0]], yerr=[[0.15],[0.15]], fmt='o', color='#ffdd00', markersize=8, capsize=4, label='NICER J0030')

ax1.set_xlabel('Radius (km)', fontsize=12)
ax1.set_ylabel('Mass (M_sun)', fontsize=12)
ax1.set_title('Neutron Star Mass-Radius Relations', fontsize=14, fontweight='bold')
ax1.legend(fontsize=8, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white', loc='upper right')
ax1.set_xlim(6, 22)
ax1.set_ylim(0.5, 2.8)
ax1.grid(True, alpha=0.15)

# --- Plot 2: P-Pdot diagram ---
ax2 = axes[0, 1]

# Generate synthetic pulsar population
rng = np.random.default_rng(42)

# Normal pulsars
N_normal = 600
P_normal = 10**(rng.normal(-0.3, 0.5, N_normal))
Pdot_normal = 10**(rng.normal(-15, 1.0, N_normal))
# Ensure physical: Pdot > 0
Pdot_normal = np.abs(Pdot_normal)

# Millisecond pulsars
N_msp = 150
P_msp = 10**(rng.normal(-2.3, 0.3, N_msp))
Pdot_msp = 10**(rng.normal(-19.5, 0.7, N_msp))
Pdot_msp = np.abs(Pdot_msp)

# Magnetars
N_mag = 30
P_mag = 10**(rng.normal(0.7, 0.2, N_mag))
Pdot_mag = 10**(rng.normal(-11, 0.5, N_mag))
Pdot_mag = np.abs(Pdot_mag)

ax2.scatter(P_normal, Pdot_normal, s=3, color='#44bbff', alpha=0.6, label='Normal pulsars')
ax2.scatter(P_msp, Pdot_msp, s=8, color='#ffa500', alpha=0.8, label='Millisecond pulsars')
ax2.scatter(P_mag, Pdot_mag, s=25, color='#ff4444', marker='*', alpha=0.9, label='Magnetars')

# Lines of constant B field
P_line = np.logspace(-3, 2, 100)
for logB, label_B in [(8, '10^8 G'), (10, '10^{10} G'), (12, '10^{12} G'), (14, '10^{14} G'), (15, '10^{15} G')]:
    B = 10**logB
    Pdot_B = (B / 3.2e19)**2 / P_line
    ax2.plot(P_line, Pdot_B, ':', color='#666', linewidth=0.8, alpha=0.6)
    # Label
    idx = len(P_line) // 3
    if Pdot_B[idx] > 1e-22 and Pdot_B[idx] < 1e-8:
        ax2.text(P_line[idx], Pdot_B[idx]*2, f'B={label_B}', color='#888', fontsize=7, rotation=-40)

# Lines of constant characteristic age
for log_tau, label_tau in [(3, '1 kyr'), (6, '1 Myr'), (9, '1 Gyr')]:
    tau = 10**log_tau * 3.156e7  # seconds
    Pdot_tau = P_line / (2 * tau)
    ax2.plot(P_line, Pdot_tau, '--', color='#444', linewidth=0.8)
    idx = 2 * len(P_line) // 3
    if Pdot_tau[idx] > 1e-22:
        ax2.text(P_line[idx]*1.2, Pdot_tau[idx], label_tau, color='#666', fontsize=8)

# Death line
P_death = np.logspace(-1, 1.5, 50)
Pdot_death = 1e-16 * (P_death / 1.0)**(-1)
ax2.plot(P_death, Pdot_death, '-', color='#ff6644', linewidth=1.5, alpha=0.5, label='Death line')

ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.set_xlabel('Period P (s)', fontsize=12)
ax2.set_ylabel('Period derivative dP/dt (s/s)', fontsize=12)
ax2.set_title('Pulsar P-Pdot Diagram', fontsize=14, fontweight='bold')
ax2.legend(fontsize=8, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white', loc='upper left')
ax2.set_xlim(1e-3, 30)
ax2.set_ylim(1e-22, 1e-9)
ax2.grid(True, alpha=0.15)

# --- Plot 3: Spin-down evolution ---
ax3 = axes[1, 0]

# Evolve pulsar spin-down for different initial B fields
t_max = 1e8  # years
dt = 1e4  # years
t_arr = np.arange(0, t_max, dt)

B_values = [1e11, 1e12, 1e13, 1e14]
B_colors = ['#44bbff', '#ffa500', '#ffdd00', '#ff4444']
P0 = 0.02  # Initial period in seconds
I_ns = 1e45 * 1e-7  # kg m^2

for B_val, col in zip(B_values, B_colors):
    P_evol = np.zeros(len(t_arr))
    P_evol[0] = P0
    Omega = 2 * np.pi / P0

for i in range(1, len(t_arr)):
        # dOmega/dt = -B^2 R^6 Omega^3 / (6 I c^3)
        R_ns = 1e4  # 10 km in m... wait, use CGS consistently
        # Use formula: Pdot = (8pi^2 B^2 R^6) / (3 I c^3 P)  (CGS)
        # In SI: Pdot = B^2 * R^6 * 8*pi^2 / (6 * I * c^3 * P) -- approximate
        # Simpler: P(t) = P0 * sqrt(1 + t/tau_sd)
        # tau_sd = 3 I c^3 / (B^2 R^6 * 4*pi^2) * P0^2
        pass

# Analytic solution: P(t) = P0 * sqrt(1 + t/tau_sd)
    R_ns = 1e4  # 10 km
    # tau_sd in seconds (CGS-like)
    B_gauss = B_val * 1e-4  # convert to Tesla for SI
    # SI formula: tau_sd = 3*I*c^3*P0^2 / (8*pi^2 * B^2 * R^6)
    tau_sd = 3 * I_ns * c_light**3 * P0**2 / (8 * np.pi**2 * (B_val*1e-4)**2 * (1e4)**6)
    tau_sd_yr = tau_sd / (3.156e7)

P_evol = P0 * np.sqrt(1 + t_arr / tau_sd_yr)
    Pdot_evol = P0**2 / (2 * tau_sd) * (1 + t_arr / tau_sd_yr)**(-0.5) / 3.156e7

# Spin-down luminosity
    Edot = 4 * np.pi**2 * I_ns * Pdot_evol / (P_evol**3 * 3.156e7)

ax3.semilogy(t_arr / 1e6, P_evol, color=col, linewidth=2.5,
                 label=f'B = 10^{{{int(np.log10(B_val))}}} G')

ax3.set_xlabel('Age (Myr)', fontsize=12)
ax3.set_ylabel('Period (s)', fontsize=12)
ax3.set_title('Pulsar Spin-Down Evolution', fontsize=14, fontweight='bold')
ax3.legend(fontsize=9, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white')
ax3.set_xlim(0, 100)
ax3.set_ylim(0.01, 100)
ax3.grid(True, alpha=0.15)

# --- Plot 4: Glitch model ---
ax4 = axes[1, 1]

# Simulate glitch events in Vela-like pulsar
t_total = 50  # years
dt_fine = 0.001  # years
t = np.arange(0, t_total, dt_fine)
P_base = 0.0893  # Vela period in seconds
Pdot_base = 1.25e-13  # Vela Pdot

# Continuous spin-down
nu = 1.0 / P_base  # frequency
nudot = -Pdot_base / P_base**2

nu_arr = nu + nudot * t * 3.156e7  # Hz

# Add glitches every ~3 years
glitch_times = [3.1, 5.8, 9.2, 12.5, 15.1, 18.7, 22.0, 25.3, 28.9, 31.5, 35.2, 38.0, 41.5, 44.8, 48.1]
glitch_sizes = rng.uniform(1e-6, 3e-6, len(glitch_times))  # Delta_nu / nu

for gt, gs in zip(glitch_times, glitch_sizes):
    mask = t >= gt
    delta_nu = gs * nu
    # Instantaneous jump
    nu_arr[mask] += delta_nu
    # Exponential recovery (partial)
    Q = 0.7  # recovery fraction
    tau_rec = 0.1  # years
    recovery = Q * delta_nu * np.exp(-(t[mask] - gt) / tau_rec)
    nu_arr[mask] -= recovery

# Plot residuals (subtract linear trend)
nu_linear = nu + nudot * t * 3.156e7
residual = (nu_arr - nu_linear) * 1e6  # micro-Hz

ax4.plot(t, residual, color='#ffa500', linewidth=1.5)

# Mark glitches
for gt in glitch_times:
    ax4.axvline(x=gt, color='#ff4444', linewidth=0.5, alpha=0.3)

ax4.set_xlabel('Time (years)', fontsize=12)
ax4.set_ylabel('Frequency residual (uHz)', fontsize=12)
ax4.set_title('Vela-like Pulsar Glitch Activity', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.15)

# Annotate
ax4.text(25, max(residual)*0.8, 'Glitches shown as\nvertical markers', color='#ff6644', fontsize=10)

# Print summary
print("=" * 65)
print("NEUTRON STAR PHYSICS: KEY RESULTS")
print("=" * 65)

print("\nTOV Equation correction factors (for M=1.4 Msun, R=12 km):")
M_ns = 1.4 * M_sun
R_ns_m = 12e3
compactness = G * M_ns / (R_ns_m * c_light**2)
print(f"  Compactness: GM/(Rc^2) = {compactness:.3f}")
print(f"  Metric factor: (1-2GM/Rc^2)^-1 = {1/(1-2*compactness):.3f}")

print("\nCharacteristic ages and B-fields:")
pulsars = [("Crab", 0.0334, 4.21e-13), ("Vela", 0.0893, 1.25e-13),
           ("B1937+21", 0.00156, 1.05e-19), ("SGR 1806-20", 7.55, 8e-11)]
for name, P, Pdot in pulsars:
    tau_c = P / (2 * Pdot) / 3.156e7  # years
    B = 3.2e19 * np.sqrt(P * Pdot)
    print(f"  {name:12s}: P = {P:.4f} s, Pdot = {Pdot:.2e}, tau_c = {tau_c:.0f} yr, B = {B:.2e} G")

print("\nNuclear saturation density: n_0 = 0.16 fm^-3 = 2.7e17 kg/m^3")
print("Maximum observed NS mass: 2.08 +/- 0.07 M_sun (J0740+6620)")
print("NICER radius: 12.4 +/- 1.0 km (J0030+0451)")

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
print("\nPlot saved to output.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter Summary

Neutron stars are described by the TOV equation, the general-relativistic generalization of hydrostatic equilibrium, where pressure itself gravitates. The nuclear equation of state at supra-nuclear densities remains uncertain, but observations of $\sim 2\,M_\odot$neutron stars and NICER radius measurements constrain the EOS significantly.

Pulsars spin down by magnetic dipole radiation with $\dot{\Omega} \propto \Omega^3$, yielding characteristic age $\tau_c = P/(2\dot{P})$ and surface field$B = 3.2 \times 10^{19}\sqrt{P\dot{P}}$ G. Magnetars are powered by magnetic field decay rather than rotation.

Glitches reveal interior superfluidity: neutron Cooper pairs form quantized vortices whose catastrophic unpinning transfers angular momentum to the crust. The cumulative glitch activity constrains the fraction of the star's moment of inertia residing in the superfluid component.

Practice Problems

Problem 1:Estimate the maximum mass of a neutron star using a simple Newtonian argument. Balance the neutron degeneracy pressure $P \sim \hbar^2 n^{5/3}/(m_n)$ against gravitational pressure. Express the Chandrasekhar-like mass limit for neutrons and compare with the TOV limit of $\sim 2\text{-}2.5\,M_\odot$.

Solution:

1. By analogy with the Chandrasekhar mass: $M_{\text{Ch}} \sim \frac{(\hbar c)^{3/2}}{G^{3/2} m_n^2}$ (replacing electron mass with nucleon mass for the particle providing degeneracy pressure is not correct -- the scaling uses the baryon mass $m_B$ throughout).

2. The Chandrasekhar mass: $M_{\text{Ch}} = \frac{5.83}{\mu_e^2}\,M_\odot$ for electron degeneracy. For neutrons: $M_{\text{OV}} \sim 5.83\,M_\odot \times (m_e/m_n)^0 = 5.73/\mu_n^2\,M_\odot$.

3. More precisely, the non-relativistic to relativistic transition gives $M_{\max} \approx 0.7\,M_\odot$ for an ideal Fermi gas of neutrons (Oppenheimer-Volkoff 1939 original result).

4. Nuclear interactions (repulsive core at short range) stiffen the EOS, increasing $M_{\max}$ to $2\text{-}2.5\,M_\odot$.

5. General relativity lowers the maximum mass (pressure gravitates in GR), counteracting the nuclear stiffening. The TOV equation: $\frac{dP}{dr} = -\frac{(P + \rho c^2)(m + 4\pi r^3 P/c^2)G}{r^2 c^2(1 - 2Gm/rc^2)}$.

6. Observational constraint: PSR J0740+6620 has $M = 2.08 \pm 0.07\,M_\odot$, ruling out soft EOS models. The maximum mass lies between 2.1 and 2.5 $M_\odot$ for most viable equations of state.

Problem 2:A pulsar has period $P = 0.033\;\text{s}$ and period derivative $\dot{P} = 4.2 \times 10^{-13}$. Calculate the surface magnetic field strength using the standard dipole formula $B = 3.2 \times 10^{19}\sqrt{P\dot{P}}\;\text{G}$.

Solution:

1. The dipole spin-down formula: $B_s = 3.2 \times 10^{19}\sqrt{P\dot{P}}\;\text{G}$ (assuming $R = 10\;\text{km}$, $I = 10^{45}\;\text{g·cm}^2$).

2. $P\dot{P} = 0.033 \times 4.2 \times 10^{-13} = 1.386 \times 10^{-14}$.

3. $\sqrt{P\dot{P}} = \sqrt{1.386 \times 10^{-14}} = 1.177 \times 10^{-7}$.

4. $B_s = 3.2 \times 10^{19} \times 1.177 \times 10^{-7} = 3.77 \times 10^{12}\;\text{G}$.

5. This is $\sim 4 \times 10^{12}\;\text{G}$, typical of the Crab pulsar. This is consistent with the strong non-thermal emission and the Crab Nebula powered by the pulsar wind.

6. Derivation: equating rotational energy loss $\dot{E} = I\Omega\dot{\Omega} = -4\pi^2 I\dot{P}/P^3$ to magnetic dipole radiation power $P_{\text{rad}} = B^2R^6\Omega^4\sin^2\alpha/(6c^3)$, and solving for $B$ with $\sin\alpha = 1$.

Problem 3:Calculate the spin-down luminosity of the pulsar in Problem 2. Compare with the observed X-ray luminosity of the Crab Nebula ($L_X \approx 2 \times 10^{38}\;\text{erg/s}$). Use moment of inertia $I = 10^{45}\;\text{g·cm}^2$.

Solution:

1. Spin-down luminosity: $\dot{E} = -I\Omega\dot{\Omega} = \frac{4\pi^2 I \dot{P}}{P^3}$.

2. $\dot{E} = \frac{4\pi^2 \times 10^{45} \times 4.2 \times 10^{-13}}{(0.033)^3}$.

3. Numerator: $4\pi^2 \times 10^{45} \times 4.2 \times 10^{-13} = 1.659 \times 10^{34}$.

4. Denominator: $(0.033)^3 = 3.594 \times 10^{-5}$.

5. $\dot{E} = \frac{1.659 \times 10^{34}}{3.594 \times 10^{-5}} = 4.6 \times 10^{38}\;\text{erg/s} \approx 1.2 \times 10^5\,L_\odot$.

6. The Crab Nebula's total luminosity ($\sim 5 \times 10^{38}\;\text{erg/s}$ across all wavelengths) is comparable to $\dot{E}$, confirming the pulsar as the energy source. The X-ray fraction $L_X/\dot{E} \approx 0.4$ represents the efficiency of converting spin-down power to synchrotron X-rays via the pulsar wind nebula.

Problem 4:Calculate the Eddington luminosity for a $1.4\,M_\odot$ neutron star. If the observed luminosity during an X-ray burst is $L = 3.5 \times 10^{38}\;\text{erg/s}$, is the source super-Eddington? Assume pure hydrogen composition.

Solution:

1. Eddington luminosity: $L_{\text{Edd}} = \frac{4\pi G M c}{\kappa}$, where $\kappa = \sigma_T / m_p = 0.40\;\text{cm}^2/\text{g}$ for Thomson scattering off hydrogen.

2. $L_{\text{Edd}} = \frac{4\pi \times 6.674 \times 10^{-8} \times 1.4 \times 1.989 \times 10^{33} \times 3 \times 10^{10}}{0.40}$.

3. Numerator: $4\pi \times 6.674 \times 10^{-8} \times 2.785 \times 10^{33} \times 3 \times 10^{10} = 7.00 \times 10^{37}$.

4. $L_{\text{Edd}} = \frac{7.00 \times 10^{37}}{0.40} = 1.75 \times 10^{38}\;\text{erg/s}$.

5. In solar units: $L_{\text{Edd}} \approx 1.26 \times 10^{38} \times (M/M_\odot) = 1.76 \times 10^{38}\;\text{erg/s}$ for $1.4\,M_\odot$.

6. Since $L_{\text{obs}} = 3.5 \times 10^{38} = 2.0\,L_{\text{Edd}}$, the source is super-Eddington. This occurs during Type I X-ray bursts (thermonuclear flashes on the NS surface), where the photosphere expands due to radiation pressure exceeding gravity -- a photospheric radius expansion (PRE) burst, used as a standard candle for distance measurements.

Problem 5:A pulsar has $P = 0.5\;\text{s}$ and $\dot{P} = 5 \times 10^{-15}$. Calculate (a) the characteristic age $\tau_c$, (b) the surface B-field, and (c) determine whether this is a normal pulsar, millisecond pulsar, or magnetar on the P-$\dot{P}$ diagram.

Solution:

1. Characteristic (spin-down) age: $\tau_c = \frac{P}{2\dot{P}} = \frac{0.5}{2 \times 5 \times 10^{-15}} = 5 \times 10^{13}\;\text{s}$.

2. Converting: $\tau_c = \frac{5 \times 10^{13}}{3.156 \times 10^7} = 1.58 \times 10^6\;\text{yr} = 1.58\;\text{Myr}$.

3. Surface B-field: $B = 3.2 \times 10^{19}\sqrt{0.5 \times 5 \times 10^{-15}} = 3.2 \times 10^{19} \times 5 \times 10^{-8} = 1.6 \times 10^{12}\;\text{G}$.

4. Classification: $P = 0.5\;\text{s}$ and $B \sim 10^{12}\;\text{G}$ place this squarely in the normal pulsar population on the $P\text{-}\dot{P}$ diagram.

5. Comparison: millisecond pulsars have $P \sim 1\text{-}30\;\text{ms}$, $B \sim 10^{8\text{-}9}\;\text{G}$; magnetars have $P \sim 2\text{-}12\;\text{s}$, $B \sim 10^{14\text{-}15}\;\text{G}$.

6. Spin-down luminosity: $\dot{E} = 4\pi^2 I\dot{P}/P^3 = 4\pi^2 \times 10^{45} \times 5 \times 10^{-15} / 0.125 = 1.58 \times 10^{33}\;\text{erg/s}$. This is ~$0.4\,L_\odot$, modest compared to the Crab but still detectable in radio and potentially in X-rays if the pulsar has a wind nebula.

← White Dwarfs Supernovae →