A targeted radio pulsar survey of redback candidates with MeerKAT

Redbacks are millisecond pulsar binaries with low mass, irradiated companions. These systems have a rich phenomenology that can be used to probe binary evolution models, pulsar wind physics, and the neutron star mass distribution. A number of high-confidence redback candidates have been identified through searches for variable optical and X-ray sources within the localisation regions of unidentified but pulsar-like Fermi -LAT gamma-ray sources. However, these candidates remain unconfirmed until pulsations are detected. As part of the TRAPUM project, we searched for radio pulsations from six of these redback candidates with MeerKAT. We discovered three new radio millisecond pulsars, PSRs J0838 − 2527, J0955 − 3947 and J2333 − 5526, confirming their redback nature. PSR J0838 − 2827 remained undetected for two years after our discovery despite repeated observations, likely due to evaporated material absorbing the radio emission for long periods of time. While, to our knowledge, this system has not undergone a transition to an accreting state, the disappearance, likely caused by extreme eclipses, illustrates the transient nature of spider pulsars and the heavy selection bias in uncovering their radio population. Radio timing enabled the detection of gamma-ray pulsations from all three pulsars, from which we obtained 15-year timing solutions. All of these sources exhibit complex orbital period variations consistent with gravitational quadrupole moment variations in the companion stars. These timing solutions also constrain the binary mass ratios, allowing us to narrow down the pulsar masses. We find that PSR J2333 − 5526 may have a neutron star mass in excess of 2 M ⊙ .


INTRODUCTION
Spider pulsars are a family of binary systems in which a millisecond pulsar (MSP) orbits with a tidally-locked, low-mass degenerate or semi-degenerate companion in a tight orbit (∼few hours up to a day; e.g.Roberts et al. 2017).Their population is generally subdivided into two main categories, black widows and redbacks, according to the mass of their companion (≪ 0.1M ⊙ and 0.2 − 0.4M ⊙ , respectively; Roberts 2012).In these systems, the two stars interact via ★ E-mail: tinn.thongmeearkom@postgrad.manchester.ac.uk an intra-binary shock between the pulsar's relativistic wind and the companion stellar wind.The very energetic pulsar is responsible for the strong irradiation of the companion whereby the side that faces the pulsar can be heated to several thousands of degrees above the opposite side, thus creating strong variability at visible wavelengths (e.g.Callanan et al. 1995;Breton et al. 2013).Furthermore, in a large fraction of systems, material blown away from the companion's surface by the pulsar wind can obstruct the radio emission from the pulsar (most likely via synchrotron absorption) and cause eclipses that can extend for a significant fraction of the orbit (e.g.Polzin et al. 2018Polzin et al. , 2019)).Spider companions therefore 'evaporate' over time, though current mass-loss rate estimates (e.g.Polzin et al. 2018) suggest that this process is not sufficient to provide a formation channel for isolated MSPs, as initially proposed by Fruchter et al. (1988).
The evolutionary history of spider pulsars leading to their current state is still poorly understood, but they most certainly experience a phase of mass transfer from the companion, which pumps angular momentum into the pulsar and where the system would be visible as an X-ray binary (Alpar et al. 1982;Radhakrishnan & Srinivasan 1982;Bhattacharya & van den Heuvel 1991).This evolutionary link was clearly evidenced with the discovery of the accreting MSPs (Wĳnands & van der Klis 1998;Chakrabarty et al. 2003) and then of the transitional millisecond pulsars in which X-ray binaries have been seen to turn into redback pulsars (Archibald et al. 2009) and vice versa (Papitto et al. 2013;Stappers et al. 2014).More recently, it became clearer that irradiation effects play a crucial role in dictating their evolution pathway, with the irradiation efficiency potentially determining if a system turns into a black widow or a redback (e.g.Chen et al. 2013;De Vito et al. 2020).
Multi-wavelength studies make it possible to constrain orbital parameters and masses of spider pulsar systems.Most notably, pulsar timing in the radio and gamma-ray regimes enables one to measure the five Keplerian parameters of the pulsar with great accuracy.Optical photometric light curves, on the other hand, can be modelled to derive the orbital inclination and, if combined with radial velocity information from spectroscopy of the companion, they can be used to infer the mass of each component (Lorimer & Kramer 2004).Evidence showing that spider systems have higher-than-average neutron star masses (1.46 ± 0.30 M ⊙ for an average neutron star; Zhang et al. 2011) has now been found, with a report that the neutron stars in redbacks have a median mass of 1.78 ± 0.09M ⊙ (Strader et al. 2019).This motivates searches for > 2M ⊙ neutron stars, which could constrain the equation of state.-some interesting examples above 2M ⊙ have possibly been found (Linares 2020;Romani et al. 2022).
Multiple approaches have been pursued to find new spider pulsars.Unguided widefield radio pulsar survey has only found a handful of the ∼ 80 known spiders thus far (the ATNF Pulsar Catalogue; Manchester et al. 2005), including the original black widow (Fruchter et al. 1988).This technique is comparatively less efficient at finding spider pulsars than other types of binary pulsars, primarily since radio eclipses and the relatively large and rapidly changing accelerations heavily hinder the detection of spiders.Instead, the bulk of the spider population has been uncovered through the targeted approach of pointing at candidates identified via multi-wavelength campaigns.By far the most successful approach to finding new Galactic spiders has been to pursue their gamma-ray signature.In 2008, the Fermi Gamma-ray Space Telescope started to operate.One of the telescope's instruments, the Large Area Telescope (LAT; Atwood et al. 2009), is especially suitable for detecting the pulsed GeV gammaray emission from millisecond pulsars (Abdo et al. 2009), which are particularly efficient high-energy emitters (Smith et al. 2023).The Fermi-LAT has detected and localised a large number of point sources.The latest catalogue, the Fermi-LAT Fourth Source Catalogue (Abdollahi et al. 2020) Data Release 4 (4FGL-DR4, Ballet et al. 2023) detecting more than 7000 sources, while the Third Fermi-LAT Gamma-ray Pulsar Catalogue contains nearly 300 confirmed gamma-ray pulsars (Smith et al. 2023).
Point-like gamma-ray sources that have not been associated with other objects at other wavelengths -of which thousands remain -provide a very effective list of locations that can be searched for energetic pulsars in the radio where they are bright enough to yield a detection over a short minutes-to-hour observation.Gamma-ray point source characteristics such as variability and spectral properties can also be used to identify sources most likely to be a pulsar (e.g.Saz Parkinson et al. 2016).Multiple surveys, initially coordinated through the Fermi Pulsar Search Consortium (PSC; Ray et al. 2012), have successfully employed this strategy, including studies conducted with the Green Bank Telescope (GBT) (Ransom et al. 2011), the Murriyang Parkes telescope (e.g.Keith et al. 2011;Camilo et al. 2015), the Arecibo telescope (Cromartie et al. 2016) and the Five-hundred-meter Aperture Spherical radio Telescope (FAST) (Wang et al. 2021).Given that an important limitation to finding pulsars is instantaneous sensitivity, the MeerKAT radio telescope has recently provided us with a new, improved tool for this quest, as was demonstrated by the first results of the TRansients and PUlsars with MeerKAT (TRAPUM) survey of unidentified Fermi-LAT sources, which discovered 9 new pulsars from a search of 79 unidentified gamma-ray sources (Clark et al. 2023b).
A similar strategy to the one laid out above is to pursue candidate spider pulsars identified through searches for periodic optical and X-ray counterparts, often confirmed with spectral velocity measurements indicating an unseen neutron-star-mass primary (e.g.Kong et al. 2012;Romani 2012).These candidates are usually identified through the monitoring of suitable Fermi gamma-ray unidentified point sources or serendipitous observations of such fields made by large-scale surveys.More than a dozen candidate redback and blackwidow binaries have been identified in this way (see Table 15 of Smith et al. 2023), but only a handful of these have been confirmed by a pulsation detection.Not including the discoveries that we present here, just three optically-identified redback candidates, PSRs J2339−0533 (Ray et al. 2020), J0212+5320 (Perez et al. 2023) and J1910−5320 (Au et al. 2023;Dodge et al. 2024) were confirmed via radio pulsations, while a further three binaries were confirmed by direct searches for gamma-ray pulsations in the Fermi-LAT data (Pletsch et al. 2012;Nieder et al. 2020b;Clark et al. 2021), using constraints from optical monitoring to reduce the orbital parameter space.The long integration time required to accumulate enough gamma-ray photons makes this latter method a challenging computational task (Nieder et al. 2020a), and in many cases the available orbital constraints are not sufficiently precise to make a gamma-ray pulsation search feasible, leaving radio searches as the only means to confirm their MSP nature.However, only two of these three gamma-ray discovered spiders were eventually detected by folding radio observations using the gamma-ray timing ephemeris (Ray et al. 2013;Corongiu et al. 2021), and one of these, PSR J1311−3430, is only very sporadically detectable.The third of these systems, PSR J1653−0158, remains undetected in radio despite efforts by many large radio telescopes (Nieder et al. 2020b), suggesting that some spider binaries may be nearly always enshrouded by evaporated material.Hence, both radio and gamma-ray searches of these candidate systems remain important complementary approaches for obtaining a more complete view of this population.
Despite the rather low yield to date, this refined strategy of targeting optically-identified candidate binaries provides a number of potential advantages over more traditional surveys.First, the optical signature provides a much stronger prior as to the existence of a spider pulsar than the gamma-ray properties.Second, the localisation of the candidate is generally much better (sub-arcsecond vs. arcminute uncertainties), thus allowing one to position the radio beam directly at the right location to avoid sensitivity losses that arise towards the edge of the beam.This is particularly crucial for MeerKAT observations, where hundreds of coherent tied-array beams with widths of a few arcsecs are required to tile the error box of an unassociated Fermi source (which typically have widths measured in arcminutes), with typical sensitivity losses of around 50% occurring midway between beams.The ensuing search for radio pulsations can also be operated on a smaller data set, which opens possibilities for longer pointings and searching at higher accelerations when processing resources are limited, as they are for regular TRAPUM observations.Finally, the optical observations also provide an approximate orbital period and reference phase for the system, which in turn can be used to schedule the radio searching to occur around the inferior conjunction of the pulsar when it is least likely to be eclipsed.
In this paper, we expand the quest to find more spider pulsars by focusing on six Fermi unidentified gamma-ray sources, which have been found to have possible counterparts with redback-like optical light curves and/or modulated X-rays.Since these systems are likely to harbor a radio pulsar, we conducted a deep survey of these targets with the MeerKAT telescope, combining observations made at Lband and Ultra High Frequency (UHF).The structure of the paper is as follows.Section 2 presents the survey strategy and targets.Section 3 reports on the discoveries, gamma-ray searches and timing.Section 4 discusses the characteristics of the discovered pulsars.Section 5 draws some conclusions and summarises the work.

SURVEY AND OBSERVATION
The MeerKAT telescope is an interferometer array that contains 64 radio dishes with a diameter of 13.5 metres.Located in the Karoo, a desert region in the Northern Cape province of South Africa, MeerKAT is the most sensitive telescope operating at centimetre wavelengths in the Southern Hemisphere.The location was chosen due to its low population density.Consequently, the site is only affected by very low levels of radio frequency interference (RFI) and is high quality for radio astronomy (Booth et al. 2009).At the time the survey was conducted, two receivers were available: an L-band receiver operating between 856-1712 MHz, and a UHF receiver operating between 544-1088 MHz (see Jonas & MeerKAT Team 2016, for more technical information).
Six targets have been observed as part of this survey (see Table 1).Multi-wavelength observations indicated they are very likely redback pulsar binaries.The Fermi-LAT detected point sources in the direction of these sources that are not clearly associated to other known objects (Abdollahi et al. 2020).Optical and X-ray searches of the few arc-minute localisation regions were conducted to shed light on their nature, and identified optical and/or X-ray counterparts with lightcurves displaying the typical orbital periodicity of redback systems (e.g.Halpern et al. 2017b;Li et al. 2018;Swihart et al. 2020).The variability in the optical data can have either two maxima per orbital cycle, due to tidal distortion of the companion star, or one maximum per cycle if the pulsar's irradiation of the inner facing hemisphere of the companion star causes a large temperature difference.The X-ray modulation, on the other hand, is caused by an intra-binary shock and displays a double-peak pattern due to the emission from particles accelerated to high Lorentz factors along the shock boundary.Mass estimates derived from the projected radial velocity curve of the companions obtained from optical spectroscopy further indicate that these systems contain an unseen primary object with a mass compatible with that of a neutron star though no pulsations have been reported from these candidates until now (e.g.Rea et al. 2017;Strader et al. 2019).
In light of the strong evidence provided by the combination of gamma-ray, X-ray and optical data, there is a very high likelihood that all six targets are redback binaries that harbour a millisecond radio pulsar.Some of these sources have been targeted several times in radio pulsation searches with other telescopes without success (e.g.Camilo et al. 2015).These non-detections could be due to the intrinsic faintness of the pulsars, but are may also be due to the eclipsing nature of these redback systems, many of which are only sporadically detectable (e.g.PSR J0212+5321 has been detected only twice despite years of searching; Perez et al. 2023).Therefore, we decided to pursue a deep radio search with MeerKAT in an attempt to detect radio pulsations from the pulsars that they may host.2.1.1 4FGL J0523.3−2527Strader et al. (2014) suggested that this LAT source is a new probable gamma-ray pulsar binary.This system has one of the largest inferred companion masses among known or candidate redbacks (Strader et al. 2019), with a companion mass between 0.8−1.3M ⊙ assuming a neutron star mass in the range 1.4−2.0M ⊙ .The companion has been classified as a late G or early K star using optical spectroscopy (Strader et al. 2014).In addition to the projected radial velocity amplitude, rotational broadening was measured from the optical spectral lines, thus enabling a direct estimate of the mass ratio ( = 0.61 ± 0.06) leading to the above companion mass estimate.The radial velocity measurement also shows evidence of eccentricity ( = 0.04) which would be the largest known in a redback system.Later, Halpern et al. (2022) reported the most luminous optical and X-ray flares seen in a non-accreting pulsar so far.Previous attempts to find radio pulsations from this source have been unsuccessful (Guillemot et al. 2012;Petrov et al. 2013).

Targeted sources/ Redback Candidates
2.1.24FGL J0838.7−28274FGL J0838.7−2827 was identified as a high-confidence MSP candidate following an optical and X-ray spectral study, but with an unknown orbital period due to the limited duration of the observations (Halpern et al. 2017a;Rea et al. 2017).The likely redback nature was confirmed later from improved photometry and spectroscopy, with the optical companion being a low-mass M dwarf star in a 5.15 hr orbit around an unseen primary with a mass consistent with a neutron star (Halpern et al. 2017b).Flaring was detected similar to that seen in transitional MSPs, though the X-ray and gamma-ray luminosities suggest it is a non-accreting MSP binary (Halpern et al. 2017a, e.g. X−ray ∼ 10 31 erg s −1 vs.   ∼ 10 33 erg s −1 ).The gamma-ray flux also shows no significant variability indicative of a transition over the 14-year Fermi-LAT data set (Ballet et al. 2023).The optical spectrum shows variable H emission, which is believed to come from the wind driven by the heated side of the companion as the emission line is too broad to be chromospheric in origin (Halpern et al. 2017b).Finally, The Australia Telescope Compact Array (ATCA) observations of the field put a 3 upper limit of ∼20 Jy at 5.5 and 9 GHz on the presence of a point source (Rea et al. 2017).

4FGL J0940.3−7610
4FGL J0940.3−7610 was proposed as a redback candidate by Swihart et al. (2021) based on the discovery of a variable optical and X-ray counterpart to the Fermi point source.The optical light curve displays both ellipsoidal variations and irradiation with amplitudes typical of redback systems, with an orbital period of 6.5 hours.Orbital phaseresolved spectroscopy and modelling of the light curve suggest a lower mass neutron star (1.2 − 1.4 M ⊙ ) and higher mass companion ( c ≳ 0.4 M ⊙ ) compared to what is typically seen in redbacks, however, without a pulsation detection or a reliable ephemeris to Table 1.List of redback candidates in our survey, with parameter constraints from optical observations from the references.Positions are from Gaia Data Release 3 (Gaia Collaboration et al. 2016, 2023). asc is the epoch of the ascending node of the pulsar in the Barycentric Dynamical Time (TDB) system. b is the orbital period.1- uncertainties on the final digits of the orbital parameters are quoted in parentheses.The DM columns represent the predicted dispersion measures along these lines-of-sight at the distances estimated from previous studies, according to the NE2001 (Cordes & Lazio 2002) and YMW16 (Yao et al. 2017) Galactic electron density models.The dagger sign ( †) represents the maximum predicted DM along the lines-of-sight from both models.DM search are the DM ranges that we used for searching.2018) via optical spectroscopy and X-ray variability.In addition, a 77 mJy pulsar-like radio counterpart was found by Frail et al. (2016) in the 150 MHz TGSS ADR catalogue (Intema et al. 2017), but was discarded as it lies outside the 3FGL Fermi localisation error (Acero et al. 2015) -this source is, in fact, coincident with the optical counterpart, which lies well within the improved 4FGL localisation region (Li et al. 2018).However, repeated radio follow-up observations by Camilo et al. (2015) with Parkes, prior to the identification as a redback candidate, detected no radio pulsations, despite this source lying within the beam.

4FGL name
2.1.54FGL J1120.0−22044FGL J1120.0−2204 was the second brightest unclassified source in the 4FGL DR3 catalogue (Abdollahi et al. 2022).The presumed optical counterpart to the Fermi source, found by Swihart et al. (2022a), features a binary companion with a mass around 0.17 M ⊙ in a sub-day orbital period ( b ∼15.1 hr) around an unseen neutron star.As suggested by these authors, this system might be at an intermediate stage leading the companion to form an exceptionally low-mass helium white dwarf orbiting a neutron star, e.g.akin to PSR J1012+5307 (Mata Sánchez et al. 2020).Notably, optical observations have revealed a lack of variability in the system's behavior, suggesting the orbit might be seen with a face-on orientation.Furthermore, 4FGL J1120.0−2204displays a lower X-ray luminosity when compared to most redback systems.While it may not fit the conventional criteria for a redback candidate, its exceptional properties have led us to include 4FGL J1120.0−2204 in our survey, recognising its importance for binary evolution and bridging the gap between spider and white dwarf pulsar systems.
No radio pulsar has been detected within this system despite extensive radio searches conducted with various telescopes, including Effelsberg (Barr et al. 2013), GMRT (Bhattacharyya et al. 2021) and unpublished PSC searches with Parkes and GBT, while Hui et al. (2015) inspected nearby radio sources from two catalogues, the Molonglo Sky Survey (Bock et al. 1999) and the NRAO VLA Sky Survey (Condon et al. 1998).

4FGL J2333.1−5527
Our final candidate, 4FGL J2333.1−5527,exhibits an orbital period of approximately 6.9 hr and features a companion in a nearly edge-on orbit (Swihart et al. 2020).The detection of a significant gamma-ray eclipse confirms this picture (Clark et al. 2023a).Optical photometry reveals a typical redback light curve, while optical spectroscopy indicates the presence of a mid-K type companion.Furthermore, the study suggests that the neutron star mass in this system likely exceeds 1.4 M ⊙ , but it remains poorly constrained.No previous radio observations of 4FGL J2333.1−5527have been published.

Survey configuration
We designed our campaign to observe each target for one hour as this integration time should provide a significant sensitivity improvement over previous efforts, yet a duration between 5 − 20% of the orbital periods means that pulsations can still be recovered from a search of the full time series, using only one or two derivatives to account for the orbital motion, i.e. the so-called acceleration (Johnston & Kulkarni 1991) and jerk searches (Andersen & Ransom 2018), respectively.Each target was observed twice at both L-band and UHF.These two bands provide complementary information; on the one hand, pulsars often display steep radio spectra which favour low-frequency observations (Jankowski et al. 2018), while on the other hand, radio eclipses tend to be shorter and less opaque at higher frequencies (e.g.Polzin et al. 2020).Similarly, repeated observations can help mitigate effects such as scintillation due to the scattering of the pulse in the interstellar medium (e.g.Camilo et al. 2015).Repeats avoid eclipses due to slightly different orbital phase coverage, and because eclipses, especially in their non-central part, vary greatly.Finally, our knowledge of the orbital periods from optical follow-up lets us plan observations when the pulsar is nearer inferior conjunction (i.e. between phases 0.5 − 1.0 in the pulsar timing convention) where it is less likely to be eclipsed.
Note that 4FGL J1120.0−2204 is an exception to the above strategy due its late addition to the list.We have only collected data from one epoch at each frequency at the time of writing this paper.
All 64 dishes were requested for our observations to maximise sensitivity.However, the number of dishes varied from one observation to another depending on the availability at that time (see Table A1  and A2).Given that the optical counterparts provide accurate source positions, we could point a single tied-array beam at the candidates without worrying about beam size.We recorded the data using 4096 frequency channels and a time resolution of 38 s and 60 s at L-band and UHF, respectively.This combination of frequency and time resolution maximises fast MSP detectability and mitigates the effect of the unknown dispersion measure (DM) along the line of sight, which can smear pulses within single channels, a particularly detrimental effect at short spin periods.

Search tools and Acceleration search
We employed three software packages: PRESTO1 , Pulsar_miner2 and Spider_twister3 .PRESTO (PulsaR Exploration and Search TOolkit) is software for analysis and pulsar searching (Ransom et al. 2003).PULSAR_MINER is a python-based pipeline that automates all the steps of a PRESTO-based acceleration/jerk search (DDplan for DM steps, de-dispersion, searching, candidate filtering and folding), while Spider_twister automatically searches in orbital phase with several trials of epoch of ascending node,  asc , to find the value that produces the highest signal-to-noise folded result using a preliminary ephemeris as a starting point.
In PRESTO the acceleration search range is parameterised by , the (dimensionless) number of Fourier bins that a signal is smeared over due to orbital motion.This is related to the pulsar's acceleration, where  is the pulsar's fundamental spin frequency,  is the speed of light,  is the observation duration, ℎ is the harmonic number (such that the fundamental harmonic has ℎ = 1) and  is the acceleration.A similar dimensionless parameter, , quantifies the frequency smearing due to jerk, which is the first time derivative of the acceleration (Andersen & Ransom 2018).
Searching a range of non-zero jerks provides additional sensitivity to short period binaries (in the range  ≈ 0.1 b to  ≈ 0.2 b ), but requires more computational power than a pure acceleration search.In our study, we performed acceleration searches up to  = 200 for all targets, and jerk searches up to  = 600 for those with no clear detections in the acceleration search.Assuming a typical MSP spin period ( = 2 ms), these searches cover accelerations of 9.26 m s −2 , 37.03 m s −2 and 333.33 m s −2 for 1 hr, 30 min and 10 min segments, respectively.We then calculated the maximum acceleration of 15 redbacks using the orbital parameters provided in Strader et al. (2019).Consequently, we found that the highest maximum acceleration from these 15 sources is 26 m s −2 while the mean and the median are 11 m s −2 and 8 m s −2 , respectively.We also limited the range of DMs to search by taking the larger of the maximum values predicted by the NE2001 (Cordes & Lazio 2002) and YMW16 (Yao et al. 2017) models along the line of sight with extra ranges (see Table 1) for each target.

RESULTS
Out of the six redback candidates we searched, we discovered a radio millisecond pulsar in three of them: PSRs J0838−2827, J0955−3947 and J2333−5526 in 4FGL J0838.7−2827,4FGL J0955.3−3949, and 4FGL J2333.1−5527,respectively.All of them have been discovered at L-band in 2021.Figure 1 shows their pulse profiles at L-band and at UHF. Accordingly, three new MSPs have spin periods of  ∼ 3.6, 2.0 and 2.1 ms, with accelerations ( ∼13.6, 10.8 and 17.0 m s −2 ) demonstrating their binary nature.Further observations (see below) confirmed that their orbital periods are consistent with the optical counterparts.The spin periods are right in the range expected for redback systems, where most pulsars in these systems are found with  ≲ 5 ms (Strader et al. 2019).
Following these discoveries we searched archival observations from Murriyang Parkes and Green Bank telescopes that were accessible to us for the three new pulsars.Knowing the DM and the pulse period, and later having a partial timing solution, enabled us to detect all of them in a number of observations, with non-detections generally believed to result from eclipses or scattering (see Section 3.4).Full details of the epochs and phase coverage for the four 1 hr observations and 5 min timing campaign observations conducted with MeerKAT as well as ancillary data obtained from other facilities are provided in Table A1.Examples of the discovery, the archival data and the detection with eclipses are shown in Figure 2, 3 and 4, respectively.

Timing
We undertook timing campaigns for each newly detected pulsar to obtain long-term phase-connected timing solutions that precisely describe their spin and orbital properties.These observations were performed at MeerKAT, using the Pulsar Timing User Supplied Equipment (PTUSE) backend (Bailes et al. 2020), with the 64 m Murriyang Parkes telescope, using the Ultra Wide-Frequency Low receiver (UWL; Hobbs et al. 2020) and the Medusa Backend (see Section 3.1 of Hobbs et al. 2020), and, for the northernmost and more elusive PSR J0838−2827, also the Nançay Radio telescope (NRT) located in France.
Initially, we adopted a pseudo-logarithmic cadence strategy using the Parkes telescope (as time to follow up these sources at MeerKAT was granted at a later time), with one or two observations taken on the first day, followed by single observations on the second, fifth, tenth days, and so on.Typically, these observations lasted 1 to 2 hr, depending on the target and, in case of a detection, provided one to four times-of-arrival (ToAs).
While PSR J0955−3949 was always detected, PSR J2333−5526 was successfully detected more sporadically and PSR J0838−2827 was never seen in our early UWL Parkes data, despite being visible in two archival observations obtained with a less sensitive backend.In October 2023, we then detected it again five times with Parkes.
A second pseudo-logarithmic campaign was hence carried out at MeerKAT on these two latter pulsars, with multiple 5 min observations taken on day 1, followed by single observations in the following days as described above.These allowed us, in the case of PSR J2333−5526, to get a good enough initial orbital solution, as described below, to properly fold previous Parkes data and, in some cases, recover the signal.PSR J0838−2827 remained undetected in our MeerKAT timing campaign until one detection was made in each of August, September and October 2023.
To follow up the sources over longer timescales, we continued our The horizontal axis displays the spin phase of each source, replicated for clarity, while the vertical axis depicts the intensity, normalised to the peak value.The grey line presents the background level, calculated as the median of the off-pulse region.Since the gamma-ray ephemerides do not perfectly fold radio data to perform an absolute alignment, we aligned the L-band (black) and the UHF (red) profiles using a cross-correlation.
timing campaign with approximately monthly observations using the Parkes telescope and NRT.We used the optical constraints given in Table 1 as starting points for solving the orbits and fitting a circular orbit model (using PRESTO's fit_circular_orbit.py) to the best-fitting spin frequency in each observation to obtain an initial estimate for the pulsar's orbital semi-major axis and to improve upon the optical estimates for  b and  asc .After obtaining an acceptable initial orbit, we used Dracula4 (Freire & Ridolfi 2018) to determine the global rotation count of these pulsars and derive fully phase-connected timing solutions.Several orbital frequency derivatives needed to be added at this stage to obtain flat timing residuals with a reduced chi-square close to one and account for variations in the orbital period, a behaviour commonly seen in redback binaries (e.g.Pletsch & Clark 2015;Deneva et al. 2016).
Phase-connected radio timing solutions were obtained in this way for two of the three pulsars, PSRs J0955−3947 and J2333−5526 (see Table 2).We used these as a basis to obtain a longer-term timing solution from the Fermi-LAT gamma-ray data in Section 3.2.PSR J0838−2827 has only been detected 14 times (see Section 3.4 and Figure 5).Due to the long time separation of PSR J0838−2827 between 2021 data and 2023 data, it was impossible to obtain a phase-connected timing solution covering all radio data.Therefore, we used only observations in 2023 for radio timing to create a starting point for gamma-ray timing (see Section 3.2) for this source.

Gamma-ray pulsation searches and timing
We used the radio timing solutions obtained in Section 3.1, to search for gamma-ray pulsations from each new pulsar in the Fermi-LAT data.We analysed SOURCE-class events using the "Pass 8" P8R3_SOURCE_V3 instrument response functions (Atwood et al. 2013;Bruel et al. 2018).We included photons according to the energy-dependent cuts that are described in Table 1 of Abdollahi et al. (2022).We used the gtsrcprob tool and the 4FGL-DR3 (Abdollahi et al. 2022) spectral and spatial models, with the gll_iem_v07.fitsGalactic interstellar (Acero et al. 2016) and iso_P8R3_SOURCE_V3_v1.txtisotropic diffuse models, to assign weights representing the probability of each photon being emitted by the pulsar (Kerr 2011), as opposed to nearby point sources or the diffuse background.
Significant gamma-ray pulsations, with weighted H-test values (Kerr 2011) of  = 71 and  = 45 were detected from PSRs J0955−3947 and J2333−5526, respectively, within the validity period of the preliminary radio timing ephemerides.However, these timing solutions did not reveal pulsations in earlier Fermi-LAT data.This is because the radio timing ephemerides for these pulsars begin at our discovery epochs, and contain several orbital frequency derivatives describing the orbital phase variations over time, and this Taylor series model extrapolates poorly to earlier times.For PSR J0838−2827, the initial radio ephemeris did not cover enough Fermi-LAT data to reveal significant pulsations, but it did provide precise measurements of the pulsar's semi-major axis, and ascending node epoch.With these values fixed, and astrometry fixed at the Gaia DR3 values for the companion star, we performed a search over the spin frequency, spin-down rate and orbital period.This led to a significant detection with  = 78, but the pulse profile was unclear in some parts of the data, again indicative of orbital phase variations not accounted for in our initial timing model.
To obtain a gamma-ray timing solution describing the entire 15year Fermi-LAT data set, we therefore first had to search over the orbital phase, period and the first orbital period derivative as a function of time, searching for coherent pulsations in overlapping 400-day intervals covering the LAT data.The results of these "sliding window" searches are shown in Figure 6, revealing shifts in the pulsar's ascending node time of over 100 s (for PSR J0955−3947) compared to a model with a constant orbital period.
We then used the resulting orbital-phase vs. time measurements as a starting point to fit a fully-coherent timing model.We used the framework introduced in Clark et al. ( 2021) in which we modeled the orbital phase variations as a stochastic Gaussian process, and jointly fit the parameters of the covariance function of this process along with the standard timing model parameters, and those of the template Table 2. Timing solutions for the newly discovered millisecond pulsars.Timing parameter values are from gamma-ray timing, described in Section 3.2 using priors on astrometry parameters from Gaia, and the radio timing value as a prior on  1 .References:  Halpern et al. (2017b),  Li et al. (2018),   pulse profile used to evaluate the likelihood of the photon-weighted pulse profile (Kerr et al. 2015).We use a Matérn covariance function, parameterised by an amplitude (ℎ), length scale (ℓ) and "smoothness" () hyperparameters.A process whose covariance function has these parameters has a smoothly broken power law power spectral density, which is flat below a cut-off frequency ∝ ℓ −1 , turning over to a power-law with spectral index  = −(2 + 1).
To make this procedure tractable, the method described in Clark et al. ( 2021) approximates the spin-phase likelihood function for each photon as a Gaussian, but we found that this approximation does not work well for fainter pulsars, or those with larger orbital phase variations.We therefore modified this procedure, using a Gibbs sampling method to deal with the multi-modal pulse profile.In this method, each photon is randomly assigned, according to its weight and an initial timing model, to either a single Gaussian pulse profile component (either one of two sharp peaks or to a bridge component connecting the two) or to the unpulsed background.The timing model (including the orbital phase variation Gaussian process) can then be fit analytically using generalised least squares, given these component assignments.The component assignments can then be updated according to a random sample drawn from the resulting posterior distribution on the timing model and template parameters, and this process is then repeated to generate a Monte-Carlo chain that approximates the full joint posterior on the timing model parame- ters, template pulse profiles, and Gaussian process hyperparameters.In the fitting, we adopted the Gaia DR3 (Gaia Collaboration et al. 2016, 2023) values and uncertainties as Gaussian priors on the astrometric parameters, and used radio timing estimates as priors on the pulsar semi-major axes.Gamma-ray timing is generally not precise enough to improve upon these measurements, with the exception of proper motion in PSR J2333−5526 and the semi-major axis of PSR J0838−2527, where the posterior uncertainties are slightly narrower.The estimated power spectral density of the orbital phase variations, and the best-fitting timing models from this procedure, are shown in Figures 6 and 7.
We checked the resulting gamma-ray timing solutions using the radio ToAs, which were generally in agreement, but revealed small residual trends (on the order of 1% to 10% of a rotation) that were consistent with residual variations in the orbital phase that are below the level of precision provided by the gamma-ray timing solution.Owing to these remaining residuals we have not attempted to phasealign the radio and gamma-ray pulse profiles.A modification of the above timing procedure to jointly fit radio ToAs and gamma-ray photon arrival times would likely resolve this.Clark et al. (2023a) detected evidence in the Fermi-LAT data for gamma-ray eclipses from both PSRs J0838−2827 and J2333−5526, but the eclipse properties were uncertain since no pulsar timing ephemeris was available to precisely determine the orbital period and phase parameters.No gamma-ray eclipse was detected from PSR J0955−3947.We confirm the detected gamma-ray eclipses in PSRs J0838−2827 and J2333−5526, and the non-detection for PSR J0955−3947, now using the precise gamma-ray timing solutions derived above to compute the orbital phases and to re-weight the photons according to their measured spin phase and best-fitting pulse profile template.For PSRs J0838−2827 and PSR J2333−5526, we find that the eclipse durations are 0.015 to 0.034 orbits, and 0.045 to 0.07 orbits, respectively (95% confidence intervals).Using the now-known mass ratio, and assuming that the companion star fills its Roche lobe, and that there is no significant gamma-ray absorption from intra-binary material outside the Roche lobe, this sets a lower bound on the binary inclination angle of  > 78 • for both pulsars.For PSR J0955−3947, the eclipse non-detection implies an inclination  < 76 • under the same assumptions.

Gamma-ray eclipses
These inclination constraints for PSRs J0955−3947 and J2333−5526 are also qualitatively consistent with the shapes of the gamma-ray pulse profiles, under the assumption that the pulsars achieved their observed spin rates through accretion from their companion stars such that their current spin axis is aligned with the orbit.For PSR J2333−5526, the two sharp peaks are very close to half a rotation apart, which is what one would expect from a gammaray pulsar viewed from close to the spin/binary equator (e.g.Watters et al. 2009).The narrower separation between peaks in the pulse profile for PSR J0955−3947 is similarly indicative of a more moderate viewing angle/inclination.Future modelling of the phase-aligned radio and gamma-ray pulse profiles (e.g.Johnson et al. 2014;Corongiu et al. 2021) may provide a quantitative test of these constraints, and also for PSR J0838−2827 whose broader pulse profile is somewhat more difficult to interpret.

Radio eclipses
We have detected radio eclipses in the data from PSRs J0838−2827 and J0955−3947.The orbital phase coverage of the various observations for all three new MSPs is shown in Figure 5 so that their eclipses can be contextualised.
Radio eclipsing is often seen in spider systems as gas outflowing from the companion star can block radio pulsations for extended parts of the orbit, generally near superior conjunction of the pulsar.Synchrotron absorption is established as the primary eclipse mechanism acting in most systems (Polzin et al. 2018) and whether the pulsar disappears completely or not, depends on the radio frequency compared to the characteristic plasma frequency of the eclipsing medium (Thompson et al. 1994).Additionally, the intervening plasma causes a local increase in the dispersion measure, which is seen as a delay curving the pulse trail.
The L-band discovery detection of PSR J0838−2827 in February 2021 displays an eclipse ingress around 20 min into the 1 hr observation followed by a total eclipse (see Figure 2).Another observation at UHF in June 2021 shows no sign of eclipse in another 1 hr long observation covering a similar range of the orbit.Since eclipses tend to last longer at lower frequencies, this clearly indicates that the medium responsible for the absorption is highly variable over time.The orbital range in which the eclipse took place, 0.7 − 0.8, is also unusual as this corresponds to the inferior conjunction of the pulsar (i.e.half an orbit away from where it normally happens) while NRT and Parkes timing observations show a typical eclipse around phase 0.4.An archival GBT observation covering orbital phases 0.5 − 0.6 also displayed a complete eclipse for part of the observation.Finally, Parkes archival data from 2017 shows a total eclipse for half of the observation in phases 0.25 − 0.40 (see Figure 3).For non-detection observations, we ran Spider_twister to check that the disappearances are due to a high orbital frequency derivative of the system with an inadequate preliminary ephemeris or actual eclipses.The latter seems to be the case as Spider_twister could not recover a pulsation from those observations.
For PSR J0955−3947, two MeerKAT observations display several short, transient eclipses lasting around 40 s near orbital phase 0.1 (see Figure 4).These broadband flux variations cannot be explained by interstellar scintillation, which for this pulsar would have a scintillation bandwidth of around 0.2 MHz according to NE2001 (Cordes & Lazio 2002).From the follow-up campaign, four of our Parkes observations starting close to orbital phase 0.4 show the pulsar coming out of the eclipse.While by design we have avoided observing near superior conjunction of the pulsar (phase 0.25), it appears quite likely that PSR J0955−3947 experiences "classic" eclipses between phases 0.1 and 0.4, with potentially short eclipses on the ingress side.
No eclipses were detected in the case of PSR J2333−5526.This is partly due to the fact that our observing scheme was successful in avoiding the superior conjunction of the pulsar (around phase 0.25) where regular eclipses normally take place.We can not rule out that eclipses can occur in this range; dedicated observations would be needed to investigate this.Note that this pulsar tends to exhibit very strong scintillation, which is expected given its low DM.We, therefore, conclude that scintillation is the most likely cause for the non-detection in some of the Murriyang Parkes observations, also owing to the fact that the telescope is less sensitive than MeerKAT.

DISCUSSION
In Section 3, we presented the discovery and timing of three new millisecond pulsars in our targeted MeerKAT searches of six highconfidence redback candidates.In this Section, we discuss possible explanations for the lack of detections from the remaining three candidates and discuss the new insights into our discoveries that our pulsation detections and timing solutions provide.

Candidates without pulsation detections
We did not detect any convincing candidate pulsar signals from 4FGL J0523.3−2527,J0940.3−7610 and J1120.0−2204.Possibly, our non-detections are simply due to intrinsic faintness.While these candidates are all found in relatively bright gamma-ray sources (Ballet et al. 2023), pulsar radio fluxes are almost entirely uncorrelated with gamma-ray fluxes (Smith et al. 2023), and so they may yet host very faint radio pulsars.A handful of "radio-quiet" gamma-ray MSPs are also known, thought to be due to emission and viewing geometries being such that their radio beams do not sweep across our line of sight (Clark et al. 2017;Smith et al. 2023).To estimate a nominal flux density upper limit for our observations, we use the pulsar radiometer equation (Lorimer & Kramer 2004), where S/N = 9 is the signal-to-noise ratio required for a confident detection,  sys = 22.5 K is the system temperature,  sky is the sky temperature taken from the 408 MHz all-sky map (Haslam et al. 1982;Remazeilles et al. 2015) scaling to the central frequency of 1284 MHz assuming a spectral index of -2.6,  = 2.8 K Jy −1 is the gain,  pol = 2 is the number of polarisations,  = 700 and 500 MHz is the effective usable bandwidth at -band and UHF,  obs is the observing time,  = 0.15 is the assumed fractional pulse width, and  ≈ 0.7 conservatively accounts for various sources of sensitivity losses in the searching process (digitisation, beamforming efficiency, incoherent harmonic summing, FFT scalloping, acceleration and DM smearing, etc.).For the 15 min, 30 min and 60 min segments that we searched, we find conservative flux density upper limits of around 40 Jy, 30 Jy and 20 Jy (L-band) and 60 Jy, 45 Jy and 30 Jy (UHF), respectively.For an assumed distance of 2.2 kpc (for 4FGL J0523.3−2527, the most distant of the unconfirmed candidates) a putative pulsar would have to have a pseudo-luminosity  1400 =  1400  2 fainter than 95% of known MSPs in the ATNF Pulsar Catalogue (Manchester et al. 2005) to escape detection in our 30 min segments.
Perhaps the most likely explanation for these non-detections is that their pulsations are nearly always rendered undetectable by eclipsing intra-binary material.While we attempted to mitigate this effect by observing these systems when the pulsar was in front of the companion (at orbital phases between 0.5-1.0)according to the published optical ephemerides, transient eclipses are often observed in other redbacks at unexpected orbital phases (e.g.Roy et al. 2014;Deneva et al. 2016), even around the pulsar's inferior conjunction (Kudale et al. 2020), implying that material may sometimes entirely enshroud Red curves and shaded regions show the PSDs of the best-fitting Matérn covariance model, and 1 and 2-sigma uncertainties estimated from the Monte-Carlo samples, respectively.Black curves and grey shaded regions show the measured PSDs and their 95% confidence intervals.The measured PSDs break to higher values than the Matérn models at high frequencies due to measurement uncertainties, and at low frequencies due to uncertainties caused by jointly fitting for the underlying orbital period and phase.Dashed and dot-dashed vertical lines show the minimum and maximum frequencies that were included in the timing model, respectively.Lower plots show the orbital phase variations as a function of time.The greyscale image shows the measured pulsation log-likelihood as a function of offsets from the pulsar's ascending node epoch, measured in overlapping 300-day intervals.Red curves show the 95% confidence interval on the deviations in the pulsar's ascending node time, obtained from the Monte-Carlo samples, with a small offset for clarity.
the binary systems.The recent discovery of PSR J0212+5320 (Perez et al. 2023), which was only detected in two out of five GBT observations targeting an optical redback candidate, is a particularly relevant example for our study.Our MeerKAT timing observations of PSR J0955−3947 (Figure 4) also exhibit short (lasting around 40 s) transient eclipses at an orbital phase close to the eclipse ingress, but still outside of the eclipsing region.
Among the sources with no pulsar identified, flaring activity has been reported in the optical light curve of 4FGL J0523.3−2527, which could be a result of magnetic field reconnection events due to intra-binary material interacting with the shocked pulsar wind (Sironi & Spitkovsky 2011;Halpern et al. 2022).While our repeated observations were intended to increase our chances of detection despite the presence of transient eclipses and interstellar scintillation (e.g.Camilo et al. 2015), at least one spider binary is detectable less than 10% of the time (Ray et al. 2013), and so further observations of these high-confidence targets remain warranted.Indeed, our discovery of PSR J0838−2527 provides another extreme case of redback elusiveness.The discovery observation (Figure 2) shows the signal disappearing due to an apparent eclipse at an orbital phase of 0.7, and following another hour-long detection at UHF in June 2021, this pulsar remained undetectable until August 2023, in spite of multiple monitoring observations.A possible low-significance detection with NRT was made half way through that period, but would have been dismissed without prior knowledge of the timing ephemeris.Observations of transient X-ray flares and optical emission lines from this system suggest that substantial material overflows the companion's Roche lobe (Halpern et al. 2017b), providing a source of eclipsing material.Another tantalising explanation for our non-detections is that this pulsar may be a transitional millisecond pulsar (tMSP) that switches from a radio millisecond pulsar state to a low-mass X-ray binary (LMXB) state, as proposed by Rea et al. (2017) based on X-ray flaring activity.However, the 14-year 4FGL-DR4 light curve shows no evidence of gamma-ray flux variability, which has been observed in the two known Galactic tMSP transitions that have occurred during the Fermi mission (Stappers et al. 2014;Johnson et al. 2015).Hence, while it might not classify as a full transitional MSP, PSR J0838−2527 might be a less extreme version in which the pulsar vanishes for extended periods due to material lying in the system but no accretion disc forms.

Orbital period variations
In Section 3.2, we applied the gamma-ray timing techniques introduced in Clark et al. (2021) to model the redback orbital period variations as stochastic Gaussian processes.Since this procedure has not yet been applied to the full population of Fermi-LAT redbacks, we do not have a lot of context in which to discuss these results.Nevertheless, the results for these three pulsars already reveal some interesting differences with respect to each other, and with PSR J2039−5617, discussed in Clark et al. (2021).The traditional explanation for the orbital period variations commonly seen in spider binaries is the Applegate mechanism (Applegate 1992), in which variations in gravitational quadrupole moment of the companion star, perhaps driven by magnetic activity, couple with the binary orbital period (e.g.Applegate & Shaham 1994;Lazaridis et al. 2011;Pletsch & Clark 2015).Under this model, the change in the companion quadrupole moment () is directly related to the change in the orbital period, via Δ/ ∝ Δ b / b (Voisin et al. 2020), which we can estimate from the derivatives of the Gaussian process described above, with Δ b / b ≲ 10 −6 for all three pulsars.Following the equations in Voisin et al. (2020), and assuming an apsidal motion constant  2 ≳ 10 −3 , we find that Δ/ ≲ 2% for each case.
The spectral indices of the power spectra of the orbital phase variations are also of potential interest, as this parameter encodes information about the mechanism driving the orbital period changes.For PSR J0955−3947, we find that the orbital phase variations are well described by a pure power-law process, with a length scale, ℓ, likely to be longer than our 15-year data set, and with Matérn degree  = 1.4 ± 0.1.This is strikingly close to, and statistically compatible with  = 1.5 (equivalently,  = −4), which is the expected spectral index for a random-walk process in the orbital period.This therefore implies that the star's quadrupole moment undergoes a random walk.In contrast, for PSR J2333−5526, a smoothly broken power law is preferred, with characteristic length scale 300 d < ℓ < 800 d, well within the frequencies probed by our data, breaking to a steep powerlaw process (with  > 5) at higher frequencies.The steeper index implies that these changes are more gradual than in PSR J0955−3947, and that there are time-varying torques (e.g.mass loss, asynchronous rotation) acting on the companion star.
In this respect, PSR J2333−5526 is much more similar to PSR J2039−5617, which has a similar characteristic length scale, and a steeply-decreasing spectrum at higher frequencies (Clark et al. 2021).The picture is less clear for PSR J0838−2827.Its lower spin frequency and smaller semi-major axis compared to the other two pulsars increase the uncertainty on the orbital phase, leading to much greater uncertainty in the hyperparameters.Here the power spectrum could either be a pure, shallow power-law like J0955−3947, or could break to a steeper power low at high frequencies, like J2333−5526.
The orbital period variations seen in PSRs J0955−3947 and J2333−5526 have a larger amplitude, and are more complicated, than those that have previously been subjected to gamma-ray timing studies (PSR J2339−0533 in Pletsch & Clark 2015 and PSR J2039−5617 in Clark et al. 2021).These pulsars required substantial improvements in our methodology to reach a satisfactory timing solution.Aside from allowing the orbital period variations in redbacks to be studied more quantitatively, our new methods also provide a possible path for bright gamma-ray redbacks to contribute to the Gammaray Pulsar Timing Array (PTA, FERMI-LAT Collaboration et al. 2022), several of which have been excluded until now due to their complicated behaviour.Both PSRs J0955−3947 and J2333−5526 have log-likelihoods above the threshold for inclusion in the PTA.
For PSR J0838−2827, we required four spin frequency derivatives to model the phase evolution over the Fermi-LAT mission.These terms are unusual for a gamma-ray MSP, where long-term red noise is all but unmeasurable in gamma-ray data for most MSPs (FERMI- LAT Collaboration et al. 2022).In one other black-widow binary, PSR J1555−2908, the presence of significant higher order spin frequency variations has been interpreted as possible evidence for a third, planetary mass object orbiting the inner binary system (Nieder et al. 2022).Further timing analyses of this pulsar as the Fermi-LAT mission continues may determine whether or not this timing noise is intrinsic to the pulsar.

Distances and energetics
The detection of radio pulsations from these MSPs provides knowledge of their DMs, which can be used to infer distances based on models of the Galactic electron density content.These can be compared to distances estimates obtained from Gaia parallax measurements (or upper limits), and constraints from previous optical modelling.
For PSRs J0838−2827 and J0955−3947 there is significant tension between the DM-inferred distances ( = 410 pc and  = 490 pc, respectively, according to YMW16) and the lower bounds from Gaia DR3 parallax ( = 0.43±0.40mas and  = 0.27±0.13mas, respectively, where the 2 upper limits correspond to minimum distances of  = 800 pc and  = 1.9 kpc, respectively).Upon closer inspection, these pulsars lie close to the edge of the Gum Nebula, perhaps indicating that the shell-shaped excess electron density added to the YMW16 model to account for this feature may be overestimated in these directions.The NE2001 model predicts larger distances for these pulsars,  = 560 pc and  = 3.3 kpc, respectively, compatible with the parallax range for PSR J0955−3947 but still not with that for J0838−2827.For PSR J2333−5526, the Gaia DR3 parallax ( = −0.37 ± 0.52, for  > 1.5 kpc) all but rules out the NE2001 distance ( = 860 pc), but the YMW16 distance  = 2.5 kpc is plausible.For further calculations, we adopt nominal distances of:  = 1.7 kpc for PSR J0838−2827 (the maximum distance found by Halpern et al. 2017a to be compatible with the observed optical flux);  = 3.5 kpc for PSR J0955−3947 (the distance obtained by Koljonen & Linares 2023 using the Gaia DR3 parallax); and  = 3.1 kpc for PSR J2333−5526 (from optical modelling by Swihart et al. 2020).
The gamma-ray timing solutions, constrained by Gaia DR3 astrometry, also provide precise values for proper motion that can be used to Doppler-correct the observed spin-down rates, and therefore the spin-down energy budget, of these pulsars for the radial acceleration induced by their transverse motion, known as the "Shklosvkii effect" (Shklovskii 1970), and for their acceleration in the Galactic potential.For PSRs J0955−3947 and J2333−5526 these corrections are small, totalling less than 10% of the observed spin-down rate, but for PSR J0838−2827 the Shklovskii component is larger, accounting for more than 20% of the observed  at  = 1.7 kpc.These spindown corrections are included in the derived parameter values given in Table 2.

Mass measurement
Our timing solutions provide precise values for the pulsars' projected semi-major axes,  1 , which were the missing pieces for pulsar mass measurements for these systems.
The pulsar mass,  psr , can be estimated from these parameters via the binary mass function, where  b is the orbital period,  1 = 2  1 / b is the pulsar's projected radial velocity,  2 is the companion radial velocity amplitude and  is the gravitational constant.The companion radial velocities, reported in Table 2, were obtained via optical spectroscopy taken when these binary systems were first identified as candidate redbacks5 .We estimated pulsar masses using a Monte-Carlo technique, in which we generate 50,000 random realisations, with values for  2 and  1 drawn from normal distributions according to our timing solutions, and with cos  drawn from a uniform distribution (representing an isotropic distribution for the binary orbital axis), discarding values lying outside the constraints from the gamma-ray eclipse observations (e.g.requiring  > 78 • for PSR J0838−2827).The results of these calculations are presented in Table 2.These values span a large range, from  psr = 0.93 ± 0.14M ⊙ for PSR J0838−2527, to  psr = 2.06 ± 0.13M ⊙ for PSR J2333−5526.
However, one significant source of systematic uncertainty in these values remains, which is the "centre-of-light" adjustment to  2 , required to correct the observed radial velocity curves to account for the effects of the large difference in surface temperature between the irradiated and un-irradiated faces of the companion star.These corrections are difficult to estimate, as they vary between different spectral lines, and depend strongly on the temperature pattern on the surface of the companion star.Swihart et al. (2020) account for this in an approximate manner for PSR J2333−5526 by producing a model of the companion star surface temperatures, and using an effective temperatureequivalent width relation for the adopted Mg b triplet to produce a "corrected" model radial velocity curve, which reduces the inferred  2 .Adopting this value, we find a lower  psr = 1.70 ± 0.06M ⊙ .However, this is at odds with the picture found by Kandel & Romani (2020) for the similarly-irradiated PSR J2215+5135, in which the Mg b lines only slightly under-estimate the centre-of-mass velocity, suggesting that the original uncorrected value may be closer to reality.Dodge et al. (2024) find almost no centre-of-light correction is required to the velocity of the Mg b line in a similarly irradiated redback, PSR J1910−5320, with an uncertainty of around 4% on the correction, suggesting a systematic uncertainty here of around 0.2-0.3M ⊙ for PSR J2333−5526.
Estimates for the centre-of-light corrections have not been made for either PSR J0838−2527 or J0955−3947.For both pulsars, radial velocities were predominantly measured from the Ca II triplet (Halpern et al. 2017a;Li et al. 2018), whose equivalent width varies little with temperature (Mallik 1997) and which will therefore track a position close to the heated face of the companion star which contributes most strongly to the observed spectra.The true centre-of-mass radial velocity amplitude will therefore be larger than observed, by up to around 20% (e.g.Romani & Shaw 2011), meaning our pulsar mass estimates are better interpreted as lower bounds for these pulsars, and the true mass could be up to 60% larger.
A full modelling treatment of the optical light curves and spectra (or radial velocities) (e.g.Linares et al. 2018;Kandel & Romani 2020;Kennedy et al. 2022) will be necessary to turn our lower limits into more precise pulsar mass estimates for these systems.For PSR J0838−2527 in particular, such an effort will likely be complicated by the strong asymmetry and variability observed by Cho et al. (2018).Nevertheless, our results are consistent with the picture from Strader et al. (2019) that redbacks tend to contain rather heavy neutron stars, and in particular the mass of PSR J2333−5526 is likely to lie towards the heavier end of this distribution.

CONCLUSION AND FUTURE WORK
We have presented the results from TRAPUM's deep survey for six redback candidates using the MeerKAT radio telescope.We found pulsations from MSPs in three of these systems, thus confirming their nature as redbacks.We detected PSRs J0955−3947 and J2333−5526 in most of our search and timing observations.The data we collected enabled us to establish phase-connected timing solutions, which enabled the detection of gamma-ray pulsations and a full 15-year gamma-ray timing solution accounting for significant stochastic variability in the orbital periods, consistent with gravitational quadrupole moment variations in the companion star.
On the contrary, PSR J0838-2827 remained undetected for two years after our discovery, except for a single faint detection with NRT, until August 2023 when we started to re-detect this source with MeerKAT and Parkes.With these new detections, we could find a local phase-connected timing solution, which enabled gammaray timing and finding gamma-ray pulsation.We should note that PSR J0838−2827 has been reported as a possible transitional millisecond pulsar candidate (Rea et al. 2017), but gamma-ray pulsations are detectable throughout 15 years of monitoring with the Fermi-LAT, with no apparent flux variations or transitions.
We encourage further monitoring of PSR J0838−2827, as its elusive behaviour, flaring X-ray activity and optical emission features might all point towards future transitional behaviour between pulsar and X-ray binary-like states.Likewise, we suggest that our inability to detect pulsations in the other three candidates is likely caused by a combination of heavy eclipses and potentially faint pulsations.Further observations might eventually unravel these pulsars.

Figure 1 .
Figure1.Pulse profiles of the three newly discovered redback pulsars at L-band and UHF.The horizontal axis displays the spin phase of each source, replicated for clarity, while the vertical axis depicts the intensity, normalised to the peak value.The grey line presents the background level, calculated as the median of the off-pulse region.Since the gamma-ray ephemerides do not perfectly fold radio data to perform an absolute alignment, we aligned the L-band (black) and the UHF (red) profiles using a cross-correlation.

Figure 2 .
Figure2.Discovery observation of PSR J0838−2827, with the folded pulse profile summed over frequency and time at the top, and the folded pulse profiles summed in frequency, but with sub-integrations running on the vertical axis at the bottom.The pulse cycle is replicated twice for clarity.In this L-band observation dedispersed at   = 80.55 pc cm −3 , the pulsations are clearly detected at the start before fading into a total eclipse around 1000 s after the start.The pulse phase is also seen shifting due to the additional electron column density from the eclipsing medium increasing the DM temporarily.

Figure 3 .
Figure3.Archival Parkes data of PSR J0838−2827 folded using the timing parameters from MeerKAT (panels are as explained in Figure2).A weak pulse can be seen later in the observation starting approximately 1500 s into the observation after coming out of an eclipse.

Figure 4 .
Figure 4.A five-minute MeerKAT L-band timing observation of PSR J0955−3947, beginning at orbital phase 0.05 (panels are as explained in Figure 2).Short eclipses can be seen throughout the observation around 100-140 s and 220-260 s.

Figure 5 .
Figure 5.The orbital phase coverage of the MeerKAT L-band, UHF, Parkes (PKS), GBT and NRT observation.Phase 0 corresponds to the pulsar at ascending node, and phase 0.25 is superior conjunction.The colours represent observations in the following bands/telescopes: UHF (red), L-band (blue), PKS (green), GBT (teal), NRT (purple), and L-band with short eclipses observed (yellow).Detections are displayed by brighter colors with dots, while non-detection observations use dimmer colour without dots.The grey region is when the eclipse clearly shows up in the observation.Eclipses in L-band and PKS archival data are shown in Figure 2 and 3.An example of short eclipses is in Figure 4.

Figure 6 .
Figure 6.Orbital phase variations over the course of the Fermi-LAT data.Upper panels show the power spectral densities (PSD) of the orbital phase variations.Red curves and shaded regions show the PSDs of the best-fitting Matérn covariance model, and 1 and 2-sigma uncertainties estimated from the Monte-Carlo samples, respectively.Black curves and grey shaded regions show the measured PSDs and their 95% confidence intervals.The measured PSDs break to higher values than the Matérn models at high frequencies due to measurement uncertainties, and at low frequencies due to uncertainties caused by jointly fitting for the underlying orbital period and phase.Dashed and dot-dashed vertical lines show the minimum and maximum frequencies that were included in the timing model, respectively.Lower plots show the orbital phase variations as a function of time.The greyscale image shows the measured pulsation log-likelihood as a function of offsets from the pulsar's ascending node epoch, measured in overlapping 300-day intervals.Red curves show the 95% confidence interval on the deviations in the pulsar's ascending node time, obtained from the Monte-Carlo samples, with a small offset for clarity.

Figure 7 .
Figure 7. Gamma-ray pulsations from the three new redback MSPs.Lower panels show the weighted gamma-ray photon phases according to the best-fitting timing solution, while upper panels show the integrated pulse profile.The best-fitting template pulse profile is shown by an orange curve, with faint black curves underneath showing templates randomly drawn from the Monte-Carlo samples to illustrate the uncertainty on the template pulse profile.

Table A1 -
continued List of observations.