The MeerKAT Pulsar Timing Array: Maps of the gravitational wave sky with the 4.5-yr data release

ABSTRACT

In an accompanying publication, the MeerKAT Pulsar Timing Array (MPTA) Collaboration reports tentative evidence for the presence of a stochastic gravitational wave background, following observations of similar signals from the European and Indian Pulsar Timing Arrays, the North American Nanohertz Observatory for Gravitational Waves, the Parkes Pulsar Timing Array, and the Chinese Pulsar Timing Array. If such a gravitational wave background signal originates from a population of inspiraling supermassive black hole binaries, the signal may be anisotropically distributed in the sky. In this paper, we evaluate the anisotropy of the MPTA signal using a spherical harmonic decomposition. We discuss complications arising from the covariance between pulsar pairs and the regularization of the Fisher matrix. Applying our method to the |$4.5 \hbox{-}\text{yr}$| data set, we obtain two forms of sky maps for the three most sensitive MPTA frequency bins between |$7 \ {\rm and} \ 21 \, {\rm nHz}$|⁠. Our ‘clean maps’ estimate the distribution of gravitational wave strain power with minimal assumptions. Our radiometer maps answer the question: Is there a statistically significant point source? We find a noteworthy hotspot in the |$7 \, \mathrm{nHz}$| clean map with a p-factor of |$p=0.015$| (not including trial factors). Future observations are required to determine if this hotspot is of astrophysical origin.

gravitational waves, methods: data analysis, methods: statistical, stars: black holes, pulsars: general

1 INTRODUCTION

Acting as cosmic clocks, millisecond pulsars emit extraordinarily regular pulses at radio frequencies that are observed with ground-based radio telescopes. The stability of millisecond pulsar rotation allows us to build accurate timing models that account for the physics both intrinsic and extrinsic to the pulsar, with the most accurate models able to predict pulse times of arrival (ToAs) to a precision of tens of nanoseconds (e.g. Agazie et al. 2023b; Reardon et al. 2024; Wang et al. 2024). Distortions in space–time caused by gravitational waves induce a temporally correlated signal in the arrival times and hence residuals of each pulsar. In a collection of pulsars, a so-called pulsar timing array (PTA), this gravitational wave signal appears as a correlation between the residuals of different pulsars depending on the separation angle between these pulsars (Hellings & Downs 1983).

PTAs are sensitive to nanohertz gravitational waves with periods of years to decades. In this frequency range, gravitational waves are likely to create a stochastic gravitational wave background from the superposition of incoherent gravitational waves radiated from many sources (Allen & Romano 1999; Rosado, Sesana & Gair 2015). The signal is expected to first appear as a time-correlated red noise process (Siemens et al. 2013), characterized by a power spectral density with a steep power-law index. However, this red process may be mimicked by pulsar noise processes (Hazboun et al. 2020; Goncharov et al. 2021, 2022; Zic et al. 2022), and so angular correlations are crucial for establishing a confident detection.

An isotropic gravitational wave background induces an angular correlation function known as the Hellings–Downs curve, which depends only on the angular separation of pairs of pulsars (Hellings & Downs 1983). The unique (nearly quadrupolar) correlations of the Hellings–Downs curve are distinct from other correlations that may be induced by systematic errors such as clock errors (which induce monopole moments) and Solar system ephemeris errors (which induce dipole moments; Tiburzi et al. 2016). A statistically significant Hellings–Downs correlation is the evidence required for gravitational wave background detection.

In 2023 June, there was a coordinated release of pulsar timing papers by members of the International Pulsar Timing Array (IPTA), with each presenting evidence for a gravitational wave background signal. The North American Nanohertz Observatory for Gravitational Waves (NANOGrav) reported evidence for a gravitational wave background with a significance of |${\sim} 3 \sigma \!-\! 4 \sigma$| (Bayesian versus frequentist methods provide slightly different levels of evidence; see Agazie et al. 2023a), while the joint publication from the European and Indian Pulsar Timing Arrays (EPTA and InPTA, respectively) reported |$3 \sigma \!-\! 4 \sigma$| (using a subset of the full data; see EPTA Collaboration & InPTA Collaboration 2023). The Parkes Pulsar Timing Array (PPTA) simultaneously found |${\sim} 2 \sigma$| evidence for a gravitational wave background signal (Reardon et al. 2023). In addition to these IPTA publications, the Chinese Pulsar Timing Array (CPTA) reported a |$4.6 \sigma$| evidence for the presence of a Hellings–Downs correlation in the data at the frequency slice of |$14 \, \mathrm{nHz}$| (Xu et al. 2023), although due to the frequency-dependent analysis, the CPTA result is not directly comparable to the results from the other PTAs.

Once a gravitational wave background has been established, one of the next key steps is to characterize how the signal is distributed in the sky. In this paper, we assess the anisotropy of the nanohertz gravitational wave background in the hope of determining its source.

The most likely source of the gravitational wave background is a population of inspiraling supermassive black hole binaries (Rajagopal & Romani 1995). Supermassive black holes are thought to be found in the centre of most massive galaxies (Kormendy & Richstone 1995), and mergers between these galaxies gravitationally bind their central black holes, leading to the formation of supermassive black hole binaries (Begelman, Blandford & Rees 1980). Circular, gravitationally driven supermassive binary black holes are expected to stall in a phenomenon known as the final parsec problem (Milosavljevic & Merritt 2003), though various solutions have been proposed including multiple-body interactions (Ryu et al. 2018; Bonetti et al. 2019), hardening from interactions with stars (Gualandris et al. 2017), and galaxy rotation (Holley-Bockelmann & Khan 2015).

A gravitational wave background arising from supermassive black hole binaries is likely to exhibit anisotropy from the finite number of merging binaries in the PTA observing band. While Mingarelli et al. (2017) estimated that the level of anisotropy due to undetected continuous wave sources is about 20 per cent, other studies (e.g. Sesana et al. 2004; Simon 2023) argue that there is a large uncertainty involved, due to our relatively poor knowledge of the supermassive binary black hole population.

Alternatively, a gravitational wave background may arise from processes in the early Universe (Maggiore 2000). Primordial gravitational waves from quantum fluctuations redshifted by inflation (Guzzetti et al. 2016; Lasky et al. 2016; Domènech 2021; Yuan & Huang 2021), phase transitions (Caprini et al. 2020; Hindmarsh et al. 2021), or cosmic strings (Hindmarsh & Kibble 1995; Saikawa 2017) are likely to give rise to an isotropic signal. Understanding the degree of anisotropy in the gravitational wave background therefore provides clues about its source. Recent analyses from NANOGrav (Agazie et al. 2023c) and the EPTA (Taylor et al. 2015) have yielded no evidence for anisotropy in the nanohertz gravitational wave background.

The remainder of this paper is organized as follows. We start by giving a brief overview of the MPTA and the data set used in this analysis in Section 2. This is then followed by an in-depth description of the methods we use for evaluating the anisotropy of the nanohertz gravitational wave background in Section 3. Our method follows techniques developed in the context of the Laser Interferometer Gravitational Wave Observatory (LIGO) and subsequently PTAs. However, we extend this framework to take into account the contribution of noise from the gravitational wave background itself. In Section 4, we present the results obtained from our analysis of the MeerKAT Pulsar Timing Array (MPTA) 4.5-yr data set. We conclude in Section 5 with a discussion of these findings and the future outlook. Lastly, the appendices of this paper contain a series of auxiliary calculations and tests documenting the performance of our data analysis pipeline.

2 THE MeerKAT PULSAR TIMING ARRAY

The MPTA uses the MeerKAT radio telescope, a 64-antenna array, located in the Northern Cape, South Africa (Jonas & MeerKAT Team 2018). Together these antennas provide a gain of |$2.8 \, \mathrm{K}\text{Jy}^{-1}$| making MeerKAT the most sensitive radio telescope in the Southern hemisphere. The MPTA is one subproject of the MeerTime large survey project (Bailes et al. 2020). For pulsar timing work the beam-formed data are coherently dedispersed (to remove the dispersive effect of the ionized interstellar medium) and folded at the topocentric period of the pulsar. They are processed to remove interference signals and calibrate polarization resulting in Stokes I pulse profiles with 32 frequency channels. In most cases the observations were fully time averaged; however, some long observations (usually observed as part of other subprojects) were subdivided to better match the integration time of standard MPTA observations (Miles et al. 2024a).

A frequency-resolved template (describing the pulse profile) is produced for each pulsar using the PulsePortraiture software (Pennucci 2019). The ToAs (of the pulses closest to the middle of the observations) are determined using the Fourier domain Monte Carlo method in the psrchive software suite (Hotan, van Straten & Manchester 2004) for each channel. The data used in this work come from a 4.5-yr period from 2019 February to 2024 August. The MPTA observes 83 pulsars with an approximately 14-d cadence. Observations were done using the L-band receiver covering a bandwidth of 856 MHz centred at |$1284 \, \mathrm{MHz}$| (see Miles et al. 2024a for more detail).

The frequency-resolved information is used to constrain models of chromatic effects, such as variations in the ionized interstellar medium, and the effect of the solar wind. Achromatic noise is modelled for both red and white noise processes. Details of this noise modelling can be found in Miles et al. (2024a). The MPTA reports a |$3.4 \sigma$| evidence for a stochastic isotropic gravitational wave background in the |${\sim} 7 \!-\!21 \, {\rm nHz}$| frequency range (Miles et al. 2024b). This analysis uses the results from the ALT analysis in Miles et al. (2024b), a conservative set of noise models where additional achromatic noise processes have been added (for all pulsars barring PSR J2129−5721). The excellent sensitivity of MeerKAT enables us to study the gravitational wave background at high frequencies and high observational cadence. High-frequency data sets are particularly useful in searches for anisotropy because the number density of binary black holes is expected to fall sharply with frequency |$p(f) \propto f^{-11/3}$| (assuming gravitationally driven, circular binaries; Peters 1964; Rajagopal & Romani 1995; Smith & Thrane 2018; Gardiner et al. 2024). Thus, at higher frequencies, we expect the background to be created by a relatively smaller number of binaries.

3 METHODS

3.1 Overview

Gravitational wave cartography was pioneered in the context of ground-based detector community (e.g. Allen & Ottewill 1997; Ballmer 2006; Mitra et al. 2008; Thrane et al. 2009; Abadie et al. 2011). In the PTA community, Mingarelli et al. (2013) derived a generalized overlap reduction function for PTA analyses, leading to the first Bayesian PTA analysis strategy towards estimating the anisotropy of a gravitational wave background by Taylor & Gair (2013), which was later generalized by Gair et al. (2014). A frequentist analysis pipeline using optimal statistics (e.g. Vigeland et al. 2018) and the likelihood maximization developed by Thrane et al. (2009) and Romano & Cornish (2017) was recently put forward by Pol, Taylor & Romano (2022).

While there are differences, these methods employ many common features: they define a basis consisting of either pixels or spherical harmonic functions (Cornish & van Haasteren 2014; Ali-Haïmoud, Smith & Mingarelli 2020, 2021; Taylor, van Haasteren & Sesana 2020; Banagiri et al. 2021; Pol et al. 2022; Agazie et al. 2023c) and estimate the gravitational wave power (proportional to strain squared) associated with each basis element by maximizing the likelihood. These are the most common choices. Other bases such as PTA eigenmaps (Cornish & van Haasteren 2014; Ali-Haïmoud et al. 2020, 2021) are also considered. Here, we follow the method outlined in Thrane et al. (2009) and Abadie et al. (2011) – and adapted for PTAs in Pol et al. (2022) and references therein – which employs a spherical harmonic basis. In this framework, the reconstructed sky map is allowed to have regions with negative gravitational wave power.

Since the stochastic background can only induce positive power, these patches of negative power represent noise fluctuations. While some authors chose to enforce positive power [e.g. using a ‘square root basis’ (Payne et al. 2020; Taylor et al. 2020; Pol et al. 2022), or through the use of priors (Taylor & Gair 2013; Gair et al. 2014)], we prefer to allow negative power as a way of visualizing the noise fluctuations. Indeed, we use the behaviour of the sky map noise fluctuations as a diagnostic to check that our pipeline produces reasonable results.

We highlight two dedicated PTA searches for anisotropy. The first one was carried out by Taylor et al. (2015) as part of the EPTA Collaboration (Kramer & Champion 2013) using the timing data from the six best-quality pulsars. The second and most recent anisotropy analysis by Agazie et al. (2023c) was presented in the paper’s series following the 15-yr NANOGrav data release. Neither found evidence for anisotropy.

In the following subsections, we introduce the formalism of our approach. Speaking broadly to give a high-level overview, we follow the spherical harmonic formalism from Mitra et al. (2008) to derive ‘clean maps’, through the regularization scheme described in Thrane et al. (2009) using the generalized overlap reduction functions from Mingarelli et al. (2013). We make some modifications to relax the assumption of the weak-signal limit.

3.2 The overlap reduction function

A PTA consisting of |$N_\mathrm{psr}$| pulsars allows for the calculation of |$N_\text{pairs} = (N_\mathrm{psr}^2-N_\mathrm{psr})/2$| correlations between distinct pulsar pairs. These correlations encode the distribution of gravitational wave power on the sky via the so-called overlap reduction function (Christensen 1992; Allen & Ottewill 1997; Finn, Larson & Romano 2009; Mingarelli et al. 2013):

$$\begin{eqnarray} \Gamma _{ab} \propto \int _{S^2} \! \mathrm{d}\hat{\boldsymbol{\Omega }} \,\, \mathcal {P}(\hat{\boldsymbol{\Omega }}) \left[\mathcal {F}_a^+(\hat{\boldsymbol{\Omega }})\mathcal {F}_b^+(\hat{\boldsymbol{\Omega }}) + \mathcal {F}_a^\times (\hat{\boldsymbol{\Omega }})\mathcal {F}_b^\times (\hat{\boldsymbol{\Omega }})\right] . \end{eqnarray}$$

(1)

Here, |$\mathcal {P}(\hat {\Omega })$| is the probability density function for gravitational wave power at different locations on the sky |$\hat{\Omega }$| as seen from the Solar system barycentre – this is the distribution that we seek to measure. The |$a,b$| indices label different pulsars. The quantity

$$\begin{eqnarray} \mathcal {F}_a^A(\hat{\boldsymbol{p}}, \hat{\boldsymbol {\Omega }}) = \frac{1}{2}\frac{\hat {p}_a^i\hat {p}_a^j}{1-\hat {\boldsymbol{\Omega }}\cdot \hat{\boldsymbol{p}}_a}\, e_{ij}^A(\hat {\boldsymbol{\Omega }}) = \frac{1}{2}\frac{\hat {p}_a^i\hat {p}_a^j}{1+\hat{\boldsymbol{k}}\cdot \hat{\boldsymbol{p}}_a}\, e_{ij}^A(\hat{\boldsymbol{\Omega }}) \end{eqnarray}$$

(2)

is the antenna factor for a gravitational wave propagating along the vector |$\hat{\boldsymbol{k}}$|⁠, i.e. the source is located at position |$\hat{\boldsymbol{\Omega} }$| with |$\hat{\boldsymbol{k}}$| = |$- \hat{\boldsymbol{\Omega} }$|⁠, with polarization state A as measured by a pulsar located in the direction |$\hat{\boldsymbol{p}}_a$|⁠. The term |$e^A_{ij}(\hat{\boldsymbol{\Omega }})$| is the polarization tensor for a gravitational wave with polarization A; |$i,j$| are indices for spatial coordinates. In our analysis, we adopt the definition of the polarization basis from Taylor et al. (2020).

3.3 Spherical harmonic decomposition

We expand the distribution of gravitational wave power in the complex spherical harmonics basis:

$$\begin{eqnarray} \mathcal {P}(\hat {\boldsymbol{\Omega }}) = \sum _{l=0}^{\ell _\text{max}}\sum _{m=-\ell }^{\ell } P_{\ell m} Y_{\ell m}(\hat {\boldsymbol{\Omega} }), \end{eqnarray}$$

(3)

where |$Y_{\ell m}(\hat {\boldsymbol{\Omega }})$| are the complex-valued spherical harmonics and |$P_{\ell m}$| the respective coefficients (Cornish 2001; Kudoh & Taruya 2005; Cornish & van Haasteren 2014). The sum truncates at |$\ell _\text{max}$|⁠, which we choose in order to match with the intrinsic resolution of our PTA. Since we have 83 pulsars in our array, a single-frequency-bin map is characterized by at most |${\lesssim} 83$| independent parameters. As we discuss below, the number is actually much smaller in practice owing to degeneracies between measurements. Thus we choose |$\ell _\text{max}=8$|⁠, which yields 81 degrees of freedom. It follows that our angular resolution is |$180^\circ /\ell _\text{max} \approx 23^\circ$|⁠.

3.4 Response matrix

Some additional bookkeeping is necessary before we proceed. We introduce the Greek indices |$\alpha$| and |$\beta$| to denote pulsar pairs |$(ab)$| so, for example, we can write |$\Gamma _\alpha$| as shorthand for |$\Gamma _{ab}$|⁠. Meanwhile, we use the indices |$\mu$| and |$\nu$| to denote the distinct spherical harmonics, i.e., for example, we write the spherical harmonic basis functions as |$Y_\mu$|⁠. A summary of indices is provided in Table 1.

Table 1.

A summary of different indices used in this paper.

Index	Represents	Maximum range
|$a,b$|	Pulsar	83
|$\alpha ,\beta$|	Pairs	3403
|$\mu ,\nu$|	Spherical harmonic	81

With this notation, we can relate the overlap reduction function for some pair |$\alpha$| to the spherical harmonic coefficients |$P_\mu$| with a matrix equation

$$\begin{eqnarray} \Gamma _\alpha = R_{\alpha \mu } P_\mu , \end{eqnarray}$$

(4)

where

$$\begin{eqnarray} R_{\alpha \mu } = \int {\rm d}\hat{ \boldsymbol{\Omega}} \left[\mathcal {F}_a^+(\hat{\boldsymbol{\Omega }})\mathcal {F}_b^+(\hat{\boldsymbol{\Omega }}) + \mathcal {F}_a^\times (\hat{\boldsymbol{\Omega }})\mathcal {F}_b^\times (\hat{\boldsymbol{\Omega }})\right] Y_\mu (\hat{\boldsymbol{\Omega }}) \end{eqnarray}$$

(5)

is referred to as the ‘response matrix’ and we adopt the Einstein convention, so that repeated indices imply summation.¹

3.5 Covariance matrices

PTA data consist of |$N_\mathrm{psr}$| vectors of pulse ToAs |$\boldsymbol{\delta t}$|⁠. By taking the expectation values of the product of different residuals, we construct covariance matrices:

$$\begin{eqnarray} \bf {C}_a = \langle \boldsymbol{\delta t}_a\boldsymbol{\delta t}_a^{\rm T}\rangle , \end{eqnarray}$$

(6)

$$\begin{eqnarray} \bf {S}_{ab} &= \langle \boldsymbol{\delta t}_a\boldsymbol{\delta t}_b^{\rm T}\rangle _{a\ne b} . \end{eqnarray}$$

(7)

Here, |$\bf{C}_a$| is the autocorrelation variance and |$\bf{S}_{ab}$| is the cross-correlation variance matrix.² The autocorrelation matrix is described by the sum of all individual noise processes determined in the single-pulsar noise analysis.

The cross-power matrix, on the other hand, depends only on the gravitational wave signal that is common to all the pulsars in the array:

$$\begin{eqnarray} \bf{S}_{ab} = \boldsymbol{F}_a \Gamma _{ab} P_\mathrm{GWB}(f) \boldsymbol{F}_b^\dagger . \end{eqnarray}$$

(8)

Here, |$P_\mathrm{GWB}(f)$| is the power spectral density of the gravitational wave signal,

$$\begin{eqnarray} P_\mathrm{GWB}(f) = \frac{A_\mathrm{GWB}^2}{12\pi ^2}\left(\frac{f}{f_\mathrm{yr}}\right)^{-\gamma _\mathrm{GWB}}, \end{eqnarray}$$

(9)

which we parametrize in terms of an amplitude |$A_\mathrm{GWB}$| and a spectral index |$\gamma _\mathrm{GWB}$|⁠. For gravitationally driven binary inspirals, a value of |$\gamma _\mathrm{GWB}=13/3$| is expected. The Fourier basis vectors |$\boldsymbol{F}_{a,b}$| convert the power spectral density from the frequency domain to the time domain.

3.6 Likelihood

Following Demorest et al. (2013), Chamberlin et al. (2015), and Pol et al. (2022), we define pulsar cross-correlations |$\rho _\alpha$| and their respective uncertainties |$\sigma _\alpha$|⁠:

$$\begin{eqnarray} \rho _\alpha \equiv \rho _{ab} = \frac{\boldsymbol{\delta t}^\dagger _a \bf {C}_a^{-1} {\hat {\bf{S}}}_{ab}\bf {C}_b^{-1}\boldsymbol{\delta t}_b}{\mathrm{tr}\left[\bf {C}_a^{-1} {\hat {\bf{S}}}_{ab}\bf {C}_b^{-1} {\hat {\bf{S}}}_{ba}\right]}, \end{eqnarray}$$

(10)

$$\begin{eqnarray} \sigma _\alpha \equiv \sigma _{ab} &= \left(\mathrm{tr}\left[\bf {C}_a^{-1}{\hat {\bf{S}}}_{ab}\bf {C}_b^{-1} {\hat {\bf{S}}}_{ba}\right]\right)^{-1/2}, \end{eqnarray}$$

(11)

where

$$\begin{eqnarray} {\hat {\bf{S}}}_{ab} = \bf {S}_{ab}/A_\mathrm{GW}^2 \Gamma _{ab} . \end{eqnarray}$$

(12)

The traces in equations (10) and (11) are taken over the ToA indices. The cross-correlations are a convenient reduced data product, which enables us to work with a simple likelihood function.

Assuming stationary Gaussian noise (and a stationary Gaussian stochastic background), the likelihood for the cross-correlations can be written as (Mitra et al. 2008; Thrane et al. 2009; Pol et al. 2022)

$$\begin{eqnarray} {\cal L}(\rho |\mathcal {P}) = \frac{\exp \left[-\frac{1}{2}\left(\rho _\alpha -R_{\alpha \mu }\mathcal {P}_\mu \right)^\dagger \Sigma _{\alpha \beta }^{-1}\left(\rho _\beta -R_{\beta \nu } \mathcal {P}_\nu \right)\right]}{\sqrt{\det (2\pi \boldsymbol {\Sigma })}} . \end{eqnarray}$$

(13)

Here, |$\boldsymbol {\Sigma }$| is the second moment³ of the cross-correlations (Allen & Romano 2023):

$$\begin{eqnarray} \Sigma _{\alpha \beta } = \mathrm{diag}(\sigma _\alpha ^2)\, \delta _{\alpha \beta } + \varsigma _{\alpha \beta }(\boldsymbol{\cal P}, A_\text{GWB}) . \end{eqnarray}$$

(14)

The diagonal of |$\boldsymbol{\Sigma}$| is just the variance associated with each cross-correlation. However, when a gravitational wave signal is present the term |$\varsigma _{\alpha \beta }$| (which includes off-diagonal components) becomes important (Allen & Romano 2023). The rather complicated expressions for |$\boldsymbol{\varsigma}$| are provided in Appendix A2.

The likelihood is maximized with the following solution (Mitra et al. 2008; Thrane et al. 2009; Pol et al. 2022):

$$\begin{eqnarray} \mathcal {P}_\mu ^{\prime } = M_{\mu \nu }^{-1} X_\nu , \end{eqnarray}$$

(15)

where

$$\begin{eqnarray} X_\nu &= R_{\mu \alpha }^\dagger \Sigma _{\alpha \beta }^{-1} \, \rho _\beta \end{eqnarray}$$

(16)

is known as the dirty map, and

$$\begin{eqnarray} M_{\mu \nu } &= R_{\mu \alpha }^\dagger \Sigma _{\alpha \beta }^{-1}R_{\beta \nu } \end{eqnarray}$$

(17)

is the Fisher matrix of the maximum likelihood estimators |$\mathcal {P}_\mu ^{\prime }$|⁠,⁴ as shown in Thrane et al. (2009) and derived in Appendix A1 in the notation used in this work, also accounting for the presence of a gravitational wave signal. The dirty map is the inverse-noise weighted representation of the sky-distributed gravitational wave power as seen through the response of the pulsars (Thrane et al. 2009; Ali-Haïmoud et al. 2021; Pol et al. 2022; Agazie et al. 2023c), i.e. it is a ‘blurred’ picture of the true gravitational wave power distribution. The maximum likelihood estimator in equation (15) is referred to as the clean map|$\boldsymbol{\mathcal {P}}^{\prime }$|⁠.⁵ Equation (15) was derived in e.g. Thrane et al. (2009) for the weak-signal limit (i.e. where |$\boldsymbol{\zeta}$| is small), but it holds also as |$\boldsymbol{\zeta}$| increases, given that the corrections to |$\boldsymbol {\Sigma }$| are correctly applied. As |$\boldsymbol{\zeta}$| is increased, the maximum likelihood solution changes from equation (15). We return to this detail later.

3.7 Regularization

The inversion of the Fisher matrix in equation (15) involves an interesting subtlety. The Fisher matrix can be diagonalized into eigenmodes, each of which represents some patch of sky (Mitra et al. 2008; Thrane et al. 2009; Ali-Haïmoud et al. 2021); see Fig. C1 in the Appendix  C. However, due to the irregular sky distribution and sensitivity of the pulsars, some eigenmodes can be measured far more precisely than others. As a result, naive inversion of the Fisher matrix produces results that are dominated by the uncertainty associated with the least well-measured modes.

We therefore regularize the Fisher matrix with a singular value decomposition scheme outlined in Thrane et al. (2009) and Abadie et al. (2011). We project out the least well-measured modes in and keep only |$N_\text{cut-off}$| eigenmodes, i.e. formally, we set the smallest eigenvalues of |$\bf{M}$| to infinity. With these poorly measured modes removed, we calculate the inverse of the regularized Fisher matrix denoted |$\tilde {M}_{\mu \nu }^{-1}$|⁠.

Choosing the number of modes to keep |$N_\text{cut-off}$| is a trade-off. By throwing out some modes, the sky map becomes biased because these rejected modes are no longer included in our sky maps. However, we obtain smaller uncertainties on the modes that we choose to keep. We experimented with different cut-offs and chose |$N_\text{cut-off}=32$| based on our ability to produce reliable reconstructions with simulated data. We discuss below how our results change for different values of |$N_\text{cut-off}$| (see Appendix  C).

3.8 Clean map

Using the regularized Fisher matrix we obtain a regularized clean map (Thrane et al. 2009),

$$\begin{eqnarray} \tilde {\mathcal {P}}^{\prime }_\mu = \tilde {M}^{-1}_{\mu \nu } X_\nu , \end{eqnarray}$$

(18)

with associated uncertainty

$$\begin{eqnarray} \tilde {\sigma }^{\mathcal {P}^{\prime }}_\mu = \sqrt{\left({\tilde {\bf{M}}}^{-1} \bf {M} {\tilde {\bf{M}}}^{-1} \right)_{\mu \mu } } = \sqrt{ {\tilde {\bf{M}}}^{-1}_{\mu \mu }}, \end{eqnarray}$$

(19)

where the tildes denote regularized quantities. The clean map signal-to-noise ratio (S/N) is given by

$$\begin{eqnarray} {\rm S/N} = \frac{\tilde{{\cal P}}^{\prime }_{\hat{\Omega }}}{\tilde{\sigma }^{\cal P^{\prime }}_{\hat{\Omega }}} . \end{eqnarray}$$

(20)

Here the subscript |$\hat{\Omega }$| indicates that we have transformed the numerator and denominator from the spherical harmonic basis to the pixel basis.

Furthermore, throughout the paper a division of non-scalar quantities (mostly in the calculation of S/Ns) denotes an element-wise division.

Well-regularized, clean, S/N maps exhibit fluctuations of |${\approx} {\pm} 3$| about a mean value, which is determined by the magnitude of the isotropic gravitational wave signal. Large positive excursions from the mean can be indications of anisotropy, whereas large negative excursions from the mean are typically evidence of a problem in the upstream noise analysis. We use equation (20) to derive the main results of this paper.

3.9 Radiometer map

The clean map |$\tilde{\cal P}$| represents our best guess for how gravitational wave power is distributed in the sky with minimal assumptions. However, we can instead assume that there is precisely one point source somewhere in the sky and ask the question: What is the inferred strain power that would be coming from this one point source as a function of sky location? The answer to this question is the radiometer map (Mitra et al. 2008; Thrane et al. 2009). While it does not produce a reliable estimate for broad-scale anisotropy, it is optimized for point source detection. The radiometer map is given by

$$\begin{eqnarray} \eta _{\hat{\Omega }} = X_{\hat{\Omega }}/M_{\hat{\Omega }\hat{\Omega }}, \end{eqnarray}$$

(21)

with associated uncertainty

$$\begin{eqnarray} \sigma ^\eta _{\hat{\Omega }} = (M_{\hat{\Omega }\hat{\Omega }})^{-1/2} . \end{eqnarray}$$

(22)

As above, the |$\hat{\Omega }$| subscripts denote that we are working on a pixel basis. We emphasize that the radiometer map does not require regularization because we invert elements of the Fisher matrix rather than the entire Fisher matrix. This is because the radiometer map does not take into account the covariance between different patches of sky. The radiometer S/N is

$$\begin{eqnarray} {\rm S/N} = \frac{\eta _{\hat{\Omega }}}{\sigma _{\hat{\Omega }}^\eta }. \end{eqnarray}$$

(23)

3.10 Significance

The final step is to quantify the statistical significance of deviations from isotropy. There are numerous ways to estimate the significance of anisotropy in sky maps, each of which answers a different question. We opt to answer the question: Is the maximum value of |${\rm S/N}$| on our sky maps consistent with what one would expect from an isotropic background? This question reflects our prior belief that the first deviation from isotropy is likely to appear as a hotspot – a patch of sky with an anomalously high S/N.

We calculate a frequentist p-value. First, following Abadie et al. (2011), we generate a distribution of dirty maps; details are provided in Appendix A1. We then construct the clean map from each of these realizations following equation (15). In order to determine the significance of the fluctuations about an isotropic background, we subtract the monopole from these maps by setting the monopole component to zero (i.e. for the clean map |$\mathcal {P}_{\ell =0,m=0}=0$|⁠, and for the radiometer |$\eta _{\ell =0,m=0}=0$|⁠) before converting back to the pixel basis. For each clean map realization, we determine the minimum and maximum |${\rm S/N}$| across the whole map and combine these values into two histograms, in order to obtain the distributions of the global minimum and maximum |${\rm S/N}$|⁠. We compare the measured maximum |${\rm S/N}$| to the expected distribution in order to obtain a p-value. A small p-value |${\ll} 0 = 1 {{\, \rm per\, cent}}$| can indicate statistically significant anisotropy.

4 RESULTS

We create both regularized clean maps and radiometer maps for three frequency bins: |$7$|⁠, |$14 $|⁠, and |$21 \,\mathrm{n}\mathrm{Hz}$| (all multiples of the |$1/t_\text{obs}$|⁠). We focus on narrow-band maps since supermassive black hole binaries evolve slowly compared to our observation time, it is unlikely that the same source appears in more than one frequency bin.⁶ Each map indicates the sky position of all MPTA pulsars with white stars. The size of each star is scaled by the inverse of the root-mean-square of the pulsar’s residuals so that bigger stars provide more information about the stochastic background than small stars.

In Fig. 1, we plot the regularized S/N clean maps, which provide minimum-assumption reconstructions of the gravitational wave sky. The p-values for maximum and minimum |${\rm S/N}$| are shown in the ‘clean map’ column of Table 2. The p-values for the minimum |${\rm S/N}$| are between 0.15 and 0.96, which is encouraging because small p-values for minimum |${\rm S/N}$| typically indicate a problem in the data analysis pipeline.⁷ While two of the maximum |${\rm S/N}$| values are typical for isotropic backgrounds, the |$7 \,\mathrm{n}\mathrm{Hz}$| maximum |${\rm S/N}$| is somewhat unusual with |$p=0.03$|⁠.

Figure 1.

Regularized clean S/N maps of the gravitational wave power across the three frequency bins. The white stars are the sky locations of the pulsars in the MPTA, where the size corresponds to the inverse of the residuals’ root-mean-square of each pulsar. All three maps are consistent with an isotropic background signal at the |${\lesssim} 2.3\sigma$| level. There is a hotspot in the |$7 \, \mathrm{nHz}$| frequency bin with a modest statistical significance of |$p=0.03$| located at RA 1^h, Dec. |$-70 \,{}^{\circ }$|⁠.

Table 2.

Quantifying radiometer map and clean map anisotropy. Frequentist p-values for the minimum and maximum |${\rm S/N}$| in each frequency bin.

	Radiometer map		Clean map
Frequency (nHz)	|$p(\mathrm{max})\, ^{a}$|	|$p(\mathrm{min}) \, ^{b}$|	|$p(\mathrm{max}) \, ^{a}$|	|$p(\mathrm{min}) \, ^{b}$|
7	0.015	0.15	0.03	0.37
14	0.72	0.78	0.62	0.89
21	0.42	0.94	0.67	0.96

^aShorthand notation for |$p({\rm S/N}_\mathrm{max}^\mathrm{means})$|⁠.

^bShorthand notation for |$p({\rm S/N}_\mathrm{min}^\mathrm{means})$|⁠.

In Fig. 2, we plot the radiometer maps, which are optimized for the detection of point sources. Comparing the radiometer to the clean maps in Fig. 1, we can see the impact of the cleaning. While there is a qualitative correspondence between hotspots on each map, the clean maps represent better estimates of the distribution of gravitational waves in the sky. The p-values for the maximum and minimum |${\rm S/N}$| are shown in the ‘radiometer map’ column of Table 2. As we saw with the clean maps, most of the p-values are |${\cal O}(1)$|⁠, suggesting that the data are consistent with an isotropic background. However, the |$7 \,\mathrm{n}\mathrm{Hz}$| maximum |${\rm S/N}$| is again somewhat unusual with |$p=0.015$|⁠.

Figure 2.

Radiometer S/N of the gravitational wave power across the three frequency bins. The white stars are the sky locations of the pulsars in the MPTA, where the size corresponds to the inverse of the residuals’ root-mean-square of each pulsar. There is a hotspot with a p-value of 0.015 in the first map. The results are consistent with isotropy at the |${\lesssim} 2.3\sigma$| level.

Further auxiliary data products characterizing the information content of the MPTA data set include the cleaned gravitational wave power sky map and a map of the PTA’s sensitivity towards gravitational waves from different areas of the sky. The respective maps are shown in Figs D1 and D2 in Appendix  D. For the interested expert reader, we also explain the mathematical details of the sensitivity estimation in Sections D1.2 and D2.

5 DISCUSSION

The presence of a hotspot in our |$7 \,\mathrm{n}\mathrm{Hz}$| maps at approximately RA 1^h, Dec. |$-70 \,{}^{\circ }$| is intriguing; however, the modest |$p=0.015$| statistical significance suggests that it could be due to a noise fluctuation. So in order to demonstrate the performance of our pipeline, we show in Appendix  B that we are able to recover injected point sources accurately, which illustrates our ability to detect point sources should they exist in the data.

The statistical significance of the hotspot is strongly affected by the PSR J2129−5721, less so by our most sensitive pulsars (see Appendix  E for details). We note that this is the pulsar that leads to the creation of the ALT model in the accompanying isotropic gravitational wave background analysis (see Miles et al. 2024b).

The hotspot is also, of course, affected by our choice of regularization cut-off (see Appendix  C for details). Reducing the number of eigenmodes in our analysis, the hotspot is still visible (and a second hotspot emerges), but the significance of these hotspots is diminished compared to our fiducial analysis. This is not unexpected: reducing |$N_\text{SpH}$| throws out some of our signal, but provides a better measurement of the remaining modes. Increasing the number of eigenmodes, the hotspot becomes indistinguishable from other fluctuations in the sky map. We interpret this to mean that the map is dominated by uncertainty in the poorly resolved eigenmodes. While we believe our choice of |$N_\text{SpH}$| is reasonable, we recommend additional work to develop a more objective method for determining |$N_\text{SpH}$|⁠, perhaps by optimizing the ability to resolve specific anisotropic signals.

If the hotspot is subsequently shown to be of astrophysical in origin, it could indicate that the gravitational wave background arises from supermassive black hole binaries. Our current angular resolution of |${\gtrsim} 23^\circ$| makes it difficult to speculate about potential host galaxies for a nearby supermassive black hole binary. Additional pulsars must be added to the MPTA network in order to improve the angular resolution of our sky maps. This would also serve to increase the overall sensitivity of the network. Since supermassive black hole binaries evolve slowly, the hotspot can be confirmed or ruled out with subsequent observations.

We have not attempted to account for cosmic variance, which means that our p-values are probably overestimating the presence of anisotropy. Fluctuations in the gravitational wave background from a statistically gravitational wave background (one that is isotropic on average) yield individual realizations with anisotropy (Allen 2023; Romano & Allen 2024). Future work is required to take into cosmic variance in our p-value calculations.

Another topic for future work is the inclusion of gravitational wave noise in our analysis, which becomes increasingly important as we depart from the low S/N regime. While our Fisher matrix includes contributions from gravitational wave noise, we do not take this into account when maximizing the likelihood function. As a result, our maximum likelihood estimators are only approximate solutions. Since our gravitational wave signal is not so large in each frequency bin, we expect the exact solutions to be qualitatively similar to our approximate solution. However, as the sensitivity of MPTA improves, it will become increasingly important to include a self-consistent treatment of gravitational wave noise.

One possible solution is a recursive approach. The initial maximum likelihood estimators for the spherical harmonics coefficients are used in the next iteration to calculate the overlap reduction function, leading to a new pulsar-pair correlation covariance matrix and an updated maximum likelihood estimate for |$\tilde {\mathcal {P}}^{\prime }$|⁠. This method is currently under development.

ACKNOWLEDGEMENTS

We thank Matthew Bailes for ongoing valuable support as chair of the MeerTime Collaboration. Furthermore, we appreciate the helpful comments on this work by Rutger van Haasteren, Huanchen Hu, Steve Taylor, Gosiya Curylo, and Levi Schult.

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory (SARAO), which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time–space reference systems department of the Paris Observatory. MeerTime observations used the PTUSE computing cluster. This cluster was funded in part by the Max-Planck-Institut für Radioastronomie (MPIfR) and the Max Planck Society (MPG). MeerTime data are stored and processed on the OzStar and Ngarrgu Tindebeek supercomputers, operated by the Swinburne University of Technology (SUT). This project utilized the MeerTime data portal, which was supported by Nick Swainston and the Astronomy Data and Computing Services (ADACS) team. We acknowledge and pay respects to the Elders and Traditional Owners of the land on which the Australian institutions stand, the Bunurong and Wurundjeri Peoples of the Kulin Nation.

We acknowledge support from the Australian Research Council (ARC) Centres of Excellence for Gravitational Wave Discovery (OzGrav) CE170100004 and CE230100016. KG, DJC, FA, MJK, and VVK acknowledge continuing valuable support from the Max Planck Society. KG acknowledges support from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the University of Bonn and University of Cologne. RSN acknowledges support from the Astronomical Society of Australia. RMS acknowledges support through ARC Future Fellowship FT190100155. FA acknowledges that part of the research activities described in this paper were carried out with the contribution of the NextGenerationEU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and Creation of Research Infrastructures, Project IR0000034 – ‘STILES – Strengthening the Italian Leadership in ELT and SKA’. Pulsar research at Jodrell Bank Centre for Astrophysics was supported by an Science and Technology Facilities Council (STFC) consolidated grant (ST/T000414/1 and ST/X001229/1). MKr acknowledges support by the CAS (Chinese Acadamy of Science)-MPG Legacy Programme. AP acknowledges financial support from the European Research Council (ERC) starting grant ‘GIGA’ (A Gamma-ray Infrastructure to Advance Gravitational Wave Astrophysics, grant agreement no. 101116134) and through the Dutch Research Council (NWO)-I Veni fellowship. JS acknowledges funding from the South African Research Chairs Initiative of The Department of Science and Technology and the National Research Foundation of South Africa. VVK acknowledges financial support from the ERC starting grant ‘COMPACT’ (Understanding gravity using a COMprehensive search for fast-spinning Pulsars And CompacT binaries, grant agreement no. 101078094). PG acknowledges support through the SUT stipend Swinburne University Postgraduate Research Award (SUPRA). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

This publication made use of open source python libraries including numpy (Harris et al. 2020) and matplotlib (Hunter 2007), as well as the python-based pulsar analysis packages, enterprise (Ellis et al. 2020) and enterprise_extensions (Taylor et al. 2021).

DATA AVAILABILITY

Data used in this analysis will be available on the MeerTime data portal https://pulsars.org.au/.

Footnotes

1

Numerically evaluating the integral in equation (1) requires a suitable choice of the discretization resolution. In this analysis, we use healpy (Górski et al. 2005; Zonca et al. 2019) to create equal-area pixels, governed by the parameter nside, which takes values of |$2^n$|⁠, where |$n \in \mathbb {N}^+$|⁠. If the sky patch size is too similar to the spherical harmonics resolution, this under-resolving strongly impairs the sky maps, i.e. they vary as a function of nside until the resolution is adequate. For the MPTA data set, we find it converges with nside|$\ge 16$|⁠.

2

In the literature, the autocovariance matrix is often also denoted as |$P_\mathrm{ a}$| (e.g. Chamberlin et al. 2015; Pol et al. 2022), but we adopt the notation of Vigeland et al. (2018) as it is the most commonly used reference, and also matches the notation of the analysis script that we use.

3

We avoid calling |$\Sigma$| a covariance matrix – which it is – to avoid confusion with |$\bf {C}$| and |$\bf {S}$|⁠, which are also covariance matrices. The former is a covariance matrix for the correlations, while the latter are covariance matrices for the residuals.

4

The apostrophe indicates the maximum likelihood estimate of a quantity.

5

Our clean maps have units of |$\text{correlation}^{2}\,\mathrm{sr}^{-1}$|⁠. The analogous maps produced by LIGO–Virgo (Abadie et al. 2011) have units of |$\text{strain}^{2}\,\mathrm{Hz}^{-1}\,\mathrm{sr}^{-1}$|⁠.

6

Assuming a circular orbit. Gravitational wave emission is no longer monochromatic over the orbital period in the presence of eccentricity (Peters & Mathews 1963).

7

As outlined in Section 3.1, negative power in the map is due to noise fluctuations and therefore should not be statistically significant.

8

These equations were simultaneously and independently derived in Gersbach et al. (2024), and were published while this paper was in preparation.

9

We set the recovery |$\log _{10}A_\mathrm{GWB}$| to different values until the |${\rm S/N}$| range of the clean map matches the |${\rm S/N}$| range of the clean map of the real data set. We find that recovering the isotropic gravitational wave background signal with a common red noise specified by |$\log _{10}A_\mathrm{GWB}= -14.0$| produced a reasonable |${\rm S/N}$| range as shown in Fig. C3. Hence, we adopted this value for the analysis across all simulations.

APPENDIX A: ADDITIONAL METHODOLOGY DETAILS: BEYOND THE WEAK SIGNAL APPROXIMATION

A1 Fisher matrix

Thrane et al. (2009) showed that in the weak signal limit, |$\bf {M} = \bf {R}^{\rm T}\boldsymbol {\Sigma }^{-1}\bf {R}$| is the Fisher matrix of the maximum likelihood estimators |$\boldsymbol{\mathcal {P}}^{\prime }$|⁠. This is appropriate in the context of audio-band detectors like LIGO because the stochastic background is extremely small compared to the instrumental noise in any given measurement. The audio-band stochastic background is only detectable by stacking a large number of measurements to dig beneath the noise (Romano & Cornish 2017).

However, the situation is different in pulsar timing where the gravitational wave signal can be comparable to the pulsar noise. In this work, we include some corrections that arise when we relax the small-signal assumption. We here show that |$\mathrm{Cov}({\boldsymbol{\mathcal {P}}^{\prime },\boldsymbol{\mathcal {P}}^{\prime }}) = \left[\bf {R}^{\rm T}\boldsymbol {\Sigma }^{-1}\bf {R}\right] = \bf {M}^{-1}$| still holds irrespective of the S/N.

We begin by simplifying the expression for the covariance as

$$\begin{eqnarray} \mathrm{Cov}({\mathcal {P}^{\prime }_\mu ,\mathcal {P}^{\prime }_\nu }) = M^{-1}_{\mu \kappa }M^{-1}_{\nu \lambda } \mathrm{Cov}({X_\kappa ,X_\lambda }) = M^{-1}_{\mu \kappa }M^{-1}_{\nu \lambda } \left[\langle X_\kappa X_\lambda \rangle -\langle X_\kappa \rangle \langle X_\lambda \rangle \right], \end{eqnarray}$$

(A1)

where we used the definition equation (15) to substitute |$\boldsymbol{\mathcal {P}}^{\prime }$|⁠. In the next step, we further break down the covariance of the dirty map X by substituting equation (16). This yields for the expectation values in equation (A1):

$$\begin{eqnarray} \langle X_\kappa \rangle &= R^{\rm T}_{\kappa \alpha }\Sigma ^{-1}_{\alpha \beta }\langle \rho _\beta \rangle , \end{eqnarray}$$

(A2)

$$\begin{eqnarray} \langle X_\kappa X_\lambda ^{\rm T} \rangle &=& R^{\rm T}_{\kappa \alpha }\Sigma ^{-1}_{\alpha \beta } R_{\lambda \delta }\Sigma ^{-1}_{\delta \epsilon }\langle \rho _\beta \rho _\epsilon ^{\rm T} \rangle \\&=& R^{\rm T}_{\kappa \alpha }\Sigma ^{-1}_{\alpha \beta } R^{\rm T}_{\lambda \delta }\Sigma ^{-1}_{\delta \epsilon } \left[\Sigma _{\beta \epsilon } + \langle \rho _\beta \rangle \langle \rho _\epsilon ^{\rm T} \rangle \right]. \end{eqnarray}$$

(A3)

Putting everything together, the components of the covariance matrix of the dirty map read

$$\begin{eqnarray} \mathrm{Cov}({X_\kappa ,X_\lambda }) &=& R^{\rm T}_{\kappa \alpha }\Sigma ^{-1}_{\alpha \beta } R_{\lambda \delta }\Sigma ^{-1}_{\delta \epsilon } \left[\Sigma _{\beta \epsilon } + \langle \rho _\beta \rangle \langle \rho _\epsilon ^{\rm T} \rangle -\langle \rho _\beta \rangle \langle \rho _\epsilon ^{\rm T} \rangle \right] \\&=& R^{\rm T}_{\kappa \alpha }\Sigma ^{-1}_{\alpha \beta } R_{\lambda \delta }\Sigma ^{-1}_{\delta \epsilon }\Sigma _{\beta \epsilon } \\&=& R^{\rm T}_{\kappa \alpha }\Sigma ^{-1}_{\alpha \beta }R_{\delta \lambda } \delta _{\delta \beta } \\&=& R^{\rm T}_{\kappa \alpha }\Sigma ^{-1}_{\alpha \beta }R_{\beta \lambda } \\&=& M_{\kappa \lambda }. \end{eqnarray}$$

(A4)

Inserting back into equation (A1), we find

$$\begin{eqnarray} \mathrm{Cov}({\mathcal {P}^{\prime }_\mu ,\mathcal {P}^{\prime }_\nu }) = M^{-1}_{\mu \kappa }M^{-1}_{\nu \lambda } \mathrm{Cov}({X_\kappa ,X_\lambda }) = M^{-1}_{\mu \kappa }\underbrace{M^{-1}_{\nu \lambda } M_{\kappa \lambda }}_{=\delta _{\kappa \nu }} = M^{-1}_{\mu \nu }. \end{eqnarray}$$

(A5)

Thus, the expression for |$\bf {M}$| as defined already in Thrane et al. (2009) and Pol et al. (2022), |$\bf {M} = \bf {R}^{\rm T}\boldsymbol {\Sigma }^{-1}\bf {R}$|⁠, holds as the Fisher matrix of the clean map, regardless of the signal strength. The correction with respect to the weak signal limit is completely absorbed in the correct calculation of the correlation covariance matrix |$\boldsymbol {\Sigma }$|⁠, i.e. the only change to the formalism lies in the calculation of |$\boldsymbol {\Sigma }$|⁠.

A2 Corrections to the sigma matrix

When we relax the weak-signal approximation, the off-diagonal entries of the cross-correlation covariance matrix are not negligible anymore and need to be included in the analysis. A general mathematical expression for the full matrix has been derived in Allen & Romano (2023). By distinguishing between different cases of pulsar-pair constituents, this general expression can be simplified and written in terms of known quantities such as the auto- and cross-correlation matrices of the pulsars. This is also the key step towards any implementation of the full matrix. In the following, we therefore derive these simplified expressions in order to bridge the gap between the work done by Allen & Romano (2023) and the code provided for the analysis done in this work.⁸

We aim to calculate the covariance of the cross-correlations, i.e.

$$\begin{eqnarray} \Sigma _{\alpha \beta } = \mathrm{Cov}({\rho _\alpha ,\rho _\beta }) = \langle \rho _\alpha \rho _\beta \rangle -\langle \rho _\alpha \rangle \langle \rho _\beta \rangle . \end{eqnarray}$$

(A6)

First, we adopt index notation, explicitly writing out the vector and matrix products in equation (10). In order to maintain a human level of readability, all pulsar indices |$a,b,c,d$| are upper indices, while the lower indices |$i,j,k,l,m,n,p,q$| run over the number of ToAs in each pulsar. We furthermore adopt the Einstein summation convention for lower indices.

In this notation, we can rewrite the pulsar cross-correlation calculation from equation (10) as

$$\begin{eqnarray} \rho _{ab} = \sigma _{ab}^2 \delta t^a_i (C^{-1})^a_{ij} \hat {S}^{ab}_jk (C^{-1})^b_{kl} \delta t^b_l, \end{eqnarray}$$

(A7)

where we used equation (11) to replace |$\left({\rm tr} [(\bf{C}_a^{-1}\hat {\bf{S}}_{ab}\bf{C}_b^{-1}\hat {\bf{S}}_{ab})]\right)^{-1} = \sigma _{ab}^2$|⁠:

$$\begin{eqnarray} \mathrm{Cov}({\rho _{ab},\rho _{cd}}) &=& \sigma _{ab}^2\sigma _{cd}^2 \left[(C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl}\right]\,\,\left[(C^{-1})_{mn}^c \hat {S}_{np}^{cd}(C^{-1})^d_{pq}\right] \underbrace{\left\lbrace \langle \delta t^a_i \delta t^b_l \delta t^c_m \delta t^d_q \rangle -\langle \delta t^a_i \delta t^b_l \rangle \langle \delta t^c_m \delta t^d_q \rangle \right\rbrace }_{{\substack{\sim \,\, \langle ab \rangle \langle cd \rangle + \langle ac \rangle \langle bd \rangle + \langle ad \rangle \langle bc \rangle \, \, -\langle ab \rangle \langle cd \rangle \\= \langle ac \rangle \langle bd \rangle + \langle ad \rangle \langle bc \rangle }}} \\&=& \sigma _{ab}^2\sigma _{cd}^2 \left[(C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl}\right]\,\,\left[(C^{-1})_{mn}^c \hat {S}_{np}^{cd}(C^{-1})^d_{pq}\right] \underbrace{\left\lbrace \langle \delta t^a_i \delta t^c_m \rangle \langle \delta t^b_l \delta t^d_q \rangle + \langle \delta t^a_i \delta t^d_q \rangle \langle \delta t^b_l \delta t^c_m \rangle \right\rbrace }_{\equiv \,\, \mathfrak {T}}. \end{eqnarray}$$

(A8)

Here, we use Isserlis’ theorem (Isserlis 1918; Wick 1950) in order to express the four-point correlation as a sum of products of two-point correlations, i.e.|$\langle abcd \rangle =\langle ab \rangle \langle cd \rangle + \langle ac \rangle \langle bd \rangle + \langle ad \rangle \langle bc \rangle$|⁠. Implicitly we introduce the shorthand notation |$\langle ab \rangle = \langle \delta t^a_i \delta t^b_l \rangle$|⁠, where

$$\begin{eqnarray} \langle aa \rangle _{ij} \equiv \langle \delta t^a_i \delta t^a_j \rangle = C^{a}_{ij}, \qquad \qquad \langle ab \rangle _{ij} \equiv \langle \delta t^a_i \delta t^b_j \rangle = S^{ab}_{ij} = A^2_\mathrm{GWB}\Gamma ^{ab}\hat {S}^{ab}_{ij}. \end{eqnarray}$$

(A9)

Depending on the pulsars making up each pair, the individual correlations in |$\mathfrak {T}$| contribute either as the autocovariance matrix |$\bf {C}_{a}$| or the cross-covariance matrix |$\bf {S}_{ab}$|⁠. In the decomposition of the cross-covariance matrix, |$\Gamma ^{ab}$| is the overlap reduction function derived from the pulsar-pair correlations, which does not necessarily have to be the Hellings–Downs correlation function, depending on the angular gravitational wave power distribution.

As pointed out above, all entries of the cross-correlation covariance matrix |$\boldsymbol {\Sigma }$| can be classified into three different cases derivation wise (both pulsars in the pairs are the same, only one matches, all are different), which turn into five cases upon numerical implementation (the one-match cases are related to each other by symmetry of |$\Gamma _{ab}$|⁠, |$\hat {\bf{S}}_{ab}$|⁠, and |$\bf{C}_a$|⁠):

• two-match case: |$a=c$|⁠, |$b=d$|⁠: |$\quad \mathfrak {T}\longrightarrow \langle aa \rangle _{im}\langle bb \rangle _{lq} + \langle ab \rangle _{iq}\langle ba \rangle _{lm} = C^a_{im}C^b_{lq} + A_\mathrm{GWB}^4{\Gamma ^{ab}}^2\hat {S}^{ab}_{iq} \hat {S}^{ab}_{lm}$|

  $$\begin{eqnarray} \qquad \qquad \qquad \mathrm{Cov}({\rho _{ab},\rho _{ab}}) &=& \sigma _{ab}^2\sigma _{ab}^2 \bigg \lbrace \underbrace{C^a_{im}(C^{-1})^a_{ij}}_{= \delta _{mj}}\hat {S}_{jk}^{ab}(C^{-1})^b_{kl}(C^{-1})_{mn}^a \hat {S}_{np}^{ab}\underbrace{(C^{-1})^b_{pq}C^b_{lq}}_{=\delta _{pl}} \\&& \qquad \qquad \qquad + A_\mathrm{GWB}^4{\Gamma ^{ab}}^2\, \hat {S}^{ab}_{iq}(C^{-1})^a_{ij}\hat {S}_{jk}^{ab}(C^{-1})^b_{kl}\hat {S}^{ab}_{lm}(C^{-1})_{mn}^a \hat {S}_{np}^{ab}(C^{-1})^b_{pq} \bigg \rbrace \\&=& \sigma _{ab}^2\sigma _{ab}^2 \left\lbrace \mathrm{tr}[\bf{C}^{-1}_a\hat {\bf{S}}_{ab}\bf{C}^{-1}_b\hat {\bf{S}}_{ba}] \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^4\Gamma _{ab}^2 \, \mathrm{tr}[\bf{C}^{-1}_a\hat {\bf{S}}_{ab}\bf{C}^{-1}_b\hat {\bf{S}}_{ba}\bf{C}^{-1}_a\hat {\bf{S}}_{ab}\bf{C}^{-1}_b\hat {\bf{S}}_{ba}] \right\rbrace , \end{eqnarray}$$

  (A10)

• one-match case

  • |$a=c$|⁠, |$b\ne d$|⁠: |$\quad \mathfrak {T}\longrightarrow \langle aa \rangle _{im}\langle bd \rangle _{lq} + \langle ad \rangle _{iq}\langle ba \rangle _{lm} = A^2_\mathrm{GWB}\Gamma ^{bd}C^a_{im}\hat {S}^{bd}_{lq} + A_\mathrm{GWB}^4\Gamma ^{ab}\Gamma ^{ad}\hat {S}^{ad}_{iq} \hat {S}^{ba}_{lm}$|

    $$\begin{eqnarray} \qquad \qquad \qquad \mathrm{Cov}({\rho _{ab},\rho _{ad}}) &=& \sigma _{ab}^2\sigma _{ad}^2 \left\lbrace A^2_\mathrm{GWB}\Gamma ^{bd}\,\, C^a_{im}(C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} \hat {S}^{bd}_{lq} (C^{-1})_{mn}^a \hat {S}_{np}^{ad}(C^{-1})^d_{pq} \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^4\Gamma ^{ab}\Gamma ^{ad}\,\,\hat {S}^{ad}_{iq} (C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} \hat {S}^{ba}_{lm} (C^{-1})_{mn}^a \hat {S}_{np}^{ad}(C^{-1})^d_{pq} \right\rbrace \\&=&\sigma _{ab}^2\sigma _{ad}^2 \left\lbrace A^2_\mathrm{GWB}\Gamma _{bd} \mathrm{tr}[\bf{C}^{-1}_a\hat {\bf{S}}_{ab}\bf{C}^{-1}_b\hat {\bf{S}}_{bd}\bf{C}^{-1}_d\hat {\bf{S}}_{da}] \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^4\Gamma ^{ab}\Gamma ^{ad} \,\,\mathrm{tr}[\bf{C}^{-1}_a \hat {\bf{S}}_{ab}\bf{C}^{-1}_b \hat {\bf{S}}_{ba} \bf{C}^{-1}_a \hat {\bf{S}}_{ad}\bf{C}^{-1}_d\hat {\bf{S}}_{da}] \right\rbrace , \end{eqnarray}$$

    (A11)
  • |$a\ne c$|⁠, |$b = d$|⁠: |$\quad \mathfrak {T}\longrightarrow \langle ac \rangle _{im}\langle bb \rangle _{lq} + \langle ab \rangle _{iq}\langle bc \rangle _{lm} = A_\mathrm{GWB}^2\Gamma ^{ac}\hat {S}^{ac}_{im}C^b_{lq} + A_\mathrm{GWB}^4\Gamma ^{ab}\Gamma ^{bc}\hat {S}^{ab}_{iq} \hat {S}^{bc}_{lm}$|

    $$\begin{eqnarray} \qquad \qquad \qquad \mathrm{Cov}({\rho _{ab},\rho _{cb}}) &=& \sigma _{ab}^2\sigma _{cb}^2 \left\lbrace A_\mathrm{GWB}^2\Gamma ^{ac} \hat {S}^{ac}_{im} (C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} (C^{-1})_{mn}^c \hat {S}_{np}^{cb}(C^{-1})^b_{pq}C^b_{lq} \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^4\Gamma ^{ab}\Gamma ^{bc} \hat {S}^{ab}_{iq} (C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} \hat {S}^{bc}_{lm} (C^{-1})_{mn}^c \hat {S}_{np}^{cb}(C^{-1})^b_{pq} \right\rbrace \\&=& \sigma _{ab}^2\sigma _{cb}^2 \left\lbrace A_\mathrm{GWB}^2\Gamma ^{ac} \mathrm{tr}[ \bf{C}^{-1}_a \hat {\bf{S}}_{ab}\bf{C}^{-1}_{b} \hat {\bf{S}}_{bc} \bf{C}^{-1}_c \hat {\bf{S}}_{ca}] \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^4\Gamma ^{ab}\Gamma ^{bc} \bf{C}^{-1}_a\hat {\bf{S}}_{ab} \bf{C}^{-1}_b\hat {\bf{S}}_{bc} \bf{C}^{-1}_c\hat {\bf{S}}_{cb} \bf{C}^{-1}_b\hat {\bf{S}}_{ba} \right\rbrace , \end{eqnarray}$$

    (A12)
  • |$a = d$|⁠, |$b\ne c$|⁠: |$\quad \mathfrak {T}\longrightarrow \langle ac \rangle _{im}\langle ba \rangle _{lq} + \langle aa \rangle _{iq}\langle bc \rangle _{lm} = A_\mathrm{GWB}^4\Gamma ^{ac}\Gamma ^{ab}\hat {S}^{ac}_{im}\hat {S}^{ba}_{lq} + A_\mathrm{GWB}^2\Gamma ^{bc} C^a_{iq} \hat {S}^{bc}_{lm}$|

    $$\begin{eqnarray} \qquad \qquad \qquad \mathrm{Cov}({\rho _{ab},\rho _{ca}}) &=& \sigma _{ab}^2\sigma _{ca}^2 \left\lbrace A_\mathrm{GWB}^4\Gamma ^{ac}\Gamma ^{ab} (C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} \hat {S}^{ba}_{lq} \hat {S}^{ac}_{im} (C^{-1})_{mn}^c \hat {S}_{np}^{ca}(C^{-1})^a_{pq} \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^2\Gamma ^{bc} (C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} \hat {S}^{bc}_{lm} (C^{-1})_{mn}^c \hat {S}_{np}^{ca}(C^{-1})^a_{pq} C^a_{iq}\right\rbrace \\&=& \sigma _{ab}^2\sigma _{ca}^2 \left\lbrace A_\mathrm{GWB}^4\Gamma ^{ac}\Gamma ^{ab} \mathrm{tr}[\bf{C}^{-1}_a \hat {\bf{S}}_{ab} \bf{C}^{-1}_b \hat {\bf{S}}_{ba} \bf{C}^{-1}_a \hat {\bf{S}}_{ac} \bf{C}^{-1}_c \hat {\bf{S}}_{ca}] \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^2\Gamma ^{bc} \mathrm{tr}[\bf{C}^{-1}_a \hat {\bf{S}}_{ab}\bf{C}^{-1}_b \hat {\bf{S}}_{bc} \bf{C}^{-1}_c \hat {\bf{S}}_{ca}] \right\rbrace , \end{eqnarray}$$

    (A13)
• no-match case: |$a \ne b \ne c \ne d$|⁠: |$\quad \mathfrak {T}\longrightarrow \langle ac \rangle _{im}\langle bd \rangle _{lq} + \langle ad \rangle _{iq}\langle bc \rangle _{lm} = A_\mathrm{GWB}^4\Gamma _{ac}\Gamma _{bd} \hat {S}^{ac}_{im}\hat {S}^{bd}_{lq} + A_\mathrm{GWB}^4\Gamma ^{ad}\Gamma ^{bc} \hat {S}^{ad}_{iq} \hat {S}^{bc}_{lm}$|

  $$\begin{eqnarray} \qquad \qquad \qquad \qquad \mathrm{Cov}({\rho _{ab},\rho _{cd}}) &=& \sigma _{ab}^2\sigma _{cd}^2 \left\lbrace A_\mathrm{GWB}^4\Gamma _{ac}\Gamma _{bd} (C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} \hat {S}^{ac}_{im} (C^{-1})_{mn}^c \hat {S}_{np}^{cd}(C^{-1})^d_{pq}\hat {S}^{bd}_{lq} \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^4\Gamma ^{ad}\Gamma ^{bc} (C^{-1})^a_{ij} \hat {S}_{jk}^{ab}(C^{-1})^b_{kl} \hat {S}^{bc}_{lm} (C^{-1})_{mn}^c \hat {S}_{np}^{cd}(C^{-1})^d_{pq} \hat {S}^{ad}_{iq} \right\rbrace \\&=& \sigma _{ab}^2\sigma _{cd}^2 \left\lbrace A_\mathrm{GWB}^4\Gamma _{ac}\Gamma _{bd} \mathrm{tr}[ \bf{C}^{-1}_a\hat {\bf{S}}_{ab} \bf{C}^{-1}_b\hat {\bf{S}}_{bd} \bf{C}^{-1}_d\hat {\bf{S}}_{dc} \bf{C}^{-1}_c\hat {\bf{S}}_{ca}] \right. \\&& \qquad \qquad \qquad \left. + A_\mathrm{GWB}^4\Gamma ^{ad}\Gamma ^{bc} \mathrm{tr}[ \bf{C}^{-1}_a\hat {\bf{S}}_{ab} \bf{C}^{-1}_b\hat {\bf{S}}_{bc} \bf{C}^{-1}_c\hat {\bf{S}}_{cd} \bf{C}^{-1}_d\hat {\bf{S}}_{da} ] \right\rbrace . \end{eqnarray}$$

  (A14)

Comparing the terms contributing to the main diagonal entries (equation A10) to the expression for the cross-correlation uncertainties |$\sigma _{\alpha \beta }$| as given in equation (11) we find that

$$\begin{eqnarray} \mathrm{Cov}({\rho _{ab},\rho _{ab}}) = \sigma _{ab}^2 + A_\mathrm{GWB}^4\Gamma _{ab}^2 \, \mathrm{tr}[\bf{C}^{-1}_a\hat {\bf{S}}_{ab}\bf{C}^{-1}_b\hat {\bf{S}}_{ba}\bf{C}^{-1}_a\hat {\bf{S}}_{ab}\bf{C}^{-1}_b\hat {\bf{S}}_{ba}]. \end{eqnarray}$$

(A15)

Following from that, |$\varsigma$| in equation (14) is given by

$$\begin{eqnarray} \boldsymbol {\varsigma } = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}A_\mathrm{GWB}^4\Gamma _{ab}^2 \, \mathrm{tr}[\bf{C}^{-1}_a\hat {S}_{ab}\bf{C}^{-1}_b\hat {S}_{ba}\bf{C}^{-1}_a\hat {S}_{ab}\bf{C}^{-1}_b\hat {S}_{ba}] &\text{if } a=c, b=d,\\\mathrm{equations \ (A11)\!-\!(A14)} & \text{otherwise.} \end{array}\right. \end{eqnarray}$$

(A16)

APPENDIX B: DEMONSTRATION WITH SIMULATED DATA

In order to test the behaviour of our analysis pipeline and its capability to identify and localize individual point sources, we create simulated data sets with a similar |${\rm S/N}$| regime to what we observe in the actual data analysis.

Each simulated data set has the same observation times, pulsar positions, and pulsar timing model as the real data set. For each MPTA pulsar, we generate simulated ToAs, which include contributions from both white noise and a gravitational wave signal. Therefore, we create a .tim file containing idealized telescope site arrival times from the original .par and .tim files using the tempo2 formideal plugin. Then we initialize a libstempo tempo2pulsar object with the original .par-file and simulated .tim file, add white noise fluctuations to the efac and equad parameters, as well as different gravitational wave signals to the flat residuals using the libstempo toolbox. The resulting data set is analysed with our standard enterprise/matlab pipeline as laid out in Section 3, assuming only white noise and a power-law gravitational wave signal with |$\log _{10}A_\mathrm{GWB}=-14.0$| and |$\gamma _\mathrm{GWB}=4.31$| to be present in the data.

We investigate the presence of an isotropic background signal as well as the presence of point sources using the following different data sets:

• White noise + isotropic gravitational wave background. We inject a common red noise signal with a power-law power spectral density characterized by the values found in Miles et al. (2024b).⁹

• White noise + point sources. For all point source (continuous wave) injections, we adopt the same recovery amplitude and spectral index as used in the white noise + isotropic gravitational wave background model. We position all sources at a distance of |$1 \,\mathrm{M}\text{pc}$| and assume them to be monochromatic with |$f_\mathrm{gw}= 1 \times 10^{-8} \,\mathrm{Hz}$|⁠. The recovered |${\rm S/N}$| range is adjusted by changing the chirp mass of the simulated source(s). Throughout the analysis, we investigate the following scenarios:

  • single point source with |$\mathcal {M}_\mathrm{c}=1 \times 10^{8} \,\rm {M_ {\odot }}$| at a position of RA 18^h, Dec. |$-45 \,{}^{\circ }$|⁠,

  • single point source with |$\mathcal {M}_\mathrm{c}=1 \times 10^{8} \,\rm {M_ {\odot }}$| at a position of RA 6^h, Dec. |$45 \,{}^{\circ }$|⁠,

  • single point source with |$\mathcal {M}_\mathrm{c}=1 \times 10^{8} \,\rm {M_ {\odot }}$| at a position of RA 6^h, Dec. |$-45 \,{}^{\circ }$|⁠,

  • two point sources |$\mathcal {M}_\mathrm{c}=10^{7.9}\rm {M_ {\odot }}$| at positions RA 6^h, Dec. |$-45 \,{}^{\circ }$| and RA 18^h, Dec. |$-45 \,{}^{\circ }$|⁠.

Although these simulations are consistent with the real MPTA data set, we emphasize that we did not include any of the various chromatic noise processes present in the observed data. Thus, we can only qualitatively compare the simulation with the actual results presented in Section 4.

B1 Isotropic background

Fig. B1 shows the clean map |${\rm S/N}$| of the first three frequency bins for the simulated data set containing the isotropic background signal. Similar to the results obtained with the full MPTA 4.5-yr data set, these maps exhibit an overall positive map mean with fluctuations around that mean.

Figure B1.

Regularized clean S/N maps of the simulated data set containing white noise and an isotropic gravitational wave background signal The white stars are the sky locations of the pulsars in the MPTA, where the size of each star is inversely proportional to the mean root-mean-square of the best-fitting timing solution.

B2 Point source(s)

If the analysis pipeline is not only able to identify but also correctly localize the point source(s), we expect both the radiometer and clean map |${\rm S/N}$| to exhibit a spatially constraint area of significantly increased |${\rm S/N}$|⁠, located at a similar position as the simulated source. As the radiometer |${\rm S/N}$| map (see Section 3.9) is ideal for identifying point sources, we present both the resulting radiometer and clean map |${\rm S/N}$| for all three data sets mentioned above, each shown in Figs B2 and B3, respectively.

Figure B2.

Radiometer maps of the first frequency bin (⁠|$7 \, \mathrm{nHz}$|⁠) for the simulated data sets containing a single gravitational wave source. The injected position of the source is indicated with a cross-hair.

Figure B3.

Regularized clean S/N maps of the first frequency bin (⁠|$7 \,\mathrm{n}\mathrm{Hz}$|⁠) for the simulated data sets containing a single gravitational wave source. The injected position of the source is indicated with a cross-hair.

Unsurprisingly, the ‘hotspot’ at the injection location is only present in the first two frequency bins, as the frequency of the simulated sources falls between the first and second frequency bins of the MPTA. Thus we present only the sky maps for the first frequency bin. In each map, the position of the injected source(s) is indicated with grey cross-hairs.

As expected, sources that are injected in the Southern hemisphere are well localized in both the radiometer and clean |${\rm S/N}$| maps, but sources injected in the Northern hemisphere are not well resolved. That is, we are able to see a signal in the data when gravitational waves are in the Northern sky, but it is difficult to differentiate an isotropic background from a point source. This is in accordance with the general understanding (Taylor & Gair 2013; Ali-Haïmoud et al. 2021; Agazie et al. 2023c) that a patch of sky populated with more pulsars is in general more sensitive to the gravitational wave signal coming from that sky region.

APPENDIX C: REGULARIZATION

Here we explore the impact of the regularization cut-off on the resulting sky distributions. The first step is to examine the ‘eigenspectrum’ in Fig. C1 showing the eigenvalues of the Fisher matrix in descending order. Larger eigenvalues correspond to modes that we measure with relatively lower uncertainty. The three modes with the highest and lowest relative uncertainty are shown in Fig. C2. As noted in the main body of the paper, the choice of regularization cut-off is a balancing act, improving sensitivity to some modes by throwing out other modes, thereby introducing a bias. If not mentioned otherwise, we employ a cut-off value of 32 throughout the analysis.

Figure C1.

Singular value spetrum of the MPTA Fisher information matrix |$\bf{M}$| (cf. equation (17)) of the pulsar-pair correlations \rho|$\rho_{\alpha \beta}$| in the first frequency bin of the full MPTA 4.5 year data set.

Figure C2.

Fisher matrix modes of the full MPTA 4.5-yr data set. Top row: Three most constrained modes. Bottom row: Three least constrained modes.

In order to illustrate how the choice of cut-off affects our results, we create sky maps for |$i_\mathrm{max}\in [1, N_\mathrm{SpH}]$| for all three frequency bins and investigate the evolution of the clean map |${\rm S/N}$| distribution with varying |$i_\mathrm{max}$|⁠. We show three exemplary maps for |$i_\mathrm{max}=\lbrace 20, 32, 54\rbrace$| of the first bin in Fig. C3. Videos showing the map evolution across all cut-offs for the three frequency bins can be found at: 10.57891/j0vh-5g31.

Fig. C3 shows how our reconstructions of an isotropic background vary with different choices of regularization cut-off. Meanwhile, in Fig. C4, we observe how our sky maps vary with regularization cut-off for a point source (continuous wave) signal. The more modes that we keep, the more pronounced and sharp the single source becomes. At the same time, the off-source sky exhibits increasingly large noise fluctuations. Artefacts from the diffraction limit are visible in the rightmost plot (⁠|$i_\mathrm{max}=54$|⁠). Encouragingly, though, the location of the hotspot is stable as we vary the cut-off. In Fig. C5, we show how the actual MPTA clean maps vary with cut-off.

Figure C3.

Visualizations of the clean map |${\rm S/N}$| across the sky of the simulated white noise + gravitational wave background data set using different regularization cut-offs.

Figure C4.

Visualizations of the clean map |${\rm S/N}$| across the sky of a simulated point source (continuous wave) signal using different regularization cut-offs.

Figure C5.

Realizations of the clean map |${\rm S/N}$| as seen by the MPTA across the sky using different regularization cut-offs.

APPENDIX D: ADDITIONAL RESULTS

In this appendix, we show additional plots, which may be of interest to expert readers.

D1 Cleaned sky maps

D1.1 Clean gravitational wave power sky map

The clean maps associated with the S/N maps are shown in Fig. D1 across the three frequency bins. We used equation (18) to derive these maps. The bright spot seen in the S/N map in Fig. 1(a) can also be seen in the corresponding clean map, Fig. 1(a).

D1.2 Clean sensitivity map

In order to visualize the sensitivity, |$\boldsymbol{\mathcal {S}}$|⁠, of the MPTA as a function of |$(\text{RA}, \text{DEC})$|⁠, we define ‘sensitivity maps’, which show the mean |${\rm S/N}$| for a source located in some direction |$\hat{\boldsymbol{n}}$| with unit amplitude:

$$\begin{eqnarray} \boldsymbol{\mathcal {P}}_{\hat {n}}(\hat{\boldsymbol{\Omega} }) \equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 & \hat {\boldsymbol{\Omega} } = \hat {\boldsymbol{n}}, \\0 & \hat {\boldsymbol{\Omega} } \ne \hat {\boldsymbol{n}}. \end{array}\right. \end{eqnarray}$$

(D1)

Formally, the sensitivity map is a plot of

$$\begin{eqnarray} \mathcal {S}_\mathrm{clean}(\hat{\boldsymbol{n}}) \equiv \frac{\boldsymbol{\mathcal {P}}_\text{eff} }{ \tilde {\boldsymbol{\sigma} }^{\boldsymbol{\mathcal {P}}^{\prime }} }(\hat{\boldsymbol{n}}) = \frac{\tilde {\bf {M}}^{-1} \bf {M} \boldsymbol{\mathcal {P}}_{\hat {n}}}{\tilde {\boldsymbol{\sigma} }^{\boldsymbol{\mathcal {P}}^{\prime }}}(\hat{\boldsymbol{n}}). \end{eqnarray}$$

(D2)

Here, |$\mathcal {P}_\text{eff}(\hat{\boldsymbol{n}})$| is the effective strain power after accounting for the loss of power due to regularization. This quantity is calculated by constructing the expected dirty map |$\boldsymbol{X}_{\hat {n}} = \bf {M}\boldsymbol{\mathcal {P}}_{\hat {n}}$| and cleaning it with the regularized inverse of the Fisher matrix. Due to the pixel-wise evaluation, the sensitivity map can be written as

$$\begin{eqnarray} \boldsymbol{\mathcal {S}}_\mathrm{clean} = \frac{\text{diag}\left(\tilde {\bf {M}}^{-1} \bf {M}\right)}{\tilde {\boldsymbol{\sigma} }^{\mathcal {P}^{\prime }}}. \end{eqnarray}$$

(D3)

The resulting sensitivity maps are shown in Fig. D2. It is easier to detect gravitational wave signals associated with high-sensitivity patches of sky. As also shown by Taylor & Gair (2013) and Agazie et al. (2023c), a PTA is expected to be most sensitive in the sky area with the highest number density of pulsars, and most sensitive pulsars. This can also be seen in the maps in Fig. D2: they exhibit interesting structures showing that MPTA is considerably more sensitive to gravitational waves from the Southern hemisphere than the Northern hemisphere. As expected, the sensitivity distribution follows the number density of pulsars as well as the location of the most sensitive pulsars in the MPTA.

Figure D1.

Clean maps of the MPTA in its first three frequency bins. The units of the colour bar are |$\text{correlation}^{2} \, \mathrm{sr}^{-1}$|⁠.

Figure D2.

Clean sensitivity map (defined in equation D2) for the MPTA in its first three frequency bins.

D2 Radiometer sensitivity maps

For completeness, we can also define the radiometer sensitivity map in analogy to the clean map sensitivity, but evaluating the radiometer |${\rm S/N}$| (cf. equation 22). The corresponding map shows the mean |${\rm S/N}$| for a source in some direction |$\hat{\boldsymbol{\Omega }}$| with unit amplitude:

$$\begin{eqnarray} \eta _{\hat{\boldsymbol{n}}} \equiv \eta _{\hat{\boldsymbol{\Omega }}}(\hat{\boldsymbol{\Omega} }|\hat{\boldsymbol{n}}) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 & \hat{\boldsymbol{\Omega} }= \hat{\boldsymbol{n}}, \\0 & \hat{\boldsymbol{\Omega} }\ne \hat{\boldsymbol{n}} \end{array}\right. \end{eqnarray}$$

(D4)

as a function of |$\hat{\boldsymbol{n}}$|⁠. Formally, the radiometer sensitivity for the pixel at position |$\hat{\boldsymbol{n}}$| is given as |$\boldsymbol{\mathcal {S}}_\mathrm{radiometer}(\hat{\boldsymbol{n}}) \equiv \eta _{\hat{n}}/\boldsymbol{\sigma} ^\eta _{\hat{\Omega }}$|⁠. The full map can be expressed as

$$\begin{eqnarray} \boldsymbol{\mathcal {S}}_\mathrm{radiometer} = \frac{\mathbb {1}}{ \boldsymbol{\sigma} ^\eta _{\hat{\Omega }} } = \frac{\mathbb {1}}{ \left(M_{\hat{\Omega }\hat{\Omega}}\right)^{-1/2} }. \end{eqnarray}$$

(D5)

These maps are shown in Fig. D3.

Figure D3.

Radiometer sensitivity map (defined in equation (D5)) for the MPTA in its first three frequency bins.

APPENDIX E: DROPOUT ANALYSIS

In this appendix, we investigate how our sky maps depend on different pulsars by performing dropout analyses. We identify the pulsars with the potentially largest influence on the sky maps based on the mean ToA uncertainty, overall timing residual mean square, and the volatility of their noise models. The three pulsars with the lowest mean ToA uncertainty are (in ascending order) J0437−4715, J2241–5236, and J1909–3744. The respective clean map |${\rm S/N}$| sky maps of the first frequency bin are shown in the top row of Fig. E1. As expected, the |${\rm S/N}$| range changes slightly upon removing each pulsar. Nevertheless, the position and size of the patches with higher and lower |${\rm S/N}$| remain similar compared to the sky map of the full data set. Thus, we conclude that our result is robust against the influence of these most precisely measured pulsars.

Figure E1.

Sky maps of the clean map |${\rm S/N}$| in the first frequency bin resulting from the dropout analysis described in Appendix  E.

Additionally, we test the influence of J2129–5721 on the sky map, due to the peculiarity of its noise model. As pointed out by Miles et al. (2024b), it shows an extremely steep-spectrum achromatic red noise. On the one hand, the same noise is not present in pulsars that are at low angular separation to it (i.e. J1909–3744 and J2241–5236). On the other hand, 16 yr of PPTA observations did not allow us to constrain any achromatic red noise in this pulsar. Thus, Miles et al. (2024b) present the most sensitive solution excluding the achromatic red noise process for this pulsar. Since J2129–5721 is located at a key angular separation from other, more sensitive pulsars in the MPTA, we also test our results against removing it from the data set. The resulting sky map in the first frequency bin is shown in the bottom row of Fig. E1. Again we find that the |${\rm S/N}$| range is decreased compared to the full data set. But in this case, the position of the |${\rm S/N}$| patches changes significantly. Most striking, the hotspot located at RA 1^h, Dec. |$-70 \,{}^{\circ }$| almost completely vanishes, while the areas around RA 14^h, Dec. |$-10 \,{}^{\circ }$| and RA 17^h, Dec. |$-45 \,{}^{\circ }$| become more prominent as a ‘triangle’ of increased |${\rm S/N}$| areas in the Western hemisphere.