2.1 Dataset

We base our results on observations of MSPs collected with the MeerKAT radio telescope obtained as part of the MeerTIME Large Survey Project (Bailes et al. Reference Bailes2020). The MeerKAT radio telescope comprises 64 13.9-m diameter antennas located in the Karoo region of South Africa. Initially, 189 MSPs were observed as part of a census of the MSP population observable with MeerKAT (Spiewak et al. Reference Spiewak2022). The purpose of the census was to characterise the Galactic field MSPs at L-band and assess the achievable timing precision. Following the conclusion of the census, a sub-sample of MSPs has been observed regularly. The pulsars and integration times were chosen such that the largest sample of MSPs could be observed to sub-microsecond rms error within $\approx 300$ hr of observing per year. The current sample includes 83 MSPs. Currently, each pulsar is observed at an approximately two-week cadence. A minimum 256 s is allocated if the required timing precision can be achieved in less time, otherwise, MSPs are observed for a duration of up to $2\,048$ s. However, there are additional observations in the data set of longer duration, that is, up to several hours, as the data set also includes observations taken as part of the MeerTIME Relativistic Binary project (Kramer et al. Reference Kramer2021).

Nearly all of the observations we consider have been taken with the L-band receiver for recording data at 856 $-$ 1 712 MHz. We have used the fold mode dataset, which folds the digitised, channelised data stream at the topocentric period of the pulsar. The resulting observations have 8 s sub-integrations with four polarisation states and 1 024 phase bins per pulse period, written in PSRFITS file format. To account for bandpass roll-off effects, 48 frequency channels (each with $\sim 0.8\, \textrm{MHz}$ bandwidth) from the top and bottom of the band have been removed. For timing studies, the 928 channels are then averaged to 16 frequency channels with a bandwidth of 54 MHz each. As MSPs are faint radio sources, mitigation of radio frequency interference (RFI) is essential. RFI excision is performed using meerguard, an extension of coastguard software package (Lazarus et al. Reference Lazarus2016) tailored for MeerKAT pulsar timing data.

Even though only 83 pulsars are currently monitored as part of the MPTA, we have studied jitter in the larger sample of 89 MSPs that were initially included in the MPTA for long-term monitoring. To measure jitter noise, we first identified the highest S/N observations taken between 2019 February and 2023 June because these would be the most jitter-dominated. Using these observations, we created a frequency, time, and polarisation-averaged pulse template. The observation archives were then fully frequency and polarisation averaged into 8 s sub-integrations. Jitter noise is known to show a frequency-dependent evolution for bright pulsars such as PSR J0437 $-$ 4715 (Kulkarni et al. Reference Kulkarni2024), but for most cases, the usage of a frequency-averaged template does not affect the jitter noise (Parthasarathy et al. 2021).

For our other assessments of the MPTA sensitivity (including the DM noise modelling presented in Section 4.1 and observing strategy study presented in Section 4.2), we base our analysis on the MPTA 4.5-yr data set. Unlike the jitter analysis, the dataset uses sub-banded observations. A detailed description of the dataset is presented in Miles et al. (Reference Miles2025a). We summarise the differences below. In the 4.5-yr data release, the observations are averaged into 32 frequency channels. Here, we have used 16-channel frequency resolution. To account for pulse profile evolution with frequency, 16-channel frequency-dependent pulse profile models (referred to as pulse portraits) were used as templates and cross-matched with the sub-banded observations to create the ToAs. We restrict our analysis to ToAs derived from pulse profiles with an S/N greater than 10. These differences are not expected to impact our assessment of sensitivity or optimisation strategies.

2.2 Sensitivity curves of PTAs

To analyse the performance of GW detectors, it is often common to produce sensitivity curves, which show the sensitivity of a detector as a function of GW frequency. We use the formalism presented in Hazboun, Romano, & Smith (Reference Hazboun, Romano and Smith2019) and the software package hasasia to calculate the sensitivity of a PTA to a GW signal. We calculated the sensitivity curves in two cases: (a) Continuous GWs from a non-evolving circular binary and (b) an SGWB. For the continuous GW signal, a matched-filtering method can be used to determine a filter function that maximises the S/N. Since the source position and polarisation are unknown, the expected S/N averaged over the sky and inclination angle is

(1) \begin{equation} \langle \rho^{2} \rangle = 2 T_\textrm{obs} \int_{0}^{f_\textrm{Nyq}} \frac{S_\textrm{det, h}(f)}{S_\textrm{det, eff}(f)} df,\end{equation}

where $S_\textrm{det, eff}$ is the effective strain-noise power spectral density of the array, which is the weighted sum of the strain-noise power spectral density of individual pulsars. $S_\textrm{det, h}$ is the strain power spectral density of a deterministic monochromatic source of GWs, and $T_\textrm{obs}$ is the time span of the PTA dataset (Hazboun et al. Reference Hazboun, Romano and Smith2019). The individual pulsar strain power spectral densities are estimated considering individual pulsar noise characteristics and the response of a pulsar to a plane GW.

For an SGWB, an optimal cross-correlation statistic is used to determine the expected S/N while considering inter-pulsar correlations. The expected S/N ratio of the optimal statistic can be expressed as

(2) \begin{equation} \langle \rho_\textrm{OS}^{2} \rangle = 2 T_\textrm{obs} \int_{0}^{f_\textrm{Nyq}} \frac{S_\textrm{stoch, h}^{2}(f)}{S_\textrm{ stoch, eff}^{2}(f)} df,\end{equation}

where $S_\textrm{stoch, eff}(f)$ is the effective strain-noise power spectral density of the entire PTA, which is the sum of strain-noise power spectral densities of pulsar pairs, including the HD correlations and $S_\textrm{stoch, h}(f)$ is the strain-noise power spectral density of an SGWB (Hazboun et al. Reference Hazboun, Romano and Smith2019). $S_\textrm{det/stoch, eff}$ can be interpreted as the detection sensitivity curve and is usually represented in terms of a dimensionless ‘characteristic strain’, $h_\textrm{c}$ such that $h_\textrm{ c}(f) \equiv \sqrt{f S_\textrm{det/stoch, eff}(f)}$ .

For a PTA, the pulsar noise models include white noise (typically defined as modifications to formal ToA uncertainties) and red noise processes. Both are identified through detailed modelling of the pulsar timing data sets. The timing model fit removes power at the frequencies corresponding to yr $^{-1}$ and 2 yr $^{-1}$ because of fitting for pulsar sky position and parallax, and at lower frequencies because of fitting for a pulsar spin frequency and frequency derivative (see Figure 5). Similarly, power at frequencies corresponding to the reciprocal of the orbital periods is absorbed for pulsars in binary systems. A GWB will induce a red noise signal and contribute to the spectra of all pulsars (i.e. GW self-noise and could manifest as CURN). The GWB amplitude that can be detected for a given configuration of the PTA at a specified S/N can be estimated using equation (2).

Here, we analyse the detection sensitivity curves for an SGWB and a single circular binary source. For GWB self-noise, we consider three scenarios: a pessimistic amplitude of $1 \times 10^{-16}$ (Sesana Reference Sesana2013), an intermediate value of $2 \times 10^{-15}$ (consistent with the amplitude of the signal inferred from the PPTA analysis Reardon et al. Reference Reardon2023a) and a larger value of $5 \times 10^{-15}$ (more consistent with that inferred from the MPTA, Miles et al. Reference Miles2025b). We summarise different datasets used in this work and the assumed GWB self-noise amplitude in each case in Table 1.

Table 1. Summary of dataset used for different studies in this work. The length of the data span, relevant section, and assumed value of GWB amplitude for each study are shown.

To account for the complex noise models in MPTA, we have added additional noise processes to hasasia, in particular, DM variations and chromatic noise following the prescription given in Agazie et al. (Reference Agazie2024). We include DM variations as a red noise signal following a power spectral density expressed as

(3) \begin{equation} P_\textrm{DM}(f) = \frac{A_\textrm{DM}^{2}}{12 \pi^{2}} \bigg(\frac{f}{f_{c}}\bigg)^{-\gamma_\textrm{DM}} \bigg(\frac{\nu}{\nu_\textrm{ref}}\bigg)^{-4}\end{equation}

where $A_\textrm{DM}$ and $\gamma_\textrm{DM}$ are the amplitude and spectral index of DM noise power law spectrum, respectively, at a reference frequency $f_{c}$ at $1\textrm{yr}^{-1}$ . Similarly, for chromatic noise addition, we assume the power spectral density of the form

(4) \begin{equation} P_\textrm{chrom}(f) = \frac{A_\textrm{chrom}^{2}}{12 \pi^{2}} \bigg(\frac{f}{f_{c}}\bigg)^{-\gamma_\textrm{chrom}} \bigg(\frac{\nu}{\nu_\textrm{ref}}\bigg)^{-2\beta}\end{equation}

where $\beta$ is the chromatic index, a free parameter encompassing the frequency-dependent multi-path propagation delays. We did not incorporate the impact of stochastic solar wind electron density variations in our analysis.

3.1 Methodology

To measure jitter noise for each pulsar in our sample, we first generated a frequency, time, and polarisation-averaged template and then calculated ToAs for frequency-averaged 8 s sub-integrations. These ToAs are then fitted using a previously derived timing model (Miles et al. Reference Miles2025a) with all parameters fixed except spin frequency to estimate the weighted-rms ( $W_\textrm{rms}$ ) of the timing residuals. We only fitted for spin frequency to absorb red noise-like signals so that the remaining scatter can be attributed to jitter noise. As we have used frequency-averaged observations, we did not fit for DM.

On short timescales, the error in ToA can be described as

(5) \begin{equation} \sigma_\textrm{tot}^{2} = \sigma_\textrm{S/N}^{2} + \sigma_{\textrm{J}}^2 + \sigma_\textrm{DISS}^{2} \end{equation}

where $\sigma_\textrm{S/N}$ , $\sigma_{\textrm{J}}$ , and $\sigma_\textrm{DISS}$ are the errors associated with the radiometer noise, jitter noise, and diffractive scintillation noise. The latter term arises from the stochasticity in the pulse broadening due to the turbulent IISM and is expected to contribute most to highly scattered pulsars that usually possess high DMs (Cordes & Shannon Reference Cordes and Shannon2010). Here we have neglected this because previous studies have suggested that for most pulsars the effect of $\sigma_\textrm{DISS}$ is found to be smaller than the $\sigma_{\textrm{J}}$ and $\sigma_\textrm{S/N}$ in a comparable sample (Wang et al. Reference Wang2024).

To measure jitter noise, we first simulated idealised ToAs with ToA uncertainties the same as the observations using TEMPO2 following the practice outlined in Parthasarathy et al. (2021), creating $1\,000$ realisations of the simulated dataset using the high flux density observation. We use simulations as they allow us to set robust limits and confidence intervals, which is especially important as we are setting jitter limits from short observations consisting of a few independent arrival times.Footnote ^a Therefore, uncertainty due to jitter noise in a sub-integration length $T_\textrm{sub}$ can be estimated to be

(6) \begin{equation} \sigma_{\mathrm{J}}^{2} (T_{\mathrm{sub}}) = \sigma_{\mathrm{obs}}^{2} (T_{\mathrm{sub}}) - \sigma_{\mathrm{sim}}^2 (T_{\mathrm{sub}}) \end{equation}

Similar to Parthasarathy et al. (2021), we report the detection of jitter noise if the variance in the observations is greater than that measured in 95% of the realisations of idealised data sets.

Jitter measurements are often reported as those that would result from one-hour duration integrations. Since jitter is expected to scale with the number of pulses being averaged as $1/\sqrt{N_p}$ , we can estimate jitter in an hour to be

(7) \begin{equation} \sigma_{\mathrm{J}} (\mathrm{hr}) = \frac{\sigma_{\mathrm{J}} (T_{\mathrm{sub}}) }{\sqrt{3\,600/T_{\mathrm{sub}}}}.\end{equation}

3.2 Results

We have measured jitter noise in 41 MSPs (20 of which are new) and report upper limits for the remaining 48 MSPs as presented in Table 2. For the brightest MSP, PSR J0437 $-$ 4715, we measured jitter noise of 52(7) ns in 1 hr. The lowest jitter level is shown by PSR J2241 $-$ 5236 at 3.5(3) ns in 1 h, consistent with the value obtained by Parthasarathy et al. (2021). Our MSP sample overlaps with previous MSP jitter studies such as (e.g. Shannon et al. Reference Shannon2014; Lam et al. Reference Lam2019; Parthasarathy et al. 2021; Wang et al. Reference Wang2024). A comparison of results is presented in Table 2. Our jitter measurements are consistent with the previous studies at an $\sim 1\sigma$ level.

Table 2. Measured jitter values and upper limits for 89 pulsars in the MPTA. For each pulsar we present the period of the pulsar (P), DM, weighted RMS of the timing residuals for the chosen observation, jitter noise in 1 hr $\sigma_{\textrm{J}}(\textrm{hr})$ , jitter noise for a single pulse $\sigma_{\textrm{J},1}$ and comparison of $\sigma_{\textrm{J}}(\textrm{hr})$ with previous studies. S14, L16, L19, P21, and W24 refer to Shannon et al. (Reference Shannon2014), Lam et al. (Reference Lam2016); Lam et al. (Reference Lam2019), Parthasarathy et al. (2021), and Wang et al. (Reference Wang2024), respectively.

Propagation of pulsar radio emission through the IISM causes flux density variations, that is, interstellar scintillation. As such, it is likely that jitter noise is the only limiting noise source for a fraction of our observations ( $F_\textrm{j}$ ), where the flux density is enhanced. We assume that an observation is jitter-limited if the ToA uncertainty ( $\sigma_\textrm{S/N}$ ), measured from a fully time and frequency-averaged observation, is less than the jitter noise in the epoch-wise observing length. In Figure 1, we plot $F_{j}$ for pulsars in which we have detected jitter with MeerKAT. We find that, for the nearest and brightest MSP, PSR J0437 $-$ 4715 (Shamohammadi et al. Reference Shamohammadi2024), all of the observations are jitter-limited, whereas for 15 pulsars (36%), none of the observations are jitter-limited.

Figure 1. Fraction of jitter-limited observations ( $F_{\textrm{J}}$ ) per pulsar in MPTA. Out of the 41 MSPs, 15 pulsars are not included in this analysis as none of the observations were jitter-limited according to our definition. The pulsars are sorted by the fraction, and PSR J0437 $-$ 4715 has 100% observations that are jitter-limited.

We have also searched for a correlation between pulse profile effective width, which is a proxy for the sharpness of the pulse profile $W_\textrm{eff}$ and jitter noise. We calculate $W_\textrm{eff}$ using pulse portraits, generated similar to those presented in Miles et al. (Reference Miles2023) using the PulsePortraiture method (Pennucci Reference Pennucci2019). We calculate $W_\textrm{eff}$ , using the relationship given in Cordes & Shannon (Reference Cordes and Shannon2010):

(8) \begin{equation} W_\textrm{eff} = \frac{\Delta t}{\sum_\textrm{j} \! (U_\textrm{j+1} - U_\textrm{j})^{2}}\end{equation}

where $U_\textrm{j}$ is the profile (normalised to peak flux density) and $\Delta t $ is the time resolution of the pulse profile. We modelled $\sigma_{\textrm{J}}$ for the MSPs using a linear relationship between $W_\textrm{eff}$ and $\sigma_{\textrm{J}}$ as shown in Figure 2. We derive a median and 2 $\sigma$ region dependence for a linear fit model using bilby (Ashton et al. Reference Ashton2019). In our model, we scaled the uncertainties by adding two parameters; one multiplied by the errors and the other added in quadrature in order to make the likelihood better describe the data, but we plot the data using the original uncertainties. We also considered a power law relationship between $W_\textrm{eff}$ and $\sigma_{\textrm{J}}$ (i.e. $\sigma_{\textrm{J}} \propto W_\textrm{eff}^\alpha)$ . We found $\alpha=0.95\pm 0.09 $ , which is consistent with a linear relationship.

Figure 2. Comparison of $W_\textrm{eff}$ and $\sigma_{\textrm{J}} (\textrm{hr})$ for the MPTA pulsars. The median relationship is $36.9(8.1) \times (W_\textrm{eff}/10^{-4} ) + 4.2(2.7)$ and the 2 $\sigma$ uncertainty regions for the likelihood fit are plotted.

Using our larger sample of pulsar jitter measurements, we searched for correlations between jitter noise and other pulsar intrinsic and extrinsic parameters, including pulse width ( $W_{50}$ ), sharpness of pulse profile ( $W_\textrm{eff}$ ), pulsar period (P), and flux density at 1 400 MHz ( $S_{1\,400}$ ) of the pulsar. We use values of $W_{50}$ reported in the ATNF pulsar catalogue (Manchester et al. Reference Manchester, Hobbs, Teoh and Hobbs2005). To be consistent with previous studies, we also converted $\sigma_J(\textrm{hr})$ to single pulse jitter $\sigma_{\textrm{J}, 1}$ , (replacing $T_\textrm{sub}$ as P in equation 7) and the jitter parameter $k_{\textrm{J}} = \sigma_{\textrm{J}, 1}/P$ , which removes bias related to pulse period. We use the Spearman correlation coefficient (R) as our test statistic and use its p-value to test for statistical significance. We summarise our findings and compare them with previous works in Table 3. Our results suggest that there is significant correlation between pulse jitter and the $W_\textrm{eff}$ of a pulsar. However, when normalising by pulse period, the correlation becomes insignificant.

Table 3. Spearman correlation coefficient between jitter noise and pulsar parameters. The first four columns denote jitter noise, pulsar parameters, correlation coefficient and associated p-value. The last column shows correlation coefficients reported in previous works and associated p-value if available. Here L+19, P+21 and W+24 refers to Lam et al. (Reference Lam2019), Parthasarathy et al. (2021) and Wang et al. (Reference Wang2024) respectively. We note that both the $W_\textrm{eff}$ and R relation used by Lam et al. (Reference Lam2019) differed from other works.

We also investigated the connection between jitter noise derived from our single-epoch observations and that derived from long-term timing observations. In the long-term pulsar noise analysis, ECORR describes the uncertainties that are correlated between contemporaneous arrival times for each observing backend (NANOGrav Collaboration et al. 2015). This term includes pulse jitter noise but could also include other noise terms correlated within an observation but uncorrelated between epochs, such as broadband RFI or IISM effects (Coles et al. Reference Coles, Rickett, Gao, Hobbs and Verbiest2010). In Figure 3, we compare the ECORR values estimated by Miles et al. (Reference Miles2025a) to the jitter noise values estimated in this work. While for some pulsars, there is the expected linear correlation between jitter noise and ECORR, for the majority, it is clear that ECORR overestimates jitter noise. Previously, Lam et al. (Reference Lam2016) reported a similar overestimation of jitter noise in the NANOGrav 9-yr dataset. If the source of excess noise measured in ECORR can be mitigated, then using the measured jitter noise values instead of ECORR will reduce the number of free parameters and could improve the sensitivity of PTAs. To demonstrate this, we compare the sensitivity of the MPTA assuming the ECORR values derived from long-term timing against the sensitivity assuming ECORR corresponding to $\sigma_{\textrm{J}}$ . We see an improvement of 8% when using jitter noise as the ECORR value for a GWB amplitude of $1 \times 10^{-16}$ , the scenario that would show the greatest improvement. For the rest of the work, we will use ECORR derived from long-term timing in our array sensitivity calculations.

Figure 3. Comparison of ECORR parameters derived in MPTA noise analysis and jitter values estimated in this work. The blue points and green points denote direct measurements of jitter noise and upper limits, respectively. Jitter measurements are systematically lower than the ECORR values. The red dashed line corresponds to a 1:1 relationship between $\sigma_{\textrm{J}} (\textrm{hr})$ and ECORR.

Jitter noise will have a larger contribution for MSPs observed with even more sensitive telescopes, such as the SKA-Mid telescope. The design sensitivity for the SKA-Mid (AA4, comprising 197 antennas) is forecast to be a factor of $4.1$ more sensitive than the MeerKAT telescope (Braun et al. Reference Braun, Bonaldi, Bourke, Keane and Wagg2019), and it will have an L-band observing system with comparable bandwidth to MeerKAT. We simulated SKA-Mid observations for our MSP sample by scaling the ToA uncertainties by their sensitivities and repeated the above procedure to calculate $F_{j}$ . For the MSPs without a $\sigma_{\textrm{J}}$ measurement, we make a prediction based on the correlation between $W_\textrm{eff}$ and $\sigma_J$ described above. Based on these relations, we simulated the SKA-Mid $\sigma_{\textrm{J}}$ using the median values for each MSP, and pessimistic (jitter noise $+ 2\sigma$ above the median) and optimistic ( $2\sigma$ below the median) values. The fraction of observations that would be jitter-limited with the SKA-Mid is displayed in Figure 4, with pulsars sorted by the ‘median’ model. In the pessimistic case, we can see spikes in $F_\textrm{j}$ for certain pulsars. We note that for these pulsars, ToA uncertainties are more Gaussian distributed than the other pulsars. Therefore, when a certain scenario is assumed, for example, pessimistic, a higher number of observations becomes jitter-limited compared to the remaining pulsars. Owing to the increase in sensitivity, more pulsars are completely jitter-limited with SKA-Mid, including one of the most precisely timed MSPs, PSR J1909 $-$ 3744.

Figure 4. Fraction of jitter-limited observations ( $F_{\textrm{J}}$ ) for each pulsar with a telescope with SKA-Mid sensitivity. Jitter noise is measured in this work for 41 MSPs. For the MSPs without jitter measurements, we have used the relation between $W_\textrm{eff}$ and $\sigma_{\textrm{J}}(\textrm{hr})$ derived from Figure 2 to estimate the jitter. Hence, three fractions are shown: one using the median relation, and the other two showing the 2 $\sigma$ region to define optimistic (green) and pessimistic (blue) scenarios for the jitter noise in these pulsars. The fractional values are sorted by the median method. Many high-precision MSPs in SKA-Mid will be completely jitter-limited if observed with the entire array.

Finally, we considered whether it is possible to improve array sensitivity by increasing the integration time for the jitter-limited pulsars and reducing the time of other pulsars such that the total observing hours remain the same. We created a jitter-based sub-array by choosing the 11 most highly jitter-limited pulsars with the MPTA, as shown in Figure 1. In particular, we considered a scenario where these pulsars were observed 10 hrs in total per observing session using 16 dishes of MeerKAT at L-band, with the remaining pulsars observed using 48 dishes. Considering a typical observation length of 256 s, this sub-array method will increase the ToA uncertainty in a jitter-limited pulsar by a factor of 1.11 while decreasing the jitter noise by a factor of 3.6. For such a jitter-based sub-array, we did not find an improvement in sensitivity relative to the standard strategy. If we compare the change in sensitivity of an array just made from the 11 jitter-limited pulsars, then observing longer increases sensitivity by 20%, but is counterbalanced by the loss in sensitivity (15%) for the remaining pulsars due to the reduction in the number of dishes and therefore telescope gain.

4.1 DM noise modelling

As the pulsar radiation traverses different lines of sight due to the motion of the pulsar, the Earth, and the IISM, it encounters different electron column densities. As a result, the DM varies with time, introducing time and frequency-correlated delays proportional to $\nu^{-2}$ , where $\nu$ is the radio frequency at which ToAs are calculated. There are two common ways to mitigate and fit for the effects of DM variations. Many PTA data analyses model the DM variations as a stochastic correlated GP with a red noise power spectral density given by equation (3). The DM GP power law amplitude and spectral index are inferred from the timing residuals with a series of Fourier basis functions. An alternative approach, in particular used by the NANOGrav collaboration, is to use a piece-wise constant function to model DM variations (Demorest et al. Reference Demorest2013). Owing to their multi-epoch observation with different frequency band receivers, they assume that DM does not change over a period of $\sim$ 14 days (although this number is tuned for each pulsar) and fit for the DM variations relative to a fiducial value over their observation time span. This method is commonly referred to as the DMX method. Here, we analyse our dataset based on both these DM variation models to assess which model provides higher sensitivity to GWs.

To compare the DM GP and DMX methods, we use the 4.5-yr MPTA data set as described above. To account for DMX, we needed to choose a DM sampling cadence, which is included in the pulsar timing ephemerides. We choose DM epochs centred around each observation of a pulsar. This is appropriate as all MPTA observations are taken in the same band, so all observations have the same frequency coverage. Using the standard timing software TEMPO2, a unique DM offset is fitted for all the ToAs in each epoch-based window using the $\nu^{-2}$ dependence while the global DM is fixed.

For all of the MSPs in the MPTA, we studied the effect of the choice of DM noise modelling methodology on the sensitivity of the array for a range of DM GP properties. For individual MSPs, we varied the DM noise amplitude (in the range $-20 \lt \log_{10} A_\textrm{ DM} \lt -10$ ) and spectral index around the value ( $\gamma_\textrm{DM} \pm 2$ ) inferred from the 4.5-yr noise analysis (Miles et al. Reference Miles2025a), and estimated the minimum SGWB amplitude that can be detected for each combination. The motivation behind this was to determine if pulsars with lower dispersion measure were more sensitive to the choice of DM modelling approach than the higher dispersion measure pulsars (or vice versa). We found that the GP method is more sensitive to an SGWB in all cases. We then compared the entire array sensitivity with respect to the choice of DM noise model. The noise sensitivity curve for both cases is shown in Figure 5. We found that the DM GP method improves the array sensitivity by 44% compared to the DMX method. We attribute this decrease in sensitivity to the large number of additional parameters added to the timing model when DMX is used. This was also observed for an analysis of a similar DM correction scheme described in Keith et al. (Reference Keith2013). However, we note that there are many caveats to choosing one DM model over the other. In the DM GP model, DM variations as a function of time, DM(t), are modelled as a finite sum of sine and cosine functions, which are truncated to lower computational costs. Therefore, the recovered amplitude and spectral index are approximations of the true solutions (van Haasteren & Vallisneri Reference van Haasteren and Vallisneri2014). Apart from that, the DM GP model assumes stationarity of the signal, whereas extreme scattering events can make this assumption invalid (Coles et al. Reference Coles2015; Kerr et al. Reference Kerr2018). The DMX model can potentially better account for this non-stationarity, whereas the DM GP model better models the long-term correlated noise due to DM.

Figure 5. Comparison of sensitivity curves for the MPTA with alternate DM noise modelling methods. The cyan curve corresponds to white noise, red noise, and chromatic noise, the purple curve additionally includes DM as a GP, whereas the crimson curve includes DMX. The DMX implementation reduces the sensitivity by 44% compared with the DM GP model.

4.2 Observing strategies with the multi-frequency receivers of MeerKAT

The flexibility of MeerKAT motivates optimising sensitivity through the choice of receiver or sub-array. In order to assess the sensitivity, it is first necessary to estimate the achievable timing precision in other observing bands. Due to the lack of timing campaigns with the UHF and S-band receivers, it was necessary to base timing precisions on simulations extrapolated from the L-band observations. Using the 16 frequency channel observations, we first created idealised ToAs with added Gaussian noise using TEMPO2. This way, the ToA errors were consistent with the original observation. In our simulations, we assumed the entire bandwidths at the UHF band and S band are free from RFI. The 544 MHz and $\sim$ 870 MHz of available bandwidth of UHF and S-band systems were divided into 10 and 16 channels, respectively, to match the channel bandwidth of the L-band system. To simulate the ToA error ( $\sigma_\textrm{ToA}$ ) for the different frequency receivers, we started with an assumed relationship between timing precision and effective width ( $W_\textrm{eff}$ ):

(9) \begin{equation} \sigma_\textrm{ToA} \propto \frac{W_\mathrm{eff}}{\textrm{S/N}} (T_\textrm{sys} + T_\textrm{sky}), \end{equation}

where $T_\textrm{sys}$ and $T_\textrm{sky}$ are the system and sky temperature, respectively. The $W_\textrm{eff}$ , S/N of observations, $T_\textrm{sys}$ and $T_\textrm{sky}$ change with frequency, so this relation can be used to extrapolate the ToA uncertainties at different frequencies. Here, we have considered $T_\textrm{sys}$ to be fixed at 18 K for L-band and 23 K for UHF band, which correspond to 1 400 and 815 MHz (centre frequency of UHF-band), respectively. For each pulsar, we calculated $T_\textrm{sky}$ based on the Haslam et al. (Reference Haslam, Salter, Stoffel and Wilson1982) 408 MHz all sky brightness temperature map for each pulsar position, and scaled it as $\nu^{-2.6}$ (Lawson et al. Reference Lawson, Mayer, Osborne and Parkinson1987) at each frequency sub-band, where $\nu$ is the radio frequency. For the S band, $T_\textrm{sky}$ will have negligible effects on sensitivity and hence can be neglected. Since the flux density of a pulsar is proportional to the observed S/N, we used the known measured spectral indices of the pulsars, reported in Gitika et al. (Reference Gitika2023), to extrapolate flux density values at different frequencies. The pulse width of a pulsar is expected to increase because of scattering and intrinsic pulse-width variations (Sieber, Reinecke, & Wielebinski Reference Sieber, Reinecke and Wielebinski1975; Slee, Bobra, & Alurkar Reference Slee, Bobra and Alurkar1987). However, Miles et al. (Reference Miles2023) showed that $W_\textrm{eff}$ can exhibit diverse and sometimes non-monotonic evolution with pulse frequency for MPTA pulsars. Using measurements of $W_\textrm{eff}$ as a function of frequency in L band, revealed that at UHF frequencies, some pulsars can have three times higher or lower $W_\textrm{eff}$ in comparison to L band. To model the frequency dependence of $W_\textrm{eff}$ , we considered power law, quadratic, and cubic polynomial relationships and selected the model that best described the L-band observations as measured by the lowest chi-squared value. Of the 83 MSPs, 62 pulsars followed a quadratic relation between $W_\textrm{eff}$ and frequency with the remaining having either power law or cubic relationships.

For each epoch, the ToA uncertainties were estimated using equation (9) and the width and flux density relations with frequency as derived above. As a result, the ToA uncertainty remains constant across all the epochs but changes across frequency. We note that this neglects the effect of scintillation modulating the S/N of the observations, and hence the effects of variable timing precision. To make the comparisons between bands fairer, we also simulated L-band observations based on these extrapolations. Between the observed L-band dataset and the simulated dataset, the sensitivity differs by just 5%. This suggests that incorporating scintillation into our simulated data set would only marginally change the results.

4.3 Confirming our UHF predictions

While we did not have sufficient observations to measure timing precision in the UHF-band, to assess if our projections were reasonable, we conducted a single epoch of observations of the MPTA pulsars in March 2024 with MeerKAT using the UHF receiving system. We investigated whether our predictions for $W_\textrm{eff}$ and flux density matched the observations. The tests were impacted by only having a single epoch. For low DM pulsars, diffractive scintillation introduces an additional systematic in the uncertainties of the flux density measurements and can make it difficult to assess how the pulse profile is evolving with frequency. For similar reasons, it was also difficult to make reliable measurements for the weaker pulsars in the MPTA sample. As a result, template creation with the PulsePortraiture method did not always converge. Instead, to estimate $W_\textrm{eff}$ for the UHF observations, we averaged the observations to 10 frequency channels and smoothed the profiles. For many low DM pulsars, we found a deviation in the projected flux density and frequency power law due to scatter broadening of the pulse profile at the lowest frequencies. In most cases, our predictions for $W_\textrm{eff}$ as a function of frequency were consistent with observations. However, for 25 MSPs, the observation followed either a broken power law or a simple power law with an amplitude and index different from those derived using L-band observations. We modified our models and simulated the dataset based on these new observations.

4.4 Sub-arrays with MeerKAT

4.4.1 Observing strategies for SGWB detection

Using the L-band, UHF-band, and S-band simulated data sets, we considered the following observing strategies and compared the resulting sensitivity:

1. 1. Observe with L-band exclusively

2. 2. Observe with UHF-band exclusively

3. 3. Observe with S-band exclusively

4. 4. Observe with UHF and L-band simultaneously with 32-antenna sub-arrays

5. 5. Observe with UHF and L-band alternatively, maintaining the same overall cadence

6. 6. Observe with the optimal strategy out of cases 1–5 for each pulsar

For case (4), we combined the simulated L-band and UHF-band data sets and multiplied the ToA uncertainties by a factor of two to account for the loss in gain due to the reduction in the number of antennas in each sub-array. For case (5), the dataset was created by discarding alternate epochs of observations from the combined L and UHF-band datasets. For case (6), we chose the best method of observation for each pulsar out of the remaining observing strategies.

For each strategy, we use hasasia to generate the sensitivity curves for individual MSPs assuming a low ( $1 \times 10^{-16}$ ), medium ( $2 \times 10^{-15}$ ) and high ( $5 \times 10^{-15}$ ) amplitude for the SGWB. Each curve assumes DM GP and ECORR model for which we have used MPTA 4.5 yr noise models given in Miles et al. (Reference Miles2025a). For the majority of MSPs, UHF-only observations are the preferred strategy, with case (5) being the second most preferred. S-band observations were the most sensitive for 13 MSPs. Many of these pulsars have significant scattering noise contributions, which dominate at lower frequencies and decrease the timing precision at low frequencies. Figure 6 shows the sensitivity curves for the entire array based on the above six observing strategies. Out of the six observing strategies, UHF exclusive observations (case 2) provide 30% improvement over L-band exclusive observations. S-band exclusive observations are a factor of $1.8$ less sensitive than L-band observations for an SGWB amplitude of $1 \times 10^{-16}$ . This is expected given the low flux density observations at the S-band and steep spectral indices. The optimal strategy is nearly identical to the UHF exclusive strategy. To understand if we could mitigate frequency-dependent noise processes in UHF and L band equivalently, we calculated the error in infinite frequency ToAs using the equation in Cordes & Shannon (Reference Cordes and Shannon2010). For the simulated dataset of PSR J0437–4715, we found the error in infinite frequency ToA in UHF-band to be 10% better than L-band. We conclude that trade-offs between timing precision (impacted by sensitivity and pulse width) and frequency coverage result in similar ability to correct for DM in UHF and L-band.

Figure 6. Sensitivity curves for MPTA observing strategies. The curves show the following strategies: L-band (L), UHF-band (UHF), S-band (S), observing with L and UHF simultaneously (L+UHF), observing with L and UHF-band alternatively (L+UHF alt) and choosing the optimal method out of the above for each pulsar (OPT). The three sub-panels show the same curves for three values of SGWB self-noise: $1 \times 10^{-16}, 2 \times 10^{-15}, 5 \times 10^{-15} $ .

Among the mixed L-band and UHF-band observations, alternate observations show higher sensitivity than simultaneous sub-array observations using half of the dishes in each array due to favourable linear ToA uncertainty scaling with gain in contrast to the square root improvement that comes from increasing observing time. This suggests that the increased ToA precision is more important than simultaneous observations, assuming propagation effects can be optimally modelled, an assumption we investigate further below. We also note that at higher GW frequencies, the slope of the sensitivity curve for case (4) is higher than for case (5) as the use of fewer antennas increases the white noise contribution. However, for practical purposes, simultaneous observations might be beneficial in characterising low-frequency effects like chromatic noise, especially if they vary on the timescale of the observing cadence.

As we increase the SGWB self-noise, the lower frequencies become GW-dominated, and the amount of improvement decreases between the UHF-only and L-band only observations. For SGWB amplitude of $2 \times 10^{-15}$ and $5 \times 10^{-15}$ , UHF is only 8% and 2% more sensitive than L-band observation, respectively. We also note these improvements assume perfect correction of chromatic effects, an assumption we test below.

4.4.2 Observing strategies for single source detection

Continuous GWs from non-evolving single SMBHBs could dominate over a stochastic background, particularly at higher frequencies ( $f \gtrsim$ a few 10 s of nHz) (Sesana, Vecchio, & Colacino Reference Sesana, Vecchio and Colacino2008; Sesana, Vecchio, & Volonteri Reference Sesana, Vecchio and Volonteri2009). We investigated whether alternate observing strategies in further observations would result in improved sensitivity by estimating the sky-averaged sensitivity of the array with five years of observations with the current MPTA observing cadence (as the MPTA has had a consistent observing strategy for five years up to the end of 2024) and then adjusting the observing strategy for the subsequent three years. We chose to modify the cadence for three years, as that marks the beginning of science verification for 64 SKA-mid dishes, which at the time of writing this paper, was expected to be combined with the 64 MeerKAT antennas in late 2028.

The strategies we considered with simulated datasets were, Case (a) 5 yr MPTA observations + 3 yr observations with the current MPTA strategy

Case (b) 5 yr MPTA observations + 3 yr observations of 20 best MSPs with a factor of four increased cadence while reducing the observation cadence of remaining MSPs, keeping the telescope time fixed

Case (c) 5 yr MPTA observations + 3 yr observations of 20 best pulsars with a quadrupled cadence

Case (d) 5 yr MPTA observations + 3 yr observations of 5 added MSPs with 4 us ToA error without any red noise.

In Figure 7, we plotted the sensitivity curves for the single source strategies. The high-cadence observations of the best pulsars, that is, case (b) provide the highest sensitivity at the higher frequencies, with a $\sim$ 14% improvement over case (a). Timing the best MSPs (case c) provides similar sensitivity towards single source detection as timing the entire array. The addition of five pulsars (case d) marginally improves the sensitivity compared to case (a), therefore, we did not include the corresponding curve in the Figure 7. We speculate these pulsars are not strongly contributing to the array’s sensitivity because of their higher ToA errors and shorter observing spans. For a more realistic scenario, when we increase the SGWB self-noise, the sensitivity to single sources remains consistent across all observing strategies.

Figure 7. Sky-averaged sensitivity of MPTA for a single continuous GW source. We consider cases where the SGWB self-noise amplitude is $1 \times 10^{-16}$ (solid), $2 \times 10^{-15}$ (dashed), and $5 \times 10^{-15}$ (dotted). We compare the sensitivities for different observing strategies: (a) L-band simulated 8 yr dataset (green) (b) L-band 5-yr dataset + 3-yr dataset with high-cadence observations of the 20 best pulsars, and low-cadence observations of the remaining pulsars (black) (c) L-band 5 yr dataset of 20 best pulsars + 3 year dataset of high-cadence observations (red). Case (a) and (b) provide similar sensitivities for SGWB and single sources in the case of higher GW self-noise, whereas case (b) provides best sensitivity for single source detection without reducing sensitivity towards an SGWB.

As the MPTA aims to optimise for both single source and SGWB detection, we need to understand how the choice of single source strategy affects SGWB sensitivity. For the same observing strategies, we also compared the SGWB sensitivity. At lower frequencies, case (b) and case (c) provide nearly identical sensitivities for low GW self-noise. Case (a) and case (b) also provide nearly identical sensitivities for GW self-noise $2 \times 10^{-15}$ and $5 \times 10^{-15}$ . We found that case (c) is less sensitive as we increase the SGWB self-noise, while the other strategies decrease the sensitivity to an SGWB by only 2–3%. This could be because the inter-pulsar correlations contribute significantly to the sensitivity of an SGWB rather than the high-cadence and high precision of the best pulsars. Therefore, out of the different observing strategies, case (b) optimises the sensitivity to single sources with minimal reduction in sensitivity to an SGWB while keeping the overall observing time constant. We note that we have not considered the ability of the array to angularly resolve a source and have used a sky-averaged approximation only. Angular resolution is improved with timing a larger number of pulsars, so there is likely a trade-off between timing fewer pulsars to higher precision and being able to localise an individual source. We defer this analysis to future work.

4.5 Challenges of timing with UHF-band

The MPTA has observed MSPs exclusively using the L-band receiver. In contrast, most of the PTAs use more than one telescope and therefore, different backends, multiple receivers, or both, which makes it important for them to address the systematic differences during data combination. In our study, we showed that observing pulsars with the UHF receiver of MeerKAT could improve the sensitivity of the MPTA more significantly than other chosen observing strategies, under ideal assumptions. While our assumptions include a variety of chromatic noise processes (both DM noise and chromatic noise potentially attributed to scattering), our modelling assumes they can be perfectly corrected. It is possible that the DM (and hence DM variations) of a pulsar are frequency-dependent as the ray bundles at different frequencies trace different lines of sight in the IISM (Cordes et al. Reference Cordes, Shannon and Stinebring2016; Lam et al. Reference Lam2020). This has been predicted to cause the dispersion measure to vary at different frequencies and has been tentatively detected in the PSR J2219+4754 and PSR J2241 $-$ 5236 (Donner et al. Reference Donner2019; Kaur et al. Reference Kaur2022). In a similar effect, multi-path propagation through non-uniform small-scale density fluctuations in the IISM causes the lower frequency emission to be highly scattered, resulting in a decorrelation in chromatic noise between different frequency components (Shannon & Cordes Reference Shannon and Cordes2017). If such decorrelated scattering noise is present in our observations, it could reduce the sensitivity of low-frequency observations. In many PTA experiments, band noise has been identified in pulsar timing data sets, modelled as additional noise unique to each observing frequency sub-band (Lazarus et al. Reference Lazarus2016; Reardon et al. Reference Reardon2023b). This could be an example of such propagation noise.

Figure 8. Sensitivity curves for individual pulsars with decorrelated chromatic noise. Top panel: For PSR J1017–7156, UHF is less sensitive than the L-band with added chromatic noise at 10% decorrelation. Bottom panel: In the case of PSR J2129–5721, UHF is less sensitive than the L-band with added chromatic noise at 50% decorrelation. The solid lines show the sensitivity curves using the fully correlated chromatic noise.

Figure 9. Comparison of MPTA 4.5-yr sensitivity with the other PTAs for an SGWB at an amplitude $2 \times 10^{-15}$ . The shaded region shows the 68% and 99.7% confidence region of SGWB amplitudes derived from merger rates from redshift-dependent galaxy mass functions, fraction of close galaxy pairs and overmassive BHs in galaxies. The top-left panel includes all models described in Sesana (Reference Sesana2013), the top-right panel includes only the best estimates for models, the bottom-left panel includes models accounting for observations of overmassive BHs and the bottom-right panel includes models allowing for a broken power law relation between stellar velocity dispersion and bulge stellar mass.

To assess the impact of such noise on the sensitivity of the MPTA, we developed a new model that includes partially correlated chromatic noise, similar to band noise modelling in Lam et al. (Reference Lam2016), Goncharov et al. (Reference Goncharov2021), Reardon et al. (Reference Reardon2023b). In this alternate chromatic noise model, one component of the noise is fully correlated in frequency and time, while another component is uncorrelated between different sub-bands (but has the same temporal correlation). The partially correlated noise was added as a power law, with power spectral density following equation (4) for ToAs, selected by frequency sub-band. Since our aim was to determine what combination of partially correlated noise would alter our UHF sensitivity projections given in Section 4.4, we modelled this additional noise assuming a range of fractional correlations. We varied the level of correlation by adjusting the amplitude of the power spectral density of the correlated components to assess the potential impact of such an effect. When we decrease the degree of correlation, the UHF-only observing strategy becomes less sensitive than the L-band. The level depends on the choice of the pulsar. Two of the most extreme cases are shown in Figure 8. The dashed curves in Figure 8 show the final L-band and UHF-band sensitivities for individual pulsars at the degree of correlation where the UHF-band becomes less sensitive than the L-band. For PSR J2129 $-$ 5721, the UHF-band becomes less sensitive than the L-band if the fraction of chromatic noise that is correlated is 50% whereas for many pulsars, even smaller levels (10%) of correlated chromatic noise decrease UHF sensitivity, which is the case for PSR J1017 $-$ 7156. This new chromatic noise model is important for understanding the effects of decorrelated chromatic noise and gives a perspective on the diminishing low-frequency effects. Using dedicated observing campaigns at low frequencies, we can determine if such processes are present in PTA data sets.