Compact binary millisecond pulsars MSPs with orbital periods lesssim1d are key to understanding binary evolution involving massive neutron stars NSs. Due to ablation of the companion by the rapidly spinning pulsar these systems are also known as spiders categorized into two main branches redbacks RBs companion mass in the range of 0.1 to 0.5 Modot and black widows BWs companion mass lesssim 0.1 Modot. We present models of low and intermediatemass Xray binaries XRBs and compare them with observations of Galactic spiders including the presenceabsence of hydrogen lines in their optical spectra. We constrain and quantify the interaction between pulsar and companion. Using MESA we create the allowed initial parameter space modeling the irradiation of the companion by the pulsar wind. For the first time in MESA we include the detailed evolution of the pulsar spin. Efficient mass accretion onto the NS 70 of mass transferred with an Xray irradiated disc followed by strong irradiation of the companion can explain most properties of the observed spiders. Our RBs continue to the BW regime connecting the two branches of spiders. Our models explain the lack of hydrogen in some observed BWs with ultralight companions. The mass required to spinup a NS to submilliseconds is high enough to collapse it into a black hole. Finally we discuss the formation of RBlike spiders with giant companions and orbital periods of several days huntsmen showing that they are unlikely to produce supermassive NSs maximum accreted mass lesssim0.5 Modot. Cannibalistic MSP binary formation depends heavily on the interplay between accretion onto the pulsar and pulsar wind irradiation. Our work supports earlier claims that RBs evolve into BWs. We also show that the fastestspinning pulsars may collapse before reaching submillisecond spin periods.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Solar and Stellar Astrophysics']
Investigating cannibalistic millisecond pulsar binaries using MESA new constraints from pulsar spin and mass evolution
{'Investigating cannibalistic millisecond pulsar binaries using MESA : new constraints from pulsar spin and mass evolution': 'Devina Misra 1 ,⋆ , Manuel Linares 1 , 2 , Claire S. Ye 3  \n- 1 Department of Physics, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway\n- 2 Departament de Física, EEBE, Universitat Politècnica de Catalunya, Av. Eduard Maristany 16, E-08019 Barcelona, Spain.\n- 3 Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, Ontario M5S 3H8, Canada  \nReceived XXYYZZ; Accepted XXYYZZ', 'ABSTRACT': 'Context. Compact binary millisecond pulsars (MSPs), with orbital periods ≲ 1d, are key to understanding binary evolution involving massive neutron stars (NSs). Due to the ablation of the companion by the rapidly spinning pulsar, these systems are also known as spiders, categorized into two main branches: redbacks (RBs; companion mass in the range of 0.1 to 0.5 M ⊙ ) and black widows (BWs; companion mass ≲ 0.1 M ⊙ ).  \nAims. We present models of low- and intermediate-mass X-ray binaries (XRBs) and compare them with observations of Galactic spiders (including the presence / absence of hydrogen lines in their optical spectra). We constrain and quantify the interaction between pulsar and companion.  \nMethods. Using MESA , we create the allowed initial parameter space, modeling the irradiation of the companion by the pulsar wind. For the first time in MESA , we also include the detailed evolution of the pulsar spin.  \nResults. E ffi cient mass accretion onto the NS (70% of the mass transferred) with an X-ray irradiated disc, followed by strong irradiation of the companion can explain most properties of the observed spiders. Our RB evolutionary tracks continue to the BW regime, connecting the two branches of spiders. Our models explain the lack of hydrogen in some observed BWs with ultra-light companions. The mass required to spin-up a NS to sub-milliseconds is high enough to collapse it into a black hole. Finally, we discuss the formation of RB-like spiders with giant companions and orbital periods of several days (huntsmen), concluding that they are unlikely to produce super-massive NSs (maximum accreted mass ≲ 0.5 M ⊙ ).  \nConclusions. Cannibalistic MSP binary formation depends heavily on the interplay between accretion onto the pulsar and pulsar wind irradiation. Our work supports earlier claims that RBs evolve into BWs. We also show that the fastest spinning pulsars may collapse before reaching sub-millisecond spin periods.  \nKey words. accretion, accretion disks - binaries: eclipsing - close - stars: neutron - pulsars: general - methods: numerical', '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 ( ∼ 10 8 G). The neutron stars (NS) in LMXBs accrete matter which 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 is also suggested to decrease with 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 the Roche-lobe overflow (RLO) has ended, the highly spinning pulsar loses energy in the form of energetic particles that irradiates 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  \nvery 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 irradiationinduced mass loss has since been observed in multiple eclipsing compact MSPs (Kluzniak et al. 1988; Ruderman et al. 1989; Stappers et al. 1996b; Polzin et al. 2020; Kumari et al. 2024).  \nDue to the cannibalistic nature of this interaction, these compact binary MSPs, with observed orbital periods P orb ≲ 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 our 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 M c): Redbacks (RBs; having companions with masses 0 . 1 M ⊙ ≲ M c ≲ 0 . 5 M ⊙ ) and black widows (BWs; having companions with masses M c < 0 . 1 M ⊙ ). The orbital period distribution is similar for both types ( ≲ 1 d).  \nFollowing classical binary evolution, the observed orbital periods of all spiders are not 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 e ff ect on the orbit (Soberman et al. 1997). Therefore, the involvement of this energetic radio  \npulsar 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 found that in order to reproduce the orbital parameters of spiders, cyclic mass transfer due to the e ff ect of X-rays is required (Benvenuto et al. 2014, 2015). Ginzburg & Quataert (2020) and Ginzburg & Quataert (2021) argued that sustained magnetic braking (that couples to stellar wind leading to angular momentum loss) is needed to reproduce spiders and ablation alone is not su ffi cient. Some spiders are also thought to come from ultra-compact X-ray binaries (UCXBs), which are a subcategory of LMXBs with hydrogen-deficient companions and orbital periods of less than an hour (see review by Nelemans &Jonker 2010).  \nSince spider binaries have passed through a long masstransfer 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 NS. 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, we investigate here the fastest spins possible for pulsars. Despite relatively large mass uncertainties, NSs in spider binaries tend towards higher masses (median M NS ∼ 1 . 8 M ⊙ ; Linares 2020) in comparison to other systems like double NSs or wide NS + white dwarf (WD) binaries (Demorest et al. 2010; Antoniadis et al. 2013; Linares 2020).  \nIn 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 e ff ect of various assumptions of binary interaction on the pulsar spin evolution and to constrain these 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, followed by the conclusions of our study in Section 5.', '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  \nof 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. (2022), by identifying failed binaries and re-running 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 10 12 G. These values correspond to the average spin parameters of observed young pulsars. We model the pulsar spin evolution by identifying di ff erent 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 a ff ected 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 ( ˙ Ω ),  \n˙ Ω = -2 µ 2 sin 2 α 3 µ 0 c 2 I Ω 3 , (1)  \nwhere µ is the pulsar magnetic moment ( µ ≈ BR 3 ), α 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 2 NS / 5 for a solid sphere, where M NS and R NS 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 di ff erence in ˙ Ω , at most (Spitkovsky 2006), hence, we do not vary this parameter in our study and fix it at an intermediate value.  \nWe evolve the isolated pulsar surface magnetic field ( B ) decay following Chattopadhyay et al. (2020),  \nB = ( B i -B min) × exp ( -t /τ d) + B min , (2)  \nwhere B i is the initial surface magnetic field, B min is the minimum magnetic field (we take B min = 10 8 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 ≳ 2 Gyr being favoured (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 disc and use the criteria by  \nDubus et al. (1999) that describes the limiting mass-transfer rate below which thermal-viscous instabilities would prevent accretion. This limit is as follows,  \n˙ M crit = 3 . 2 × 10 -9 GLYPH<18> M NS 1 . 4 M ⊙ GLYPH<19> 1 / 2 GLYPH<18> M c 1 M ⊙ GLYPH<19> -1 / 5 GLYPH<18> P orb 1 d GLYPH<19> 7 / 5 M ⊙ yr -1 . (3)  \nThe change in the NS mass ( ˙ M 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 β a ff ects the amount of accretion, lower values of β result in faster spin-up of the pulsar compared to higher values. Additionally, the value of β a ff ects the orbital evolution following the equation for the relative change in orbital separation (˙ a / a ; Tauris & van den Heuvel 2006), assuming other changes in orbital angular momentum are represented by ˙ J orb,  \n˙ a a = ˙ J orb J orb -2 GLYPH<18> 1 + ( β -1) M c M NS -β M c 2( M c + M NS) GLYPH<19> ˙ M c M c , (4)  \nwhere, M c and M NS are companion and NS masses, respectively. Other processes to determine ˙ J orb (apart from mass transfer / loss) are magnetic braking (Rappaport et al. 1983), gravitational wave radiation (Peters 1964), and spin-orbit coupling.  \nTheoretical 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 e ff ect (Illarionov & Sunyaev 1975) and thermal-viscous instabilities in the accretion disc (Pringle 1981; Mineshige 1993; van Paradijs 1996); all e ff ects 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 disc instabilities (explained further below in Section 2.3) and explore three values of β (0.0, 0.3, and 0.7) to study its e ff ect 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 disc instabilities, since these have large uncertainties and there might still be ine ffi cient transfer of angular momentum due to some unforeseen physical process.  \nThe 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,  \nB = ( B i -B min) × exp ( -∆ M NS / M d) + B min . (5)  \nIn the above equation, ∆ M NS is the accreted matter and M d is the magnetic field mass decay scale (0.025, although other values for M d 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,  \n˙ J = κ ˙ M NS R 2 NS ω di ff , (6)  \nwhere κ is an e ffi ciency factor and ω di ff is the angular velocity of the incoming material (Kiel et al. 2008; Chattopadhyay et al. 2020) that is defined by,  \nω di ff = s GM NS r 3 mag -s GM NS r 3 cor . (7)  \nIn the above equation, r cor is the co-rotation radius, which is defined as r cor = ( GM NS / Ω 2 ) 1 / 3 and r mag is the magnetospheric radius, which is defined as follows (Kiel et al. 2008; Chattopadhyay et al. 2020),  \nr mag = r Alfv'en 2 = 2 π 2 G µ 2 0 ! 1 / 7 R 6 NS ˙ M NS M 1 / 2 NS ! 2 / 7 B 4 / 7 , (8)  \nwhere, r Alfv'en is the Alfv'en radius and G is the gravitational constant. The two radii, r cor and r mag, together with the light-cylinder radius, defined as r lc = c / Ω , are associated with the NS and used to study the interacting binary phase. Using the three characteristic radii ( r mag, r cor, and r lc), we define below the three main regimes of accretion:  \n2.3.1. Accretion regime, r mag < ( r cor , r lc)  \nIf 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 disc extends down to the NS surface. Hence, in Equation 6, R NS 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 ˙ M NS = ϵ ˙ M c = (1 -β ) ˙ M c, ϵ being the accretion e ffi ciency. We also assume a κ of unity, similar to Chattopadhyay et al. (2020).", '2.3.2. Propeller regime, r cor ≤ r mag < r lc': '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 ˙ M NS = β ˙ M 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 for 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  \ncondition for the onset of the pulsar wind phase is that the masstransfer rate has decreased below the criteria for thermal-viscous instabilities and can no longer sustain a persistent accretion disc (see Equation 3), under the e ff ect of X-ray irradiation during the RLO phase.  \nOnce 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,  \nL pulsar = 4 π 2 I ˙ P spin P 3 spin . (9)  \nIn the above equation, the initial values of the NS spin ( P spin) and change in NS spin with time ( ˙ P 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 ( ˙ M c), defined by Stevens et al. (1992),  \n˙ M c = -f pulsar × L pulsar 2 v 2 esc GLYPH<18> R c a GLYPH<19> 2 , (10)  \nwhere f pulsar is the e ffi ciency that describes the transfer of energy from the incoming spin-down luminosity to the generated wind loss. It also stands for any geometrical beaming e ff ects that might be present during the process. The other parameters in the above equation include L pulsar as the spin-down luminosity, v esc as the escape velocity from the companion star surface, R c as the companion radius, and a as the binary separation.  \nThe 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 catalogue of BWs, RBs, and huntsmen (including candidates). A total of 57 spiders are used, 52 of which are catalogued in Linares & Kachelrieß (2021, and references therein). The more recent ones are from various other works in the literature, and an updated catalogue 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 M c , 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 catalogue have radio pulses been confirmed and M c , 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.1. Mass-accretion efficiency': "We investigate the formation of di ff erent 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 di ff erent values of f pulsar. The tracks correspond to initial NS masses of 1.3 M ⊙ and have final NS spins  \n≲ 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 e ffi ciencies. Similar behavior is seen for huntsmen spiders ( P orb > 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 f pulsar, more and more BWs and RBs' observed orbital properties are reproduced. The reason for the di ff erence 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, which 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.  \nWe find that all the MESA tracks that pass through the RB region of the observed parameter space, end up in the BW region (see Figure 1). The only exception is for one track corresponding to f pulsar = 0 . 5, that ends with M c of around 0.15 M ⊙ (discussed below in Section 3.4). Hence, to study the e ff ect of the assumed β on pulsar spin-up, we compare only the observed BWs to the final pulsar spins from the grid in Figure 2, corresponding to f pulsar = 0 . 0 and β = 0.0, 0.3, and 0.7. The figure shows the di ff erent binaries with di ff erent initial NS masses, 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. NSs with final masses above 2.5 M ⊙ are considered to have collapsed into BHsand would no longer be observed as radio pulsars. The value β = 0 . 0 leads to a majority of the NSs gaining spins in the subMSP range ( ≲ 1 ms) and having final NS masses ≳ 2 . 5 M ⊙ .  \nWith a β of 0 . 7 (or an accretion e ffi ciency of 0 . 3), the minimum spin period that is reached in all the grids is 1.14 ms, consistent with the fastest spinning Galactic BW observed so far (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.2 ms, which is more than twice the minimum observed BW spin. Looking at the case of β = 0 . 3 in Figure 2, all the observed spins can be reproduced. Comparing the minimum simulated 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 ⊙ ; Romani et al. 2022). There are again some binary tracks that ended with final spins below 1 ms, in the sub-MSP regime (all with M i NS > 1 . 3 M ⊙ ). The NSs with M i NS = 1 . 3 M ⊙ gained a maximum of 1 . 06 M ⊙ and did not spinup to 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 a 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, 1.8 and 2.0 M ⊙ (see Figure 2), 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 -3 spin , and not a ff ecting the orbital evolution as much.  \nFig. 1: Orbital period evolution versus companion mass of the simulated binaries, for β = 0 . 3, having an initial NS of mass 1 . 3 M ⊙ and f pulsar values of 0.0 (top-left), 0.05 (top-right), 0.1 (bottom left), and 0.5 (bottom-right). M c is the mass of the companion to the NS and P orb is the orbital period. The colorbar shows the changing surface abundance of the companion. We also compare our tracks to estimated M c , min and observed P orb 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.  \n<!-- image -->  \nIncreased wind mass loss from the companion would lead to increased orbital widening during the radio pulsar phase.", '3.2. Orbital evolution and surface hydrogen abundance': 'The top left panel in Figure 1 shows the evolutionary tracks with no pulsar wind irradiation ( f pulsar = 0 . 0). The colorbar 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, PSRJ1311-3430 (Pletsch et al. 2012; Romani et al. 2012) and PSRJ1653-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 P orb ≲ 2 hr and M c ≲ 0 . 015 M ⊙  \n(Romani et al. 2016). They are reproduced regardless of the value of f pulsar. Their evolutionary tracks also pass through the ultra-compact X-ray binary region (UCXB), with orbital periods less than 1 hr, 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.  \nWe can understand their evolution better in Figure 3, where we show tracks that pass through the tidarren region (solid gray line, with an M i c = 1 . 5 M ⊙ and P i 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 i NS = 1 . 3 M ⊙ . We also show a binary that goes well into the UCXB region during its  \nFig. 2: Distribution of the final NS spins ( P f spin ) in converging binaries, for β = 0 . 0 (top), β = 0 . 3 (middle), and β = 0 . 7 (bottom). We group the binaries by the initial NS mass ( M i NS ), and if the final mass was above 2 . 5 M ⊙ (hatched histograms, with faded colours) or below it (step histograms, with brighter colours). We also compare to the observed spins of BWs (shown by the gray shaded area).  \n<!-- image -->  \nevolution, to reach the overall minimum orbital period (19 min) in the grid slice (with an M i c = 1 . 0 M ⊙ and P i orb = 2.6 d). With increasing orbital periods, the binaries go from starting RLO with the companion on the MS (at P orb = 0 . 5 d, also called case A RLO), to starting RLO when the companion has exhausted core-H (also called case B RLO, at P orb = 1 . 1 d). This transition occurs at the bifurcation period ( P bif), because companions in wider orbits have more time to evolve and move out of the MS. In the presented example, P bif is between P orb = 0 . 93 d (with P i orb = 2 . 4 d and M i c = 1 . 5 M ⊙ ) and P orb = 1 . 1 d (with P i orb = 2 . 6 d and M i c = 1 . 5 M ⊙ ). We term systems just below P bif as nearP bif binaries, since similar studies have shown that the formation of UCXBs from LMXBs happens just below P bif (Podsiadlowski et al. 2002, 2003; Gossage et al. 2023).  \nFor 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 (corresponding to the black curve in Figure 3) proceeds on a nuclear timescale during the evolution of the stellar MS, any  \nFig. 3: Evolution of orbital period versus companion mass for selected binaries from our simulated grids ( β = 0 . 3, f pular = 0 . 0, and M i NS = 1 . 3 M ⊙ ). M c is the mass of the companion to the NS and P orb is the orbital period. The three di ff erent kinds of binaries that we discuss in the text are studied here, case A ( P i orb = 1 . 0 d), nearP bif ( P i orb = 2 . 4 d), and diverging ( P i 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 f pulsar = 0 . 0. We show the start of the first RLO phase as circles and the formation of a fully convective star as squares (which happens only for two of the tracks shown). We also compare our tracks to estimated M c , min and observed P orb 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 2FGLJ0846.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.  \n<!-- image -->  \nmass loss via winds is not fast enough to dominate over magnetic braking. Additionally, in compact binaries ( ≲ 1 d), angular momentum losses due to gravitational wave radiation are also present 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 with M i c = 1 . 5 M ⊙ and P i 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 or-  \nhich is the drop in orbital period seen in the figure). The gravitational wave radiation sustains this second RLO phase and the orbit eventually expands.  \nFor case B RLO (corresponding to M i c = 1 . 5 M ⊙ and P i 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 e ff ect of magnetic braking is weaker in comparison, we get diverging orbits. 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 systems with diverging orbits are NS / He-WDs in wide orbits. For nearP bif systems, the companion becomes fully convective at a significantly lower mass ( M c ≲ 0 . 1 M ⊙ , see gray and orange lines in Figure 3) than case A systems. They have a delayed turno ff of the first RLO phase, during which the magnetic braking e ff ect only increases as the orbit shrinks ( ∝ P -3 orb ; Rappaport et al. 1983). The reason for the delayed formation of a fully convective structure is the first RLO phase; as the companion loses mass, its central temperature decreases due to loss in internal pressure. The case A binary has an RLO phase that started earlier compared to the nearP bif binary, providing enough time to cool internal structure of the star and making it fully convective while the companion is still on the MS. Additionally, the nearP bif binary also has strong stellar winds that, along with the magnetic braking driven RLO, strip the outer layers of the star. The more there is a delay in magnetic braking turno ff , the further the orbit contracts into the UCXB region; the binary that reaches the minimum P orb (19 min; Figure 3) in the grid slice (with an M i c = 1 . 0 M ⊙ and P i orb = 2 . 6 d), never forms a fully convective companion.  \nThere is another source in the same observed parameter space as the two tidarrens (see Figure 3), for which the surface Habundance 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 f pulsar, the e ff ect 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 P bif. As there is no irradiation included in their evolution, these systems could depict some of the nonirradiated RBs observed. However, 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 f pulsar to form RBs in the subsequent sections.', '3.3. Effect of pulsar wind irradiation': 'Figure 1 shows the evolutionary grids with an M i NS = 1 . 3 M ⊙ and varying pulsar wind e ffi ciency ( f pulsar). 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 f pulsar = 0 . 0 and 0.5 and initial NS masses 1 . 3 M ⊙ . Many of the binaries shown in Figure 4 (around 50% in simulated spiders) end up su ff ering a dynamical instability (shown by diamond symbols in the figure); the companion exceeds its Roche lobe to fill out the gravitational equipotential surface passing through the sec-  \nLagrange point (Misra et al. 2020). This is a dynamically unstable process and leads to the binary merging.  \nThe companion at the end of the first RLO phase is a fully convective star. Convective envelopes are known to expand rapidly on mass loss (Hjellming & Webbink 1987), driving sudden high mass-transfer rates during the second RLO phase and a dynamical instability. This occurs with and without pulsar wind irradiation. Irradiation removes some matter from the surface of the companion star and causes rapid stellar expansion. Increased f pulsar 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 f pulsar). Despite this, we include these dynamically unstable systems in our study since by the time the instability onset happens, the average MSP binary is already ∼ 10 1 0 yrs old and could have been observed during the preceding MSP phase. Hence, there is a reasonable chance to observe these binaries as BWbinaries within a Hubble time (1 . 4 × 10 10 yrs). For a 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.  \nGoing from f pulsar = 0 . 0 to f pulsar = 0 . 5 in the converging systems (Figure 4), the final orbital periods widen significantly and this e ff ect depends on the initial properties of binary. The widest orbits are seen with M i c = 1 . 0 M ⊙ , 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 5 shows the orbital evolution corresponding to three M i c (1.0, 1.5, and 2.5 M ⊙ ) and f pulsar = 0 . 0 to f pulsar = 0 . 5, for various values of P i orb . For case A RLO binaries ( P i orb = 1 . 0 d; shown by black lines in the left panel), we see increased orbital widening with f pulsar = 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 f pulsar. Going from M i c = 1 . 0 M ⊙ to 1 . 5 M ⊙ (for f pulsar = 0 . 5) in right and left panels, respectively, 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 i c = 4 . 0 M ⊙ , the accreted matter decreases down to about 0 . 004 M ⊙ as the binary su ff ers a dynamical instability soon after the first RLO-onset. This is significantly less than 1 . 0 M ⊙ accreted with M i c = 3 . 0 M ⊙ . 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.  \nNext we look at the tracks for P i orb = 2 . 4 d in Figure 5, corresponding to the gray curves (left panel). For M i c = 1 . 5 M ⊙ , the e ff ect of f pulsar = 0 . 5 is significant on the orbital evolution. The nearP bif binary maintains a wide orbit (around 0.2 d) crossing the RB regime. Going to M i c = 2 . 5 M ⊙ and P i orb = 2 . 0 d (gray curves in right panel), the nearP bif tracks no longer pass through the RB region with the heavier companion (irrespective of f pulsar). For intermediate-mass companions ( ≳ 2 . 0 M ⊙ ), the initial angular momentum loss by RLO dominates over magnetic braking. 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 5) and magnetic braking starts to dom-  \nFig. 4: Allowed initial parameter space to form spider MSPs, for an initial NS of mass 1 . 3 M ⊙ and f pulsar values of 0.0 (left) and 0.5 (right). M i c is the initial mass of the companion to the NS and P i orb is the initial orbital period. The colorbar shows the final orbital periods ( P f 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 (triangle symbols) or ended with one of the criteria for successful termination (star symbols). The diverging symbols are shown by gray crosses. We also show the nearP bif 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 J09520607 (with observed spin period of 1.41 ms).  \n<!-- image -->  \nfollowing 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 nearP bif 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).  \nWe can determine also the evolutionary history of PSRJ0952-0607, the fastest spinning Galactic BW so far (with spin 1.41 ms), from our simulations, considering β = 0 . 3 and f pulsar = 0 . 5. To account for uncertainties during mass transfer and in our calculation of P spin, we consider a 10% margin and take all binary tracks with P spin ≲ 1 . 55 ms. The tracks are shown in Figure 6. We compare the simulated tracks to the estimated orbital parameters, M NS = 2 . 35 + 0 . 17 -0 . 17 M ⊙ , M c = 0 . 032 + 0 . 002 -0 . 002 M ⊙ , and P orb = 6 . 42 hr (Romani et al. 2022), and select the best matching track (within the errors). The best matching track started with M i NS = 1 . 3 M ⊙ , M i c = 2 . 0 M ⊙ and P i 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.', '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).  \nThey evolve from a case B RLO phase in LMXBs. In terms of companion masses, they are similar to RBs, with the di ff erence 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 f pulsar. Contrary to RBs and BWs, with a high f pulsar the orbits do not widen with increasing f pulsar. Figure 5 (left panel, with M i c = 1 . 5 M ⊙ and P i orb = 2 . 6 d, shown by blue curves) shows that wider orbits correspond to f pulsar = 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. The resulting e ff ect is a decreased expansion of the orbits.  \nWe 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. Fig 7 shows the final spin period distributions 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 a NS with an initial mass of 1 . 3 M ⊙ to spin-up to 2.66 ms. With a β = 0 . 3, only NSs with initial masses ≳ 1 . 5 M ⊙ are able to reproduce 1FGL J1417.7440.  \nFig. 5: Orbital period evolution versus companion mass of selected binaries from our simulated grid shown in Figure 1. M c is the mass of the companion to the NS and P orb is the orbital period. The figure shows the e ff ect of pulsar wind irradiation, comparing f pulsar = 0 . 0 (dashed lines) and 0.5 (solid lines). Left panel shows the tracks with M i c = 1 . 5 M ⊙ , P i orb = 1 . 0, 2.4, and 2.6 d. Right panel shows two initial companion masses, M i c = 1 . 0 (with P i orb = 1 d) and 2 . 5 M ⊙ ( P i orb = 2 . 0 and 2.2 d). We compare our tracks to M c , min and P orb 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 J102347.67 + 003841.2 ( M c = 0 . 22 M ⊙ and P orb = 4 . 75 hr; Shahbaz et al. 2019; Thorstensen & Armstrong 2005), IGR J18245-2452 ( M c , min = 0 . 17 M ⊙ and P orb = 11 . 03 hr; Papitto et al. 2013), and XSS J12270-4859 ( M c , min = 0 . 15 M ⊙ and P orb = 6 . 91 hr; Roy et al. 2015; de Martino et al. 2015; Bassa et al. 2014).  \n<!-- image -->  \nTable 1: Observed and estimated properties of the two observed huntsmen spiders. M c , min is the minimum mass, calculated assuming a pulsar mass of 1 . 4 M ⊙ and inclination of 90°, P orb is the orbital period, P spin is the observed pulsar spin, M NS is the estimated pulsar mass, M c is the companion mass observed from Fermi-LAT, and e is the observed eccentricity.  \n.  \n.  \nStrader et al. (2015),  \n2  \nSwihart et al. (2018),  \n3  \nAssuming β = 0 . 0, and looking at the final mass distributions of the NSs in Figure 8, we see the maximum final NS masses that we can get in huntsmen spiders. We also compare our distributions to the estimated masses of the pulsars in 1FGL J1417.7440 (1 . 62 + 0 . 43 -0 . 17 ; Swihart et al. 2018) and 2FGL J0846.0 + 2820 (1 . 96 + 0 . 41 -0 . 41 M ⊙ ; Swihart et al. 2017). Even if the binaries had β = 0 . 0, for a huntsman binary to have a NS with mass greater than 2.0 M ⊙ , it should have M i NS > 1 . 3 M ⊙ (with 1.3 M ⊙ , the maximum mass reached is 1.85 M ⊙ ). Therefore, with β > 0 . 0, the initial mass should be even higher to get a high mass NS at the end of RLO.  \nThe observed and estimated properties of the two currently known huntsmen spiders are described in Table 1. The spin of 1FGLJ1417.7-440 (2.66 ms) can be used to constrain its evo-  \nlutionary 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 Figures 7 and 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 di ff erent values of f pulsar. 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 . 07 -0 . 03 M ⊙ 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 bi-  \nFig. 6: Orbital period versus companion mass for diverging binary tracks from our simulated grids that have similar NS spin as the observed PSR J0952-0607 (1.41 ms) within a 10 % margin (dashed gray curves). M c is the mass of the companion to the NS and P orb is the orbital period. All the tracks have β = 0 . 3 and f pulsar = 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.  \n<!-- image -->  \nFig. 7: Distribution of the final NS spins ( P f spin ) in diverging binaries, for β = 0 . 0 (top), β = 0 . 3 (middle), and β = 0 . 7 (bottom). We group the binaries by the initial NS mass ( M i NS ) and also compare the spins to the observed spin of the huntsmen source 1FGLJ1417.7-440 (2.66 ms, shown by the vertical gray dashed lines; Camilo et al. 2016).  \n<!-- image -->  \ntrack that best matches the source's observed orbital properties (within the errors) corresponds to β = 0 . 0 and f pulsar = 0 . 1, and with the initial parameters, M i c = 2 . 0 M ⊙ , M i NS = 1 . 3 M ⊙ ,  \nFig. 8: Distribution of the final NS masses ( M f NS ) in diverging binaries (for β = 0 . 0), grouping the binaries by the initial NS mass ( M i NS ). We also compare to the estimated pulsar masses of the two huntsmen sources 1FGL J1417.7-440 (1 . 62 + 0 . 43 -0 . 17 M ⊙ , shown by the vertical gray dashed lines and gray shaded area for the uncertainties; Swihart et al. 2018) and 2FGL J0846.0 + 2820 (1 . 96 + 0 . 41 -0 . 41 M ⊙ , shown by the vertical blue dashed lines and blue shaded area for the uncertainties; Swihart et al. 2017).  \n<!-- image -->  \nand P i 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.  \nAs for the source, 2FGL J0846.0 + 2820, since we do not know its spin, there is a wider range of binary tracks which could be used to explain it and almost all f pulsar values can reproduce its orbital parameters (see Figure 1). Compared to 1FGL J1417.7440 (and other spider systems), 2FGL J0846.0 + 2820 has a nonnegligible 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 e ffi ciently circularized the orbits (Verbunt & Phinney 1995). This observed eccentricity could be the result of the interaction with a circumbinary disc (Antoniadis 2014) or ine ffi cient tides (such as in radiative envelopes, which would be the case for MS stars; Hurley et al. 2002).  \nAlthough 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 f pulsar = 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 colorbar 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  \nFig. 9: Allowed initial parameter space to form huntsmen spiders ( f pulsar = 0 . 0), for an initial NS of mass 1 . 3 M ⊙ (left) and 2 . 0 M ⊙ (right). M i c is the initial mass of the companion to the NS and P i orb is the initial orbital period. The colorbar shows the final NS spins ( P f spin ), the gray stars denoting 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.  \n<!-- image -->  \nFig. 10: Orbital period versus companion mass for diverging binary tracks from our simulated grids that have similar NS spin as the observed huntsman source 1FGL J1417.7-440 (2.66 ms) within a 20 % margin (dashed gray curves). M c is the mass of the companion to the NS and P orb is the orbital period. We further show the binary that best matches the sources in its orbital evolution in blue curves, having β = 0 . 0 , f pulsar = 0 . 1, M i c = 2 . 0 M ⊙ , M i NS = 1 . 3 M ⊙ , and P i orb = 2 . 6 d.  \n<!-- image -->  \n2 . 0 M ⊙ and 2.5-3.5 d to being centered around 2 . 5 M ⊙ and 1.63.0 d.  \nThe boundary of the huntsmen parameter space depends on the initial q and P i 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 e ffi cient NS spin-up. If it is too massive, the RLO phase is very short as the star evolves quick. Additionally, 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 i NS = 2 . 0 M ⊙ ), wealso see an isolated binary with M i c = 4 . 5 M ⊙ and P i orb = 1 . 8 d that has a NS with MSP-like spin. The final binary is a wide HeWD / 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 do not find a similar binary with a lower initial NS mass.", '3.5. Transitional MSPs': '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 J102347.67 + 003841.2 (Archibald et al. 2009), IGR J18245-2452 (Papitto et al. 2013), and XSSJ12270-4859 (Bassa et al. 2014). The companion  \nmasses for all of them are in the range 0.15 to 0.2 M ⊙ and orbital periods 4 to 11 hr (Papitto & de Martino 2022, and references therein), which put them in the RB regime of the observed parameter space. Figure 5 shows the three tMSPs along with some of the tracks from our simulated binaries. Two of these lie the nearP bif tracks for M i c = 1 . 5 M ⊙ (gray lines in the figure), where the tracks di ff erentiate for f pulsar = 0 . 0 and 0.5. The point of di ff erentiation 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 disc instabilities (Lasota 2001) which are not yet fully understood.', '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 e ffi ciency. For BWs, they found f pulsar ≲ 0 . 08 to be the required irradiation e ffi ciency, while for RBs, it was f pulsar ≳ 0 . 08. They explored mainly case A RLO systems, M i c around 1 . 0 M ⊙ , M i NS around 1 . 4 M ⊙ , and P i orb from 0.08 to 1.4 d, and found increased orbital expansion with increasing f pulsar. Qualitatively, we also see similar behavior in our results (see Section 3.3). Particularly for BWs, looking at similar systems in Figure 1, for f pular = 0 . 0 to 0.05, we cover the lower-left region of the observed BW space. The degree to which the orbit expands for f pular ≳ 0 . 1 is lesser than Chen et al. (2013).  \nFor RBs, Chen et al. (2013) found them to be produced from case A binaries ( ∼ 1 d) with strong irradiation ( f pulsar ≳ 0 . 08). In our work, the majority of RBs are produced at much wider initial orbits (nearP bif at RLO onset, starting at around 2.0 to 2.4 d at ZAMS), along with high irradiation. With f pulsar = 0 . 5, we also see rapid expansion of the case A binaries (solid black curve and M i c = 1 . 5 in the right panel of Figure 5). However, this expansion in orbits is not enough to bring them completely into the RB regime. Benvenuto et al. (2014) and Benvenuto et al. (2015) studied the e ff ect 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 both BWs and RBs in our simulated tracks, by the time the spider phase begins the binaries are ≳ 10 9 yr old, as most of the prior evolution happened during the MS phase, which is long for low-mass stars.  \nThe need for a persistent magnetic braking was also suggested by Ginzburg & Quataert (2020) and Ginzburg & Quataert (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 f pulsar = 0 . 5 in right  \npanel of Figure 5). With increasing companion masses, we cover more of the observed BW parameter space. The MSP source PSRJ1953 + 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 catalogue. Still, Pan et al. (2023) suggested that this source was formed from a similar process as UCXBs (with P i orb values near P bif), which is similar to our results, albeit at di ff erent metallicities.', '4.2. Accretion efficiency and non-conservative evolution': 'As we saw in Section 2.3, in the literature non-conservative mass-accretion e ffi ciencies ( β ≳ 0 . 5) are favoured 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 e ffi ciencies are better suited ( β ≲ 0 . 3). We also find that with a β of 0.7 (see Figure 2), we can reproduce the observed spin and estimated pulsar mass of PSR J1738 + 0333, 5.85 ms and 1 . 47 + 0 . 07 -0 . 06 M ⊙ , respectively. 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 e ff ects 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 e ff ectiveness, we assume a low e ffi ciency of loss of angular momentum during this phase. We also assumed an full e ffi ciency for transfer of angular momentum from the accreted matter during RLO (see Section 2.3.1).  \nAs for accretion disc instabilities, we include a model that describes the presence of a accretion disc above a critical masstransfer rate, below which an X-ray irradiated disc would be unstable (see Section 2.3). Since we have already included ine ffi cient mass accretion due to these two processes, the massaccretion e ffi ciency ( ϵ = 1 -β ) which 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 ine ffi cient accretion, since we use simplistic approximations. However, the binary may evolve differently due to a high value of fixed β throughout compared to a low β with restricted accretion at some instances. To test this, we ran a MESA model without the X-ray irradiated accretion disc assumptions and with β = 0 . 7 and f pulsar = 0 . 0, with M i c = 1 . 5 M ⊙ , M i NS = 1 . 3 M ⊙ , and P i orb = 1 . 0 d. We compared it to runs with same initial parameters, with β = 0 . 3 and the X-ray irradiated disc assumption. We found that without the restriction on accretion, the NS spun up to 0.67 ms, while with the restriction it spun up to 2.9 ms. Assuming a fixed β might result in overestimation of highly spinning pulsars, as most of these systems undergo case A RLO that is a long-lasting phase with moderately high mass-transfer rate (10 -10 M ⊙ yr -1 to 10 -8 M ⊙ yr -1 ). Despite this di ff erence, our results do not contradict the LMXB observations supporting highly non-conservative behavior.', '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 pul-  \nsar 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 ( ˙ Ω ) with the angular frequency ( Ω ), which is ˙ Ω ∝ -Ω 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 10 6 (Hobbs et al. 2004; Biryukov et al. 2012; Zhang & Xie 2012).  \nStudies 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 f pulsar = 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.  \nAccording 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 f spin ≤ 30 . 0 ms). 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. Haskell et al. (2018) found that masses greater than 2 . 0 M ⊙ have an upper limiting frequency of around 810 Hz (lower P spin 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 i NS = 1 . 3 M ⊙ , it is at 1.38 ms and for M i NS = 2 . 0 M ⊙ , it is at 2.66 ms. For a NS lighter than M i NS = 1 . 3 M ⊙ , 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-ofstate 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, subMSPs may not exist if those NSs have collapsed into BH.', '4.4. Uncertainties in RB 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  \n<!-- image -->  \nFig. 11: Accreted matter versus final NS spins ( P f spin ) in our simulated grid with β = 0 . 3. We show the di ff erent initial masses as di ff erent colors (shown in the legend) and distinguish between the NSs that had M f NS ≥ 2 . 5 M ⊙ (crosses) and those where M f NS < 2 . 5 M ⊙ (circles). We also compare our results to the equilibrium spin period calculations from Tauris et al. (2012), done for M i NS = 1 . 3 M ⊙ (solid blue curve) and 2 . 0 M ⊙ (solid red curve). The gray dashed line shows the observed spin for PSRJ0952-0607 (spin period of 1.41 ms).  \n<!-- image -->  \nfrom LMXBs starting from close to P bif (Podsiadlowski et al. 2002; van der Sluys et al. 2005b,a; Istrate et al. 2014). 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 nearP bif) and companion masses from which UCXBs can be produced (van der Sluys et al. 2005b,a; Istrate et al. 2014). Even though strong irradiation widens nearP bif orbits and some of them might no longer be UCXBs, some UCXBs can still be reproduced (see f pulsar = 0 . 5 in Figure 1).  \nThe dependence of P bif on magnetic braking assumptions and the e ff ects 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. 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 M c ∼ 0 . 2 M ⊙ and the companion has lost its radiative core. Gossage et al. (2023) compared single star and LMXB observations and found saturated magnetic braking to best reproduce the observations. Still the standard magnetic braking prescription can reproduce LMXBs with periods less than 1 hr to around 10 d, better than some other prescriptions (Echeveste et al. 2024; Gossage et al. 2023). Carrying out a population synthesis study  \nof UCXBs would also tell us what fraction of LMXBs in a stellar population start RLO close to their P bif. 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 us in understanding the formation of spiders.', '5. Conclusions': "Compact MSP binaries (or spider binaries) are good laboratories for studying binary evolution. We calculated the NS spin evolution and identified binaries that formed MSPs (using the criteria for NS spins: P f pular ≤ 30 ms). Assuming a model for an X-ray irradiated accretion disc and three di ff erent values of the accretion e ffi ciency (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 di ff erent values for the e ffi ciency 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. Comparing our simulated grids to the observed spider populations, we arrive at the following conclusions:  \n- -With the assumption of unstable accretion disc due to X-ray irradiation during the mass-transfer phase, we need an accretion e ffi ciency of 70% to reproduce the observed spins and pulsar masses of BW systems. High accretion e ffi ciency ( > 70%) also allows for the formation of the huntsmen spiders from a binary with a canonical NS.\n- -For BW formation, a high pulsar wind e ffi ciency ( f pulsar = 0 . 5) are required to reproduce the full observed parameter space, although lower e ffi ciencies (0 . 0 < f pulsar < 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 . 01M ⊙ companions) require no irradiation. We also predict that J0636 + 5128 (similar to its neighbours) 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 nearP bif and pass through an UCXB phase. Their evolution is dependent on the magnetic braking prescription used.\n- -With high pulsar wind irradiation ( f pulsar = 0 . 5), the orbits of the nearP bif binaries expand drastically, covering most of the RB region. This is true for the companion masses ≲ 1 . 5 M ⊙ , the nearP bif 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 initial parameter space can reproduce the observations.\n- -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 Xray irradiated accretion disc) a NS in diverging orbits with initial mass 1 . 3 M ⊙ gains a maximum of 0 . 55 M ⊙ . If a huntsman pulsar is discovered to have a pulsar mass ≳ 2 . 0 M ⊙ , most likely it was formed more massive than 1 . 3 M ⊙ .  \n- -We identified the prior evolution of the pulsars: BW PSRJ0952-0607 and huntsman 1FGLJ1417.7-440, by matching their observed spins and orbital properties to binary tracks. For PSR J0952-0607, high irradiation ( f pulsar = 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 ( f pulsar = 0 . 1) matches the observations well, though irradiation is not required to reproduce the orbital properties of these diverging NS binaries.\n- -NSs in LMXBs / 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 i NS = 1 . 3 M ⊙ ). 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 i 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 PSRJ0952-0607 (with observed spin 1.41 ms and pulsar mass of 2 . 35 M ⊙ ). Hence, fast spinning NSs are key in identifying the boundary between NSs and BHs.\n- -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, the surviving binaries for which would be isolated spinning NSs. These could explain the observations of isolated MSPs.  \nSince this study involved running detailed MESA simulations with a fixed set of initial parameters, we do 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 a ff ecting the evolution prior to the XRB phase (like the supernova mechanism) and more directly involved physical processes (like magnetic braking and pulsar wind irradiation). Since LMXBs have very unstable accretion discs causing outbursts that would a ff ect 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 e ff ects that would cause non-conservative mass transfer, these instabilities and their e ff ects 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).  \nAcknowledgements. The authors thank Thomas Tauris and Hai-Liang Chen 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': "Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 848 Abolmasov, P., Biryukov, A., & Popov, S. B. 2024, Galaxies, 12, 7  \n- Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Nature, 300, 728\n- Antoniadis, J. 2014, ApJ, 797, L24\n- Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448\n- Antoniadis, J., Tauris, T. M., Ozel, F., et al. 2016, arXiv e-prints, arXiv:1605.01665  \nAntoniadis, J., van Kerkwijk, M. H., Koester, D., et al. 2012, MNRAS, 423, 3316 Archibald, A. M., Bogdanov, S., Patruno, A., et al. 2015, ApJ, 807, 62 Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Science, 324, 1411 Archibald, R. F., Gotthelf, E. V., Ferdman, R. D., et al. 2016, ApJ, 819, L16 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 Bassa, C. G., Patruno, A., Hessels, J. W. T., et al. 2014, MNRAS, 441, 1825 Bassa, C. G., Pleunis, Z., Hessels, J. W. T., et al. 2017, ApJ, 846, L20 Bejger, M., Blaschke, D., Haensel, P., Zdunik, J. L., & Fortin, M. 2017, A&A,  \n- 600, A39  \nBenvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2014, ApJ, 786, L7  \nBenvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2015, ApJ, 798, 44  \n- Bhattacharya, D. & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1\n- Bhattacharyya, S. & Chakrabarty, D. 2017, ApJ, 835, 4\n- Bildsten, L. 1998, ApJ, 501, L89\n- Biryukov, A., Beskin, G., & Karpov, S. 2012, MNRAS, 420, 103  \nBreton, R. P., van Kerkwijk, M. H., Roberts, M. S. E., et al. 2013, ApJ, 769, 108  \nCamilo, F., Reynolds, J. E., Ransom, S. M., et al. 2016, ApJ, 820, 6  \n- Chattopadhyay, D., Stevenson, S., Hurley, J. R., Rossi, L. J., & Flynn, C. 2020, MNRAS, 494, 1587  \nChen, H.-L., Chen, X., Tauris, T. M., & Han, Z. 2013, ApJ, 775, 27  \n- Chen, W.-C. & Podsiadlowski, P. 2016, ApJ, 830, 131\n- Cumming, A., Zweibel, E., & Bildsten, L. 2001, ApJ, 557, 958  \nde Martino, D., Papitto, A., Belloni, T., et al. 2015, MNRAS, 454, 2190  \nDemorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels,  \nJ. W. T. 2010, Nature, 467, 1081  \nDeng, Z.-L., Li, X.-D., Gao, Z.-F., & Shao, Y. 2021, ApJ, 909, 174  \nDraghis, P. & Romani, R. W. 2018, ApJ, 862, L6  \n- Dubus, G., Lasota, J.-P., Hameury, J.-M., & Charles, P. 1999, MNRAS, 303, 139 Echeveste, M., Novarino, M. L., Benvenuto, O. G., & De Vito, M. A. 2024, MNRAS\n- El-Badry, K., Conroy, C., Fuller, J., et al. 2022, MNRAS, 517, 4916\n- Ergma, E. & Fedorova, A. V. 1992, A&A, 265, 65\n- Ertan, U. 2015, arXiv e-prints, arXiv:1504.03996\n- Ertan, Ü. 2018, MNRAS, 479, L12\n- Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2023, ApJS, 264, 45\n- Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2022, arXiv e-prints, arXiv:2202.05892\n- Fruchter, A. S., Stinebring, D. R., & Taylor, J. H. 1988, Nature, 333, 237\n- Gerlach, U. H. 1968, Phys. Rev., 172, 1325\n- Ginzburg, S. & Quataert, E. 2020, MNRAS, 495, 3656\n- Ginzburg, S. & Quataert, E. 2021, MNRAS, 500, 1592\n- Glendenning, N. K. & Kettner, C. 2000, A&A, 353, L9\n- Gossage, S., Kalogera, V., & Sun, M. 2023, ApJ, 950, 27\n- Hamil, O., Stone, J. R., Urbanec, M., & Urbancová, G. 2015, Phys. Rev. D, 91, 063007\n- Haskell, B., Zdunik, J. L., Fortin, M., et al. 2018, A&A, 620, A69  \nHessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901  \n- Hjellming, M. S. & Webbink, R. F. 1987, ApJ, 318, 794\n- Hobbs, G., Lyne, A. G., Kramer, M., Martin, C. E., & Jordan, C. 2004, MNRAS, 353, 1311  \nHurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897  \nIllarionov, A. F. & Sunyaev, R. A. 1975, A&A, 39, 185  \n- Istrate, A. G., Tauris, T. M., & Langer, N. 2014, A&A, 571, A45\n- Jacoby, B. A., Hotan, A., Bailes, M., Ord, S., & Kulkarni, S. R. 2005, ApJ, 629, L113\n- Jia, K. & Li, X. D. 2014, ApJ, 791, 127\n- Jia, K. & Li, X.-D. 2015, ApJ, 814, 74\n- Johnston, S. & Karastergiou, A. 2017, MNRAS, 467, 3493  \nKaplan, D. L., Bhalerao, V. B., van Kerkwijk, M. H., et al. 2013, ApJ, 765, 158  \nKiel, P. D., Hurley, J. R., Bailes, M., & Murray, J. R. 2008, MNRAS, 388, 393  \nKluzniak, W., Ruderman, M., Shaham, J., & Tavani, M. 1988, Nature, 334, 225  \nKonar, S. & Bhattacharya, D. 1997, MNRAS, 284, 311  \n- Konar, S. & Bhattacharya, D. 1999a, MNRAS, 303, 588\n- Konar, S. & Bhattacharya, D. 1999b, MNRAS, 308, 795\n- Kong, A. K. H., Jin, R., Yen, T. C., et al. 2014, ApJ, 794, L22  \nKumari, S., Bhattacharyya, B., Sharan, R., et al. 2024, ApJ, 961, 155  \nLasota, J.-P. 2001, New A Rev., 45, 449  \n- Li, L., Han, Z., & Zhang, F. 2004, MNRAS, 355, 1383\n- Linares, M. 2020, in Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. Palermo, 23\n- Linares, M. & Kachelrieß, M. 2021, J. Cosmology Astropart. Phys., 2021, 030 Linares, M., Shahbaz, T., & Casares, J. 2018, ApJ, 859, 54\n- Liu, X.-J., You, Z.-Q., Chen, Z.-C., et al. 2024, ApJ, 962, 80\n- Ma, B. & Li, X.-D. 2009, ApJ, 691, 1611  \nManchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005a, AJ, 129, 1993  \nManchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005b, VizieR Online  \nData Catalog, VII  \n/  \n245  \n- Mineshige, S. 1993, Ap&SS, 210, 83\n- Misra, D., Fragos, T., Tauris, T. M., Zapartas, E., & Aguilera-Dena, D. R. 2020, A&A, 642, A174\n- Nedreaas, I. F. 2024, Master's thesis, Norwegian University of Science and Technology\n- Nelemans, G. & Jonker, P. G. 2010, New A Rev., 54, 87\n- Nelemans, G., Jonker, P. G., & Steeghs, D. 2006, MNRAS, 370, 255\n- Nelson, L. A., Rappaport, S. A., & Joss, P. C. 1986, ApJ, 304, 231\n- Osłowski, S., Bulik, T., Gondek-Rosi'nska, D., & Belczy'nski, K. 2011, MNRAS, 413, 461\n- Pallanca, C., Mignani, R. P., Dalessandro, E., et al. 2012, ApJ, 755, 180\n- Pan, Z., Lu, J. G., Jiang, P., et al. 2023, Nature, 620, 961\n- Papaloizou, J. & Pringle, J. E. 1978, MNRAS, 184, 501\n- Papitto, A. & de Martino, D. 2022, in Astrophysics and Space Science Library,  \nVol. 465, Astrophysics and Space Science Library, ed. S. Bhattacharyya,  \nA. Papitto, & D. Bhattacharya, 157-200  \n- Papitto, A., de Martino, D., Belloni, T. M., et al. 2015, MNRAS, 449, L26  \nPapitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517  \n- Patruno, A., Haskell, B., & D'Angelo, C. 2012, ApJ, 746, 9\n- Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3\n- Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4\n- Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15\n- Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34  \nPaxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10  \nPeters, P. C. 1964, Physical Review, 136, 1224  \nPletsch, H. J., Guillemot, L., Fehrmann, H., et al. 2012, Science, 338, 1314  \nPodsiadlowski, P., Han, Z., & Rappaport, S. 2003, MNRAS, 340, 1214  \nPodsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107  \nPolzin, E. J., Breton, R. P., Bhattacharyya, B., et al. 2020, MNRAS, 494, 2948  \nPringle, J. E. 1981, ARA&A, 19, 137  \nPylyser, E. H. P. & Savonije, G. J. 1989, A&A, 208, 52  \nRaaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, ApJ, 918, L29  \nRadhakrishnan, V. & Srinivasan, G. 1982, Current Science, 51, 1096  \nRappaport, S., Verbunt, F., & Joss, P. C. 1983, ApJ, 275, 713  \nRitter, H. & King, A. R. 2001, in Astronomical Society of the Pacific Conference  \nSeries,  \nVol. 229, Evolution of Binary and Multiple Star Systems, ed.  \nP. Podsiadlowski, S. Rappaport, A. R. King, F. D'Antona, & L. Burderi, 423  \nRoberts,  \nM. S. E. 2013, in Neutron Stars and Pulsars: Challenges and  \nOpportunities after 80 years, ed. J. van Leeuwen, Vol. 291, 127-132  \nRomani, R. W., Filippenko, A. V., & Cenko, S. B. 2014, ApJ, 793, L20  \n- Romani, R. W., Filippenko, A. V., Silverman, J. M., et al. 2012, ApJ, 760, L36\n- Romani, R. W., Graham, M. L., Filippenko, A. V., & Zheng, W. 2016, ApJ, 833, 138\n- Romani, R. W., Kandel, D., Filippenko, A. V., Brink, T. G., & Zheng, W. 2022, ApJ, 934, L17\n- Roy, J., Ray, P. S., Bhattacharyya, B., et al. 2015, ApJ, 800, L12  \nRuderman, M., Shaham, J., Tavani, M., & Eichler, D. 1989, ApJ, 343, 292  \nRuderman, M. A. & Sutherland, P. G. 1975, ApJ, 196, 51  \n- Shahbaz, T., Linares, M., Rodríguez-Gil, P., & Casares, J. 2019, MNRAS, 488, 198  \nSkumanich, A. 1972, ApJ, 171, 565  \n- Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, A&A, 327,", '620': '- Spitkovsky, A. 2006, ApJ, 648, L51\n- Spruit, H. C. & Ritter, H. 1983, A&A, 124, 267\n- Stappers, B. W., Bailes, M., Lyne, A. G., et al. 1996a, ApJ, 465, L119\n- Stappers, B. W., Bessell, M. S., & Bailes, M. 1996b, ApJ, 473, L119\n- Stappers, B. W., van Kerkwijk, M. H., Lane, B., & Kulkarni, S. R. 1999, ApJ, 510, L45\n- Stevens, I. R., Rees, M. J., & Podsiadlowski, P. 1992, MNRAS, 254, 19P\n- Stovall, K., Lynch, R. S., Ransom, S. M., et al. 2014, ApJ, 791, 67\n- Strader, J., Chomiuk, L., Cheung, C. C., et al. 2015, ApJ, 804, L12\n- Strader, J., Swihart, S., Chomiuk, L., et al. 2019, ApJ, 872, 42\n- Sturrock, P. A. 1971, ApJ, 164, 529\n- Swihart, S. J., Strader, J., Johnson, T. J., et al. 2017, ApJ, 851, 31\n- Swihart, S. J., Strader, J., Shishkovsky, L., et al. 2018, ApJ, 866, 83\n- Tauris, T. M., Langer, N., & Kramer, M. 2012, MNRAS, 425, 1601\n- Tauris, T. M. & Savonije, G. J. 1999, A&A, 350, 928\n- Tauris, T. M. & van den Heuvel, E. P. J. 2006, in Compact stellar X-ray sources, Vol. 39, 623-665\n- Thorstensen, J. R. & Armstrong, E. 2005, AJ, 130, 759  \nTurchetta, M., Linares, M., Koljonen, K., & Sen, B. 2023, MNRAS, 525, 2565  \nvan der Sluys, M. V., Verbunt, F., & Pols, O. R. 2005a, A&A, 431, 647  \nvan der Sluys, M. V., Verbunt, F., & Pols, O. R. 2005b, A&A, 440, 973 van Paradijs, J. 1996, ApJ, 464, L139 Verbunt, F. & Phinney, E. S. 1995, A&A, 296, 709 Verbunt, F. & Zwaan, C. 1981, A&A, 100, L7 Wagoner, R. V. 1984, ApJ, 278, 345 Wang, B., Chen, W.-C., Liu, D.-D., et al. 2021, MNRAS, 506, 4654 Wijnands, R. & van der Klis, M. 1998, Nature, 394, 344 Woan, G., Pitkin, M. D., Haskell, B., Jones, D. I., & Lasky, P. D. 2018, ApJ, 863, L40 Ye, C. S., Kremer, K., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2019, ApJ, 877, 122 Zhang, C. M. 1998, A&A, 330, 195 Zhang, C. M. & Kojima, Y. 2006, MNRAS, 366, 137 Zhang, C. M., Wang, J., Zhao, Y. H., et al. 2011, A&A, 527, A83 Zhang, S.-N. & Xie, Y. 2012, ApJ, 761, 102', 'Appendix A: Numerical stability of MESA tracks': '<!-- image -->  \nFig. A.1: Stellar luminosity ( L ) versus e ff ective temperature ( T e ff ) 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 di ff erent 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 (that ended without any numerical instability) in solid orange.  \n<!-- image -->  \nFragos et al. (2022, 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 re-run 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 re-run κ max value used by Fragos et al. (2022) is 0.5 cm 2 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  \nin 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 ⊙ , initial NS mass of 2.0 M ⊙ , and initial orbital period of 1.2 d, that failed due to numerical instability. We also show a neighbouring 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 behaviour to their close neighbours, we re-run the failed binaries in two batches, first we used a κ max of 10 cm 2 g -1 . Then we run the binaries that failed again, with a κ max of 5 cm 2 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.', '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 e ffi ciencies for our study.  \nFig. B.1: Orbital period evolution versus companion mass of the simulated binaries, for β = 0 . 7, having an initial NS of mass 1 . 3 M ⊙ and f pulsar values of 0.0 (top-left), 0.05 (top-right), 0.1 (bottom left), and 0.5 (bottom-right). M c is the mass of the companion to the NS and P orb is the orbital period. We also compare our tracks to estimated M c , min and observed P orb for BWs (black crosses) and RBs (red crosses), and a huntsmen spider (green cross), 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 initial binaries have a MS star (range of 1.0 to 4.5 M ⊙ ) and initial periods of 0.6 to 4.0 d.  \n<!-- image -->'}