astrophysics - Pulsars: How do astronomers measure minute changes in period (~picoseconds per year)?

Let us suppose that the pulsar is spinning down at a uniform rate. So it has a period $P$ and a rate of change of period $dP/dt$ that is positive and constant (in practice there are also second, third, fourth etc. derivatives to worry about, but this doesn't change the principle of my answer).

Now let's assume you can measure the period very accurately - say you look at the pulsar today and measure its radio signals for a few hours, do a Fourier transform of the signal and get a nice big peak with a period of 0.1 seconds (for example).

With that period, you can "fold" the data to create an average pulse profile. This pulse profile can then be cross-correlated with subsequent measurements of the pulse to determine an offset between the predicted time of "phase zero" in the profile, calculated using the 0.1 s period, and the actual time of phase zero. This is often called an "O-C" curve or a residuals curve.

If you have the correct period and $dP/dt=0$, then the residuals will scatter randomly around zero with no trend as you perform later and later observations (see plot (a) from Lorimer & Kramer 2005, The Handbook of Pulsar Astronomy). If the initial period was in error, then the residuals would immediately begin to depart from zero on a linear trend.

If however, you have the period correct, but $dP/dt$ is positive, then the residuals curve will be in the form of a parabola (see plot (b)).

If you have second, third etc. derivatives in the period, then this will affect the shape of the residuals curve correspondingly.

The residuals curve is modelled to estimate the size of the derivatives of $P$. The reason that $dP/dt$ can be measured so precisely is that pulsars spin fast and have repeatable pulse shapes, so changes in the phase of the pulse quickly become apparent and can be tracked over many years.

Mathematically it works something like this. The phase $phi(t)$ is given by
$$phi(t) simeq phi_0 + 2pi frac{Delta t}{nP} - frac{2pi}{2}frac{(Delta t)^2}{nP^2} frac{dP}{dt} + ...,$$
where $phi_0$ is an arbitrary phase zero, $Delta t$ is the time between the first and last observation and $n$ is the integer number of full turns the pulsar has made during that time. If the period is approximately correct, then $n = int(Delta t/P)$.

The "residual curve" would be given by
$$phi_0 - 2pifrac{Delta t}{nP} -phi(t) simeq frac{2pi}{2}frac{(Delta t)^2}{nP^2} frac{dP}{dt} + ...,$$

For example, if the period of a $P sim 0.01$ second pulsar changed by a picosecond in a year, then there would be an accumulated residual of almost $10^{-4}$ seconds after 1 year of observation. Depending on how "sharp" the pulse is, then this shift of about 1% in the phase of the pulse might be detectable.

Perhaps needless to say, but there are a host of small effects and corrections to make in order to get this very high precision timing. You need to know exactly how the Earth is moving in its orbit. The proper motion of the pulsar on the sky also has an effect. These and more can be found in Lorimer and Kramer's book, but there is also a summary here.