The MeerTime Pulsar Timing Array -- A Census of Emission Properties and Timing Potential

MeerTime is a five-year Large Survey Project to time pulsars with MeerKAT, the 64-dish South African precursor to the Square Kilometre Array. The science goals for the programme include timing millisecond pulsars (MSPs) to high precision (<1 $\mu$s) to study the Galactic MSP population and to contribute to global efforts to detect nanohertz gravitational waves with the International Pulsar Timing Array (IPTA). In order to plan for the remainder of the programme and to use the allocated time most efficiently, we have conducted an initial census with the MeerKAT"L-band"receiver of 189 MSPs visible to MeerKAT and here present their dispersion measures, polarization profiles, polarization fractions, rotation measures, flux density measurements, spectral indices, and timing potential. As all of these observations are taken with the same instrument (which uses coherent dedispersion, interferometric polarization calibration techniques, and a uniform flux scale), they present an excellent resource for population studies. We used wideband pulse portraits as timing standards for each MSP and demonstrated that the MeerTime Pulsar Timing Array (MPTA) can already contribute significantly to the IPTA as it currently achieves better than 1 $\mu$s timing accuracy on 89 MSPs (observed with fortnightly cadence). By the conclusion of the initial five-year MeerTime programme in July 2024, the MPTA will be extremely significant in global efforts to detect the gravitational wave background with a contribution to the detection statistic comparable to other long-standing timing programmes.


INTRODUCTION
Millisecond pulsars (MSPs) are the sub-population of pulsars that are characterised by short spin periods ( < 50 ms) and low magnetic field strengths ( surf < 10 10 G) . The first binary pulsar, PSR B1913+16 (J1915+1606; Hulse & Taylor 1975), * E-mail: renee.spiewak@manchester.ac.uk (RS) This definition of an MSP is relatively simple (c.f. Lee et al., 2012) but useful for this work and in relation to the other MeerTime themes described in section 1. possessed an unusually low spin period ( = 59 ms) and magnetic field strength ( surf = 2 × 10 10 G) compared to the bulk of the population known at the time, which, apart from the young Crab and Vela pulsars, consisted mainly of solitary pulsars with spin periods of typically 0.1-2 s and magnetic fields of order 10 12 G. Prior to the discovery of PSR B1913+16, Bisnovatyi-Kogan & Komberg (1974) had predicted that mass accretion onto neutron stars (NSs) may decrease their magnetic fields, and PSR B1913+16 appeared to agree with 2 R. Spiewak et al. their prediction as at least one of the NSs in this system must have accreted matter during a previous stage of mass transfer. When the first MSP, PSR B1937+21 (J1939+2134), was discovered with = 1.56 ms and surf = 4 × 10 8 G (Backer et al., 1982), it was soon proposed that a (now missing) companion had been responsible for its spin-up and the low field made spin-up to its millisecond period possible (Alpar et al., 1982;Radhakrishnan & Srinivasan, 1982). As more recycled pulsars were discovered, a self-consistent model began to emerge (e.g., Bailes, 1989) in which most pulsars lost their binary companions in the initial explosion due to impulsive kicks and mass loss at the time of the first explosion, but those that remained bound had a chance for a second life upon accreting material from their evolved companions that not only reduced their magnetic field strengths but made it possible for them to be spun-up to millisecond periods. NSs with lower-mass companions accrete material for longer than those with high-mass companions because their evolutionary timescales are longer, and they are able to be spun up to very short periods ( ∼ 1.5 ms) with an associated reduction of their magnetic fields to ∼ 10 8 G. MSPs with heavier white dwarf companions tend to spin more slowly ( ∼ 10 -16 ms) with some notable exceptions, such as PSR J1614−2230 (Crawford et al., 2006;Demorest et al., 2010), and those with NS companions slower still ( ∼ 16 -60 ms). For an early paper describing these models, see Bhattacharya & van den Heuvel (1991), and, for a more recent discussion of their binary properties, Tauris et al. (2017). Support for these models (see, e.g., Phinney, 1992) comes in part from the large fraction of MSPs observed to have binary companions compared to the fraction of "normal" pulsars in binaries (∼ 74% of MSPs are in binaries but only ∼ 2% of longer-period pulsars are in binaries; from the Australia Telescope National Facility (ATNF) Pulsar Catalogue, v.1.64, Manchester et al. 2005).
Studies of the properties of MSPs and the Galactic MSP population overall have often been limited or biased by small sample sizes (e.g., Johnston & Bailes, 1991) or by instrumentation that prohibits studies of the polarization or novel features of the pulse profile, which are probably determined by the pulsar's magnetosphere and spin period. See, for example, Kramer et al. (1998), Xilouris et al. (1998), and Kramer et al. (1999) for a comprehensive study of 27 MSPs, or the review by Kramer & Xilouris (2000), and references therein. The majority of MSPs have been discovered by quasi-uniform surveys by large-aperture telescopes such as the Parkes 64-m radio telescope (also known as Murriyang; e.g., Manchester et al. 2001;Keith et al. 2010), the Arecibo 305-m dish (e.g., Cordes et al., 2006), and the Green Bank 100-m Telescope (GBT; e.g., Boyles et al. 2013;Stovall et al. 2014). More recently, other instruments like the Low-Frequency Array (LOFAR) and Giant Metrewave Radio Telescope (GMRT) have discovered pulsars in a combination of quasi-uniform surveys and targeted observations of gamma-ray point sources with pulsar characteristics (see, e.g., Ray et al., 2012;Bhattacharyya et al., 2016;Sanidas et al., 2019). In many cases, the only published pulse profiles are those formed using survey data and thus limited in resolution due to survey data rates as well as the use of incoherent dispersion correction. It is not uncommon for pulsars with large dispersion measures (DMs) to have narrow features smeared beyond recognition using the discovery hardware which cannot a priori know the DM of the pulsars.
There is also some inconsistency in reported MSP flux densities in the ATNF Pulsar Catalogue and the literature. This is for a variety of reasons that include uncertainties in the pulsar position within the primary beam during initial timing, pulsar scintillation that makes pulsars more likely to have been brighter in their discovery observations, and design limitations that make absolute flux calibration a secondary consideration in the construction of some pulsar survey hardware. For example, the flux density of PSR J1843−1113 was first determined by Hobbs et al. (2004) to be 1400 = 0.1 mJy, but was revised to 0.3 mJy by Levin et al. (2013), and then to 0.6 mJy by Desvignes et al. (2016).
Analyses of MSP polarization profiles and their evolution with radio frequency, which could potentially lead to a more consistent theory of the pulsar emission mechanism, have only been performed on samples of a few to ∼ 20 MSPs (e.g., Rankin et al., 2017, and references therein), and the field could significantly benefit from a large comprehensive analysis derived from high signal-to-noise ratio ( / ) observations. Additional areas of pulsar astronomy that would benefit from such a dataset include probing the Galactic magnetic field with rotation measure (RM) measurements (e.g., Han et al., 2018;Sobey et al., 2019), improving Galactic electron density models with new independent distance measurements (see, e.g., Matthews et al., 2016;Yao et al., 2017), and constraining models of the NS equation of state through MSP mass measurements (Özel & Freire, 2016).
Estimates of the total MSP (Levin et al., 2013) and binary NS-NS populations (e.g., Burgay et al., 2003) rely upon our understanding of the MSP luminosity function, which in turn depends upon measured MSP distances, flux densities and their pulse shapes. Hence an authoritative survey of MSPs with the same flux density scale, coherent dedispersion to ascertain the true pulse shape, and ultimately pulsar parallaxes from timing is very valuable when we want to estimate their underlying populations and ascertain potential yields for future pulsar surveys and timing programmes such as those planned for the Square Kilometre Array (SKA; e.g., Levin et al., 2018).
The new precursor to the SKA in South Africa, MeerKAT (Meer Karoo Array Telescope), presents an opportunity to revisit the fluxes, spectral indices, pulse shapes, polarimetry, and timing potential of the MSPs visible to it. The MeerTime project (Bailes et al., 2020) is a MeerKAT Large Survey Project and has four major themes. These include the Thousand Pulsar Array (Johnston et al., 2020), the Relativistic Binary programme (Kramer et al., 2021), the Globular Cluster pulsar timing programme (Ridolfi et al., 2021), and the MeerTime Pulsar Timing Array (MPTA), which studies MSPs with the primary aim of detecting gravitational waves (GWs). Parthasarathy et al. (2021) analysed the pulse-to-pulse vari-ability ("jitter") of 29 of these MSPs, placing limits on their timing precision. In this work, we present a census of 189 MSPs visible to MeerKAT at ∼ 1400 MHz frequencies, providing high-time-resolution, polarization-and flux-calibrated profiles, measured DMs, polarization fractions, RMs, subbanded flux densities, and spectral indices, as well as timing potential for the MPTA. The first results of the MPTA will be presented by Miles et al. (2022, in prep.), including a full release of timing residuals and refined ephemerides, and closely followed by an analysis of red noise and GW search (Middleton et al., 2022, in prep.. The 64 13.5-m offset Gregorian MeerKAT dishes have a high aperture efficiency and low system temperature that achieves a system equivalent flux density of just ≈ 7.0 Jy at ∼ 1400 MHz frequencies with the full array (Jonas & MeerKAT Team, 2016). With an antenna gain roughly 4 times that of Parkes, a lower system temperature (∼ 18 K), and a larger bandwidth (856 MHz of the MeerKAT L-band receiver compared with the 340 MHz of the Parkes multibeam receiver that was used for much of the last decade for the Parkes Pulsar Timing Array (PPTA; Manchester et al., 2013)), MeerKAT has become the best telescope in the southern hemisphere for a large-scale MSP observing programme and population study (Bailes et al., 2020). Compared to, for example, the Fivehundred-meter Aperture Synthesis Telescope (FAST), which can observe down to Dec > −15 degrees, MeerKAT has a lower sensitivity but greater sky coverage (Dec < +40 degrees, which includes the MSP-rich Galactic centre) and significantly higher slew rates (60 and 120 deg min −1 in elevation and azimuth respectively). The 500-m diameter FAST has much superior gain but can only observe about one source every ten minutes due to reconfiguration overheads. The ability to quickly move between sources with minimal overhead is one of the advantages of large-N (number), small-D (diameter) arrays like MeerKAT.
The timing precision and stability of MSPs can additionally be leveraged to search for low-frequency GWs through observations of arrays of MSPs like a galactic-scale interferometer, termed Pulsar Timing Arrays (PTAs; Hellings & Downs, 1983). The current global effort, the International Pulsar Timing Array (IPTA), places stringent constraints on the amplitude of a GW background but no detection has been made (Perera et al., 2019), although there are tantalising indications of similar spectral slopes in the red noise of the timing residuals of many of the MSPs with no spatial correlations indicative of a GW background (e.g., Arzoumanian et al., 2020;Chen et al., 2021;Goncharov et al., 2021). While the sensitivity of current PTAs to a GW signal scales favourably with time, it can be greatly accelerated by adding more pulsars to the array and increasing the cadence (Siemens et al., 2013). As Parkes is the largest telescope in the southern hemisphere currently contributing to PTA efforts, MeerTime can contribute significantly to the IPTA by ) observing the current IPTA pulsars with higher precision, ) observing additional pulsars not currently included in the IPTA due to low precision, and ) helping offset the recent loss of the Arecibo 305-m telescope.
Throughout this paper, we define MSPs as having rotational periods < 50 ms and period derivatives < 2 × 10 −17 s s −1 . In section 2, we discuss the sources included in this census and the observational parameters. We further discuss the polarization properties in section 3, the flux density measurements in section 4, and the preliminary timing results in section 5. We conclude with some predictions for the SKA in section 6. Information on accessing FITS-format files of integrated pulse profiles as well as the full table of results is provided at the end of this manuscript.

METHODS
All observations included in this paper were performed with the MeerKAT L-band receiver (856-1712 MHz) with the Pulsar Timing User Supplied Equipement (PTUSE) backend that uses coherent dedispersion. A full description of the observing system can be found in Bailes et al. (2020). Typical observations used 54-61 dishes resulting in a system equivalent flux density of 6.3-7.2 Jy near 1400 MHz. The system temperature varies as a function of frequency and is slightly worse near the band edges but difficult to measure in the areas of the band permanently affected by radio frequency interference (RFI) (Bailes et al., 2020;Geyer et al., 2021). While the MeerKAT Ultra-High Frequency (UHF) receiver at 544-1088 MHz is now fully commissioned for pulsar timing, time constraints did not allow for a full MSP census with the UHF receiver. Using the results of this L-band census, a smaller census of low-DM, steep-spectrum sources with the UHF receiver may be feasible, and such sources may be observed with the UHF receiver for the MPTA. The same applies for high-DM, flat-spectrum sources with the MeerKAT S-band receiver in future.
As the original network interface card for the first PTUSE machine was unable to ingest the full bandwidth in "1K mode" (using 1024 channels across the 856-MHz band), observations between mid-April 2019 and February 2020 were restricted to 928 frequency channels (dropping 48 channels from the top and bottom of the band where the roll-off adversely affects sensitivity in any case), reducing the recorded bandwidth to 775.75 MHz. In order to maintain consistency throughout our analyses, we therefore chose to reduce the observations with 856 MHz of recorded bandwidth down to the same 775.75 MHz. The loss in sensitivity is nearly negligible because of the sharp roll-off at the edges of the receiver band; comparisons of / measurements of a sample of 20 pulsars (250 observations) show a median sensitivity loss of only 2.2% (mean of 4.5%) . Almost all of our pulsars have negative spectral indices (typically between −1 and −3), which makes the loss of the lower part of the band worse than the top of it. Broader bands also make determination of DMs more accurate, which in turn aids in precision timing. Since 20% of the observations in our sample of 250 showed an increase of up to ≈ 5% in / after trimming the band edges, which is likely an artefact of band-pass calibration over-weighting edge channels.
late Feb 2020, all of our observations have used the entire 856 MHz of bandwidth.
Our source list was derived from the ATNF pulsar catalogue (using v.1.63;Manchester et al., 2005), limited to radio-loud MSPs not associated with globular clusters (GCs) and visible to MeerKAT (Dec ≤ +30 degrees to avoid severe time restrictions). We further restricted the list to exclude sources with positional uncertainties greater than 2 arcsec (near the resolution of the tied array beam) or without complete timing solutions (ephemerides) available to us, and we removed 12 Northern sources (Dec > 0 degrees) with known low flux densities ( 1400 < 0.1 mJy), under the assumption that these could be observed with greater efficiency by Northern telescopes and were unlikely to efficiently contribute to the timing array . Finally, some sources were removed from our list (after one or more initial observations) due to insufficient precision in ephemeris parameters (leading to significant "smearing" in individual 8-second sub-integrations) or due to low / in 1024-second observations . In total, 33 sources were excluded from our original list. To this list, we added 4 unpublished sources (or sources with currently unpublished timing solutions) matching the same sample limits; these MSPs were discovered in the Einstein@Home reprocessing of the Parkes Multibeam Pulsar Survey (Knispel et al., 2013) and the Green Bank North Celestial Cap survey (Stovall et al., 2014). We list our sources with basic parameters, including number of observations, and references in Table 1, where the DMs indicated are fitted from the brightest epoch per source . We also plot the sources in a -diagram in Figure 1, where is the rotation period of the pulsar and the dot denotes its time derivative. We note that, although initial ephemerides from published sources, the ATNF pulsar catalogue, or private communication were used as a starting point, we regularly refined ephemerides to ensure high-quality observations. These updated ephemerides will be presented in future work (Miles et al., 2022, in prep.).
There remain 10 pulsars in this study with Dec > 0 degrees and 1400 < 0.1 mJy as none of these had flux densities in v.1.63 of the ATNF pulsar catalogue but the / was sufficient in 1024-second observations for inclusion.
Some of the 33 sources that were excluded from the list of census sources were later observed to good precision after, for example, improving the ephemeris, but there were insufficient usable observations during the selected time range for this census We chose to provide DMs from single observations to demonstrate the precision achievable with the MeerKAT L-band receiver. Long-term timedependent variations, where significant, are accounted for in our ephemerides using "DM1" and "DM2" parameters. We include the epoch for each measured DM, and further information, in a supplementary table.  The 189 pulsars included in the MeerTime census with basic parameters, measured RMs with 1-uncertainties, percentages of linear polarization ( / ), circular polarization ( / ) and absolute circular polarization (| |/ ), median ToA uncertainties normalised to 256 seconds (see section 5 for details), flux densities at 1400 MHz, spectral indices from the power law model fit, and references (initial publication). The number of observations included in the analysis for each pulsar is given in the second column. Values less than 1-for constrained-positive parameters are given as 2-upper limits. The pulsars included in the regular observing programme (88 out of the 89) are indicated with asterisks. RMs measured from summed observations (not the mean of RMs from multiple observations) are indicated with a dagger ( †).    Science observations for MeerTime commenced on 2019 February 12, and observations for this census were of a high priority, beginning with the lowest-declination sources that are inaccessible from telescopes in the Northern hemisphere. The main motivation for the census was to quickly establish which pulsars should be included in the PTA programme. For each source, we planned a minimum of 6 observations in order to avoid random scintillation extrema biasing our results when deriving their timing potential and flux densities; due to delays and data quality issues, 4 sources have only 5 observations, and 1 was observed only 4 times. Sources that were found to have good timing precision (< 1 μs uncertainty in < 2000 s) were added to a list to be observed with a ∼ 2 week cadence for the MPTA, with source observing times set to achieve sub-μs timing precision on average (with a minimum observing time of 256 s to ensure a large fraction of telescope time is spent on-source). This list is updated regularly, to add new pulsars found to have high timing precision and to remove those found not to contribute significantly to the PTA. There are currently 89 that are in the MPTA, and these, excluding one which was included in this analysis due to issues with the ephemeris, are denoted by an asterisk in Table 1 below. In total, we analysed 3966 observations of 189 MeerTime MSPs from February 2019 to February 2021. This number includes observations performed for the MeerTime Relativistic Binaries sub-project (Kramer et al., 2021), with which there is considerable overlap of sources; data are shared between MeerTime sub-projects to maximise the science output of the entire Large Survey Project.
The processing of all observations for this project used (Hotan et al., 2004;van Straten et al., 2012), in a pipeline developed by Parthasarathy et al. (2021) which includes RFI excision with a modified version of the C -G software (Lazarus et al., 2020) that is optimised for the RFI present at MeerKAT (see, e.g., Bailes et al., 2020). Observations taken prior to mid-April 2020 were polarizationcalibrated with the tool using post-facto derived calibration solutions provided by the South African Radio Astronomy Observatory (SARAO), whereas later observations are polarization-calibrated at the telescope prior to the data being received by the PTUSE computers (see Serylak et al., 2021, for a full descrption of both calibration methods). All polarization values and derivatives follow standard pulsar conventions as described in van Straten et al. (2010).
In selecting pulsars for the regular timing programme, southern sources were prioritised over northern sources. Some northern sources with very high timing precision, such as PSR J1713+0747, have been included.
Due to various possible issues with the data, some observations were only usable for one part of this analysis or another. For example, an observation with failed polarization calibration might still be useful for the timing analysis, which only used Stokes I. The number of such observations is ∼ 200. The numbers of observations listed in Table 1 are the total unique observations used in any part of this analysis.
The observed polarization position angles are measured from celestial North and increasing towards East, and the IEEE definition of circular polarization is used. Stokes V is defined as left-circular minus right-circular.

POLARIZATION
For each pulsar, we produce a high / profile by summing the brightest individual observations partially averaged in frequency to 232 channels (≈ 3.34 MHz channel widths), and full time integration after summing. We show the polarization profiles in Figures 7-27, integrated to the centre frequency (after correcting for Faraday rotation as described below). All calculations of linear polarization, L, have been de-biased with a procedure modified from Everett & Weisberg (2001), where, for every bin across the profile, where meas is the measured linear polarization ( 2 = 2 + 2 ) and is the root-mean-square of the total intensity profile in the off-pulse region . We determine the uncertainties on the polarization position angle (PA, ) following Everett & Weisberg (2001): for each bin across the profile, if 0 = true, / > 10, the uncertainty is , = /2 ; otherwise, we determine the uncertainty numerically from the probability distribution (Naghizadeh-Khouei & Clarke, 1993) where = 0 √ 2 cos 2( − true ) and erf ( ) is the Gaussian error function. The uncertainty, , can therefore be determined by setting true = 0 and adjusting the bounds of integration of the probability distribution to satisfy When analysing our pulse profiles with , it is necessary to consider the noise baseline in each polarization channel, especially with complicated profiles with emission across the entire phase range. Accurate calibration will set the noise values appropriately, but, when calculating / or polarization fractions, will "subtract the baseline" by default. For all but one pulsar in our sample, we use the "minimum" baseline estimator , which finds a continuous region with the minimum mean by smoothing the profile with a boxcar with custom widths varying from 1% to 90% (where the region width was determined manually for each pulsar). One pulsar for which the "minimum" baseline estimator To make the brightest profiles, we summed observations in descending order of / until the next observation satisfied: / obs < / max /3 and either / obs < 5 or the number of observations summed was greater than 7.
Note that Everett & Weisberg (2001) apply a lower limit of 0 to each true,i , but we have chosen not to do so, just as the values of Stokes I can be negative.
was not optimal is PSR J1421−4409 (Spiewak et al., 2020), which instead used the "normal" estimator . Note that this pulsar still shows apparent over-polarization at pulse phases ∼ 0.6-1.0 (shown in the central panels of Fig. 13), as there is significant emission in both Q and U throughout that range. This pulsar may be similar to PSR J0218+4232 with underlying unpulsed emission (Navarro et al., 1995), which would contradict the general assumption that all polarization channels are consistent with noise over some phase range (the off-pulse baseline). Further study of this pulsar with MeerKAT could reveal the nature of this emission.
We measure significant fractional linear polarization ( / ; from fully time-and frequency-integrated profiles centred at 1284 MHz, as shown in Fig. 7-27) for 180 pulsars in our sample (> 1-significance), and significant absolute fractional circular polarization (| |/ ) for 119. Our full results are listed in Table 1.
We calculate the RM for each pulsar from the individual observations (fully integrated in time, with 29 frequency channels across the 775.75 MHz band) or, for faint pulsars, from the summed observations (also fully time-integrated, with 29 frequency channels) and list all in Table 1. We use to determine RMs using a brute force linear polarization maximisation which is then refined using an iterative fit to the measured PA as a function of frequency. After calculating RM from the individual observations, we then correct for the contributions of the ionosphere using the F R software (Sotomayor-Beltran et al., 2013) with IGSG maps . For each pulsar with ≥ 10 corrected RMs, we remove outliers using a basic drop-out algorithm comparing the standard deviation and median-absolute-deviation from the median. Outliers are removed manually for pulsars with fewer observations. We then report the mean of the corrected RM values with the standard deviation as the 1-uncertainty. For faint pulsars, because we measure the RM from a summed profile, we report the RMs and uncertainties as given by after correcting by the mean of the ionospheric contributions per pulsar (the uncertainties on the ionospheric corrections are summed in quadrature with the uncertainties given by ). No RMs are given for pulsars with linear polarization fractions less than ≈ 0.05 (∼ 1.5-).
In total, we measure 171 non-zero RMs (at the 1-level; 163 with 3-measurements), increasing the number of MSPs with known RMs by 88. To verify the reliability of our measurements, we compare the RMs in the ATNF pulsar catalogue (v1.64) for the 83 MSPs with previously published values and calculate the significance of the difference of the values: = (RM psrcat − RM census )/ √︃ 2 psrcat + 2 census . Nearly 50% of the pulsars have RMs that are 1-consistent with the previously published values, and 85% of the pairs are We found that the "normal" baseline estimator performed worse than the "minimum" estimator with custom widths for most pulsars with 95% duty cycles.
IGSG maps for the ionospheric correction were downloaded from NASA's CDDIS database at https://cddis.nasa.gov/archive/slr/ data/npt_crd/lageos1 (Noll, 2010). consistent to 2.6-. Given that pulsar RMs are known to vary with time due to inaccurate models of the ionospheric Faraday rotation (e.g., Yan et al., 2011), that not all values reported in the ATNF pulsar catalogue are corrected for the ionosphere, and that measurements at different frequencies can result in statistically inconsistent RMs (see, e.g., the position angle fits in Dai et al. 2015), differences between our RMs and published values at this level are not unexpected. However, two pulsars have significantly different values in our census and the literature: PSR J1502−6752 has a published value of −225(2) rad m −2 from Keith et al. (2012), while our measurement is 232.8(2) rad m −2 ; and PSR J0154+1833 has a published value of 21.6(1) rad m −2 (Martinez et al., 2019), while our measurement is −21.9(6) rad m −2 . Given the consistency of our other RM measurements, we are inclined to believe the Keith et al. (2012) and Martinez et al. (2019) values suffer from a sign convention error.
We show in Figure 2 the positions of the MSPs in our sample, coloured by their measured RM (black marks indicate the positions of pulsars without measured RMs). Overall, the RMs for our sample do not indicate strong (average) magnetic fields along the lines of sight, with 90% of values in the range −2.8 -2.3 μG (full range −6.7 -5.0 μG).

FLUX DENSITIES AND SPECTRAL INDICES
Ideally, pulsar flux densities are obtained using an injected and pulsed broadband calibration signal that can be measured against a source of known flux density to derive a flux scale. While attempting to derive a polarisation calibration solution using the MeerKAT L-Band receiver pulsed cal, we found that the degree of polarisation of the cal often exceeds 105%. This is probably due to the strength of the cal being comparable to the system equivalent flux density and driving the system into a non-linear regime. This made this conventional method unsuitable for absolute flux calibration.
An alternative calibration method involves using the radiometer equation and assuming that the rms of the pulsar baseline region is well-defined by the system equivalent flux density plus the galactic sky temperature divided by the antenna gain. In the presence of unrecognised sources of noise due to RFI, this method can underestimate radio fluxes and must be approached with caution. Nevertheless, as shown in Bailes et al. (2020), there are many parts of the MeerKAT L-band (856-1712 MHz) that are almost always devoid of interference, and we found that many of our high DM (DM 100 pc cm −3 ) sources gave very reliable / and hence flux values from epoch to epoch across much of the band. To assess the reliability of this method, we computed the modulation indices of the derived flux densities for the five pulsars with the largest DMs in our sample and found them to be between 7-15%, which means that the reliability of the flux densities from epoch to epoch is probably no worse than this figure.
Absolute flux densities for each observation were thus calculated using the radiometer equation and assuming the Figure 2: The positions of the MSPs in this sample on an Aitoff projection, with coloured circles indicating sources with measured RMs (including those formally consistent with zero) and black 'x' markers indicating sources without measured RMs, as described in the text. Haslam et al. (1981) sky temperatures scale as the observing frequency to the power −2.55 (Lawson et al., 1987). To determine mean pulsar flux densities the integrated flux above a baseline is often computed by firstly converting the counts to Jy in a calibration procedure. The mean pulsed flux is then just the integral of the flux above a defined baseline divided by the number of bins. Millisecond pulsars can have profiles that are very complex and this makes baseline estimation difficult if using a simple boxcar, especially if it extends into the wings of the pulse profile. For every MSP, the ideal boxcar width will differ, or even require multiple disjoint sections, so we used an interactive tool to define the baseline through visual inspection for each pulsar from high / observations to help determine the rms of the noise in the baseline. Software was written to cross-correlate the template with each profile, and the rms residual of the phase-aligned profiles in the baseline region was evaluated.
There are sections of the MeerKAT band that are rarely affected by RFI (see Bailes et al. 2020), and these were used to determine the flux density in those bands. We found that many of the high-DM MSPs that would be expected to have consistent epoch to epoch flux densities often exhibited dips in their mean spectra where RFI is regularly present in the fourth of the 8 sub-bands. As an example, for PSR J1017−7156, the off-pulse rms was roughly 5-21% higher in this band with a mean 10% higher after our attempts at RFI removal with C G . Scaling the flux densities of all channels by the ratio of their off-pulse rms counts to that of the median off-pulse rms counts led to mean spectra largely free of this dip in channel 4 and was performed on all the data. After integrating all of the observations for each pulsar, the vast majority of our pulsars were shown to have a smooth power law variation of mean flux density with frequency as is often seen in the population between 1-1.7 GHz, especially those at high dispersion ( 100 pc cm −3 ) where strong scintillation events across ∼100-MHz bands are not seen (Dai et al. 2015; see also Fig. 3, showing PSRs J1017−7156, J1600−3053, and J2241−5236, with DMs of 94.2, 52.3, and 11.4 pc cm −3 , respectively). To derive spectral indices and normalise our flux density measurements to a standard frequency (1400 MHz), we fitted a power-law model, = 1400 ( /1400 MHz) , to the mean flux density in each of 8 sub-bands across the 775.75 MHz band , and our results are listed in Table 1. In Figure 3, we show an example of the flux density measurements and best-fitting models for 3 pulsars with varying DM values, selected for comparison with the study by Dai et al. (2015) on pulsars in the PPTA (Kerr et al., 2020).
We refer the reader to Gitika et al. (in prep.) for a thorough analysis of the flux density measurements from these and more recent MeerTime MSP observations, and here we provide a simple comparison of our results with overlapping samples from Dai et al. (2015) and Alam et al. (2021a) to verify the accuracy of our method. Note that the Dai et al. (2015) spectral indices were derived from flux density measurements in 3 bands centred at 732, 1369, and 3100 MHz, and therefore they may be sensitive to, e.g., spectral turn-over, which is possibly seen in their analysis of PSR J1600−3053. With this caveat, we find good agreement between our 1400 values and spectral indices and those of Dai et al. (2015), with PSR J1744−1134 showing the greatest disparities in both 1400 and spectral index (although the values are consistent at the 2-level and this pulsar has a very small DM which makes it prone to large flux density variations). The flux density of this pulsar in the 97-Mhz sub-band centred at 1429 MHz varied from 0.14 to 14 mJy with a mean of 3.7 mJy and a standard deviation of 4.0 mJy. Thus its difference in mean flux density is not that worrying.
Recently, the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) collaboration presented a comprehensive study of flux densities for its sample with which we have significant overlap (Alam et al., 2021a). Of the 15 non-PPTA pulsars for which we have common data, 12 had flux densities within ∼30% of our values and those with the largest deviations (differing by roughly a factor of 2 from our values) have low DMs. Alam et al. (2021a) pointed out that the median flux they obtained for PSR J2317+1439 (0.45 +0.35 −0.19 mJy) differed from the catalogue flux density of 4(1) mJy significantly . Our mean value 0.60(6) mJy (median of 0.37 mJy) is much closer to the Alam et al. (2021a) result, and we note that this pulsar scintillates strongly. The pulsar We performed the power-law fitting in log space, weighting each channel mean by the standard deviation of the observed values.
The ATNF catalogue flux density for PSR J2317+1439 was from Kramer et al. (1998) and represents the mean of their measured distribution. The median quoted in that work was 2 mJy, hinting at significant scintillation. Figure 4: Histogram of spectral indices. In green, we show the spectral indices for pulsars with DM > 100 pc cm −3 (and therefore low scintillation) and those with > 10 observations to reduce the impact of scintillation on the measurements. We show the remaining spectral indices in our sample in the blue histogram.
with the most observations (134 with reliable fluxes) was PSR J1909−3744. Alam et al. (2021a) measured a median flux density at 1400 MHz of 1.03 mJy, matching our median value precisely, which gives us great confidence in our relative flux density scales. This pulsar experiences extreme scintillation events, and our measured mean flux density was 1.80(9) mJy. Dai et al. (2015) recorded a mean flux density for PSR J1909−3744 of 2.5(2) mJy at 1400 MHz. Note also that different observing strategies can lead to biases in observed flux density distributions. MeerKAT observations use a queue system with minimal intervention (and NANOGrav observations are run with minimal intervention), while observations for the PPTA are done manually, and observers are encouraged to switch sources when the source / is low and to repeat observations when the / is high. This leads to better uncertainties on the times of arrival (ToAs) but biases the flux density distributions to higher values, which is consistent with our results.
Overall, we find a mean spectral index of −1.92(6), with a standard deviation of 0.6, for our sample of 189 pulsars, and we show the full distribution of spectral indices in Figure 4. This mean is consistent at the 2-level with that found by Dai et al. (2015), who found a value of −1.81(1), as well as previous results for MSP studies (Toscano et al., 1998;Kramer et al., 1999) and studies of normal pulsars (Jankowski et al., 2018). The typical uncertainties on the spectral indices for pulsars with > 10 observations are ∼ 0.1, with 95% of the uncertainties < 0.4. These results could be improved in future using the MeerKAT UHF and S-band receivers to increase the frequency coverage.
For all pulsars in this study, we list basic parameters in Table 1, and provide our measured median ToA uncertainties. The templates used for our standard timing procedure are produced from summed, temporally-integrated observations with 232 frequency channels across the 775.75-MHz band. From these high-/ observations, we then derive pulse "portraits" using the P P software (Pennucci, 2019) and produce noise-free templates with 8 sub-bands for timing . We use the "Fourier Domain with Markov chain Monte Carlo" (FDM) fitting algorithm from in to produce ToAs from archives partially-integrated to 256second sub-integrations and 8 sub-bands across the 775.75-MHz band (∼ 97-MHz sub-bands). We compute "typical" ToA uncertainties per pulsar by first normalising the measured uncertainties to the nominal minimum observing time, 256 s, then calculating the error on the weighted mean ToA for each sub-integration as where i represents the frequency channels, and finally returning the median of the resulting sub values as the typical ToA uncertainty. We note that this method approximates the ToA uncertainty achievable by proper wide-band timing using a 2-dimensional template (i.e., returning a single ToA per sub-integration). We summarise the resulting median ToA uncertainties in Figure 5, showing the cumulative histogram of the median uncertainties. Note that the median ToA uncertainties should not be confused with the ultimate fitted rms residuals, as they do not factor in pulse jitter nor the deleterious effects of red noise.

Timing Forecast
In assembling a large dataset of MSP data across a minimum timespan of six months, we are able to determine the timing potential for each source and inform future studies. The potential for MeerTime to contribute to global efforts to detect GWs with the IPTA can be estimated from these data. As mentioned in section 2, we continue to regularly observe pulsars that have low ToA uncertainties, < 1 μs, calculated using frequency-scrunched templates and observations with integration times, 256 < obs < 2000 s.
To compare the sensitivity of PTA observations with MeerKAT to other current programmes, and its potential impact on global efforts, we consider the detection of a GW background (assuming a strain spectrum of ℎ ( ) = True wide-band timing with the Pennucci (2019) method is left to future work.
Note that frequency-dependent profile evolution and scintillation can result in poor ToA fits, and therefore higher uncertainties, using frequencyscrunched templates than when using frequency-resolved templates. This effect, and the differing datasets used in this Census and the MPTA programme, leads to different numbers of pulsars with sub-microsecond timing precision. 1.9 × 10 −15 ( /1 yr) −2/3 ; e.g., Siemens et al. 2013) through correlations between pulsars in the array (Hellings & Downs, 1983). The background is chosen to be the amplitude of the apparent common red noise found in analysis of the NANOGrav 12.5-year data release (Alam et al., 2021b). We can approximate the sensitivity of a given PTA to both the auto-correlated signal (the "pulsar term") and the cross-correlated signal (the "Earth term"), following Siemens et al. (2013) . They construct matched detection statistics in the Fourier domain (see their Equation 17) and calculate how the signal to noise ratio of the detection statistic depends on the noise properties of the array.
The expected signal-to-noise ratio in the autocorrelations is where ( ) is the power spectral density of the GWB in frequency channel , where = 1, is a Fourier series starting 1/ at the reciprocal of the common effective observing eff .
( ) = ( ) + , ( ) + , is the total power (the sum of the GWB and the red noise ( , ) and the white noise ( , ) in the same channel for pulsar .
The expected signal-to-noise ratio in the cross correlations is found as the S/N of a matched filter with the Hellings-Downs correlation function Shannon et al. (2013) used the same approach to place the bound on a gravitational wave background. 14 R. Spiewak et al. where are the angular correlation coefficients (Hellings & Downs, 1983) = 3 4 1 − cos( ) log 1 − cos( ) We use the published source lists, median pulsar timing precisions, timespans, cadences, and observational parameters from the European Pulsar Timing Array (EPTA; Desvignes et al., 2016), NANOGrav (Alam et al., 2021b), and the PPTA (Kerr et al., 2020) to form a representative sample of the IPTA. For sources in multiple PTAs, the data set with the longest time span is used. We show how the / for this background increases with time in Figure 6.
For the NANOGrav sensitivity curves, we have assumed that the Arecibo timing programme ceased in August 2020 and that the pulsars are not observed elsewhere, which demonstrates the importance of Arecibo to international pulsar timing efforts. The modelling includes red noise in the sensitivity estimates, for pulsars where red noise has previously been detected (Arzoumanian et al., 2020;Goncharov et al., 2021). Undetected red noise is likely to be present in some of the pulsars timed by MeerKAT, which will affect our sensitivity estimates.
These calculations imply that the MPTA will be comparable in sensitivity to the PPTA by 2023, NANOGrav by 2024, and the EPTA by ≈ 2025. In less than four years, the MPTA will be contributing to the IPTA effort. As the MeerKAT timing array programme has observed a large number of pulsars from its start, the sensitivity to the cross correlated signal (dashed lines) is always comparable to the auto-correlated signal, similar to the EPTA with its 42 pulsars. In contrast, the PPTA project includes only 26 pulsars, and NANOGrav started out with a smaller number of pulsars before gradually increasing in size starting in 2009, and, in the most recent data release (12.5-yr), they listed 47 pulsars (Alam et al., 2021b). The MPTA could be expanded in the future with timing of additional pulsars in other bands or timing of MSPs discovered by MeerKAT and other facilities. Future work by Middleton et al. (2022, in prep.) will examine the sensitivity of the MPTA and potential improvements.

SUMMARY
In this work, we demonstrated the outstanding potential of MeerKAT for precision pulsar timing of MSPs with this initial census, analysing polarization properties, flux densities and spectral indices, and timing precision. We presented a dataset of nearly 4000 observations of 189 pulsars, spanning 24 months, including polarization-calibrated profiles and high-precision ToAs. Our timing simulations offer simple predictions of future results, and we predict that the MPTA will be a major contributor to global PTA efforts within the next 5 years. More work is forthcoming, including the first Figure 6: Comparison of PTA sensitivity to the crosscorrelated component of the GW background (solid lines) and the uncorrelated signal (dashed lines), using the current MeerTime timing programme with 89 pulsars (blue lines), the EPTA programme with 42 pulsars (orange lines), the NANOGrav programme with 47 sources (red lines), the PPTA with 26 (brown lines), or the IPTA as the union of the four (black lines).
MPTA data release with full timing and noise analyses (Miles et al., 2022, in prep.), an analysis of flux density measurements for MeerTime MSPs (Gitika et al., in prep.), and a study of the sensitivity of the MPTA and potential for improvements (Middleton et al., 2022, in prep.).

DATA AVAILABILITY
We have made available all 189 polarization profiles, as shown in Figures 7-27, in PSRFITS format with full phase resolution (1024 bins), 4 polarization channels, and 8 frequency channels across the 775.75-MHz band, as well as the table containing the analysis results. The table, in csv format with 189 rows, contains mean flux densities (and standard deviations) in 8 frequency bands per pulsar, spectral indices and fit 1400 values with uncertainties, median timing uncertainties as described in section 5, polarization fractions and rotation measures, dispersion measures and reference epochs, and basic pulsar properties (rotation properties and positions) for easy reference. The dataset will be made available prior to publication of this manuscript, under the DOI 10.5281/zenodo.5347875.