<figure><img src="image_1.png" /></figure>

HAL
open science

# The TRAPUM L-band survey for pulsars in Fermi-LAT gamma-ray sources

C.J. Clark, R.P. Breton, E.D. Barr, M. Burgay, T. Thongmeearkom, L.

Nieder, S. Buchner, B. Stappers, M. Kramer, W. Becker, et al.

► To cite this version:

C.J. Clark, R.P. Breton, E.D. Barr, M. Burgay, T. Thongmeearkom, et al. The TRAPUM L-band survey for pulsars in Fermi-LAT gamma-ray sources. Monthly Notices of the Royal Astronomical Society, 2023, 519 (4), pp.5590-5606. 10.1093/mnras/stac3742. hal-03925287

HAL Id: hal-03925287

https://hal.science/hal-03925287

Submitted on 22 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The TRAPUM L-band survey for pulsars in Fermi-LAT gamma-ray sources

C. J. Clark<sup>1,2,3★</sup> R. P. Breton<sup>3</sup> E. D. Barr<sup>4</sup> M. Burgay<sup>5</sup> T. Thongmeearkom<sup>3</sup> L. Nieder<sup>1,2</sup>

S. Buchner<sup>6</sup> B. Stappers<sup>3</sup> M. Kramer<sup>4,3</sup> W. Becker<sup>7,4</sup> M. Mayer<sup>7</sup> A. Phosrisom<sup>3</sup>

A. Ashok<sup>1,2</sup> M. C. Bezuidenhout<sup>3</sup> F. Calore<sup>8</sup> I. Cognard<sup>9,10</sup> P. C. C. Freire<sup>4</sup> M. Geyer<sup>6</sup>

J.-M. Grießmeier,<sup>9,10</sup> R. Karuppusamy<sup>4</sup> L. Levin<sup>3</sup> P. V. Padmanabh<sup>4,1,2</sup> A. Possenti,<sup>5</sup>

S. Ransom<sup>11</sup> M. Serylak<sup>12,13</sup> V. Venkatraman Krishnan<sup>4</sup> L. Vleeschower<sup>3</sup> J. Behrend,<sup>4</sup>

D. J. Champion<sup>4</sup> W. Chen<sup>4</sup> D. Horn<sup>6</sup> E. F. Keane<sup>14</sup> L. Künkel<sup>15</sup> Y. Men<sup>4</sup> A. Ridolfi<sup>5,4</sup>

V. S. Dhillon<sup>16,17</sup> T. R. Marsh<sup>18</sup> and M. A. Papa<sup>1,2</sup>

*Affiliations are listed at the end of the paper*

Accepted 2022 December 13. Received 2022 December 13; in original form 2022 October 14

## ABSTRACT

More than 100 millisecond pulsars (MSPs) have been discovered in radio observations of gamma-ray sources detected by the Fermi Large Area Telescope (LAT), but hundreds of pulsar-like sources remain unidentified. Here, we present the first results from the targeted survey of Fermi-LAT sources being performed by the Transients and Pulsars with MeerKAT (TRAPUM) Large Survey Project. We observed 79 sources identified as possible gamma-ray pulsar candidates by a Random Forest classification of unassociated sources from the 4FGL catalogue. Each source was observed for 10 min on two separate epochs using MeerKAT’s L-band receiver (856–1712 MHz), with typical pulsed flux density sensitivities of ~ 100 µJy. Nine new MSPs were discovered, eight of which are in binary systems, including two eclipsing redbacks and one system, PSR J1526–2744, that appears to have a white dwarf companion in an unusually compact 5 h orbit. We obtained phase-connected timing solutions for two of these MSPs, enabling the detection of gamma-ray pulsations in the Fermi-LAT data. A follow-up search for continuous gravitational waves from PSR J1526–2744 in Advanced LIGO data using the resulting Fermi-LAT timing ephemeris yielded no detection, but sets an upper limit on the neutron star ellipticity of 2.45 × 10⁻⁸. We also detected X-ray emission from the redback PSR J1803–6707 in data from the first eROSITA all-sky survey, likely due to emission from an intrabinary shock.

**Key words:** binaries: general – pulsars: general – pulsars: individual: J1036–4353, J1526–2744, J1803–6707 – gamma rays: stars.

## 1 INTRODUCTION

Observations by the Large Area Telescope (LAT; Atwood et al. 2009) onboard the Fermi Gamma-ray Space Telescope have led to the detection of gamma-ray pulsations from nearly 300 pulsars<fnref n="1" /> (The Fermi-LAT Collaboration, in preparation). These fall into two main classes: canonical pulsars that are still rapidly spinning-down from their initial rotation periods; and millisecond pulsars (MSPs) that have been spun-up (or ‘recycled’ Bhattacharya & van den Heuvel 1991) to rates of hundreds of rotations per second by accretion from an orbiting star (Smarr & Blandford 1976; Alpar et al. 1982).

Both pulsar classes have characteristic gamma-ray emission properties (curved spectra and low flux variability over time) that are distinct from those of other gamma-ray emitting objects (e.g. Ackermann et al. 2012). These properties can be used to identify promising

pulsar-like gamma-ray sources, which can then be targeted with radio telescopes to detect pulsations and confirm their nature (Ray et al. 2012). The few-arcmin localization regions of unassociated Fermi-LAT sources enable targeted, long and repeated observations of promising sources, and a higher detection efficiency than can be achieved when surveying a broad region of the sky.

This method has proven highly successful at discovering new MSPs; more than a quarter of the 400 MSPs known in the Galactic field have been discovered in Fermi-LAT sources<fnref n="2" /> (e.g. Cognard et al. 2011; Keith et al. 2011; Ransom et al. 2011; Kerr et al. 2012; Barr et al. 2013; Camilo et al. 2015; Cromartie et al. 2016) in a global effort coordinated by the Fermi Pulsar Search Consortium (Ray et al. 2012). These searches have been particularly effective in finding exotic but elusive ‘spider’ binary pulsars (‘black widows’ and ‘redbacks’, Roberts 2012) whose long radio eclipses due to diffuse intrabinary material make them easily missed in single-pass

<footnote>* E-mail: colin.clark@aei.mpg.de</footnote>

<footnote n="1">¹http://tinyurl.com/fermipulsars</footnote>

<footnote n="2">²http://astro.phys.wvu.edu/GalacticMSPs/</footnote>
untargeted surveys. Around two-thirds of the known spider binaries in the Galactic field were found by targeting *Fermi* sources².

Many new MSPs remain to be found amongst the LAT sources: pulsars make up around 6 per cent of the identified or associated<fnref n="3" /> sources in the recent 12-yr iteration (Abdollahi et al. 2022, hereafter 4FGL DR3) of the *Fermi*-LAT Fourth Source Catalog (Abdollahi et al. 2020, hereafter 4FGL), while over 2000 sources remain unassociated. Observed gamma-ray and radio fluxes from pulsars are not strongly correlated with one another (Abdo et al. 2013; The Fermi-LAT Collaboration, in preparation), and so prospects remain high for detecting radio pulsars even in faint new gamma-ray sources that have only recently been detected thanks to the exposure that has accumulated during *Fermi*'s ongoing all-sky survey. This also means that new radio MSPs discovered within faint new *Fermi*-LAT sources can still be bright enough to be valuable astrophysical tools, and indeed several new MSPs found using this method have been added to pulsar timing array (PTA) projects aiming to detect gravitational waves (e.g. Spiewak et al. 2022).

The potential importance of detecting new MSPs in *Fermi*-LAT sources is illustrated by the fact that several previous discoveries found in this manner now mark the extreme edges of the MSP population, and are therefore the best current probes for several fundamental astrophysics questions. These include: the fastest known Galactic MSP (PSR J0952−0607, Bassa et al. 2017); the pulsar binary system with the shortest known orbital period (PSR J1653−0158, Nieder et al. 2020b); and a group of massive black-widow MSPs that probe the maximum neutron star mass (Romani et al. 2022).

One crucial benefit of finding a new MSP within a gamma-ray source is that gamma-ray pulsations can often be detected and timed directly in the *Fermi*-LAT data. Initial timing solutions can be refined and extrapolated backwards to the start of the LAT data (which currently spans more than 14 yr), providing long and precise rotational ephemerides without the need for lengthy radio timing campaigns. Illustrating the potential scientific benefits of this capability, gamma-ray timing of the recently discovered black-widow PSR J1555−2908 (Ray et al. 2022) may have revealed a second, planetary mass object in a long period orbit around the inner binary system (Nieder et al. 2022). The LAT data have also even been recently exploited to build a gamma-ray PTA (The Fermi-LAT Collaboration 2022), whose sensitivity to a stochastic gravitational wave background may reach that of current radio PTAs within a decade. While searches for pulsations in the *Fermi*-LAT data itself can reveal new MSPs (e.g. Pletsch et al. 2012; Clark et al. 2018), trials factors and computational costs limit these searches to the brighter pulsars. The detection of gamma-ray pulsations from MSPs in binary systems (which most are) is also impossible without prior orbital constraints (Nieder et al. 2020a). Radio surveys and initial timing therefore remain critical for expanding the population of Galactic MSPs.

The efforts of searching for new MSPs in *Fermi*-LAT sources have recently been bolstered by new radio telescopes. These bring capabilities of observing in new parameter spaces (e.g. at low radio frequencies with LOFAR and GMRT, Bassa et al. 2017; Pleunis

et al. 2017; Bhattacharyya et al. 2021) or with greater sensitivity (e.g. FAST; Wang et al. 2021).

In this paper, we present the first results from one such new survey of unassociated *Fermi*-LAT sources, using the MeerKAT radio telescope (Jonas 2009; Jonas & MeerKAT Team 2016). The full MeerKAT array is around five times more sensitive than the Murriyang Parkes telescope (Bailes et al. 2020), the next most sensitive radio telescope in the Southern Hemisphere. The Transients and Pulsars with MeerKAT (TRAPUM) project is a large survey project using MeerKAT to search for new pulsars (Stappers & Kramer 2016). All TRAPUM observations target sky locations in which pulsars are particularly likely to lie: globular clusters (Ridolfi et al. 2021); nearby galaxies (Carli et al. 2022); supernova remnants, pulsar wind nebulae and other TeV sources; and GeV gamma-ray sources. A separate dedicated *L*-band survey of the Galactic plane (MMGPS-L) is also ongoing using the same instrumentation and processing infrastructure (Kramer et al. 2016, Padmanabh et al., in preparation). To date, the TRAPUM and MMGPS-L searches have discovered more than 150 new pulsars,<fnref n="4" /> the majority of which are MSPs. Here, we present the first results from TRAPUM's survey of unassociated *Fermi*-LAT sources, which led to 9 of these MSP discoveries.

The paper is organized as follows: Section 2 describes the survey setup (recording and processing infrastructure, target selection and observation strategy); Section 3 presents the new discoveries and subsequent investigation (localization, timing and multiwavelength follow-ups) and an estimate of the survey's sensitivity; and finally a brief discussion and conclusions follow in Sections 4 and 5.

## 2 SURVEY PROPERTIES

### 2.1 MeerKAT and the TRAPUM processing infrastructure

MeerKAT is a radio interferometer located in the Karoo, South Africa, consisting of 64 antennas with 13.5-m effective diameter. Here, we give a brief description of MeerKAT and the TRAPUM infrastructure used to perform our pulsar search observations. For a full technical description of the instrument, we refer the reader to Jonas & MeerKAT Team (2016). At the time of data taking (between 2020 June and 2021 February), two receivers were available: the *L*-band receiver operating between 856 and 1712 MHz and the Ultra High Frequency (UHF) receiver between 544 and 1088 MHz. This survey was conducted at *L* band, but follow-up observations of new pulsars were also made at UHF.

At the *L*-band centre frequency, a coherent tied-array beam produced using all 64 MeerKAT antennas has a typical full width at half-maximum of the order of a few arcseconds. Furthermore, the data rate of complex voltages from the antennas is too high to record while observing. Efficiently searching the several-arcmin localization region of an unassociated *Fermi*-LAT source in a single pointing therefore relies on the ability to form and record the Stokes intensities from a large number of coherent beams simultaneously. This capability is provided by the Filterbanking BeamFormer User Supplied Equipment (FBFUSE), a 32-node, GPU-based, software beamformer developed by the Max Planck Institute for Radio Astronomy (Barr 2018; Chen et al. 2021). FBFUSE coherently sums the channelized complex voltages in real time, using sky position-dependent complex weights computed by the purpose-built Mosaic software<fnref n="5" /> from a delay model provided by the MeerKAT

<footnote n="3">³Gamma-ray sources are only deemed 'identified' if they have pulsed gamma-ray emission, gamma-ray variability that correlates with that seen in other wavelengths, or are resolved and have an angular extent consistent with that of a known source seen at other wavelengths. If a known source is likely to be the source of the gamma-ray emission, but these conditions are not met, then the gamma-ray source is described as 'associated' but not identified.</footnote>

<footnote n="4">⁴http://trapum.org/discoveries/</footnote>

<footnote n="5">⁵https://github.com/wchenastro/Mosaic</footnote>
Science Data Pipeline. FBFUSE can produce channelized time series for up to 864 coherent beams, as well as one incoherent beam produced by summing the (real-valued) Stokes intensities from each antenna. FBFUSE’s beamforming algorithm requires a multiple of four antennas, and not all antennas are available for all observations. As such, we used either 56 or 60 antennas during our observations, depending on availability.

The channelized data from each coherent beam are then recorded on to a distributed 3.5-PiB file system accessible from the Accelerated Pulsar Search User Supplied Equipment (APSUSE) instrument, a second, 60-node computer cluster, with two NVIDIA GeForce GTX 1080 Ti GPUs per node, on which the pulsar search takes place. The number of beams that can be stored is limited by the data rate at which APSUSE can record, and so down-sampling in time from the native data rate is necessary for a large number of beams to be recorded. All of our observations used the 4096-channel MeerKAT F-engine channeliser mode, but beamformed spectra were down-sampled in time by a factor of 16 from the native time resolution (4096/856 MHz = 4.785 µs) to give 76-µs time resolution. With these time and frequency resolutions, up to 288 coherent beams could be formed and recorded. Additionally, downsampling in frequency by a factor of 2 allowed for 480 coherent beams with 2048 frequency channels. Both of these recording modes were used in our survey.

The full set of filterbank files for 480 coherent beams constitute 46 TiB of data per hour of observation, allowing at most 73 h of data in this format to be stored for processing on APSUSE. It was therefore necessary to process these data quickly, identify promising candidates for further follow-up, and delete the raw data to ensure there was sufficient storage space for other TRAPUM projects to continue observing. Only the filterbank files for beams in which promising candidates were identified were retained for later use.

## 2.2 Observing strategy

The very high gain and low system temperature of the MeerKAT array ($G = 2.8$ K Jy$^{-1}$ for the full array and $T_{sys} = 18$ K, Bailes et al. 2020) enable the detection of pulsars with low flux densities, even with short observation lengths. This motivated a strategy involving short pointings towards as many sources as possible. We chose 10-min observations, for which we still obtain flux density limit estimates of around 95µJy that compare favourably against previous surveys of Southern sources with longer observations (see Sections 3.4 and 4).

Previous surveys for MSPs within unidentified *Fermi*-LAT sources have revealed the importance of observing sources more than once to mitigate non-detections due to scintillation (Camilo et al. 2015) or unfavourable orbital phases (e.g. due to spider eclipses or ‘jerk’ effects Andersen & Ransom 2018). In this paper, we describe the first two passes of this survey, both performed at L band. At least two further passes are planned for each source at the UHF band, where the lower frequency will provide additional sensitivity to pulsars with steep spectra, but where propagation effects and dispersive smearing are larger.

During each observation, FBFUSE can be configured to form coherent beams that are distributed within the primary field of view either at pre-specified locations, or automatically using an optimal hexagonal tiling pattern. In this latter mode, used for all of our survey observations, the beam spacing is defined by an overlap parameter, which is the fractional sensitivity level of a 2D Gaussian function fit to the tied-array beam response, simulated by Mosaic (see Chen et al. 2021) at the centre of the frequency band, at the point mid-way between two neighbouring beams. The simple Gaussian model of

<figure><img src="image_3.png" /><figcaption>**Figure 1.** An example of the coherent beam tiling pattern, from the second observation of 4FGL J1623.9–6936. The greyscale image shows the fractional sensitivity according to the Mosaic beam simulation, maximized over neighbouring coherent tied-array beams, reaching 100 per cent at the centre of each beam. The solid red ellipse is the 95 per cent confidence region for the *Fermi*-LAT source. The dashed green circle shows the minimum region that we aimed to cover, and has a radius corresponding to the semimajor axis of an approximate 99 per cent confidence region. The blue star shows the location of PSR J1623−6936 discovered in this observation. The positional uncertainty (obtained in Section 3.2) is smaller than this marker.</figcaption></figure>

the earlier versions of Mosaic, used for all our survey observations, tended to overestimate the true overlap when using the full MeerKAT array, and so e.g. points mid-way between beams in a tiling with an intended 50 per cent overlap actually only achieved ~40 per cent sensitivity.

In the first pass of all our targeted sources, we used the 4096-channel mode, with a maximum of 288 coherent beams. The beam tiling patterns for each source were configured with a desired overlap of 50 per cent. We prioritized observing sources with larger positional uncertainties when they were at lower elevations, where the coherent tied-array beams covered a larger solid angle, to ensure that the full localization region could be covered. For the sources with the largest uncertainty regions, observations were scheduled to ensure a sufficiently low elevation that the beam tiling covered a circle with a radius at least as large as the semimajor axis of the elliptical 95 per cent confidence region (hereafter $r_{95}$) from 4FGL. No significant sensitivity penalty due to ground spillover is incurred for elevations above 30°. For better-localized sources, the outer beams were well outside the *Fermi* source region.

Further development of FBFUSE after our survey began provided the capability to alter the tiling overlap between sources within one observing block. For 75 per cent of the observations in the second pass, we used the 2048-channel, 480-beam mode and adjusted the tiling overlaps for each source to maximize sensitivity while ensuring that the coherent beams covered a circular region with radius at least $r = \sqrt{\log(0.01)/\log(0.05)}r_{95}$, the semimajor axis of an approximate 99 per cent-confidence ellipse. An example of the resulting tiling pattern is shown in Fig. 1. All sources had $r_{95} < 5$ arcmin, such that the tiled region covered a very small patch around the centre of the ~1 deg primary beam, meaning no significant sensitivity loss occurs for coherent beams that are not located at the central pointing position.
## 2.3 Search pipeline

The filterbank data from all coherent tied-array beams and incoherent beams were searched by a dedicated pipeline built around the peasoup<fnref n="6" /> GPU-accelerated pulsar search code, which performs an FFT-based acceleration search via time-domain resampling with incoherent harmonic summing (Barr 2020, and described in detail in Morello et al. 2019).

The data were de-dispersed at trial dispersion measures (DMs) up to 300 pc cm⁻³ with spacing ΔDM = 0.06 pc cm⁻³. This range covers the maximum DMs predicted by the Galactic electron density model of Yao, Manchester & Wang (2017, hereafter YMW16) in the direction of any of our target sources, and 70 of our 79 target sources have maximum predicted DMs that are less than half of this range (i.e. below 150 pc cm⁻³).

The acceleration search covered a range of $a \pm 50$ m s⁻², which is slightly less than twice as large as the maximum acceleration seen from any known fully recycled binary MSP (26.1 m s⁻² from PSR J0024−7204V, Ridolfi et al. 2016) in the ATNF Pulsar Catalogue (v1.65, Manchester et al. 2005). Higher accelerations have been observed from relativistic binaries, such as pulsars with massive white dwarf companions or double neutron star systems, but these are all mildly recycled systems, which typically do not have sufficient spin-down power to emit gamma-ray pulsations, although the Double Pulsar is one notable exception (Guillemot et al. 2013). At each trial DM, the spacing between acceleration trials was chosen according to the scheme described in Morello et al. (2019), which ensures that signals lying between two acceleration trials suffer a total smearing that is not more than 10 per cent larger than the other (unavoidable) smearing effects due to intrachannel dispersion smearing and finite time resolution.

After performing the acceleration search, a clustering algorithm was used to search for clusters of candidates from neighbouring coherent beams with similar periods and DMs. Candidates with periods close to known radio-frequency interference (RFI) signals were excluded. This clustering step also aims to distinguish between astrophysical and terrestrial RFI signals, based on how the signal-to-noise ratio (S/N) drops off in neighbouring beams as a function of the angular offset from the beam in which the signal is detected most strongly. With the exception of side-lobe detections of very bright pulsars, signals from astrophysical point sources should only be detected in single beams, or in a small number of beams close to their sky positions. Candidates are therefore rejected if they are detected in many non-neighbouring coherent beams with S/Ns that drop off too slowly with angular offset, as this is indicative of a terrestrial interference signal.

The data from the beam containing the strongest detection of each candidate cluster were then folded using the PulsarX software,<fnref n="7" /> and folded candidates were scored by the Pulsar Image-based Classification System (PICS) machine-learning classifier (Zhu et al. 2014). Candidates surpassing a conservatively low PICS-score threshold of 10 per cent, typically a few hundred per pointing, were then grouped for visual inspection. As mentioned previously, the full set of raw data could not be stored indefinitely, except for data from beams containing high-confidence pulsar candidates. Reduced data products, i.e. candidate lists, and folded data for every candidate, were retained from all observations.

## 2.4 Target selection

We built a list of observing targets by identifying pulsar-like unassociated Fermi-LAT gamma-ray sources within the 4FGL catalogue.

Several studies have had success in using machine-learning classification techniques to identify pulsar candidates from the population of unassociated sources detected by the Fermi LAT (Lee et al. 2012; Saz Parkinson et al. 2016; Luo et al. 2020; Finke, Krämer & Manconi 2021). These methods rely on the fact that gamma-ray pulsars have characteristic spectral properties that distinguish them from other gamma-ray emitting source classes: (1) gamma-ray pulsar spectra typically have significant curvature, in that they deviate from a simple power-law spectrum due to a sub-exponential cutoff at photon energies above a few GeV; and (2) gamma-ray pulsars are very stable emitters on long time-scales. The Random Forest algorithm (Breiman 2001) has been shown to perform well for the purpose of classifying Fermi-LAT sources (Saz Parkinson et al. 2016; Luo et al. 2020), and so we employed this method to rank sources from 4FGL.

We used five parameters from 4FGL for the classification: PLEC_SigCurve, the significance of the log-likelihood improvement when fitting the source spectrum with a (curved) sub-exponentially cutoff power-law spectrum typical of pulsars, compared to a simple power-law; Variability2 Index, the chi-squared value of the energy flux measured in 2-month time bins; Signif_Avg, the point-source significance (brighter sources can have higher curvature and variability significance, including this parameter in the ranking accounts for this); PLEC_Index and PLEC_Expfactor, the photon index and the pre-factor *a* in the exponential cutoff term $\exp(-aE^b)$, where $b = 2/3$ was used for all unassociated sources and all but six bright pulsars in 4FGL. These final two parameters encode the energy at which the gamma-ray spectrum peaks, which for pulsars typically lies between 0.5 and 4 GeV, but which can be at much higher energies for active galactic nuclei (AGNs) with curved spectra.

We combined the gamma-ray source classes listed in 4FGL into three broad categories: AGNs, pulsars, or 'other' (which contains e.g. supernova remnants and pulsar wind nebulae). The classifier was trained to identify these classes using sources with highly likely or confirmed associations listed in 4FGL. To evaluate the classifier performance, we removed 33 per cent of the 4FGL sources, chosen at random, trained the classifier on the remaining population, and compared the classifier results to the known association classes. For identifying pulsars, the classifier had a precision of 82 per cent (i.e. 82 per cent of sources predicted to be pulsars were in fact pulsars) and a recall of 71 per cent (i.e. 71 per cent of pulsars in the sample were correctly identified as such).

For each unassociated source in 4FGL the Random Forest algorithm estimates the probability, $P(psr)$, of this source being a pulsar. We used this list of unassociated sources, ranked by their predicted pulsar probability, and made further cuts to reduce the number of sources to observe. First, only sources with declinations below +20° were included. Next, we removed sources with Galactic latitudes within $|b| < 10°$ of the Galactic plane. This is because gamma-ray pulsars close to the Galactic plane tend to be slowly spinning young pulsars, which have much narrower radio beams, and are therefore more likely to be radio-quiet gamma-ray pulsars undetectable to our survey. The Galactic plane is also being surveyed with MeerKAT with similar sensitivity as part of the MMGPS-L survey up to $|b| < 5°$. Our survey was planned before the MMGPS-L survey region had been finalized, and hence our Galactic latitude cut aimed to avoid redundant observations. We then removed sources whose

<footnote n="6">⁶https://github.com/ewanbarr/peasoup</footnote>

<footnote n="7">⁷https://github.com/ypmen/PulsarX</footnote>
4FGL 95 per cent confidence regions had semimajor axes larger than 5 arcmin, which is the largest region that could be covered in a single pointing using 288 coherent beams overlapping at 50 per cent sensitivity. This semimajor axis cut removes around 40 sources that pass our other cuts, but several of these will be covered in future UHF observations, where the coherent beams are wider.

Finally, we removed any sources with *P*(psr) < 12 per cent, which leaves 79 sources to search, while retaining 95 per cent of the total sum of all *P*(psr) values. Summing the classifier probabilities gives a (crude) prediction that 38 pulsars (not necessarily all detectable in radio surveys) should exist within these sources, with only 2 expected pulsars in the remaining sources that did not pass the probability threshold. The sources that were searched, along with the pulsar probabilities predicted by the classifier, are listed in Table 1.

## 3 RESULTS

From our two-pass survey, nine candidate signals were identified in the final visual inspection step as likely being new millisecond pulsars, and all of these were later confirmed by additional detections (see Section 3.1). The pulse profiles for the newly discovered MSPs are shown in Fig. 2.

Our survey also independently detected PSR J0312−0921, a black-widow MSP discovered recently in a GBT observation of 4FGL J0312.1−0921 (Tabassum et al., in preparation). Shortly after our first observations, another two target sources, 4FGL J0802.1−5612 and 4FGL J1231.6−5116, were identified as young gamma-ray pulsars by the direct detection of gamma-ray pulsations in the *Fermi*-LAT data by *Einstein@Home*.<fnref n="8" /> Neither pulsar was detected in our survey, but this is not surprising, as only a very small fraction of young pulsars discovered in gamma-ray pulsation searches have been detected in radio observations (Wu et al. 2018).

### 3.1 Follow-up observations

To confirm the pulsar nature of the detected candidates, we performed dedicated follow-up observations with MeerKAT at both *L* band and UHF, and checked for archival data from previous search observations of the *Fermi*-LAT sources in which they were found. Three pulsars (J1757−6032, J1803−6707, and J1823−3543) were re-detected in archival data from the Parkes radio telescope.

The five high-confidence pulsar candidates detected in the first pass were removed from the scheduled second survey pass, so that they could be observed in a dedicated set of confirmation observations along with candidates from the second pass. In these confirmation observations, we employed a very dense tiling (with 90 per cent overlap) with a smaller number of beams around the location of the coherent beam in which the candidate was detected most strongly in the initial survey observations. This dense tiling ensured high sensitivity for re-detecting these candidates, while also enabling us to more precisely localize each pulsar using the method described below (see Section 3.2).

Eight of the pulsar candidates (excluding PSR J1036−4353, discussed below) were re-detected in these confirmation observations. One pulsar, PSR J1709−033 was only seen in the UHF observation, most likely due to unfavourable scintillation that was seen during observations in which this pulsar was detected.

One pulsar, PSR J1036−4353 was not included in the confirmation observations, as it had not been immediately identified as a candidate due to a bug in the folding pipeline that caused it to be folded with the wrong acceleration sign. Instead, it was identified at a later date when the folded archives were corrected. It was confirmed in a later UHF observation as part of the next stage of this survey, which will be presented in a future paper.

### 3.2 Localization

Following the detection of a pulsar in our survey, the sky position could be estimated to a much higher precision than the size of a coherent beam by triangulation using the measured S/Ns in neighbouring beams. The method with which we localized candidates from our survey is described in Bezuidenhout et al. (in preparation), based on the concept outlined in Obrocka, Stappers & Wilkinson (2015), which works as follows. Given the model of the coherent beam point spread function (PSF) provided by Mosaic, one can compute the expected ratio of the S/Ns that would be recovered in two neighbouring beams for a pulsar at any given point nearby. This expected ratio will match the observed ratio, within uncertainties, for a strip of positions between the two beams. These strips can then be computed for each pair of beams, and the location of the candidate can be inferred from where these strips all cross one another.

This procedure is implemented by the *SeeKAT* package.<fnref n="9" /> We used *SeeKAT* on each of the confirmation observations, which had the most dense beam tilings. The best-fitting positions and uncertainties were in good agreement across the three observations (one at *L* band, two at UHF), and so we combined these results by summing together the log-likelihood surfaces from each observation.

The best-fitting positions and uncertainties for each pulsar are given in Table 2, and an example of the localization log-likelihood surface is illustrated in Fig. 3. The *SeeKAT* positional uncertainties for new pulsar candidates are typically of the order of a few arc-seconds, sufficiently precise to enable the identification of potential counterparts in multiwavelength catalogues.

### 3.3 Timing

Follow-up timing campaigns have begun for all newly detected pulsars. Initially, all timing observations were performed at other telescopes. The pulsars visible from the Northern hemisphere are followed up mainly using the Nançay (for PSRs J1526−2744 and J1823−3543) and Effelsberg (for PSRs J1709−0333 and J1906−1754) telescopes at *L* band, while all other pulsars (and, initially, also those observed at other telescopes) are followed up at the Murriyang Parkes telescope, using the Ultra-wide-band Low (UWL) receiver (Hobbs et al. 2020), covering a frequency range from 0.7 to 4 GHz. Depending on the discovery S/N at MeerKAT or, when available, the S/N at the improved position, we follow-up the TRAPUM discoveries for 1–2 h. Observations are carried out in search mode. Whenever possible the pulsars have been observed with a pseudo-logarithmic cadence to help achieve phase coherence in our timing analysis. Phase-connected timing solutions for two MSPs, described in detail below, were obtained using the *Dracula* algorithm<fnref n="10" /> (Freire & Ridolfi 2018); this was necessary given the sparsity of the timing data in both cases.

<footnote n="8">⁸https://einsteinathome.org/gammaraypulsar/FGRP1_discoveries.html</footnote>

<footnote n="9">⁹https://github.com/BezuidenhoutMC/SeeKAT</footnote>

<footnote n="10">¹⁰https://github.com/pfreire163/Dracula</footnote>
<table><caption>**Table 1.** List of survey observations of Fermi-LAT sources. P(psr) is the probability predicted by the Random Forest classifier of each source being a gamma-ray pulsar (see Section 2.4). Sopt, S50, and S95 are sensitivity estimates (see Section 3.4) from the radiometer equation for a 500-Hz pulsar, with DM = 100 pc cm⁻³ and a 15 per cent duty cycle, based on optimistic (no losses due to finite frequency, acceleration, and DM trial spacings, and for a pulsar lying at the centre of a coherent beam), realistic (median losses), and conservative (95th centile losses) assumptions, respectively.</caption><thead><tr><th>4FGL source</th><th>P(psr)</th><th>Epoch 1</th><th>S<sub>opt</sub> (μJy)</th><th>S<sub>50</sub> (μJy)</th><th>S<sub>95</sub> (μJy)</th><th>Epoch 2</th><th>S<sub>opt</sub> (μJy)</th><th>S<sub>50</sub> (μJy)</th><th>S<sub>95</sub> (μJy)</th></tr></thead><tbody><tr><td>J0048.6-6347</td><td>0.53</td><td>59034.2131</td><td>53</td><td>114</td><td>135</td><td>59196.6117</td><td>50</td><td>88</td><td>97</td></tr><tr><td>J0139.5-2228</td><td>0.29</td><td>59123.8042</td><td>53</td><td>110</td><td>142</td><td>59196.6191</td><td>50</td><td>94</td><td>112</td></tr><tr><td>J0251.1-1830</td><td>0.13</td><td>59034.206</td><td>52</td><td>112</td><td>136</td><td>59202.7276</td><td>49</td><td>78</td><td>89</td></tr><tr><td>J0312.1-0921<sup>b</sup></td><td>0.94</td><td>59034.1976</td><td>52</td><td>112</td><td>132</td><td>59202.7204</td><td>49</td><td>86</td><td>95</td></tr><tr><td>J0414.7-4300</td><td>0.45</td><td>59034.1435</td><td>52</td><td>113</td><td>136</td><td>59202.6874</td><td>49</td><td>92</td><td>106</td></tr><tr><td>J0529.9-0224</td><td>0.3</td><td>59065.4275</td><td>49</td><td>120</td><td>144</td><td>59188.1007</td><td>51</td><td>97</td><td>117</td></tr><tr><td>J0540.0-7552<sup>c</sup></td><td>0.19</td><td>59034.221</td><td>52</td><td>113</td><td>138</td><td>59188.1084</td><td>50</td><td>89</td><td>105</td></tr><tr><td>J0657.4-4658</td><td>0.26</td><td>59065.4594</td><td>50</td><td>121</td><td>142</td><td>59188.1538</td><td>50</td><td>96</td><td>115</td></tr><tr><td>J0712.0-6431</td><td>0.17</td><td>59065.5499</td><td>50</td><td>122</td><td>189</td><td>59188.161</td><td>49</td><td>99</td><td>118</td></tr><tr><td>J0802.1-5612<sup>b</sup></td><td>0.91</td><td>59065.444</td><td>52</td><td>128</td><td>148</td><td>59182.2759</td><td>54</td><td>94</td><td>108</td></tr><tr><td>J0940.3-7610<sup>c</sup></td><td>0.94</td><td>59065.4674</td><td>51</td><td>126</td><td>144</td><td>59188.1155</td><td>50</td><td>81</td><td>93</td></tr><tr><td>J0953.6-1509</td><td>0.89</td><td>59065.4514</td><td>49</td><td>118</td><td>137</td><td>59182.2831</td><td>53</td><td>94</td><td>112</td></tr><tr><td>J1036.6-4349<sup>a</sup></td><td>0.45</td><td>59065.4908</td><td>53</td><td>130</td><td>148</td><td>59182.3887</td><td>54</td><td>94</td><td>112</td></tr><tr><td>J1106.7-1742</td><td>0.57</td><td>59065.498</td><td>50</td><td>122</td><td>140</td><td>59182.3508</td><td>53</td><td>94</td><td>113</td></tr><tr><td>J1120.0-2204<sup>c</sup></td><td>0.75</td><td>59065.5058</td><td>52</td><td>126</td><td>143</td><td>59188.1761</td><td>49</td><td>64</td><td>74</td></tr><tr><td>J1126.0-5007</td><td>0.21</td><td>59065.4357</td><td>50</td><td>119</td><td>224</td><td>59182.412</td><td>54</td><td>95</td><td>111</td></tr><tr><td>J1204.5-5032</td><td>0.5</td><td>59020.8389</td><td>50</td><td>113</td><td>182</td><td>59188.0925</td><td>51</td><td>99</td><td>115</td></tr><tr><td>J1207.4-4536</td><td>0.7</td><td>59065.4835</td><td>52</td><td>127</td><td>152</td><td>59250.2806</td><td>54</td><td>103</td><td>126</td></tr><tr><td>J1213.9-4416</td><td>0.23</td><td>59065.5129</td><td>51</td><td>122</td><td>140</td><td>59188.1316</td><td>50</td><td>87</td><td>99</td></tr><tr><td>J1231.6-5116<sup>b</sup></td><td>0.96</td><td>59065.528</td><td>53</td><td>130</td><td>152</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1303.1-4714</td><td>0.15</td><td>59065.4207</td><td>67</td><td>159</td><td>248</td><td>59188.1244</td><td>51</td><td>103</td><td>121</td></tr><tr><td>J1345.9-2612</td><td>0.34</td><td>59065.4747</td><td>49</td><td>116</td><td>133</td><td>59188.1389</td><td>51</td><td>94</td><td>113</td></tr><tr><td>J1400.0-2415</td><td>0.91</td><td>59020.7561</td><td>50</td><td>117</td><td>134</td><td>59188.1461</td><td>51</td><td>87</td><td>100</td></tr><tr><td>J1416.7-5023</td><td>0.55</td><td>59020.7634</td><td>53</td><td>125</td><td>147</td><td>59188.1686</td><td>54</td><td>95</td><td>108</td></tr><tr><td>J1450.8-1424</td><td>0.14</td><td>59065.5427</td><td>52</td><td>130</td><td>158</td><td>59188.1981</td><td>51</td><td>96</td><td>111</td></tr><tr><td>J1458.8-2120</td><td>0.46</td><td>59065.5578</td><td>51</td><td>122</td><td>139</td><td>59188.1909</td><td>51</td><td>83</td><td>95</td></tr><tr><td>J1513.7-1519</td><td>0.19</td><td>59020.7276</td><td>51</td><td>120</td><td>143</td><td>59188.2053</td><td>52</td><td>91</td><td>103</td></tr><tr><td>J1517.7-4446</td><td>0.24</td><td>59065.5205</td><td>57</td><td>143</td><td>168</td><td>59188.1837</td><td>55</td><td>109</td><td>130</td></tr><tr><td>J1526.6-2743<sup>a</sup></td><td>0.68</td><td>59020.7346</td><td>51</td><td>120</td><td>140</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1526.6-3810</td><td>0.23</td><td>59065.5353</td><td>57</td><td>140</td><td>168</td><td>59188.2124</td><td>55</td><td>110</td><td>132</td></tr><tr><td>J1539.4-3323</td><td>0.9</td><td>59020.7418</td><td>52</td><td>116</td><td>140</td><td>59182.3427</td><td>57</td><td>100</td><td>116</td></tr><tr><td>J1543.6-0244</td><td>0.83</td><td>59182.3063</td><td>54</td><td>94</td><td>109</td><td>59250.3519</td><td>55</td><td>104</td><td>125</td></tr><tr><td>J1544.2-2554</td><td>0.59</td><td>59020.749</td><td>51</td><td>120</td><td>138</td><td>59196.4976</td><td>52</td><td>86</td><td>97</td></tr><tr><td>J1612.1+1407</td><td>0.84</td><td>59020.7713</td><td>52</td><td>123</td><td>144</td><td>59196.4842</td><td>66</td><td>114</td><td>129</td></tr><tr><td>J1622.2-7202</td><td>0.37</td><td>59065.5655</td><td>53</td><td>130</td><td>153</td><td>59182.2916</td><td>54</td><td>95</td><td>112</td></tr><tr><td>J1623.9-6936<sup>a</sup></td><td>0.17</td><td>59020.7123</td><td>51</td><td>119</td><td>137</td><td>59196.5798</td><td>51</td><td>83</td><td>94</td></tr><tr><td>J1630.1-1049</td><td>0.57</td><td>59020.8087</td><td>51</td><td>122</td><td>166</td><td>59182.2992</td><td>56</td><td>98</td><td>114</td></tr><tr><td>J1646.7-2154</td><td>0.17</td><td>59020.7198</td><td>53</td><td>127</td><td>175</td><td>59182.3136</td><td>58</td><td>102</td><td>270</td></tr><tr><td>J1656.4-0410</td><td>0.18</td><td>59123.674</td><td>134</td><td>271</td><td>337</td><td>59182.3591</td><td>55</td><td>97</td><td>112</td></tr><tr><td>J1659.0-0140</td><td>0.61</td><td>59123.7308</td><td>57</td><td>116</td><td>146</td><td>59182.3735</td><td>56</td><td>99</td><td>114</td></tr><tr><td>J1709.4-0328<sup>a</sup></td><td>0.13</td><td>59123.7082</td><td>57</td><td>117</td><td>147</td><td>59196.5426</td><td>53</td><td>98</td><td>112</td></tr><tr><td>J1711.9-1922</td><td>0.35</td><td>59020.8238</td><td>54</td><td>127</td><td>154</td><td>59196.5498</td><td>54</td><td>90</td><td>101</td></tr><tr><td>J1717.5-5804</td><td>0.66</td><td>59020.7791</td><td>54</td><td>126</td><td>147</td><td>59196.603</td><td>55</td><td>92</td><td>103</td></tr><tr><td>J1720.6+0708</td><td>0.6</td><td>59020.7944</td><td>53</td><td>124</td><td>144</td><td>59196.521</td><td>54</td><td>97</td><td>112</td></tr><tr><td>J1722.8-0418</td><td>0.29</td><td>59123.7452</td><td>58</td><td>118</td><td>147</td><td>59196.557</td><td>53</td><td>85</td><td>96</td></tr><tr><td>J1727.4+0326</td><td>0.13</td><td>59182.3663</td><td>57</td><td>99</td><td>115</td><td>59250.3446</td><td>57</td><td>119</td><td>138</td></tr><tr><td>J1730.4-0359</td><td>0.86</td><td>59123.738</td><td>58</td><td>118</td><td>150</td><td>59196.5641</td><td>54</td><td>82</td><td>95</td></tr><tr><td>J1735.3-0717</td><td>0.32</td><td>59123.7816</td><td>59</td><td>119</td><td>147</td><td>59196.5956</td><td>55</td><td>99</td><td>112</td></tr><tr><td>J1747.6-0324</td><td>0.46</td><td>59123.7673</td><td>59</td><td>121</td><td>152</td><td>59182.3969</td><td>59</td><td>103</td><td>122</td></tr><tr><td>J1749.8-0303</td><td>0.22</td><td>59020.8015</td><td>55</td><td>128</td><td>146</td><td>59196.5884</td><td>56</td><td>95</td><td>105</td></tr><tr><td>J1757.7-6032<sup>a</sup></td><td>0.84</td><td>59020.8314</td><td>52</td><td>121</td><td>140</td><td>59182.321</td><td>56</td><td>100</td><td>120</td></tr><tr><td>J1803.1-6708<sup>a</sup></td><td>0.23</td><td>59020.8161</td><td>51</td><td>118</td><td>138</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1813.7-6846</td><td>0.12</td><td>59123.6858</td><td>56</td><td>115</td><td>146</td><td>59182.3281</td><td>55</td><td>98</td><td>118</td></tr><tr><td>J1816.4-6405</td><td>0.47</td><td>59123.6931</td><td>57</td><td>116</td><td>148</td><td>59182.3353</td><td>56</td><td>99</td><td>120</td></tr><tr><td>J1816.7+1749</td><td>0.2</td><td>59123.7009</td><td>58</td><td>118</td><td>147</td><td>59196.5063</td><td>53</td><td>94</td><td>105</td></tr><tr><td>J1818.6+1316</td><td>0.77</td><td>59123.7601</td><td>57</td><td>118</td><td>147</td><td>59196.5282</td><td>53</td><td>92</td><td>104</td></tr><tr><td>J1822.9-4718</td><td>0.2</td><td>59123.7159</td><td>56</td><td>117</td><td>145</td><td>59202.6346</td><td>55</td><td>94</td><td>104</td></tr><tr><td>J1823.8-3544<sup>a</sup></td><td>0.13</td><td>59034.0991</td><td>58</td><td>125</td><td>150</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1824.2-5427</td><td>0.61</td><td>59123.6786</td><td>55</td><td>114</td><td>143</td><td>59182.3811</td><td>55</td><td>96</td><td>111</td></tr><tr><td>J1827.5+1141</td><td>0.96</td><td>59123.7744</td><td>59</td><td>121</td><td>151</td><td>59196.5354</td><td>55</td><td>92</td><td>103</td></tr></tbody></table>
<table><caption>**Table 1** – continued</caption><thead><tr><th>4FGL source</th><th>P(psr)</th><th>Epoch 1</th><th>S<sub>opt</sub> (μJy)</th><th>S<sub>50</sub> (μJy)</th><th>S<sub>95</sub> (μJy)</th><th>Epoch 2</th><th>S<sub>opt</sub> (μJy)</th><th>S<sub>50</sub> (μJy)</th><th>S<sub>95</sub> (μJy)</th></tr></thead><tbody><tr><td>J1831.1-6503</td><td>0.67</td><td>59123.7231</td><td>55</td><td>114</td><td>143</td><td>59202.6716</td><td>52</td><td>73</td><td>85</td></tr><tr><td>J1845.8-2521</td><td>0.85</td><td>59123.8123</td><td>61</td><td>124</td><td>155</td><td>59202.6571</td><td>56</td><td>96</td><td>107</td></tr><tr><td>J1858.3-5424<sup>a</sup></td><td>0.71</td><td>59034.1361</td><td>54</td><td>118</td><td>140</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1906.0-1718</td><td>0.33</td><td>59034.1066</td><td>57</td><td>123</td><td>150</td><td>59182.4041</td><td>57</td><td>101</td><td>123</td></tr><tr><td>J1906.4-1757<sup>a</sup></td><td>0.14</td><td>59123.7962</td><td>58</td><td>120</td><td>151</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1913.4-1526</td><td>0.51</td><td>59123.789</td><td>58</td><td>120</td><td>150</td><td>59202.6428</td><td>54</td><td>95</td><td>105</td></tr><tr><td>J1916.8-3025</td><td>0.59</td><td>59123.7525</td><td>57</td><td>117</td><td>147</td><td>59202.6643</td><td>53</td><td>72</td><td>84</td></tr><tr><td>J1924.8-1035</td><td>0.78</td><td>59202.65</td><td>54</td><td>83</td><td>95</td><td>59250.31</td><td>59</td><td>102</td><td>128</td></tr><tr><td>J1947.6-1121</td><td>0.18</td><td>59034.1285</td><td>55</td><td>118</td><td>142</td><td>59202.6791</td><td>52</td><td>89</td><td>100</td></tr><tr><td>J2026.3+1431</td><td>0.19</td><td>59196.5136</td><td>51</td><td>90</td><td>100</td><td>59250.3282</td><td>55</td><td>110</td><td>130</td></tr><tr><td>J2043.9-4802</td><td>0.87</td><td>59034.1757</td><td>53</td><td>116</td><td>143</td><td>59202.7498</td><td>51</td><td>77</td><td>88</td></tr><tr><td>J2112.5-3043</td><td>0.96</td><td>59034.183</td><td>53</td><td>115</td><td>138</td><td>59202.6952</td><td>50</td><td>60</td><td>67</td></tr><tr><td>J2121.8-3412</td><td>0.23</td><td>59034.2283</td><td>53</td><td>114</td><td>141</td><td>59202.757</td><td>50</td><td>94</td><td>111</td></tr><tr><td>J2133.1-6432</td><td>0.77</td><td>59034.1903</td><td>53</td><td>115</td><td>137</td><td>59202.7354</td><td>50</td><td>81</td><td>91</td></tr><tr><td>J2201.0-6928</td><td>0.41</td><td>59034.1214</td><td>53</td><td>113</td><td>133</td><td>59202.7426</td><td>50</td><td>86</td><td>95</td></tr><tr><td>J2212.4+0708</td><td>0.92</td><td>59034.1595</td><td>53</td><td>113</td><td>133</td><td>59202.7041</td><td>50</td><td>87</td><td>99</td></tr><tr><td>J2219.7-6837</td><td>0.52</td><td>59034.1138</td><td>53</td><td>114</td><td>135</td><td>59196.4899</td><td>50</td><td>78</td><td>90</td></tr><tr><td>J2241.4-8327</td><td>0.33</td><td>59034.152</td><td>53</td><td>113</td><td>132</td><td>59202.7119</td><td>50</td><td>87</td><td>96</td></tr><tr><td>J2355.5-6614</td><td>0.24</td><td>59034.1679</td><td>53</td><td>115</td><td>137</td><td>59196.5724</td><td>50</td><td>91</td><td>104</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Pulse profiles of the newly discovered millisecond pulsars. Black and red curves show the profiles from L-band and UHF observations, respectively. For clarity, two identical pulses are shown, and pulses have been arbitrarily normalized and phase shifted to peak at phase 0.5.</figcaption></figure>

Obtaining a phase-connected timing solution for the remaining seven pulsars has required a dedicated timing campaign using the Pulsar Timing User Supplied Equipment (PTUSE) system developed for the MeerTime project (Bailes et al. 2020), which additionally provides coherent de-dispersion and full Stokes polarization information. This MeerKAT timing campaign, along with follow-up searches for gamma-ray pulsations from these seven pulsars, will be presented in a dedicated paper (Burgay et al., in preparation). For these seven pulsars, we give a preliminary timing solution in Table 2, obtained using PRESTO’s fitCircular_orbit.py routine to fit a sinusoidal modulation to the observed barycentric spin periods from multiple observations without requiring phase-alignment across

observations. We also plot the orbital properties of the 8 newly detected binary MSPs (excluding PSR J1709-0333 which appears to be an isolated pulsar) according to these ephemerides, in comparison to the overall MSP population, in Fig. 4.

### 3.3.1 PSR J1526-2744

Using the SeeKAT position as a starting point, we obtained a phase-connected timing solution using TOAs produced from the original L-band observation, the two UHF observations, as well as several Nançay and Parkes observations taken in 2021 as part of our dedicated timing campaign. This timing solution provided precise constraints on the pulsar’s orbital semimajor axis and time of ascending node, as well as a refined position, but the spin-down rate was not significantly measurable and was highly correlated with the declination.

We used the radio timing solution to search for gamma-ray pulsations in the photon arrival times measured by the Fermi LAT. For this, we used SOURCE-class gamma-ray photons from within a 3° region-of-interest around the SeeKAT position, with energies greater than 100 MeV and with a zenith angle below 105° according to the ‘Pass 8’ P8R3_SOURCE_V2 (Atwood et al. 2013; Bruel et al. 2018) instrument response functions.<fnref n="11" /> To increase sensitivity to faint pulsations, we used gtsrcprob to compute photon probability weights (Kerr 2011), using the spectral and spatial parameters for nearby sources from the 4FGL DR2 catalogue, as well as the gll_iem_v07.fits Galactic and iso_P8R3_SOURCE_V2_v1.txt isotropic diffusion emission models.

Folding the gamma-ray data over the validity interval of the radio timing solution did not yield a significant detection. To search for gamma-ray pulsations in earlier data it was necessary to search over a 5D parameter space (α, δ, ν, v̇, and $P_{\text{orb}}$). The pulsar’s orbital semimajor axis and time of ascending node were already constrained

<footnote><sup>a</sup>New MSP discovered in this work.</footnote>

<footnote><sup>b</sup>Identified as a gamma-ray pulsar by other groups after this work began.</footnote>

<footnote><sup>c</sup>Identified as a likely pulsar binary system through optical/X-ray observations, but without pulsation detections.</footnote>

<footnote n="11">¹¹See https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_essentials.html</footnote>
<table><caption>**Table 2.** Properties of newly discovered millisecond pulsars without phase-connected timing solutions. Initial orbital solutions are obtained from fitting sinusoids to observed barycentric spin periods as a function of time. These models are approximate, and so we only provide limited precision on orbital parameters and spin frequency. Positional parameters (RA and Dec.) are estimated from the measured S/N in neighbouring coherent beams using SeeKAT (see Section 3.2), except for J1036−4353 whose position is from Gaia DR3.</caption><thead><tr><th>Parameter</th><th>PSR J1036-4353</th><th>PSR J1623-6939</th><th>PSR J1709-0333</th><th>PSR J1757-6032</th></tr></thead><tbody><tr><td>RA, α (J2000)</td><td>10<sup>h</sup>36<sup>m</sup>30<sup>s</sup>21513(1)</td><td>16<sup>h</sup>23<sup>m</sup>51<sup>s</sup>41<sup>°</sup>33<sup>′</sup>33<sup>″</sup></td><td>17<sup>h</sup>09<sup>m</sup>32<sup>s</sup>79<sup>°</sup>18<sup>′</sup>18<sup>″</sup></td><td>17<sup>h</sup>57<sup>m</sup>45<sup>s</sup>53<sup>°</sup>22<sup>′</sup>24<sup>″</sup></td></tr><tr><td>Dec., δ (J2000)</td><td>-43°53′08″7252(2)</td><td>-69°36′48″3<sup>+</sup>2<sup>′′</sup>2<sup>″</sup></td><td>-03°33′17″7<sup>+</sup>6<sup>′′</sup>0<sup>″</sup></td><td>-60°32′10″7<sup>+</sup>2<sup>′′</sup>6<sup>″</sup></td></tr><tr><td>Dispersion measure, DM (pc cm<sup>-3</sup>)</td><td>61.1</td><td>46.4</td><td>35.7</td><td>62.9</td></tr><tr><td>Spin frequency, ν (Hz)</td><td>595.2</td><td>415.0</td><td>283.8</td><td>343.3</td></tr><tr><td>Orbital Period, P<sub>orb</sub> (d)</td><td>0.259</td><td>11.01</td><td>—</td><td>6.28</td></tr><tr><td>Projected semimajor axis, x (lt-s)</td><td>0.665</td><td>6.73</td><td>—</td><td>9.62</td></tr><tr><td>Epoch of ascending node, T<sub>asc</sub> (MJD)</td><td>59536.31</td><td>59192.91</td><td>—</td><td>59183.40</td></tr><tr><td>Minimum companion mass, M<sub>c, min</sub> (M<sub>☉</sub>)</td><td>0.23</td><td>0.19</td><td>—</td><td>0.43</td></tr><tr><td>Distance (from YMW16), d (kpc)</td><td>0.40</td><td>1.3</td><td>0.21</td><td>3.5</td></tr><tr><td>Parameter</td><td>PSR J1823-3543</td><td>PSR J1858-5422</td><td>PSR J1906-1754</td><td>—</td></tr><tr><td>RA, α (J2000)</td><td>18<sup>h</sup>23<sup>m</sup>43<sup>s</sup>06<sup>°</sup>14<sup>′</sup>14<sup>″</sup></td><td>18<sup>h</sup>58<sup>m</sup>07<sup>s</sup>92<sup>°</sup>26<sup>′</sup>26<sup>″</sup></td><td>19<sup>h</sup>06<sup>m</sup>14<sup>s</sup>94<sup>°</sup>11<sup>′</sup>11<sup>″</sup></td><td>—</td></tr><tr><td>Dec., δ (J2000)</td><td>-35°43′40″8<sup>+</sup>1<sup>′′</sup>4<sup>″</sup></td><td>-54°22′14″6<sup>+</sup>3<sup>′′</sup>8<sup>″</sup></td><td>-17°54′33″7<sup>+</sup>1<sup>′′</sup>8<sup>″</sup></td><td>—</td></tr><tr><td>Dispersion measure, DM (pc cm<sup>-3</sup>)</td><td>81.7</td><td>30.8</td><td>98.1</td><td>—</td></tr><tr><td>Spin frequency, ν (Hz)</td><td>421.4</td><td>424.5</td><td>347.7</td><td>—</td></tr><tr><td>Orbital period, P<sub>orb</sub> (d)</td><td>144.5</td><td>2.58</td><td>6.49</td><td>—</td></tr><tr><td>Projected semimajor axis, x (lt-s)</td><td>51.8</td><td>1.68</td><td>1.35</td><td>—</td></tr><tr><td>Epoch of ascending node, T<sub>asc</sub> (MJD)</td><td>59091.36</td><td>59564.98</td><td>59304.97</td><td>—</td></tr><tr><td>Minimum companion mass, M<sub>c,min</sub> (M<sub>☉</sub>)</td><td>0.27</td><td>0.12</td><td>0.05</td><td>—</td></tr><tr><td>Distance (from YMW16), d (kpc)</td><td>3.7</td><td>1.2</td><td>6.8</td><td>—</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** An example of pulsar localization using SeeKAT, from the three confirmation observations of PSR J1526−2744. The dashed black contours show the localization probabilities in three separate observations. The colour scale and the solid red contour lines show the joint likelihood (i.e. the product of the likelihoods from all three observations). Contour lines are at approximate 1σ and 2σ levels. The green star is at the Fermi timing position for this pulsar.</figcaption></figure>

precisely enough by the initial radio ephemeris that only one trial was required in these dimensions. The search was performed using the weighted H-test (de Jager, Raubenheimer & Swanepoel 1989; Kerr 2011), a statistic that normally performs an optimal incoherent sum of the Fourier power in the first 20 harmonics. In this case,

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Orbital properties of known binary MSPs. Each companion type, according to the ATNF Pulsar Catalogue (Manchester et al. 2005), is denoted by a different marker as noted in the legend. We mark two pulsars categorized as having He white dwarf companions in the ATNF Pulsar Catalogue, PSRs J1518+0204C and J1653−0158, as black widows following Pallanca et al. (2014) and Nieder et al. (2020b), respectively. Filled markers show MSPs that were found by targeting an unidentified Fermi-LAT gamma-ray source. The new MSPs from this work are highlighted with blue stars and labelled.</figcaption></figure>

we only summed over 3 harmonics, as detecting power in higher harmonics requires increased search grid density in each dimension, but gamma-ray pulsars have most power in lower harmonics. This search detected a significant pulsed signal with $H = 108.9$,
<table><caption>**Table 3.** Timing solution for PSR J1526−2744. Timing parameters are obtained from the gamma-ray timing analysis, with the exception of the DM measurement which is from the original MeerKAT discovery observation. Parameter values are in the Barycentric Dynamical Time (TDB) scale.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2"><strong>Timing parameters</strong></td></tr><tr><td>Data span (MJD)</td><td>54681–59476</td></tr><tr><td>Reference epoch (MJD)</td><td>59355.468037</td></tr><tr><td>RA, α (J2000)</td><td>15<sup>h</sup>26<sup>m</sup>45.103(2)</td></tr><tr><td>Dec., δ (J2000)</td><td>−27°44′05″91(8)</td></tr><tr><td>Dispersion measure, DM (pc cm<sup>−3</sup>)</td><td>30.95(3)</td></tr><tr><td>Spin frequency, ν (Hz)</td><td>401.7446020975(3)</td></tr><tr><td>Spin-down rate, ṝ (Hz s<sup>−1</sup>)</td><td>−5.71(1) × 10<sup>−16</sup></td></tr><tr><td>Orbital period, P<sub>orb</sub> (d)</td><td>0.2028108285(7)</td></tr><tr><td>Projected semi-major axis, x (lt s)</td><td>0.22410(3)</td></tr><tr><td>Epoch of ascending node, T<sub>asc</sub> (MJD)</td><td>59303.20598(1)</td></tr><tr><td colspan="2"><strong>Derived parameters</strong></td></tr><tr><td>Spin period, p (ms)</td><td>2.489143587192(2)</td></tr><tr><td>Spin period derivative, ṡ</td><td>3.537(6) × 10<sup>−21</sup></td></tr><tr><td>Spin-down power, Ė (erg s<sup>−1</sup>)</td><td>9.1 × 10<sup>33</sup></td></tr><tr><td>Surface magnetic field strength (G)</td><td>9.5 × 10<sup>7</sup></td></tr><tr><td>Light-cylinder magnetic field strength (G)</td><td>5.7 × 10<sup>4</sup></td></tr><tr><td>Minimum companion mass, M<sub>min</sub> (M<sub>☉</sub>)</td><td>0.083</td></tr><tr><td>Distance (from YMW16), d (kpc)</td><td>1.3</td></tr></tbody></table>

and with phase-connected pulsations visible from the start of the LAT data.

Following this detection, we derived a precise 13-yr gamma-ray timing solution by varying the timing parameters to maximize the unbinned Poisson log-likelihood of the weighted photon phases (Abdo et al. 2013) using a template pulse profile consisting of two wrapped Gaussian peaks whose parameters were also free to vary in the fit. The best-fitting parameter values and uncertainties are given in Table 3. The gamma-ray photon phases according to the best-fitting timing model and template pulse are shown in Fig. 5. We also tested for but did not significantly detect proper motion ($|\mu| < 49$ mas yr$^{-1}$) or eccentricity ($e < 8 \times 10^{-4}$).

The nature of the companion star in this system is currently unclear. The short orbital period (4.9 h) and low minimum companion mass (0.08 M<sub>☉</sub>) suggest a heavy black-widow or light redback companion, but no radio eclipses typical of these systems have been seen. Many black-widow binaries have detectable optical counterparts, and so we searched for an optical counterpart to this system using the ULTRACAM (Dhillon et al. 2007) high-speed multiband imager on the 3.5-m New Technology Telescope (NTT) at ESO La Silla. The longest observation lasted 3.5 h, covering orbital phases between 0.1 and 0.8. We did not find any counterpart at the pulsar’s timing position, with $3\sigma$ magnitude upper limits of $i_s = 23.2$, $g_s = 24.0$, and $u_s = 23.0$ in the deepest 5-min stacked image, which was obtained close to the companion’s superior conjunction, where a heated black-widow companion would appear at its brightest. It therefore seems likely that the companion is a light-weight ($M_{min} = 0.083$ M<sub>☉</sub>) white dwarf. If confirmed, e.g. through the detection of a non-variable optical counterpart below our ULTRACAM detection threshold, this would be the shortest orbital period of any known fully recycled Galactic MSP-WD system, and the third shortest orbit of any Galactic PSR-WD binary system, after the relativistic binary PSRs J0348+0432 (Antoniadis et al. 2013) and J1141-6545 (Kaspi et al. 2000).

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Gamma-ray pulsations from PSR J1526−2744 in Fermi-LAT data. The lower panel shows each photon’s phase according to the best-fitting timing solution. Each photon is assigned a weight, computed from its incidence angle and energy, and illustrated by greyscale representing the probability that it was emitted by the targeted source as opposed to a foreground background source. The upper panel shows the integrated pulse profile, computed by summing the photon probability weights in each rotational phase bin. The dashed blue horizontal line shows the estimated background rate. In the upper panel, the opaque orange curve shows the best-fitting template pulse profile used to estimate the log-likelihood of the gamma-ray pulsations. The transparent black curves show templates with parameter values randomly drawn from our MCMC samples to illustrate the uncertainty on the pulse profile.</figcaption></figure>

For the observed DM = 30.95(3) pc cm$^{-3}$, YMW16 predicts a distance of 1.3 kpc. The gamma-ray energy flux from this source above 100 MeV is $G_{100} = 2.5 \pm 0.4 \times 10^{-12}$ erg cm$^{-2}$ s$^{-1}$, corresponding to a gamma-ray luminosity of $L_{\gamma} = 5.3 \times 10^{32}$ erg s$^{-1}$. This can be compared to the pulsar’s spin-down power $\dot{E} = 4\pi^2 I \nu \dot{\nu} = 9.1 \times 10^{33}$ erg s$^{-1}$, for an assumed pulsar moment-of-inertia $I = 10^{45}$ g cm$^2$, corresponding to an efficiency $\eta = L_{\gamma}/\dot{E} \approx 0.06$, which is typical of gamma-ray MSPs (Abdo et al. 2013).

Using the timing solution obtained above, we carried out a search for continuous gravitational waves from this pulsar. Continuous gravitational wave emission is expected due to a non-axisymmetric deformation (ϵ) of the neutron star and the dominant signal frequency is twice its rotational frequency. Our search was therefore targeted at a frequency $f \approx 803.5$ Hz and spin-down $\dot{f} \approx -1.14 \times 10^{-15}$ Hz/s and directed at the timing position.
An *a priori* estimate of the detectable strain amplitude using Advanced LIGO data yields a value about an order of magnitude larger than the spin-down upper limit amplitude (Aasi et al. 2014),

$$h_0^{sd} = \left( \frac{5}{2} \frac{GI}{c^3} \frac{|\dot{v}|}{d^2 v} \right)^{1/2}, \qquad (1)$$

where we again assume $I = 10^{45}$ g cm$^2$. This is the strain amplitude of the signal assuming that all the rotational kinetic energy lost by the pulsar ($\dot{E}$) was converted into gravitational waves. We therefore do not expect a detection, but, in the spirit of leaving no stones unturned (e.g. Abbott et al. 2018, 2019; Nieder et al. 2019, 2020b; Abbott et al. 2022a,b), we carried out the search.

We used all of the Advanced LIGO data from the Hanford and Livingston detectors collected during the O1, O2, and O3 runs (Abbott et al. 2021). The gamma-ray timing solution allows us to coherently combine these data sets using a single template, thus achieving the maximum possible sensitivity this search could have to date. Loud detector glitches in the data were removed through gating (Steltner, Papa & Eggenstein 2022). The frequency range relevant to this search is free of known lines in the detectors. The multidetector, matched-filtering *F*-statistic (Cutler & Schutz 2005) was used for the analysis. Our search result has a *p*-value of 0.29, estimated using off-source data. This result implies a non-detection and based on it we set a 95 per cent confidence upper limit of $1.25 \times 10^{-26}$ on the intrinsic gravitational wave amplitude. Our upper limit is a factor of ≈20 larger than the spin-down upper limit, and the uncertainty on this upper limit is not more than 8 per cent including calibration uncertainties. Translating the amplitude upper limit into an upper limit on the ellipticity of the source, we constrain the ellipticity of J1526−2744 to be <2.45 × 10$^{-8}$, which is close to the minimum ellipticity proposed for millisecond pulsars by Woan et al. (2018).

We also searched in a band of $f$ and $\dot{f}$ accounting for mismatches between the phase of the gravitational wave signal and the phase locked to the electromagnetic observations. Such mismatches could result from a differential rotation between the parts of the star emitting the gravitational wave and the electromagnetic pulsations or if the star was biaxial and consequently freely precessing. We searched in a bandwidth of 0.4 per cent of the spin parameter values following Abbott et al. (2020), and estimated the *p* values of the results as done in Ashok et al. (2021). The results from this band search were also consistent with expectations from Gaussian noise.

### 3.3.2 PSR J1803−6707

PSR J1803−6707 was the first pulsar to be found in our survey, and was quickly confirmed by a detection in 1 h of archival Parkes data from 2015. Detection of the pulsar in this observation required a significant jerk term, indicating that the pulsar was likely to be in a short-period binary system. Pulsations were only detected for ~30 min in one of two 1 h dedicated follow-up observations with Parkes in 2020 December, indicative of a wide eclipse typical of black-widow or redback systems.

Establishing an orbital timing solution (given in Table 4) revealed an ~9 h orbit and a minimum companion mass of 0.26 M$_{\odot}$ (assuming a pulsar mass of 1.4 M$_{\odot}$), as well as significant variations in the orbital period requiring several orbital frequency derivative terms to describe the orbital phase throughout the 1 yr of phase-connected timing data. These features are characteristic of redback binary systems (e.g. Deneva et al. 2016).

We used the UHF MeerKAT observations to estimate the sky position of this pulsar using *SeeKAT*. The refined position is

<table><caption>**Table 4.** Radio timing solution for PSR J1803−6707. Parameter values are in the TDB scale.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2"><strong>Gaia astrometry</strong></td></tr><tr><td>RA, α (J2000)</td><td>18h03m04s235339(9)</td></tr><tr><td>Dec., δ (J2000)</td><td>-67°07′36″1576(2)</td></tr><tr><td>Proper motion in RA, μαcos δ (mas yr<sup>-1</sup>)</td><td>-8.4(2)</td></tr><tr><td>Proper motion in Dec. μδ (mas yr<sup>-1</sup>)</td><td>-6.5(2)</td></tr><tr><td>Parallax, ω (mas)</td><td>0.18(17)</td></tr><tr><td>Position reference epoch (MJD)</td><td>57388.0</td></tr><tr><td colspan="2"><strong>Timing parameters</strong></td></tr><tr><td>Data span (MJD)</td><td>59020–59583</td></tr><tr><td>Reference epoch (MJD)</td><td>59364.893677</td></tr><tr><td>Dispersion measure, DM (pc cm<sup>-3</sup>)</td><td>38.382(3)</td></tr><tr><td>Spin frequency, ν (Hz)</td><td>468.46771214886(7)</td></tr><tr><td>Spin-down rate, ṽ (Hz s<sup>-1</sup>)</td><td>-4.01(2) × 10<sup>-15</sup></td></tr><tr><td>Orbital period, P<sub>orb</sub> (d)</td><td>0.38047324(8)</td></tr><tr><td>Projected semimajor axis, x (lt s)</td><td>1.061910(3)</td></tr><tr><td>Epoch of ascending node, T<sub>asc</sub> (MJD)</td><td>59020.99710(1)</td></tr><tr><td>1st orbital frequency derivative, f<sub>orb</sub><sup>(1)</sup> (Hz s<sup>-1</sup>)</td><td>1.05(9) × 10<sup>-19</sup></td></tr><tr><td>2nd orbital frequency derivative, f<sub>orb</sub><sup>(2)</sup> (Hz s<sup>-1</sup>)</td><td>-1.5(1) × 10<sup>-24</sup></td></tr><tr><td>3rd orbital frequency derivative, f<sub>orb</sub><sup>(3)</sup> (Hz s<sup>-1</sup>)</td><td>1.03(7) × 10<sup>-31</sup></td></tr><tr><td>4th orbital frequency derivative, f<sub>orb</sub><sup>(4)</sup> (Hz s<sup>-1</sup>)</td><td>-3.3(2) × 10<sup>-40</sup></td></tr><tr><td colspan="2"><strong>Derived parameters</strong></td></tr><tr><td>Spin period, p (ms)</td><td>2.1346188308539(3)</td></tr><tr><td>Spin period derivative, ṡ</td><td>1.828(7) × 10<sup>-20</sup></td></tr><tr><td>Spin-down power, Ẽ (erg s<sup>-1</sup>)</td><td>7.4 × 10<sup>34</sup></td></tr><tr><td>Surface magnetic field strength (G)</td><td>2 × 10<sup>8</sup></td></tr><tr><td>Light-cylinder magnetic field strength (G)</td><td>1.9 × 10<sup>5</sup></td></tr><tr><td>Minimum companion mass, M<sub>min</sub> (M<sub>☉</sub>)</td><td>0.26</td></tr><tr><td>Distance (from YMW16), d (kpc)</td><td>1.4</td></tr></tbody></table>

coincident with the location of a star in the Gaia DR3 catalogue (Gaia Collaboration 2016, 2022). We also observed this star with ULTRACAM, which revealed optical variability with the same 9-h periodicity, confirming that this is indeed the optical counterpart. The optical light curve, shown in Fig. 6, varies by around 1.5 mag, with a single peak indicative of significant heating via irradiation from the pulsar (similar to e.g. PSR J2215+5135, Schroeder & Halpern 2014) but the counterpart is detectable at all orbital phases, as is typical for nearby redback companions which tend to have hot surfaces (T > 4000 K) even on the non-irradiated side. Light-curve modelling to estimate properties of the companion such as its temperature, irradiation, radius, and the binary inclination angle will be performed in a dedicated follow-up project (Phosrisom et al., in preparation).

As with PSR J1526−2744, we folded the gamma-ray data using the radio ephemeris to check for gamma-ray pulsations. Significant pulsations are detected within the 1-yr interval in which the radio ephemeris is valid, with a weighted *H*-test of *H* = 51, corresponding to a 5.6σ detection using the false-alarm probability calibration from Bruel (2019). These pulsations are shown in Fig. 7. However, these pulsations quickly disappear when extrapolating to earlier data. This is not unexpected, as Taylor-series models for orbital period variations lack predictive power, and the spin-down rate is not measured precisely enough in the radio data to fold the 13 yr of Fermi-LAT data. We have not been successful in extending the ephemeris by timing the gamma-ray data, likely as a result of the pulsar's faint gamma-ray flux and the large amplitude of the orbital phase variations. A
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Optical light curves for PSRs J1036−4353 and J1803−6707 obtained using ULTRACAM+NTT. Observations are folded on the orbital period using the radio timing ephemerides. The pulsar's superior and inferior conjunctions are marked with dotted and dashed vertical lines, and labelled SC and IC, respectively. The two peaks in the light curve for PSR J1036−4353 are due to ellipsoidal deformation of the companion star in the pulsar's gravitational field, while the single peak in the light curve for PSR J1803−6707 indicates strong heating from the pulsar.</figcaption></figure>

longer radio timing baseline will hopefully solve this problem by providing a more precise spin-down rate measurement, and a longer phase-connected radio ephemeris that will fold more *Fermi*-LAT data to provide better statistics with which to build a gamma-ray pulse profile template to search with. Nevertheless, this detection confirms the association between 4FGL J1803.1−6708 and PSR J1803−6707. The gamma-ray energy flux from this source above 100 MeV is $G_{100} = 5.0 \pm 0.5 \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$, corresponding to an efficiency $\eta = 0.017$ for the assumed YMW16 distance of 1.4 kpc (for DM = 38.382 pc cm$^{-3}$), which is again fairly typical of gamma-ray MSPs (Abdo et al. 2013).

Using data from the first *eROSITA* all-sky survey (eRASS1; Predehl et al. 2021), the X-ray counterpart of PSR J1803−6707 was independently detected (positional match within ~5 arcsec) in a pilot search for likely X-ray counterparts of unassociated *Fermi*-LAT sources. Fig. 8 depicts the X-ray counterpart as seen in eRASS1. The formal detection significance in the 0.2–2.3 keV bandwidth is $4.6\sigma$ for a vignetting corrected exposure time of 167.63s and a count rate of $(7.2 \pm 2.4) \times 10^{-2} \text{ ct s}^{-1}$ obtained from all seven telescope modules. The source is not detected in the 2.3–8.0 keV bandpass. The detected X-ray photons do not support a detailed spectral analysis. *eROSITA*'s temporal resolution is 50 ms which prevents the detection of a periodicity at the millisecond level as well.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Gamma-ray pulsations from PSR J1803−6707 in *Fermi*-LAT data. The lower panel shows each photon arrival during the radio timing solution's validity interval. Each photon is assigned a weight, computed from its incidence angle and energy, and illustrated by greyscale representing the probability that it was emitted by the targeted source as opposed to a foreground background source. The upper panel shows the integrated pulse profile, computed by summing the photon probability weights in each rotational phase bin. The dashed blue horizontal line shows the estimated background rate.</figcaption></figure>

Assuming a power-law spectrum with a photon index of 2.0 and an absorption column density of $6 \times 10^{20} \text{ cm}^{-2}$ (HI4PI Collaboration 2016), the *eROSITA* counting rate implies an unabsorbed X-ray flux of $F_X = (1.5 \pm 0.5) \times 10^{-13} \text{ erg s}^{-1} \text{ cm}^{-2}$ in the 0.2–10 keV band. For the Yao et al. (2017) distance, we obtain an X-ray luminosity of $L_X \sim 3.5 \times 10^{31} \text{ erg s}^{-1}$ and an X-ray efficiency of $L_X/\dot{E} \sim 4.7 \times 10^{-4}$, a typical value for X-ray detected MSPs (see e.g. Becker & Truemper 1997). A significant contribution by the $g \sim 20$ optical counterpart to the detected X-ray flux is not expected. The X-ray-to-optical flux ratio of the system is around 1, much greater than that expected for coronal X-ray emission from stars (Krautter et al. 1999).

It is interesting to note that, considering the $\gamma$-ray flux of the likely counterpart, the $\gamma$-ray to X-ray flux ratio of this source is comparatively small at $F_\gamma/F_X \sim 35$. This value is located at the lower end of the observed distribution for high-energy pulsars, which extends over a range of around $10–10^4$ (Marelli, De Luca & Caraveo 2011; Berteaud et al. 2021). It is therefore likely that the X-ray flux contains a contribution from, or is dominated by, an intrabinary shock, as is commonly seen in redback binary systems (e.g. Roberts
<figure><img src="image_10.png" /><figcaption>**Figure 8.** eRASS1 image in the 0.2–2.3 keV band, with the 95 per cent positional error ellipse of 4FGL J1803.1–6708 indicated in blue, and the GAIA position overlaid on the eROSITA counterpart of PSR J1803–6707 marked in red. The image has been spatially binned to a pixel size of 10 arcsec and smoothed with a 15 arcsec kernel. The color bar indicates the count rate in units of $10^{-3}$ ct s$^{-1}$ pix$^{-1}$.</figcaption></figure>

et al. 2015). Follow-up observations are planned to further investigate the X-ray emission from this system.

### 3.3.3 PSR J1036–4353

This pulsar is both the fastest-spinning (with $v \approx 595$ Hz), and the most highly-accelerated (with $a \approx 15$ m s$^{-2}$) pulsar found by our survey. Follow-up observations have revealed this to also be a redback system. A preliminary orbital solution revealed an ~6.2 h orbit and a minimum companion mass of 0.23 $M_{\odot}$. Several non-detections have occurred during our timing campaign in the half of the orbit in which the pulsar is behind its companion, likely due to eclipses by material from the redback companion. A single Gaia source lies within its SeeKAT localization region. Like PSR J1803–6707, we observed this source with ULTRACAM, revealing a double-peaked light curve, shown in Fig. 6, very typical of redback binary systems. Modelling of this data will also be presented in a future paper (Phosrisom et al., in preparation).

### 3.3.4 The other new MSPs

Of the remaining MSPs, one (PSR J1709–0333) appears to be isolated, while three (PSRs J1623–6939, J1823–3543, and J1858–5422) have companion masses and orbital periods that are consistent with helium white dwarf (He–WD) companions, but the remaining two warrant further mention. The first, PSR J1906–1754 lies at an unusual position in the orbital parameter space shown in Fig. 4, with a minimum companion mass ($M_c > 0.05 M_{\odot}$) much lower than that of typical white dwarf companions, but with an orbital period ($P_{orb} = 6.5$ d) much longer than a typical black-widow system. Of the MSP binaries with $P_{orb} > 1$ d, only PSR J1737–0811 has a similarly low minimum companion mass ($M_c > 0.06 M_{\odot}$, Boyles et al. 2013). The low minimum mass could indicate a very face-on binary inclination angle, with $i < 14.5$ deg required for $M_c = 0.2 M_{\odot}$. However, such orbits are a priori unlikely, as

only around 3 per cent of orbits should have a lower inclination, assuming a random distribution of viewing angles. The final pulsar, PSR J1757–6032, has a larger minimum companion mass ($M_c > 0.45 M_{\odot}$) that suggests it is perhaps more likely to have a CO–WD companion, similar to that of the relativistic binary PSR J1614–2230 (Demorest et al. 2010; Tauris, Langer & Kramer 2011).

## 3.4 Sensitivity

As our observations were not flux calibrated, we estimated the flux density thresholds ($S$) above which a pulsar should have been detectable to our survey. To estimate this, we used the pulsar radiometer equation (Lorimer & Kramer 2004),

$$S = \frac{\rho(T_{\text{sys}} + T_{\text{sky}})}{\beta G n_{\text{ant}} \sqrt{n_{\text{pol}} T_{\text{obs}} B}} \sqrt{\frac{w}{1-w}}, \quad (2)$$

where we have assumed: $G = 0.044$ K Jy$^{-1}$ is the gain per antenna (Bailes et al. 2020); $\rho = 9$ was the S/N threshold required for a candidate to be folded for visual inspection; $n_{\text{ant}}$ is the number of antennas; $n_{\text{pol}} = 2$ polarizations are summed; $T_{\text{obs}}$ is the exposure time; $B = 700$ MHz is the estimated useable bandwidth after RFI excision; and $T_{\text{sys}} = 22.5$ K is the combination of the receiver temperature, atmospheric, and ground spillover contributions (Ridolfi et al. 2021). The sky temperature, $T_{\text{sky}}$, towards each source is taken from the 408 MHz all-sky map of Haslam et al. (1982), reprocessed by Remazeilles et al. (2015), and scaled to the central frequency of 1284 MHz according to an assumed spectral index of $-2.6$. The resulting flux-density limits are given in Table 1. Estimated sensitivity also depends strongly on the assumed pulse duty cycle, $w$. For the new MSPs, this varies from 0.05 (PSR J1803–6707) to 0.35 (PSR J1906–1754); we have assumed $w = 0.15$ for our estimated sensitivities.

The remaining factor in equation (2), $\beta$, accounts for various losses incurred during observing and searching. This includes constant fractional losses due to the 8-bit voltage digitization (5 per cent) and beamforming efficiency (5 per cent); sensitivity losses due to using incoherent harmonic summing (15 per cent for our assumed 15 per cent duty cycle, Morello et al. 2020); finite time resolution and intrachannel dispersion smearing (both <1 per cent). In addition to these fractional losses, there are random losses due to the location of a signal between FFT bins (averaging 8 per cent), DM trials and acceleration trials. We estimate these losses for an assumed $v = 500$ Hz, DM = 100 pc cm$^{-3}$ pulsar, via a Monte Carlo procedure – generating signals with random offsets from the nearest search trial and evaluating the S/N losses. The final loss factor to consider is that due to the angular offset between a pulsar and the centre of the nearest coherent beam. To evaluate this, for each source we drew random locations from the 4FGL localization probability densities and evaluated the sensitivity at that location of the nearest coherent beam according to the simulated PSF model from Mosaic. Combining all of these losses provided a Monte Carlo distribution for $\beta$.

For each source we quote three estimated sensitivities: $S_{\text{opt}}$ is the optimum sensitivity, assuming only the constant losses described above, and therefore estimates a fundamental flux density limit for our survey; while $S_{50}$ and $S_{95}$ are the flux density limits obtained using equation (2) using the median and 95th centile values, respectively, of the $\beta$ values obtained from our Monte Carlo estimates. Average sensitivities were $S_{\text{opt}} \approx 55$ µJy, $S_{50} = 120 (90)$ µJy, and $S_{95} = 150 (110)$ µJy for the first (second) pass, respectively. The flux density threshold for the incoherent beam (for pulsars lying outside
the region tiled by coherent beams) is approximately $\sqrt{n_{\text{ant}}}$ ≈ 8 larger
than that for a coherent beam.

4 DISCUSSION

In this paper, we have presented the first MeerKAT survey for new
radio pulsars in unassociated *Fermi*-LAT sources. While our strategy
of targeting pulsar-like gamma-ray sources is certainly not novel,
the capabilities of this next-generation radio telescope do lead to
significant advantages over previous surveys.

Foremost among these is the extremely high sensitivity of the full MeerKAT array, surpassed only by the Arecibo (Cromartie et al. 2016) and FAST (Wang et al. 2021) telescopes, but unprecedented in the Southern hemisphere. In Section 3.4, we estimated typical 95 per cent flux density upper limits of ~100 µJy, which can be compared to the nominal ~200 µJy sensitivity (not including many of the loss terms that we consider) that was achieved with typical hour-long pointings in the similar survey of Fermi-LAT sources performed at Parkes by Camilo et al. (2015).

Of the newly discovered pulsars, only PSR J1803−6707 is
bright enough that it could perhaps have been discovered in hour-
long Parkes observations, but such a discovery would have been
complicated by jerk effects from its short orbital period and high
acceleration. All of the new pulsars have eventually been detected in
dedicated Parkes observations, but with low S/Ns that would be hard
to detect without prior knowledge of the DM and spin period.

A key factor for our survey, as Cromartie et al. (2016) also discussed for their survey of *Fermi*-LAT sources performed with the Arecibo telescope, is that these competitive flux-density thresholds are achieved with only short 10-min observations, bringing several benefits. Firstly, these short observation times allow for more sources to be observed within a given observing time budget, while MeerKAT’s rapid (~30 s) slew time ensures that costly overheads are not incurred by doing so. Secondly, short observations enable the search processing to be performed quickly (the computing cost of an acceleration search scales with at least $T^3 \log(T)$ for integration time $T$), allowing storage space to be freed quickly enough that observation scheduling is not limited by this factor, and allowing us to search up to high accelerations. Thirdly, short observation times enable the detection of very short-period binary MSPs. For very short orbits, the assumptions that go into an acceleration search (that the orbital motion within the observation can be approximated with a constant acceleration) break down as higher order ‘jerk’ terms become significant, with sensitivity to binary MSPs only maintained for observations lasting less than 10 per cent of an orbit (Johnston & Kulkarni 1991). Our short 10-min observation strategy mitigates the worst of the jerk effects, but still means that sensitivity is lost for binaries with periods shorter than ~100 min. Only three Galactic field MSPs with shorter orbital periods than this are known, although all three are gamma-ray MSPs (Pletsch et al. 2012; Stovall et al. 2014; Nieder et al. 2020b).

Additionally, short observations mean that multiple passes can be
performed to minimize missed discoveries due to the time-varying
effects that can contribute to a pulsar's detectability in a given
observation. Interstellar scintillation introduces time- and frequency-
dependent variations in the observed spectrum, which can lead to a
pulsar being undetectable over large frequency ranges (Camilo et al.
2015). Indeed PSR J1709−0333 remained undetected in our L-band
confirmation observation, despite a dense beam tiling that covered its
now-known position. Our UHF confirmation observations reveal the
reason for this – in both observations, separated by 70 min, the pulsar
was only detected in ~20 per cent (~100 MHz) of the bandwidth,

indicative of scintillation. Subsequent detections with the Parkes
UWL and the Effelsberg L-band receivers revealed similar behaviour.

MSPs in black-widow and redback binaries can also be unde-
tectable for large fractions of an orbit, usually (but not exclusively)
around the pulsar's superior conjunction, as a result of diffuse intra-
binary plasma dispersing, scattering and absorbing radio pulsations.
The jerk terms that limit the sensitivity of an acceleration search to
short period binaries also have an orbital phase dependence (Johnston
& Kulkarni 1991), and so sensitivity also depends on the orbital
phase at which observations take place. Our short exposure, two-pass
survey partially mitigates these effects. The case of PSR J1526−2744
further illustrates this: we were unable to detect this pulsar in
dedicated follow-up observations with Parkes and Nançay prior to
obtaining an orbital solution without performing jerk searches, a
computationally expensive technique (Andersen & Ransom 2018)
that has only recently started being employed in radio surveys of
Fermi sources.

One additional time-dependent effect that is specific to this survey
is the fact that the shape of a coherent tied-array beam on the
sky depends on the elevation at which a source is observed, with
beams being more elongated for sources at lower elevations due to
the smaller projected baselines between antennas. The beam tiling
pattern used to cover a given source region therefore depends on the
exact sidereal time and array configuration.

Indeed, four of the nine pulsars discovered here were only detected
in the second observation of their respective *Fermi* sources. This
was partially due to the improved flexibility of FBFUSE’s tiling
patterns that was developed between the first and second passes,
and our use of a coarser frequency resolution but larger number
of coherent beams, which allowed us to cover a larger solid angle
around each source with a more sensitive tiling pattern. Two of these
pulsars, PSRs J1709−0333 and J1036−4353, lay just outside the
corresponding *Fermi*-LAT source 95 per cent localization regions
that we aimed to tile with coherent beams in the first pass, but were
detected in the second pass where we targeted a larger nominal
99 per cent confidence region. The other two pulsars that were only
detected in the second pass, PSRs J1623−6936 and J1757−6032,
lay in a less sensitive location in the coherent beam tiling pattern
in their first observation, mid-way between three neighbouring
beams, where sensitivity was ~50 per cent of that at the centre
of a beam, but were at a more favourable location in the second
observation.

One final advantage to our survey, over projects using single-
dish telescopes, is the rapid and precise localization that can be
obtained from a multibeam detection using *SeeKAT*. This has two
significant scientific benefits: rapid localizations enable immediate
multiwavelength follow-up and catalogue searches; and precise
knowledge of a pulsar’s location greatly assists in obtaining a phase-
connected timing solution to fully exploit the scientific potential of a
new MSP discovery. In the absence of an interferometric localization,
it often takes lengthy timing campaigns to reveal the location of a
pulsar, as astrometric parameters can be highly degenerate with a
pulsar’s spin-down rate (and sometimes orbital parameters) until a
data set spanning several months has been obtained. The *SeeKAT*
positions for PSRs J1803−6707 and J1036−4353, were precise
enough to unambiguously link these pulsars to their *Gaia* coun-
terparts, establishing these as redback binary systems, and providing
sub-millisecond astrometric uncertainties that are sufficiently
precise that these parameters no longer needed to be fit for when
building a timing solution. For PSR J1526−2744, the ~2 arcsec
positional uncertainty represents a reduction of a factor of ~10 000
in the number of sky positions that had to be searched to detect
gamma-ray pulsations, greatly decreasing the required computational effort.

With only two pointings towards our targeted sources there is still a high chance that some otherwise detectable pulsars among our targets may have been missed due to scintillation or eclipsing. We have therefore planned another two observations towards the sources observed here using the UHF receiver. All 9 pulsars discovered here were also detected at a higher S/N in our UHF observations, illustrating the sensitivity gains that can be made by observing at lower frequencies to exploit typical MSP spectra that decrease with frequency (Frail et al. 2016; Jankowski et al. 2018), and MeerKAT's extremely low RFI environment at UHF (Bailes et al. 2020).

It is important to consider what prospects remain for further detections among Fermi-LAT sources. To obtain a crude estimate for the number of pulsars remaining in this sample, we can simply sum the probabilities of each source being a gamma-ray pulsar that the Random Forest classifier computed in Section 2.4, obtaining $\sum P(\text{psr}) = 38.4$. Subtracting from this sum the 9 new pulsars discovered here, plus the 3 others mentioned in Section 3 that were recently discovered by other projects, gives ~26 pulsars remaining in this sample, and around 190 pulsars in the full 4FGL catalogue. However, this is certainly an overestimate, as many of these sources have been searched already with other telescopes, and therefore our target list represents the remainder from a larger, unbiased sample from which more easily detected pulsars have already been removed. There are hints of this effect in our results: five of the nine new pulsars had fairly low pulsar probabilities $P(\text{psr}) < 0.25$, perhaps because the more promising sources have already been surveyed extensively by other telescopes, while the lower probability targets tend to be fainter, newer sources many of which have not been searched before. This also suggests that further searches of a larger number of apparently less promising Fermi-LAT sources below the $P(\text{psr})$ threshold that we used here may yet bear fruit.

Furthermore, even if this is an accurate estimate of the number of gamma-ray pulsars in our sample, it does not mean that there are 26 detectable radio pulsars in this sample. A large fraction of the young gamma-ray pulsars detected by the Fermi-LAT remain undetected in radio observations (and indeed the two recently discovered young gamma-ray pulsars that were covered in our survey were not detected), presumably due to the viewing angle missing the radio beam. While we deliberately avoided sources at low Galactic latitudes in order to reduce the chance of young pulsars entering our sample, we note that there are 17 sources in our list at lower Galactic latitudes than 4FGL J0802.1-5612, one of the two recently discovered young pulsars in our sample. Although MSPs tend to have wider radio beams, and hence are less likely to be 'radio-quiet', a handful of gamma-ray MSPs have now been discovered in the Fermi-LAT data, but have remained undetected in deep radio follow-up searches (Clark et al. 2018; Nieder et al. 2020b), and similar objects may exist within our targets.

Nevertheless, it is certain that there are still several detectable MSPs lurking in our target list, and iterative 4FGL releases (Ballet et al. 2020; Abdollahi et al. 2022) bring new, albeit fainter, unassociated sources to target. Indeed, since our survey began, likely MSP binaries have been discovered in optical and X-ray searches of three of our target sources, 4FGL J0540.0-7552 (Strader et al. 2021), 4FGL J0940.3-7610 (Swihart et al. 2021), and 4FGL J1120.0-2204 (Swihart et al. 2022), while a fourth, 4FGL J1702.7-5655, has recently been identified through the detection of gamma-ray eclipses (Corbet et al. 2022). 4FGL J1120.0-2204 in particular appears to have a white dwarf companion in a fairly long (15 h) orbit, but has remained undetected in many previous long radio observations,

as well as in our 10-min pointings. Its non-detection is therefore unlikely to be explained by orbital or eclipsing effects, suggesting that deeper searches, in addition to our continuing 'shallow' survey, may be necessary to uncover the nature of many of the remaining unidentified but pulsar-like Fermi-LAT sources.

## 5 CONCLUSIONS AND FUTURE WORK

We have presented the first results from TRAPUM's survey for new pulsars in unassociated Fermi-LAT sources using the MeerKAT radio telescope, discovering nine new MSPs, of which eight are in binary systems. Our results continue the trend of short-period binaries being discovered at a far higher rate within Fermi-LAT sources than in untargeted surveys, with the discovery of two new redback binaries (PSRs J1036-4353 and J1803-6707) with optically bright companion stars and radio eclipses, and a third (PSR J1526-2744) that is possibly the most compact known MSP-WD binary system in the Galactic field. Two other MSPs have preliminary orbital solutions that mark their companions as outliers among the known MSP-binary population, PSR J1757-6032 appears to have a less typical CO-WD companion, while PSR J1906-1754 has a much longer orbital period than any known black widow, but a minimum companion mass that is far lighter than that of a typical WD companion.

We obtained phase-connected timing solutions for two of the new MSPs, PSRs J1526-2744 and J1803-6707, using timing observations at Parkes and Nançay, that enabled us to detect gamma-ray pulsations in the Fermi-LAT data. Orbital period variations in PSR J1803-6707 precluded further extrapolation of the timing solution, but for PSR J1526-2744, gamma-ray timing provides a full 14-yr timing solution. This in turn enabled us to search for continuous gravitational waves from this pulsar in the Advanced LIGO O1, O2, and O3 data, although none were detected and the strain upper limit remains well above the pulsar's spin-down luminosity budget. A dedicated timing campaign is underway at MeerKAT to obtain phase-connected timing solutions, and afterwards gamma-ray pulsations, from the remaining seven discoveries (Burgay et al., in preparation), as are optical follow-up observations and modelling of the new redback binaries (Phosrisom et al., in preparation).

Our results also emphasize the promise of continued radio surveys of Fermi-LAT sources: several high-confidence pulsar candidates still remain within our target list but have eluded detection in ours and earlier surveys, perhaps due to scintillation, eclipses, or simple intrinsic faintness. Indeed, our TRAPUM survey continues to make discoveries beyond the first results presented here. Two further passes of the sources surveyed here are currently underway using MeerKAT's UHF receiver, from which we have discovered 8 additional MSPs already.<fnref n="12" /> We have also performed several hour-long observations of a group of high confidence redback candidates that were previously identified in optical observations, detecting a further three new MSPs that will be presented in a dedicated paper (Thongmeearkom et al., in preparation). We will continue our initial survey strategy of observing a large number of sources with short (10-min) observations, expanding our target list to include more sources from the most recent 4FGL DR3 catalogue, but deeper observations may be required to fully explore the most promising gamma-ray sources.

<footnote n="12">¹²http://trapum.org/discoveries/</footnote>
## ACKNOWLEDGEMENTS

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. We thank staff at SARAO for their help with observations and commissioning. TRAPUM observations used the FBFUSE and AP-SUSE computing clusters for data acquisition, storage and analysis. These clusters were funded and installed by the Max-Planck-Institut für Radioastronomie (MPIfR) and the Max-Planck-Gesellschaft. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Parkes radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42) which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Wiradjuri people as the traditional owners of the Observatory site. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS) and Université d'Orléans. It is partially supported by the Region Centre Val de Loire in France. Partly based on observations with the 100-m telescope of the MPIfR at Effelsberg.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Etudes Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

eROSITA is the primary instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from IKI and the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nürnberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the ESASS/NRTA software system developed by the German eROSITA consortium.

This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO Laboratory and Advanced

LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. The construction and operation of KAGRA are funded by Ministry of Education, Culture, Sports, Science and Technology (MEXT), and Japan Society for the Promotion of Science (JSPS), National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea, Academia Sinica (AS) and the Ministry of Science and Technology (MoST) in Taiwan.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Based on observations collected at the European Southern Observatory under ESO programmes 105.20RJ.001 and 105.20RJ.002.

CJC and RPB acknowledge support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 715051; Spiders). BWS and MCB acknowledge funding from the ERC under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 694745). MB, AP, and AR gratefully acknowledge financial support by the research grant 'iPeska' (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016. EDB, CJC, DJC, WC, MK, VVK, LN, PVP, and AR acknowledge continuing valuable support from the Max-Planck Society. LV acknowledges financial support from the Dean's Doctoral Scholar Award from the University of Manchester. VSD, TRM, and ULTRACAM acknowledge the support of the STFC. SMR is a CIFAR Fellow and is supported by the NSF Physics Frontiers Center awards 1430284 and 2020265.

We would like to thank Matthew Kerr, Guillem Martí-Devesa, and David Smith for reviewing this paper on behalf of the Fermi-LAT collaboration.

## DATA AVAILABILITY

TRAPUM data products are available upon reasonable request to the TRAPUM collaboration. The Fermi-LAT data are available from the Fermi Science Support Center (http://fermi.gsfc.nasa.gov/ssc).

## REFERENCES

* Aasi J. et al., 2014, ApJ, 785, 119

* Abbott B. P. et al., 2018, Phys. Rev. Lett., 120, 031104

* Abbott B. P. et al., 2019, Phys. Rev. D, 99, 122002

* Abbott R. et al., 2020, ApJ, 902, L21

* Abbott R. et al., 2021, SoftwareX, 13, 100658

* Abbott R. et al., 2022a, ApJ, 932, 133

* Abbott R. et al., 2022b, ApJ, 935, 1

* Abdo A. A. et al., 2013, ApJS, 208, 17

* Abdollahi S. et al., 2020, ApJS, 247, 33

* Abdollahi S. et al., 2022, ApJS, 260, 53
Ackermann M. et al., 2012, ApJ, 753, 83

Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982, Nature, 300, 728

Andersen B. C., Ransom S. M., 2018, ApJ, 863, L13

Antoniadis J. et al., 2013, Science, 340, 448

Ashok A. et al., 2021, ApJ, 923, 85

Atwood W. B. et al., 2009, ApJ, 697, 1071

Atwood W. et al., 2013, in Brandt T. J., Omodei N., Wilson-Hodge C. et al., Proc. 4th Fermi Symposium 2012. eConf C121028. Monterey, California, p. 8

Bailes M. et al., 2020, PASA, 37, e028

Ballet J., Burnett T. H., Digel S. W., Lott B., 2020, preprint (arXiv:2005.11208)

Barr E. D., 2018, in Weltevrede P., Perera B. B. P., Preston L. L., Sanidas S., eds, Pulsar Astrophysics the Next Fifty Years, Vol. 337. Jodrell Bank Observatory, UK, p. 175

Barr E., 2020, Astrophysics Source Code Library, record ascl:2001.014

Barr E. D. et al., 2013, MNRAS, 429, 1633

Bassa C. G. et al., 2017, ApJ, 846, L20

Becker W., Truemper J., 1997, A&A, 326, 682

Berteaud J., Calore F., Clavel M., Serpico P. D., Dubus G., Petrucci P.-O., 2021, Phys. Rev. D, 104, 043007

Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203, 1

Bhattacharyya B. et al., 2021, ApJ, 910, 160

Boyles J. et al., 2013, ApJ, 763, 80

Breiman L., 2001, Mach. Learn., 45, 5

Bruel P., 2019, A&A, 622, A108

Bruel P., Burnett T. H., Digel S. W., Johannesson G., Omodei N., Wood M., 2018, preprint (arXiv:1810.11394)

Camilo F. et al., 2015, ApJ, 810, 85

Carli E. et al., 2022, MNRAS, 517, 5406

Chen W., Barr E., Karuppusamy R., Kramer M., Stappers B., 2021, J. Astron. Instrum., 10, 2150013

Clark C. J. et al., 2018, Sci. Adv., 4, eaao7228

Cognard I. et al., 2011, ApJ, 732, 47

Corbet R. H. D. et al., 2022, ApJ, 935, 2

Cromartie H. T. et al., 2016, ApJ, 819, 34

Cutler C., Schutz B. F., 2005, Phys. Rev. D, 72, 063006

Demorest P. B., Pennucci T., Ransom S. M., Roberts M. S. E., Hessels J. W. T., 2010, Nature, 467, 1081

Deneva J. S. et al., 2016, ApJ, 823, 105

de Jager O. C., Raubenheimer B. C., Swanepoel J. W. H., 1989, A&A, 221, 180

Dhillon V. S. et al., 2007, MNRAS, 378, 825

Finke T., Krämer M., Manconi S., 2021, MNRAS, 507, 4061

Frail D. A., Jagannathan P., Mooley K. P., Intema H. T., 2016, ApJ, 829, 119

Freire P. C. C., Ridolfi A., 2018, MNRAS, 476, 4794

Gaia Collaboration, 2016, A&A, 595, A1

Gaia Collaboration, 2022, preprint (arXiv:2208.00211)

Guillemot L. et al., 2013, ApJ, 768, 169

Haslam C. G. T., Salter C. J., Stoffel H., Wilson W. E., 1982, A&AS, 47, 1

HI4PI Collaboration, 2016, A&A, 594, A116

Hobbs G. et al., 2020, PASA, 37, e012

Jankowski F., van Straten W., Keane E. F., Bailes M., Barr E. D., Johnston S., Kerr M., 2018, MNRAS, 473, 4436

Johnston H. M., Kulkarni S. R., 1991, ApJ, 368, 504

Jonas J. L., 2009, IEEE Proc., 97, 1522

Jonas J., MeerKAT Team, 2016, in Taylor R., Camilo F., Leeuw L., Moodley K., eds, Proc. MeerKAT Science: On the Pathway to the SKA – PoS(MeerKAT2016). SISSA, Trieste, PoS#1

Kaspi V. et al., 2000, ApJ, 543, 321

Keith M. J. et al., 2011, MNRAS, 414, 1292

Kerr M., 2011, ApJ, 732, 38

Kerr M. et al., 2012, ApJ, 748, L2

Kramer M. et al., 2016, in Taylor R., Camilo F., Leeuw L., Moodley K., eds, Proc. MeerKAT Science: On the Pathway to the SKA – PoS(MeerKAT2016). SISSA, Trieste, PoS#3

Krautter J. et al., 1999, A&A, 350, 743

Lee K. J., Guillemot L., Yue Y. L., Kramer M., Champion D. J., 2012, MNRAS, 424, 2832

Lorimer D. R., Kramer M., 2004, Handbook of Pulsar Astronomy, Vol. 4. Cambridge Univ. Press, Cambridge

Luo S., Leung A. P., Hui C. Y., Li K. L., 2020, MNRAS, 492, 5377

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

Marelli M., De Luca A., Caraveo P. A., 2011, ApJ, 733, 82

Morello V. et al., 2019, MNRAS, 483, 3673

Morello V., Barr E. D., Stappers B. W., Keane E. F., Lyne A. G., 2020, MNRAS, 497, 4654

Nieder L. et al., 2019, ApJ, 883, 42

Nieder L., Allen B., Clark C. J., Pletsch H. J., 2020a, ApJ, 901, 156

Nieder L. et al., 2020b, ApJ, 902, L46

Nieder L., Kerr M., Clark C. J., Bruel P., Cromartie H. T., Ransom S. M., Ray P. S., 2022, ApJ, 931, L3

Obrocka M., Stappers B., Wilkinson P., 2015, A&A, 579, A69

Pallanca C., Ransom S. M., Ferraro F. R., Dalessandro E., Lanzoni B., Hessels J. W. T., Stairs I., Freire P. C. C., 2014, ApJ, 795, 29

Pletsch H. J. et al., 2012, Science, 338, 1314

Pleunis Z. et al., 2017, ApJ, 846, L19

Predehl P. et al., 2021, A&A, 647, A1

Ransom S. M. et al., 2011, ApJ, 727, L16

Ray P. S. et al., 2012, in Morselli A., ed., Proc. 2011 Fermi Symposium, eConf C110509

Ray P. S. et al., 2022, ApJ, 927, 216

Remazeilles M., Dickinson C., Banday A. J., Bigot-Sazy M. A., Ghosh T., 2015, MNRAS, 451, 4311

Ridolfi A. et al., 2016, MNRAS, 462, 2918

Ridolfi A. et al., 2021, MNRAS, 504, 1407

Roberts M. S. E., 2012, in van Leeuwen J., ed., Proc. IAU Symp. S291: Neutron Stars and Pulsars: Challenges and Opportunities After 80 Years. Cambridge Univ. Press, Cambridge, p. 127

Roberts M. S. E., McLaughlin M. A., Gentile P. A., Ray P. S., Ransom S. M., Hessels J. W. T., 2015, in Fukazawa Y., Tanaka Y., Itoh R., eds, Proc. 5th Fermi Symposium, eConf C14102.1. p. 62

Romani R. W., Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2022, ApJ, 934, L17

Saz Parkinson P. M., Xu H., Yu P. L. H., Salvetti D., Marelli M., Falcone A. D., 2016, ApJ, 820, 8

Schroeder J., Halpern J., 2014, ApJ, 793, 78

Smarr L. L., Blandford R., 1976, ApJ, 207, 574

Spiewak R. et al., 2022, PASA, 39, e027

Stappers B., Kramer M., 2016, in Taylor R., Camilo F., Leeuw L., Moodley K., eds, Proc. MeerKAT Science: On the Pathway to the SKA – PoS(MeerKAT 2016). SISSA, Trieste, PoS#9

Steltner B., Papa M. A., Eggenstein H.-B., 2022, Phys. Rev. D, 105, 022005

Stovall K. et al., 2014, ApJ, 791, 67

Strader J. et al., 2021, ApJ, 917, 69

Swihart S. J., Strader J., Aydi E., Chomiuk L., Dage K. C., Shishkovsky L., 2021, ApJ, 909, 185

Swihart S. J., Strader J., Aydi E., Chomiuk L., Dage K. C., Kawash A., Sokolovsky K. V., Ferrara E. C., 2022, ApJ, 926, 201

Tauris T. M., Langer N., Kramer M., 2011, MNRAS, 416, 2130

The Fermi-LAT Collaboration, 2022, Science, 376, 521

Wang P. et al., 2021, Sci. China Phys. Mech. Astron., 64, 129562

Woan G., Pitkin M. D., Haskell B., Jones D. I., Lasky P. D., 2018, ApJ, 863, L40

Wu J. et al., 2018, ApJ, 854, 99

Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29

Zhu W. W. et al., 2014, ApJ, 781, 117

<footnote n="1">¹Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hannover, Germany</footnote>

<footnote n="2">²Leibniz Universität Hannover, D-30167 Hannover, Germany</footnote>

<footnote n="3">³Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="4">⁴Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>
⁵INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius, CA, Italy

⁶South African Radio Astronomy Observatory, 2 Fir Street, Black River Park, Observatory 7925, South Africa

⁷Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstraße, D-85748 Garching, Germany

⁸Laboratoire d'Annecy-le-Vieux de Physique Théorique (LAPTh), CNRS, USMB, F-74940 Annecy, France

⁹LPC2E – Université d'Orléans / CNRS, F-45071 Orléans cedex 2, France

¹⁰Observatoire Radioastronomique de Nançay (ORN), Observatoire de Paris, Université PSL, Univ Orléans, CNRS, F-18330 Nançay, France

¹¹National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

¹²SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK

¹³Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa

¹⁴School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 E8C0, Ireland

¹⁵Fakultät für Physik, Universität Bielefeld, Postfach 100131, D-33501 Bielefeld, Germany

¹⁶Department of Physics and Astronomy, University of Sheffield, Sheffield S3 7RH, UK

¹⁷Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

¹⁸Department of Physics, University of Warwick, Coventry CV4 7AL, UK

This paper has been typeset from a LaTeX file prepared by the author.