Numerical Pulsar Timing Model

• The Numerical Pulsar Timing Model is a quantitative framework that decomposes deterministic spin/orbital components and stochastic propagation effects to predict pulse arrival times.
• It employs multi-frequency least squares fitting to isolate dispersion, scattering, and white noise, enhancing the precision of astrophysical parameter estimates.
• The model mitigates interstellar scattering, radiometer noise, and intrinsic jitter, providing robust error budgeting crucial for pulsar timing arrays and gravitational wave research.

A numerical pulsar timing model is a quantitative framework for modeling, error budgeting, and mitigating the various stochastic and deterministic effects that perturb the observed pulse times of arrival (TOAs) from their ideal, infinite-frequency values. Such a model explicitly includes multi-frequency plasma propagation effects (e.g., dispersion, scattering), radiometer noise, intrinsic pulse phase jitter, diffractive scintillation, and achromatic stochastic contributions, providing formulae that allow the separation and minimization of statistical and systematic uncertainties. These models are foundational for precision pulsar timing applications including pulsar timing arrays (PTAs) targeting nanohertz gravitational wave detection, astrophysical parameter inference, and tests of fundamental physics.

1. Timing Equation and Model Structure

The core of the numerical timing model is an additive decomposition of the measured TOA at radio frequency ν, transformed to the solar system barycenter, as

$$t(\nu) = t_{(\infty)} + t_{\text{DM}}(\nu) + t_C(\nu) + \epsilon(\nu)$$

where:

• $t_{(\infty)}$: Infinite-frequency TOA, incorporating deterministic components from pulsar spin and orbital motion, plus intrinsic red ("achromatic") timing noise.
• $t_\text{DM}(\nu)$: Dispersive delay, scaling as $\nu^{-2}$ (Eq. 5: $t_\text{DM} = a_\text{DM} \text{DM}/\nu^2$, with $a_\text{DM} \simeq 4.15$ ms for GHz and pc cm⁻³ units).
• $t_C(\nu)$: Strongly chromatic delay from interstellar scattering, typically $\propto \nu^{-4.4}$ for a Kolmogorov spectrum (Eq. 6).
• $\epsilon(\nu)$: Wideband white-noise term encapsulating radiometer noise, pulse-phase jitter, and the effect of finite diffractive scintles.

In practical timing analysis, TOAs are collected over broadband, multi-channel observations, and a least-squares fitting procedure is employed across frequency to determine $t_{(\infty)}$, the dispersion measure (DM) correction, and potentially an additional chromatic coefficient for scattering.

2. White-Noise Error Components

The total white-noise TOA uncertainty ($\Delta t_\text{WHITE}$) arises via quadrature sum (Eq. 4):

$$\Delta t_\text{WHITE} = \sqrt{ (\Delta t_\text{RN})^2 + (\Delta t_\text{J})^2 + (\Delta t_\text{DISS})^2 }$$

with:

• Radiometer noise: For $N$ averaged pulses of effective width $W_\text{eff}$,

$$\Delta t_\text{RN} \simeq \frac{W_\text{eff}}{\text{SNR}} = \frac{W_\text{eff}}{\text{SNR}_1 \sqrt{N}}$$

or

$$\Delta t_\text{RN} \simeq \frac{S_\text{sys}}{S_\text{peak}} \frac{1}{ \sqrt{2B N} }$$

where $B$ is the observing bandwidth, $S_\text{sys}$ the system flux density, and $S_\text{peak}$ the pulsar's peak flux density.

• Intrinsic pulse jitter: For a Gaussian pulse with single-pulse phase standard deviation $f_J$, width $W_i$, intensity modulation index $m_I$,

$$\Delta t_\text{J} \simeq \frac{f_J W_i \sqrt{1 + m_I^2}} {2\sqrt{2N_i \ln 2}}$$

where $N_i$ is the number of statistically independent pulses. The jitter term often dominates at high SNR; observed values of $f_J$ are $0.3$–$0.5$.

• Scintillation white noise: The number of independent scintles in the observation $N_\text{iss}$ sets the white-noise floor due to the finite number of sampled interstellar diffractive patterns:

$$\Delta t_\text{DISS} \simeq \frac{\tau_d}{\sqrt{N_\text{iss}}}$$

where $\tau_d$ is the diffractive pulse broadening timescale.

3. Chromatic Propagation Effects and Error Modeling

Propagation through the ionized ISM introduces two dominant chromatic delays:

• Dispersion: $t_\text{DM}(\nu) = a_\text{DM} \text{DM}/\nu^2$, with DM varying in time and stochasticity imparted by ISM turbulence (Eq. 5).
• Scattering: $t_C(\nu) \simeq a_C \nu^{-X}$, with $X \approx 4.4$ for a Kolmogorov spectrum (Eq. 6), causing deterministic pulse broadening and stochastic arrival-time jitter via variable ISM conditions.

At sufficiently high DM or low $\nu$, $t_C$ may dominate over the dispersive term.

Additional chromaticity may arise from refractive angle-of-arrival scintillation or other refractive delays, which are often strongly frequency-dependent.

4. Multi-Frequency Fitting and Parameter Estimation

At each observing epoch, the set of multi-channel TOAs is modeled as:

$$t(\nu) = t_{(\infty)} + a_\text{DM} \nu^{-2} + a_C \nu^{-X} + \epsilon(\nu)$$

A design matrix is constructed with columns for $1$, $\nu^{-2}$, and $\nu^{-X}$, and a weighted least squares fit is performed. The frequency leverage is quantified via weighted sums

$$s_p = \sum_k w_k\, \nu_k^{-p} \quad (p = 0,2,4,\ldots)$$

The uncertainty in $t_{(\infty)}$ in a two-parameter ($t_{(\infty)}, a_\text{DM}$) fit is

$$\sigma_{t_{(\infty)}}^2 = \frac{s_4}{s_0 s_4 - s_2^2}$$

Neglecting significant scattering leads to systematic biases in $t_{(\infty)}$, characterized by a residual

$$\delta t_{(\infty)} \simeq - a_C \cdot R_{t_{(\infty)}}$$

where $R_{t_{(\infty)}}$ quantifies the bandpass lever arm and $X$ index dependence. Aggressive mitigation fits for additional chromatic parameters but increases statistical errors per additional parameter.

5. Impact and Mitigation of Interstellar Scattering

Scattering introduces two effects:

• Deterministic pulse broadening: Arrival time is shifted by the mean of the pulse broadening function $\tau_{d, \mathrm{mean}}\propto \nu^{-4.4}$ (Kolmogorov spectrum).
• Diffractive scintillation: Adds a white noise term due to the finite sampling of scintles:

$$\Delta t_\text{DISS} = \frac{\tau_d}{\sqrt{N_\text{iss}}}$$

Scattering bias can be reduced by:

• Observing at higher frequencies where $\tau_d$ is reduced.
• Applying dynamic spectrum analysis to estimate $\Delta\nu_d$ (the scintillation bandwidth), using the relation $2\pi \tau_d \Delta\nu_d \simeq C_1$ (with $C_1 \sim 1$).
• Subtracting $\tau_d$ from measured TOAs on a per-epoch basis, though geometric and anisotropic uncertainties limit this in practice.

6. Precision Limits: DM, Bandwidth, and Red Noise

The DM determines the feasibility of high-precision timing at a given frequency via

$t_\text{DM} = a_\text{DM}~\frac{\text{DM}}{\nu^2}$

$$\delta t_\text{DM} = a_\text{DM}~\frac{\delta \text{DM}}{\nu^2}$$

Sub-microsecond precision at $1.4$ GHz is only attainable for $\text{DM} \lesssim 30$ pc cm⁻³ unless scattering is aggressively mitigated. The separation of dispersion and scattering delays requires fractional bandwidth coverage of nearly an octave.

After all propagation effects are removed, the ultimate limit is set by "red" timing noise (intrinsic stochasticity in the spin rate). Unlike plasma effects, this is achromatic, cannot be reduced via frequency leverage, and will set a floor to $t_{(\infty)}$ measurement precision.

7. Synthesis: Full Measurement Model and Application

The overall TOA perturbation is written as

$$\Delta t_\text{TOA} = \Delta t_\text{WHITE} + \Delta t_\text{SLOW}$$

$$\Delta t_\text{WHITE} = \Delta t_\text{J} + \Delta t_\text{RN} + \Delta t_\text{DISS}$$

$$\Delta t_\text{SLOW} = \Delta t_\text{DM} + \Delta t_\text{PBF} + \text{(other chromatic terms)}$$

Operationally, multi-frequency least squares fits are performed to estimate $t_{(\infty)}$, DM, and scattering parameters at each epoch; joint time series are then used in further “multi-epoch fitting” to infer parameters of astrophysical interest (e.g., ephemerides, binary motion, GW signatures).

The model defines a robust error budget and provides formulaic prescriptions [Eqs. (2), (3), (4), (5), (6), (7), (10), (11)] to quantify all major sources of timing error.

Application of this measurement model is critical for PTAs in the detection of gravitational waves, where residual systematic errors from either scattering or incompletely characterized red noise can obscure or mimic GW signatures at the $\sim$100 ns precision target.

8. Summary Table: Key Timing Model Components

Term / Effect	Frequency Scaling	Formula / Reference
Dispersion delay	$\nu^{-2}$	$t_\text{DM}(\nu)$, Eq. (5)
Scattering delay	$\nu^{-4.4}$	$t_C(\nu)$, Eq. (6)
Radiometer noise	Weakly chromatic	$\Delta t_\text{RN}$, Eq. (2)/(2a)
Pulse-phase jitter	Weakly chromatic	$\Delta t_\text{J}$, Eq. (3)
Diffractive scintle noise	$\propto \tau_d/\sqrt{N_\text{iss}}$	Eq. (7)
Achromatic red noise	None (spin)	Intrinsic; sets ultimate floor

9. Significance and Limitations

The measurement model unifies deterministic spin/orbital timing with stochastic and frequency-dependent plasma and emission effects, enabling aggressive mitigation strategies for propagation errors. It demonstrates that, for high-SNR millisecond pulsars, precision is fundamentally limited by irreducible interstellar scattering or intrinsic pulse jitter, not radiometer noise. The model also quantifies the conditions under which bandwidth increases yield diminishing returns or even degrade timing if chromatic effects are not simultaneously fitted and subtracted.

At current and next-generation sensitivity levels (Arecibo, FAST, SKA), robust error budgeting and frequency leverage will determine which pulsars and which DMs are admissible for sub-microsecond or sub-100 ns timing, an essential prerequisite for nanohertz gravitational wave detection and high-precision neutron star physics.

References: All claims, equations, and methodology above are taken verbatim from "A Measurement Model for Precision Pulsar Timing" (Cordes et al., 2010).