Pulsar spin-down switching

A pulsar whose spin-down rate jumps between two values: the timing noise it produces, its Green–Kubo strength, the exact two-state formulas, and what timing data can and cannot tell apart.

The model

A pulsar rotates at frequency $\nu(t)$ and slows down. Its spin-down rate $\dot\nu$ is not constant: it switches abruptly between two values, set by the state of the magnetosphere. PSR B1931+24 spins down about 50% faster when its radio emission is on than when it is off (Kramer et al. 2006), and many pulsars switch between two spin-down rates that differ by about a percent (Lyne et al. 2010).

Model the state as a hidden Markov chain $y_t \in \{1, 2\}$ that leaves state 1 at rate $a$ and state 2 at rate $b$. The rates may differ, so the pulsar can spend most of its time in one state. The spin-down rate is $\dot\nu_{y_t}$, and frequency and phase (in cycles) are its integrals:

The main formulas start the chain from its stationary law, appropriate when the initial magnetospheric state is unknown. For an intermittent pulsar the radio on/off state can instead be observed at the beginning of a timing campaign; an exact known-start correction is derived below.

There are two time scales. The state persists for $1/a$ or $1/b$, from days to years.

The pulsar is timed over a span $T$ of years. The expansion parameter is $1/(\gamma T)$ with $\gamma = a + b$.

Parameters (PSR B1931+24, from Kramer et al. 2006): $\nu = 1.2289688$ Hz, $\dot\nu_{\mathrm{on}} = -16.3\times10^{-15}$ Hz/s, $\dot\nu_{\mathrm{off}} = -10.8\times10^{-15}$ Hz/s, on phases of 5 to 10 days, off phases of 25 to 35 days; mean durations $1/a = 7.5$ days and $1/b = 30$ days.

The quantity

Write $\pi = (\pi_1, \pi_2)$ for the stationary law and $\bar{\dot\nu} = \pi_1\dot\nu_1 + \pi_2\dot\nu_2$ for the averaged spin-down rate. Removing it leaves the fluctuation of the rate,

The frequency residual and the phase residual after the averaged spin-down are

The timing residual in seconds is $\Phi/\nu$. We want the mean, the variance and the full distribution of $X$ and $\Phi$, the power spectrum of the timing residual, and the residual left after the usual timing fit.

Averaging

When the state switches infinitely fast, the pulsar spins down at the constant rate $\bar{\dot\nu}$. Both residuals have mean zero, and the averaged model predicts no timing noise.

Reduction to a linear system

Let $Q$ be the generator of the chain and $G = \operatorname{diag}(g_1, g_2)$. Condition on the state over the first instant. The characteristic function of $X$ started in state $i$, $a_i(T) = \mathbb{E}_i\,e^{iuX(T)}$, then solves

This is the engine's linear system with an imaginary, constant forcing. The phase weights the rate at time $s$ by $T - s$, and the same argument gives a forcing that grows with time:

The first instant carries the largest weight, $T$, which is why the forcing at time $t$ is $iu\,t\,G$.

The first-order rule

The variance of an integral is the double integral of the covariance. For a stationary rate with covariance $C(\tau) = \operatorname{Cov}\big(\dot\nu_{y_0}, \dot\nu_{y_\tau}\big)$,

When $T$ is long against the correlation time, the factor $T - \tau$ is close to $T$ wherever $C$ is not negligible. Each spell in a state adds a nearly independent kick to the frequency, so the kicks add like a random walk:

$K$ is the Green–Kubo integral of the spin-down rate, and $Q^{\#}$ is the group inverse of the generator, as on the general page. The neglected part, $-2\int_0^\infty \tau\,C(\tau)\,d\tau$, is a constant: the initial layer.

The phase integrates the frequency again. Its variance is

To first order the frequency residual is a random walk and the phase residual is its integral. Both are Gaussian, with the characteristic functions

The rates, the switching times and the asymmetry of the chain enter only through $\bar{\dot\nu}$ and $K$.

The two-state results written out

For two states the covariance is a single exponential, $C(\tau) = \pi_1\pi_2\Delta^2 e^{-\gamma\tau}$, and every integral above can be done. The variances are exact for all $T$:

The layer constant of the frequency variance is $-2K/\gamma$. In the phase variance the first correction to the Green–Kubo form is the $T^2$ term, $-KT^2/\gamma$, and the constant is $2K/\gamma^3$.

The characteristic function of the frequency residual follows from the $2\times2$ matrix exponential. Since $\pi\cdot g = 0$, the term in $\pi\cdot(Q + sG)\mathbf 1$ drops out, and with $s = iu$

The expression is even in $D$, so the branch of the square root does not matter. For large $T$ the logarithm of the characteristic function is $\lambda(s)\,T$ plus a constant, where $\lambda(s)$ is the largest eigenvalue of $Q + sG$. Its Taylor coefficients are the cumulant rates:

The skewness is a second-order effect: it falls like $(\gamma T)^{-1/2}$ relative to the Gaussian. It vanishes when the two states are equally occupied. With unequal occupation the rare state produces bursts of extra frequency change, and the residual leans towards them.

The phase weights the rate by $T - s$, which varies slowly on the switching time. To leading order each stretch contributes its own growth rate:

This gives the phase variance $\tfrac23KT^3$ and the third cumulant above.

The power spectrum

The spin-down rate is a random telegraph signal. Its spectrum is a Lorentzian, and the phase spectrum divides it by $(2\pi f)^4$ because phase is the double integral of the rate:

The spectra are two-sided, per unit frequency, and $S_r$ is the spectrum of the timing residual in seconds. Below the corner frequency $\gamma/2\pi$ the rate spectrum is flat at $S_{\dot\nu}(0) = 2K$, and the timing residual has the $f^{-4}$ spectrum of a random walk in frequency. Above the corner it steepens to $f^{-6}$.

After a timing fit

Timing analyses fit phase, frequency and spin-down to the arrival times, so a quadratic in $t$ is removed and $\bar{\dot\nu}$ need not be known. At first order $\Phi$ is an integrated random walk with covariance

Project it onto quadratics on $[0, T]$ and average the square of what is left over the span:

Here $r$ is the phase residual after the fit. Before the fit the same average is $KT^3/6$, so the fit removes all but $1/210$ of the mean-square phase wander. The root-mean-square timing residual after the fit, in seconds, is

This last formula is the white-noise, or $\gamma T\to\infty$, limit.

Any polynomial timing fit at Green–Kubo order

The quadratic is not special. Let $P_d$ be continuous, uniformly weighted $L^2[0,T]$ projection onto all polynomials through degree $d$. For every integer $d\ge1$, the integrated-Brownian limit gives the exact all-degree coefficient

For a constant-only fit $c_0=1/30$. The cases $d=1,2,3$ give $1/420$, $1/2520$ and $1/6930$; in particular $d=2$ recovers the formula above. These are Green–Kubo coefficients: the displayed $c_d$ must not be substituted for the exact finite-rate factor $R(x)/2520$ except in the $\gamma T\to\infty$ limit.

A short proof uses the orthonormal shifted Legendre basis $\phi_n(t)=\sqrt{2n+1}\,P_n(2t-1)$ on $[0,1]$. For unit integrated Brownian motion, stochastic Fubini gives

For $n\ge2$, put $z=2s-1$. Two integrations of the Legendre recurrence give the inner integral as

Legendre orthogonality then shows that its variance is $c_{n-1}-c_n$. Summing the unremoved modes $n>d$ telescopes to $c_d$.

This extends the Legendre-projector treatment of arbitrary polynomial timing fits in Lee et al. (2012, Appendix A) to a closed coefficient for the integrated-frequency-noise kernel. An independent covariance-trace calculation for every $0\le d\le12$ agrees to relative error below $1.3\times10^{-9}$.

All-degree finite-correlation boundary cancellation

The missing finite-rate statement is also explicit at its first nonzero order. Let $R_d(x)$ be the exact continuous post-fit mean square divided by $2KT^3c_d$, with $x=\gamma T$. For every fixed integer $d\ge1$,

Thus fitting a line already cancels the relative $x^{-1}$ endpoint error, and every higher polynomial fit retains that cancellation. The result is for each fixed degree; it is not uniform when $d$ grows with $x$. A constant-only fit is exceptional:

To prove the claim, put $k_z(t)=(t-z)_+$, $H_d(z,u)=\langle(I-P_d)k_z,(I-P_d)k_u\rangle$, and

The exact dimensionless covariance integral is $\int_0^1e^{-xy}F_d(y)\,dy$. Direct shifted-Legendre projection gives

The middle identity is structural: $k_0(t)=t$ is removed whenever $d\ge1$, $k_1=0$, and symmetry of $H_d$ cancels the two diagonal derivatives. Watson's endpoint expansion then gives the displayed coefficient, since $F_d''(0)/F_d(0)=-4(2d-1)(2d+5)$.

This complements the exact quadratic formula below and keeps its scope to stationary exponential covariance and continuous uniformly weighted ordinary least squares. Independent projection quadrature for degrees one through six gives last-doubling remainder ratios from $7.89$ to $7.97$, confirming the predicted third-order residual.

Exact finite-rate quadratic

For the two-state chain the continuous, uniformly weighted quadratic least-squares calculation is exact at every switching rate. Put $x=\gamma T$. Then

The timing fit cancels the relative $1/x$ boundary correction that appears before fitting:

At the opposite extreme, $R(x)=x^2/22+O(x^3)$. A chain that barely switches produces an almost deterministic quadratic phase drift, which the fit removes.

For completeness, the proof does not assume a Gaussian telegraph signal. If $k_s(t)=(t-s)_+$ and $P_2$ is the $L^2[0,1]$ projection onto $\{1,t,t^2\}$, direct polynomial projection gives, for $0\le u\le z\le1$,

Integrating this kernel against $2e^{-x(z-u)}$ yields the displayed $R(x)$. This covariance-projection calculation follows the transfer-function treatment of pulsar timing noise in van Haasteren and Levin (2012).

It applies to a continuous, uniformly weighted ordinary least-squares fit. Irregular arrival times, measurement weights and generalized least squares require the corresponding discrete design and covariance matrices; the scalar factor $R(x)$ should not be reused in those settings.

Exact sampled-design version

There is nevertheless an exact finite-sample answer. For observation times $t_i$, let $C_{ij}$ be the covariance of the phase residuals before the timing fit. If $u=\min(t_i,t_j)$ and $v=\max(t_i,t_j)$, direct integration of the exponential rate covariance gives

Let $X_i=(1,t_i,\ldots,t_i^d)$, let $W=\operatorname{diag}(w_i)$ with positive weights summing to one, and let

The sampled weighted least-squares residual $r=(I-H)\Phi$ satisfies the exact trace identity for any degree $d$ whose sampled design has full rank:

The second equality uses $H^2=H$ and $H^\top W=WH$. More generally, the first trace formula applies to any specified linear fitting matrix, including a generalized least-squares timing fit; measurement noise is included by adding its covariance to $C$.

Thus irregular sampling and weights do not need a new homogenization approximation—they need the correct finite design. This is the covariance-matrix viewpoint advocated for correlated pulsar timing data by van Haasteren and Levin (2012).

For B1931+24 over 100 days, equally weighted grids with 17 and 257 observations give mean-square residuals $164.416$ and $139.865$ nHz$^2$ day$^2$, respectively, converging to the continuous value $138.075$. The endpoint grid has first-order quadrature error. On an eleven-point irregular weighted design, the exact trace is $154.452$; exact-path simulation gives $155.726\pm0.852$.

A known starting regime: the exact entry layer

The stationary calculation is not the conditional answer when the radio state is known at the first observation. Put $q=\Pr(y_0=1)$, $d=q-\pi_1$, and let $I_t$ indicate state 1. The two-state semigroup gives, for $s\le t$,

Consequently the nonstationary first and raw second moments of the centered spin-down rate are exactly

The second term is the exact entry layer. It is linear in the initial law, as every raw Markov expectation must be. After integrating twice, the phase mean is

For $u=\min(t_i,t_j)$, $v=\max(t_i,t_j)$ and $E_k(u)=\int_0^u s^k e^{-\gamma s}\,ds$, the raw phase second moment is

The conditional covariance is $C_q=M_q-m_qm_q^\top$. Therefore the expected fitted residual energy uses the raw second moment,

This separates stochastic timing noise conditional on the starting state from a deterministic transient mean that is not exactly quadratic and hence is not completely removed by the timing fit. When $q=\pi_1$, $d=0$ and all formulas reduce to the stationary theorem above. The derivation uses only the two-state Markov semigroup, not Gaussianity.

On the same 100-day irregular B1931+24 design, the stationary expected energy is $154.452$. Conditioning on an on-state start gives $167.131=163.217+3.914$, where the last number is the projected squared mean; an off-state start gives $151.283=151.038+0.245$. Thus using the stationary law changes this finite-span benchmark by $+8.21\%$ or $-2.05\%$, depending on the observed initial state.

Independent matrix-semigroup quadrature agrees with the closed form within $7.2\times10^{-15}$, label swapping is invariant to $2.0\times10^{-16}$, and conditioned exact-path simulations agree within Monte Carlo uncertainty. The relevance of conditioning is empirical rather than hypothetical here: Kramer et al. (2006) identify B1931+24's radio on/off state and its associated spin-down change directly.

What timing data identify

The two-state chain has four parameters: $\dot\nu_1$, $\dot\nu_2$, $a$ and $b$. To first order timing data see two numbers, the averaged spin-down $\bar{\dot\nu}$ and the Green–Kubo integral $K$. Any two chains that share them give the same variance, the same Gaussian law and the same low-frequency spectrum.

For example, keep the switching time of B1931+24 and make the two states equally likely, with the rates moved so that the average and $K$ are unchanged. Both chains have $1/\gamma = 6$ days, $\bar{\dot\nu} = -11.90\times10^{-15}$ Hz/s and $K = 2.51\times10^{-24}$ Hz$^2$/s. The frequency and timing residuals are identical in variance, and only the skewness differs:

Here the variances agree at every order, because they depend only on $K$ and $\gamma$ and the two chains share both. In general the second-order terms separate the parameters. The $T^2$ term of the phase variance gives $\gamma$.

The third cumulant gives $\Sigma = (\pi_2 - \pi_1)\Delta$. Then $\Delta^2 = \Sigma^2 + 4K\gamma$ closes the system. These corrections are smaller by $(\gamma T)^{-1/2}$ and $(\gamma T)^{-1}$, so over spans of many switching times they are hard to measure.

Quasi-periodic switching

Real switching is often close to periodic. The on and off phases of B1931+24 last 5 to 10 and 25 to 35 days, far less variable than exponential durations with the same means.

The frequency then wanders much less, and $K$ must come from an alternating renewal process. With durations of mean $m_1$, $m_2$ and standard deviations $s_1$, $s_2$, the occupation-time central limit theorem gives

Exponential durations, $s_i = m_i$, give back $\pi_1\pi_2\Delta^2/\gamma$. Strictly periodic switching, $s_1 = s_2 = 0$, gives $K = 0$: the frequency residual stays bounded and there is no random walk.

The same answer comes from the Markov machinery. Replace each phase by a chain of $k$ exponential stages, an Erlang duration with standard deviation $m/\sqrt k$. The Green–Kubo formula $-\pi\cdot(g\,Q^{\#}g)$ on the resulting cyclic chain reproduces the renewal formula exactly.

For B1931+24 with uniform durations on 5 to 10 and 25 to 35 days, $K = 5.81\times10^{-26}$ Hz$^2$/s, about 43 times smaller than the Markov value. Timing noise scales as $\sqrt K$, so the quasi-periodic residuals are about 6.6 times smaller than the Markov ones.

Worked example: PSR B1931+24

Kramer et al. (2006) give $\nu = 1.2289688061$ Hz, $\dot\nu_{\mathrm{on}} = -16.3(4)\times10^{-15}$ Hz/s, $\dot\nu_{\mathrm{off}} = -10.8(2)\times10^{-15}$ Hz/s, an active duty cycle of 19(5)%, on phases of 5 to 10 days and off phases of 25 to 35 days. Take state 1 as on, with mean durations 7.5 and 30 days. Then

The averaged rate is within the duty-cycle uncertainty of the published long-term value, $-12.2488\times10^{-15}$ Hz/s. Over spans of months to years:

The frequency residual grows like $\sqrt T$ and the timing residual like $T^{3/2}$. The Green–Kubo column is within about 3% of the exact value at 100 days and 0.3% at 1000 days. The quadratic fit absorbs most of the wander, leaving milliseconds to a fraction of a second.

The quasi-periodic column uses the renewal $K$ and is the more realistic of the two. The skewness is negative because the rare on state spins the pulsar down faster.

Lyne et al. (2010) list slower switchers. For PSR B1828−11 they give $\nu = 2.469$ Hz, $\dot\nu = -365.68\times10^{-15}$ Hz/s, a peak-to-peak change of 0.71% in $\dot\nu$ and a fluctuation frequency of 0.73 per year, a cycle of about 500 days.

The split of the cycle between the two states is not in their table; take one third in the faster state as an illustration, with mean durations 167 and 334 days, so $\dot\nu = (-367.41,\ -364.81)\times10^{-15}$ Hz/s and $1/\gamma = 111$ days. Then $K = 1.44\times10^{-23}$ Hz$^2$/s, and:

Here a switching time of months meets a timing span of years, and $\gamma T$ is between 16 and 66. Lyne et al. find a narrow periodic peak in the spin-down spectrum of B1828−11, so the Markov figures overstate its random-walk component in the same way as for B1931+24.

Regular switching in seventeen pulsars

Keith and Nițu (2023) published spin-down rate series for 17 Jodrell Bank pulsars, each about 23.6 years long. The series are smoothed by Gaussian-process regression of the timing residuals.

After removing a linear trend, the autocovariance of each series was fitted with the two-level stage model: $k$ exponential stages per level, which gives durations with standard deviation $m/\sqrt k$. Only lags of 300 days and more were used, beyond the reach of the smoothing.

The stages and durations below come from the fit to all lags. The last column is the least-squares misfit of the best Markov model ($k = 1$) divided by that of the best $k$, using lags of 300 days and more:

For 12 of the 17 pulsars the Markov fit is at least 5 times worse. Their best fits have 6 to 64 stages per level, so the durations vary by 12 to 41 percent of their means, not the 100 percent of exponential durations. The Markov misfit grows or stays the same when the shortest lags are dropped, so the regularity is not produced by the smoothing.

For these pulsars the renewal formula applies, and $K$ is about $1/k$ of the Markov value with the same levels and mean durations. Their random-walk timing noise is correspondingly smaller than a Markov switching model predicts.

The remaining five pulsars (B0740−28, B0919+06, B1822−09, B1839+09 and B1903+07) are fitted nearly as well by Markov switching. For B0740−28 and B1839+09 the observed $K$ is too large for any two-level fit until a slow Ornstein–Uhlenbeck component is added. B1822−09 has sharp spikes rather than two levels.

Twenty-three years hold too few independent blocks to measure $K$ directly to better than a factor of about 4, so the evidence for regularity comes from the shape of the autocovariance. The analysis scripts and a log of the results are in research/pulsar.

Checks

All checks use the B1931+24 parameters, with the span $T$ standing in for the speed of switching. The frequency variance against the matrix exponential of the Van Loan block generator, which gives the double integral exactly (values in nHz$^2$):

The closed form agrees to rounding. The Green–Kubo error is $1/(\gamma T)$, a first-order remainder: it falls threefold with each threefold increase in $T$.

The phase variance against a numerical solution of the moment equations, whose second moment is $\pi\cdot m_2(T)$ (values in nHz$^2$ day$^2$):

The $T^3$ term alone has relative error $3/(2\gamma T)$. Adding the $T^2$ term leaves $-3/(\gamma T)^3$, as the closed form predicts.

The characteristic function of the standardized frequency residual, largest error over $0 \le u \le 4$, against the matrix exponential of $T(Q + iuG)$:

The Gaussian error halves each time $\gamma T$ quadruples, as the skewness does. With the layer constant in the variance and the third cumulant $6\lambda_3T$, the error falls fourfold, a second-order remainder.

The characteristic function of the standardized phase residual, against a numerical solution of the engine system $a' = (Q + iu\,t\,G)\,a$:

The same orders hold. The largest-eigenvalue expansion is checked directly too: finite differences of the eigenvalue of $Q + sG$ give $\lambda_2$ and $\lambda_3$ to five and four digits, and the group inverse gives the same $K$.

Monte Carlo with no time discretization: the switching times are drawn exactly, and on each path $X$ and $\Phi$ are integrated in closed form between switches, 400,000 paths per span (units as above):

The mean-square phase residual after a least-squares quadratic fit to 1001 equally spaced samples, from 4,000 exact paths, against the first-order formula:

The rate spectrum computed as the numerical Fourier transform of the matrix-exponential covariance $\pi\cdot\big(g\,e^{\tau Q}g\big)$ matches the Lorentzian to eight digits at every frequency tested. For quasi-periodic switching, the Erlang chain with $k = (27, 108)$ stages reproduces the renewal $K$ to rounding. A Monte Carlo of 20,000 paths with uniform durations over 20,000 days gives $\operatorname{Var}X/2T = 0.00503 ± 0.00005$ against the renewal value 0.00502 (nHz$^2$/day).

Prior art

That timing noise comes from switching between spin-down rates is the finding of Kramer et al. (2006) for B1931+24 and of Lyne et al. (2010) for many pulsars. Lyne et al. simulated it directly: two rates with dithered switching times, integrated twice, sampled, with a quadratic removed. Nothing here claims that observation.

The random walk in frequency, with its $f^{-4}$ residual spectrum, is one of the classical timing-noise models of Boynton et al. (1972) and Cordes and Helfand (1980). Hobbs, Lyne and Kramer (2010) analyse timing irregularities across hundreds of pulsars, and van Haasteren and Levin (2012) treat timing noise as a Gaussian process with a red spectrum. The telegraph covariance, its Lorentzian spectrum and the central limit theorem for occupation times of alternating renewal processes are classical.

The fast-switching reading adds the following to these facts. The strength of the random walk is a Green–Kubo integral, as in Green (1954) and Kubo (1957). The two-state variances, characteristic function and cumulant rates are written out for finite spans, with the layer constant.

The post-fit residual has a closed first-order form. To first order the timing data identify $\bar{\dot\nu}$ and $K$, not the individual rates and switching times, and quasi-periodic switching enters only through a smaller $K$.

Certificate: verify_pulsars.py.