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

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.

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 1 , a sizable population of compact binary millisecond pulsar systems with semi-degenerate companion stars has emerged 2 . These are often split into two main classes: black widows with companion stars with masses below M c ≲ 0.05 solar masses (M ⊙ ); and redbacks with companion masses 0.1 M ⊙ ≲ M c ≲ 0.5 M ⊙ . 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 3,4 . 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 Article https://doi.org/10.1038/s41550-022-01874-x 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 ≈ 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 × 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 σ γ = 0.03σ T (ref. 19 ), where σ 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 τ γ = σ γ I −1 ΔDM ≈ 6 × 10 −8 I −1 (ΔD M/1 pc cm −3 ), where I is the ionization fraction.
Polzin et al. 4 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 ΔDM ≈ 0.1 pc cm −3 at a phase 0.025 orbits after conjunction, corresponding to τ γ ≈ 6 × 10 −9 I −1 . Assuming an isotropic wind, with density ρ decreasing with radius r from the companion star as ρ(r) ∝ 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. 20 ).
Of course, this model for the companion wind is overly simplistic: spider companions have non-isotropic swept-back winds 11 that often vary with time 13 , 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 ℒ 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 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. [5][6][7] and there are hints that spider pulsars may be systematically heavier than other species of binary neutron star 8,9 . Several classes of theoretical neutron star equation-of-state (see ref. 10 and references therein) models predict maximum masses close to 2 M ⊙ , 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 ∝ 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 11 . In several redbacks, this temperature pattern is even seen to vary over time (for example refs. [12][13][14]. 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) 15 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 16 , 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. 17 found evidence in the Fermi-LAT data, later confirmed by Kennedy et al. 18 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 × 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 Article https://doi.org/10.1038/s41550-022-01874-x 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 ⊙ . 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 ⊙ .
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° ≲ i ≲ 67° (refs. 23,24 ). When combined with optical spectroscopy results 5 , this corresponds to an extremely high mass of M psr = 2.4 ± 0.1 M ⊙ , 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 10 . Our detection of a gamma-ray eclipse, however, requires a much higher i, >84.1°. This lower bound on the inclination corresponds to M psr = 1.81 ± 0.07 M ⊙ , 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   The dotted and solid black curves show the empirical survival function (that is, the fraction of pulsars that survive a given log ℒ 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.
Article https://doi.org/10.1038/s41550-022-01874-x Guillemot et al. 28 and Johnson et al. 26 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. 23 and Draghis et al. 24 , 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 Contour lines are drawn at log ℒ corresponding to false-alarm probabilities of 5% (yellow) and 1% (red). These levels are also marked on the colour bar, along with the log ℒ corresponding to a false-alarm probability of 0.1% (dark red), although this level is never reached. The position of the maximum likelihood ( ℒ max ) is marked by a cyan cross, and the corresponding log ℒ 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 ℒ 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 ℒ max found were greater than 40%.

Table 1 | Constraints for pulsars with detected eclipses
Pulsar Class θ 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.
Article https://doi.org/10.1038/s41550-022-01874-x 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 11 or diffusion on the stellar surface causing heat to 'leak' over the terminator 29 -a possibility that Reynolds et al. 23 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 11 , 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 30 , 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 31 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 7 and Kandel and Romani 32 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][34][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. Parameters are the same as in Table 1, but inclination upper limits i max and mass lower limits M min psr and M min c 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.

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 37,38 . Sensitive unbinned-likelihood-based methods for detecting eclipses (for example ref. 39 ) 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 40 . 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 41 . For this, we used the 10 yr incremental version (Data Release 2) of the Fermi-LAT Fourth Source Catalog (4FGL; 42,43 ) (https://fermi.gsfc. nasa.gov/ssc/data/access/lat/10yr_catalog/) 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. 41 . 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 44 . For each pulsar, we computed the orbital phase at which each photon was emitted according to these ephemerides using the TEMPO2 software 45 . 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. 18 and ref. 39 . 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 θ of the orbital period. Within the eclipse, we assumed a constant flux level, which we parameterized with α, the fractional flux level within the eclipse relative to the overall average flux. The increase in log ℒ for such an eclipse, compared with the null hypothesis of photons being uniformly distributed in orbital phase, is: where w i is the photon probability weight for the ith photon; Θ (Θ ) refers to photons with orbital phases inside (outside) the eclipse; and η Θ (ηΘ ) 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 39 , and folded these on the orbital period to compute η Θ and ηΘ . 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 α = 0, testing for θ ∈ [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 psr /M c = 3.93 ± 0.06, which would eclipse for 8.4% of an orbit if the companion filled its Roche lobe and was observed at i = 90°.
For pulsars in which significant eclipses were detected we also tested alternative eclipse models with 0 < α < 1 and curved (rather than sharp) ingresses and egresses. No significant log ℒ 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 39 , 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 ϕ at which the photon was emitted by applying Bayes' theorem:
Here P(S|ϕ) is now the posterior probability of the photon having been emitted by the pulsar, given its rotational phase; p(ϕ|S) is the phase distribution of photons emitted by the pulsar (that is, the pulsar's pulse profile, which we hereafter denote as f(ϕ)); and p(ϕ|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(ϕ|B) = 1. Rewriting equation (2) with these values gives us the reweighting equation: Hereafter w is the prior weights, and w′ is the posterior weights. For phases within peaks of the pulse profile, where f(ϕ) > 1, these posterior weights are always greater than the prior weights, and for phases outside of peaks, where f(ϕ) < 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 ℒ.
To obtain the pulse profile models, f(ϕ), we fitted a set of wrapped Gaussian functions to the prior-weighted photon phases using the Article https://doi.org/10.1038/s41550-022-01874-x maximum-likelihood method described by Abdo et al. 44 . The number of Gaussian functions used to model each pulse profile was chosen to minimize the Bayesian information criterion 46 .
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 ℒ values, with δ log ℒ increasing by at least 2.6 for these pulsars. Only the eclipse from PSR J1555−2908 was undetected using the prior weights, with δ log ℒ = 3.35 compared with δ log ℒ = 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 eclipse 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 20 . 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. 17,18 ), 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 θ 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 θ 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 ≈ 5 × 10 −5 .
The δ log ℒ 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 δ log ℒ 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 ≈ 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 δ log ℒ 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. 47 ) or from the amplitude of the 'ellipsoidal' component of an observed optical light curve 48 . 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 48 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 49 , 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, 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 Article https://doi.org/10.1038/s41550-022-01874-x 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° 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. 47 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. 24 ), 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° (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 13,50 . 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. 50 fixed the inclination to i = 76°, 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. 50 found that f RL ≈ 0.85 was compatible with multiple light curves despite long-term variability. From optical spectroscopy, Miraval Zanon et al. 51 found an observed companion radial velocity amplitude of 343.3 ± 4.4 km s −1 . Using these values, we found that the observed eclipse duration required an inclination greater than 80. 9° (compare with 80.4° assuming f RL = 1 in Table 1). Heating corrections reduced the estimated centre-of-mass velocity to 298.7 ± 7.7 km s −1 , for a larger mass ratio (and hence lower minimum inclination of i = 80.1°) but a lower pulsar mass, M psr ≈ 1.1 M ⊙ . 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. 52 in a targeted radio search of a steep-spectrum radio continuum source identified within a pulsar-like gamma-ray source by Frail et al. 53 . Modelling of both optical photometry and spectroscopy by Kennedy et al. 30 , using a model that includes the possibility of heat diffusion across the terminator, revealed the companion's projected radial velocity to be 397 ± 2 km s −1 , and indicated a high binary inclination of i > 75°, giving a maximum pulsar mass of 1.82 M ⊙ . 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° < i < 86.2°, for a pulsar mass 1.58 M ⊙ < M psr < 1.71 M ⊙ 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. 47 . 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. 47 estimated a radius that corresponds to f RL ≈ 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° and a lower pulsar mass range of 1.68 M ⊙ < M psr < 2.11 M ⊙ For PSR J2129−0429, Bellm et al. 54 measured a projected companion radial velocity amplitude of K 2 = 250 ± 4 km s −1 (for mass ratio q = 3.93 ± 0.06), and inferred a filling factor of f RL = 0.82 ± 0.03 and an inclination i > 68° 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° < i < 78.3°, consistent with the range allowed by optical modelling. This corresponds to a pulsar mass range of 1.61 M ⊙ < M psr < 1.88 M ⊙ at 95% confidence. Al Noori et al. 55 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 21,[56][57][58][59][60][61][62] . 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 14,22,63,64 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 Article https://doi.org/10.1038/s41550-022-01874-x 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 δ log ℒ, 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 δ log ℒ ≈ 9, as opposed to δ log ℒ ≈ 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 δ log ℒ = 10.8 and δ log ℒ = 11.3, corresponding to trial-corrected P FA of 3 × 10 −3 and 1 × 10 −3 , respectively. We show the δ log ℒ 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 δ log ℒ values (δ log ℒ = 11.8 and δ log ℒ = 14.9, respectively) at orbital phases that are ΔT asc ≈ 300 s and ΔT asc ≈ 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. 60 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 ≤ 90°, as systems at i are indistinguishable from those at 180° − 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 ≡ M psr /M 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 redbacks 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 redbacks. 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 45 , used to compute rotational and orbital photon phases, is available via Bitbucket at https://bitbucket.org/psrsoft/tempo2. Icarus 48 , used to commpute expected eclipse durations to derive pulsar mass estimates, is available via GitHub at https://github.com/bretonr/ Icarus. Godot 39 , used to compute Fermi-LAT exposure, is available via GitHub at https://github.com/kerrm/godot. PINT 65  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 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