Open Access
Issue
A&A
Volume 693, January 2025
Article Number A314
Number of page(s) 17
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202452035
Published online 29 January 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The formation of radio millisecond pulsars (MSPs) has long been connected to low-mass X-ray binaries (LMXBs; Alpar et al. 1982; Radhakrishnan & Srinivasan 1982; Bhattacharya & van den Heuvel 1991). Observed MSPs have short spin periods (≲30 ms) and low surface magnetic fields (∼108 G). The neutron stars (NSs) in LMXBs accrete matter that spins them up to millisecond periods (Alpar et al. 1982). This phase is observed as an X-ray phase and is also called the recycling phase due to the pulsar spin-up. Observations of accreting X-ray MSPs confirmed the link between MSPs and accreting LMXBs (Wijnands & van der Klis 1998; Archibald et al. 2009). The NS surface magnetic field has also been suggested to decrease with the accreted mass due to Ohmic decay, explaining the observed magnetic field strengths of MSPs (Konar & Bhattacharya 1997; Zhang 1998; Konar & Bhattacharya 1999a,b; Cumming et al. 2001). Compact binary MSPs are a sub-category of MSP binaries with compact orbits and eclipses in their radio signals. In compact binary MSPs, after the mass-transfer phase or Roche-lobe overflow (RLO) has ended, the highly spinning pulsar loses energy in the form of energetic particles that irradiate the outer surface of the companion, leading to mass loss (Archibald et al. 2009). Initial observations of eclipsing MSPs showed that some of them have very low mass companions (0.02 to 0.05 M; Fruchter et al. 1988; Stappers et al. 1996a, 1999), which could be explained by mass loss due to pulsar wind irradiation. Evidence of irradiation-induced mass loss has since been observed in multiple eclipsing compact MSPs (Kluzniak et al. 1988; Ruderman et al. 1989a,b; Stappers et al. 1996b; Polzin et al. 2020; Kumari et al. 2024).

Due to the cannibalistic nature of this interaction, these compact binary MSPs, with observed orbital periods Porb ≲ 1 d (Pallanca et al. 2012; Romani et al. 2012; Kaplan et al. 2013; Breton et al. 2013; Roberts 2013), are also called spider pulsars. So far, about 60 spiders have been observed and studied in the Galaxy (for instance, see Abdo et al. 2009; Pletsch et al. 2012; Pallanca et al. 2012; Romani et al. 2012; Kaplan et al. 2013; Breton et al. 2013). There are two main categories of spiders (Roberts 2013) based on the masses of the companions to the pulsar (denoted as Mc): redbacks (RBs), which have companions with masses 0.1 M ≲ Mc ≲ 0.5 M, and black widows (BWs), whose companions have masses of Mc < 0.1 M. The orbital period distribution is similar for both types (≲1 d).

Following classical binary evolution, the observed orbital periods of all spiders have not been explained, as they fall short by about an order of magnitude (Podsiadlowski et al. 2002). Fast winds that leave a binary isotropically and without interacting with anything else have a widening effect on the binary orbit (Soberman et al. 1997). Therefore, the involvement of this energetic radio pulsar wind has been invoked to explain the orbital properties of recycled pulsars (Ergma & Fedorova 1992; Podsiadlowski et al. 2002; Benvenuto et al. 2014, 2015; Chen et al. 2013). There are a few ideas that have been proposed to explain spider MSPs. Chen et al. (2013) proposed strong and weak irradiation producing RBs and BWs, respectively. Another group studied the formation of spiders by including X-ray irradiation feedback, and they found that in order to reproduce the orbital parameters of spiders, cyclic mass transfer due to the effect of X-rays is required (Benvenuto et al. 2014, 2015). Ginzburg & Quataert (2020) and Ginzburg & Quataert (2021) argue that sustained magnetic braking (which couples to stellar wind and leads to angular momentum loss) is needed to reproduce spiders and that ablation alone is not sufficient. Some spiders are also thought to come from ultra-compact X-ray binaries (UCXBs), which are a sub-category of LMXBs with hydrogen-deficient companions and orbital periods of less than an hour (see review by Nelemans & Jonker 2010).

Since spider binaries have passed through a long mass-transfer phase, it has been suggested that they host some of the most massive (Antoniadis et al. 2016; Linares et al. 2018; Strader et al. 2019; Linares 2020) and fastest spinning NSs. The fastest spinning pulsar known so far spins at a frequency of 716 Hz (PSR J1748–2446ad spinning at about 1.4 ms; Hessels et al. 2006) in the globular cluster Terzan 5. Various reasons have been discussed to explain the lack of pulsars spinning faster than 1 ms, for example an additional torque (from gravitational wave emission) that acts at high frequencies and spins down the pulsar (Papaloizou & Pringle 1978; Wagoner 1984; Bildsten 1998; Bhattacharyya & Chakrabarty 2017). Since pulsar spins depend on the accretion phase, in this work we investigate the fastest spins possible for pulsars. Despite relatively large mass uncertainties, NSs in spider binaries tend toward higher masses (median MNS ∼ 1.8 M; Linares 2020) in comparison to other systems, such as double NSs or wide NS plus white dwarf (WD) binaries (Demorest et al. 2010; Antoniadis et al. 2013; Linares 2020).

In this study, we carry out detailed binary evolution simulations, calculating the evolution of the pulsar spin and taking into account isolated and mass-transfer evolution. The main aim is to investigate the effect of various assumptions of binary interaction on the pulsar spin evolution and to constrain them when comparing with observed spider properties. In Section 2, we describe the initial parameters of the simulated tracks, the added physics of pulsar spin evolution, and the treatment of the pulsar wind irradiation. Section 3 presents our simulated grids in the context of the observed population of spiders. In Section 4, we compare our results with those in the literature and discuss the caveats involved. This is followed by the conclusions of our study in Section 5.

2. Methods

2.1. Initial pulsar properties and binary physics

We use the detailed stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019) to carry out our calculations for the evolution of spiders and included pulsar spin evolution. One of the benefits of using a detailed stellar code like MESA is that compared to rapid population synthesis codes, detailed codes provide a more self-consistent calculation of the mass-transfer phase along with a physically motivated estimation of the stability of binary interaction. All the initial binaries in this study have a low- to intermediate-mass main-sequence star (MS; range of 1.0 to 4.5 M) and a point mass NS (range of 1.3 to 2.0 M). The initial orbital periods range from 0.6 to 4 d. The simulations are run at Solar metallicity (0.0142; Asplund et al. 2009). Since most of the observed spider binaries with dynamical mass measurements are within the Milky Way (Linares 2020), we are focusing on Galactic spiders for our study. The physics in the MESA grids follows the description in Fragos et al. (2023). We treat our numerically failed runs similar to Fragos et al. (2023), by identifying failed binaries and rerunning them with changed parameters that will help computational stability and reach an acceptable termination condition. We discuss our treatment of the numerically failed binaries in Appendix A.

2.2. Isolated pulsar phase

As initial pulsar spin parameters, we take a spin of 1 s and a spin period derivative of 10−15 s/s. The corresponding magnetic field strength is 1012 G. These values correspond to the average spin parameters of observed young pulsars. We model the pulsar spin evolution by identifying different stages of the evolution and applying the corresponding equations, as described in the following sections. Pulsars in detached binaries, where the binary components are orbiting each other without transferring mass, are essentially assumed to spin down as isolated pulsars as their spins and magnetic fields are not affected by the companion. Since pulsars are rapidly spinning dipoles, they lose energy as they spin down in the form of magnetic dipole radiation, provided the dipole has an inclination. Hence, the spin-down of a pulsar results in the following decrease in its angular frequency ( Ω ˙ $ \dot{\Omega} $):

Ω ˙ = 2 μ 2 sin 2 α 3 μ 0 c 2 I Ω 3 , $$ \begin{aligned} \dot{\Omega } = -\frac{2\mu ^{2}\sin ^{2}{\alpha }}{3\mu _{0}c^{2}I}\Omega ^{3}, \end{aligned} $$(1)

where μ is the pulsar magnetic moment (μ ≈ BR3), α is the angle between the rotation axis and magnetic axis, μ0 is the permeability of free space (it is unity in CGS units), I is the pulsar moment of inertia (approximated as 2 M NS R NS 2 / 5 $ 2M_{\mathrm{NS}}R_{\mathrm{NS}}^{2}/5 $ for a solid sphere, where MNS and RNS are the NS mass and radius), and c is the speed of light. We assumed α = 45°, which is a simplistic approximation, as the inclination would decay with time (Spitkovsky 2006; Johnston & Karastergiou 2017). The full range of α values would give an order of magnitude difference in Ω ˙ $ \dot{\Omega} $, at most (Spitkovsky 2006); hence, we did not vary this parameter in our study and fixed it at an intermediate value.

We evolved the isolated pulsar surface magnetic field (B) decay following Chattopadhyay et al. (2020):

B = ( B i B min ) × exp ( t / τ d ) + B min , $$ \begin{aligned} B = (B^{i} - B_{\rm min})\times \exp {(-t/\tau _{\rm d})} + B_{\rm min}, \end{aligned} $$(2)

where Bi is the initial surface magnetic field, Bmin is the minimum magnetic field (we take Bmin = 108 G; Zhang & Kojima 2006; Osłowski et al. 2011), t is time, and τd is the magnetic field decay timescale (we take τd = 3 Gyr; Ye et al. 2019). The superscripts i and f stand for initial and final values. Previous studies have explored τd values ranging from a few Myr to a few Gyr, with decay timescales ≳1 Gyr being favored (Kiel et al. 2008; Ye et al. 2019; Chattopadhyay et al. 2020).

2.3. Interacting binary phase

Once the companion to the pulsar fills its Roche lobe and the mass-transfer phase begins, accreted matter will carry angular momentum and spin up the pulsar. We consider mass accretion when there is a persistent accretion disk and use the criteria by Dubus et al. (1999) that describes the limiting mass-transfer rate below which thermal-viscous instabilities would prevent accretion. This limit is as follows:

M ˙ crit = 3.2 × 10 9 ( M NS 1.4 M ) 1 / 2 ( M c 1 M ) 1 / 5 ( P orb 1 d ) 7 / 5 M yr 1 . $$ \begin{aligned} \dot{M}_{\rm crit} = 3.2\times 10^{-9} \bigg (\frac{M_{\rm NS}}{1.4\,\mathrm{M}_\odot }\bigg )^{1/2} \bigg (\frac{M_{\rm c}}{1\,\mathrm{M}_\odot }\bigg )^{-1/5} \bigg (\frac{P_{\rm orb}}{{1\,\mathrm d}} \bigg )^{7/5} {\mathrm{M_\odot \,yr^{-1}}}. \end{aligned} $$(3)

The change in the NS mass (NS) is determined by the efficiency of accretion (ϵ), which is assumed to be equal to 1 − β, where β is the fraction of mass transferred by the companion that is lost from the vicinity of the NS as a fast isotropic wind, taking away angular momentum from the NS. The rest of the material (1 − β) is accreted by the NS, until the mass-transfer rate reaches the Eddington limit, beyond which all accretion is limited. Since the term β affects the amount of accretion, lower values of β result in a faster spin-up of the pulsar compared to higher values. Additionally, the value of β affects the orbital evolution following the equation for the relative change in orbital separation (ȧ/a; Tauris & van den Heuvel 2006), assuming other changes in orbital angular momentum are represented by J ˙ orb $ \dot{J}_{\mathrm{orb}} $:

a ˙ a = J ˙ orb J orb 2 ( 1 + ( β 1 ) M c M NS β M c 2 ( M c + M NS ) ) M ˙ c M c , $$ \begin{aligned} \frac{\dot{a}}{a} = \frac{\dot{J}_{\rm orb}}{J_{\rm orb}} -2\,\bigg (1 + (\beta - 1) \frac{M_{\rm c}}{M_{\rm NS}} - \frac{\beta M_{\rm c}}{2 (M_{\rm c} + M_{\rm NS})}\bigg ) \frac{\dot{M}_{\rm c}}{M_{\rm c}}, \end{aligned} $$(4)

where, Mc and MNS are companion and NS masses, respectively. Other processes to determine J ˙ orb $ \dot{J}_{\mathrm{orb}} $ (apart from mass transfer/loss) are magnetic braking (Rappaport et al. 1983), gravitational wave radiation (Peters 1964), and spin-orbit coupling.

Theoretical and observational studies of various types of binaries have shown that the parameter β could be more than 0.5 in LMXBs despite the evolution being mostly sub-Eddington (Tauris & Savonije 1999; Ritter & King 2001; Podsiadlowski et al. 2002; Jacoby et al. 2005; Zhang et al. 2011; Antoniadis et al. 2012; Tauris et al. 2012). Some possible reasons for this are the propeller effect (Illarionov & Sunyaev 1975) and thermal-viscous instabilities in the accretion disk (Pringle 1981; Mineshige 1993; van Paradijs 1996); all effects that would prevent accretion. Studies of accreting pulsars use values in the range of 0.5 to 0.7 (Podsiadlowski et al. 2002; Chen et al. 2013; Jia & Li 2014, 2015; Misra et al. 2020; Wang et al. 2021). We add a prescription for thermal-viscous accretion disk instabilities (explained further below in Section 2.3) and explore three values of β (0.0, 0.3, and 0.7) to study its effect on the pulsar spin evolution. This fractional value β is also a proxy for the change in angular momentum when mass is accreted. We include β even though we have a separate treatment for the propeller phase (discussed below in Section 2.3.2) and the aforementioned disk instabilities, since these have large uncertainties and there might still be inefficient transfer of angular momentum due to some unforeseen physical process.

The accreted matter suppresses the magnetic field due to rapid Ohmic decay (Konar & Bhattacharya 1997, 1999a,b), which we include as an exponential decay with accreted mass (Osłowski et al. 2011; Chattopadhyay et al. 2020), as follows:

B = ( B i B min ) × exp ( Δ M NS / M d ) + B min . $$ \begin{aligned} B = (B^{i} - B_{\rm min})\times \exp {(-\Delta M_{\rm NS}/M_{\rm d})} + B_{\rm min}. \end{aligned} $$(5)

In the above equation, ΔMNS is the accreted matter and Md is the magnetic field mass decay scale (0.025, although other values for Md have been explored in the range of 0.0033 to 0.05; Osłowski et al. 2011; Chattopadhyay et al. 2020). The accreted matter also causes a change in the angular momentum of the NS, which can be estimated as

J ˙ = κ M ˙ NS R NS 2 ω diff , $$ \begin{aligned} \dot{J} = \kappa \dot{M}_{\rm NS} R^{2}_{\rm NS} \omega _{\rm diff}, \end{aligned} $$(6)

where κ is an efficiency factor and ωdiff is the angular velocity of the incoming material (Kiel et al. 2008; Chattopadhyay et al. 2020) that is defined by

ω diff = G M NS r mag 3 G M NS r cor 3 · $$ \begin{aligned} \omega _{\rm diff} = \sqrt{\frac{GM_{\rm NS}}{r^{3}_{\rm mag}}} - \sqrt{\frac{GM_{\rm NS}}{r^{3}_{\rm cor}}}\cdot \end{aligned} $$(7)

In the above equation, rcor is the co-rotation radius, which is defined as rcor = (GMNS2)1/3 and rmag is the magnetospheric radius, which is defined as follows (Kiel et al. 2008; Chattopadhyay et al. 2020):

r mag = r Alfv e ´ n 2 = ( 2 π 2 G μ 0 2 ) 1 / 7 ( R NS 6 M ˙ NS M NS 1 / 2 ) 2 / 7 B 4 / 7 , $$ \begin{aligned} r_{\rm mag} = \frac{r_{\rm Alfv\acute{e}n}}{2} = \Bigg (\frac{2\pi ^{2}}{G\mu ^{2}_{0}}\Bigg )^{1/7} \Bigg (\frac{R^{6}_{\rm NS}}{\dot{M}_{\rm NS}M^{1/2}_{\rm NS}}\Bigg )^{2/7} B^{4/7}, \end{aligned} $$(8)

where, rAlfvén is the Alfvén radius and G is the gravitational constant. The two radii, rcor and rmag, together with the light-cylinder radius, defined as rlc = c/Ω, are associated with the NS and used to study the interacting binary phase. Using the three characteristic radii (rmag, rcor, and rlc), we define below the three main regimes of accretion:

2.3.1. Accretion regime, rmag < (rcor, rlc)

If the mass-transfer rate is high enough that the magnetospheric radius is below the co-rotation radius (following Equation 8), we assume that accretion can proceed unimpeded. Since the magnetic field is weak, we assume that the inner accretion disk extends down to the NS surface. Hence, in Equation (6), RNS is equal to the NS radius, which we take as 12.5 km (for example, Raaijmakers et al. 2021). Since, the increase in angular momentum comes from accreted material, we take NS = ϵ = (1−β)c, ϵ being the accretion efficiency. We also assume a κ of unity, similar to Chattopadhyay et al. (2020).

2.3.2. Propeller regime, rcor ≤ rmag < rlc

While no mass is thought to have been accreted during the propeller phase (Illarionov & Sunyaev 1975), it is uncertain how much material is propelled away. PSR J1023+0038 appears to be undergoing spin-down in a weak-propeller state, while still showing pulsed X-rays that signify accretion (Papitto et al. 2015; Archibald et al. 2015; Ertan 2018). Following Equation (7), the change in angular momentum during the propeller phase is negative and the result is a spin down of the pulsar. We assume 10% of mass is accreted out of the mass transferred by the companion (β = 0.9), while the rest is driven away from the NS, slowing it down. To account for spin down due to matter propelled away from the pulsar, we take NS = βṀc in Equation (6). Additionally, we assume a κ of 0.001 in the Equation (6), due to the uncertain nature of this phase.

2.3.3. Ejection regime and pulsar wind phase

If the magnetic field is strong enough or the mass-accretion rate is low enough, the magnetosphere extends far outside of the light cylinder radius. Any accretion or transfer of angular momentum is prevented and the pulsar spin (and the corresponding pulsar magnetic field) is evolved as that of an isolated pulsar (see Section 2.2). We use two criteria to determine when the ejection regime and the consequent pulsar wind phase begins. First, the pulsar has accreted enough matter to spin up to MSP-like spins (we use 30 ms as the identifying upper limit). The second condition for the onset of the pulsar wind phase is that the mass-transfer rate has decreased below the criteria for thermal-viscous instabilities and can no longer sustain a persistent accretion disk (see Equation 3), under the effect of X-ray irradiation during the RLO phase.

Once the ejection phase begins and as the pulsar spins down, the pulsar emits an energetic wind with a certain luminosity, known as the spin-down luminosity, which is defined using the spin properties of the pulsar as follows:

L pulsar = 4 π 2 I P ˙ spin P spin 3 · $$ \begin{aligned} L_{\rm pulsar} = \frac{4\pi ^{2} I \dot{P}_{\rm spin}}{{P}_{\rm spin}^{3}}\cdot \end{aligned} $$(9)

As this energetic wind hits the companion, it ablates the outer surface of the star, which, when observed, signifies the presence of a spider pulsar (Kluzniak et al. 1988; Ruderman et al. 1989b,a; Stappers et al. 1996b). In the above equation, the initial values of the NS spin (Pspin) and change in NS spin with time (spin) are determined by the preceding accretion phase (see Section 2.3.1) and follow the spin evolution of a detached pulsar (see Section 2.2). The irradiation from the pulsar wind drives a mass loss from the companion surface (c,irr) defined by Stevens et al. (1992):

M ˙ c , irr = f pulsar × L pulsar 2 v esc 2 ( R c a ) 2 , $$ \begin{aligned} \dot{M}_{\rm c, irr} = -f_{\rm pulsar}\times \frac{L_{\rm pulsar}}{2v^{2}_{\rm esc}}\bigg (\frac{R_{\rm c}}{a}\bigg )^{2}, \end{aligned} $$(10)

where fpulsar is the efficiency that describes the transfer of energy from the incoming spin-down luminosity to the generated wind loss. It also stands for any geometrical beaming effects that might be present during the process. The other parameters in the above equation include Lpulsar as the spin-down luminosity, vesc as the escape velocity from the companion star surface, Rc as the companion radius, and a as the binary separation. The effect of this increased mass loss from the companion is described by the fast mode or Jean’s mode of mass transfer (Jeans 1924; Soberman et al. 1997). This mode of mass transfer assumes that a small fraction of the stellar mass is lost in the form of fast, spherically symmetric winds that leave the binary without interacting with any binary components, escaping with the specific angular momentum of the mass-losing star. The orbital angular momentum (Jorb) is defined as

J orb = μ G a ( M c + M NS ) , $$ \begin{aligned} J_{\rm orb} = \mu \sqrt{Ga(M_{\rm c}+M_{\rm NS})}, \end{aligned} $$(11)

where μ is the reduced mass of the binary. The change in orbital angular momentum per unit of reduced mass (J/μ) is negligible and the orbit expands accommodate the change in total mass (Huang 1963). As the effect of irradiation will not necessarily be spherically symmetric, we have the efficiency term fpulsar as described above (Stevens et al. 1992; Chen et al. 2013).

The definition of the three regimes of mass transfer described above based on the relative radii of the magnetosphere, co-rotation radius, and light cylinder radius, are only approximate estimations. There is work showing that the boundaries between these regimes are bit more complicated (Ertan 2015, 2018). However, since they are highly uncertain, we assume the more straightforward definitions described.

2.4. Observations used for comparisons to simulations

To compare our simulations to observations we use a catalog of BWs, RBs, and huntsmen (including candidates). A total of 57 spiders are used, 52 of which are cataloged in Linares & Kachelrieß (2021, and references therein). The more recent ones are from various other works in the literature, and an updated catalog has been prepared by Nedreaas (2024, and references therein). The reported values for the companion masses are in terms of the minimum mass (denoted as Mc, min), calculated assuming a pulsar mass of 1.4 M and inclination of 90° from observations of radio pulses. Only for one of the huntsmen source in the catalog have radio pulses been confirmed and Mc, min estimated (1FGL J1417.7–440, Camilo et al. 2016). While the other huntsmen, 2FGL J0846.0+2820, only has Fermi Large Area Telescope (LAT) observations (Swihart et al. 2017).

3. Results

3.1. Mass-accretion efficiency

We investigate the formation of different branches of compact MSP binaries using detailed evolutionary sequences. Figures 1 and B.1 show the evolutionary tracks with β of 0.3 and 0.7, respectively, each with different values of fpulsar. The tracks correspond to initial NS masses of 1.3 M and have final NS spins ≲30 ms. While β = 0.7 reproduces the very compact spiders, it cannot explain a large part of the binary parameter space where these systems are observed. This behavior is irrespective of the pulsar wind efficiencies. Similar behavior is seen for huntsmen spiders (Porb > 1.0 d), as β = 0.7 does not reproduce the small population of observed huntsmen spiders. With β of 0.3, a large part of the observed parameter space is explained and on increasing fpulsar, more and more BWs and RBs’ observed orbital properties are reproduced. The reason for the difference between β = 0.3 and β = 0.7 can be understood by looking at Equation (10), which shows that the mass lost due to pulsar wind irradiation is directly proportional to the pulsar spin-down luminosity and in turn depends on the pulsar spin at the end of the RLO phase (Equation (9)). Since the preceding accretion phase dictates the spin-up of the pulsar, we can compare the resulting pulsar spins in our grids with observations to see why more accretion helps.

thumbnail Fig. 1.

Orbital period evolution versus companion mass of the simulated binaries for β = 0.3 with an initial NS of mass 1.3 M and fpulsar values of 0.0 (top left), 0.05 (top right), 0.1 (bottom left), and 0.5 (bottom right). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. The color bar shows the changing surface abundance of the companion. We also compare our tracks to estimated Mc, min and observed Porb for BWs (black crosses) and RBs (red crosses), and huntsmen spiders (green crosses), including whether H-lines have been detected (filled circles) or not (filled triangles) or if it is unknown (crosses), from Nedreaas (2024, and references therein). The green star is the huntsman candidate 2FGL J0846.0+2820, observed by Fermi-LAT (Swihart et al. 2017) (surface H is unknown). The three tidarrens (J1311–3430, J1653–0159, and J0636+5128) are also shown. The initial binaries have a MS star (range of 1.0 to 4.5 M) and initial periods of 0.6 to 4.0 d.

We compared the observed BWs and RBs to the pulsar spin evolution from the grid in Figure 2, corresponding to fpulsar = 0.0 and β = 0.0, 0.3, and 0.7. The figure shows the binaries for two initial NS masses 1.3 M and 2.0 M, and also separates out the cases where the final NS mass exceeds 2.5 M, which is an approximate upper limit for the NS mass. The NSs with final masses above 2.5 M are considered to have collapsed into black holes (BHs) and would no longer be observed as radio pulsars. The value β = 0.0 leads to a majority of the NSs with M NS i > 1.3 M $ M^{i}_{\mathrm{NS}} > 1.3\,\mathrm{M}_\odot $ gaining spins in the sub-MSP range (≲1 ms) and having final NS masses ≳2.5 M. With a β of 0.7 (or an accretion efficiency of 0.3), the fastest spinning Galactic BW observed so far can be reproduced (PSR J0952–0607, with a spin period of 1.41 ms; Bassa et al. 2017). However, these short spins correspond to the case where the final NS mass exceeds the maximum NS mass limit. Removing these, the minimum spin period for β = 0.7 is 3.8 ms, which is almost thrice the minimum observed BW spin. Looking at the case of β = 0.3 in Figure 2, most of the observed spins can be reproduced. Comparing the minimum final spin period for β = 0.3 (1.38 ms) and the corresponding final pulsar mass (2.36 M) with PSR J0952–0607, we find that our simulations can reproduce well this extreme BW pulsar (its estimated pulsar mass is 2 . 35 0.17 + 0.17 M $ 2.35^{+0.17}_{-0.17}\,\mathrm{M}_\odot $; Romani et al. 2022).

thumbnail Fig. 2.

Evolution of the NS spins (Pspin) versus the orbital periods (Porb) in converging binaries for β = 0.0 (top), β = 0.3 (middle), and β = 0.7 (bottom). We distinguish the binaries if the final NS mass was above 2.5 M (faint gray curves) or below it, which are further grouped by two initial NS masses ( M NS i $ M^{i}_{\mathrm{NS}} $) 1.3 M (blue curves) and 2.0 M (orange curves). We also compare to the observed spins of BWs and RBs, shown by the black and red circles, respectively.

In Figure 2, there are some binary tracks that ended with final spins below 1 ms, in the sub-MSP regime (all with M NS i > 1.3 M $ M^{i}_{\mathrm{NS}} > 1.3\,\mathrm{M}_\odot $). The NSs with M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $ gained a maximum of 1.06 M and did not spin up to the sub-MSP regime (final spin was 1.37 ms). All the binaries where the final NS mass exceeded the maximum NS mass limit, were spun up to sub-milliseconds. Hence, the mass required to spin up an NS to sub-millisecond range is high enough to collapse it into a BH. Additionally, the sub-millisecond regime is mainly dominated by initially more massive NSs (majority of the M NS i = 2.0 M $ M^{i}_{\mathrm{NS}} = 2.0\,\mathrm{M}_\odot $ tracks in Figure 2 appear to be above the maximum NS mass limit) while the initially less massive NSs are not always able to accrete enough matter to enter the sub-millisecond regime. This is because more massive NSs leads to more a stable and longer lasting mass transfer phase (Misra et al. 2020). Comparing the final spin distributions between β = 0.3 and β = 0.7, we can see that β = 0.7 results in slower spinning NS, which would then have lower spin-down luminosities as L pulsar P spin 3 $ L_{\mathrm{pulsar}} \propto P^{-3}_{\mathrm{spin}} $, and not affecting the orbital evolution as much. Increased wind mass loss from the companion would lead to increased orbital widening during the radio pulsar phase. With β = 0.0, while the entire observed spider parameter space is reproduced, the spin evolution extends down to 1 ms. The value of β that is most consistent with the observed minimum pulsar spin is 0.3. Additionally, β = 0.3 leaves some room for uncertainty in accretion and in our implementation of the X-ray irradiated accretion disk instabilities. With future observations of spiders, any spins approaching 1 ms would signify more conservative nature of RLO at play. Hence, β ≲ 0.3 is required to explain the observed spins of spider pulsars.

3.2. Orbital evolution and surface hydrogen abundance

The top left panel in Figure 1 shows the evolutionary tracks with no pulsar wind irradiation (fpulsar = 0.0). The color bar corresponds to the surface H abundance of the companion, which can be compared to optical spectroscopic observations of spiders. We mark in Figure 1, spiders where H lines have not been detected in their optical spectra (filled triangles) and those where H lines are detected (filled circles). Two of the most compact BWs, PSR J1311–3430 (Pletsch et al. 2012; Romani et al. 2012) and PSR J1653–0159 (Kong et al. 2014; Romani et al. 2014), with negligible surface H abundance, are easily reproduced by our evolutionary tracks. These sources also form a sub-category of BWs, called tidarrens, with Porb ≲ 2 h and Mc ≲ 0.015 M (Romani et al. 2016). They are reproduced regardless of the value of fpulsar. Their evolutionary tracks also pass through the ultra-compact X-ray binary region (UCXB), with orbital periods less than one hour, with the minimum orbital period being 30.8 min. Most UCXBs are also considered to be H-deficient due to their small orbits (Nelson et al. 1986; Nelemans et al. 2006), similar to our binary evolution tracks for tidarrens. Hence, tidarrens (the most compact BWs with the lightest companions) passed through a UCXB phase in their previous evolution. Additionally, we also see these systems reproduced in Figure B.1 (with β = 0.7), since not a lot of accretion is needed to form these sources.

We can understand their evolution better in Figure 3 (top left panel) where we show tracks that pass through the tidarren region (solid gray line, with an M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $ and P orb i = 2.4 $ P^{i}_{\mathrm{orb}} = 2.4 $ d) and compare these with other tracks with initially narrower and wider orbits. All the tracks shown in the figure have an M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $. We also show a binary that goes well into the UCXB region during its evolution, to reach the overall minimum orbital period (19 min) in the grid slice (with an M c i = 1.0 M $ M^{i}_{\mathrm{c}} = 1.0\,\mathrm{M}_\odot $ and P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d). With increasing orbital periods, the binaries go from starting RLO with the companion on the MS (at Porb = 0.5 d, also called case A RLO), to starting RLO when the companion has exhausted core-H (also called case B RLO, at Porb = 1.1 d). We also see this delay in RLO in bottom right in Figure 3, showing the time evolution of the mass-loss rates of the companions. This transition occurs at the bifurcation period (Pbif), because companions in wider orbits have more time to evolve and move out of the MS. In the presented example, Pbif is between Porb = 0.93 d (with P orb i = 2.4 $ P^{i}_{\mathrm{orb}} = 2.4 $ d and M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $) and Porb = 1.1 d (with P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d and M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $). We term systems just below Pbif as near-Pbif binaries; similar studies have shown that the formation of UCXBs from LMXBs happens just below Pbif (Podsiadlowski et al. 2002, 2003; He et al. 2019; Gossage et al. 2023).

thumbnail Fig. 3.

Evolution of selected binaries from our simulated grids for β = 0.3, fpular = 0.0, and M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $. The three different kinds of binaries that we discuss in the text are studied here: case A ( P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $ d), near-Pbif ( P orb i = 2.4 $ P^{i}_{\mathrm{orb}} = 2.4 $ d), and case B or diverging ( P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d). Also shown is the binary that goes down to the orbital period of 19 min (orange curve) for comparison, which is the minimum in the grid with fpulsar = 0.0. Top left shows the orbital period versus companion mass, where Mc is the mass of the companion to the NS and Porb is the orbital period. The formation of a fully convective star and the start of the ejection phase are shown as squares. We also compare our tracks to estimated Mc, min and observed Porb for BWs (faint black crosses), RBs (faint red crosses), and huntsmen spiders (faint green crosses), from Nedreaas (2024, and references therein). The green star is the huntsman candidate 2FGL J0846.0+2820, observed by Fermi-LAT (Swihart et al. 2017). The tidarrens, namely, J1311–3430, J1653–0159, and J0636+5128, are highlighted as pink, blue, and black crosses, respectively. Top right, bottom left, and bottom right show the time evolution of the NS spin, companion mass loss rate, and NS mass, respectively, for the same binaries. The time axis is zoomed in to focus on the spin-up phase.

For binary evolution governed only by conservative mass transfer, the orbit would expand as long as the mass is being transferred from the less massive to the more massive star, while it would contract if the companion star is more massive. With the involvement of magnetic braking, this evolution is not as simple. In NS binaries with MS companions having masses around 1 M, orbital angular momentum evolution is dominated by magnetic braking, which contracts the orbits. Since case A RLO proceeds on a nuclear timescale during the evolution of the stellar MS, any mass loss via winds is not fast enough to dominate over magnetic braking. We can see this nuclear timescale RLO phase in Figure 3 (bottom left panel). In compact binaries with orbital periods ≲1 d, angular momentum losses due to gravitational wave radiation are also significant and further contract the orbit. Once the companion in case A binaries decreases to about 0.2 to 0.3 M, it develops a fully convective structure and magnetic braking stops operating (Rappaport et al. 1983; Spruit & Ritter 1983; Pylyser & Savonije 1989; Podsiadlowski et al. 2002). This is most prominent for the binary in Figure 3 (top left panel) with M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $ and P orb i = 1 $ P^{i}_{\mathrm{orb}} = 1 $ d; there is a decrease in the orbit without any change in mass in the companion mass. The companion detaches from its Roche lobe as the main process driving RLO (the first RLO phase for this binary) has stopped operating, until it fills its Roche lobe again due to gravitational wave radiation continuing to contract the orbit (which is the drop in orbital period seen in the figure). The gravitational wave radiation sustains this second RLO phase and the orbit eventually expands. Correspondingly, we also see this break in the mass-transfer rate in the same figure after which the mass-transfer rate increases again for the second RLO phase. The overall effect of this RLO phase is seen in the top and bottom left panels of Figure 3, case A NSs undergo accretion much more (∼0.85 M) and have an earlier spin-up than other binaries.

For case B RLO (corresponding to M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $ and P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d in Figure 3), strong stellar winds take away stellar angular momentum. Since angular momentum lost due to stellar winds expands the orbits and the effect of magnetic braking is weaker in comparison, we get diverging orbits. The companion has almost an order of magnitude higher mass-loss rate than case A. The binary is at the very end of the companion MS at the onset of RLO (central H abundance less than 0.1); hence, there is a small contraction in the orbit before it expands during the giant phase. The resulting binaries from these diverging systems are NS/He-WDs in wide orbits. Since this is a shorter-lived phase than case A, the NS does not accrete as much for the same initial donor mass (∼0.4 M).

The mass-transfer phase for near-Pbif systems is similar to case B, though the outcomes differ. The companions in near-Pbif binaries have not fully left the MS at the start of RLO, central abundance is close to about 0.1 and magnetic braking is still operating. The stellar radius is about 1.5 times that of the case A companion at the onset of RLO, making the MS winds almost an order of magnitude stronger (∼10−12 M yr−1). The companion becomes fully convective at a significantly lower mass than case A systems (Mc ≲ 0.1 M, see gray and orange lines in the top left panel of Figure 3). The magnetic braking effect increases as the orbit shrinks ( P orb 3 $ {\propto}P_{\mathrm{orb}}^{-3} $; Rappaport et al. 1983). The reason for the fully convective structure forming after the star has lost most of its mass is the duration of the first RLO phase. As the companion loses mass, its central temperature decreases due to loss in internal pressure. The case A binary has a RLO phase that started earlier and lasted longer compared to the near-Pbif binary (seen also in Figure 3, bottom left panel) providing enough time to cool internal structure of the star and making it fully convective while the companion is still on the MS. The near-Pbif binary has stellar winds that, along with the magnetic braking driven RLO, strip the outer layers of the star over a much shorter duration than case A. The shorter the first RLO phase is, the further the orbit contracts into the UCXB region.

The binary that reaches the minimum Porb (19 min; Figure 3 top left panel) in the grid slice (with an M c i = 1.0 M $ M^{i}_{\mathrm{c}} = 1.0\,\mathrm{M}_\odot $ and P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d) forms a fully convective companion after the companion reaches ≲0.02 M. The NS present for this binary also accretes the least in our selected binaries, and the NS spin-up occurs the latest. Toward the end of the second RLO phase, while the companion star is detaching from its Roche lobe, orbital angular momentum losses due to gravitational radiation are still operating (strong for periods ≲1 d). This causes the companion star to fill its Roche lobe again. As the companion is fully convective, it responds by expanding (Hjellming & Webbink 1987) further extracting angular momentum from the orbit (due to spin-orbit coupling). As both angular momentum losses (due to RLO and spin-orbit coupling) reduce the binary separation, the Roche-lobe radius is proportionally affected and the companion overfills its Roche lobe, increasing the mass-transfer rate to ≲10−2 M yr−1 and triggering a dynamical instability.

There is another source in the same observed parameter space as the two tidarrens (see Figure 1), for which the surface H abundance is unknown: PSR J0636+5128 (Stovall et al. 2014; Draghis & Romani 2018). Since it is in the same region of the parameter space as PSR J1311–3430 and PSR J1653–0159 and lies on similar evolutionary tracks, it is very likely this source had a similar formation history and will be H deficient on its surface. Since tidarrens can be produced with any value of fpulsar, the effect of irradiation on their evolution is not significant. They reach their extremely low masses due to the extended RLO phase along with stellar winds during the giant phase, during which they lose their entire H envelope. As we saw in Figure 3, the lower-right part of the RB population is also reproduced. These systems are the very same that pass through the UCXB region and lead to the tidarrens, and correspond to systems that start their RLO just below Pbif. As there is no irradiation included in their evolution, these systems could depict some of the non-irradiated RBs observed. Since all currently observed BWs and about half of the currently observed RBs are strongly irradiated (Turchetta et al. 2023), we investigate higher values of fpulsar to form RBs in the subsequent sections.

3.3. Effect of pulsar wind irradiation

Figure 1 shows the evolutionary grids with an M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $ and varying pulsar wind efficiency (fpulsar) and we see the need for strong irradiation in order to explain all of the observed spiders. We also show the initial parameter space for spiders in Figure 4, their dynamical end, as well as the final orbital periods in the color bar, for fpulsar = 0.0 and 0.5 and initial NS masses 1.3 M. Many of the binaries shown in Figure 4 (around 50 to 60% in all simulated spiders) end up suffering a dynamical instability (shown by triangle symbols in the figure) divided into two categories; those that become unstable during the first RLO phase and those that become unstable toward the end of the second RLO phase. For these binaries in the second category, the NS still manages to accrete enough material to spin it up to milliseconds during the preceding case A RLO. Hence, these could be seen as spider binaries until the onset of instability, after which they would be seen as isolated spinning pulsars. Figure 5 shows the total mass-loss rates and the mass-transfer rates of the companion versus time of three types of interacting binaries (with pulsar wind irradiation) we have studied so far: case A, near-Pbif, and case B binaries. The binaries have initial orbital properties as M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $ and M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, and three values for P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $, 2.4, and 2.6 d. Going from low irradiation (fpulsar = 0.05) to high (fpulsar = 0.5), we see the mass-loss rates spike to very high values for case A and near-Pbif binaries, which could instigate an instability.

thumbnail Fig. 4.

Allowed initial parameter space to form spider MSPs for an initial NS of mass 1.3 M and fpulsar values of 0.0 (left) and 0.5 (right). The term M c i $ M^{i}_{\mathrm{c}} $ is the initial mass of the companion to the NS, and P orb i $ P^{i}_{\mathrm{orb}} $ is the initial orbital period. The color bar shows the final orbital periods ( P orb f $ P^{f}_{\mathrm{orb}} $) and the gray symbols did not attain MSP-like spins. We show the end states of the tracks, if they encountered a dynamical instability (with gray triangles or with colorful triangles depending on if onset of instability occurred during the first RLO phase or the second, respectively) or ended with one of the criteria for successful termination (star symbols). The diverging symbols are shown by gray crosses. We also show the near-Pbif binaries that passed through the RB part of the observed parameter space enclosed by dashed orange lines. The magenta triangle shows the initial parameters of the binary that reproduces PSR J0952–0607 (with observed spin period of 1.41 ms).

thumbnail Fig. 5.

Evolution of the mass-loss rates (c, black curves) and the mass-transfer rates (transfer, blue curves) of the companion star. The initial binary masses are M c i = 1.0 M $ M^{i}_{\mathrm{c}} = 1.0\,\mathrm{M}_\odot $ and M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, with β = 0.3. We present the three cases of interacting binaries studied, case A (top row; with P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $ d), near-Pbif (middle row; with P orb i = 2.4 $ P^{i}_{\mathrm{orb}} = 2.4 $ d), and case B (bottom row; with P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d), each for three values of fpulsar = 0.05 (left column), 0.1 (middle column), and 0.5 (right column).

We discussed the occurrence of dynamical instability toward the end of the second RLO phase in the previous section. The presence of irradiation would only increase the occurrence of this instability. Irradiation removes matter from the surface of the companion star, causing rapid stellar expansion. Increased fpulsar results in a higher number of dynamical unstable systems (goes from 50% to 62% in simulated spiders when going from 0.0 to 0.5 in fpulsar). We include these dynamically unstable systems in our study since by the time the instability onset happens, the average MSP binary is already ∼1010 yr old and could have been observed during the preceding MSP phase. Hence, there is a reasonable chance to observe these binaries as BW binaries within a Hubble time (1.4 × 1010 yr). For this part of the parameter space the companion star would be destroyed due the instability. Therefore, it could provide an additional channel to the population of Galactic isolated MSPs (20% of the total Galactic MSP population) that is observed (Manchester et al. 2005a), with previous RLO spinning up the NS and subsequent merging destroying the companion. The idea that isolated MSPs could originate from BWs irradiating away their companions has been previously suggested both theoretically (van den Heuvel & van Paradijs 1988) and through observations of MSPs with planets/asteroid belts (Wang et al. 2009; Shannon et al. 2013).

Going from fpulsar = 0.0 to fpulsar = 0.5 in the converging systems (Figure 4), the final orbital periods widen significantly and this effect depends on the initial properties of binary. The widest orbits are seen with M c i = 1.0 M $ M^{i}_{\mathrm{c}} = 1.0\,\mathrm{M}_\odot $, where the final orbital periods are ≳10 d. For higher companion masses, the final orbital periods lie in the range of 0.1 to 1 d. Figure 6 shows the orbital evolution corresponding to three M c i $ M^{i}_{\mathrm{c}} $ (1.0, 1.5, and 2.5 M) and fpulsar = 0.0 to fpulsar = 0.5, for various values of P orb i $ P^{i}_{\mathrm{orb}} $. For case A RLO binaries ( P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $ d; shown by black lines), we see increased orbital widening with fpulsar = 0.5. The companions to the NS become fully convective around 0.2 M, after which the pulsar wind is turned on. The orbit responds to the mass loss by expanding depending on fpulsar. Going from M c i = 1.0 M $ M^{i}_{\mathrm{c}} = 1.0\,\mathrm{M}_\odot $ to 1.5 M (for fpulsar = 0.5) the expanding orbits cover a larger part of the observed BW space. This is due the fact that all of these binaries are undergoing RLO on a nuclear timescale. For a fixed NS mass and orbital period, more massive companions (or higher q) lead to more accretion onto the NS and a higher spin-down luminosity (and correspondingly high orbital expansion). However, with an M c i = 4.0 M $ M^{i}_{\mathrm{c}} = 4.0\,\mathrm{M}_\odot $, the accreted matter decreases down to about 0.004 M as the binary suffers a dynamical instability soon after the first RLO-onset. This is significantly less than 1.0 M accreted with M c i = 3.0 M $ M^{i}_{\mathrm{c}} = 3.0\,\mathrm{M}_\odot $. Hence, there is a limit to how high in companion mass we can go with a fixed initial NS mass and reproduce BWs. We can see this clearly by looking at Figure 4. The effect of irradiation on RLO can be seen in Figure 5 (top row) with increasing mass loss after the first RLO phase with increasing fpulsar.

thumbnail Fig. 6.

Orbital period evolution versus companion mass of selected binaries from our simulated grid shown in Figure 1. The term Mc is the mass of the companion to the NS, and Porb is the orbital period. The figure shows the effect of pulsar wind irradiation, comparing fpulsar = 0.0 (dashed lines) and 0.5 (solid lines). The left panel shows the tracks with M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $, P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $, 2.4, and 2.6 d. The right panel shows two initial companion masses, M c i = 1.0 $ M^{i}_{\mathrm{c}} = 1.0 $ (with P orb i = 1 $ P^{i}_{\mathrm{orb}} = 1 $ d) and 2.5 M ( P orb i = 2.0 $ P^{i}_{\mathrm{orb}} = 2.0 $ and 2.2 d). We compared our tracks to Mc, min and Porb for BWs (faint black crosses), RBs (faint red crosses), and huntsmen spiders (faint green crosses), from Nedreaas (2024, and references therein). The green star is the huntsman candidate 2FGL J0846.0+2820, observed by Fermi-LAT (Swihart et al. 2017). The magenta stars show the three transitional MSPs (tMSPs) observed, PSR J1023+0038 (Mc = 0.22 M and Porb = 4.75 h; Shahbaz et al. 2019; Thorstensen & Armstrong 2005), IGR J18245–2452 (Mc, min = 0.17 M and Porb = 11.03 h; Papitto et al. 2013), and XSS J12270–4859 (Mc, min = 0.15 M and Porb = 6.91 h; Roy et al. 2015; de Martino et al. 2015; Bassa et al. 2014).

Next we look at the tracks for P orb i = 2.4 $ P^{i}_{\mathrm{orb}} = 2.4 $ d in Figure 6, corresponding to the gray curves (left panel). For M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $, the effect of fpulsar = 0.5 is significant on the orbital evolution. The near-Pbif binary maintains a relatively wide orbit (around 0.2 d) crossing the RB regime. With increasing fpulsar, the duration of the mass-transfer phase is shortened as majority of mass loss is due to irradiation (see Figure 5). The radio pulsar wind will start as long as the two conditions described in Section 2.3.3 are met. Hence, strong irradiation disrupts the RLO phase as the companion detaches more easily from its Roche lobe due to higher mass loss. This does not significantly reduce the number of spiders binaries produced because most of the accretion happens early on in the evolution. In Figure 3, we see the NSs spin-up to milliseconds toward the start of RLO, for all the binaries presented. Going to M c i = 2.5 M $ M^{i}_{\mathrm{c}} = 2.5\,\mathrm{M}_\odot $ and P orb i = 2.0 $ P^{i}_{\mathrm{orb}} = 2.0 $ d (gray curves in right panel), the near-Pbif tracks no longer pass through the RB region with the heavier companion (irrespective of fpulsar). For intermediate-mass companions (≳2.0 M), the initial angular momentum loss by RLO dominates, as magnetic braking is not operating. For very short durations, the mass-transfer rate could even breach the Eddington limit, which results in any excess material being lost from the system taking away even more angular momentum. This continues until the star reaches around 1.0 M (and 1 d in Figure 6) and magnetic braking turns on. The following evolution is governed by magnetic braking losses and behaves similar to case A RLO, during its first RLO phase. The fast initial mass loss causes the initially 2.5 M star to become fully convective around 0.2 M, after which it follows a UCXB-like evolution. Hence, initial companion masses ≳2.0 M in near-Pbif binaries do not reproduce RBs (also shown in Figure 4). If a heavy NS is present (∼2.0 M), it could be argued that correspondingly heavier companions (≳2.0 M) could form RBs. However, NSs with masses 1.8 M and above very easily cross the maximum NS mass threshold and very few binaries in those grids would be viable spiders (see Section 3.1).

We can determine also the evolutionary history of PSR J0952–0607, the fastest spinning Galactic BW so far (with spin 1.41 ms), from our simulations, considering β = 0.3 and fpulsar = 0.5. To account for uncertainties during mass transfer and in our calculation of Pspin, we consider a 10% margin and take all binary tracks with Pspin ≲ 1.55 ms. The tracks are shown in Figure 7. We compare the simulated tracks to the estimated orbital parameters, M NS = 2 . 35 0.17 + 0.17 M $ M_{\mathrm{NS}} = 2.35^{+0.17}_{-0.17}\,\mathrm{M}_\odot $, M c = 0 . 032 0.002 + 0.002 M $ M_{\mathrm{c}} = 0.032^{+0.002}_{-0.002}\,\mathrm{M}_\odot $, and Porb = 6.42 h (Romani et al. 2022), and select the best matching track (within the errors). The best matching track started with M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, M c i = 2.0 M $ M^{i}_{\mathrm{c}} = 2.0\,\mathrm{M}_\odot $ and P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $ d (also shown in Figure 4). The binary underwent case A RLO and ended due to dynamically unstable RLO from a fully convective companion; the final NS mass was 2.35 M and final NS spin 1.2 ms. Hence, strong irradiation of the companion explains this source.

thumbnail Fig. 7.

Orbital period versus companion mass for diverging binary tracks from our simulated grids that have a similar NS spin as the observed PSR J0952–0607 (1.41 ms) within a 10% margin (dashed gray curves). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. All the tracks shown have β = 0.3 and fpulsar = 0.5. We further show the binary that best matches the sources in its orbital evolution in blue curve, which has initial companion mass of 2.0 M, initial NS mass of 1.3 M, and initial orbital period of 1.0 d.

3.4. Huntsman spiders

Huntsman spiders are considered to comprise of a red giant and a pulsar in a widening orbit (Strader et al. 2015, 2019). They evolve from a case B RLO phase in LMXBs. In terms of companion masses, they are similar to RBs, with the difference being that their orbits are an order of magnitude wider than RBs. Looking at Figure 1, we see that the binary evolution tracks reproduce the huntsmen regime irrespective of the value of fpulsar. Contrary to RBs and BWs, with a high fpulsar the orbits do not widen with increasing fpulsar. Figure 6 (left panel, with M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $ and P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d, shown by blue curves) shows that wider orbits correspond to fpulsar = 0.0. Since these binaries are undergoing case B RLO, the companions to the NS are giant stars that have developed a deep convective envelope. Convective envelopes expand rapidly on mass loss and drive a high mass-transfer rate. When the radio pulsar wind turns on and irradiates the companion, it further increases the instantaneous mass loss from the companion. The companion detaches from its Roche lobe, pausing the RLO phase. We see the shortened RLO phase in Figure 5 (bottom row). The resulting effect is a decreased expansion of the orbits.

We also find more conservative mass-accretion rates better reproducing the only current huntsman system with an observed spin period (2.66 ms for the source 1FGL J1417.7–440; Camilo et al. 2016). Table 1 describes the various properties of this huntsman source, as well as the other huntsman known so far, 2FGL J0846.0+2820. Figure 8 shows the final spin period distributions versus the final NS masses of binary tracks that end up in the huntsmen regime, corresponding to the three β values investigated in our study (0.0, 0.3, and 0.7). Compared to BWs and RBs, more conservative mass transfer (β < 0.3) is required for an NS with an initial mass of 1.3 M to spin up to 2.66 ms. With a β = 0.7, the minimum spin period reached is 5.83 ms. With a β = 0.3, only NSs with initial masses ≳1.5 M are able to reproduce the pulsar spin of 1FGL J1417.7–440. We also compare our simulations to the estimated masses of the pulsars in 1FGL J1417.7–440 ( 1 . 62 0.17 + 0.43 $ 1.62^{+0.43}_{-0.17} $; Swihart et al. 2018) and 2FGL J0846.0+2820 ( 1 . 96 0.41 + 0.41 M $ 1.96^{+0.41}_{-0.41}\,\mathrm{M}_\odot $; Swihart et al. 2017). If the binaries had fully conservative RLO (β = 0.0), for a huntsman to have an NS with mass greater than 2.0 M, it should have M NS i > 1.3 M $ M^{i}_{\mathrm{NS}} > 1.3\,\mathrm{M}_\odot $ (with 1.3 M, the maximum mass reached is 1.85 M). Therefore, a huntsman with a confirmed high mass pulsar would have started RLO with a pulsar heavier than the canonical mass of 1.3 M.

thumbnail Fig. 8.

Final NS spins ( P spin f $ P^{f}_{\mathrm{spin}} $) versus final NS masses ( M NS f $ M^{f}_{\mathrm{NS}} $) in diverging binaries for β = 0.0 (blue colored symbols), β = 0.3 (orange colored symbols), and β = 0.7 (green colored symbols). We group the binaries by the initial NS mass ( M NS i $ M^{i}_{\mathrm{NS}} $) corresponding to different symbols (see legend). We also compare to the observed spin and estimated pulsar mass of the huntsman source 1FGL J1417.7–440 (2.66 ms and 1 . 62 0.17 + 0.43 M $ 1.62^{+0.43}_{-0.17}\,\mathrm{M}_\odot $, shown by the gray dashed lines and gray shaded area for the uncertainties; Camilo et al. 2016; Swihart et al. 2018) and the estimated pulsar mass of the huntsman candidate 2FGL J0846.0+2820 ( 1 . 96 0.41 + 0.41 M $ 1.96^{+0.41}_{-0.41}\,\mathrm{M}_\odot $, shown by the vertical pink dashed lines and pink shaded area for the uncertainties; Swihart et al. 2017).

The observed and estimated properties of the two currently known huntsmen spiders are described in Table 1. The spin of 1FGL J1417.7–440 (2.66 ms) can be used to constrain its evolutionary history. To reproduce the observed spin and estimated current pulsar mass with an initially canonical NS, we need a fairly conservative RLO phase (β ≲ 0.3; see Figure 8). We also use an upper limit of 2.926 ms and a lower limit of 2.394 ms on the final spin to better identify the matching binaries, to account for 20% uncertainty in our calculations for the spin. Though pulsar irradiation is not necessary to reproduce the orbital properties of huntsman spiders, since it does have an effect on the orbit, we compare the different values of fpulsar. In Figure 10 we show the binary tracks that satisfy our selection criteria. We also compare these to the estimated orbital properties of 1FGL J1417.7–440 (mass of 0 . 28 0.03 + 0.07 M $ 0.28^{+0.07}_{-0.03}\,\mathrm{M}_\odot $ and orbital period 5.37 d; Swihart et al. 2018; Strader et al. 2015), and see what tracks pass through the observed properties. We find that the binary track that best matches the source’s observed orbital properties (within the errors) corresponds to β = 0.0 and fpulsar = 0.1, and with the initial parameters, M c i = 2.0 M $ M^{i}_{\mathrm{c}} = 2.0\,\mathrm{M}_\odot $, M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, and P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d. The final spin of the pulsar is 2.73 ms, which is within 3% of the observed spin. This binary ends as a He-WD with mass 0.25 M and a 1.73 M NS in an orbit of 6.46 d.

Table 1.

Observed and estimated properties of the two observed huntsmen spiders.

As for the source 2FGL J0846.0+2820, since we do not know its spin, there is a wider range of binary tracks that could be used to explain it, and almost all fpulsar values can reproduce its orbital parameters (see Figure 1). Compared to 1FGL J1417.7–440 (and other spider systems), 2FGL J0846.0+2820 has a non-negligible eccentricity in its orbit (0.06; Swihart et al. 2017). This contradicts the general assumption that these systems are old enough for tidal forces to have efficiently circularized the orbits (Verbunt & Phinney 1995). This observed eccentricity could be the result of the interaction with a circumbinary disk (Antoniadis 2014) or inefficient tides (such as in radiative envelopes, which would be the case for MS stars; Hurley et al. 2002).

Although we can reproduce the orbital parameters observed, we also need to study which diverging LMXB/IMXB systems can spin up NSs to MSP spins. Figure 9 shows the allowed initial parameter space for the spin-up of NSs in diverging LMXBs/IMXBs. The presented figure shows grids of initial pulsar masses of 1.3 and 2.0 M, for β = 0.3 and fpulsar = 0.0. We present the figures with β = 0.3 instead of 0.0 as we leave room for uncertainties during the RLO phase that were not accounted for in our simulations. The color bar highlights the region of the initial parameter space where the NS was able to spin up to ≤30 ms. Not all diverging binaries will be able to spin up to MSP range. As we go from a less massive NS to a more massive NS, the peak of the spin distribution (systems with the maximum spin-up due to accretion) moves to higher companion masses and lower orbital periods, from centered around 2.0 M and 2.5–3.5 d to being centered around 2.5 M and 1.6–3.0 d.

thumbnail Fig. 9.

Allowed initial parameter space to form huntsmen spiders (fpulsar = 0.0) for an initial NS of mass 1.3 M (left) and 2.0 M (right). The term M c i $ M^{i}_{\mathrm{c}} $ is the initial mass of the companion to the NS, and P orb i $ P^{i}_{\mathrm{orb}} $ is the initial orbital period. The color bar shows the final NS spins ( P spin f $ P^{f}_{\mathrm{spin}} $), and the gray stars denote the systems that formed diverging orbits but did not spin up to the MSP range. We show the end states of the tracks, if they encountered a dynamical instability (triangle symbols) or ended with one of the criteria for successful termination (star symbols). Converging tracks are shown by gray crosses.

thumbnail Fig. 10.

Orbital period versus companion mass for diverging binary tracks from our simulated grids that have a similar NS spin as the observed huntsman source 1FGL J1417.7–440 (2.66 ms) within a 20% margin (dashed gray curves). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. We further show the binary that best matches the sources in its orbital evolution in blue curves, where β = 0.0, fpulsar = 0.1, M c i = 2.0 M $ M^{i}_{\mathrm{c}} = 2.0\,\mathrm{M}_\odot $, M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, and P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d.

The boundary of the huntsmen parameter space depends on the initial q and P orb i $ P^{i}_{\mathrm{orb}} $ of the binary (see Figure 9). The range of acceptable q is almost constant for both NS masses and corresponds to regions of high accreted mass. If the companion to the NS is not massive enough, the mass-transfer rate will be too low for efficient NS spin-up. If it is too massive, the RLO phase is very short as the star evolves quick. If the RLO phase is super-Eddington, accretion would be further limited. The top boundary of the parameter space is limited due to the evolutionary state of the companion. Wider orbits than ∼3.5 d results in the companion filling its Roche lobe after it has developed a deep convective envelope (late case B RLO), making the RLO phase short. In Figure 9 (right panel, with M NS i = 2.0 M $ M^{i}_{\mathrm{NS}} = 2.0\,\mathrm{M}_\odot $), we also see an isolated binary with M c i = 4.5 M $ M^{i}_{\mathrm{c}} = 4.5\,\mathrm{M}_\odot $ and P orb i = 1.8 $ P^{i}_{\mathrm{orb}} = 1.8 $ d that has an NS with MSP-like spin. The final binary is a wide He-WD/NS binary (final WD mass 0.31 M and final period 33 d) with a final NS mass of 2.26 M. The neighboring binaries do not drive high enough mass-transfer to spin up their NS, and we did not find a similar binary with a lower initial NS mass.

3.5. Transitional millisecond pulsars

There is also a sub-category of compact MSPs that have been observed, binaries that switch from the rotation-powered state (where the spin-down of the pulsar dominates the evolution) to accretion-powered state (mass-transfer from the companion dominates) and vice-versa (see review by Papitto & de Martino 2022). These binaries are called transitional MSPs (tMSPs), and so far there are 3 observed, PSR J1023+0038 (Archibald et al. 2009), IGR J18245–2452 (Papitto et al. 2013), and XSS J12270–4859 (Bassa et al. 2014). The companion masses for all of them are in the range 0.15 to 0.2 M and orbital periods 4 to 11 h (Papitto & de Martino 2022, and references therein), which put them in the RB regime of the observed parameter space. Figure 6 shows the three tMSPs along with some of the tracks from our simulated binaries. Two of these lie the near-Pbif tracks for M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $ (gray lines in the figure), where the tracks differentiate for fpulsar = 0.0 and 0.5. The point of differentiation of the two tracks is the turning on of the pulsar wind, and it explains the transitioning of the MSP binaries from accretion-powered to rotation powered state. Some other factors that should be considered to explain tMSPs are accretion disk instabilities (Lasota 2001), which are not yet fully understood.

4. Discussion

4.1. Spider binary formation

We can compare our simulated spiders to those from previous works. While we find the formation of RBs and BWs to be linked, Chen et al. (2013) found RBs and BWs to come from two distinct populations of LMXBs, depending on the irradiation efficiency. For BWs, they found fpulsar ≲ 0.08 to be the required irradiation efficiency, while for RBs, it was fpulsar ≳ 0.08. They explored mainly case A RLO systems, M c i $ M^{i}_{\mathrm{c}} $ around 1.0 M, M NS i $ M^{i}_{\mathrm{NS}} $ around 1.4 M, and P orb i $ P^{i}_{\mathrm{orb}} $ from 0.08 to 1.4 d, and found increased orbital expansion with increasing fpulsar. Qualitatively, we also see similar behavior in our results (see Section 3.3). Particularly for BWs, looking at similar systems in Figure 1, for fpular = 0.0 to 0.05, we cover the lower-left region of the observed BW space. The degree to which the orbit expands for fpular ≳ 0.1 is lesser than Chen et al. (2013). Using different criteria for the start of the radio pulsar phase than Chen et al. (2013), we see the phase start at similar evolutionary stages for similar systems (at donor mass of 0.2 − 0.3 M for the case A XRBs).

For RBs, Chen et al. (2013) found them to be produced from case A binaries (∼1 d) with strong irradiation (fpulsar ≳ 0.08). In our work, the majority of RBs are produced at much wider initial orbits (near-Pbif at RLO onset, starting at around 2.0 to 2.4 d at ZAMS), along with high irradiation. With fpulsar = 0.5, we also see rapid expansion of the case A binaries (solid black curve and M c i = 1.5 $ M^{i}_{\mathrm{c}} = 1.5 $ in the right panel of Figure 6). However, this expansion in orbits is not enough to bring them completely into the RB regime. Benvenuto et al. (2014, 2015) studied the effect of X-rays on the companion surface, as well as pulsar wind irradiation. They found that all BWs should have passed through a RB stage previously. We find that while not all BWs passed through the RB region, most binary tracks passing RB tracks ended up in the BW region of the observed parameter space, and the BWs that passed through the RB region, were H-surface deficient and in narrow orbits. The number of observed BWs to RBs observed is similar (approximately 63% and 37%, respectively), which implies their evolutionary lifetimes are comparable. For the BWs and RBs in our simulated tracks, by the time the spider phase begins the binaries are ≳109 yr old, as most of the prior evolution happened during the MS phase, which is long for low-mass stars.

The need for a persistent magnetic braking was also suggested by Ginzburg & Quataert (2020, 2021), where they found that magnetic braking must remain active even when the companion has developed a fully convective structure and the pulsar wind has turned on. While the authors focused their works on BWs, they suggested that the pulsar wind is too weak to remove enough mass from the pulsar companions. They looked at companions with initial masses of 1.0 M. Their conclusion is similar to our results, as we also find the binary tracks not reproducing all of the observed BWs for initial companion masses of 1.0 M (see fpulsar = 0.5 in right panel of Figure 6). With increasing companion masses, we cover more of the observed BW parameter space. The MSP source PSR J1953+1844 (in the globular cluster M71) that was recently discovered, seems to be the evolutionary link between tidarrens and RBs (Pan et al. 2023). It has an orbital period of 53.3 min and an estimated companion mass of 0.07 M. Since it is associated with a globular cluster (and would hence have a lower metallicity) and our work seeks to study Galactic spiders, we do not include it in our observed sources catalog. Still, Pan et al. (2023) suggested that this source was formed from a similar process as UCXBs (with P orb i $ P^{i}_{\mathrm{orb}} $ values near Pbif), which is similar to our results, albeit at different metallicities.

4.2. Accretion efficiency and non-conservative evolution

As we saw in Section 2.3, in the literature non-conservative mass-accretion efficiencies (β ≳ 0.5) are favored since they can reproduce the LMXB observations. We find that in order to reproduce the observed properties of spider binaries and the observed spins, more conservative efficiencies are better suited (β ≲ 0.3). An assumption of non-conservative mass accretion throughout the entire evolution is proxy for any instabilities that might be affecting accretion (Illarionov & Sunyaev 1975; Pringle 1981; van Paradijs 1996, for instance, propeller effects or thermal-viscous instabilities). We include the propeller phase in the NS spin evolution (see Section 2.3.2) by removing 10% of the transferred mass due to the propeller action. Since this phase is highly uncertain even in its effectiveness, we assume a low efficiency of loss of angular momentum during this phase. We also assumed an full efficiency for transfer of angular momentum from the accreted matter during RLO (see Section 2.3.1). Our results are similar to those by Kar et al. (2024), who found that lower β values are required to reproduce accreting MSPs.

As for accretion disk instabilities, we include a model that describes the presence of a accretion disk above a critical mass-transfer rate, below which an X-ray irradiated disk would be unstable (see Section 2.3). Since we have already included inefficient mass accretion due to these two processes, the mass-accretion efficiency (ϵ = 1 − β) that is set to a fixed value during the entire binary evolution accounts for any additional causes of mass or angular momentum loss. It also accounts for any uncertainties in our models of inefficient accretion, since we use simplistic approximations. A binary would evolve similarly due to a high value of fixed β compared to a low β with restricted accretion at some instances. To test this, we ran MESA models without the X-ray irradiated accretion disk assumptions and with β = 0.7 and 0.3, and fpulsar = 0.0. The initial properties were M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $, M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, and P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $ d. We compared it to runs with same initial parameters and with the X-ray irradiated disk assumption, and found no significant differences in the evolution.

4.3. Pulsar spin evolution and sub-millisecond pulsars

There is a lot of uncertainty in our understanding of pulsar spin evolution, including pulsar magnetic field evolution and the pulsar braking index (see review by Abolmasov et al. 2024). The spin down for an isolated pulsar follows a power law that dictates the rate of change in the angular frequency ( Ω ˙ $ \dot{\Omega} $) with the angular frequency (Ω), which is Ω ˙ Ω 3 $ \dot{\Omega}\propto -\Omega^{3} $. The range of values generally explored for the braking index are from 2.5 to 6.0 (Archibald et al. 2016; Woan et al. 2018), with values closer to 3.0 considered more canonical (Sturrock 1971; Ruderman & Sutherland 1975; Manchester et al. 2005b; Hamil et al. 2015). Still there are observations of the braking index in the wide range of 10−6 to 106 (Hobbs et al. 2004; Biryukov et al. 2012; Zhang & Xie 2012).

Studies have suggested that MSPs are close to their spin equilibrium, determined by the accretion disk and magnetosphere interaction (Patruno et al. 2012). Tauris et al. (2012) studied the formation of MSPs and derived the equilibrium spin depending on the mass accreted onto the pulsar. They suggested that since most LMXBs undergo non-conservative mass transfer, the MSP can be spun to ≳1.4 ms. Similarly, in our results, we find a minimum spin of around 1.38 ms (for β = 0.3; see Section 3.1), after disregarding systems where the NS crossed the maximum NS mass limit (taken to be 2.5 M). Figure 11 shows the accreted matter onto the NS versus the final NS spin from our simulated grid (for fpulsar = 0.0). The correlation is compared to the equilibrium spin period calculations from Tauris et al. (2012); our simulated NSs require more accreted mass for the same spin-up in the MSP regime. Reasons for discrepancy could be existing uncertainties in spin evolution and accretion physics. We separate out the cases where the final NS mass was greater than 2.5 M (shown as stars in the figure). We find a correlation between the NS spin and the accreted mass, with sub-millisecond spins produced only by NSs beyond the 2.5 M mass limit.

thumbnail Fig. 11.

Accreted mass versus final NS spins ( P spin f $ P^{f}_{\mathrm{spin}} $) in our simulated grid with β = 0.3. We show the different initial masses as different colors (shown in the legend) and distinguish between the NSs that had M NS f 2.5 M $ M^{f}_{\mathrm{NS}}\ge 2.5\,\mathrm{M}_\odot $ (crosses) and those where M NS f < 2.5 M $ M^{f}_{\mathrm{NS}} < 2.5\,\mathrm{M}_\odot $ (circles). We also compare our results to the equilibrium spin period calculations from Tauris et al. (2012), done for M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $ (light blue curve) and 2.0 M (dark orange curve). The black dotted line shows the observed spin for PSR J0952–0607 (spin period of 1.41 ms) and the gray dotted line shows the upper limit for MSPs considered (30 ms).

According to Figure 11, NSs in LMXBs and IMXBs accrete between 0.2 and 1.0 M to spin up to millisecond periods ( 1.38 ms P spin f 30.0 ms $ 1.38\,\mathrm{ms}\leq P^{f}_{\mathrm{spin}}\leq 30.0\,\mathrm{ms} $). The NSs with initially higher masses (≳1.8 M) accrete more than lighter NSs and hence spin up to a higher degree. Naturally, some of these heavy accretors could exceed the NS mass limit and enter the sub-MSP regime. Accretion induced collapse (AIC) of massive NSs was also studied by Chen et al. (2023) where they found inefficient RLO required initial NS masses of at least 2.16 M in order for AIC to occure. Haskell et al. (2018) found that masses greater than 2.0 M have an upper limiting frequency of around 810 Hz (lower Pspin limit of 1.23 ms) depending on the NS equation of state, while lighter NSs could spin faster (up to 1200 Hz or 0.83 ms), inconsistent with observations. The spin at which an accreting NS would collapse into a BH also depends on the initial mass: for M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, it is at 1.38 ms and for M NS i = 2.0 M $ M^{i}_{\mathrm{NS}} = 2.0\,\mathrm{M}_\odot $, it is at 2.66 ms. For an NS lighter than M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, the switch from NS to BH could occur at a spin lower than 1.38 ms. While our resulting spins are not as small as those by Haskell et al. (2018), we see a similar trend. For masses ≳2.0 M the equation-of-state could change, leading to a collapse of the NS into a BH, or the formation of a hybrid star with sub-MSP spin (Gerlach 1968; Glendenning & Kettner 2000; Bejger et al. 2017). Liu et al. (2024) accounted for selection bias due to pulse broadening and suggested chances for Square Kilometer Array to observe sub-millisecond pulsar spins would be negligible. Since there are no current observations of any such sub-MSP star, sub-MSPs may not exist if those NSs have collapsed into BH.

4.4. Uncertainties in redback formation

We find an evolutionary connection between UCXBs, RBs, and tiderrans (see Section 3). Our results for the formation of UCXBs are similar to other studies where they find UCXBs originating from LMXBs starting from close to Pbif (Podsiadlowski et al. 2002; van der Sluys et al. 2005a,b; Istrate et al. 2014; He et al. 2019). The idea that the Tidarren PSR J1311–3430 went through a UCXB phase was also hinted at by Chen et al. (2013) and. There is the fine-tuning problem in the formation of UCXBs. With the standard magnetic braking prescription (Skumanich 1972; Verbunt & Zwaan 1981; Rappaport et al. 1983) often used (also used in our study), there is a narrow range of orbital periods (that is near-Pbif) and companion masses from which UCXBs can be produced (van der Sluys et al. 2005a,b; Istrate et al. 2014). Even though strong irradiation widens near-Pbif orbits and some of them might no longer be UCXBs, some UCXBs can still be reproduced (see fpulsar = 0.5 in Figure 1).

The dependence of Pbif on magnetic braking assumptions and the effects on LMXB evolution have been extensively studied (Podsiadlowski et al. 2002; Li et al. 2004; Ma & Li 2009; Chen & Podsiadlowski 2016; Deng et al. 2021; El-Badry et al. 2022; Echeveste et al. 2024), and it has been suggested that a stronger magnetic braking process is needed to form UCXBs (and consequently RBs) from a larger initial parameter space. Chen et al. (2021) tested the magnetic braking prescription suggested by Van & Ivanova (2019) in the formation of binary MSPs and UCXBs. They found that while efficient magnetic braking does increase the allowed initial parameter space to form UCXBs, wide-orbit binary MSPs are not reproduced. In regards to spiders, Ginzburg & Quataert (2020) and Ginzburg & Quataert (2021) focused on reproducing the observed BW properties and proposed a sustained magnetic braking prescription that will drive RLO; magnetic braking will not stop operating at Mc ∼ 0.2 M and the companion has lost its radiative core. Gossage et al. (2023) compared MESA simulations with single star and LMXB observations and also found saturated magnetic braking to best reproduce both types of sources. Still the standard magnetic braking prescription can reproduce LMXBs with periods less than one hour to around ten days, better than some other proposed prescriptions (Echeveste et al. 2024; Gossage et al. 2023). Carrying out a population synthesis study of UCXBs would also tell us what fraction of LMXBs in a stellar population start RLO close to their Pbif. Hence, to disentangle the various formation channels and types of magnetic braking prescriptions, a combined investigation of detailed simulations and population synthesis analysis would be required, which would further help in understanding the formation of spiders.

5. Conclusions

Compact MSP binaries (or spider binaries) are good laboratories for studying binary evolution. In this work, we calculated the NS spin evolution and identified binaries that formed MSPs (using the criteria for NS spins: P pular f 30 $ P^{f}_{\mathrm{pular}}\le 30 $ ms). Assuming a model for an X-ray irradiated accretion disk and three different values of the accretion efficiency (0.0, 0.3, and 0.7), we ran simulation grids spanning four values of initial NS mass (1.3 M, 1.5 M, 1.8 M, and 2.0 M), a range of initial companion masses (1.0 M to 4.5 M), and initial orbital periods (0.6 d to 4.0 d). Following the RLO phase, as the pulsar spins down, it loses rotational energy in the form of spin-down luminosity, ablating the outer envelope of the companion star and causing changes in the orbital evolution. We used four different values for the efficiency of this process (0.0, 0.05, 0.1, and 0.5) and presented the allowed initial parameter space to produce BWs, RBs, and huntsmen spiders. After comparing our simulated grids to the observed spider populations, we arrived at the following conclusions:

  • With the assumption of an unstable accretion disk due to X-ray irradiation during the mass-transfer phase, an accretion efficiency of at least 70% is needed in order to reproduce the observed spins and pulsar masses of BW systems. High accretion efficiency (> 70%) also allows for the formation of the huntsmen spiders from a binary with a canonical NS in a widening orbit during RLO.

  • For BW formation, a high pulsar wind efficiency (fpulsar = 0.5) is required to reproduce the full observed parameter space, although lower efficiencies (0.0 < fpulsar < 0.5) can still reproduce part of the systems. However, the tidarrens (J1311–3430, J1653–0159, and J0636+5128; BWs with orbital periods ≲1.0 d and 0.01 M companions) require no irradiation. We also predict that J0636+5128 (similar to its neighbors) should not have any hydrogen on its companion’s surface, as the companion loses its hydrogen envelope through stable mass transfer. The binary tracks passing through the tidarren region also pass through the RB region, connecting BWs and RBs as part of the same initial population. These binaries start RLO at near-Pbif and pass through an UCXB phase. Their evolution is dependent on the magnetic braking prescription used.

  • With high pulsar wind irradiation (fpulsar = 0.5), the orbits of the near-Pbif binaries expand drastically, covering most of the RB region. This is true for the companion masses ≲1.5 M. The near-Pbif binaries are the main formation channel for RBs. Hence, RBs are subject to the fine-tuning problem seen in UCXBs, where only a narrow range of values for the initial orbital parameter space, with standard magnetic braking, can reproduce the observations.

  • We identified the initial parameter space for diverging NS binaries to spin up as MSPs. These systems would form the small population of huntsmen spiders that have been observed. With fully conservative accretion (assuming an X-ray irradiated accretion disk), an NS in diverging orbits with an initial mass of 1.3 M gains a maximum of 0.55 M. If a huntsman pulsar is discovered to have a pulsar mass ≳2.0 M, it was most likely more massive than 1.3 M when formed.

  • We identified the prior evolution of the pulsars BW PSR J0952–0607 and huntsman 1FGL J1417.7–440 by matching their observed spins and orbital properties to binary tracks. For PSR J0952–0607, high irradiation (fpulsar = 0.5) after case A RLO is required to reproduce its estimated orbital properties. For 1FGL J1417.7–440, the decreased expansion of the orbit due to pulsar wind irradiation (fpulsar = 0.1) matches the observations well, though irradiation is not required to reproduce the orbital properties of these diverging NS binaries.

  • The NSs in LMXBs and IMXBs accrete between 0.2 and 1.0 M to spin up to MSP spins. For NSs to attain sub-MSP spins, the mass accreted is large enough to collapse it into a BH (≳1.0 M for M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $). Additionally, the spin at which an accreting NS would enter the sub-MSP regime (and collapse) increases with increasing mass, going from 1.38 ms to 2.66 ms when M NS i $ M^{i}_{\mathrm{NS}} $ goes from 1.3 M to 2.0 M. This could explain why sub-MSPs are not observed and why the fastest spinning BW spins at 1.4 ms. With β = 0.3, we can reproduce the fastest spinning BW currently observed PSR J0952–0607 (observed spin of 1.41 ms and pulsar mass of 2.35 M). Hence, fast spinning NSs are key in identifying the boundary between NSs and BHs.

  • Around 50 to 60% of our simulated spiders ended up becoming dynamically unstable due to mass transfer initiated from a fully convective companion. These binaries would go through a common envelope phase, and the surviving binaries would be isolated spinning NSs. These binaries could explain the observations of isolated MSPs.

Since this study involved running detailed MESA simulations with a fixed set of initial parameters, we did not know the formation probability of these types of systems and how it would compare to the observed populations. In order to calculate the formation rates of spiders, we would need to carry out a population synthesis study. Using large-scale parameter studies, we can better constrain certain physical parameters governing binary evolution. These include processes affecting the evolution prior to the XRB phase (e.g., the supernova mechanism) and more directly involved physical processes (e.g., magnetic braking and pulsar wind irradiation). Since LMXBs have very unstable accretion disks that cause outbursts that would affect pulsar spin-up (for instance, due to X-ray irradiation during RLO; Dubus et al. 1999), unsteady or transient accretion must be further investigated in detail. While we do include a model to account for these effects that would cause non-conservative mass transfer, these instabilities and their effects on the observed populations are highly uncertain and require further investigation. Finally, this study was conducted at solar metallicity because of the substantial population of Galactic spiders that have been observed. However, other metallicities should also be investigated to account for the spiders that have been observed in globular clusters with metallicities below solar (e.g., PSR J1953+1844 in the globular cluster M71 at 1/10th solar metallicity; Pan et al. 2023).

Acknowledgments

The authors thank the anonymous referee for their constructive comments that helped improve the manuscript. The authors thank Thomas Tauris, Hai-Liang Chen, and Raphaël Mignon-Risse for their useful discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101002352, PI: M. Linares). C.S.Y. acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) DIS-2022-568580.

References

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 848 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abolmasov, P., Biryukov, A., & Popov, S. B. 2024, Galaxies, 12, 7 [NASA ADS] [CrossRef] [Google Scholar]
  3. Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Nature, 300, 728 [NASA ADS] [CrossRef] [Google Scholar]
  4. Antoniadis, J. 2014, ApJ, 797, L24 [NASA ADS] [CrossRef] [Google Scholar]
  5. Antoniadis, J., van Kerkwijk, M. H., Koester, D., et al. 2012, MNRAS, 423, 3316 [NASA ADS] [CrossRef] [Google Scholar]
  6. Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448 [Google Scholar]
  7. Antoniadis, J., Tauris, T. M., Ozel, F., et al. 2016, ArXiv e-prints [arXiv:1605.01665] [Google Scholar]
  8. Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Science, 324, 1411 [Google Scholar]
  9. Archibald, A. M., Bogdanov, S., Patruno, A., et al. 2015, ApJ, 807, 62 [NASA ADS] [CrossRef] [Google Scholar]
  10. Archibald, R. F., Gotthelf, E. V., Ferdman, R. D., et al. 2016, ApJ, 819, L16 [NASA ADS] [CrossRef] [Google Scholar]
  11. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bassa, C. G., Patruno, A., Hessels, J. W. T., et al. 2014, MNRAS, 441, 1825 [Google Scholar]
  13. Bassa, C. G., Pleunis, Z., Hessels, J. W. T., et al. 2017, ApJ, 846, L20 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bejger, M., Blaschke, D., Haensel, P., Zdunik, J. L., & Fortin, M. 2017, A&A, 600, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Benvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2014, ApJ, 786, L7 [NASA ADS] [CrossRef] [Google Scholar]
  16. Benvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2015, ApJ, 798, 44 [Google Scholar]
  17. Bhattacharya, D., & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1 [Google Scholar]
  18. Bhattacharyya, S., & Chakrabarty, D. 2017, ApJ, 835, 4 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bildsten, L. 1998, ApJ, 501, L89 [Google Scholar]
  20. Biryukov, A., Beskin, G., & Karpov, S. 2012, MNRAS, 420, 103 [NASA ADS] [CrossRef] [Google Scholar]
  21. Breton, R. P., van Kerkwijk, M. H., Roberts, M. S. E., et al. 2013, ApJ, 769, 108 [Google Scholar]
  22. Camilo, F., Reynolds, J. E., Ransom, S. M., et al. 2016, ApJ, 820, 6 [NASA ADS] [CrossRef] [Google Scholar]
  23. Chattopadhyay, D., Stevenson, S., Hurley, J. R., Rossi, L. J., & Flynn, C. 2020, MNRAS, 494, 1587 [NASA ADS] [CrossRef] [Google Scholar]
  24. Chen, W.-C., & Podsiadlowski, P. 2016, ApJ, 830, 131 [CrossRef] [Google Scholar]
  25. Chen, H.-L., Chen, X., Tauris, T. M., & Han, Z. 2013, ApJ, 775, 27 [NASA ADS] [CrossRef] [Google Scholar]
  26. Chen, H.-L., Tauris, T. M., Han, Z., & Chen, X. 2021, MNRAS, 503, 3540 [NASA ADS] [CrossRef] [Google Scholar]
  27. Chen, H.-L., Tauris, T. M., Chen, X., & Han, Z. 2023, ApJ, 951, 91 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cumming, A., Zweibel, E., & Bildsten, L. 2001, ApJ, 557, 958 [NASA ADS] [CrossRef] [Google Scholar]
  29. de Martino, D., Papitto, A., Belloni, T., et al. 2015, MNRAS, 454, 2190 [NASA ADS] [CrossRef] [Google Scholar]
  30. Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081 [Google Scholar]
  31. Deng, Z.-L., Li, X.-D., Gao, Z.-F., & Shao, Y. 2021, ApJ, 909, 174 [NASA ADS] [CrossRef] [Google Scholar]
  32. Draghis, P., & Romani, R. W. 2018, ApJ, 862, L6 [NASA ADS] [CrossRef] [Google Scholar]
  33. Dubus, G., Lasota, J.-P., Hameury, J.-M., & Charles, P. 1999, MNRAS, 303, 139 [Google Scholar]
  34. Echeveste, M., Novarino, M. L., Benvenuto, O. G., & De Vito, M. A. 2024, MNRAS, 530, 4277 [NASA ADS] [CrossRef] [Google Scholar]
  35. El-Badry, K., Conroy, C., Fuller, J., et al. 2022, MNRAS, 517, 4916 [NASA ADS] [CrossRef] [Google Scholar]
  36. Ergma, E., & Fedorova, A. V. 1992, A&A, 265, 65 [NASA ADS] [Google Scholar]
  37. Ertan, U. 2015, ArXiv e-prints [arXiv:1504.03996] [Google Scholar]
  38. Ertan, Ü. 2018, MNRAS, 479, L12 [NASA ADS] [CrossRef] [Google Scholar]
  39. Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2023, ApJS, 264, 45 [NASA ADS] [CrossRef] [Google Scholar]
  40. Fruchter, A. S., Stinebring, D. R., & Taylor, J. H. 1988, Nature, 333, 237 [Google Scholar]
  41. Gerlach, U. H. 1968, Phys. Rev., 172, 1325 [NASA ADS] [CrossRef] [Google Scholar]
  42. Ginzburg, S., & Quataert, E. 2020, MNRAS, 495, 3656 [Google Scholar]
  43. Ginzburg, S., & Quataert, E. 2021, MNRAS, 500, 1592 [Google Scholar]
  44. Glendenning, N. K., & Kettner, C. 2000, A&A, 353, L9 [NASA ADS] [Google Scholar]
  45. Gossage, S., Kalogera, V., & Sun, M. 2023, ApJ, 950, 27 [NASA ADS] [CrossRef] [Google Scholar]
  46. Hamil, O., Stone, J. R., Urbanec, M., & Urbancová, G. 2015, Phys. Rev. D, 91, 063007 [NASA ADS] [CrossRef] [Google Scholar]
  47. Haskell, B., Zdunik, J. L., Fortin, M., et al. 2018, A&A, 620, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. He, X., Meng, X.-C., & Chen, H.-L. 2019, RAA, 19, 110 [NASA ADS] [Google Scholar]
  49. Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901 [CrossRef] [PubMed] [Google Scholar]
  50. Hjellming, M. S., & Webbink, R. F. 1987, ApJ, 318, 794 [Google Scholar]
  51. Hobbs, G., Lyne, A. G., Kramer, M., Martin, C. E., & Jordan, C. 2004, MNRAS, 353, 1311 [NASA ADS] [CrossRef] [Google Scholar]
  52. Huang, S.-S. 1963, ApJ, 138, 471 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  54. Illarionov, A. F., & Sunyaev, R. A. 1975, A&A, 39, 185 [NASA ADS] [Google Scholar]
  55. Istrate, A. G., Tauris, T. M., & Langer, N. 2014, A&A, 571, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Jacoby, B. A., Hotan, A., Bailes, M., Ord, S., & Kulkarni, S. R. 2005, ApJ, 629, L113 [Google Scholar]
  57. Jeans, J. H. 1924, MNRAS, 85, 2 [NASA ADS] [CrossRef] [Google Scholar]
  58. Jia, K., & Li, X. D. 2014, ApJ, 791, 127 [NASA ADS] [CrossRef] [Google Scholar]
  59. Jia, K., & Li, X.-D. 2015, ApJ, 814, 74 [NASA ADS] [CrossRef] [Google Scholar]
  60. Johnston, S., & Karastergiou, A. 2017, MNRAS, 467, 3493 [CrossRef] [Google Scholar]
  61. Kaplan, D. L., Bhalerao, V. B., van Kerkwijk, M. H., et al. 2013, ApJ, 765, 158 [NASA ADS] [CrossRef] [Google Scholar]
  62. Kar, A., Ojha, P., & Bhattacharyya, S. 2024, MNRAS, 535, 344 [NASA ADS] [CrossRef] [Google Scholar]
  63. Kiel, P. D., Hurley, J. R., Bailes, M., & Murray, J. R. 2008, MNRAS, 388, 393 [NASA ADS] [CrossRef] [Google Scholar]
  64. Kluzniak, W., Ruderman, M., Shaham, J., & Tavani, M. 1988, Nature, 334, 225 [NASA ADS] [CrossRef] [Google Scholar]
  65. Konar, S., & Bhattacharya, D. 1997, MNRAS, 284, 311 [NASA ADS] [CrossRef] [Google Scholar]
  66. Konar, S., & Bhattacharya, D. 1999a, MNRAS, 303, 588 [NASA ADS] [CrossRef] [Google Scholar]
  67. Konar, S., & Bhattacharya, D. 1999b, MNRAS, 308, 795 [NASA ADS] [CrossRef] [Google Scholar]
  68. Kong, A. K. H., Jin, R., Yen, T. C., et al. 2014, ApJ, 794, L22 [NASA ADS] [CrossRef] [Google Scholar]
  69. Kumari, S., Bhattacharyya, B., Sharan, R., et al. 2024, ApJ, 961, 155 [NASA ADS] [CrossRef] [Google Scholar]
  70. Lasota, J.-P. 2001, New Astron. Rev., 45, 449 [Google Scholar]
  71. Li, L., Han, Z., & Zhang, F. 2004, MNRAS, 355, 1383 [NASA ADS] [CrossRef] [Google Scholar]
  72. Linares, M. 2020, Multifrequency Behaviour of High Energy Cosmic Sources – XIII, 3–8 June 2019, Palermo, 23 [Google Scholar]
  73. Linares, M., & Kachelrieß, M. 2021, JCAP, 2021, 030 [CrossRef] [Google Scholar]
  74. Linares, M., Shahbaz, T., & Casares, J. 2018, ApJ, 859, 54 [Google Scholar]
  75. Liu, X.-J., You, Z.-Q., Chen, Z.-C., et al. 2024, ApJ, 962, 80 [NASA ADS] [CrossRef] [Google Scholar]
  76. Ma, B., & Li, X.-D. 2009, ApJ, 691, 1611 [Google Scholar]
  77. Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005a, AJ, 129, 1993 [Google Scholar]
  78. Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005b, VizieR Online Data Catalog: VII/245 [Google Scholar]
  79. Mineshige, S. 1993, Ap&SS, 210, 83 [NASA ADS] [CrossRef] [Google Scholar]
  80. Misra, D., Fragos, T., Tauris, T. M., Zapartas, E., & Aguilera-Dena, D. R. 2020, A&A, 642, A174 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Nedreaas, I. F. 2024, Master’s Thesis, Norwegian University of Science and Technology, Norway [Google Scholar]
  82. Nelemans, G., & Jonker, P. G. 2010, New Astron. Rev., 54, 87 [Google Scholar]
  83. Nelemans, G., Jonker, P. G., & Steeghs, D. 2006, MNRAS, 370, 255 [Google Scholar]
  84. Nelson, L. A., Rappaport, S. A., & Joss, P. C. 1986, ApJ, 304, 231 [Google Scholar]
  85. Osłowski, S., Bulik, T., Gondek-Rosińska, D., & Belczyński, K. 2011, MNRAS, 413, 461 [CrossRef] [Google Scholar]
  86. Pallanca, C., Mignani, R. P., Dalessandro, E., et al. 2012, ApJ, 755, 180 [NASA ADS] [CrossRef] [Google Scholar]
  87. Pan, Z., Lu, J. G., Jiang, P., et al. 2023, Nature, 620, 961 [CrossRef] [Google Scholar]
  88. Papaloizou, J., & Pringle, J. E. 1978, MNRAS, 184, 501 [NASA ADS] [CrossRef] [Google Scholar]
  89. Papitto, A., & de Martino, D. 2022, in Astrophysics and Space Science Library, eds. S. Bhattacharyya, A. Papitto, & D. Bhattacharya, 465, 157 [NASA ADS] [CrossRef] [Google Scholar]
  90. Papitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517 [NASA ADS] [CrossRef] [Google Scholar]
  91. Papitto, A., de Martino, D., Belloni, T. M., et al. 2015, MNRAS, 449, L26 [NASA ADS] [CrossRef] [Google Scholar]
  92. Patruno, A., Haskell, B., & D’Angelo, C. 2012, ApJ, 746, 9 [NASA ADS] [CrossRef] [Google Scholar]
  93. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  94. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  95. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  96. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  97. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  98. Peters, P. C. 1964, Phys. Rev., 136, 1224 [Google Scholar]
  99. Pletsch, H. J., Guillemot, L., Fehrmann, H., et al. 2012, Science, 338, 1314 [CrossRef] [Google Scholar]
  100. Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107 [NASA ADS] [CrossRef] [Google Scholar]
  101. Podsiadlowski, P., Han, Z., & Rappaport, S. 2003, MNRAS, 340, 1214 [NASA ADS] [CrossRef] [Google Scholar]
  102. Polzin, E. J., Breton, R. P., Bhattacharyya, B., et al. 2020, MNRAS, 494, 2948 [Google Scholar]
  103. Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
  104. Pylyser, E. H. P., & Savonije, G. J. 1989, A&A, 208, 52 [Google Scholar]
  105. Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, ApJ, 918, L29 [NASA ADS] [CrossRef] [Google Scholar]
  106. Radhakrishnan, V., & Srinivasan, G. 1982, Curr. Sci., 51, 1096 [NASA ADS] [Google Scholar]
  107. Rappaport, S., Verbunt, F., & Joss, P. C. 1983, ApJ, 275, 713 [Google Scholar]
  108. Ritter, H., & King, A. R. 2001, in Evolution of Binary and Multiple Star Systems, eds. P. Podsiadlowski, S. Rappaport, A. R. King, F. D’Antona, & L. Burderi, ASP Conf. Ser., 229, 423 [Google Scholar]
  109. Roberts, M. S. E. 2013, in Neutron Stars and Pulsars: Challenges and Opportunities after 80 Years, ed. J. van Leeuwen, 291, 127 [NASA ADS] [Google Scholar]
  110. Romani, R. W., Filippenko, A. V., Silverman, J. M., et al. 2012, ApJ, 760, L36 [NASA ADS] [CrossRef] [Google Scholar]
  111. Romani, R. W., Filippenko, A. V., & Cenko, S. B. 2014, ApJ, 793, L20 [NASA ADS] [CrossRef] [Google Scholar]
  112. Romani, R. W., Graham, M. L., Filippenko, A. V., & Zheng, W. 2016, ApJ, 833, 138 [NASA ADS] [CrossRef] [Google Scholar]
  113. Romani, R. W., Kandel, D., Filippenko, A. V., Brink, T. G., & Zheng, W. 2022, ApJ, 934, L17 [NASA ADS] [CrossRef] [Google Scholar]
  114. Roy, J., Ray, P. S., Bhattacharyya, B., et al. 2015, ApJ, 800, L12 [Google Scholar]
  115. Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51 [Google Scholar]
  116. Ruderman, M., Shaham, J., & Tavani, M. 1989a, ApJ, 336, 507 [NASA ADS] [CrossRef] [Google Scholar]
  117. Ruderman, M., Shaham, J., Tavani, M., & Eichler, D. 1989b, ApJ, 343, 292 [NASA ADS] [CrossRef] [Google Scholar]
  118. Shahbaz, T., Linares, M., Rodríguez-Gil, P., & Casares, J. 2019, MNRAS, 488, 198 [NASA ADS] [CrossRef] [Google Scholar]
  119. Shannon, R. M., Cordes, J. M., Metcalfe, T. S., et al. 2013, ApJ, 766, 5 [NASA ADS] [CrossRef] [Google Scholar]
  120. Skumanich, A. 1972, ApJ, 171, 565 [Google Scholar]
  121. Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, A&A, 327, 620 [NASA ADS] [Google Scholar]
  122. Spitkovsky, A. 2006, ApJ, 648, L51 [Google Scholar]
  123. Spruit, H. C., & Ritter, H. 1983, A&A, 124, 267 [NASA ADS] [Google Scholar]
  124. Stappers, B. W., Bailes, M., Lyne, A. G., et al. 1996a, ApJ, 465, L119 [NASA ADS] [CrossRef] [Google Scholar]
  125. Stappers, B. W., Bessell, M. S., & Bailes, M. 1996b, ApJ, 473, L119 [NASA ADS] [CrossRef] [Google Scholar]
  126. Stappers, B. W., van Kerkwijk, M. H., Lane, B., & Kulkarni, S. R. 1999, ApJ, 510, L45 [NASA ADS] [CrossRef] [Google Scholar]
  127. Stevens, I. R., Rees, M. J., & Podsiadlowski, P. 1992, MNRAS, 254, 19P [CrossRef] [Google Scholar]
  128. Stovall, K., Lynch, R. S., Ransom, S. M., et al. 2014, ApJ, 791, 67 [CrossRef] [Google Scholar]
  129. Strader, J., Chomiuk, L., Cheung, C. C., et al. 2015, ApJ, 804, L12 [NASA ADS] [CrossRef] [Google Scholar]
  130. Strader, J., Swihart, S., Chomiuk, L., et al. 2019, ApJ, 872, 42 [Google Scholar]
  131. Sturrock, P. A. 1971, ApJ, 164, 529 [Google Scholar]
  132. Swihart, S. J., Strader, J., Johnson, T. J., et al. 2017, ApJ, 851, 31 [NASA ADS] [CrossRef] [Google Scholar]
  133. Swihart, S. J., Strader, J., Shishkovsky, L., et al. 2018, ApJ, 866, 83 [NASA ADS] [CrossRef] [Google Scholar]
  134. Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928 [NASA ADS] [Google Scholar]
  135. Tauris, T. M., & van den Heuvel, E. P. J. 2006, Compact Stellar X-ray Sources, 39, 623 [NASA ADS] [CrossRef] [Google Scholar]
  136. Tauris, T. M., Langer, N., & Kramer, M. 2012, MNRAS, 425, 1601 [NASA ADS] [CrossRef] [Google Scholar]
  137. Thorstensen, J. R., & Armstrong, E. 2005, AJ, 130, 759 [NASA ADS] [CrossRef] [Google Scholar]
  138. Turchetta, M., Linares, M., Koljonen, K., & Sen, B. 2023, MNRAS, 525, 2565 [CrossRef] [Google Scholar]
  139. van den Heuvel, E. P. J., & van Paradijs, J. 1988, Nature, 334, 227 [NASA ADS] [CrossRef] [Google Scholar]
  140. van der Sluys, M. V., Verbunt, F., & Pols, O. R. 2005a, A&A, 440, 973 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. van der Sluys, M. V., Verbunt, F., & Pols, O. R. 2005b, A&A, 431, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  142. van Paradijs, J. 1996, ApJ, 464, L139 [NASA ADS] [CrossRef] [Google Scholar]
  143. Van, K. X., & Ivanova, N. 2019, ApJ, 886, L31 [NASA ADS] [CrossRef] [Google Scholar]
  144. Verbunt, F., & Phinney, E. S. 1995, A&A, 296, 709 [NASA ADS] [Google Scholar]
  145. Verbunt, F., & Zwaan, C. 1981, A&A, 100, L7 [NASA ADS] [Google Scholar]
  146. Wagoner, R. V. 1984, ApJ, 278, 345 [CrossRef] [Google Scholar]
  147. Wang, Z., Archibald, A. M., Thorstensen, J. R., et al. 2009, ApJ, 703, 2017 [NASA ADS] [CrossRef] [Google Scholar]
  148. Wang, B., Chen, W.-C., Liu, D.-D., et al. 2021, MNRAS, 506, 4654 [NASA ADS] [CrossRef] [Google Scholar]
  149. Wijnands, R., & van der Klis, M. 1998, Nature, 394, 344 [CrossRef] [Google Scholar]
  150. Woan, G., Pitkin, M. D., Haskell, B., Jones, D. I., & Lasky, P. D. 2018, ApJ, 863, L40 [Google Scholar]
  151. Ye, C. S., Kremer, K., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2019, ApJ, 877, 122 [NASA ADS] [CrossRef] [Google Scholar]
  152. Zhang, C. M. 1998, A&A, 330, 195 [NASA ADS] [Google Scholar]
  153. Zhang, C. M., & Kojima, Y. 2006, MNRAS, 366, 137 [NASA ADS] [CrossRef] [Google Scholar]
  154. Zhang, S.-N., & Xie, Y. 2012, ApJ, 761, 102 [NASA ADS] [CrossRef] [Google Scholar]
  155. Zhang, C. M., Wang, J., Zhao, Y. H., et al. 2011, A&A, 527, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Appendix A: Numerical stability of MESA tracks

Fragos et al. (2023, Section 5.2) list a wide range of criteria for successful termination of the MESA track, (1) the age of the companion to the NS exceeds the Hubble age (13.8 Gyr), (2) the companion to the NS becomes a WD, and (3) the binary enters a common envelope. The binaries that did not satisfy the criteria for successful termination (mentioned in Section 2.1), were rerun with a limitation on the maximum radiative opacity (κmax) of the companion. The initial failure rate was 16.5% in the entire initial simulated grid. Since the stellar opacity influences the internal stellar structure, it prevents the star properties from fluctuating too drastically for MESA to handle. The rerun κmax value used by Fragos et al. (2023) is 0.5 cm2 g−1, which brought down our failure rate to 0.27%. However, we find this value to be too restrictive for some stars in our initial parameter space, resulting in a significant deviation in the resulting tracks from the successful tracks in close vicinity. This can be seen in Figure A.1, where we show a binary (with β = 0.0) having an initial companion mass of 4.5 M, an initial NS mass of 2.0 M, and an initial orbital period of 1.2 d, which failed due to numerical instability. We also show a neighboring binary that avoided numerical instabilities with 4.5 M, initial NS mass of 2.0 M, and initial orbital period of 1.4 d. Decreasing the value of κmax to 0.5 results in the stellar properties deviating significantly from previous reruns, from the very beginning of the evolution. The binary that previously should have undergone a dynamical instability with a strict κmax does not have a high mass-transfer rate and ends as a WD. Hence, to keep the tracks as similar to the original run as possible (while avoiding numerical instabilities) as well as similar in behavior to their close neighbors, we reran the failed binaries in two batches. First we used a κmax of 10 cm2 g−1, and then we ran the binaries that failed again with a κmax of 5 cm2 g−1. Our final failure rate was 9.19%, which is a lot higher than the case of a strict κmax; however, we compromised by checking the final properties of the binaries, and if they seemed to have an MSP-like spin for the NS (≤30 ms), we considered the binary successful.

thumbnail Fig. A.1.

Stellar luminosity (L) versus effective temperature (Teff) of a binary from the simulated grids with an initial companion mass of 4.5 M, initial NS mass of 2.0 M, and initial orbital period of 1.4 d for different maximum opacity values (κmax). The solid blue curve is the failed binary. The solid green, dashed red, and dashed purple curves are with κmax value of 10.0, 5.0, and 0.5, respectively. We compare the tracks to a neighboring track (without any numerical instability) in solid orange.

Appendix B: Low accretion efficiency grid

We present the evolutionary grids with β = 0.7 in Figure B.1. The figure also shows the observed population of spiders, BWs, RBs, and huntsmen. As can be seen from the figures, most of the tracks miss the spider populations, and hence we look to more conservative mass-accretion efficiencies for our study.

thumbnail Fig. B.1.

Orbital period evolution versus companion mass of the simulated binaries for β = 0.7 and with an initial NS of mass 1.3 M and fpulsar values of 0.0 (top left), 0.05 (top right), 0.1 (bottom left), and 0.5 (bottom right). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. We also compare our tracks to estimated Mc, min and observed Porb for spiders (Section 2.4). The initial binaries have a MS star (range of 1.0 to 4.5 M) and initial periods of 0.6 to 4.0 d.

All Tables

Table 1.

Observed and estimated properties of the two observed huntsmen spiders.

All Figures

thumbnail Fig. 1.

Orbital period evolution versus companion mass of the simulated binaries for β = 0.3 with an initial NS of mass 1.3 M and fpulsar values of 0.0 (top left), 0.05 (top right), 0.1 (bottom left), and 0.5 (bottom right). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. The color bar shows the changing surface abundance of the companion. We also compare our tracks to estimated Mc, min and observed Porb for BWs (black crosses) and RBs (red crosses), and huntsmen spiders (green crosses), including whether H-lines have been detected (filled circles) or not (filled triangles) or if it is unknown (crosses), from Nedreaas (2024, and references therein). The green star is the huntsman candidate 2FGL J0846.0+2820, observed by Fermi-LAT (Swihart et al. 2017) (surface H is unknown). The three tidarrens (J1311–3430, J1653–0159, and J0636+5128) are also shown. The initial binaries have a MS star (range of 1.0 to 4.5 M) and initial periods of 0.6 to 4.0 d.

In the text
thumbnail Fig. 2.

Evolution of the NS spins (Pspin) versus the orbital periods (Porb) in converging binaries for β = 0.0 (top), β = 0.3 (middle), and β = 0.7 (bottom). We distinguish the binaries if the final NS mass was above 2.5 M (faint gray curves) or below it, which are further grouped by two initial NS masses ( M NS i $ M^{i}_{\mathrm{NS}} $) 1.3 M (blue curves) and 2.0 M (orange curves). We also compare to the observed spins of BWs and RBs, shown by the black and red circles, respectively.

In the text
thumbnail Fig. 3.

Evolution of selected binaries from our simulated grids for β = 0.3, fpular = 0.0, and M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $. The three different kinds of binaries that we discuss in the text are studied here: case A ( P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $ d), near-Pbif ( P orb i = 2.4 $ P^{i}_{\mathrm{orb}} = 2.4 $ d), and case B or diverging ( P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d). Also shown is the binary that goes down to the orbital period of 19 min (orange curve) for comparison, which is the minimum in the grid with fpulsar = 0.0. Top left shows the orbital period versus companion mass, where Mc is the mass of the companion to the NS and Porb is the orbital period. The formation of a fully convective star and the start of the ejection phase are shown as squares. We also compare our tracks to estimated Mc, min and observed Porb for BWs (faint black crosses), RBs (faint red crosses), and huntsmen spiders (faint green crosses), from Nedreaas (2024, and references therein). The green star is the huntsman candidate 2FGL J0846.0+2820, observed by Fermi-LAT (Swihart et al. 2017). The tidarrens, namely, J1311–3430, J1653–0159, and J0636+5128, are highlighted as pink, blue, and black crosses, respectively. Top right, bottom left, and bottom right show the time evolution of the NS spin, companion mass loss rate, and NS mass, respectively, for the same binaries. The time axis is zoomed in to focus on the spin-up phase.

In the text
thumbnail Fig. 4.

Allowed initial parameter space to form spider MSPs for an initial NS of mass 1.3 M and fpulsar values of 0.0 (left) and 0.5 (right). The term M c i $ M^{i}_{\mathrm{c}} $ is the initial mass of the companion to the NS, and P orb i $ P^{i}_{\mathrm{orb}} $ is the initial orbital period. The color bar shows the final orbital periods ( P orb f $ P^{f}_{\mathrm{orb}} $) and the gray symbols did not attain MSP-like spins. We show the end states of the tracks, if they encountered a dynamical instability (with gray triangles or with colorful triangles depending on if onset of instability occurred during the first RLO phase or the second, respectively) or ended with one of the criteria for successful termination (star symbols). The diverging symbols are shown by gray crosses. We also show the near-Pbif binaries that passed through the RB part of the observed parameter space enclosed by dashed orange lines. The magenta triangle shows the initial parameters of the binary that reproduces PSR J0952–0607 (with observed spin period of 1.41 ms).

In the text
thumbnail Fig. 5.

Evolution of the mass-loss rates (c, black curves) and the mass-transfer rates (transfer, blue curves) of the companion star. The initial binary masses are M c i = 1.0 M $ M^{i}_{\mathrm{c}} = 1.0\,\mathrm{M}_\odot $ and M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, with β = 0.3. We present the three cases of interacting binaries studied, case A (top row; with P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $ d), near-Pbif (middle row; with P orb i = 2.4 $ P^{i}_{\mathrm{orb}} = 2.4 $ d), and case B (bottom row; with P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d), each for three values of fpulsar = 0.05 (left column), 0.1 (middle column), and 0.5 (right column).

In the text
thumbnail Fig. 6.

Orbital period evolution versus companion mass of selected binaries from our simulated grid shown in Figure 1. The term Mc is the mass of the companion to the NS, and Porb is the orbital period. The figure shows the effect of pulsar wind irradiation, comparing fpulsar = 0.0 (dashed lines) and 0.5 (solid lines). The left panel shows the tracks with M c i = 1.5 M $ M^{i}_{\mathrm{c}} = 1.5\,\mathrm{M}_\odot $, P orb i = 1.0 $ P^{i}_{\mathrm{orb}} = 1.0 $, 2.4, and 2.6 d. The right panel shows two initial companion masses, M c i = 1.0 $ M^{i}_{\mathrm{c}} = 1.0 $ (with P orb i = 1 $ P^{i}_{\mathrm{orb}} = 1 $ d) and 2.5 M ( P orb i = 2.0 $ P^{i}_{\mathrm{orb}} = 2.0 $ and 2.2 d). We compared our tracks to Mc, min and Porb for BWs (faint black crosses), RBs (faint red crosses), and huntsmen spiders (faint green crosses), from Nedreaas (2024, and references therein). The green star is the huntsman candidate 2FGL J0846.0+2820, observed by Fermi-LAT (Swihart et al. 2017). The magenta stars show the three transitional MSPs (tMSPs) observed, PSR J1023+0038 (Mc = 0.22 M and Porb = 4.75 h; Shahbaz et al. 2019; Thorstensen & Armstrong 2005), IGR J18245–2452 (Mc, min = 0.17 M and Porb = 11.03 h; Papitto et al. 2013), and XSS J12270–4859 (Mc, min = 0.15 M and Porb = 6.91 h; Roy et al. 2015; de Martino et al. 2015; Bassa et al. 2014).

In the text
thumbnail Fig. 7.

Orbital period versus companion mass for diverging binary tracks from our simulated grids that have a similar NS spin as the observed PSR J0952–0607 (1.41 ms) within a 10% margin (dashed gray curves). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. All the tracks shown have β = 0.3 and fpulsar = 0.5. We further show the binary that best matches the sources in its orbital evolution in blue curve, which has initial companion mass of 2.0 M, initial NS mass of 1.3 M, and initial orbital period of 1.0 d.

In the text
thumbnail Fig. 8.

Final NS spins ( P spin f $ P^{f}_{\mathrm{spin}} $) versus final NS masses ( M NS f $ M^{f}_{\mathrm{NS}} $) in diverging binaries for β = 0.0 (blue colored symbols), β = 0.3 (orange colored symbols), and β = 0.7 (green colored symbols). We group the binaries by the initial NS mass ( M NS i $ M^{i}_{\mathrm{NS}} $) corresponding to different symbols (see legend). We also compare to the observed spin and estimated pulsar mass of the huntsman source 1FGL J1417.7–440 (2.66 ms and 1 . 62 0.17 + 0.43 M $ 1.62^{+0.43}_{-0.17}\,\mathrm{M}_\odot $, shown by the gray dashed lines and gray shaded area for the uncertainties; Camilo et al. 2016; Swihart et al. 2018) and the estimated pulsar mass of the huntsman candidate 2FGL J0846.0+2820 ( 1 . 96 0.41 + 0.41 M $ 1.96^{+0.41}_{-0.41}\,\mathrm{M}_\odot $, shown by the vertical pink dashed lines and pink shaded area for the uncertainties; Swihart et al. 2017).

In the text
thumbnail Fig. 9.

Allowed initial parameter space to form huntsmen spiders (fpulsar = 0.0) for an initial NS of mass 1.3 M (left) and 2.0 M (right). The term M c i $ M^{i}_{\mathrm{c}} $ is the initial mass of the companion to the NS, and P orb i $ P^{i}_{\mathrm{orb}} $ is the initial orbital period. The color bar shows the final NS spins ( P spin f $ P^{f}_{\mathrm{spin}} $), and the gray stars denote the systems that formed diverging orbits but did not spin up to the MSP range. We show the end states of the tracks, if they encountered a dynamical instability (triangle symbols) or ended with one of the criteria for successful termination (star symbols). Converging tracks are shown by gray crosses.

In the text
thumbnail Fig. 10.

Orbital period versus companion mass for diverging binary tracks from our simulated grids that have a similar NS spin as the observed huntsman source 1FGL J1417.7–440 (2.66 ms) within a 20% margin (dashed gray curves). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. We further show the binary that best matches the sources in its orbital evolution in blue curves, where β = 0.0, fpulsar = 0.1, M c i = 2.0 M $ M^{i}_{\mathrm{c}} = 2.0\,\mathrm{M}_\odot $, M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $, and P orb i = 2.6 $ P^{i}_{\mathrm{orb}} = 2.6 $ d.

In the text
thumbnail Fig. 11.

Accreted mass versus final NS spins ( P spin f $ P^{f}_{\mathrm{spin}} $) in our simulated grid with β = 0.3. We show the different initial masses as different colors (shown in the legend) and distinguish between the NSs that had M NS f 2.5 M $ M^{f}_{\mathrm{NS}}\ge 2.5\,\mathrm{M}_\odot $ (crosses) and those where M NS f < 2.5 M $ M^{f}_{\mathrm{NS}} < 2.5\,\mathrm{M}_\odot $ (circles). We also compare our results to the equilibrium spin period calculations from Tauris et al. (2012), done for M NS i = 1.3 M $ M^{i}_{\mathrm{NS}} = 1.3\,\mathrm{M}_\odot $ (light blue curve) and 2.0 M (dark orange curve). The black dotted line shows the observed spin for PSR J0952–0607 (spin period of 1.41 ms) and the gray dotted line shows the upper limit for MSPs considered (30 ms).

In the text
thumbnail Fig. A.1.

Stellar luminosity (L) versus effective temperature (Teff) of a binary from the simulated grids with an initial companion mass of 4.5 M, initial NS mass of 2.0 M, and initial orbital period of 1.4 d for different maximum opacity values (κmax). The solid blue curve is the failed binary. The solid green, dashed red, and dashed purple curves are with κmax value of 10.0, 5.0, and 0.5, respectively. We compare the tracks to a neighboring track (without any numerical instability) in solid orange.

In the text
thumbnail Fig. B.1.

Orbital period evolution versus companion mass of the simulated binaries for β = 0.7 and with an initial NS of mass 1.3 M and fpulsar values of 0.0 (top left), 0.05 (top right), 0.1 (bottom left), and 0.5 (bottom right). The term Mc is the mass of the companion to the NS, and Porb is the orbital period. We also compare our tracks to estimated Mc, min and observed Porb for spiders (Section 2.4). The initial binaries have a MS star (range of 1.0 to 4.5 M) and initial periods of 0.6 to 4.0 d.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.