# Neutron star mass estimates from gamma-ray eclipses in spider millisecond pulsar binaries

C.J. Clark, M. Kerr, E.D. Barr, B. Bhattacharyya, R.P. Breton, P. Bruel, F.

Camilo, W. Chen, Ismaël Cognard, H.T. Cromartie, et al.

► To cite this version:

C.J. Clark, M. Kerr, E.D. Barr, B. Bhattacharyya, R.P. Breton, et al. Neutron star mass estimates from gamma-ray eclipses in spider millisecond pulsar binaries. Nature Astron., 2023, 7 (4), pp.451-462.
10.1038/s41550-022-01874-x. hal-03960032

HAL Id: hal-03960032

https://hal.science/hal-03960032

Submitted on 28 Jan 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Neutron star mass estimates from gamma-ray eclipses in spider millisecond pulsar binaries

Received: 21 January 2022

Accepted: 1 December 2022

Published online: 26 January 2023

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

Check for updates

C. J. Clark Ⓔ¹,²,³, M. Kerr Ⓔ⁴, E. D. Barr⁵, B. Bhattacharyya⁶, R. P. Breton Ⓔ³, P. Bruel Ⓔ⁷, F. Camilo⁸, W. Chen Ⓔ⁵, I. Cognard⁹,¹⁰, H. T. Cromartie¹¹, J. Deneva⁴,¹², V. S. Dhillon Ⓔ¹³,¹⁴, L. Guillemot⁹,¹⁰, M. R. Kennedy Ⓔ³,¹⁵, M. Kramer Ⓔ³,⁵, A. G. Lyne³, D. Mata Sánchez³,¹⁴,¹⁶, L. Nieder Ⓔ¹,², C. Phillips¹⁷, S. M. Ransom Ⓔ¹⁸, P. S. Ray Ⓔ⁴, M. S. E. Roberts Ⓔ¹⁹, J. Roy⁶, D. A. Smith Ⓔ²⁰, R. Spiewak Ⓔ³,²¹,²², B. W. Stappers³, S. Tabassum Ⓔ²³,²⁴, G. Theureau⁹,¹⁰,²⁵ & G. Voisin Ⓔ³,²⁵

Reliable neutron star mass measurements are key to determining the equation of state of cold nuclear matter, but such measurements are rare. Black widows and redbacks are compact binaries consisting of millisecond pulsars and semi-degenerate companion stars. Spectroscopy of the optically bright companions can determine their radial velocities, providing inclination-dependent pulsar mass estimates. Although inclinations can be inferred from subtle features in optical light curves, such estimates may be systematically biased due to incomplete heating models and poorly understood variability. Using data from the Fermi Large Area Telescope, we have searched for gamma-ray eclipses from 49 spider systems, discovering significant eclipses in 7 systems, including the prototypical black widow PSR B1957+20. Gamma-ray eclipses require direct occultation of the pulsar by the companion, and so the detection, or significant exclusion, of a gamma-ray eclipse strictly limits the binary inclination angle, providing new robust, model-independent pulsar mass constraints. For PSR B1957+20, the eclipse implies a much lighter pulsar (1.81 ± 0.07 solar masses) than inferred from optical light curve modelling.

Since the discovery of the first black widow pulsar, B1957+20, in 1988¹, a sizable population of compact binary millisecond pulsar systems with semi-degenerate companion stars has emerged². These are often split into two main classes: black widows with companion stars with masses below $M_c \lesssim 0.05$ solar masses ($M_\odot$); and redbacks with companion masses $0.1 M_\odot \lesssim M_c \lesssim 0.5 M_\odot$. The characteristic signatures of black widow and redback systems are periodic disappearances of radio pulsations, often lasting for a large fraction of the orbital period. Despite being referred to as 'eclipses', these events are too long to be caused by occultations by the companion star, but are instead explained by dispersion, scattering and absorption of radio emission by diffuse

intra-binary material³,⁴. This material is thought to have been ablated from the companion star's outer envelope by the intense pulsar wind. The spider nicknames come from this destructive behaviour, by analogy with arachnid species whose females have a (perhaps unfair) reputation for killing their lighter mates.

A key motivation for finding and studying new spider pulsars is that they are one of the few types of pulsar binary system from which neutron star mass estimates can be obtained. This is because these systems have optically bright companion stars, whose radial velocities can be measured via optical spectroscopy. Dividing these by the pulsar radial velocities (measured by pulsar timing) provides binary
mass ratio measurements, with which inclination-dependent mass estimates can be obtained by solving the binary mass function (see 'Pulsar mass constraints' in the Methods). Large neutron star masses have been inferred in this way from individual spider systems (for example refs. ⁵⁻⁷) and there are hints that spider pulsars may be systematically heavier than other species of binary neutron star⁸,⁹. Several classes of theoretical neutron star equation-of-state (see ref. ¹⁰ and references therein) models predict maximum masses close to $2 M_\odot$, and so precise measurements of neutron star masses close to or above this level could have substantial implications for fundamental nuclear physics.

However, spider pulsar mass estimates via radial velocity measurements depend strongly on the estimated binary inclination angle, *i*, with the inferred mass $M_{psr} \propto 1/\sin^3 i$. The ability to accurately measure inclination angles for spider systems is therefore crucial if their masses are to be used to probe the nuclear equations of state.

Binary inclination angles in spider systems are commonly estimated by modelling their optical light curves, which exhibit inclination-dependent features due to tidal deformation of the companion star in the pulsar's gravitational field, and heating by the pulsar. However, these models are sensitive to the exact temperature pattern on the companion star's surface, which often deviates significantly from that predicted by simple models in which the pulsar wind directly heats the inner face of the companion star, for example due to heating contributions from an intra-binary shock between the pulsar and stellar winds¹¹. In several redbacks, this temperature pattern is even seen to vary over time (for example refs. ¹²⁻¹⁴). Owing to the $\sin^3 i$ scaling, a systematic error in the estimated inclination angle due to an incomplete heating model can lead to a large bias in the resulting pulsar mass estimate.

Millisecond pulsars (MSPs) also emit gamma-ray pulsations, as revealed by the Large Area Telescope (LAT)¹⁵ onboard the Fermi Gamma-ray Space Telescope. Gamma rays are particularly helpful in discovering and studying spider pulsars because, unlike radio waves, they are not absorbed in the diffuse intra-binary material. Partially as a result of radio searches repeatedly targeting unidentified Fermi-LAT sources¹⁶, the number of known Galactic spider systems has increased tenfold since Fermi's launch.

The LAT data also offer a new and independent means to constrain binary inclination angles and pulsar masses in spider systems by enabling searches for, and studies of, gamma-ray eclipses. Strader et al.¹⁷ found evidence in the Fermi-LAT data, later confirmed by Kennedy et al.¹⁸ using a longer dataset, for the first such gamma-ray eclipse in a candidate accreting transitional MSP.

In this Article we present a systematic search for gamma-ray eclipses from 49 confirmed/candidate spider pulsar systems. We first searched for eclipses in 42 Fermi-LAT detected confirmed spider pulsars, and found significant gamma-ray eclipses from five pulsars: PSRs B1957+20, J1048+2339, J1555–2908, J1816+4510 and J2129–0429. This number of detected eclipses is consistent with the number we would expect to observe from the tested population, assuming randomly distributed orbital axes and nearly Roche-lobe filling companion stars. The orbital gamma-ray light curves for these systems are shown in Fig. 1, and the results of Monte Carlo simulations used to estimate their significance (see 'Significance calibration via Monte Carlo simulations' in the Methods) are shown in Fig. 2. Of these, the eclipse in PSR J1555–2908 has the lowest significance, but still has a false-alarm probability of $4 \times 10^{-5}$ after accounting for the trials factor introduced by testing a range of possible eclipse widths. The most significant eclipse, from PSR J2129–0429, represents a deficit of no more than 20 (weighted) photons over the 11.4 yr of LAT data considered here.

For 32 of the pulsars without detected eclipses the gamma-ray data significantly exclude otherwise feasible eclipses above a certain duration, but the faintest five systems yield no such constraints. In one of these, PSR J0251+2606, there is marginal evidence for an eclipse, with a false-alarm probability of 0.002 (see 'Significance calibration via

Monte Carlo simulations' in the Methods). Given the number of pulsars included in our search, this is around a factor of ten lower than expected for the largest outlier ($1/n_{psr} \approx 0.02$). If this is indeed an eclipsing system, another 7 yr of accumulated Fermi-LAT data will be required to reach the same significance as the eclipse in PSR J1555–2908.

We also searched for gamma-ray eclipses from seven likely redback systems that were first identified from the discovery of periodic optical and/or X-ray sources within pulsar-like Fermi-LAT sources, but either currently lack radio or gamma-ray pulsation detections to confirm their nature or have only recently been detected as pulsars. These systems do not yet have precise orbital ephemerides from pulsar timing, and so a search over orbital period and phase is required, which introduces a large trials factor and therefore greatly reduces sensitivity. Nevertheless, from two of these systems, PSRs J0838–2827 and J2333–5526, we found evidence for eclipses with trials-corrected false-alarm probabilities below $3 \times 10^{-3}$, which are also shown in Fig. 1. The search results for these two systems are shown in Fig. 3.

We consider it unlikely that the gamma-ray eclipses are caused by the same mechanism as radio eclipses—that is, absorption by diffuse material evaporated from the companion star. At LAT photon energies, the primary interaction between gamma rays and matter is through pair production, which for hydrogen gas has cross-section $\sigma_v = 0.03\sigma_T$ (ref. ¹⁹), where $\sigma_T$ is the Thomson cross-section. We can estimate the electron column depth in the radio eclipse region using measurements of the excess radio dispersion measure (DM). In terms of the radio DM excess, the gamma-ray optical depth is $\tau_v = \sigma_v r^{-1} \Delta DM \approx 6 \times 10^{-8} r^{-1} (\Delta M/1 \text{ pc cm}^{-3})$, where *r* is the ionization fraction.

Polzin et al.⁴ have studied the radio eclipses in PSRs B1957+20 and J1816+4510. The larger DM excess was seen in J1816+4510 with a value of $\Delta DM \approx 0.1 \text{ pc cm}^{-3}$ at a phase 0.025 orbits after conjunction, corresponding to $\tau_v \approx 6 \times 10^{-9} r^{-1}$. Assuming an isotropic wind, with density $\rho$ decreasing with radius *r* from the companion star as $\rho(r) \propto r^{-2}$, the optical depth may be around 100 times higher for orbital phases half way between conjunction and our measured eclipse egress. This is still several orders of magnitude too low to explain the observed eclipses unless the ionization fraction is extremely low, which seems highly unlikely given the intense environment. Similar values were obtained for B1957+20 and J1048+2339 (ref. ²⁰).

Of course, this model for the companion wind is overly simplistic: spider companions have non-isotropic swept-back winds¹¹ that often vary with time¹³, and the wind density profile and degree of ionization are not yet known. If the stellar wind contributed meaningfully to the gamma-ray optical depth, we could hope to see signs of this in the Fermi-LAT data, for example from gradual ingresses/egresses due to tenuous intra-binary material that increases in density towards the companion star or variability with time in the eclipse properties. Unfortunately, the data are not sensitive to these effects due to the very low number of expected missing photons within the relevant orbital phases. Nevertheless, for all detected eclipses, the inferred fluxes within the eclipse regions are consistent with zero, and sudden (rather than gradual) ingresses or egresses are statistically preferred. The eclipse durations and depths also do not seem to vary over time, at least on the long timescales that we are sensitive to: eclipse widths measured from the first and second halves of the Fermi-LAT data are consistent within their 1σ uncertainties, and the eclipse log L increase approximately linearly with accumulated exposure. The observed eclipse durations are also consistent when measured in different energy bands (above or below 1 GeV).

The observed eclipses are also short enough to be caused by companion stars that fill some or all of their Roche lobes (Table 1). The longest eclipse is observed in PSR J1048+2339, lasting for 6–12% of the orbital period, while the maximum eclipse duration expected for a Roche-lobe filling companion in this system is 8%. Interestingly, optical spectroscopy has revealed emission from matter close to the L1 Lagrange point in this system, and emission lines are seen in the
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Gamma-ray orbital light curves of seven eclipsing spider pulsars. The red dashed lines show the estimated background level. Phase 0 corresponds to the pulsar's ascending node. The phase of the pulsar's superior conjunction, where eclipses would be expected to occur, has been placed at the centre of a</figcaption></figure>

phase bin, and is shown at the centre of the plots for emphasis. Bin widths were chosen to be close to the best-fitting eclipse duration. Bin heights show the sum of the photon weights in each orbital phase bin, and error bars show the corresponding $1\sigma$ Poisson uncertainties.

<figure><img src="image_5.png" /><figcaption>**Fig. 2** | Results of Monte Carlo simulations used to calibrate eclipse false-alarm probabilities. Vertical lines show the measured log likelihood (log L) values, maximized over eclipse widths, for each pulsar. Those for pulsars with significant eclipses are marked by coloured lines. The coloured curves show the false-alarm probabilities from simulations using the distributions of photon weights from each of the five eclipsing pulsars. Horizontal dashed lines show the corresponding false-alarm probabilities according to the Monte Carlo calibration. The dotted and solid black curves show the empirical survival function (that is, the fraction of pulsars that survive a given log L threshold) for the real population of spiders studied here, before and after removing the five pulsars with significant eclipses, respectively. The diagonal dashed line is an extrapolation of the fit to the simulated false-alarm probability curves used to estimate the false-alarm probabilities for the most significant eclipses.</figcaption></figure>

spectra of several other redbacks, including PSR J0838−2527 (ref. 21) in which we detected an eclipse with a shorter duration (≤2% of an orbit; Supplementary Table 1). This suggests that some degree of overflowing material may be common in redback systems, but the observed eclipse durations do not provide evidence for gamma-ray absorption from this material at present.

Our observations are therefore all consistent with eclipses that are solely due to occultations of the pulsar by the companion star. Under this simpler assumption, the detection of a gamma-ray eclipse and the measurement of its duration, or the significant non-detection of an eclipse, provide a robust constraint on the binary inclination. For spider systems whose companion radial velocity curves have been measured through optical spectroscopy, these inclination limits in turn constrain the pulsar masses (see 'Pulsar mass constraints' in the Methods). For eclipsing spider systems, the minimum eclipse duration provides a lower limit for the inclination, and hence an upper limit on the pulsar mass. By the same logic, we can obtain upper limits on the inclination and lower limits on the pulsar mass for systems that are not eclipsing. We list these pulsar mass constraints for eclipsing and non-eclipsing systems in Tables 1 and 2, respectively, and illustrate these results in Fig. 4. One of the eclipsing pulsars, PSR J1816+4510, has a mass upper limit larger than $2 M_\odot$. Of the non-eclipsing pulsars, the extremely compact black widow binary PSR J1653−0158 (ref. 22) has the largest minimum mass at $1.76 M_\odot$.

The resulting inclination limits also provide crucial independent tests that can validate or refute multiwavelength models, including optical (for example refs. 23,24) and X-ray light curve models (for example ref. 25), and radio and gamma-ray pulse profile models (for example ref. 26), all of which have *i* as a free parameter.

For four of the five eclipsing pulsars, our inferred inclination constraints are consistent with existing optical modelling results (or no inclination constraints from optical modelling exist in the literature, see 'Optical constraints for eclipsing pulsars' in the Methods), but this is not the case for PSR B1957+20. Modelling of photometric observations of B1957+20 yields inclination estimates of 63° $\lesssim i \lesssim 67°$ (refs. 23,24). When combined with optical spectroscopy results³, this corresponds to an extremely high mass of $M_{\text{psr}} = 2.4 \pm 0.1 M_\odot$, higher than that of any other known neutron star. This is at odds with most equation-of-state models, which predict lower maximum neutron star masses¹⁰. Our detection of a gamma-ray eclipse, however, requires a much higher $i$, >84.1°. This lower bound on the inclination corresponds to $M_{\text{psr}} = 1.81 \pm 0.07 M_\odot$, with the uncertainty now dominated by the radial velocity measurement and centre-of-mass correction. This mass is more consistent with the most massive neutron stars found by more robust pulsar timing studies (such as ref. 27 and references therein). The nearly edge-on inclination is also more consistent with the estimates by
<figure><img src="image_6.png" /></figure>

**Fig. 3 |** Search results for the two candidate redback systems in which eclipses are detected. Left: eclipse log *L* as a function of the orbital parameters, maximized over eclipse durations. The blue crosshairs denote the 2σ ranges around the orbital ephemeris from optical observations (see references in the Supplementary Information). Contour lines are drawn at log *L* corresponding to false-alarm probabilities of 5% (yellow) and 1% (red). These levels are also marked on the colour bar, along with the log *L* corresponding to a false-alarm probability of 0.1% (dark red), although this level is never reached. The position of the

maximum likelihood ($L_{max}$) is marked by a cyan cross, and the corresponding log *L* is also marked in cyan on the colour bar. Right: results of the Monte Carlo simulations used to calibrate these false-alarm probabilities. The 5%, 1% and 0.1% levels are marked in the same colours used in the left panels, with the log $L_{max}$ found in the search and corresponding false-alarm probability marked in cyan. Top: PSR J0838−2827. Bottom: PSR J2333−5526. In all other candidate systems, the false-alarm probabilities for the log $L_{max}$ found were greater than 40%.

<table><caption>**Table 1 |** Constraints for pulsars with detected eclipses</caption><thead><tr><th>Pulsar</th><th>Class</th><th>δ log L</th><th>P<sub>FA</sub></th><th>θ<sup>min</sup></th><th>θ<sup>max</sup></th><th>K<sub>c</sub> (kms<sup>-1</sup>)</th><th>i<sup>min</sup>(°)</th><th>M<sub>psr</sub> (M<sub>☉</sub>)</th><th>M<sub>c</sub> (M<sub>☉</sub>)</th><th>Reference</th></tr></thead><tbody><tr><td>B1957+20</td><td>Black widow</td><td>12.63</td><td>2×10<sup>-6</sup></td><td>0.007</td><td>0.011</td><td>353.0±4.0</td><td>84.1</td><td>1.67–1.94</td><td>0.025–0.027</td><td>5</td></tr><tr><td>J1048+2339</td><td>Redback</td><td>13.28</td><td>2×10<sup>-6</sup></td><td>0.058</td><td>0.120</td><td>343.3±4.4</td><td>80.4</td><td>1.44–1.72</td><td>0.31–0.35</td><td>51</td></tr><tr><td>J1555–2908</td><td>Black widow</td><td>10.07</td><td>4×10<sup>-5</sup></td><td>0.023</td><td>0.040</td><td>397.0±2.0</td><td>83.1</td><td>1.58–1.71</td><td>0.057–0.060</td><td>30</td></tr><tr><td>J1816+4510</td><td>Redback</td><td>17.32</td><td>&lt;1×10<sup>-7</sup></td><td>0.014</td><td>0.019</td><td>343.0±7.0</td><td>79.0</td><td>1.64–2.17</td><td>0.18–0.22</td><td>47</td></tr><tr><td>J2129–0429</td><td>Redback</td><td>17.67</td><td>&lt;1×10<sup>-7</sup></td><td>0.030</td><td>0.036</td><td>250.3±4.3</td><td>76.3</td><td>1.48–1.93</td><td>0.39–0.47</td><td>54</td></tr></tbody></table>

$θ^{min}$ and $θ^{max}$ are the minimum and maximum eclipse durations (in orbits) at 95% confidence. $i^{min}$ is the limiting inclination at which the minimum eclipse duration would be reached for a Roche-lobe filling companion, assuming the 2σ lower limit on the companion radial velocity amplitude $K_c$, which is taken from the references listed in the final column, with 1σ uncertainties quoted. $M_{psr}$ and $M_c$ give the (conservative) range of pulsar and companion masses that are allowed by the eclipse detection. The minimum masses were found by assuming $i=90°$ and the 2σ lower limit on $K_c$, and require companions to substantially underfill their Roche lobes (with the exception of PSR J1048+2339). The maximum masses were found by assuming the 2σ upper limit on $K_c$, and the minimum inclination required to produce the minimum eclipse duration with the conservative assumption of a Roche-lobe filling companion.

Guillemot et al.²⁸ and Johnson et al.²⁶ obtained from modelling the radio and gamma-ray pulse profiles, implying that their assumption that the pulsar’s spin becomes aligned with the orbit during recycling is correct.

How then do we interpret the light curve models of Reynolds et al.²³ and Draghis et al.²⁴, which consistently estimate far lower inclinations than we find here? Compared with models with intermediate inclinations, a model with nearly edge-on inclination will predict less
<table><caption>Table 2 | Constraints for pulsars without detected eclipses</caption><thead><tr><th>Pulsar</th><th>Class</th><th>δ log L</th><th>θ<sup>max</sup></th><th>K<sub>c</sub> (kms<sup>-1</sup>)</th><th>r<sup>max</sup>(s)</th><th>M<sup>min</sup><sub>psr</sub> (M<sub>☉</sub>)</th><th>M<sup>min</sup><sub>c</sub> (M<sub>☉</sub>)</th><th>Reference</th></tr></thead><tbody><tr><td>J0952-0607</td><td>Black widow</td><td>0.36</td><td>0.010</td><td>376.1±5.1</td><td>86.3</td><td>1.40</td><td>0.020</td><td>72</td></tr><tr><td>J1023+0038</td><td>Redback</td><td>0.00</td><td>0.007</td><td>295.0±3.0</td><td>81.9</td><td>0.65</td><td>0.09</td><td>73</td></tr><tr><td>J1227-4853</td><td>Redback</td><td>0.00</td><td>0.002</td><td>294.4±4.0</td><td>81.1</td><td>1.01</td><td>0.18</td><td>73</td></tr><tr><td>J1301+0833</td><td>Black widow</td><td>0.31</td><td>0.013</td><td>284.9±4.4</td><td>85.7</td><td>0.63</td><td>0.014</td><td>74</td></tr><tr><td>J1311-3430</td><td>Black widow</td><td>0.00</td><td>0.001</td><td>633.9±5.3</td><td>86.9</td><td>1.66</td><td>0.009</td><td>6</td></tr><tr><td>J1431-4715</td><td>Redback</td><td>2.54</td><td>0.066</td><td>278.0±3.0</td><td>90</td><td>1.13</td><td>0.11</td><td>8</td></tr><tr><td>J1622-0315</td><td>Redback</td><td>0.05</td><td>0.009</td><td>423.0±8.0</td><td>83.4</td><td>1.33</td><td>0.10</td><td>8</td></tr><tr><td>J1628-3205</td><td>Redback</td><td>2.72</td><td>0.009</td><td>358.0±10.0</td><td>82.2</td><td>1.09</td><td>0.14</td><td>8</td></tr><tr><td>J1653-0158</td><td>Black widow</td><td>0.74</td><td>0.000</td><td>700.2±7.9</td><td>86.7</td><td>1.76</td><td>0.012</td><td>22</td></tr><tr><td>J1810+1744</td><td>Black widow</td><td>0.00</td><td>0.001</td><td>462.3±2.2</td><td>84.7</td><td>1.59</td><td>0.049</td><td>75</td></tr><tr><td>J2039-5617</td><td>Redback</td><td>0.40</td><td>0.001</td><td>327.2±5.0</td><td>81.7</td><td>1.02</td><td>0.14</td><td>8,14</td></tr><tr><td>J2215+5135</td><td>Redback</td><td>0.00</td><td>0.001</td><td>412.3±5.0</td><td>81.6</td><td>1.58</td><td>0.23</td><td>7</td></tr><tr><td>J2339-0533</td><td>Redback</td><td>0.00</td><td>0.000</td><td>377.6±17.7</td><td>81.1</td><td>1.21</td><td>0.24</td><td>76</td></tr></tbody></table>

Parameters are the same as in Table 1, but inclination upper limits *r*<sup>max</sup> and mass lower limits *M*<sup>min</sup><sub>psr</sub> and *M*<sup>min</sup><sub>c</sub> were computed for the maximum eclipse duration, and assuming that the companion is 50% Roche-lobe filling, to obtain lower limits on the component masses.

<figure><img src="image_7.png" /><figcaption>**Fig. 4 | Neutron star mass constraints for gamma-ray detected spider MSPs using constraints obtained from the detection or exclusion of gamma-ray eclipses.** The five pulsars with detected eclipses are highlighted in bold. The two additional eclipsing systems, PSRs J0838–2827 and J2333–5526, are excluded from this plot as their mass ratios are not yet known from pulsar timing and so their masses cannot yet be estimated. The colour indicates the subclass of spider system: black widows are shown in black, redbacks are shown in red. For pulsars with no detected eclipses we show lower limits on the pulsar mass indicated with arrows of arbitrary length.</figcaption></figure>

flux at minimum (for the same stellar temperature model) as less of the heated face is visible when the companion is in front of the pulsar. The irradiation pattern must therefore extend further around the back side of the companion star than predicted by a direct-heating model to increase the minimum flux to match the observed photometry at these orbital phases. Such excess heat could be caused by redirection of heating flux by an intra-binary shock wrapping around the companion star<sup>11</sup> or diffusion on the stellar surface causing heat to ‘leak’ over the terminator<sup>29</sup>—a possibility that Reynolds et al.<sup>23</sup> noted in their

original modelling of this system. Revision of the optical modelling for this pulsar, using extended heating models such as these, will be required to resolve the tension with the inclination range inferred from our eclipse detection. All but one of the other eclipsing systems are redbacks, whose companions tend to have smaller temperature differences between the heated and unheated sides, and intra-binary shocks that wrap around the pulsar, rather than the companion star<sup>11</sup>, making these effects less strong for these systems. Optical observations of the remaining black widow, PSR J1555–2908, have been investigated with a model that takes heat diffusion into account<sup>30</sup>, resulting in inclination constraints that are consistent with our eclipse detection.

We have not found a case in which previous optical modelling suggested a high inclination (and therefore a low pulsar mass) that is now ruled out by the non-detection of a gamma-ray eclipse. For PSR J2215+5135, Romani, Graham, Filippenko and Kerr<sup>31</sup> inferred a high inclination from optical modelling, and even found marginal evidence for a low-significance gamma-ray eclipse in the Fermi-LAT data. However, eclipses in this system lasting longer than 0.1% of an orbit are strongly ruled out by our longer dataset, showing that this earlier hint was likely to be a chance false-alarm, and indeed more recent modelling by Linares, Shahbaz and Casares<sup>7</sup> and Kandel and Romani<sup>32</sup> found lower, non-eclipsing inclinations.

Finally, we note that pulsars that are eclipsed by their companion stars also necessarily pass in front of the heated face of the companion star half an orbit later. As neutron stars are very small in size compared with their companion stars, but have intense gravitational fields, they will act as gravitational lenses, magnifying the optical flux from the companion star (for example, refs. 33–35). The exact degree of the magnification depends only on the pulsar mass and the orbital separation. The detection of this gravitational lensing would therefore provide an independent measurement of the neutron star mass. Unfortunately, the magnification due to lensing is expected to be on the order of $10^{-3}$ mag (ref. 33). Effects of this level can be dwarfed by both short- and long-timescale variability on the order of 0.1 mag (refs. 13,36), as well as by systematic uncertainties in the underlying light curve due to incomplete heating models. Detecting the lensing effect will therefore require extremely sensitive optical photometry, stacked over several orbits to average out variability, and careful modelling to disentangle this effect from underlying heating effects.
# Methods

## Gamma-ray observations

For each system in our sample, we analysed 11.4 yr of observations taken by the Fermi LAT. We selected SOURCE-class photons detected with reconstructed energies 50 MeV < E < 300 GeV, and with reconstructed directions from within a 3° region of interest around each pulsar, according to the P8R3_SOURCE_V2 instrument response functions³⁷,³⁸.

Sensitive unbinned-likelihood-based methods for detecting eclipses (for example ref. ³⁹) account for each photon individually, and therefore must account for the relative probability of each photon having been emitted by the target source, as opposed to by a fore/ background source. This is achieved by weighting the contribution of each photon to the relevant statistic⁴⁰. Computing these weights requires an accurate spectral and spatial model of the emission from the target pulsar and all fore/background sources in the region of interest⁴¹. For this, we used the 10 yr incremental version (Data Release 2) of the Fermi-LAT Fourth Source Catalog (4FGL,⁴²,⁴³) (https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr catalogue/) and the gll_iem_v07.fits Galactic diffuse emission and iso_P8R3_SOURCE_V3_v1.txt isotropic diffuse emission models to describe the diffuse background emission. The parameters of the spectra of the target pulsars were then refined such that the resulting photon weights maximized the significances of their gamma-ray pulsations, as described in ref. ⁴¹. These photon weights made use of the 'PSF' event types (https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone Data/LAT_DP.html) to benefit from the narrower point-spread function for well-reconstructed photon events.

The required timing ephemerides for each spider pulsar were compiled as part of an upcoming third iteration of the Fermi-LAT Pulsar Catalogue⁴⁴. For each pulsar, we computed the orbital phase at which each photon was emitted according to these ephemerides using the TEMPO2 software⁴⁵. The validity of the orbital ephemeris was verified by the presence of gamma-ray pulsations throughout the data.

## Test statistic for eclipse detection

To test for possible eclipses we adopted the unbinned-likelihood estimation methods described in ref. ¹⁸ and ref. ³⁹. Under this model, we assumed that the eclipse has sharp in/egresses, is centred on the pulsar's superior conjunction, and lasts for a fraction $\theta$ of the orbital period. Within the eclipse, we assumed a constant flux level, which we parameterized with $\alpha$, the fractional flux level within the eclipse relative to the overall average flux. The increase in $\log \mathcal{L}$ for such an eclipse, compared with the null hypothesis of photons being uniformly distributed in orbital phase, is:

$$ \delta \log \mathcal{L}(\alpha, \theta) = \sum_{i \in \Theta} \log(w_i \alpha + 1 - w_i) + \sum_{i \in \bar{\Theta}} \log\left(w_i \frac{1-\alpha\theta}{1-\theta} + 1 - w_i\right) - \left(\alpha\eta_{\Theta} + \frac{1-\alpha\theta}{1-\theta}\eta_{\bar{\Theta}} - 1\right) \sum_i w_i \quad (1) $$

where $w_i$ is the photon probability weight for the ith photon; $\Theta(\bar{\Theta})$ refers to photons with orbital phases inside (outside) the eclipse; and $\eta_\theta(\eta_{\bar{\Theta}})$ denotes the fractional exposure inside (outside) the eclipse. The last term in equation (1) accounts for variations in the exposure as a function of orbital phase. We computed the exposure for each pulsar in 30 s time intervals over the Fermi mission using godot³⁹, and folded these on the orbital period to compute $\eta_\theta$ and $\eta_{\bar{\Theta}}$. After several years of observations, corresponding to several thousands of orbits of each pulsar system included here, the exposure is usually very evenly distributed across all orbital phases and hence exposure variations typically have very little effect on the resulting likelihood calculation; we correct for this effect nevertheless.

For each pulsar, we tested the hypothesis of a complete eclipse of the gamma-ray emission, corresponding to $\alpha = 0$, testing for $\theta \in [0, 0.2)$ with fine spacing. The upper bound on this range was more than twice as large as the maximum possible eclipse duration for our studied population, assuming that companion stars do not overflow their Roche lobes. The pulsar with the smallest mass ratio in our population is PSR J2129-0429, with $q = M_{\text{psr}}/M_c = 3.93 \pm 0.06$, which would eclipse for 8.4% of an orbit if the companion filled its Roche lobe and was observed at $i = 90^\circ$.

For pulsars in which significant eclipses were detected we also tested alternative eclipse models with $0 < \alpha < 1$ and curved (rather than sharp) ingresses and egresses. No significant $\log \mathcal{L}$ improvements were observed.

## Posterior photon weights

A significant improvement in sensitivity when searching for eclipses could be obtained by incorporating into our analysis the fact that the gamma-ray emission is pulsed; that is, that gamma-ray photons observed at pulse phases that fall within a peak in the gamma-ray pulse profile are more likely to have originated from the pulsar than from the background. To make use of this knowledge, we used the photon reweighting method of Kerr³⁹, which we briefly describe here.

A photon weight, $w$, computed as above using the spectral and spatial model of the region of interest, is our best estimate for the probability of that photon having been emitted by the pulsar, before including knowledge of the pulsar rotational phase at which the photon was emitted. We can denote this as a prior probability $P(S) = w$, where $S$ denotes the binary statement that the photon was emitted by our target source. The probability for the opposite case, $B$, where the photon is emitted by a background source, is then $P(B) = 1 - w$. The reweighting method updates our knowledge of the probability of the photon being emitted by the target source on the basis of the rotational phase $\phi$ at which the photon was emitted by applying Bayes' theorem:

$$ P(S|\phi) = \frac{p(\phi|S)P(S)}{p(\phi|S)P(S) + p(\phi|B)P(B)} \quad (2) $$

Here $P(S|\phi)$ is now the posterior probability of the photon having been emitted by the pulsar, given its rotational phase; $p(\phi|S)$ is the phase distribution of photons emitted by the pulsar (that is, the pulsar's pulse profile, which we hereafter denote as $f(\phi)$); and $p(\phi|B)$ is the phase distribution of background photons, which we could safely assume to be uniform when folding on the millisecond pulse periods of the pulsars included here, hence $p(\phi|B) = 1$. Rewriting equation (2) with these values gives us the reweighting equation:

$$ P(S|\phi) = w' = \frac{w f(\phi)}{w f(\phi) + 1 - w} \quad (3) $$

Hereafter $w$ is the prior weights, and $w'$ is the posterior weights. For phases within peaks of the pulse profile, where $f(\phi) > 1$, these posterior weights are always greater than the prior weights, and for phases outside of peaks, where $f(\phi) < 1$, the posterior weights are always lower. Thus, photons within pulse peaks are upweighted, while the rest are downweighted. When searching for eclipses, the posterior weights help to increase the detection statistic values for true eclipses by downweighting the detrimental effect of photons that by chance have high weights, and fall within the eclipse region, but whose rotational phases do not lie within a pulse peak and are therefore less likely to have been emitted by the pulsar than initially predicted by the prior weight. Similarly, photons lying outside the eclipse region but within a pulse peak, and therefore more likely to have been emitted by the pulsar, have a larger positive contribution to the eclipse $\log \mathcal{L}$.

To obtain the pulse profile models, $f(\phi)$, we fitted a set of wrapped Gaussian functions to the prior-weighted photon phases using the
maximum-likelihood method described by Abdo et al.⁴⁴. The number of Gaussian functions used to model each pulse profile was chosen to minimize the Bayesian information criterion⁴⁶.

We initially performed our search using the prior weights, but changed to using the posterior weights after finding that they significantly improved the sensitivity to eclipses. Of the four significant eclipses that were found using the posterior weights for the eclipse search (prior weights were used for PSR J1048+2339 as discussed below), three were originally significantly detected with the prior weights, but the posterior weights gave significantly larger log L values, with $\delta \log L$ increasing by at least 2.6 for these pulsars. Only the eclipse from PSR J1555-2908 was undetected using the prior weights, with $\delta \log L = 3.35$ compared with $\delta \log L = 10.07$ with posterior weights, probably owing to its weak overall flux but very narrow pulse peaks.

For one pulsar in which a significant eclipse is found, PSR J1048+2339, the radio timing ephemeris only covers a shorter 3 yr portion of the LAT mission, with variations in the orbital period and a low photon flux preventing the generation of a full timing ephemeris using the LAT data. The radio timing ephemeris for this pulsar contains several orbital frequency derivatives to model these variations, but this ephemeris becomes highly uncertain when extrapolating outside the time interval in which it was derived. For our study, we removed these orbital frequency derivatives from the ephemeris and computed orbital phases assuming a constant orbital period. During the radio timing interval, these orbital period variations cause orbital phase shifts of up to $-10^{-3}$ orbits²⁰. This is around 2% of the duration of the eclipse detected in this system, and therefore we did not expect this additional source of uncertainty to substantially affect our results. For this pulsar, given that pulsations were not observed outside the period covered by the radio ephemeris, we used $w$, rather than $w'$, when searching for eclipses.

We included the two gamma-ray detected transitional MSPs, PSRs J1023+0038 and J1227-4853 in our study, classifying these as redbacks, as they seem to be very similar to this class when in their non-accreting state. We note that the source of their increased gamma-ray flux during the accreting states is unclear, but we assumed that it also originates close to the neutron star (as seems to be the case for the gamma-ray eclipsing transitional MSP candidate 4FGL J0427.8-6704; refs. ¹⁷,¹⁸), and included data from both the accreting- and non-accreting states in our analysis. Pulsations are not detected from these pulsars in their accreting states, and so for these we again used the prior weights, rather than the posterior weights.

## Significance calibration via Monte Carlo simulations

The search over the eclipse width $\theta$ introduced an unknown number of independent trials to our search. We therefore calibrated false-alarm probabilities ($P_{FA}$) via Monte Carlo analysis. For each pulsar, we took the observed set of posterior weights, randomly sampled orbital phases from a uniform distribution, computed the log-likelihood of equation (1) for the same set of $\theta$ values used in the eclipse search and stored the maximum value, iterating $10^7$ times.

Figure 2 shows the results of the Monte Carlo simulations that we used to calibrate the statistical significances of these eclipses. Of the five eclipsing pulsars, the eclipse in PSR J1555-2908 has the lowest significance, but still has $P_{FA} \approx 5 \times 10^{-5}$.

The $\delta \log L$ values observed from PSRs J1816+4510 and J2129-0429 are larger than any obtained in our simulations. To estimate their false-alarm probabilities, we therefore performed a simple linear fit to the observed $\delta \log L$ versus $\log(P_{FA})$ curves for these pulsars, and extrapolated to the observed values.

In Fig. 2, we also show the empirical survival function; that is, the fraction of pulsars whose measured eclipse log-likelihoods would survive a given threshold. If there were no eclipses in our dataset, then the set of measured log-likelihood values would be drawn from the null-hypothesis distribution and this empirical survival function

curve would closely follow the simulated curves. The ratio between the empirical and simulated curves at the highest measured log-likelihood value illustrates the significance of the largest outlier, given the number of pulsars included in the sample.

If we removed the five eclipsing systems, then the empirical survival function curve did closely follow the simulated null-hypothesis curve, and only started to deviate for the final pulsar, PSR J0251+2606, which had a false-alarm probability of around 0.2%. With $n_{psr} = 37$ pulsars remaining in this sample, the largest outlier should have a false-alarm probability of around $1/n_{psr} = 1/37 \approx 2.7\%$. This pulsar therefore has an eclipse log-likelihood value that has an estimated false-alarm rate around ten times lower than expected for the largest outlier from our study, given the number of pulsars included. This could be viewed as marginal evidence for an eclipse, with all other measured $\delta \log L$ values being consistent with the null hypothesis.

## Pulsar mass constraints

The significant detection or exclusion of a gamma-ray eclipse provides a constraint on the binary inclination angle that depends on the angular size of the companion star as seen from the pulsar. The angular size of the companion star's Roche lobe only depends on the binary mass ratio, $q$, and hence it is convenient to parameterize the size of the companion star by $q$, and its Roche-lobe filling factor $f_{RL}$ (which we define as the radius of the star along the binary separation vector divided by the Roche lobe radius in the same direction). These parameters can be constrained by optical observations. The mass ratio was derived from measurements of the pulsar and companion projected radial velocity amplitudes ($K_{psr}$ and $K_c$), measured via pulsar timing and optical spectroscopy, respectively, with $q = M_{psr}/M_c = K_c/K_{psr}$. The Roche-lobe filling factor can be estimated from rotational broadening or surface gravity measurements via optical spectroscopy (for example ref. ⁴⁷) or from the amplitude of the 'ellipsoidal' component of an observed optical light curve⁴⁸. However, this parameter is often correlated with the estimated inclination, and so previous estimates of $f_{RL}$ are not necessarily consistent with new inclination constraints from an eclipse detection or exclusion.

To compute expected eclipse durations, we generated model stars using the Icarus⁴⁸ binary modelling software. Icarus assumes that the surface of the star follows an equipotential contour within its Roche lobe, and therefore the simulated surface accounts for the non-spherical shape of the star due to tidal and rotational deformation. For a given binary mass ratio and Roche-lobe filling factor we could then compute the range of orbital phases at which the pulsar is eclipsed by the model star, when viewed from a given inclination. We assumed that the pulsar is effectively a point-source of gamma-ray emission, since gamma-ray emission is thought to either be produced inside, or just outside, the pulsar's light cylinder⁴⁹, which is thousands of times smaller than the orbital separation in a spider binary. As we do not detect gradual in/egresses in the eclipses, and since the density profile of the outer envelope of the companion star is unknown, we assumed that any line-of-sight crossing the photosphere will be fully eclipsed.

The pulsar and companion masses can be estimated, as a function of inclination, from the binary mass function,

$$M_{psr} = \frac{K_c^2 P_{orb}(1 + 1/q)^2}{2\pi G \sin^3 i}, \qquad (4)$$

$$M_c = \frac{K_{psr}^3 P_{orb}(1 + q)^2}{2\pi G \sin^3 i}. \qquad (5)$$

Tables 1 and 2 list our resulting constraints on the inclination and component masses for eclipsing and non-eclipsing systems, respectively, with existing companion radial velocity measurements. When an eclipse is detected, we assumed that the companion fills its Roche
lobe to obtain a lower bound on the inclination, and hence a conservatively high upper bound on the pulsar and companion masses, while assuming that $i = 90^\circ$ provided a strict lower limit on the masses with no assumption on the filling factor. For systems without detected eclipses, we assumed a low $f_{RL} = 0.5$ to obtain an upper bound on the inclination, and lower bound on the component masses. This limit was based on the low filling factor for PSR J1816+4510 estimated by Kaplan et al.⁴⁷ using the surface gravity determined by optical spectroscopy. Optical models for black widow and redback systems tend to have higher estimated filling factors (for example ref. ²⁴), and so we adopted this value as a conservative estimate.

Where possible, we took radial velocity amplitudes that had been corrected for heating effects that shift the centre of light away from the companion's centre of mass. These corrections tended to increase the inferred $K_c$, and hence increase the pulsar mass estimate. Three redback pulsars in our list do not have published centre-of-light-corrected radial velocity amplitudes (PSRs J1431−4715, J1622−0315 and J1625−3205), and so the mass limits for these may be slightly underestimated. However, all three have optical light curves that suggest very little heating effect is present, so the required corrections are likely to be small for these systems.

All systems studied here that do not have radial velocity measurements are black widows (which tend to be fainter at optical wavelengths, and hence often inaccessible to spectroscopic studies). For these, we assumed a typical mass ratio of $q = 70$, and list the resulting inclination constraints in Supplementary Table 2. Although the inclination constraints vary slowly with $q$ at values typical for black widows, the component masses do depend strongly on the assumed value of $q$, and so we do not list mass constraints here.

## Optical constraints for eclipsing pulsars

Previous results from optical observations and modelling of PSR B1957+20 are discussed in the main text. In the paragraphs below we discuss the existing multiwavelength observations for the other four eclipsing pulsars. Where previous works provide constraints on the Roche-lobe filling factor, we used these constraints to obtain larger (but less robust) lower limits on the inclination angle (and therefore tighter upper bounds on the pulsar masses) than were obtained by assuming $f_{RL} = 1$ in the previous section. We also used these estimates to obtain upper limits on the inclination angle, rather than simply assuming $i < 90^\circ$ (which may imply a very low filling factor).

Optical observations of PSR J1048+2339 show significant long-term variability, with the optical maximum varying by up to one magnitude¹³,⁵⁰. Such variability cannot yet be taken into account by precise light curve modelling, and so no measurement of the inclination angle from optical modelling exists in the literature. When modelling their optical observations of this pulsar, Yap et al.⁵⁰ fixed the inclination to $i = 76^\circ$, the maximum value that was compatible with the lack of an observed X-ray eclipse. However, while a thermal X-ray component from the neutron star surface would indeed be eclipsed at higher inclinations, X-ray emission in redbacks tends to be dominated by emission from an extended intra-binary shock, and so the lack of an X-ray eclipse does not necessarily preclude a higher inclination. Yap et al.⁵⁰ found that $f_{RL} \approx 0.85$ was compatible with multiple light curves despite long-term variability. From optical spectroscopy, Miraval Zanon et al.⁵¹ found an observed companion radial velocity amplitude of $343.3 \pm 4.4$ km s⁻¹. Using these values, we found that the observed eclipse duration required an inclination greater than 80.⁹° (compare with 80.4° assuming $f_{RL} = 1$ in Table 1). Heating corrections reduced the estimated centre-of-mass velocity to $298.7 \pm 7.7$ km s⁻¹, for a larger mass ratio (and hence lower minimum inclination of $i = 80.1^\circ$) but a lower pulsar mass, $M_{psr} \approx 1.1 M_\odot$. We used the uncorrected value of $K_c$ in Table 1 to obtain a conservative bound on the pulsar mass. As noted in the main text, eclipses longer than 8% of an orbit were consistent with the data, but would require the companion star to be overflowing its Roche lobe.

PSR J1555−2908 is a black widow pulsar that was recently discovered by Ray et al.⁵² in a targeted radio search of a steep-spectrum radio continuum source identified within a pulsar-like gamma-ray source by Frail et al.⁵³. Modelling of both optical photometry and spectroscopy by Kennedy et al.³⁰, using a model that includes the possibility of heat diffusion across the terminator, revealed the companion's projected radial velocity to be $397 \pm 2$ km s⁻¹, and indicated a high binary inclination of $i > 75^\circ$, giving a maximum pulsar mass of $1.82 M_\odot$. The Roche-lobe filling factor was found to be high, $f_{RL} > 0.93$. The duration of the gamma-ray eclipse observed here requires an inclination $83.0^\circ < i < 86.2^\circ$, for a pulsar mass $1.58 M_\odot < M_{psr} < 1.71 M_\odot$ with the uncertainty dominated by that of the companion's radial velocity.

Optical spectroscopy of PSR J1816+4510 has been modelled by Kaplan et al.⁴⁷. They found that this system is perhaps more similar to a white-dwarf companion than a normal redback, owing to its extremely high temperature, but due to the presence of radio eclipses that are not otherwise seen in pulsar–white-dwarf binaries we categorize it here as the latter. Detailed modelling of optical photometry to determine the inclination or Roche-lobe filling factor has not been performed, but from spectroscopic models Kaplan et al.⁴⁷ estimated a radius that corresponds to $f_{RL} \approx 0.5$, which is much smaller than observed in other redbacks, motivating our use of this value as a low estimate for $f_{RL}$ in Table 2. Adopting this value instead of $f_{RL} = 1$ results in a higher minimum inclination of $82.6^\circ$ and a lower pulsar mass range of $1.68 M_\odot < M_{psr} < 2.11 M_\odot$.

For PSR J2129−0429, Bellm et al.⁵⁴ measured a projected companion radial velocity amplitude of $K_2 = 250 \pm 4$ km s⁻¹ (for mass ratio $q = 3.93 \pm 0.06$), and inferred a filling factor of $f_{RL} = 0.82 \pm 0.03$ and an inclination $i > 68^\circ$ from optical light curve modelling. With these values of $q$ and $f_{RL}$, the duration of the gamma-ray eclipse requires an inclination between $76.6^\circ < i < 78.3^\circ$, consistent with the range allowed by optical modelling. This corresponds to a pulsar mass range of $1.61 M_\odot < M_{psr} < 1.88 M_\odot$ at 95% confidence. Al Noori et al.⁵⁵ also found dips in the XMM-Newton light curve for PSR J2129−0429, consistent with a thermal X-ray component from the neutron star surface being eclipsed by the companion.

## Searching for eclipses in recently discovered redbacks and candidates

Seven Fermi-LAT sources have been found to contain periodic optical and X-ray sources that are almost certainly spider binary systems²¹,⁵⁶⁻⁶². Shortly before submitting this Article, millisecond radio pulsations were detected from three of these objects (4FGL J0838.7–2827, 4FGL J0955.3–3949 and 4FGL J2333.1–5527, hereafter PSRs J0838–2827, J0955–3949 and J2333–5526, respectively) by the TRAPUM collaboration (http://trapum.org/discoveries.html), but a full timing solution is not yet available for them, and gamma-ray pulsations have not yet been detected. Pulsations have not yet been detected at any wavelength from the remaining four of these systems. A further four similar systems¹⁴,²²,⁶³,⁶⁴ were initially identified in the same way, but have since been confirmed as spiders through radio or gamma-ray pulsation discoveries, and hence are already included in our search.

To search for gamma-ray eclipses in these systems, we prepared Fermi-LAT datasets in the same way as for the confirmed spider cases, but included 12.4 yr of data. As these datasets were not bound to the validity period of a pulsar timing ephemeris, we included this extra year of data to allow for stronger detections to partially mitigate the large trials factor (see below). For these datasets we used gtsrcprob to compute the photon weights, rather than optimizing these to maximize the pulsation significance (since this optimization is not possible here without a gamma-ray pulsation detection), and used the prior probability weights, since posterior weights are unavailable in the absence of pulsations.

Unlike the majority of systems studied here, where the pulsar's timing ephemeris provides precise orbital period and phase measurements over the Fermi-LAT dataset, for these systems we have only
imprecise orbital phase information from optical light curves and radial velocity curves. We therefore had to also search small ranges of orbital phases (parameterized by each pulsar's ascending node epoch, $T_{asc}$ and orbital period, $P_{orb}$). We chose the search ranges to be ±3× the published uncertainty on these parameters. Our step sizes in each parameter were chosen such that the maximum offset on the orbital phase for each photon would be 0.001 orbits.

This searching introduced yet more trials in our search, reducing our sensitivity. For each of these sources, we again calibrated our search significances via Monte Carlo simulations. We took the observed set of photon phases, assigned randomly generated arrival times evenly distributed throughout the Fermi mission interval, performed the search over orbital period and phase as above, and took the maximum resulting $\delta \log L$, iterating 1,000 times (or 10,000 times if the resulting false-alarm probability was low).

Sensitivity was greatly reduced in these searches, with $P_{FA} = 0.01$ corresponding to $\delta \log L \approx 9$, as opposed to $\delta \log L \approx 4.5$ when the pulsar's precise orbital ephemeris is known. Nevertheless, there is evidence of eclipses in two systems, PSRs J0838−2827 and J2333−5526, with $\delta \log L = 10.8$ and $\delta \log L = 11.3$, corresponding to trial-corrected $P_{FA}$ of $3 \times 10^{-3}$ and $1 \times 10^{-3}$, respectively. We show the $\delta \log L$ over the searched parameter space for these two systems and the corresponding Monte Carlo calibrations in Fig. 3.

For both PSR J0838−2827 and PSR J2333−5526, although significant detections were made within the searched parameter space, the eclipse likelihoods were still high towards the borders of the searched regions. We therefore searched outside this region and found higher $\delta \log L$ values ($\delta \log L = 11.8$ and $\delta \log L = 14.9$, respectively) at orbital phases that are $\Delta T_{asc} \approx 300$ s and $\Delta T_{asc} \approx 620$ s later than that predicted by the ephemerides obtained from fitting the companions' radial velocity curves, corresponding to -4σ and -11σ deviations, respectively. Such offsets may be caused by the irradiation of the companion star by the pulsar, which can cause the radial velocity curve to depart slightly from a simple sinusoid, although Swihart et al.⁶⁰ did not find any evidence for this effect in their modelling of PSR J2333−5526. The recent detection of radio pulsations from these pulsars will probably clarify these tensions by providing precise orbital ephemerides.

## Expected number of eclipsing spiders

Although our population is small, we could also use the number of detected eclipses to probe whether or not the population of Fermi-LAT detected spiders are viewed from randomly distributed inclination angles. This is not necessarily expected; as many of these sources have been discovered by targeting Fermi-LAT sources, the population may be biased towards those that are bright gamma-ray emitters, and gamma-ray emission models predict that MSPs are brightest around their rotational equator, which should in turn be aligned with the orbital plane during recycling. This would manifest in our population as a greater than expected number of eclipsing pulsars. Alternatively, if we observed a smaller number of eclipses than expected, this would be evidence that the companion stars in these systems tend to fill only a small fraction of their Roche lobes.

Under the assumption of randomly distributed orbital axes, the binary inclination angles will be drawn from a probability distribution $p(i) = \sin(i)$, which we adopted as a prior. We restricted $i$ to $i \le 90^\circ$, as systems at $i$ are indistinguishable from those at $180^\circ - i$ because the orbital direction cannot be determined. The prior probability of a pulsar being eclipsed is therefore the integral of this prior over inclination angles greater than the minimum inclination at which a pulsar would be eclipsed, $i_{ecl}$; that is, $P(i > i_{ecl}) = \cos(i_{ecl})$. Assuming Roche-lobe filling companions, for typical black widow and redback mass ratios ($q \equiv M_{\text{psr}}/M_{\text{c}}$) of $q = 70$ and $q = 5$, respectively, this gives a prior probability for a black widow or redback being eclipsed of 10% or 23%, respectively. The probability of observing a certain number of eclipses from the studied population follows a binomial distribution with these

success factors. From the 28 black widows and 16 rebacks in our sample (including the candidates discussed in the previous section) that are bright enough for us to significantly detect or rule out an eclipse, we found 2 eclipsing black widows and 5 eclipsing rebacks. The binomial probabilities for these samples are 24% and 16% respectively, entirely consistent with our assumptions of randomly distributed inclinations and Roche-lobe filling companions, and indeed 7 eclipses is the second most likely number of eclipses to observe from the combined population.

## Data availability

The Fermi-LAT data are available from the Fermi Science Support Center at http://fermi.gsfc.nasa.gov/ssc. Ephemerides and folded Fermi-LAT datasets, including prior and posterior photon weights for systems with detected eclipses, are available via Zenodo at https://doi.org/10.5281/zenodo.7133502. Ephemerides and folded datasets for other pulsars included in this study may contain unpublished information about unrelated scientific results; these are available from the corresponding author upon request.

## Code availability

The Fermitools, including gtsrcprob, used to analyse Fermi-LAT data are available at https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/. TEMPO⁴⁵, used to compute rotational and orbital photon phases, is available via Bitbucket at https://bitbucket.org/psrsoft/tempo2. Icarus⁴⁸, used to compute expected eclipse durations to derive pulsar mass estimates, is available via GitHub at https://github.com/bretonr/Icarus. Godot³⁹, used to compute Fermi-LAT exposure, is available via GitHub at https://github.com/kerrm/godot. PINT⁶⁵, used to evaluate template pulse profiles to derive posterior weights, is available via GitHub at https://github.com/nanograv/PINT. The scripts used to perform the eclipse searches and false-alarm calibrations are available via Zenodo at https://doi.org/10.5281/zenodo.7133502. These scripts made extensive use of the open source NumPy⁶⁶, SciPy⁶⁷, Matplotlib⁶⁸, Astropy⁶⁹,⁷⁰ and PyCUDA⁷¹ Python packages.

## References

1. Fruchter, A. S., Stinebring, D. R. & Taylor, J. H. A millisecond pulsar in an eclipsing binary. *Nature* **333**, 237–239 (1988).

2. Roberts, M. S. E. (2013) Surrounded by spiders! New black widows and rebacks in the Galactic field. In van Leeuwen, J. (ed) *Proc. International Astronomical Union* vol. 8. Cambridge University Press, Cambridge, pp 127–132

3. Stappers, B. W. et al. Probing the eclipse region of a binary millisecond pulsar. *Astrophys. J. Lett.* **465**, L119 (1996).

4. Polzin, E. J. et al. Study of spider pulsar binary eclipses and discovery of an eclipse mechanism transition. *Mon. Not. R. Astron. Soc.* **494**, 2948–2968 (2020).

5. van Kerkwijk, M. H., Breton, R. P. & Kulkarni, S. R. Evidence for a massive neutron star from a radial-velocity study of the companion to the black-widow pulsar PSR B1957+20. *Astrophys. J.* **728**, 95 (2011).

6. Romani, R. W., Filippenko, A. V. & Cenko, S. B. A spectroscopic study of the extreme black widow PSR J1311-3430. *Astrophys. J.* **804**, 115 (2015).

7. Linares, M., Shahbaz, T. & Casares, J. Peering into the dark side: magnesium lines establish a massive neutron star in PSR J2215+5135. *Astrophys. J.* **859**, 54 (2018).

8. Strader, J. et al. Optical spectroscopy and demographics of redback millisecond pulsars binaries. *Astrophys. J.* **872**, 42 (2019).

9. Linares, M. Super-massive neutron stars and compact binary millisecond pulsars. *Proc. Sci. MULTIF2019*, 023 (2020).

10. Özel, F. & Freire, P. Masses, radii, and the equation of state of neutron stars. *Annu. Rev. Astron. Astrophys.* **54**, 401–440 (2016).
11. Romani, R. W. & Sanchez, N. Intra-binary shock heating of black widow companions. Astrophys. J. 828, 7 (2016).

12. van Staden, A. D. & Antoniadis, J. An active, asynchronous companion to a redback millisecond pulsar. Astrophys. J. Lett. 833, L12 (2016).

13. Cho, P. B., Halpern, J. P. & Bogdanov, S. Variable heating and flaring of three redback millisecond pulsar companions. Astrophys. J. 866, 71 (2018).

14. Clark, C. J. et al. Einstein@Home discovery of the gamma-ray millisecond pulsar PSR J2039-5617 confirms its predicted redback nature. Mon. Not. R. Astron. Soc. 502, 915-934 (2021).

15. Atwood, W. B. et al. The Large Area Telescope on the Fermi gamma-ray space telescope mission. Astrophys. J. 697, 1071-1102 (2009).

16. Ray, P. S. et al. Radio searches of Fermi LAT sources and blind search pulsars: the Fermi Pulsar Search Consortium. In Morselli, A. (ed) Proc. 2011 Fermi Symposium eConf C110509 (2012).

17. Strader, J. et al. A new γ-ray loud, eclipsing low-mass X-ray binary. Astrophys. J. 831, 89 (2016).

18. Kennedy, M. R. et al. Optical, X-ray, and γ-ray observations of the candidate transitional millisecond pulsar 4FGL J0427.8-6704. Mon. Not. R. Astron. Soc. 494, 3912-3926 (2020).

19. Petrosian, V., McTiernan, J. M. & Marschhauser, H. Gamma-ray emission and electron acceleration in solar flares. Astrophys. J. 434, 747 (1994).

20. Deneva, J. S. et al. Multiwavelength observations of the redback millisecond pulsar J1048+2339. Astrophys. J. 823, 105 (2016).

21. Halpern, J. P., Strader, J. & Li, M. A likely redback millisecond pulsar counterpart of 3FGL J0838.8-2829. Astrophys. J. 844, 150 (2017).

22. Nieder, L. et al. Discovery of a gamma-ray black widow pulsar by GPU-accelerated Einstein@Home. Astrophys. J. Lett. 902, L46 (2020).

23. Reynolds, M. T. et al. The light curve of the companion to PSR B1957+20. Mon. Not. R. Astron. Soc. 379, 1117-1122 (2007).

24. Draghis, P. et al. Multiband optical light curves of black-widow pulsars. Astrophys. J. 883, 108 (2019).

25. Wadiasingh, Z. et al. Constraining relativistic bow shock properties in rotation-powered millisecond pulsar binaries. Astrophys. J. 839, 80 (2017).

26. Johnson, T. J. et al. Constraints on the emission geometries and spin evolution of gamma-ray millisecond pulsars. Astrophys. J. Suppl. Ser. 213, 6 (2014).

27. Cromartie, H. T. et al. Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar. Nat. Astron. 4, 72-79 (2019).

28. Guillemot, L. et al. Pulsed gamma rays from the original millisecond and black widow pulsars: a case for caustic radio emission? Astrophys. J. 744, 33 (2012).

29. Voisin, G. et al. A model for redistributing heat over the surface of irradiated spider companions. Mon. Not. R. Astron. Soc. 499, 1758-1768 (2020).

30. Kennedy, M. R. et al. Measuring the mass of the black widow PSR J1555-2908. Mon. Not. R. Astron. Soc. 512, 3001-3014 (2022).

31. Romani, R. W., Graham, M. L., Filippenko, A. V. & Kerr, M. Keck spectroscopy of millisecond pulsar J2215+5135: a moderate-$M_{NS}$, high-inclination binary. Astrophys. J. Lett. 809, L10 (2015).

32. Kandel, D. & Romani, R. W. Atmospheric circulation on black widow companions. Astrophys. J. 892, 101 (2020).

33. Marsh, T. R. Gravitational lensing in eclipsing binary stars. Mon. Not. R. Astron. Soc. 324, 547-552 (2001).

34. Sahu, K. C. & Gilliland, R. L. Near-field microlensing and its effects on stellar transit observations by Kepler. Astrophys. J. 584, 1042-1052 (2003).

35. Beskin, G. M. & Tuntsov, A. V. Detection of compact objects by means of gravitational lensing in binary systems. Astron. Astrophys. 394, 489-503 (2002).

36. Romani, R. W. et al. PSR J1311-3430: a heavyweight neutron star with a flyweight helium companion. Astrophys. J. Lett. 760, L36 (2012).

37. Atwood, W. et al. Pass 8: toward the full realization of the Fermi-LAT scientific potential. In Brandt, T. J. et al. (eds) Proc. 4th Fermi Symposium eConf, C121028, pp 8-13 (2012).

38. Bruel, P. et al. Fermi-LAT improved Pass 8 event selection. Preprint at https://arxiv.org/abs/1810.11394 (2018).

39. Kerr, M. Multiscale time- and frequency-domain likelihood analysis with photon weights. Astrophys. J. 885, 92 (2019).

40. Kerr, M. Improving sensitivity to weak pulsations with photon probability weighting. Astrophys. J. 732, 38 (2011).

41. Bruel, P. Extending the event-weighted pulsation search to very faint gamma-ray sources. Astron. Astrophys. 622, A108 (2019).

42. Abdollahi, S. et al. Fermi Large Area Telescope fourth source catalog. Astrophys. J. Suppl. Ser. 247, 33 (2020).

43. Ballet, J., Burnett, T. H., Digel, S. W. & Lott, B. Fermi Large Area Telescope fourth source catalog data release 2. Preprint at https://arxiv.org/abs/2005.11208 (2020).

44. Abdo, A. A. et al. The second Fermi Large Area Telescope catalog of gamma-ray pulsars. Astrophys. J. Suppl. Ser. 208, 17 (2013).

45. Edwards, R. T., Hobbs, G. B. & Manchester, R. N. TEMPO2, a new pulsar timing package - II. The timing model and precision estimates. Mon. Not. R. Astron. Soc. 372, 1549-1574 (2006).

46. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 6, 461-464 (1978).

47. Kaplan, D. L. et al. A metal-rich low-gravity companion to a massive millisecond pulsar. Astrophys. J. 765, 158 (2013).

48. Breton, R. P., Rappaport, S. A., van Kerkwijk, M. H. & Carter, J. A. KOI 1224: a fourth bloated hot white dwarf companion found with Kepler. Astrophys. J. 748, 115 (2012).

49. Kalapotharakos, C., Harding, A. K., Kazanas, D. & Wadiasingh, Z. A fundamental plane for gamma-ray pulsars. Astrophys. J. Lett. 883, L4 (2019).

50. Yap, Y. X. et al. Face changing companion of the redback millisecond pulsar PSR J1048+2339. Astron. Astrophys. 621, L9 (2019).

51. Miraval Zanon, A. et al. Evidence of intra-binary shock emission from the redback pulsar PSR J1048+2339. Astron. Astrophys. 649, A120 (2021).

52. Ray, P. S. et al. Discovery, timing, and multiwavelength observations of the black widow millisecond pulsar PSR J1555-2908. Astrophys. J. 927, 216 (2022).

53. Frail, D. A. et al. An image-based search for pulsars among Fermi unassociated LAT sources. Mon. Not. R. Astron. Soc. 475, 942-954 (2018).

54. Bellm, E. C. et al. Properties and evolution of the redback millisecond pulsar binary PSR J2129-0429. Astrophys. J. 816, 74 (2016).

55. Al Noori, H. et al. X-Ray and optical studies of the redback system PSR J2129-0429. Astrophys. J. 861, 89 (2018).

56. Strader, J. et al. 1FGL J0523.5-2529: a new probable gamma-ray pulsar binary. Astrophys. J. Lett. 788, L27 (2014).

57. Li, K. L. et al. Discovery of a redback millisecond pulsar candidate: 3FGL J0212.1+5320. Astrophys. J. 833, 143 (2016).

58. Linares, M. et al. A millisecond pulsar candidate in a 21-h orbit: 3FGL J0212.1+5320. Mon. Not. R. Astron. Soc. 465, 4602-4610 (2017).

59. Li, K. L. et al. Multiwavelength observations of a new redback millisecond pulsar candidate: 3FGL J0954.8-3948. Astrophys. J. 863, 194 (2018).
60. Swihart, S. J. et al. A new likely redback millisecond pulsar binary with a massive neutron star: 4FGL J2333.1-5527. Astrophys. J. 892, 21 (2020).

61. Swihart, S. J. et al. Discovery of a new redback millisecond pulsar candidate: 4FGL J0940.3-7610. Astrophys. J. 909, 185 (2021).

62. Li, K. L., Jane Yap, Y. X., Hui, C. Y. & Kong, A. K. H. Revealing a new black widow binary 4FGL J0336.0+7502. Astrophys. J. 911, 92 (2021).

63. Pletsch, H. J. et al. Binary millisecond pulsar discovery via gamma-ray pulsations. Science 338, 1314-1317 (2012).

64. Ray, P. S. et al. Radio discovery of and gamma-ray pulsations from PSR J2339-0533. Res. Not. Am. Astron. Soc. 4, 37 (2020).

65. Luo, J. et al. PINT: a modern software package for pulsar timing. Astrophys. J. 911, 45 (2021).

66. Harris, C. R. et al. Array programming with NumPy. Nature 585, 357-362 (2020).

67. Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261-272 (2020).

68. Hunter, J. D. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 9, 90-95 (2007).

69. Astropy Collaboration Astropy: a community Python package for astronomy. Astron. Astrophys. 558, A33 (2013).

70. Astropy Collaboration The Astropy Project: building an open-science project and status of the v2.0 core package. Astron. J. 156, 123 (2018).

71. Klöckner, A. et al. PyCUDA and PyOpenCL: a scripting-based approach to GPU run-time code generation. Parallel Comput. 38, 157-174 (2012).

72. Romani, R. W. et al. PSR J0952-0607: the fastest and heaviest known galactic neutron star. Astrophys. J. Lett. 934, L17 (2022).

73. Stringer, J. G. et al. Optical photometry of two transitional millisecond pulsars in the radio pulsar state. Mon. Not. R. Astron. Soc. 507, 2174-2191 (2021).

74. Romani, R. W., Graham, M. L., Filippenko, A. V. & Zheng, W. PSR J1301+0833: a kinematic study of a black-widow pulsar. Astrophys. J. 833, 138 (2016).

75. Romani, R. W. et al. PSR J1810+1744: companion darkening and a precise high neutron star mass. Astrophys. J. Lett. 908, L46 (2021).

76. Romani, R. W. & Shaw, M. S. The orbit and companion of probable γ-ray pulsar J2339-0533. Astrophys. J. Lett. 743, L26 (2011).

## Acknowledgements

C.J.C. thanks B. Allen for useful discussions that led to the use of posterior weights that increased the significances of the detected eclipses. We thank S. Digel, T. Johnson, M. Pesce-Rollins, D. Thompson and Z. Wadiasingh for carefully reviewing the manuscript on behalf of the Fermi-LAT collaboration. C.J.C., R.P.B, M.R.K., D.M.S. and G.V. acknowledge support from the ERC under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 715051; Spiders). This work was supported by the Max-Planck-Gesellschaft (MPG). B.B. acknowledges the support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0700. H.T.C. is a Hubble Fellowship Program Einstein Postdoctoral Fellow. Support for H.T.C. was provided by NASA through the NASA Hubble Fellowship Program grant no. HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract no. NAS5-26555. V.S.D. was supported by the STFC. M.R.K. acknowledges support from the Irish Research Council in the form of a Government of Ireland Postdoctoral Fellowship (GOIPD/2021/670: Invisible Monsters). S.M.R. is a CIFAR Fellow and is supported by the NSF Physics Frontiers Center award no. 1430284.

D.M.S. also acknowledges the Fondo Europeo de Desarrollo Regional (FEDER) and the Canary Islands government for the financial support received in the form of grant no. PROID2020010104. Pulsar research at Jodrell Bank Centre for Astrophysics and access to the Lovell telescope is supported by a consolidated grant from the UK Science and Technology Facilities Council (STFC). Work at the Naval Research Laboratory was supported by the NASA Fermi programme. 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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Etudes Spatiales in France. This work performed in part under DOE contract no. DE-AC02-76SF00515.

## Author contributions

C.J.C. performed the eclipse search analyses and wrote the manuscript. M. Kerr assisted the analysis and exposure calculations. P.B. and L.G. produced the Fermi-LAT datasets and photon probability weights for each pulsar. R.P.B., V.S.D., M.R.K., D.M.S. and G.V. contributed to the interpretation of optical modelling and discussion of gravitational lensing. M. Kerr, B.B., R.P.B, V.S.D., M.R.K., L.N., M.S.E.R. and D.M.S. reviewed the manuscript and provided comments. All remaining authors contributed pulsar timing ephemerides required to phase-fold the Fermi-LAT data.

## Funding

Open access funding provided by Max Planck Society.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-022-01874-x.

**Correspondence and requests for materials** should be addressed to C. J. Clark.

**Peer review information** Nature Astronomy thanks Adam Ingram and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless

indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023

<footnote n="1">¹Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Hannover, Germany. ²Leibniz Universität Hannover, Hannover, Germany.</footnote>

<footnote n="3">³Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Manchester, UK. ⁴Space Science Division, Naval Research Laboratory, Washington, DC, USA. ⁵Max-Planck-Institut für Radioastronomie, Bonn, Germany. ⁶National Centre for Radio Astrophysics, Tata Institute of Fundamental Research, Pune, India. ⁷Laboratoire Leprince-Ringuet, École Polytechnique, CNRS/IN2P3, Palaiseau, France. ⁸South African Radio Astronomy Observatory, Cape Town, South Africa. ⁹Laboratoire de Physique et Chimie de l'Environnement et de l'Espace-Université d'Orléans, CNRS, Orléans, France. ¹⁰Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, Nançay, France. ¹¹Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY, USA. ¹²College of Science, George Mason University, Fairfax, VA, USA. ¹³Department of Physics and Astronomy, University of Sheffield, Sheffield, UK. ¹⁴Instituto de Astrofísica de Canarias, La Laguna, Spain. ¹⁵Department of Physics, University College Cork, Cork, Ireland. ¹⁶Departamento de Astrofísica, Universidad de La Laguna, La Laguna, Spain. ¹⁷University of Virginia, Charlottesville, VA, USA. ¹⁸National Radio Astronomy Observatory, Socorro, NM, USA. ¹⁹Eureka Scientific, Oakland, CA, USA. ²⁰Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, Pessac, France. ²¹ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, Victoria, Australia. ²²Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, Victoria, Australia. ²³New York University Abu Dhabi, Abu Dhabi, United Arab Emirates. ²⁴Department of Physics and Astronomy, West Virginia University, Morgantown, WV, USA. ²⁵Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, Meudon, France. ✉e-mail: colin.clark@aei.mpg.de</footnote>