1 of 42

i do not know what it is about you that closes

and opens; only something in me understands

the voice of your eyes is deeper than all roses

-- E. E. Cummings

Challenges of pulsar timing array (PTA) data analysis

Golam Shaifullah

2 of 42

It’s not a bird, it’s not a plane, definitely not LGM

Pulsars are giant flywheels in space, their compact masses give rise to incredibly stable rotation.

On each rotation, the pulsar beam produces a ‘pulse’ at Earth, and the photons in that pulse can be assigned a time-of-arrival (TOA).

2

Discovered by

Prof. Bell-Burnell

/ 26

3 of 42

Models, models, models

TOAs can be predicted using a model with the following (sets of) parameters:

astrometric,
pulsar rotation and
binary (when applicable).

Apart from these pulsar emission is affected by:

Dispersive delays due to the intervening ionised plasma
Red noise (low frequency) processes

3

/ 26

4 of 42

Pulsar timing

However, once we have estimates of those parameters, we can predict very precisely when the next pulse will arrive. Or the one after 20 million rotations.

When pulses are averaged this precision quickly tends to tens of microseconds to hundreds of nanoseconds.

4

Time tagged to a precision of picoseconds

Observed pulse train

Template & polynomial model

Timing Residuals

/ 26

5 of 42

Millisecond pulsars as stable clocks

5

See Shannon et al (2016), Lam et al (2018) & others

Adapted from Hartnett & Luiten, 2011

psrqpy, (Pitkin, 2011)

/ 26

6 of 42

A quick word from our (astrophysical) sources

6

Burke-Spolaor A&A Rev. (2019)

/ 26

7 of 42

All of the light we cannot see

7

Izquierdo-Villalba et al (2024),

Curylo et al (2023), Bromm & Loeb 2003, Cole et al 2000, Benson (2012)

Manzini, MSc thesis, 2023

/ 26

8 of 42

The true picture & a wider landscape

8

Ellis et al (2023; 2308.08546)

/ 26

9 of 42

Pulsar timing arrays

GWs are expected to induce timing residuals on the order of a few tens of nanoseconds.

TOA stability scales with number of rotations averaged - use millisecond pulsars (MSPs)!

Single pulsars are ‘jittery’ and affected by noise, use an array of MSPs

9

Ferranti, MSc thesis, 2023

/ 26

10 of 42

FreqBayes^TM pulsar timing:

Observe a pulsar
De-disperse
Stack
Average
Make a template
Cross-correlate
Line up your TOAs
Repeat for another 20 - 100 sources
Sprinkle post-docs for flavour
Bake for ~30 years, turning it over once or twice a decade.

10

Figure from Verbiest & Shaifullah, 2018, CQG

/ 26

PTA data creation starts with observations of pulsars using radio telescopes. Each observation lasts for a total integration time of between 30 to 60 minutes for the EPTA. The observation or integration lengths are typically chosen to meet a predefined signal to noise ratio for the pulsar profile. Each integration is recorded as a series of smaller sub-integrations, which are typically tens of seconds long. Each sub-integration is further made up of individual pulsar rotations stacked or ‘folded’ to increase the S/N of the pulse-profile (as sketched in the box marked ‘Folding’).

When averaging across wide bandwidths, it is also necessary to correct for the interstellar dispersion of the radio pulses, which I have sketched here before the folding step, but in principle can be performed after the folding as well. Given the ever increasing bandwidths of pulsar receivers, this is an extremely important step to recover the profile of the pulsar and has direct consequences on the achievable TOA precision -- see e.g., coherent and incoherent dedispersion in the Handbook of pulsar astronomy (Lorimer & Kramer, 2005).

Once the observations have been averaged in time (for instance the final sketch on the top right); we can further integrate along the frequency axis to produce a time and frequency averaged profile. Typically, we retain only the total intensity (Stokes I) component at this stage, but the total observed band can be split into multiple sub-bands to produce multiple profiles per observation, referred to the individual sub-bands.

Once a number of observations have become available we can also construct a template for the most regular features of the pulse-profile. Although the pulse profile is the detector response as a function of the pulsar’s rotational phase, and small differences are expected for each rotation, these differences get averaged out over multiple pulsar rotations leading to a very characteristic average profile for each pulsar. See e.g., Shannon & Cordes (2011), Oslowski et al (), etc for a discussion of the time-scales at which the per rotation variations (or single pulse variations) occur and the limits these place on the achievable TOA precision or the timing RMS. However it is possible to show that the average profile is in fact highly repeatable, e.g., as shown in Liu et al (2011).

11 of 42

What is the signal PTAs are looking for?

11

Hellings & Downs, 1983

Also see Romano & Allen, arxiv:2308.05847

/ 26

12 of 42

EPTA data combination 2.0

5+1 telescopes + LOFAR (next gen).
Augmented with the EPTA-DR1 (Desvignes et al, 2016).
New data spanning ~10 years (a total of 24.5)
Double the observing bandwidth
Coherent dedispersion
Two to ten times greater observing cadence
Significantly boost timing sensitivity
Combine only 25 out of 42 sources

12

/ 26

13 of 42

EPTA + InPTA

New data spanning ~3 years
Overlap of 10 sources
Significantly improves DM sensitivity
Simultaneous low and high frequency!

13

/ 26

14 of 42

14

EPTA DR2 - Paper I A&A , 2023, doi: 10.1051/0004-6361/202346841

Zenodo: https://doi.org/10.5281/zenodo.8164424

Gitlab: https://epta.pages.in2p3.fr/epta-dr2

/ 48

15 of 42

The detection statistic and search algorithm

15

pulsar index

We assume that noise is Gaussian: the likelihood function (likelihood of the signal with given parameters) is

P(δ⃗t, 𝜽⃗) = {1/√(2π)ⁿdet(C) } exp( - ½ (δ⃗t - s⃗)^T C^-1(δ⃗t - s⃗))

δt - concatenated residuals from all pulsars in the array: total size n
s - is a model of deterministic signals (e.g. - GW signals from individually resolvable SMBHBs)
C is the noise variance-covariance matrix (size n ⨯ n );

C_𝜶i,𝜷j = C^WN δ_𝜶𝜷δ_ij + C^RN _ij δ_𝜶𝜷 + C^DM _ij δ_𝜶𝜷 + C^GW_ij δ_𝜶𝜷 + …

white

noise

red (spin)

noise

dispersion

noise

stochastic

GW

noise

toa index

/ 26

16 of 42

The design

Rows correspond to each observation epoch.
Columns correspond to the parameters of the pulsar timing model (e.g., pulse period, spin-down rate, astrometric parameters, binary parameters, etc.).

The matrix transforms the model parameters into predicted changes in pulse arrival times. Mathematically, the residuals r can be expressed as:

r=t−M⋅p

r is the vector of timing residuals,
t is the vector of observed TOAs,
M is the design matrix,
p is the vector of pulsar timing model parameters.

16

1	t1	f_π(t1)	f_μα(t1)	f_μδ(t1)		f_Pb(t1)
1	t2	f_π(t2)	f_μα(t2)	f_μδ(t2)		f_Pb(t2)
1	t3	f_π(t3)	f_μα(t3)	f_μδ(t3)		f_Pb(t3)
...	...	...	...	...	...	...
1	ti	f_π(ti)	f_μα(ti)	f_μδ(ti)		f_Pb(ti)
...	...	...	...	...	...	...
...	...	...	...	...	...	...
1	tn	f_π(tn)	f_μα(tn)	f_μδ(tn)		f_Pb(tn)

M =

/ 26

17 of 42

PTA noise sources

17

Figure adapted from Verbiest & Shaifullah, 2018, CQG

/ 26

PTA data creation starts with observations of pulsars using radio telescopes. Each observation lasts for a total integration time of between 30 to 60 minutes for the EPTA. The observation or integration lengths are typically chosen to meet a predefined signal to noise ratio for the pulsar profile. Each integration is recorded as a series of smaller sub-integrations, which are typically tens of seconds long. Each sub-integration is further made up of individual pulsar rotations stacked or ‘folded’ to increase the S/N of the pulse-profile (as sketched in the box marked ‘Folding’).

When averaging across wide bandwidths, it is also necessary to correct for the interstellar dispersion of the radio pulses, which I have sketched here before the folding step, but in principle can be performed after the folding as well. Given the ever increasing bandwidths of pulsar receivers, this is an extremely important step to recover the profile of the pulsar and has direct consequences on the achievable TOA precision -- see e.g., coherent and incoherent dedispersion in the Handbook of pulsar astronomy (Lorimer & Kramer, 2005).

Once the observations have been averaged in time (for instance the final sketch on the top right); we can further integrate along the frequency axis to produce a time and frequency averaged profile. Typically, we retain only the total intensity (Stokes I) component at this stage, but the total observed band can be split into multiple sub-bands to produce multiple profiles per observation, referred to the individual sub-bands.

Once a number of observations have become available we can also construct a template for the most regular features of the pulse-profile. Although the pulse profile is the detector response as a function of the pulsar’s rotational phase, and small differences are expected for each rotation, these differences get averaged out over multiple pulsar rotations leading to a very characteristic average profile for each pulsar. See e.g., Shannon & Cordes (2011), Oslowski et al (), etc for a discussion of the time-scales at which the per rotation variations (or single pulse variations) occur and the limits these place on the achievable TOA precision or the timing RMS. However it is possible to show that the average profile is in fact highly repeatable, e.g., as shown in Liu et al (2011).

18 of 42

Noise models & their validity

18

EPTA DR2 - Paper II A&A , 2023, doi: 10.1051/0004-6361/202346842

/ 26

19 of 42

19

https://www.epta.eu.org/epta-dr2.html

/ 48

20 of 42

PTAs inching up to the GWB

On June 29, 2023 4 PTAs announced evidence for an HD correlated process in their data.

The significance ranges from ~2 to 4.6σ; below the 5σ detection threshold.

Further this amplitude is loud (~2-3 x 10^-15) and the spectrum is flat (~3).

20

/ 26

21 of 42

Do the PTAs agree?

21

The IPTA collaboration, 2023

/ 26

22 of 42

Extended EPTA DR2 to 60 psrs

23 of 42

Data

24 of 42

EPTA - latest results & the near future

Future:

New data up to 2024 January
Expanding up to 60 best pulsars (~100 timed)
Will commit data to IPTA for as many as feasible.
Highest cadence (~3-7 days)
Longest PTA dataset (~27 years)
With InPTA, sensitive from 350 MHz up to 5 GHz
Adding LOFAR and NENUFAR to go down to 20MHz.
7 operating telescopes - 25 MHz to 100 GHz

24

Sensitivity forecast EPTA 2024 (60 PSRs; 27 years)

/ 26

25 of 42

Noise analysis - focus

26 of 42

26

MPTA

/ 48

27 of 42

The IPTA DR3

27

Add data from 5+ PTAs:

EPTA (+ LOFAR, NENUFAR)
NANOGrav (+ CHIME)
PPTA
InPTA
MPTA (MeerKAT)

2+ independent data combination pipelines
121 pulsars, down to <100 ns for a few pulsars
Greater sky coverage!
More pulsar pairs for angular correlation searches.
Lots of TOAs
Loads of compute
Currently working on “Early Data Release” (eDR3), which includes the 20 best/longest-timed pulsars
20 pulsars have been combined, first noise runs too!

PSR J1909-3744 combination with PINT pipeline

PSR J1909-3744 combination with Tempo2 pipeline

Fig: D. Good

Fig: G. Shaifullah

/ 26

28 of 42

IPTA DR3 dimensions

28

PTA	Dataset	PSRs	Tspan (years)	f_GW,low (nHz)	f_radio (MHz)
EPTA	DR2 / DR3	25 / +35	24.5	1.29	283 - 5107
EPTA	LOFAR + NENUFAR	17	9.6	-	30 - 190
NANOGrav	15-yr	68	15.9	1.99	302 - 3988
NANOGrav	CHIME	11	2.5		400 - 800
PPTA	DR3	24	18.1	1.75	704 - 4032
InPTA	DR1	15	3.5	9.05	300 - 1460
MeerKAT	DR2	88	4.5	7.04	856 - 1412
IPTA	DR3	121	~25/40	1.29/0.79	30 - 5107

/ 26

29 of 42

IPTA DR3 dimensions

29

In total 121 pulsars in full DR3;

The biggest / most sensitive PTA dataset ever made !!

The EPTA + InPTA result - a loud background
SMBHB generated backgrounds
Comparisons with Semi-Analytical Models
Stellar hardening?
Biased by cosmic variance?
Cosmic Strings
Curvature perturbations
Challenging the ultralight dark matter paradigm

/ 26

36 of 42

The astrophysical implications

36

… but a biased consideration of the uncertainties of the Hellings & Downs curve might.

The EPTA + InPTA result - a loud background
SMBHB generated backgrounds
Comparisons with Semi-Analytical Models
Stellar hardening?
Biased by cosmic variance?
Cosmic Strings
Curvature perturbations
Challenging the ultralight dark matter paradigm

/ 26

37 of 42

The astrophysical implications

37

The EPTA + InPTA result - a loud background
SMBHB generated backgrounds
Comparisons with Semi-Analytical Models
Stellar hardening?
Biased by cosmic variance?
Cosmic Strings
Curvature perturbations
Challenging the ultralight dark matter paradigm

/ 26

38 of 42

The astrophysical implications

38

The EPTA + InPTA result - a loud background
SMBHB generated backgrounds
Comparisons with Semi-Analytical Models
Stellar hardening?
Biased by cosmic variance?
Cosmic Strings
Curvature perturbations
Challenging the ultralight dark matter paradigm

/ 26

39 of 42

The astrophysical implications

39

The EPTA + InPTA result - a loud background
SMBHB generated backgrounds
Comparisons with Semi-Analytical Models
Stellar hardening?
Biased by cosmic variance?
Cosmic Strings
Curvature perturbations
Challenging the ultralight dark matter paradigm

/ 26

40 of 42

Nihil Ultra?

40

J2145-0750

J1713+0747

J0218+4232

The one-dimensional marginalized Bayesian posteriors for the amplitude of a low-frequency spectral process with a fixed spectral index of γ = 13/3.

Wang, GMS et al, A&A (2022)

/ 26

41 of 42

Outliers (Fumagalli, GMS + submitted)

41

/ 26

42 of 42

Testing biases in the MCMC algorithms

42

Red Noise models

DM Noise models

Samajdar, GMS et al (2022), arXiv:2205.04332

/ 26