Welcome!

This is the beginning of your Camp Pulsar interactive tutorials. We've created tutorials covering the entire PTA process, from an initial solution for a pulsar through to astrophysical interpretation of a gravitaitonal wave result.

We encourage you to team up with other students, especially those from other institutions or PTAs, as you work through these tutorials. There is a choose your own adventure aspect to this -- there's no way we can unpack every subtlety of PTA science in a few short days. Therefore, you have options: you can either work through all the main path tutorials to get a high level picture of how PTA science works, or you can camp on the tutorials that are of most interest to you.

Throughout the week, we'll have "counselors" available to help answer questions. These are postdocs and faculty members who work on these topics.

So, without futher ado, let's start building a PTA.

Timing pulsars:

This tutorial is adapted from existing NANOGrav Timing Tutorials created by Sofia Sosa and Michael Lam.

An introduction

Before we can do precision pulsar timing array science, we have to develop an initial timing solution. In fact, the 0th step in this process is to search for pulsars. But in the interest of time, we will assume that has already been done for us.

When a new pulsar is discovered in a pulsar search, we generally have minimal information about it -- just its general posiion on the sky, its spin period, and its dispersion measure.

All of the rest of the pulsar details we'll find by "solving" the pulsar. This is an iterative process that continues at some level forever, as we get more observations of the pulsar and refine our timing solution. But, the initial pulsar solving proces is perhaps the most satisfying activity in pulsar astronomy, so we'll focus in on it.

Unfortunately, our actual PTA pulsars tend to be a bit complicated, so we will start with solving a fake pulsar then move on to refining a timing solution for a real PTA pulsar.

What is pulsar timing?

Pulsar timing precisely predicts when pulses from a pulsar will arrive at our telescope. We define a "timing model" that includes spin parameters of the pulsar (how the pulsar spins), astrometric parameters of the pulsar (how the pulsar moves across the sky), binary parameters (where applicable -- how the pulsar orbits its companion), and how the interstellar medium affects the pulse as it propagates from the pulsar to Earth.

The cartoon below attempts to illustrated all of these features: pulsar spin, binary motion, motion throughout the galaxy of the pulsar, motion of Earth around the Sun and through the galaxy, and the interstellar medium. The relative motion of the Earth and the pulsar is illustrated with the big black arrows. The two different colored (red and blue) curves demonstrate the effects of the interstellar medium: a single broadband pulse then spreads out in frequencies due to dispersion, such that longer wavelengths (red) arrive later than shorter wavelengths (blue).

If our timing model isn't accurate, we won't be able to predict the arrival of pulses. This is illustrated in the figure below. A good model (like the top panel) predicts the pulse arrival times exactly, while a bad model (like the middle panel) will not. Our mission is to improve our timing model, by adding and modifying parameters until our predictions match the more complex physical reality of the pulsar (as in the bottom panel).

The difference between the observed time of arrival (TOA) of a pulse and the predicted arrival time is known as a timing residual, and plots of residuals help us visualize whether or not we have a good model or not.

In this exercise, you will develop an initial pulsar timing solution, and along the way discover how incorrect or missing parameters change the shape of your residual plot. This exercise is written in PINT, but the same techniques can be used in TEMPO2 or TEMPO -- it's just a matter of using the right syntax. We've chosen to present the PINT process to keep the tutorials consistently in a Python framework, but depending on what PTA you're working with and your goal, TEMPO2 may be more appropriate in the long-run.

What information do we have about a newly discovered pulsar?

When we first discover a pulsar, we only have a few pieces of information about it (the spin period, the DM, and the position). This is enough information to observe the pulsar reliably, but it's not enough fit parameters to adequately predict times of arrival. We use this very basic set of information to observe the pulsar enought to fit other parameters.

Suppose we have just discovered a pulsar, J1910-3020, with a spin period of $P = 15 ms$. This is not a real pulsar -- this is a simple simulated pulsar designed to help us learn the process of timing pulsars. But, today, we're going to practice fitting it like it's real new pulsar.

The activity below is not an exact analogy to this process -- it's very oversimplified because we're generating simulated TOAs as we go along. However, solving a real pulsar is similar process (and could very easily eat an entire week).

# This cell imports the software we'll use in our analysis
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from io import StringIO
import pint.fitter
from pint.models import get_model
from pint.toa import get_TOAs
from pint.residuals import Residuals
from pint.simulation import make_fake_toas_uniform
import pint.logging
# Advanced User Note: if you would like to see the DEBUG and INFO statements associated with PINT functions, 
# comment out the following line. These are not necessary for learners, but can be useful for experts.
pint.logging.setup("ERROR")

Now, we'll begin to use PINT to analyze this dataset.

Ordinarily, our first task is to load in the relevant timing model and TOAs. However, in the interest of simplicity, we're going to gradually build a very simple model for this theoretical pulsar and simulate a few TOAs along the way.

We'll start with three functions:

• make_model() will generate a simple PINT timing model object. This will contain the same information as a timing par file and can be used to save a par file directly from our notebook. By default, the only quantities in the model are the pulsar name and an epoch -- that means you'll need to add in a multi-line string containing each additional parameter you intend to use each time you add a parameter.

• make_TOAs() will take a model object and simulate TOAs using PINT.

• do_timing() will take the TOAs and the model, time the pulsar, and by default produce a plot of the residuals. It returns the "best fit" model if fit=True is provided as an argument.

# Helpful Functions for us to use.
def make_model(string):
    par_sim = """MODE 1
PSR      J1910-3020
PEPOCH    50000.0
"""
    par_sim += string
    return get_model(StringIO(par_sim))
def make_TOAs(model):
    t = make_fake_toas_uniform(startMJD=53000,
                               endMJD=56000,
                               ntoas=200,
                               fuzz = 1 * u.d,
                               model=model,
                               error=0.1 * u.us,
                               add_noise=True,
                               name="fake",
                               include_bipm=True
                              )
    return t
def do_timing(t, m, fit=False, plot=True):
    f = pint.fitter.Fitter.auto(t, m)
    if fit:
        f.fit_toas()
    rs = f.resids.time_resids.to(u.us).value
    xt = t.get_mjds()
    err = t.get_errors().to(u.us).value
    if plot:
        plt.figure()
        plt.errorbar(xt, rs, yerr=err, fmt=".")
        plt.xlabel("MJD")
        plt.ylabel(r"Residual ($\mu$s)")
        plt.grid()
        plt.show()
    if fit:
        return f.model

Our next step is to set up our model and simulate TOAs.

According to our discovery observations, we know the following information about J1910-3020:

• RA: 19:10:00
• Dec: -30:20:00
• Spin Period: 15 ms
• Dispersion Measure: 10.3 pc/cm$^3$

Instead of fitting for period, it is generally better to fit for spin frequency, because it is more fundamental. Therefore, we'll use $f = 66.66\bar{6} \text{Hz}$ instead of $P = 15 \text{ms}$.

Run the following cell to generate a small set of TOAs, with this model:

m = make_model('''
F0       66.66666667
RAJ 19:20:00
DECJ -30:20:00
DM 10.31
''') 
t = make_TOAs(m)

This is not the way we usually create pulsar models and TOAs (even fake TOAs)!

Instead, we usually save pulsar parameters in a text file called a "par" file (short for parameters) and TOAs in a text file called a "tim" file. We won't be using this file-based framework in this tutorial, but we will in the future, so you can try your hand at loading par and tim files now.

In the cells below, try using PINT to load the same fake TOAs and fake model with the functions get_model and get_TOAs A typical call looks like my_timing_model = get_model(path-to-your-file) or my_toas = get_TOAs(path-to-your-file).

Our fake discovery parameters are saved in a file called fakepsr_discovery.par and our fake TOAs are saved in a file called faketoas.tim.

m_file = 
t_file = 

Initial Fit and Understanding F0

Let's begin by assessing how our initial solution looks using the function do_timing defined above.

At this time, scroll back up to the cell containing our function definitions and read through do_timing and answer the following questions.

Beginning Users: What input does this function take? What output does it have? What are the optional arguments for this function?

Advanced Users: What steps does this function include? How do these steps correlate with other pulsar timing software you've used?

Everyone: run do_timing in the cell below. Before you do so, ask yourself what a "good" result would be. How would you know if the timing model is a good fit to the TOAs?

You may record your answers in the cell below. Select "Markdown" instead of "Code" from the dropdown box at the top of the page to type text freely in a cell.

do_timing(t,m)

What do you notice about the residuals? Do you think this is a reasonable fit? Why or why not? Again, record your answer in the cell below.

Now, what would happen if we were to slightly change the value for our parameters?

In particular, what do you think would happen if (without changing anything about the TOAs), we changed the spin frequency of the pulsar by 10 nHz. This corresponds to adding a 1 to the final digit in F0.

Let's try it: Make a new model that has the same RA, Dec, and DM, but has a 10 nHz larger spin frequency. Repeat the timing with this new model and the same TOAs. How does your result change?

m_10nHz = 

do_timing(t,m_10nHz)

What does the residuals plot look like? Does this suggest a good or poor fit? Why do you think this is happening?

Now, repeat the above exercise, but add 1 nHz to the original value instead of 10 nHz.

m_1nHz = 

do_timing(t,m_1nHz)

How is this result similar to the result when you added 10 nHz? How is it different?

Often, our first value for the spin period isn't quite right, and we need to fit our timing model to the data. Let's learn how to fit the do that. To do fit a parameter, we add a "1" after the value of the parameter in our model or par file. This is a very old convention, which tells pulsar timing software that the parameter should be fit for. A "0" or nothing after the value of a parameter means means to hold that value fixed.

If your model or par file has the line \ F0 66.666666671\ it will hold the value of the spin frequency constant at 66.666666671. If it has the line\ F0 66.666666671 1\ it will fit for the spin frequency using an initial value of 66.666666671.

For PINT at least, the number of spaces before your 1 is not important.

Our simple timing function can do the fit if we provide the argument we can supply the fit=True argument to the do_timing() function, like below. Note also now we will grab the resulting model so we can examine its components.

m_1nHz_fit  = 

result = do_timing(t, m_1nHz_fit, fit=True)

How does that look? Is it a better or worse fit than before?

We would also like to know what value the fitter decided on for F0. In PINT, we can view the entire model by typing the name of the model (in this case result) or the value for a specific parameter by typing the name of the model . the name of teh parameter, e.g., result.F0.

In the cells below, try printing the total model and the value for F0 after the fit.

Understanding F1

We know that since pulsars are emitting radiation, they are losing energy and thus are slowing down. The spin frequency derivative, $\dot{f}$, is written in the par file as F1. The units on $\dot{f}$ are typically Hertz per seconds, that is, how many Hertz the pulsar is slowing down every second. As you might expect, it's a very tiny number.

Our initial model didn't contain any F1 parameter, but it's likely that if this were a real, newly-discovered pulsar, fitting F1 would be a critically important step.

Let's create a new set of simulated TOAs that include F1. NOTE: this cell will redefine t, m, and result. Be careful!

m = make_model('''
F0       66.66666667
F1       -1e-19
RAJ 19:20:00
DECJ -30:20:00
DM 10.31
''') 
t = make_TOAs(m)
result = do_timing(t, m)

Now suppose we have a slightly wrong initial guess for F1. What do you think will happen?

Let's changing F1 by $0.1\times 10^{-19}$ and see what happens.

m_f1change = 

result = do_timing(t, m_f1change)

What does your result look like? Does an incorrect $\dot{f}$ have the same effect as an incorrect $f$?

Next, try fitting for F0 and F1 at the same time (both with slightly wrong starting values), and print the values of the parameters. How do they compare to the true values? Is the fit good?

m_f0f1fit = 

result = do_timing(t, m_f0f1fit, fit=True)

Understanding Astrometry

We've now understood spin parameters: how fast is the pulsar spinning and how fast is it slowing down. Next, we'll consider astrometric parameters -- parameters that have to do with where the pulsar is and where the Earth is.

Understanding Position

At present, we're assuming we know perfectly where the pulsar is located. What if there is a small discrepancy in our position?

Try definining a model with a 1 milliarcsecond error in declination. (00:00:0.001).

m_dec = 

do_timing(t,m_dec)

What does your residual plot look like? Is this a good fit? How do you know?

You can already guess what's next: start with a milliarcsecond deviation and then fit.

m_dec = 
do_timing(t, m_dec, fit = True)

How do your results look now?

We'll skip over explicitly examining what an error in RA looks like -- it's very similar to an error in Dec.

Understanding Parallax

The pulsar is at some distance from us, which produces a parallax. The units on the parallax are in milliarcseconds (mas). A pulsar with a parallax of 1 mas is at a distance of 1 kpc (kiloparsec). Let's put a pulsar at 0.5 kpc away, so it will have a parallax of 2.0 mas.

PX 2.0

Then let's see what happens when the pulsar is plotted with the wrong parallax, let's say 2.5 mas.

Instead of making two separate cells, let's condense our work into one cell. Now just call make_model() with the appropriate values, then make_TOAs(), then make_model() again with the slightly adjusted parallax before plotting the residuals with do_timing().

An important note: Quite often, we're not able to measure parallax for a pulsar in a statistically significant way. Different PTAs have different conventions regarding this issue -- make sure to check whether you should be fitting for parallax when you start a new IPTA timing project.

m_px = 
t_px = 


do_timing(t_px, m_px)

What do you notice about the fit? What is the timescale of the variation?

It turns out that the "why" of this timescale is a little bit more complicated than for an incorrect pulsar position. It has to do with the fact that the pulsar emission arriving at the solar system is not a plane wave but there is some curvature to it.

Yet again, let's fit all of the parameters. Then print out the parallax value. How does the precision on this value compare to some of the others you've examined?

Understanding Proper Motion

Let's add one more set of astrometric terms, the proper motion in the RA and DEC directions. These are given in units of milliarcseconds per year (mas/yr). In the model, we write PMRA and PMDEC. Set the true PMRA to 50.0 and the PMDEC to 0.0 just for ease. Then simulate your TOAs, make a new model file with PMRA set to 50.1, and plot the results.

What do you notice about the residuals now? What is happening to the timescale of the variation? Any ideas as to why the structure of the residuals looks the way it does?

m_pm = 
do_timing(t_pm, m_pm)

Once again, let's do the fit! Check to see if the values of the proper motion make sense.

Understanding Binary Parameters

In the interest of time, this tutorial is not going to cover fitting for binary parameters. However, many (more than half) PTA pulsars are in binary systems, so most par files you encounter will include a binary model. There are several different parameterizations -- different ways to write down the relevant orbital parameters. This will generally include some way to parameterize the binary period, the semi-major axis, the eccentricity, and the $T0$ value for each pulsar. Depending on the properties of the system and how well modeled it is, the binary model may also include the inclination angle, the mass of the other object, and other post-Keplerian parameters, things that are inherently relativistic.

One of the best quick reference guides for binary pulsar models is (oddly) found on the TEMPO sourceforge page here: https://tempo.sourceforge.net/ref_man_sections/binary.txt

The most common models for PTA uses are ELL1, ELL1H, DD, DDK, and BT. Each has its own benefits and drawbacks; ELL1 is a probably starting point for most systems. It's appropriate for low-eccentricity systems and uses the Laplace parameters in place of e and $\Omega$.

TEMPO2 also has an option called T2 that will allow the computer to help you choose a binary model -- if you're working on an IPTA project, first check if this is an appropriate model.

Understanding Dispersion Measure

For most initial pulsar solutions (and for some applications longer-term), the value of the dispersion measure is treated as fixed. For PTA science, this doesn't work. Time-dependent variations in the interstellar medium are an important source of noise in PTA data, and mitigating these effects is an active area of study.

There are many different solutions. At present, the IPTA DR3 process is using a hybrid approach, where DM and its first and second derivatives are fit during the timing process. Later, during the noise analysis process, a Gaussian Process model is created for the DM time series.

We'll talk about this more in the Noise Tutorial.

Saving our final solution

One final, simple task we might want to undertake is saving our new par and tim files. Fortunately, PINT has built-in functions to allow this.

To write out the par file, we type result.write_parfile('desired_file_name.par'), where result is the name of our fitted model. To write out the tim file, we type t.write_TOA_file('desired_file_name.tim') where t is the name of our TOA object.

Let's practice below:

Optional: Derived Parameters

There are several things we can determine from a basic pulsar timing solution. If you'd like to proceed with the PTA-specific tutorials, you can skip this section with no consequneces. However, if you'd like, we can callculate a few properties of our simulated pulsar.

First, to make things easy, let's convert from $f$ and $\dot{f}$ to $P$ and $\dot{P}$. We have by definition that the spin period and frequency are related by $P = 1/f$. Using some calculus, then we have that

$$\frac{dP}{df} = -\frac{1}{f^2}$$

and therefore

$$\dot{P} = -\frac{\dot{f}}{f^2}$$

.a

Use the value of $P$ given above (or re-calculate it) and calculate $\dot{P}$. Note that the units of $\dot{P}$ are seconds per seconds.

Let's calculate the amount of energy lost by the system, given by (see Lorimer & Kramer 2005)

$$\dot{E} = 4 \pi^2 I \dot{P} P^{-3} \simeq 3.95 \times 10^{31}~\mathrm{erg~s^{-1}}\left(\frac{\dot{P}}{10^{-15}}\right)\left(\frac{P}{\mathrm{s}}\right)^{-3}$$

How does this compare to the luminosity of the Sun ($4.0 \times 10^{33}~\mathrm{ergs~s^{-1}}$)?

Now let's calculate the surface magnetic field of the pulsar, given by

$$ B_s = 3.2 \times 10^{19}~\mathrm{G} \sqrt{P \dot{P}} \simeq 10^{12}~\mathrm{G}~\left(\frac{\dot{P}}{10^{-15}}\right)^{1/2} \left(\frac{P}{\mathrm{s}}\right)^{1/2} $$

For reference, Earth's magnetic field is usually around $\lesssim 1~\mathrm{G}$!