pulsar_timing.tex

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{subcaption}
\usepackage[makestderr]{pythontex}

% Theorem environments
\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{example}{Example}
\newtheorem{remark}{Remark}

\title{Pulsar Timing: From Spin-down to Gravitational Wave Detection\\
\large A Comprehensive Analysis of Pulsar Astrophysics}
\author{High-Energy Astrophysics Division\\Computational Science Templates}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This comprehensive analysis explores pulsar timing theory and applications, from basic spin-down physics to gravitational wave detection via pulsar timing arrays. We derive the fundamental pulsar equations including period derivatives, characteristic ages, and magnetic field strengths. The analysis covers the P-$\dot{P}$ diagram for pulsar classification, timing residuals analysis, binary pulsar systems, and the search for nanohertz gravitational waves. We examine synthetic pulsar populations representing normal pulsars, millisecond pulsars, and magnetars, and explore their evolutionary relationships.
\end{abstract}

\section{Introduction}

Pulsars are rapidly rotating, highly magnetized neutron stars that emit beams of electromagnetic radiation. Discovered in 1967 by Jocelyn Bell Burnell and Antony Hewish, pulsars have become invaluable tools for precision astrophysics, from testing general relativity to detecting gravitational waves.

\begin{definition}[Pulsar]
A pulsar is a rotating neutron star with a dipolar magnetic field misaligned with its rotation axis. The lighthouse effect produces periodic pulses as the beam sweeps past Earth with each rotation.
\end{definition}

\section{Theoretical Framework}

\subsection{Spin-down Physics}

Pulsars lose rotational energy through magnetic dipole radiation and particle winds:

\begin{theorem}[Spin-down Luminosity]
The rate of rotational energy loss is:
\begin{equation}
\dot{E} = -\frac{d}{dt}\left(\frac{1}{2}I\Omega^2\right) = I\Omega\dot{\Omega} = \frac{4\pi^2 I \dot{P}}{P^3}
\end{equation}
where $I \approx 10^{45}$ g cm$^2$ is the neutron star moment of inertia, $P$ is the period, and $\dot{P}$ is the period derivative.
\end{theorem}

\subsection{Characteristic Age}

The characteristic age assumes constant braking index $n=3$ (pure magnetic dipole):

\begin{equation}
\tau_c = \frac{P}{(n-1)\dot{P}} = \frac{P}{2\dot{P}}
\end{equation}

\begin{remark}[Age Limitations]
The characteristic age equals the true age only if:
\begin{itemize}
    \item Birth period $P_0 \ll P$ (current period)
    \item Braking index $n = 3$ (magnetic dipole)
    \item Magnetic field is constant
\end{itemize}
For young pulsars, $\tau_c$ often overestimates the true age.
\end{remark}

\subsection{Magnetic Field Strength}

Equating spin-down luminosity to magnetic dipole radiation:

\begin{theorem}[Surface Magnetic Field]
The equatorial surface field strength is:
\begin{equation}
B = \sqrt{\frac{3c^3 I P\dot{P}}{8\pi^2 R^6 \sin^2\alpha}} \approx 3.2 \times 10^{19} \sqrt{P\dot{P}} \text{ G}
\end{equation}
where $R \approx 10$ km is the neutron star radius and $\alpha$ is the inclination angle.
\end{theorem}

\subsection{Braking Index}

The braking index describes the spin-down mechanism:
\begin{equation}
n = \frac{\Omega\ddot{\Omega}}{\dot{\Omega}^2} = 2 - \frac{P\ddot{P}}{\dot{P}^2}
\end{equation}

\begin{itemize}
    \item $n = 3$: Pure magnetic dipole radiation
    \item $n = 1$: Particle wind dominated
    \item $n = 5$: Gravitational wave emission
\end{itemize}

\section{Computational Analysis}

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

np.random.seed(42)

# Physical constants
c = 2.998e10  # cm/s
I = 1e45  # Moment of inertia (g cm^2)
R_ns = 1e6  # Neutron star radius (cm)

# Generate synthetic pulsar population
# Normal pulsars (canonical)
n_normal = 200
P_normal = 10**np.random.uniform(-0.3, 1.0, n_normal)  # 0.5 - 10 s
Pdot_normal = 10**np.random.uniform(-16, -13, n_normal)

# Millisecond pulsars (recycled)
n_msp = 80
P_msp = 10**np.random.uniform(-2.5, -1.3, n_msp)  # 3 - 50 ms
Pdot_msp = 10**np.random.uniform(-21, -18, n_msp)

# Magnetars (Anomalous X-ray Pulsars / Soft Gamma Repeaters)
n_mag = 30
P_mag = 10**np.random.uniform(0.0, 1.1, n_mag)  # 1 - 12 s
Pdot_mag = 10**np.random.uniform(-12, -10, n_mag)

# Combine all populations
populations = {
    'Normal': {'P': P_normal, 'Pdot': Pdot_normal, 'color': 'blue', 'marker': 'o', 'size': 8},
    'MSP': {'P': P_msp, 'Pdot': Pdot_msp, 'color': 'red', 'marker': 's', 'size': 8},
    'Magnetar': {'P': P_mag, 'Pdot': Pdot_mag, 'color': 'purple', 'marker': '^', 'size': 12}
}

# Calculate derived quantities for all pulsars
all_P = np.concatenate([P_normal, P_msp, P_mag])
all_Pdot = np.concatenate([Pdot_normal, Pdot_msp, Pdot_mag])
all_tau = all_P / (2 * all_Pdot) / 3.154e7 / 1e6  # Myr
all_B = 3.2e19 * np.sqrt(all_P * all_Pdot)  # Gauss
all_Edot = 4 * np.pi**2 * I * all_Pdot / all_P**3  # erg/s

# Famous pulsars for reference
famous_pulsars = {
    'Crab (B0531+21)': {'P': 0.0334, 'Pdot': 4.21e-13},
    'Vela (B0833-45)': {'P': 0.0893, 'Pdot': 1.25e-13},
    'PSR B1937+21': {'P': 0.00156, 'Pdot': 1.05e-19},
    'PSR B1913+16': {'P': 0.0590, 'Pdot': 8.63e-18},
    'SGR 1806-20': {'P': 7.55, 'Pdot': 8.3e-11}
}

# Calculate properties for famous pulsars
for name, props in famous_pulsars.items():
    props['tau'] = props['P'] / (2 * props['Pdot']) / 3.154e7
    props['B'] = 3.2e19 * np.sqrt(props['P'] * props['Pdot'])
    props['Edot'] = 4 * np.pi**2 * I * props['Pdot'] / props['P']**3

# Generate timing residuals
def generate_timing_residuals(n_epochs, rms_noise, has_gw=False):
    """Generate synthetic timing residuals"""
    t = np.linspace(0, 10, n_epochs)  # years
    residuals = np.random.normal(0, rms_noise, n_epochs)

    if has_gw:
        # Add GW-induced correlation (simplified red noise)
        gw_signal = 100e-9 * np.sin(2*np.pi*t/5)  # 5-year period
        residuals += gw_signal

    return t, residuals

# Pulsar timing array (PTA) concept
n_pta_pulsars = 20
pta_rms = np.random.uniform(50, 500, n_pta_pulsars)  # ns

# Create comprehensive figure
fig = plt.figure(figsize=(14, 16))

# Plot 1: P-Pdot diagram
ax1 = fig.add_subplot(3, 3, 1)
for name, pop in populations.items():
    ax1.scatter(pop['P'], pop['Pdot'], c=pop['color'], s=pop['size'],
               alpha=0.5, marker=pop['marker'], label=name)

# Add famous pulsars
for name, props in famous_pulsars.items():
    ax1.plot(props['P'], props['Pdot'], 'k*', markersize=10)

# Add constant age lines
P_line = np.logspace(-3, 1.5, 100)
for age_yr, style in [(1e3, '--'), (1e6, '-'), (1e9, ':')]:
    Pdot_line = P_line / (2 * age_yr * 3.154e7)
    label_age = f'{age_yr/1e6:.0f} Myr' if age_yr >= 1e6 else f'{age_yr/1e3:.0f} kyr'
    ax1.loglog(P_line, Pdot_line, 'g', alpha=0.4, linewidth=1, linestyle=style)

# Add constant B lines
for B_val in [1e10, 1e12, 1e14]:
    Pdot_B = (B_val / 3.2e19)**2 / P_line
    ax1.loglog(P_line, Pdot_B, 'orange', alpha=0.4, linewidth=1, linestyle=':')

# Death line
P_death = np.logspace(-3, 1, 100)
Pdot_death = (P_death / 0.17)**2 * 1e-16
ax1.loglog(P_death, Pdot_death, 'k--', alpha=0.7, linewidth=2, label='Death Line')

ax1.set_xlabel('Period (s)')
ax1.set_ylabel('Period Derivative (s/s)')
ax1.set_title('Pulsar P-$\\dot{P}$ Diagram')
ax1.legend(fontsize=7, loc='lower right')
ax1.set_xlim([1e-3, 20])
ax1.set_ylim([1e-22, 1e-9])
ax1.grid(True, alpha=0.3)

# Plot 2: Period histogram by population
ax2 = fig.add_subplot(3, 3, 2)
bins = np.linspace(-3, 1.5, 30)
for name, pop in populations.items():
    ax2.hist(np.log10(pop['P']), bins=bins, alpha=0.6, label=name, color=pop['color'])
ax2.set_xlabel('$\\log_{10}(P$/s)')
ax2.set_ylabel('Count')
ax2.set_title('Period Distribution')
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)

# Plot 3: Magnetic field distribution
ax3 = fig.add_subplot(3, 3, 3)
ax3.hist(np.log10(all_B), bins=30, alpha=0.7, color='green', edgecolor='black')
ax3.axvline(x=8, color='r', linestyle='--', alpha=0.7, label='MSP')
ax3.axvline(x=12, color='b', linestyle='--', alpha=0.7, label='Normal')
ax3.axvline(x=15, color='purple', linestyle='--', alpha=0.7, label='Magnetar')
ax3.set_xlabel('$\\log_{10}(B$/G)')
ax3.set_ylabel('Count')
ax3.set_title('Magnetic Field Distribution')
ax3.legend(fontsize=7)
ax3.grid(True, alpha=0.3)

# Plot 4: Spin-down luminosity vs age
ax4 = fig.add_subplot(3, 3, 4)
for name, pop in populations.items():
    tau_i = pop['P'] / (2 * pop['Pdot']) / 3.154e7 / 1e6  # Myr
    Edot_i = 4 * np.pi**2 * I * pop['Pdot'] / pop['P']**3
    ax4.scatter(tau_i, Edot_i, c=pop['color'], s=pop['size'],
               alpha=0.5, marker=pop['marker'], label=name)
ax4.set_xscale('log')
ax4.set_yscale('log')
ax4.set_xlabel('Characteristic Age (Myr)')
ax4.set_ylabel('$\\dot{E}$ (erg/s)')
ax4.set_title('Spin-down Luminosity Evolution')
ax4.legend(fontsize=7)
ax4.grid(True, alpha=0.3)

# Plot 5: Timing residuals example
ax5 = fig.add_subplot(3, 3, 5)
t_res, residuals = generate_timing_residuals(100, 100e-9, has_gw=False)
ax5.errorbar(t_res, residuals*1e6, yerr=100e-3, fmt='b.', capsize=2, alpha=0.7)
ax5.axhline(y=0, color='r', linestyle='-', alpha=0.5)
ax5.set_xlabel('Time (years)')
ax5.set_ylabel('Residual ($\\mu$s)')
ax5.set_title('Timing Residuals (MSP)')
ax5.grid(True, alpha=0.3)

# Plot 6: PTA sensitivity curve
ax6 = fig.add_subplot(3, 3, 6)
f_gw = np.logspace(-9, -7, 100)  # Hz
# Characteristic strain sensitivity (simplified)
h_c = 1e-14 * (f_gw / 1e-8)**(-2/3)  # nHz regime
ax6.loglog(f_gw * 1e9, h_c, 'b-', linewidth=2, label='PTA Sensitivity')
ax6.fill_between(f_gw * 1e9, h_c, 1e-10, alpha=0.3)
# SMBHB background
h_gw = 2e-15 * (f_gw / 1e-8)**(-2/3)
ax6.loglog(f_gw * 1e9, h_gw, 'r--', linewidth=2, label='SMBHB Background')
ax6.set_xlabel('Frequency (nHz)')
ax6.set_ylabel('Characteristic Strain $h_c$')
ax6.set_title('PTA GW Sensitivity')
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3)
ax6.set_xlim([1, 100])

# Plot 7: Binary pulsar orbital decay
ax7 = fig.add_subplot(3, 3, 7)
# PSR B1913+16 orbital decay
T_orbit = 7.75  # hours
Pb_dot_obs = -2.423e-12  # s/s (observed)
Pb_dot_gr = -2.403e-12  # s/s (GR prediction)

years = np.linspace(1974, 2024, 100)
cumulative_shift = (years - 1974) * Pb_dot_obs * 3.154e7 / T_orbit / 3600  # orbits

ax7.plot(years, -cumulative_shift, 'b-', linewidth=2, label='GR Prediction')
# Add data points
obs_years = np.array([1975, 1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2020])
obs_shift = (obs_years - 1974) * Pb_dot_obs * 3.154e7 / T_orbit / 3600
ax7.plot(obs_years, -obs_shift + np.random.normal(0, 0.1, len(obs_years)),
         'ro', markersize=6, label='Observations')
ax7.set_xlabel('Year')
ax7.set_ylabel('Cumulative Periastron Shift (s)')
ax7.set_title('PSR B1913+16 Orbital Decay')
ax7.legend(fontsize=8)
ax7.grid(True, alpha=0.3)

# Plot 8: Glitch behavior
ax8 = fig.add_subplot(3, 3, 8)
t_glitch = np.linspace(0, 100, 1000)  # days
nu_0 = 30.0  # Hz (spin frequency)
nu_dot = -3.7e-10  # Hz/s (spin-down)

nu = nu_0 + nu_dot * t_glitch * 86400
# Add glitch at day 50
glitch_time = 50
glitch_idx = np.argmin(np.abs(t_glitch - glitch_time))
delta_nu = 1e-6  # Hz (glitch amplitude)
nu[glitch_idx:] += delta_nu

# Recovery
tau_d = 5  # days (recovery timescale)
recovery = delta_nu * 0.1 * (1 - np.exp(-(t_glitch[glitch_idx:] - glitch_time)/tau_d))
nu[glitch_idx:] -= recovery

ax8.plot(t_glitch, (nu - nu_0)*1e6, 'b-', linewidth=1.5)
ax8.axvline(x=glitch_time, color='r', linestyle='--', alpha=0.7, label='Glitch')
ax8.set_xlabel('Time (days)')
ax8.set_ylabel('$\\Delta\\nu$ ($\\mu$Hz)')
ax8.set_title('Pulsar Glitch Event')
ax8.legend(fontsize=8)
ax8.grid(True, alpha=0.3)

# Plot 9: Dispersion measure effect
ax9 = fig.add_subplot(3, 3, 9)
freq = np.linspace(0.3, 3, 100)  # GHz
DM = 50  # pc/cm^3
delay = 4.149 * DM / freq**2  # ms

ax9.plot(freq, delay, 'b-', linewidth=2)
ax9.set_xlabel('Frequency (GHz)')
ax9.set_ylabel('Dispersion Delay (ms)')
ax9.set_title(f'Dispersion (DM = {DM} pc/cm$^3$)')
ax9.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('pulsar_timing_plot.pdf', bbox_inches='tight', dpi=150)
print(r'\begin{center}')
print(r'\includegraphics[width=\textwidth]{pulsar_timing_plot.pdf}')
print(r'\end{center}')
plt.close()

# Reference pulsar: Crab
crab = famous_pulsars['Crab (B0531+21)']
\end{pycode}

\section{Results and Analysis}

\subsection{Pulsar Population Statistics}

\begin{pycode}
# Generate population statistics table
print(r'\begin{table}[h]')
print(r'\centering')
print(r'\caption{Pulsar Population Statistics}')
print(r'\begin{tabular}{lcccc}')
print(r'\toprule')
print(r'Population & Count & $\langle P \rangle$ (s) & $\langle B \rangle$ (G) & $\langle \tau_c \rangle$ (Myr) \\')
print(r'\midrule')
for name, pop in populations.items():
    P = pop['P']
    tau_i = P / (2 * pop['Pdot']) / 3.154e7 / 1e6
    B_i = 3.2e19 * np.sqrt(P * pop['Pdot'])
    print(f"{name} & {len(P)} & {np.median(P):.3f} & {np.median(B_i):.2e} & {np.median(tau_i):.1f} \\\\")
print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')
\end{pycode}

\subsection{Famous Pulsars}

\begin{pycode}
# Famous pulsars table
print(r'\begin{table}[h]')
print(r'\centering')
print(r'\caption{Properties of Notable Pulsars}')
print(r'\begin{tabular}{lcccc}')
print(r'\toprule')
print(r'Pulsar & $P$ (ms) & $\dot{P}$ (s/s) & $B$ (G) & $\tau_c$ (yr) \\')
print(r'\midrule')
for name, props in famous_pulsars.items():
    print(f"{name} & {props['P']*1000:.2f} & {props['Pdot']:.2e} & {props['B']:.2e} & {props['tau']:.0f} \\\\")
print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')
\end{pycode}

\begin{example}[The Crab Pulsar]
The Crab pulsar (PSR B0531+21) is a young pulsar in the Crab Nebula supernova remnant:
\begin{itemize}
    \item Period: $P = $ \py{f"{crab['P']*1000:.1f}"} ms
    \item Period derivative: $\dot{P} = $ \py{f"{crab['Pdot']:.2e}"} s/s
    \item Characteristic age: $\tau_c = $ \py{f"{crab['tau']:.0f}"} years
    \item True age (SN 1054): 970 years
    \item Surface magnetic field: $B = $ \py{f"{crab['B']:.2e}"} G
    \item Spin-down luminosity: $\dot{E} = $ \py{f"{crab['Edot']:.2e}"} erg/s
\end{itemize}
\end{example}

\section{Pulsar Evolution}

\subsection{Evolutionary Tracks}

Pulsars evolve through the P-$\dot{P}$ diagram:
\begin{enumerate}
    \item \textbf{Birth}: Upper left, short period, high $\dot{P}$
    \item \textbf{Spin-down}: Move right as period increases
    \item \textbf{Death line}: Stop emitting when $B/P^2$ drops below threshold
    \item \textbf{Recycling}: Binary accretion spins up to MSP phase
\end{enumerate}

\subsection{Millisecond Pulsar Formation}

\begin{remark}[Recycling Scenario]
MSPs are ``recycled'' pulsars:
\begin{itemize}
    \item Originally normal pulsars in binary systems
    \item Accretion from companion spins them up
    \item Magnetic field buried by accreted material
    \item Final state: $P \sim 1-10$ ms, $B \sim 10^8-10^9$ G
\end{itemize}
\end{remark}

\section{Pulsar Timing Arrays}

\subsection{Gravitational Wave Detection}

Pulsar timing arrays use the correlated timing residuals of many MSPs to detect low-frequency gravitational waves:

\begin{definition}[Hellings-Downs Curve]
The angular correlation of timing residuals between pulsar pairs due to a stochastic GW background:
\begin{equation}
\Gamma(\theta) = \frac{3}{2}x\ln x - \frac{x}{4} + \frac{1}{2} + \frac{1}{2}\delta_{12}
\end{equation}
where $x = (1 - \cos\theta)/2$ and $\theta$ is the angular separation.
\end{definition}

\subsection{Sources in PTA Band}

\begin{itemize}
    \item Supermassive black hole binaries (SMBHB)
    \item Cosmic strings
    \item Primordial gravitational waves
    \item Individual continuous sources
\end{itemize}

\section{Binary Pulsars}

\subsection{Tests of General Relativity}

Binary pulsars provide precision tests of GR through post-Keplerian parameters:
\begin{itemize}
    \item Periastron advance: $\dot{\omega}$
    \item Gravitational redshift: $\gamma$
    \item Orbital decay: $\dot{P}_b$
    \item Shapiro delay: $r$, $s$
\end{itemize}

\begin{remark}[Hulse-Taylor Pulsar]
PSR B1913+16 provided the first indirect evidence for gravitational waves through its orbital decay, matching GR predictions to 0.3\%.
\end{remark}

\section{Timing Phenomena}

\subsection{Glitches}

Sudden spin-up events caused by transfer of angular momentum from the superfluid interior:
\begin{equation}
\frac{\Delta\nu}{\nu} \sim 10^{-9} - 10^{-6}
\end{equation}

\subsection{Dispersion}

Interstellar plasma delays lower frequencies:
\begin{equation}
\Delta t = 4.149 \times 10^3 \cdot \text{DM} \cdot \left(\frac{1}{\nu_1^2} - \frac{1}{\nu_2^2}\right) \text{ s}
\end{equation}
where DM is the dispersion measure in pc cm$^{-3}$ and $\nu$ is in MHz.

\section{Limitations and Extensions}

\subsection{Model Limitations}
\begin{enumerate}
    \item \textbf{Vacuum dipole}: Real magnetospheres are plasma-filled
    \item \textbf{Constant B}: Field may decay over time
    \item \textbf{Point dipole}: Higher multipoles exist
    \item \textbf{Rigid rotation}: Differential rotation possible
\end{enumerate}

\subsection{Possible Extensions}
\begin{itemize}
    \item Full magnetosphere models (force-free, MHD)
    \item Timing noise characterization
    \item Pulsar wind nebula evolution
    \item Equation of state constraints from masses
\end{itemize}

\section{Conclusion}

This analysis demonstrates the rich physics of pulsar timing:
\begin{itemize}
    \item P-$\dot{P}$ diagram reveals distinct populations and evolution
    \item Characteristic ages provide evolutionary timescales
    \item MSPs are precision clocks with $\sigma_\text{rms} \lesssim 100$ ns
    \item PTAs are sensitive to nanohertz gravitational waves
    \item Binary pulsars test GR with unprecedented precision
\end{itemize}

\section*{Further Reading}
\begin{itemize}
    \item Lorimer, D. R. \& Kramer, M. (2012). \textit{Handbook of Pulsar Astronomy}. Cambridge University Press.
    \item Manchester, R. N. et al. (2005). The Australia Telescope National Facility Pulsar Catalogue. \textit{AJ}, 129, 1993.
    \item Hobbs, G. \& Dai, S. (2017). A review of pulsar timing array gravitational wave research. \textit{National Science Review}, 4, 707.
\end{itemize}

\end{document}