# Radio and gamma-ray timing of TRAPUM L-band Fermi pulsar survey discoveries

M Burgay, L Nieder, C.J Clark, P.C.C Freire, S Buchner, T Thongmeearkom, J.D Turner, E Carli, Ismaël Cognard, Jean-Mathias Griessmeier, et al.

► To cite this version:

M Burgay, L Nieder, C.J Clark, P.C.C Freire, S Buchner, et al. Radio and gamma-ray timing of TRAPUM L-band Fermi pulsar survey discoveries. Astronomy & Astrophysics - A&A, 2024, 691, pp.A315. 10.1051/0004-6361/202451530. hal-04805239

HAL Id: hal-04805239

https://hal.science/hal-04805239v1

Submitted on 28 Nov 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radio and gamma-ray timing of TRAPUM L-band Fermi pulsar survey discoveries

M. Burgay<sup>1,*</sup> ◦, L. Nieder<sup>2,3</sup>, C. J. Clark<sup>2,3</sup> ◦, P. C. C. Freire<sup>4</sup> ◦, S. Buchner<sup>5</sup>, T. Thongmeearkom<sup>6,7</sup> ◦,
J. D. Turner<sup>6</sup> ◦, E. Carli<sup>6</sup> ◦, I. Cognard<sup>8,9</sup>, J.-M. Grießmeier<sup>8,9</sup> ◦, R. Karuppusamy<sup>4</sup>, M. C. i Bernadich<sup>4,1</sup> ◦,
A. Possenti<sup>1</sup> ◦, V. Venkatraman Krishnan<sup>4</sup> ◦, R. P. Breton<sup>6</sup> ◦, E. D. Barr<sup>4</sup>, B. W. Stappers<sup>6</sup>, M. Kramer<sup>4,6</sup>,
L. Levin<sup>6</sup> ◦, S. M. Ransom<sup>10</sup>, and P. V. Padmanabh<sup>2,3</sup> ◦

Received 16 July 2024 / Accepted 26 September 2024

## ABSTRACT

This paper presents the results of a joint radio and gamma-ray timing campaign on the nine millisecond pulsars (MSPs) discovered as part of the L-band targeted survey of Fermi-LAT sources performed in the context of the Transients and Pulsars with MeerKAT (TRAPUM) Large Survey Project. Out of these pulsars, eight are members of binary systems; of these eight, two exhibit extended eclipses of the radio emission. Using an initial radio timing solution, pulsations were found in the gamma rays for six of the targets. For these sources, a joint timing analysis of radio times of arrival and gamma-ray photons was performed, using a newly developed code that optimises the parameters through a Markov chain Monte Carlo (MCMC) technique. This approach has allowed us to precisely measure both the short- and long-term timing parameters. This study includes a proper motion measurement for four pulsars, which a gamma ray-only analysis would not have been sensitive to, despite the 15-year span of Fermi data.

**Key words.** binaries: eclipsing – binaries: general – pulsars: general – gamma rays: stars

## 1. Introduction

Pulsars and, in particular, the fastest rotating millisecond pulsars (MSPs), are excellent laboratories for astrophysics and fundamental physics. Thanks to their clock-like nature, for instance, they can be used to test the predictions of relativistic gravity theories (e.g. Freire & Wex 2024; Kramer et al. 2021), or as detectors of gravitational waves in the nano-hertz regime (e.g. Hobbs et al. 2010; Agazie et al. 2023b; EPTA Collaboration & InPTA Collaboration 2023; Reardon et al. 2023; FERMI-LAT Collaboration 2022), or as probes of the equation of state (EoS) of nuclear matter (e.g. Hu et al. 2020) and other phenomena (Lorimer & Kramer 2004).

Over one quarter of the presently known population of Galactic MSPs (see the ATNF Pulsar Catalogue PSRCat<fnref n="1" />, Manchester et al. 2005, and the GalacticMSP catalogue<fnref n="2" />, including unpublished discoveries) have been discovered targeting unassociated gamma-ray point sources detected by the Large Area Telescope (LAT) instrument on board the Fermi satellite (Smith et al. 2023, hereafter 3PC). These include some of the

most extreme and/or elusive pulsars, such as the second tightest orbit binary MSP J1653−0158 (Nieder et al. 2020b), which was only recently surpassed by J1953+1844, discovered in the M71 globular cluster (Pan et al. 2023); the fastest spinning MSP in the Galactic field (Bassa et al. 2017); and over two-thirds of all Galactic eclipsing binaries, namely, the so-called ‘spider’ pulsars (Roberts 2012). The latter are categorised as either ‘black widows’ (BWs), with short duration eclipses and companion stars with masses below 0.1 $M_\odot$, or ‘redbacks’ (RBs), with eclipses covering large fractions of the orbital period and larger companion masses. Some RBs, namely, the so-called transitional MSPs, represent the missing link between accreting low-mass X-ray binaries and fully recycled MSPs (e.g. Archibald et al. 2009).

The full scientific potential of MSPs as science probes can only be extracted through long-term follow-up campaigns aimed at finding a ‘coherent timing solution’ that would provide very-high-accuracy astrometric, spin and (in the case of binary systems) orbital parameters. Typically, at least one year is needed to get an accurate position and ensure no covariance with the spin-period derivative measurement (a crucial parameter in estimating pulsar energetics). Typically, several years are needed to get a proper motion measurement.

<footnote n="1">¹ INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius, (CA), Italy</footnote>

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

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

<footnote n="4">⁴ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="5">⁵ South African Radio Astronomy Observatory, 2 Fir Street, Black River Park Observatory, Cape Town 7925, South Africa</footnote>

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

<footnote n="7">⁷ National Astronomical Research Institute of Thailand, Don Kaeo, Mae Rim, Chiang Mai 50180, Thailand</footnote>

<footnote n="8">⁸ LPC2E – Université d'Orléans/CNRS, 45071 Orléans cedex 2, France</footnote>

<footnote n="9">⁹ Observatoire Radioastronomique de Nançay (ORN), Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France</footnote>

<footnote n="10">¹⁰ National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA</footnote>

<footnote>* Corresponding author; marta.burgay@inaf.it</footnote>

<footnote n="1">¹ https://www.atnf.csiro.au/people/pulsar/psrcat/</footnote>

<footnote n="2">² http://astro.phys.wvu.edu/GalacticMSPs/</footnote>
<table><caption>**Table 1.** Observational parameters for the timing observations.</caption><thead><tr><th>Telescope</th><th>Freq (GHz)</th><th>BW (MHz)</th><th>No. chan</th><th>T<sub>samp</sub> (µs)</th><th>Npol</th><th>Nbit</th><th>Backend</th><th>Dedisp</th></tr></thead><tbody><tr><td rowspan="2">MeerKAT</td><td>1.3</td><td>856</td><td>1024</td><td>9.6</td><td>4</td><td>8</td><td>PTUSE</td><td>Coh</td></tr><tr><td>0.8</td><td>544</td><td>1024</td><td>7.5</td><td>4</td><td>8</td><td>PTUSE</td><td>Coh</td></tr><tr><td>Parkes</td><td>2.4</td><td>3328</td><td>3328</td><td>64</td><td>1</td><td>8</td><td>Medusa</td><td>Coh</td></tr><tr><td>Nançay</td><td>1.5</td><td>512</td><td>1024</td><td>64</td><td>1</td><td>4</td><td>Nuppi</td><td>Inc</td></tr><tr><td rowspan="2">Effelsberg</td><td>1.5</td><td>250</td><td>1024</td><td>82</td><td>1</td><td>32</td><td>PSRIX</td><td>Inc</td></tr><tr><td>1.5</td><td>250</td><td>512</td><td>51</td><td>1</td><td>8</td><td>EDD</td><td>Inc</td></tr></tbody></table>

**Notes.** Left to right the columns report: telescope used, central observing frequency, useable observing bandwidth, number of frequency channels, sampling time, number of polarisations (1=total intensity, 4=full Stokes), number of digitising bins, backend used and de-dispersion type (coherent or incoherent)

In the case of Fermi-selected targets, once an initial timing solution is obtained from multiple radio observations, if high-energy pulsations are found and the signal-to-noise ratio (S/N) is high enough, gamma-ray photons can be used to extend the timing solution to cover the total time span of the Fermi mission (15 years at the time of writing). This enables the immediate measurement of timing parameters that require longer time spans, minimising the required radio follow-up time. Timing solutions covering long time spans also have the advantage of making multi-wavelength and multi-messenger searches for pulsar counterparts easier, as they may provide timing solutions valid for the epoch of other instruments' observations (including e.g. the long time spans covered by the LIGO and VIRGO detectors for gravitational waves; see e.g. Abbott et al. 2017). Besides the advantages for timing purposes, the availability of information on both radio and gamma-ray pulsations allows for a comparison between emission geometries (and eventually emission mechanisms) in these two very different energy ranges (e.g. Venter et al. 2012; Pétri & Mitra 2021).

In this paper we describe a multi-telescope and multi-wavelength timing campaign to precisely characterise the pulsar discoveries made as part of the TRAPUM shallow L-band survey of Fermi-LAT unassociated point sources (Clark et al. 2023). The paper is organised as follows. In Sect. 2, we describe the multi-telescope radio timing campaign. In Sect. 3 we report the gamma-ray data, searches, and discoveries. In Sect. 4, we describe a new method to jointly time pulse times of arrival obtained from the radio observations and gamma-ray photons. In Sect. 5, we present the multi-wavelength timing results and a brief comparison of radio and gamma-ray properties. In Sect. 6, our basic polarimetric results are shown. Finally, in Sect. 7, we present our summary and conclusions.

## 2. Radio timing campaign

Some information on the initial phases of the timing campaign on discoveries made during TRAPUM's L-band survey of unassociated Fermi-LAT sources was already reported in our previous paper describing the survey (Clark et al. 2023). We presented an initial coherent timing solution for two pulsars, and basic spin, astrometric, and orbital parameters for the remaining seven. In the current work, we describe the full timing campaign in more detail.

In order to minimise telescope usage at MeerKAT, where no open time for time-domain experiments was available until September 2023, a multi-telescope timing campaign was set up to follow up TRAPUM discoveries starting in March 2021. For the brightest pulsars, with an S/N of the integrated pulse profile

above ~20 in the 10 minute discovery observation at MeerKAT, we obtained time at the Murriyang 64-m telescope (Parkes) in Australia. There we exploited the 3 gigahertz (GHz) bandwidth of the Ultra-Wideband Low receiver (UWL, Hobbs et al. 2020). For pulsars with a declination above ~ - 35 degrees, we also obtained time at the Effelsberg 100-m radio telescope (Germany) and the NRT/Nançay decimetric radio telescope (France). Depending on the spectral properties of each pulsar and on their exact sky location, this usually allowed us to get a good detection with one- to two-hour pointings.

Starting in November 2021, 53 hours were also granted to us to follow up TRAPUM Fermi discoveries at MeerKAT. This time was used to get initial parameters for our new pulsars, in order to better assess the suitability of less sensitive telescopes for their subsequent follow-up, and to obtain full timing solutions for the fainter and/or tighter orbit ones, for which the smearing of the pulse period due to binary motion over one to two hours would prohibit the detection of pulsations in an acceleration search.

Details on the observational set-up at each telescope (central frequency, bandwidth, frequency and time resolution, backend, de-dispersion type) are reported in Table 1.

Whenever possible, observations were executed with a pseudo-logarithmic cadence (a few observations on day 1 and single observations on days 2, 5, 10, etc.), followed by monthly observations, in order to efficiently achieve timing coherence.

The basic steps to obtain a timing solution for each pulsar are briefly described below. First, the position of the pulsar was refined through a triangulation method based on the comparison of the S/N of the integrated pulse profiles of different detections (obtained in multiple neighbouring coherent tied-array beams used in the survey's discovery observations), given a model of the beam point spread function. This method made use of the SeeKAT software described in Bezuidenhout et al. (2023) and allowed us to achieve positional accuracy typically better than 5 arcseconds. The improved position was used in all the follow-up timing observations and allowed us to promptly investigate multi-wavelength counterparts of our discoveries (e.g. Dodge et al. 2024 and Phosrisom et al. in prep.). For short-period binary pulsars, we used this interferometric position to search for possible optical counterparts in the Gaia DR3 catalogue (Gaia Collaboration 2016, 2023), in which case the Gaia astrometric solution could be used as a highly precise, independent prior for subsequent analyses.

Subsequently, for each observation, the most prominent radio frequency interference (RFI) signals were removed using the routine rfifind from the PRESTO<fnref n="3" /> software package

<footnote n="3">³ https://github.com/scottransom/presto</footnote>
<table><caption>**Table 2.** Main parameters of the radio timing campaign.</caption><thead><tr><th rowspan="3">PSR</th><th rowspan="3">Data span<br>(MJD)</th><th colspan="4">N. obs</th><th colspan="4">tobs<br>(h)</th><th colspan="4">N. ToA</th></tr><tr><th colspan="4"></th><th colspan="4"></th><th colspan="4"></th></tr><tr><th>MK</th><th>PKS</th><th>NRT</th><th>EFF</th><th>MK</th><th>PKS</th><th>NRT</th><th>EFF</th><th>MK</th><th>PKS</th><th>NRT</th><th>EFF</th></tr></thead><tbody><tr><td>J1036-4353</td><td>59536-60015</td><td>20</td><td>15</td><td>-</td><td>-</td><td>1.7</td><td>29.8</td><td>-</td><td>-</td><td>99</td><td>81</td><td>-</td><td>-</td></tr><tr><td>J1526-2744</td><td>59304-59432</td><td>2</td><td>6</td><td>7</td><td>-</td><td>0.3</td><td>9.3</td><td>6.3</td><td>-</td><td>16</td><td>28</td><td>27</td><td>-</td></tr><tr><td>J1623-6936</td><td>59250-59683</td><td>23</td><td>14</td><td>-</td><td>-</td><td>2.2</td><td>14.6</td><td>-</td><td>-</td><td>140</td><td>38</td><td>-</td><td>-</td></tr><tr><td>J1709-0333</td><td>59304-60084</td><td>20</td><td>4</td><td>-</td><td>2</td><td>1.8</td><td>6.1</td><td>-</td><td>3.4</td><td>111</td><td>32</td><td>-</td><td>8</td></tr><tr><td>J1757-6032</td><td>59250-59725</td><td>29</td><td>12</td><td>-</td><td>-</td><td>4.1</td><td>13.9</td><td>-</td><td>-</td><td>274</td><td>36</td><td>-</td><td>-</td></tr><tr><td>J1803-6707</td><td>59197-59624</td><td>2</td><td>19</td><td>-</td><td>-</td><td>1.4</td><td>11.0</td><td>-</td><td>-</td><td>56</td><td>148</td><td>-</td><td>-</td></tr><tr><td>J1823-3544</td><td>59250-59684</td><td>19</td><td>9</td><td>10</td><td>-</td><td>1.8</td><td>9.7</td><td>7.4</td><td>-</td><td>192</td><td>32</td><td>35</td><td>-</td></tr><tr><td>J1858-5422</td><td>59250-59759</td><td>37</td><td>6</td><td>-</td><td>-</td><td>6.0</td><td>6.1</td><td>-</td><td>-</td><td>265</td><td>10</td><td>-</td><td>-</td></tr><tr><td>J1906-1754</td><td>59250-60148</td><td>43</td><td>5</td><td>-</td><td>20</td><td>4.0</td><td>6.1</td><td>-</td><td>48.9</td><td>208</td><td>5</td><td>-</td><td>72</td></tr></tbody></table>

**Notes.** Left to right the columns report: Pulsar name, data span of the radio timing solution and, for each telescope (left to right: MeerKAT, Parkes, Nançay, Effelsberg), number, total duration of the observations used to obtain the timing solution and the number of ToAs. Observations with no detections are not included in the total counts. Discovery observations have been included, whenever it has been possible to use them for timing.

(Ransom 2011). The same software was used to de-disperse the cleaned data (with *prepdata*) and either perform a blind acceleration search on the resulting file (through *accelsearch*, Ransom et al. 2002), or, once a preliminary ephemeris was available, directly fold it (*prepfold*).

In the case of a binary pulsar, the initial orbital solution was found, using the first few observations, by fitting an ellipse in the acceleration-period plane (following Freire et al. 2001, especially useful in case of detections widely separated in time) and later by fitting a sinusoidal modulation to the barycentric spin periods (using *fit circular orbit*.py from PRESTO).

Multiple times of arrival (ToAs; typically from one to a dozen, depending on the S/N) were extracted from each observation using PRESTO’s *get_TOAs.py* by cross-correlating a profile template with each folded observation. The template initially used was the highest S/N profile for each pulsar. Once a coherent timing solution was achieved, multiple observations were added together in phase, resulting in profiles from which we obtained analytic, noiseless templates by fitting multiple Gaussian components.

The software *Tempo2*<fnref n="4" /> (Hobbs et al. 2006) was used to obtain a coherent timing solution by minimising the root mean square (rms) of the ToAs residuals in an iterative way, with increasingly more precise ephemerides being used to refine both the template and the ToAs (two iterations were used in most cases). Wherever needed, constant time offsets (or ‘jumps’) were also fitted, together with the pulsar parameters, to account for observatory-related time delays between the different telescopes, and/or observing systems.

In some cases, when the coverage of ToAs in time, and/or the orbital coverage, was too sparse, we achieved timing coherence using a new method. This algorithm utilises the software *Dracula*<fnref n="5" /> (Freire & Ridolfi 2018), is implemented through *Tempo* (Nice et al. 2015) and is able to determine the exact number of pulsar rotations also in the longer observational gaps. This was particularly useful for wider orbit, fainter pulsars.

In Table 2 we list, for each new pulsar, the span of radio data used to obtain the timing solution, the number and total duration of the timing observations at each telescope included in the analysis, and the total number of ToAs obtained.

The time span of each data set (also including, whenever possible, the survey discovery observation) was determined by the number of observations needed to either get a coherent timing solution, including the measurement of the pulsar spin period first derivative or to obtain a good enough ephemeris, which would enable the detection of gamma-ray pulsations, as described in the next section.

### 3. Gamma-ray observations and pulsation searches

For each newly discovered MSP, upon obtaining an initial phase-connected radio timing solution, we performed a search for gamma-ray pulsations in the *Fermi*-LAT data, and used these pulsations, when detected, to improve the timing ephemerides (see Sect. 4). For these analyses, we used *Fermi*-LAT data according to the ‘Pass 8’ P8R3_SOURCE_V3 instrument response functions (Atwood et al. 2013; Bruel et al. 2018), with the same energy-dependent photon selections that were used to produce the 4FGL-DR3 source catalogue (Abdollahi et al. 2022). Gamma-ray pulsation searching and timing relies on photon probability weights, derived from a spectral and spatial model of the gamma-ray sky, to disentangle fore- and background events from photons emitted by the target pulsar (Bickel et al. 2008; Kerr 2011; Bruel 2019). These weights were computed using the *gtsrcprob* routine, using the positions and spectra of sources from the 4FGL-DR3 (Abdollahi et al. 2022), as well as the *gll_iem_v07.fits* and *iso_P8R3_SOURCE_V3_v1.txt* Galactic and isotropic diffuse emission models.

The distribution of photon weights, *w*, for a source can be used to estimate the detectability of gamma-ray pulsations from a newly discovered pulsar, prior to searching. The significance of gamma-ray pulsations is typically quantified by the $H_M$-test (Kerr 2011), which incoherently sums the weighted Fourier power up to $M = 20$ harmonics. The expected value for $H_M$ is proportional to $W^2 = \sum_i w_i^2$ (Nieder et al. 2020a), with a prefactor that depends on the pulse profile shape, but which is typically of order unity. The $W^2$ quantity can also guide a minimum cutoff ($w_{min}$) on the photon weights, used to remove the large number of low-weight photons that do not contribute meaningfully to the pulsation detection to speed up computations during searching or timing. We chose $w_{min}$ values for which 99% of the full $W^2$ was retained, typically corresponding to $w_{min} \approx 0.01-0.03$.

<footnote n="4">⁴ https://bitbucket.org/psrsoft/tempo2/src/master/</footnote>

<footnote n="5">⁵ https://github.com/pfreire163/Dracula</footnote>
For four pulsars, PSRs J1623-6936, J1709-0333, J1757-6032, and J1858-5422, the radio ephemerides obtained in Sect. 2 were sufficiently accurate that a single fold of the Fermi-LAT data revealed pulsations reaching back to the beginning of the Fermi mission in 2008, albeit with some phase drifts requiring small adjustments to the spin frequency, $f$, and its first derivative, $\dot{f}$, which were consistent with the radio-derived uncertainties on these parameters. For a fifth MSP, PSR J1803-6707, whose gamma-ray pulsations were presented in Clark et al. (2023), significant pulsations were detected within the epochs covered by the radio timing solution but unpredictable orbital period variations prevented an extrapolation to earlier data. We address the gamma-ray timing for this pulsar in Sect. 4.2.

For one of these pulsars, the isolated MSP J1709-0333, pulsations were initially detected only faintly (with $H = 44$, after tweaking $\dot{f}$ slightly), and the resulting pulse profile contains a significant unpulsed component above the estimated background level, indicating that photon weights may be overestimated. In addition, the timing position for this pulsar lies 5 arcmin from the corresponding 4FGL-DR3 source, 4FGL J1709.3-0328, well outside the 4 arcmin 95% localisation ellipse. Taken together, these features suggest that the pulsar and background gamma-ray spectra are mis-modelled in this region, with some background flux perhaps being misattributed to 4FGL J1709.3-0328 in the model. This is perhaps unsurprising as, at Galactic coordinates $(l, b) = (17.4^{\circ}, +20.7^{\circ})$, this pulsar lies in a particularly complex region of the diffuse gamma-ray sky, along the edge of the northern Fermi Bubble (Su et al. 2010; Ackermann et al. 2014) and in a region with a significant dust excess (Acero et al. 2016). To address this, we recomputed the photon weights with gtsrcprob using the pulsar timing position instead of the 4FGL position, and then applied an energy-dependent re-weighting to maximise the $H$-test, similar to the “model weights” method developed by Bruel (2019). For this, we used the photon re-weighting equation derived by Kerr (2019), $w' = sw/(sw+1-w)$, multiplying the source flux with a power-law spectral function $s(E) = A(E/1000\,\text{MeV})^{\gamma}$, and varied the overall scale factor $A$ and spectral index $\gamma$ to maximise the re-weighted $H$-test. We found that $A = 0.3$ and $\gamma = 0.7$ results in the maximum $H = 76$, indicating that the low energy spectrum of the pulsar/background is over/under-estimated, respectively, or some combination thereof, in the 4FGL-DR3 model.

For the other MSPs the validity intervals of the radio solutions were too short for significant gamma-ray pulsations to be detectable within their validity period, and timing uncertainties were large enough that phase connection was not maintained when extrapolating to earlier gamma-ray data. In these cases it was necessary to perform a ‘brute force’ search over the uncertain parameters ($f, \dot{f}, P_b, \alpha$ and $\delta$). To save computing resources and given that most pulsars have the most Fourier power in the lowest harmonics (Pletsch & Clark 2014), we chose to truncate $H_M$ at $M = 3$. Significant detections have been claimed when $H_M > 25$ from a single trial (Smith et al. 2019), but much larger values are required to overcome the ‘trial factors’ that arise when a range of parameters must be searched over (e.g. Nieder et al. 2019).

As already reported in Clark et al. (2023), PSR J1526-2744 was detected in this way, but gamma-ray pulsations were not detected from the remaining three MSPs, PSRs J1036-4353, J1823-3543, and J1906-1754, despite these efforts. For the first two of these pulsars, this is most likely explained by the low gamma-ray flux. For PSR J1036-4353, in particular, the radio timing ephemeris covers around 500 days, and contains four

orbital frequency derivatives describing short-timescale variations in the orbital period. However, $W^2 = 134$ for this pulsar, meaning that we would not expect to detect pulsations on timescales shorter than around 3 years, but the orbital phase variations mean that the timing model becomes uncertain on far shorter timescales, preventing a detection of gamma-ray pulsations from this pulsar. PSR J1823-3543 is even fainter, with $W^2 = 67$, and we would require a timing model covering nearly all of the Fermi-LAT data to detect it. We suspect that its signal is buried in the noise of the large number of trials that were required in our search. The non-detection of PSR J1906-1754 is harder to explain (see Sect. 5.3), as with $W^2 = 111$ we would expect to detect pulsations even in a 5D search around the radio ephemeris.

# 4. Joint radio and gamma-ray timing

Radio and gamma-ray data are very different in various aspects, which directly affects the timing accuracy. The gamma-ray data are sparse, sometimes with millions of rotations between two individual gamma-ray photons from the pulsar, but cover the full 15 years between the launch of the Fermi satellite and today. Such data are well suited to measure pulsar parameters for which the precision depends directly on the data time span, namely, precision improves more rapidly than with the observing span $t^{-1/2}$, such as the spin-frequency derivative, $\dot{f}$, the orbital period, $P_b$, or proper motion. The radio data give precise pulse arrival times, but have gaps and usually (unless archival data are found at the pulsar position) only span the time between discovery and the most recent timing observation. Due to their higher S/N, the radio data typically measure parameters, such as the sky position or the binary projected semi-major axis, with a higher precision.

Gamma-ray pulsations have been detected from six of the pulsars presented here; thus, using both data sets should increase sensitivity in the timing analysis for these. For several of these pulsars, (e.g. PSR J1623-6936), the typical iterative timing approach of alternating radio and gamma-ray timing analyses did not converge on a consistent sky position, because of the relatively short time span of the radio data. To fix this, we used a timing approach in which we fit the pulsar parameters utilising both data jointly.

## 4.1. Joint timing technique

Our joint radio and gamma-ray timing analysis code<fnref n="6" /> uses the pulsar timing software PINT<fnref n="7" /> (Luo et al. 2021). With PINT the radio and gamma-ray data are folded individually for a given set of pulsar parameters. For these, a joint log-likelihood function, $\log \mathcal{L}_j$ is constructed from the two respective standard timing statistics – the radio $\chi_r^2$ statistic and the gamma-ray template-based log-likelihood $\log \mathcal{L}_g$ – as

$$ \log \mathcal{L}_j = \log \mathcal{L}_p - 0.5\chi_r^2 + \log \mathcal{L}_g, \quad (1) $$

where $\log \mathcal{L}_p$ represents the log-priors. The radio $\chi_r^2$ statistic is computed from the ToA residuals and uncertainties. $\log \mathcal{L}_g$ is the commonly used ‘unbinned’ log-likelihood statistic (Kerr et al. 2015). This parameter gives the likelihood that the distribution of pulsar phases at the emission times of the gamma-ray photons is described by a given template pulse profile, with the template consisting of wrapped Gaussian peaks (Ray et al. 2011).

<footnote n="6">⁶ https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/jrag-timing.py</footnote>

<footnote n="7">⁷ https://github.com/nanograv/PINT</footnote>
We performed the full joint radio and gamma-ray timing analysis on four binary pulsars and one isolated pulsar. Using both radio and gamma-ray data, we fit for the astrometry, spin, and orbital parameters, as well as marginalising over the gamma-ray pulse-profile template (Nieder et al. 2019), using the Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler emcee (Foreman-Mackey et al. 2013) to explore the parameter space. The timing solutions are presented in Tables 3, 4, and 5, together with the radio-only solutions of the three pulsars with no gamma-ray pulsed counterpart.

Solving the problem of the non-convergence seen for some of the pulsars was not the only achievement of the joint radio and gamma-ray timing approach. The higher sensitivity enabled the measurement of the pulsars' proper motion, which was not significant in either data sets individually (illustrated in Fig. 1), making it possible to account for the kinematic contributions to the observed spin-down rate to more accurately estimate the energetics of these pulsars (see Sect. 5). The benefit of using radio and gamma-ray data jointly is demonstrated in the panels which show sky position parameters (RAJ, DECJ) versus the proper motion (PMRA, PMDEC). While the former are more accurately measured due to the higher S/N of the radio data, the proper motion measurement benefits from the longer time span of the gamma-ray data set.

## 4.2. Gamma-ray timing of the redback PSR J1803–6707

Timing analyses of RBs are significantly complicated by long-term variations in their orbital periods, believed to be due to changes in the gravitational quadrupole moment of their companion stars (Applegate & Shaham 1994; Lazaridis et al. 2011). If they are not correctly accounted for, these lead to gamma-ray pulses apparently disappearing from data outside the range of validity of the radio ephemerides, which cannot be deterministically extrapolated over longer time spans. This effect led to us being initially unable to obtain a long-term gamma-ray timing solution for PSR J1803–6707 beyond the interval covered by the radio timing in Clark et al. (2023). Since then, we have refined our methods for timing gamma-ray RBs, allowing us to resolve this. These methods, also used in Thongmeearkom et al. (2024), build upon the procedure developed in Clark et al. (2021), in which we treat the orbital phase variations as a stationary Gaussian process.

Starting from the radio ephemeris presented in Clark et al. (2023), we first searched for pulsations throughout the gamma-ray data, using a 500-day sliding window to strike a balance between including enough data to detect significant pulsations, while also being robust against short-timescale variations in the orbital phase. This search spanned small ranges in epoch of ascending node, $T_{asc}$, orbital period $P_b$, and the first derivative of the orbital frequency, $FB1$. This resulted in a connected ‘track’ of detections dating back to the start of the Fermi-LAT data, to which we could fit an initial Gaussian process model describing the orbital phase changes. The results of this sliding window approach are illustrated in Fig. 2. We used this same sliding window approach to try to find gamma-ray pulsations from PSR J1036–4353 but did not detect significant pulsations, which is unsurprising given the low weighted-photon count from this pulsar.

Using this as starting point, we then performed a full gamma-ray timing analysis in which we jointly fit the parameters of the template pulse profile and the timing model (astrometry, spin and orbital parameters), alongside the orbital phase variation

<table><caption>**Table 3.** Timing parameters for the RBs J1036–4353 and J1803–6707.</caption><thead><tr><th colspan="3">Timing parameters</th></tr><tr><th></th><th>J1036–4353</th><th>J1803–6707</th></tr></thead><tbody><tr><td>α (h:m:s)</td><td>10:36:30.215127(14)</td><td>18:03:04.235339(10)</td></tr><tr><td>δ (deg:'.''')</td><td>-43:53:08.7252(3)</td><td>-67:07:36.15763(15)</td></tr><tr><td>μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-11.6(3)</td><td>-8.43(13)</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>2.9(3)</td><td>-6.3(1)</td></tr><tr><td>ω (mas)</td><td>0.4(4)</td><td>0.2(3)</td></tr><tr><td>PosEpoch (MJD)</td><td>57388</td><td>57388</td></tr><tr><td>f (Hz)</td><td>595.1998208151(4)</td><td>468.46771214902(10)</td></tr><tr><td>f&#775; (Hz s<sup>-1</sup>)</td><td>-2.26(3)×10<sup>-15</sup></td><td>-4.0434(6)×10<sup>-15</sup></td></tr><tr><td>PEpoch (MJD)</td><td>59646.78438</td><td>59364.893677</td></tr><tr><td>DM (pc cm<sup>-3</sup>)</td><td>61.119(3)</td><td>38.382(3)</td></tr><tr><td>Binary</td><td>BTX</td><td>ELL1</td></tr><tr><td>P<sub>b</sub> (d)</td><td>0.259621257(6)</td><td>0.38047333(11)</td></tr><tr><td>x (lt-s)</td><td>0.664789(3)</td><td>1.061914(2)</td></tr><tr><td>T<sub>asc</sub> (MJD)</td><td>59536.3052068(8)</td><td>59020.9969(5)</td></tr><tr><td>FB1 (Hz s<sup>-1</sup>)</td><td>-2.3(3)×10<sup>-18</sup></td><td>—</td></tr><tr><td>FB2 (Hz s<sup>-2</sup>)</td><td>7.9(6)×10<sup>-25</sup></td><td>—</td></tr><tr><td>FB3 (Hz s<sup>-3</sup>)</td><td>-1.12(7)×10<sup>-31</sup></td><td>—</td></tr><tr><td>FB4 (Hz s<sup>-4</sup>)</td><td>6.0(4)×10<sup>-39</sup></td><td>—</td></tr><tr><td>RMS (μs)</td><td>14.803</td><td>11.547</td></tr><tr><th colspan="3">Hyperparameters of orbital-phase-covariance function</th></tr><tr><td>h (s)</td><td>—</td><td>17<sup>+7</sup><sub>-10</sub></td></tr><tr><td>ℓ (d)</td><td>—</td><td>>560</td></tr><tr><td>ν</td><td>—</td><td>1.84<sup>+0.44</sup><sub>-0.41</sub></td></tr><tr><th colspan="3">Derived parameters</th></tr><tr><td>P (ms)</td><td>1.6801080326781(12)</td><td>2.1346188308481(12)</td></tr><tr><td>P&#775;<sub>int</sub> (s s<sup>-1</sup>)</td><td>5.50×10<sup>-21</sup></td><td>1.84 × 10<sup>-20</sup></td></tr><tr><td>E&#775; (erg s<sup>-1</sup>)</td><td>4.58×10<sup>34</sup></td><td>7.47 × 10<sup>34</sup></td></tr><tr><td>B<sub>S</sub> (G)</td><td>9.73×10<sup>7</sup></td><td>2.01×10<sup>8</sup></td></tr><tr><td>M<sub>2,min</sub> (M<sub>☉</sub>)</td><td>0.23</td><td>0.29</td></tr><tr><td>d<sub>NE2001</sub> (kpc)</td><td>2.1</td><td>1.2</td></tr><tr><td>d<sub>YMW16</sub> (kpc)</td><td>0.4</td><td>1.4</td></tr><tr><td>μ (mas yr<sup>-1</sup>)</td><td>12.0(3)</td><td>10.6(3)</td></tr><tr><td>ν<sub>⊥</sub> (km s<sup>-1</sup>)</td><td>119</td><td>60</td></tr><tr><td>S<sub>1300</sub> (mJy)</td><td>0.14</td><td>0.17</td></tr></tbody></table>

**Notes.** The top part of the table reports the measured timing parameters. From top to bottom we list: right ascension (α), declination (δ), proper motion in α and δ, parallax (ω), epoch of the reported position, spin frequency (f), spin frequency derivative, epoch of the reported f, dispersion measure (DM taken from the best detection obtained at 800 MHz at MeerKAT), binary model used (Hobbs et al. 2006), orbital period ($P_b$), projected semi-major axis (x), time of the ascending node ($T_{asc}$), four orbital frequency derivatives and the Tempo2 rms of the timing residuals. Position, proper motion and parallax are taken from Gaia DR3 (Gaia Collaboration 2023). For each parameter, the (Gaia or timing) one-sigma error on the last quoted digit(s) is reported in parenthesis. The hyperparameters of the orbital-phase-covariance function are the amplitude, h, the length scale, ℓ, and the Matérn degree, ν. Their values are reported with 95% confidence intervals. The bottom part of the table lists the derived parameters. Top to bottom: spin period, intrinsic spin period derivative $P_{int}$ corrected for the Shklovskii and Galactic potential effects, spin-down power and surface magnetic field (calculated using the corrected $P_{int}$), minimum companion mass (assuming a neutron star mass of 1.4 M<sub>☉</sub>), DM-derived distances (using the NE2001 and YMW16 models), total proper motion and transverse velocity, using the NE2001 DM-derived distance and estimate of the 1300-MHz flux density. The full ephemeris for J1803–6707 is available as supplementary material on Zenodo (see Data availability section).

Gaussian process and three ‘hyperparameters’ describing its covariance function. We adopted the commonly used Matérn
<table><caption>**Table 4.** Measured and derived timing parameters for the non-eclipsing binary pulsars.</caption><thead><tr><th></th><th colspan="6">Timing parameters</th></tr><tr><th></th><th>J1526-2744</th><th>J1623-6936</th><th>J1757-6032</th><th>J1823-3543</th><th>J1858-5422</th><th>J1906-1754</th></tr></thead><tbody><tr><td>α (h:m:s)</td><td>15:26:45.096(2)</td><td>16:23:51.3868(6)</td><td>17:57:45.4474(3)</td><td>18:23:42.9989(4)</td><td>18:58:07.7664(3)</td><td>19:06:14.7894(4)</td></tr><tr><td>δ (deg:′:″)</td><td>-27:44:06.02(9)</td><td>-69:36:49.254(3)</td><td>-60:32:12.300(4)</td><td>-35:43:40.879(13)</td><td>-54:22:15.527(74)</td><td>-17:54:34.31(4)</td></tr><tr><td>μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-12(4)</td><td>-5.2(11)</td><td>-2.3(14)</td><td>–</td><td>-8.7(14)</td><td>–</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-4(14)</td><td>-7.1(13)</td><td>-3.0(23)</td><td>–</td><td>5(3)</td><td>–</td></tr><tr><td>PosEpoch (MJD)</td><td>59355.468037</td><td>59332</td><td>59400</td><td>59091</td><td>59563.356205</td><td>59700</td></tr><tr><td>f (Hz)</td><td>401.7446020972(2)</td><td>414.95962743561(13)</td><td>343.34786701089(10)</td><td>421.3525986326(14)</td><td>424.53897607040(8)</td><td>347.66288115516(17)</td></tr><tr><td>ḟ (Hz s<sup>-1</sup>)</td><td>-5.724(10)×10<sup>-16</sup></td><td>-1.5598(7)×10<sup>-15</sup></td><td>-3.529(6)×10<sup>-16</sup></td><td>-3.736(5)×10<sup>-15</sup></td><td>-7.420(5)×10<sup>-16</sup></td><td>-4.73(11)×10<sup>-16</sup></td></tr><tr><td>PEpoch (MJD)</td><td>59355.468037</td><td>59332</td><td>59400</td><td>59091</td><td>59563.356205</td><td>59700</td></tr><tr><td>DM (pc cm<sup>-3</sup>)</td><td>30.953(3)</td><td>46.414(3)</td><td>62.918(3)</td><td>81.686(3)</td><td>30.819(3)</td><td>98.090(3)</td></tr><tr><td>Binary</td><td>ELL1</td><td>ELL1</td><td>ELL1</td><td>ELL1</td><td>ELL1</td><td>BT</td></tr><tr><td>P<sub>b</sub> (d)</td><td>0.2028108285(3)</td><td>11.01366617(2)</td><td>6.280139735(6)</td><td>144.568639(3)</td><td>2.581951179(5)</td><td>6.4910850(3)</td></tr><tr><td>x (lt-s)</td><td>0.224080(9)</td><td>6.728198(4)</td><td>9.612624(2)</td><td>51.733306(3)</td><td>1.6889079(11)</td><td>1.346404(7)</td></tr><tr><td>T<sub>asc</sub> (MJD)</td><td>59303.2059794(10)</td><td>59192.9162236(12)</td><td>59183.4010969(6)</td><td>59091.325976(8)</td><td>59564.9829160(6)</td><td>59700.935268(5)</td></tr><tr><td>EPS1</td><td>–</td><td>2.28(14)×10<sup>-5</sup></td><td>1.1(6)×10<sup>-6</sup></td><td>3.991(7)×10<sup>-5</sup></td><td>–</td><td>–</td></tr><tr><td>EPS2</td><td>–</td><td>4.54(9)×10<sup>-5</sup></td><td>-1.3(5)×10<sup>-6</sup></td><td>2.6631(10)×10<sup>-4</sup></td><td>–</td><td>–</td></tr><tr><td>RMS (μs)</td><td>47.180</td><td>34.153</td><td>24.597</td><td>16.975</td><td>21.312</td><td>32.570</td></tr><tr><th></th><th colspan="6">Derived parameters</th></tr><tr><td>P (ms)</td><td>2.4891435871939(19)</td><td>2.4098729945848(8)</td><td>2.9124980699779(7)</td><td>2.373309202899(8)</td><td>2.3554963298210(9)</td><td>2.8763496312214(15)</td></tr><tr><td>Ṗ<sub>int</sub> (s s<sup>-1</sup>)</td><td>3.88×10<sup>-21</sup></td><td>8.57×10<sup>-21</sup></td><td>2.24×10<sup>-21</sup></td><td>2.08×10<sup>-20</sup></td><td>3.70×10<sup>-21</sup></td><td>3.61×10<sup>-21</sup></td></tr><tr><td>Ė (erg s<sup>-1</sup>)</td><td>6.48×10<sup>32</sup></td><td>2.42×10<sup>34</sup></td><td>3.58×10<sup>33</sup></td><td>6.15×10<sup>34</sup></td><td>1.12×10<sup>34</sup></td><td>5.98×10<sup>33</sup></td></tr><tr><td>B<sub>S</sub> (G)</td><td>7.83×10<sup>7</sup></td><td>1.45×10<sup>8</sup></td><td>8.17×10<sup>7</sup></td><td>2.25×10<sup>8</sup></td><td>9.44×10<sup>7</sup></td><td>1.03×10<sup>8</sup></td></tr><tr><td>M<sub>2,min</sub> (M<sub>☉</sub>)</td><td>0.08</td><td>0.19</td><td>0.42</td><td>0.26</td><td>0.12</td><td>0.05</td></tr><tr><td>e</td><td>&lt;2.2×10<sup>-4</sup></td><td>5.00(13)×10<sup>-5</sup></td><td>1.5(8)×10<sup>-6</sup></td><td>2.6928(10)×10<sup>-4</sup></td><td>&lt;9.7×10<sup>-6</sup></td><td>–</td></tr><tr><td>ω (degrees)</td><td>–</td><td>26.0(20)</td><td>44(25)</td><td>171.477(16)</td><td>–</td><td>–</td></tr><tr><td>d<sub>NE2001</sub> (kpc)</td><td>1.2</td><td>1.3</td><td>1.8</td><td>2.1</td><td>1.0</td><td>2.9</td></tr><tr><td>d<sub>YMW16</sub> (kpc)</td><td>1.3</td><td>1.3</td><td>3.5</td><td>3.7</td><td>1.2</td><td>6.8</td></tr><tr><td>μ (mas yr<sup>-1</sup>)</td><td>9(7)</td><td>8.6(18)</td><td>8(4)</td><td>–</td><td>9(3)</td><td>–</td></tr><tr><td>v<sub>⊥</sub> (km s<sup>-1</sup>)</td><td>51</td><td>53</td><td>65</td><td>–</td><td>45</td><td>–</td></tr><tr><td>S<sub>1300</sub> (mJy)</td><td>0.07</td><td>0.10</td><td>0.06</td><td>0.17</td><td>0.06</td><td>0.12</td></tr></tbody></table>

**Notes.** For each parameter, the Tempo2 one-sigma error on the last quoted digit(s) is reported in parenthesis. The parameters listed (except for the non measured parallax and orbital frequency derivatives) are the same as in Table 3. For two pulsars the ELL1 binary model (Lange et al. 2001) was used to fit orbits with low eccentricities (e); the reported values for e and the longitude of periastron ω have been derived from the EPS1 and EPS2 Laplace parameters given by the model, the limits are 95% upper limits. Pulsars whose names are in bold have been jointly timed in radio and gamma-rays. The reported spin period derivative (and derived parameters) have been corrected for the Shklovskii effect whenever a significant total proper motion μ was measured.

<table><caption>**Table 5.** Timing parameters for the isolated MSP J1709-0333.</caption><thead><tr><th></th><th colspan="2">Timing parameters</th></tr><tr><th></th><th>J1709-0333</th><th></th></tr></thead><tbody><tr><td>α (h:m:s)</td><td>17:09:32.65211(12)</td><td></td></tr><tr><td>δ (deg:′,″)</td><td>-03:33:19.250(5)</td><td></td></tr><tr><td>f (Hz)</td><td>283.76747691122(6)</td><td></td></tr><tr><td>ḟ (Hz s<sup>-1</sup>)</td><td>-1.510(5)×10<sup>-16</sup></td><td></td></tr><tr><td>PEpoch (MJD)</td><td>59694</td><td></td></tr><tr><td>DM (pc cm<sup>-3</sup>)</td><td>35.791(6)</td><td></td></tr><tr><td>RMS (μs)</td><td>25.851</td><td></td></tr><tr><th></th><th colspan="2">Derived parameters</th></tr><tr><td>P (ms)</td><td>3.5240120216913(8)</td><td></td></tr><tr><td>Ṗ<sub>int</sub> (s s<sup>-1</sup>)</td><td>1.89×10<sup>-21</sup></td><td></td></tr><tr><td>Ė (erg s<sup>-1</sup>)</td><td>1.70×10<sup>33</sup></td><td></td></tr><tr><td>B<sub>S</sub> (G)</td><td>8.25×10<sup>7</sup></td><td></td></tr><tr><td>d<sub>NE2001</sub> (kpc)</td><td>1.4</td><td></td></tr><tr><td>d<sub>YMW16</sub> (kpc)</td><td>0.2</td><td></td></tr><tr><td>S<sub>1300</sub> (mJy)</td><td>0.05</td><td></td></tr></tbody></table>

**Notes.** For each parameter, the Tempo2 one-sigma error on the last quoted digit(s) is reported in parenthesis. The parameters listed (except for the missing proper motion, compatible with zero, at 1-sigma, and the orbital ephemeris) are the same as in the previous table.

covariance function, which is parameterised by an amplitude, *h*, length scale, *ℓ*, and a smoothness parameter, *v*. This covariance function corresponds to a noise power spectral density that follows a smoothly broken power law, that is flat below a corner frequency $f_c = \sqrt{\nu/\sqrt{2\pi}\ell}$, breaking to a power-law at higher frequencies with a spectral index of $\alpha = -(2\nu + 1)$. As described in Thongmeearkom et al. (2024) we used a Gibbs sampling approach to deal with the multi-component likelihood function for each photon. This approach also allows radio ToAs to be incorporated in the sampling. Details on these methods and a public release of the code will be presented in future works (in prep.).

Since PSR J1803-6707 has a bright optical counterpart (see Clark et al. 2023), we used Gaia DR3 measurements as priors on the astrometric parameters (position, proper motion and parallax). In this case, the gamma-ray timing is not sensitive enough to improve upon these. The results of this timing procedure are shown in Fig. 2, and the full gamma-ray timing solution is provided as supplementary material on Zenodo (see Data availability section) in two versions ('ORBIFUNC' and 'ORBWAVES') that can be used with Tempo2 and PINT, respectively<fnref n="8" />.

<footnote n="8">⁸ The ORBIFUNC parameters (for Tempo2) are a time-domain interpolating function describing the orbital phase variations over time, while the ORBWAVES parameters (for PINT) are the Fourier-domain representation of the same function.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Corner plot for the astrometry parameters for PSR J1858−5422, showing resulting posterior distributions for radio-only (red), gamma-ray- only (blue) and joint-radio-and-gamma-ray (black) timing analyses. The 1D posterior distributions on the diagonal are given in arbitrary units. The high S/N data allow us to pin down the positional parameters more precisely, while the proper motion parameters are measured more accurately by the long-term gamma-ray data. The sensitivity improvement of the joint analysis allows the precise measurement of all four parameters. The complete plot with posterior distributions for all pulsar parameters is shown in Fig. A.5.</figcaption></figure>

## 5. Timing results

The results of the timing campaign are summarised in Tables 3, 4, and 5, where the best set of measured and derived timing parameters are reported for the two RBs, the other binary pulsars and the isolated MSP J1709−0333, respectively. Pulsars whose timing parameters have been obtained from the joint analysis of the radio ToAs and gamma-ray photons are in boldface font. Barycentric Dynamical Time (TDB) scale and the DE421 JPL Solar system ephemeris were used for the timing analysis of all pulsars.

For PSR J1036−4353, the lack of a gamma-ray pulsation detection and the relatively short interval spanned by the radio data make the Gaussian-process treatment (described in Sect. 4.2) unnecessary. Instead, we include four orbital frequency derivatives in the timing model to obtain flat timing residuals over the radio data span.

When reporting parameters derived from the spin-down rate of the pulsar (e.g. the surface magnetic field B and the spin-down energy $\dot{E}$), one needs to take into account that the observed first derivative of the pulsar spin period $\dot{P}_{\text{obs}}$ is contaminated by the Galactic acceleration of the pulsar relative to the Earth and by its

proper motion. The intrinsic value $\dot{P}_{\text{int}}$ can be written as:

$$ \dot{P}_{\text{int}} = \dot{P}_{\text{obs}} - \dot{P}_{\text{Gal}} - \dot{P}_{\text{Shk}}, \quad (2) $$

where $\dot{P}_{\text{Gal}}$ is the contribution due to the Galactic gravitational acceleration and $\dot{P}_{\text{Shk}}$ is the Shklovskii effect (Shklovskii 1970), a purely kinematic contribution to the observed period derivative, which we could compute for those pulsars for which a significant measurement of the proper motion was available either from the Gaia catalogue or as a result of our joint timing. The Galactic contribution $\dot{P}_{\text{Gal}} = P \frac{a_{\text{Gal}}}{c}$ (where P is the spin period of the pulsar and $a_{\text{Gal}}$ is the acceleration in the Galactic potential) is calculated using the model from McMillan (2017). The Shklovskii effect is expressed as:

$$ \dot{P}_{\text{Shk}} = \frac{P \mu^2 d}{c} = \frac{P v_{\perp}^2}{c d}, \quad (3) $$

where $\mu$ is the total proper motion, d the distance to the pulsar, and $v_{\perp}$ the transverse velocity.

The distance estimates reported in Tables 3, 4, and 5 are derived from the dispersion measure (DM) of the pulsars
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Gamma ray pulsations and orbital phase variations over the course of the Fermi-LAT data for PSR J1803−6707. The panels on the left show the weighted gamma-ray photon phases for the highest-likelihood timing solution (lower panel) and the integrated pulse profile (upper panel). The orange and the 100 faint black curves represent the associated highest-likelihood pulse-profile template and randomly drawn samples from the Monte Carlo analysis to visualise the uncertainty on the pulse-profile template. The bottom-right panel shows the orbital phase variations as a function of time. The grey-scale image shows the log-likelihood for offsets from the pulsar's $T_{\text{asc}}$ as measured in overlapping 800-day windows. The red curves represent the 95% confidence interval on those deviations, obtained from Monte Carlo timing analysis, with a 5 s-offset for clarity. The top-right panel shows the power spectral density of the orbital phase variations. The red shaded regions illustrate the 68% and 95% confidence intervals of the Matérn model. The black curve and grey shaded region show the estimated power spectrum and its 95% confidence interval from the joint radio and gamma-ray fitting, while the magenta curve, mostly hidden behind the black one, shows the power spectrum estimated from gamma-ray timing alone. The magenta curve is added to emphasise the increased sensitivity at higher frequencies caused by utilising the radio data.</figcaption></figure>

and from two different models (NE2001 and YMW16) of the distribution of free electrons in the Galactic interstellar medium (Cordes & Lazio 2002; Yao et al. 2017). For the two RBs we can compare the DM-derived distances (considering a 30% error on their estimates) with the 2 sigma lower limit derived from the Gaia parallax measurement for the companion stars. For J1036−4353, the Gaia distance (>0.9 kpc) is incompatible with the YMW16 value (0.4 kpc), while the lower limit for J1803−6707 (1.25 kpc) is in line with the DM-derived values from both models (1.2 kpc and 1.4 kpc for NE2001 and YMW16, respectively). For the calculation of the transverse velocity $v_{\perp}$ and of $\dot{P}_{\text{int}}$ we hence chose to use the distance derived from the NE2001 model in all cases.

Since timing data were not flux calibrated, an estimate of the flux density has been obtained using the radiometer equation modified for pulsars (see e.g. Dewey et al. 1985). Figure 3 shows our nine MSPs as red stars in the period-period derivative diagram. As already noted in Ray et al. (2012) for the entire MSP population discovered in Fermi sources (drawn as black points in the diagram), the pulsars in our small sample appear to have shorter spin periods with respect to the general population of Galactic MSPs (plotted as grey circles and defined, here, as pulsars with spin period smaller than 10 ms). Our pulsars have a median *P* of 2.4 ms, 1.2 ms shorter than that of the Galactic population (and 0.4 ms smaller than that of the MSPs previously discovered targeting gamma-ray point sources).
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Period-period derivative diagram for Galactic millisecond pulsars. Grey circles show pulsars from PSRCat 1.70, while smaller black points denote MSPs found targeting Fermi unassociated point sources. Red stars are the MSPs from this paper. The *Ṗ* used is, whenever possible, the one corrected for the Shklovskii effect. Equal spin-down energy lines are drawn as dotted and equal surface magnetic field lines are dashed.

In the following, we highlight some interesting features of
our timed targets.

5.1. PSR J1757−6032 and the detection of its Shapiro delay

This pulsar is in a binary system with an orbital period of 6.28 d with a very small orbital eccentricity and a relatively large projected semi-major axis of 9.6 light seconds. If we assume a pulsar mass of 1.4 M☉, then the minimum companion mass is 0.42 M☉. This is significantly larger than the prediction of Tauris & Savonije (1999) for the mass of a helium white dwarf (He WD) in a 6.3-d orbit around a MSP, which is about 0.24 M☉; this implies that the companion is likely a CO white dwarf.

The basic characteristics of the system resemble
PSR J1614−2230 (Demorest et al. 2010) and a few other
systems like PSR J1933−6211 (see detailed discussions in
Geyer et al. 2023 and Shamohammadi et al. 2023). Compared
to the other known pulsars with similarly high companion
masses, these systems have 1 to 2 orders of magnitude faster
spin periods, lower eccentricities (three orders of magnitude,
on average) and smaller magnetic fields (up to three orders of
magnitude).

These characteristics can be explained if the system origi-
nated in Case A Roche lobe overflow (RLO), namely, the pulsar
started accreting mass from its companion when the latter was
still in the main sequence (Tauris et al. 2011). In this model the
unusually long accretion episode results in the fast spin, small B-
fields and unusually small orbital eccentricities, even compared
to those of MSP-He WD systems.

One of the reasons why these systems are interesting is
that the measurement of the Shapiro delay (Shapiro 1964)

in PSR J1614−2230 showed that the pulsar is quite massive
(Demorest et al. 2010), a measurement that is important for the
study of the EoS of nuclear matter at densities above that of the
atomic nucleus (see e.g. Özel & Freire 2016). The latest mass
measurements for this system (e.g. Shamohammadi et al. 2023;
Agazie et al. 2023a) confirm a pulsar mass around 1.94 M⊙.

When studying the evolution of the system, Tauris et al. (2011) suggested that, despite its length, the accretion episode is not responsible for the high mass of PSR J1614−2230, the latter is more likely a birth property. This was to some extent confirmed recently by the measurement of the mass of PSR J1933−6211: despite its strong evolutionary similarities with PSR J1614−2230 it has a mass around 1.4 M☉ (Geyer et al. 2023). However, the precision of the mass of PSR J1933−6211 still needs to be improved, and other similar systems, like PSR J1101−6424, do not yet provide precise measurements (Shamohammadi et al. 2023). In order to start answering the question of whether neutron star (NS) masses are determined from their birth or by subsequent evolution, we need to measure more pulsar masses from this type (Case A RLO) of system. These are especially promising for detecting the role of accretion in the mass of the pulsar, not only because in this case the accretion episode is exceptionally long, but also because the relatively large companion masses and short spin periods (with resulting potentially high timing precision) make the measurement of the Shapiro delay easier. Furthermore, if accretion adds a significant amount of mass to the pulsar, we might find additional massive pulsars in this class of binaries, which might introduce additional constraints on the EoS.

Interestingly, it might be possible to make a mass measure-
ment for PSR J1757-6032. We have a significant detection of
the orthometric amplitude of the Shapiro delay, h₃ = 3.4 ± 0.8 µs
(estimated assuming an inclination of 87.5 degrees, see below),
however, fitting for both parameters of the Shapiro delay (h₃,
ş, Freire & Wex 2010) yields no significant detection at present,
which precludes precise constraints on the masses. The lack of
precision of the Shapiro delay is caused by the insufficient pre-
cision of the radio timing (this is significantly worse for the
gamma-ray timing) and the currently limited number of obser-
vations near superior conjunction.

To better quantify the detection, we made a quality-of-fit (χ²) map of the M₂-cos i plane (mass of the companion vs cosine of the orbital inclination), restricting the pulsar masses to the causality limit at the upper end (3.2 M☉) and a value of 1.17 M☉ at the lower end, the lowest measured NS mass (Martinez et al. 2015). We chose cos i because this quantity has a uniform probability for randomly aligned orbits, constituting a reasonable prior. From this χ² map, we then derived a Bayesian joint posterior probability density function (pdf) and marginalised pdfs using the procedure described by Splaver et al. (2002). The results are displayed in Fig. 4. The main conclusion to take from this figure is that, for the whole range of NS masses we considered, the orbital inclination is always high, with a median value of 87.5±1.4 deg (68.3% confidence level) and 87.5±10.0 deg (95.4% confidence level). The range of companion masses (from 0.39 to about 0.75 M☉) results from this fact and the limits we set on the range of NS masses.

This high inclination indicates good prospects for the mea-
surement of the masses of the components of this system in the
near future via the two Shapiro delay parameters, provided that
(a) the timing precision improves and that (b) a long MeerKAT
observation can be made near superior conjunction as part of a
dense orbital campaign.
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Companion mass and inclination constraints for the PSR J1757−6032 binary system. The main panel shows the $M_2-\cos i$ plane, (with $M_2$ the mass of the companion star and $i$ the orbital inclination). The grey zone is excluded by the requirement that the pulsar mass must be positive. The gray scale shows the joint posterior probability density function (pdf), with white representing low or no probability. This pdf is limited, at the lower and upper sides, by pulsar masses of 1.17 and 3.2 $M_\odot$ (see text for justification); these are indicated by the blue lines. The red lines show the constraints derived from the median and $\pm1-\sigma$ measurements of $h_3$, measured assuming $i = 87.5$ deg (see text). The side panels include marginalized posterior probability distribution functions for $\cos i$ (top) and (to the right) $M_2$. The solid vertical lines indicate the medians of the distributions, the dotted lines indicate the 1-, 2-, and 3-$\sigma$ equivalent percentiles around the medians. Despite the lack of precise mass measurements, we can see that nearly edge-on orbital inclinations (i.e. with $\cos i$ near zero) are strongly preferred for the whole range of pulsar masses.

## 5.2. Redback pulsars J1036−4353 and J1803−6707

Figure 5 shows, for the two RBs J1036−4353 and J1803−6707, the orbital light curve obtained by summing all of the Parkes data. A very extended eclipse is visible between orbital longitudes ~20 deg and 160 deg (orbital phases 0.05 to 0.45) for J1036−4353 and between longitudes ~10 deg and 170 deg for J1803−6707.

Thanks to the UWL receiver at Parkes, the availability of data over a 3.3 GHz-wide frequency band can in principle allow for a phase-dependent study of the eclipse duration; this would provide an interesting insight to supplement our understanding of the nature and geometry of the eclipsing material (see e.g. Polzin et al. 2020 and Kansabanik et al. 2021). For PSR J1036−4353, however, the S/N of the summed data is not sufficient to create bright enough phase-binned orbital light curves to reliably study these effects over more than ~1 GHz in bandwidth. Figure 6 shows the orbital light curve of the pulsed flux, obtained by selecting $\pm 0.1$ in rotational phase around the peak, for the bottom part of the band (from 704 to 1728 MHz) split into four sub-bands with increasing bandwidth (to account for the steepness of the pulsar spectrum). No clear evidence for frequency evolution is visible.

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

**Fig. 5.** Waterfall plots showing, for the two RB pulsars J1036−4353 and J1803−6707, orbital longitude (with respect to the ascending node) versus pulse phase. The brighter the colour-scale, the stronger the signal. Orbital longitude ranges where the plots show a solid colour have no associated data, as observations were planned to avoid the eclipses. Data are from the Parkes telescope.

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

**Fig. 6.** Amplitude of the pulsed radio emission as a function of the longitude from the ascending node for PSR J1036−4353 in four different frequency sub-bands. The amplitude is measured over a pulse phase range of 0.2 around the peak. The central frequency and bandwidth of each sub-band are reported in the bottom-right corner. Points between longitudes 22° and 130° have not been plotted either because no data at all (between 34° and 94°) or not enough data had been recorded, resulting in a much higher level of noise.

For J1803−6707, for which the S/N is high enough to produce orbital light curves in different sub-bands across the whole UWL band, the combination of a non-uniform orbital coverage and of strong interstellar scintillation, prevents us from measuring the frequency dependence of the eclipse duration. What may appear, in Fig. 7, as a frequency-dependent longitude of eclipse ingress, for instance (where the eclipse seems to begin later at intermediate frequencies), is actually due to the fact that the orbital longitudes below ~20 deg were observed only once (on MJD 59243). On that occasion, scintillation made the signal prominent almost exclusively at frequencies around 1.7 GHz (4th panel). The same kind of effect is clearly visible in the 6th panel of Fig. 7, displaying 256 MHz of band centred around 2.1 GHz: here we see the effect of a very bright scintile present, at around 2 GHz, in one (of only three) observations covering the orbital longitudes between 270 deg and 320 deg. The phase
<figure><img src="image_8.png" /><figcaption>Fig. 7. Waterfall plots representing the orbital light curves (longitude with respect to the ascending node vs pulse phase) for PSR J1803-6707 in different frequency sub-bands. The central frequency and bandwidth are marked on the top of each panel. Only a pulse phase range of 0.4 around the peak is plotted.</figcaption></figure>

variations visible in the lowest frequency sub-band plotted in the leftmost box of Fig. 7 are likely the effect of dispersion measure delays caused by matter surrounding the system also away from the eclipse region; scintillation and profile evolution effects, however, can also contribute to create apparent pulse shape variations across the orbit, given the non-uniform orbital coverage.

### 5.3. The other binary pulsars

As already noted by Clark et al. (2023) in the discovery paper, given its minimum companion mass and orbital period, PSR J1526-2744 could also be a “spider” pulsar: either a relatively heavy BW, or a high-inclination RB (a large, albeit statistically unlikely, inclination would be required to both explain the small minimum companion mass with respect to the bulk of the RB population and the lack of visible eclipses in the system). The lack of a bright optical counterpart, also reported in Clark et al. (2023), seems, however, to favour the hypothesis that the companion star is a light (or on a highly inclined-orbit) He WD.

Despite having a minimum companion mass also in the range typical of BWs, PSR J1906-1754 has a much larger orbital period of 6.49 d and shows no eclipses. If the companion is a He WD, then the Tauris & Savonije (1999) relation predicts a mass of around $0.24 M_\odot$. If we assume a pulsar mass of $1.35 M_\odot$, this results in an orbital inclination of ~13 deg, which is a priori unlikely but should occur occasionally among MSP-He WD systems. Given the alignment of the spin with the orbital angular momentum one should expect for MSP-WD binaries, this would imply we are observing the pulsar from close to its spin axis. Two predictions stem from this: the first is that the radio pulse should be very broad, which is indeed observed: in the plot in Sect. 6, we see that this object has by far the widest pulse profile among all discoveries, with observable emission between spin phases of 0.2 to 1.0. Another consequence would be fainter observable gamma-ray emission, which is expected to be preferentially emitted around the pulsar's spin equator (e.g. Kalapotharakos et al. 2023), as well as a weaker modulation of the gamma-ray profile; this could be the reason why no gamma-ray pulsations have been detected for this system.

This speculative picture would be strengthened if the observed gamma-ray flux suggested a particularly low efficiency for this system, but this is not the case here. Given the rather uncertain DM-model distances for this pulsar ($d_{YMW16} = 6.8$ kpc vs. $d_{NE2001} = 2.9$ kpc), the energy flux above 100 MeV 4FGL J1906.4-1757 ($G_{100} = 3.4 \times 10^{-12}$ erg s$^{-1}$, Ballet et al. 2023) implies gamma-ray efficiencies ($\eta = 4\pi G_{100}d^2/\dot{E}$) of 300% (YMW16) or 60% (NE2001). While apparent efficiencies greater than 100% are possible if emission is strongly beamed towards us, the YMW16 value is far higher than any seen in a gamma-ray MSP (Smith et al. 2023), while the NE2001 value is more plausible, but still on the higher end of the population.

Alternatively, the non-detection of gamma-ray pulsations from this pulsar, and also from PSR J1823-3543, could be due to the pulsar not actually being the counterpart of the associated gamma-ray source. Kerby et al. (2023) searched for X-ray counterparts to our pulsars in Swift-XRT data, finding, for 4FGL J1906.4-1757 and 4FGL J1823.8-3544, alternative X-ray sources whose positions are inconsistent with our pulsar timing positions, but are still within the Fermi-LAT localisation regions. For 4FGL J1906.4-1757 this X-ray source was deemed likely to be a blazar, which could provide an alternative explanation for the gamma-ray source. However, this picture requires the serendipitous discovery of two unrelated MSPs within the targeted 4FGL gamma-ray localisation region, a scenario which is rather unlikely. From the GBNCC survey discoveries (McEwen et al. 2020), Bhattacharyya et al. (2022) estimated a density of 0.002 MSPs per square degree that could be discovered by the GBNCC survey. This value suggests the chance of us finding a single unrelated MSP among the 79 Fermi-LAT sources (with typical semi-major axes below 5') included in our survey was just 0.3%, although this number could be somewhat higher given the relative sensitivity of our searches.

Three pulsars, J1623-6939, J1823-3543, and J1858-5422, appear to have typical He WD companions, along with about half of the binary MSPs in the Galactic field. The first two have longer orbital periods (11 and 144 d, respectively) and a measurable eccentricity, while for the third, with an orbital period of 2.6 d, only an upper limit on the eccentricity ($e < 6 \times 10^{-6}$) is measurable. This is in line with the eccentricity-orbital period relation of Phinney & Kulkarni (1994) valid for NS-WD binary
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Phase-aligned radio (red curves) and gamma-ray (grey histograms) pulse profiles. The gamma-ray background level, estimated from the distribution of photon weights as $b = \sum_i w_i (1 - w_i) / n_{\text{bins}}$, is shown by dashed horizontal black lines. Radio pulse profiles are in arbitrary flux density units, scaled to match the same pulse amplitude as the gamma-ray peak, and with the baseline flux matching the gamma-ray background level.</figcaption></figure>

systems formed via stable mass transfer through Roche lobe overflow.

## 5.4. Radio-gamma ray phase alignment

The relative phase offset between radio and gamma-ray pulses encodes information about a pulsar's emission and viewing geometry (e.g. Johnson et al. 2014; Smith et al. 2023). For binary MSPs, assumed to have been spun-up via accretion from the companion star, and therefore spin-aligned to the orbit, this can also constrain the binary inclination angle (e.g. Corongiu et al. 2021). Gamma-ray pulsar emission models predict a dependency between the phase lag ($\delta_{g-r}$) between the radio and gamma-ray pulses and the separation ($\Delta$) between the two main peaks in the gamma-ray pulse profile, with the largest $\Delta$ values expected to occur for pulsars whose radio and gamma-ray peaks are closely aligned with each other. Figure 10 in 3PC demonstrates that this correlation is significantly detected in the known population of gamma-ray pulsars.

A full modelling treatment is beyond the scope of this work, but we show the phase-aligned radio and gamma-ray pulse profiles for the 6 gamma-ray detected pulsars in Fig. 8, with the gamma-ray photon phases rotated to align with the fiducial zero-phase points of the radio pulse profiles. We also reproduce the $\delta_{g-r}-\Delta$ from the 3PC plot here in Fig. 9, along with $\delta_{g-r}$ and $\Delta$ (for pulsars with more than one gamma-ray peak in their pulse profile) for the MSPs studied here extracted from the gamma-ray pulse profile templates used for timing in Sect. 3. We followed the same criteria used by 3PC to determine the first gamma-ray peak from which to measure $\delta_{g-r}$. For PSRs J1526-2744 and J1858-5422 this determination is possibly ambiguous as these pulsars exhibit two peaks almost half a rotation apart; however, with our chosen definition, their $\delta_{g-r}-\Delta$ values lie within the bulk of the known population. Four of these MSPs have gamma-ray

peaks that appear to arrive before the radio pulse ($\delta_{g-r} > 0.5$), a phenomenon often seen in MSPs but almost never in non-recycled pulsars.

## 6. Polarimetric results

In this section, we briefly describe the main polarimetric properties of our MSPs. Stokes parameters have been recorded only at MeerKAT through the Pulsar Timing User Supplied Equipment (PTUSE; Bailes et al. 2020) as the long integrations (hence, the large data size) were already prohibitive for the other telescopes. We therefore have polarimetry on all but two pulsars, J1526-2744 and J1803-6707, timed at Nançay and Parkes (MeerKAT data on these sources have been obtained, as part of the Fermi L-band Shallow survey, with the Accelerated Pulsar Search User Supplied Equipment (APSUSE) recording total intensity information only).

Polarisation calibration at MeerKAT is performed on-line as described in Serylak et al. (2021). All PTUSE data were then folded off-line using DSPSR<fnref n="9" /> (van Straten & Bailes 2011) and the best ephemeris for each pulsar. The pulsar software PSRCHIVE<fnref n="10" /> (Hotan et al. 2004) was used for all the subsequent steps: the routine *pac* was used to correct for parallactic angle variations and *pazi* was used to remove RFI. For each pulsar, cleaned, calibrated data were summed in phase with *psradd*, in order to obtain a high S/N profile. Finally *rmfit* was run to maximise the S/N of the linear polarisation component of the summed profiles and obtain a measurement of the rotation measure of each pulsar.

Figure 10 shows the calibrated, cleaned, de-rotated pulse profiles for the six pulsars (all except the isolated J1709-0333,

<footnote n="9">⁹ http://dspsr.sourceforge.net/</footnote>

<footnote n="10">¹⁰ http://psrarchive.sourceforge.net/</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Gamma-ray peak separations ($\Delta$) vs. radio/gamma-ray phase-lag $\delta_{g-r}$ for gamma-ray pulsars from 3PC (coloured markers) and the MSPs studied here (black stars). For two pulsars, the definition of the primary gamma-ray peak is potentially ambiguous, and so for these we show the alternative values with empty markers, linked by dashed lines. Pulsars with only one gamma-ray peak are shown below $\Delta = 0$ with a random scatter in the y-axis for clarity.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Polarisation profiles of the six Fermi MSPs with significant polarisation detections. The total intensity is shown in black, while the total linear and circular polarisation are shown in red and blue, respectively. The top panels for each pulsar show the polarisation position angle (plotted for the profile bins where the linear polarisation has a S/N larger than 3).</figcaption></figure>

which shows no detectable polarisation) obtained by summing all the observations collected with the MeerKAT L-band receiver (1.3 GHz). The black line is the total intensity profile, while red and blue are, respectively, the linear and circular polarisation profiles. The top part of each panel shows the variation of the polarisation position angle across the pulse rotational phase.

<table><caption>**Table 6.** Polarisation properties for six TRAPUM Fermi MSPs.</caption><thead><tr><th>Pulsar</th><th>%L</th><th>%V</th><th>%|V|</th><th>RM<br>(rad m<sup>-2</sup>)</th></tr></thead><tbody><tr><td>J1036-4353</td><td>11.1(13)</td><td>1.7(12)</td><td>1.7(12)</td><td>-18(8)</td></tr><tr><td>J1623-6936</td><td>36.9(16)</td><td>6.7(9)</td><td>7.7(9)</td><td>75(10)</td></tr><tr><td>J1757-6032</td><td>65.6(17)</td><td>-0.6(11)</td><td>5.1(10)</td><td>99(4)</td></tr><tr><td>J1823-3543</td><td>42.1(9)</td><td>-7.9(8)</td><td>9.1(8)</td><td>62(4)</td></tr><tr><td>J1858-5422</td><td>67.3(22)</td><td>0.0(17)</td><td>1.4(17)</td><td>38(11)</td></tr><tr><td>J1906-1754</td><td>30.4(11)</td><td>-2.5(8)</td><td>7.3(8)</td><td>-26(12)</td></tr></tbody></table>

**Notes.** Left to right we report: pulsar name, percentage of linear polarisation, percentage of circular polarisation and of its absolute value, and rotation measure.

Table 6 lists, for the six pulsars of Fig. 10, the percentage of linear (L) and circular (V) polarisation and its absolute value (|V|), and the rotation measure obtained from rmfit. The percentages of L, V, and |V| were calculated following Tiburzi et al. (2013).

## 7. Summary and conclusions

In this paper, we present the results of a multi-telescope, multi-wavelength timing campaign on the nine MSPs (eight binaries and one isolated) discovered as part of the TRAPUM Fermi Shallow survey at L-band (Clark et al. 2023). Preliminary ephemerides obtained from radio observations allowed us to detect gamma-ray pulsations for six of our targets. Thanks to this, we were able to extend their timing solutions over the 15-year Fermi-LAT data span. This was also possible for the RB pulsar J1803-6707, in spite of the presence of large orbital variations that are typical of RBs and that usually prevent the extrapolation of timing results beyond the radio data coverage; this was possible thanks to a new method (described in Sect. 4.2 and initially given in Thongmeearkom et al. 2024). For the three pulsars for which no gamma-ray pulsed emission was found, a fully coherent timing solution was achieved with a typical longer radio follow-up campaign.

To take advantage of the full potential of both the radio and gamma-ray data, for the six sources with pulsations detected in both bands, we developed a new joint timing technique (described in Sect. 4). It allowed us to obtain precise measurements for the long-term parameters constrained primarily by the gamma-ray photons and, in particular, to significantly measure the proper motion of four MSPs. The availability of well constrained, long- and short-term timing parameters is important not only to better characterise our targets, but also to extend their study at other wavelengths and through messengers other than the electromagnetic. Our timing results, for instance, are being currently used for a search of continuous gravitational waves over the entire LIGO data span (Ashok et al. in prep.).

The results of the (joint) timing campaign (presented in Sect. 5) show that our sample of pulsars differs from the previously known Galactic population (and, marginally, also from the earlier subset of Galactic MSPs discovered in association with Fermi sources); these differences are mainly seen in the distribution of the spin periods, peaking at significantly lower values, in our case. We note that the difference in spin-period distribution does not correspond to a difference in spin-down energy distribution. It is therefore unlikely that the smaller spin periods of our targets are related to a difference in the selected gamma-ray targets. It may rather simply depend on the higher
sensitivity of the MeerKAT telescope (allowing us to find fast-
spinning sources also with larger duty cycles) and on the conse-
quently shorter observing time required (allowing us to minimise
the pulse smearing induced by the pulsar accelerating in a binary
system). Shorter spin periods imply, potentially, smaller errors
in the ToAs; shorter period pulsars have, therefore, a higher
potential for being good candidates for high precision experi-
ments. One example is the pulsar timing array experiments for
the detection of nano-hertz gravitational waves, where MeerKAT
promises to be a powerful asset.

Individual sources of interest reported in this paper include
the two eclipsing RBs J1036-4353 and J1803-6708, whose
timing results are being used to complement optical observa-
tions (Phosrisom et al. in prep.) and will allow us to constrain
the masses of the NSs in the systems, adding valuable insights
into the NS mass distribution and potentially supplementing our
understanding of the EoS of nuclear matter. Similarly, albeit
through radio timing-only, the high-mass binary containing PSR
J1757-6032, has also demonstrated potential for precise NS
mass measurements through the detection (currently with lim-
ited significance) of the Shapiro delay relativistic effect. This
pulsar appears to be one of the few known objects to have been
formed through a Case A Roche lobe overflow: an evolutionary
path possibly leading to the formation of massive NSs. Further
dedicated observations will be needed to fully characterise this
system.

Building on our results and on the methods developed to
fully exploit the potential of our multi-wavelength data set, we
have now expanded our searches of *Fermi* targets both in terms
of frequency (observing the same sources covered in our L-
band survey also at 800 MHz; Thongmeearkom et al. in prep.)
and source sample (both new gamma-ray-only-selected targets
and spider candidates selected through multi-wavelength infor-
mation; Thongmeearkom et al. 2024). The preliminary results
once again confirm the high scientific potential of pulsar experi-
ments targeting *Fermi* unidentified point-like sources, especially
when coupled to such a sensitive and versatile instrument as the
MeerKAT telescope.

Data availability

The full ORBIFUNC and ORBWAVES ephemeris for PSR
J1803-6708 can be downloaded from Zenodo at the following
link: https://zenodo.org/records/13862732

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. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time space reference systems department of the Paris Observatory. TRAPUM observations use the FBFUSE and APSUSE computing clusters for data acquisition, storage and analysis. These clusters were funded, installed and operated by the Max-Planck-Institut für Radioastronomie and the Max-Planck-Gesellschaft. PTUSE was developed with support from the Australian SKA Office and Swinburne University of Technology. Murriyang, 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. This publication is partly based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. 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. This work was supported in part by the "Italian Ministry of Foreign Affairs and International Cooperation", grant number ZA23GR03, under the project "RADIOMAP-Science and technology pathways to MeerKAT+": the Italian and South African synergy" MBu and AP acknowledge the use of resources

from the research grant "iPeska" (P.I. Andrea Possenti) funded under the INAF
national call Prin-SKA/CTA (approved with the Presidential Decree 70/2016),
and the use of the TRG computer cluster at INAF - Cagliari, funded by the
Autonomous Region of Sardinia (Regional Law 7 August 2007 n. 7, year 2015,
"Highly qualified human capital"; P.I. Marta Burgay). R.P.B. acknowledges sup-
port from the European Research Council (ERC) under the European Union's
Horizon 2020 research and innovation programme (grant agreement No. 715051;
Spiders) The National Radio Astronomy Observatory is a facility of the National
Science Foundation operated under cooperative agreement by Associated Uni-
versities, Inc. SMR is a CIFAR Fellow and is supported by the NSF Physics
Frontiers Center award 2020265. 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 scien-
tific data analysis. These include the National Aeronautics and Space Adminis-
tration 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. We gratefully acknowledge Thankful H. Cro-
martie, Matthew Kerr, David A. Smith and David J. Thompson for reviewing
this manuscript on behalf of the Fermi-LAT collaboration.

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 839, 12

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53

Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 223, 26

Ackermann, M., Albert, A., Atwood, W. B., et al. 2014, ApJ, 793, 64

Agazie, G., Alam, M. F., Anumarlapudi, A., et al. 2023a, ApJ, 951, L9

Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023b, ApJ, 951, L8

Applegate, J. H., & Shaham, J. 1994, ApJ, 436, 312

Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Science, 324, 1411

Atwood, W., Albert, A., Baldini, L., et al. 2013, in Proceedings of the 4th Fermi Symposium, Monterey, California, eds. T. J. Brandt, N. Omodei, C. Wilson-Hodge, et al., 2012, 8

Bailes, M., Jameson, A., Abbate, F., et al. 2020, PASA, 37, e028

Ballet, J., Bruel, P., Burnett, T. H., Lott, B., & The Fermi-LAT collaboration 2023, ArXiv e-prints [arXiv:2307.12546]

Bassa, C. G., Pleunis, Z., Hessels, I. W. T., et al. 2017, ApJ, 846, L20

Bezuidenhout, M. C., Clark, C. J., Breton, R. P., et al. 2023, RAS Tech. Instrum., 2, 114

Bhattacharyya, B., Roy, J., Freire, P. C. C., et al. 2022, ApJ, 933, 159

Bickel, P., Kleijn, B., & Rice, J. 2008, ApJ, 685, 384

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

Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, ArXiv e-prints
[arXiv:1810.11394]

Clark, C. J., Nieder, L., Voisin, G., et al. 2021, MNRAS, 502, 915

Clark, C. J., Breton, R. P., Barr, E. D., et al. 2023, MNRAS, 519, 5590

Cordes, J. M., & Lazio, T. J. W. 2002, ArXiv e-prints [arXiv:astro-ph/0207156]

Corongiu, A., Mignani, R. P., Seyffert, A. S., et al. 2021, MNRAS, 502, 935

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

Dewey, R. J., Taylor, J. H., Weisberg, J. M., & Stokes, G. H. 1985, ApJ, 294, L25

Dodge, O. G., Breton, R. P., Clark, C. J., et al. 2024, MNRAS, 528, 4337

EPTA Collaboration, & InPTA Collaboration (Antoniadis, J., et al.) 2023, A&A, 678, A50

FERMI-LAT Collaboration (Ajello, M., et al.) 2022, Science, 376, 521

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

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

Freire, P. C. C., & Wex, N. 2010, MNRAS, 409, 199

Freire, P. C. C., & Wex, N. 2024, Living Reviews in Relativity, 27, 5

Freire, P. C., Kramer, M., & Lyne, A. G. 2001, MNRAS, 322, 885

Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Geyer, M., Venkatraman Krishnan, V., Freire, P. C. C., et al. 2023, A&A, 674, A169

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

Hobbs, G., Archibald, A., Arzoumanian, Z., et al. 2010, CQG, 27, 084013
Hobbs, G., Manchester, R. N., Dunning, A., et al. 2020, PASA, 37, e012

Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302

Hu, H., Kramer, M., Wex, N., Champion, D. J., & Kehl, M. S. 2020, MNRAS, 497, 3118

Johnson, T. J., Venter, C., Harding, A. K., et al. 2014, ApJS, 213, 6

Kalapotharakos, C., Wadiasingh, Z., Harding, A. K., & Kazanas, D. 2023, ApJ, 954, 204

Kansabanik, D., Bhattacharyya, B., Roy, J., & Stappers, B. 2021, ApJ, 920, 58

Kerby, S., Falcone, A. D., & Ray, P. S. 2023, ApJ, 954, 94

Kerr, M. 2011, ApJ, 732, 38

Kerr, M. 2019, ApJ, 885, 92

Kerr, M., Ray, P. S., Johnston, S., Shannon, R. M., & Camilo, F. 2015, ApJ, 814, 128

Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2021, Phys. Rev. X, 11, 041050

Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326, 274

Lazaridis, K., Verbiest, J. P. W., Tauris, T. M., et al. 2011, MNRAS, 414, 3134

Lorimer, D. R., & Kramer, M. 2004, Handbook of Pulsar Astronomy (Cambridge University Press), 4

Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45

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

Martinez, J. G., Stovall, K., Freire, P. C. C., et al. 2015, ApJ, 812, 143

McEwen, A. E., Spiewak, R., Swiggum, J. K., et al. 2020, ApJ, 892, 76

McMillan, P. J. 2017, MNRAS, 465, 76

Nice, D., Demorest, P., Stairs, I., et al. 2015, Astrophysics Source Code Library [record ascl:1509.002]

Nieder, L., Clark, C. J., Bassa, C. G., et al. 2019, ApJ, 883, 42

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

Nieder, L., Clark, C. J., Kandel, D., et al. 2020b, ApJ, 902, L46

Özel, F., & Freire, P. 2016, ARA&A, 54, 401

Pan, Z., Lu, J. G., Jiang, P., et al. 2023, Nature, 620, 961

Pétri, J., & Mitra, D. 2021, A&A, 654, A106

Phinney, E. S., & Kulkarni, S. R. 1994, ARA&A, 32, 591

Pletsch, H. J., & Clark, C. J. 2014, ApJ, 795, 75

Polzin, E. J., Breton, R. P., Bhattacharyya, B., et al. 2020, MNRAS, 494, 2948

Ransom, S. 2011, Astrophysics Source Code Library [record ascl:1107.017]

Ransom, S. M., Eikenberry, S. S., & Middleditch, J. 2002, AJ, 124, 1788

Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17

Ray, P. S., Abdo, A. A., Parent, D., et al. 2012, ArXiv e-prints [arXiv:1205.3089]

Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6

Roberts, M. S. E. 2012, Proc. Int. Astron. Union, 8, 127

Serylak, M., Johnston, S., Kramer, M., et al. 2021, MNRAS, 505, 4483

Shamohammadi, M., Bailes, M., Freire, P. C. C., et al. 2023, MNRAS, 520, 1789

Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789

Shklovskii, I. S. 1970, Sov. Ast., 13, 562

Smith, D. A., Bruel, P., Cognard, I., et al. 2019, ApJ, 871, 78

Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191

Splaver, E. M., Nice, D. J., Arzoumanian, Z., et al. 2002, ApJ, 581, 509

Su, M., Slatyer, T. R., & Finkbeiner, D. P. 2010, ApJ, 724, 1044

Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928

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

Thongmeearkom, T., Clark, C. J., Breton, R. P., et al. 2024, MNRAS, 530, 4676

Tiburzi, C., Johnston, S., Bailes, M., et al. 2013, MNRAS, 436, 3557

van Straten, W., & Bailes, M. 2011, PASA, 28, 1

Venter, C., Johnson, T. J., & Harding, A. K. 2012, ApJ, 744, 34

Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29
## Appendix A: Additional figures

<figure><img src="image_12.png" /><figcaption>Fig. A.1. Corner plot for all the jointly fitted parameters for J1526-2744 (colours as in Fig. 1).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. A.2. Corner plot for all the jointly fitted parameters for J1623-6936 (colours as in Fig. 1).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. A.3. Corner plot for all the jointly fitted parameters for J1709-0333 (colours as in Fig. 1).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. A.4. Corner plot for all the jointly fitted parameters for J1757-6032 (colours as in Fig. 1).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. A.5. Corner plot for all the jointly fitted parameters for J1858-5422 (colours as in Fig. 1).</figcaption></figure>