© ESO 2020

1. Introduction

Ultra-luminous X-ray sources (or ULXs) are extra-galactic sources first discovered by Fabbiano (1989) that have been observed to have X-ray luminosities of L_X ≳ 10³⁹ erg s⁻¹. X-ray variability suggests that ULXs are binaries with a non-degenerate star, referred to as a donor, transferring mass onto a compact object, the accretor (Basko & Sunyaev 1976; Heil et al. 2009). These objects often dominate the total X-ray emission of their host galaxy. They are generally too bright to be low-mass X-ray binaries (LMXBs; for which L_X ≲ 10³⁷ erg s⁻¹), and too dim when compared to active galactic nuclei (AGNs; for which L_X ≳ 10⁴¹ erg s⁻¹) with their position being off-centre in their host galaxy (see Kaaret et al. 2017, for a review on the observational properties of ULXs).

Eddington luminosity is the limit above which any accretion onto the compact object is stopped by outgoing radiation. In general, ULX luminosities far exceed the Eddington limit of stellar-mass black holes (which is of the order of 10³⁹ erg s⁻¹) provided the emission is assumed to be isotropic. Initial suggestions as to the nature of the accretor include intermediate-mass black holes (IMBHs; with masses ≳100 M_⊙) accreting at sub-Eddington rates (Colbert & Mushotzky 1999; Ebisuzaki et al. 2001; Miller 2006; Maccarone et al. 2007), or stellar-mass black holes (SMBHs) accreting at super-Eddington rates. For the latter, in order not to exceed the Eddington limit locally while still producing the observed luminosity, magnetic fields and geometric X-ray beaming effects have been proposed as an explanation (e.g. King et al. 2001; Christodoulou et al. 2014).

Since the first gravitational wave event, GW150914, detected by LIGO (Abbott et al. 2016), there has been vast renewed interest in studying the formation of double compact objects. Many of the proposed formation channels for these double compact objects predict that the binary will go through one or more phases of super-Eddington accretion onto a compact object. Marchant et al. (2017) explored how some observed ULXs might indeed be progenitors of coalescing double compact object binaries, specifically through the chemically homogeneous channel. In this type of binary evolution, the star would avoid large expansion of the envelope, and therefore the binary would avoid a common-envelope (CE) phase (Mandel & de Mink 2016). Finke & Razzaque (2017) compared observation trends in the number of ULXs per unit of star-formation rate and as a function of metallicity of the host galaxy to the merger rate of binary BHs. These latter authors found that the majority of ULXs could be progenitors of binary BH mergers. Therefore, studying the formation of ULXs could shed light on merging binary compact objects, and vice versa.

In recent decades, debate over the nature of ULXs has focused on whether they are SMBHs or IMBHs, and neutron stars (NSs) have not been considered. However, theoretical studies have extended ULX binary models to include NSs (e.g. King et al. 2001; King 2009). That is until Bachetti et al. (2014) reported the first ever observations of X-ray pulsations in the M82 galaxy with a period of about 1.37 s and a 2.52-day sinusoidal modulation in the ULX X-2. Detection of X-ray pulsations in ULXs suggests the presence of NSs as accretors instead of BHs, at least in a certain fraction of the ULX population where such pulses have been observed. This is because X-ray pulsations are characteristic of accreting NSs with radiation being emitted along their magnetic poles as they rotate about their axes. Therefore, the question of accreting NSs in binaries with relatively massive donors was raised. This discovery was followed by the detection of X-ray pulsations in several other ULXs (Israel et al. 2017a,b; Carpano et al. 2018; Brightman et al. 2018; Sathyaprakash et al. 2019; Zhang et al. 2019; Rodríguez Castillo et al. 2020).

The two most prominent questions regarding the nature of NS ULXs refer to the emission mechanism and the formation channel. For the first question, a number of mechanisms are invoked to explain super-Eddington luminosities. With the apparent extremely high mass-transfer rate, a lot of the transferred mass could be blown away by strong radiation outflows (Shakura & Sunyaev 1973). Begelman et al. (2006) applied this idea to SS433 where observations of massive outflows suggest that the source is a ULX seen from the side. The absorption features associated with these outflows have been observed in some but not all ULXs; see for example Holmberg IX X-1 and NGC 1313 X-1 (Walton et al. 2012). However, this does not invalidate the theory, as beamed X-ray emission would not be visible unless the observer had a direct line of sight to the accreting compact object, down the collimating funnel, in which direction the outflows are limited. Shakura & Sunyaev (1973), Begelman (2002), and Poutanen et al. (2007) explored the idea of the presence of strong optically thick outflows that blow away some part of the disc from where radiation can escape. Shakura & Sunyaev (1973) and Poutanen et al. (2007) also suggested the formation of geometrically thick accretion discs. This structure would cause the emission to be beamed, and therefore the observed isotropic-equivalent luminosity would be much higher than the intrinsic one. King et al. (2001) postulated that based on the assumption of mild beaming, intermediate- and high-mass X-ray binaries (IMXBs and HMXBs) undergoing mass transfer on a thermal timescale would be the best candidates for ULXs. The effect of beaming on the emission has, in general, been explored (King 2009; King et al. 2017; Wiktorowicz et al. 2019). In addition to that, NSs with strong magnetic fields (around 10¹⁴ G) reduce the electron-scattering cross-section and could anchor the infalling matter to accretion columns above the magnetic poles and thereby produce sufficiently high luminosities (Basko & Sunyaev 1976; Brightman et al. 2018). Takahashi & Ohsuga (2017) carried out general relativistic radiation magnetohydrodynamics simulations of super-Eddington accretion onto a non-rotating, magnetised NS and found a spin-up rate of ∼ − 10⁻¹¹ s s⁻¹, which is consistent with observations. In contrast, King & Lasota (2019) suggested that the observed ULX properties are explained by NSs with normal magnetic fields and not by the presence of magnetars. For the specific case of M82 X-2, Lyutikov (2014) suggested the presence of an optically thick accretion disc that acts as a curtain and shields some of the outgoing radiation, thus allowing for super-Eddington luminosities driven by Roche-lobe overflow (RLO). Finally, Weng & Zhang (2011) studied accretion discs in weakly magnetised NSs, finding that at super-Eddington rates, the magnetic field has little effect on the accretion disc.

Despite a lot of research done in the field, the question surrounding how ULXs attain super-Eddington luminosities is still open and an active field of research. Regarding a possible formation channel for NS ULXs, a lot of work has been carried out to address evolution and mass transfer in X-ray binaries. Tauris & Savonije (1999) performed non-conservative mass-transfer calculations of low-mass X-ray binaries with a 1.3 M_⊙ NS and found that the binaries can undergo very high mass-transfer rates (super-Eddington by a factor of ∼10⁴) for donors with deep convective envelopes (mass range 1.6–2.0 M_⊙). An important X-ray source in the context of high mass-transfer rates is Cygnus X-2, which is an X-ray binary containing a NS. Cygnus X-2 has a mass ratio of ≈0.34 (=M_donor/M_acc), so for an estimated accretor mass of about 1.78 M_⊙ the donor has a mass of 0.6 M_⊙ (Casares et al. 1998; Orosz & Kuulkers 1999). King & Ritter (1999) argued that the donor star lost a lot of mass (∼3.0 M_⊙) in an intense mass-transfer phase on a thermal timescale. King & Ritter (1999), Podsiadlowski & Rappaport (2000), Kolb et al. (2000) showed that Cygnus X-2 observations can be explained by case B mass transfer from a donor of mass 3.5 M_⊙; that is, the progenitor for the source was an IMXB.

Contrary to the general understanding of the stability of mass-transfer in X-ray binaries, there is similar evidence in the literature that IMXBs can undergo stable mass transfer with a NS and avoid CE. Tauris et al. (2000) carried out numerical calculations of IMXBs with 2.0–6.0 M_⊙ donor and 1.3 M_⊙ accretor masses using an updated version of the Eggleton code (Eggleton 1971, 1972; Pols et al. 1998). The authors studied the initial parameter space for producing binary millisecond pulsars with a heavy carbon–oxygen (CO) white dwarf companion, and demonstrated for the first time the full stability of IMXBs. King (2009) stated that NSs would have lower accretion rates for the same mass-transfer rates because NSs would have stronger beaming as their Eddington limit would be lower than that of BHs. Calculations similar to those of Tauris et al. (2000) were carried out by Shao & Li (2012) using the Eggleton code. These latter authors studied the initial parameter space for binary pulsars while considering orbital angular momentum losses from gravitational wave radiation, magnetic braking, and mass lost from the system. Tauris et al. (2017) compared the orbital evolution between IMXBs and HMXBs in studying the connection between ULXs and double NS systems and concluded that the orbital period evolution of IMXBs makes them more likely to be NS ULXs than HMXBs.

Fragos et al. (2015) studied the origin of the NS ULX M82 X-2 specifically, by combining parametric population synthesis calculations (using BSE; Hurley et al. 2002) with detailed binary evolution calculations (using MESA; Paxton et al. 2011, 2013). Assuming highly non-conservative mass transfer and that a significant fraction of the mass lost from the binary carries the specific orbital angular momentum of the donor star, these latter authors found that the most probable parameters to form a NS ULX are donors with initial masses in the range of 3.0–8.0 M_⊙, and initial orbital periods of 1.0–3.0 days. Shao & Li (2015) suggested that NS ULXs, in general, might represent a higher contribution to the general ULX population than BH ULXs, and a significant portion of those would be IMXBs. Similarly, Wiktorowicz et al. (2015, 2017), using parametric population synthesis calculations, showed that ULXs are more likely to be NS ULXs than BH ULXs, especially in solar metallicity environments. These latter authors found a typical NS ULX to have a ∼1.0 M_⊙ red giant. However, they found that extreme NS ULXs (L_X ≳ 10⁴² erg s⁻¹) typically have evolved, low-mass, stripped helium-star donors (∼2.0 M_⊙). Wiktorowicz et al. (2019) found that most BH ULXs emit X-rays isotropically, while NS ULXs are generally beamed. Therefore, even though NS ULXs might be intrinsically more numerous than BH ULXs, these latter authors predict that BH accretors dominate the observed ULX population.

Another discussed property of pulsating ULXs is the very high NS spin-up rate (i.e. the rate of change of spin period) in some observed pulsars. NGC 5907 ULX1 shows a change from a spin period of 1.43 s to 1.13 s in about 10 years (Israel et al. 2017b) with an inferred spin-up rate of −8 ± 0.1 × 10⁻¹⁰ s s⁻¹. The spin period of NGC 300 ULX1 went from 31.71 s to 31.54 s in 4 days (Carpano et al. 2018) with a spin-up rate ∼ − 5.56 × 10⁻⁷ s s⁻¹, which is the highest rate observed so far for a NS ULX. However, the rate of change of spin period given by

$$\begin{array}{c}\hfill \dot{P}=\dot{(\frac{1}{\nu })}=-\frac{\dot{\nu }}{{\nu }^2},\end{array}$$(1)

where ν is the NS spin frequency, does not directly reflect the mass accretion rate. Rather, the rate of change of ν is roughly proportional to the mass accretion rate (for instance see Eq. (2) in King et al. 2017). Nevertheless, we continue to mention Ṗ values when talking about observations (Table 1) as they are the ones most often reported in the literature.

High rates of frequency increase could suggest efficient spin up from very high accretion rates onto the NS. This is because high mass accretion provides the NS with enough torque to spin it up to such high rates (Ritter & King 2001; Tauris & van den Heuvel 2006). However, there is a caveat in calculating accretion rates from spin-up frequency as they might not represent the secular average accretion rate over an evolutionary timescale of the donor star. In some cases, the extremely high spin-up rate would grossly overestimate the amount of matter accreted by the NS. For example, M82 X-2 showed a high spin-up rate of −2 × 10⁻¹⁰ s s⁻¹ when X-ray pulses were first discovered (Bachetti et al. 2014). However, later, Bachetti et al. (2020) observed an average spin down of −5 × 10⁻¹¹ Hz s⁻¹ over a period of two years. These latter authors suggested that the source might be close to spin equilibrium and is alternating between phases of spin up and spin down. King et al. (2017) estimated the spin-up timescale for three pulsating ULXs and concluded that we observe them close to equilibrium. Chashkina et al. (2019) performed semi-analytical calculations for accretion onto a magnetised NS using the NS being close to spin equilibrium as a boundary condition. The reason for this is that even though $\dot{\nu }$ gives an estimate of the instantaneous accretion rate, it should not be compared directly to long-term average estimates that binary evolution models are giving, as many of these systems might be close to their spin equilibrium.

In our work, we do not investigate X-ray pulses and the super-Eddington emission mechanism. Instead, we focus on ULX formation and long-term evolution. Therefore, in the entirety of this study, we refer to LMXBs/IMXBs with NS accretors that drive super-Eddington mass-transfer rates as NS ULXs. Pulsating ULXs are a subset of NS ULXs as pulsations are not a necessary outcome of super-Eddington mass transfer but rather they are a product of the presence of a relatively strongly magnetised NS. In this paper, we investigate how these NS ULXs could be formed and try to explain the physical properties involved using numerical computations. In doing so, we study how the stability of binaries is affected by spin-orbit coupling, and by a higher accretor mass (2.0 M_⊙ instead of a typical NS mass of ∼1.3 M_⊙). We also assume that there is no precession of the accretion disc or absorption of X-ray flux by optically thick material around the source that might cause the luminosity to vary. In the following section (Sect. 2), we summarise the properties of the currently observed pulsating ULX sample. In Sect. 3, we discuss the numerical methods and physics employed for the simulations, while in Sect. 4 we present the results from our simulations, highlighting the allowed initial parameter space for NS ULX formation and the properties of the formed population. Section 5 discusses the observed NS ULXs in the context of our results and how the angular momentum exchange between spin and orbit affects the result through tides. Finally, we end with concluding remarks in Sect. 6.

2. Currently observed pulsating ULX sample

In this section, we discuss some of the observed and predicted parameters for the NS ULXs present so far in the literature. In most cases, the observables are not well constrained, and therefore large uncertainties are involved and assumptions about the physical properties have been made. In the entirety of this paper we refer to the mass of the donor as M_donor, that of the accretor as M_acc, the orbital separation as a, and the orbital period as P_orb. For our work, the most important pulsating ULX is M82 X-2 as it has the most well-constrained observational parameters, even though we comment on other pulsating ULXs observed as well. We summarise all the relevant properties of the currently known sample of pulsating ULXs in Table 1.

2.1. M82 X-2

This source is one of the most studied NS ULXs. It is observed in the core of M82 galaxy and was discovered to show X-ray pulsations by Bachetti et al. (2014). It has a peak X-ray luminosity of 1.8 × 10⁴⁰ erg s⁻¹. Using the observed orbital period of 2.52 days and inclination of < 60°, the binary mass function for M82 X-2 was calculated as 2.1 M_⊙. Assuming that the NS mass is 1.4 M_⊙, the donor mass was estimated to be ≳5.2 M_⊙.

2.2. NGC 7793 P13

This ULX source appears to be a HMXB containing an accreting NS in galaxy NGC 7793, with a luminosity of about 5.0 × 10³⁹ erg s⁻¹. Earlier, Motch et al. (2011) identified the donor in the then-not-discovered ULX system as a late B-type super-giant star with a mass of 18.0 M_⊙ <  M_donor <  23.0 M_⊙. Motch et al. (2014) estimated the orbital period to be 64 days using optical modulation in observations. The NS was found to be present after pulsations were detected from the source (Israel et al. 2017a; Fürst et al. 2016).

2.3. NGC 5907 ULX1

Discovered to have X-ray pulsations by Israel et al. (2017b), this ULX is observed in the galaxy NGC 5907. Its X-ray luminosity, at around 10⁴¹ erg s⁻¹, makes it the most luminous ULX discovered to date. This source does not have as well-constrained parameters as M82 X-2. However, some estimates could be made. The authors derived constraints on the orbital period and arrive at a value of $5.3_{-0.9}^{+2.0}$ days using a projected semi-major axis of $asini=2.5_{-0.8}^{+4.3}$ lt-s (light-seconds). They also compared HMXBs, LMXBs, and IMXBs as possible explanations and found them all to be consistent with the upper limits derived from HST images (Sutton et al. 2013).

2.4. NGC 300 ULX1

This X-ray source, located in the galaxy NGC 300, was recently discovered to show pulsations (Carpano et al. 2018). The authors assumed an orbital period in the range 1–3 days and used a projected semi-major axis of ≲4.6 × 10⁻⁵ lt-s and arrived at a binary mass function of 0.8 × 10⁻³ M_⊙. The X-ray luminosity for ULX1 was found to be 4.7 × 10³⁹ erg s⁻¹. Heida et al. (2019) identified the donor as a red super-giant star (≳8.0 M_⊙). They reported that the observed effective temperature and luminosity of the donor is consistent with evolutionary tracks of singles stars of initial masses of 8.0–10.0 M_⊙. Vasilopoulos et al. (2019) observed that the flux of the source decreased by a factor of ∼50 over a few months, while its spin-up remained constant at 4.0 × 10⁻¹⁰ Hz s⁻¹. This could result from absorption from optically thick material nearby (perhaps from outflows) or from a precessing accretion disc.

2.5. M51 ULX-8

Brightman et al. (2018) identified M51 ULX-8 as a NS ULX on the basis of the detection of cyclotron resonance scattering features instead of direct observations of X-ray pulses. The cyclotron features observed can be translated into magnetic field strengths, which in this case correspond to that of a highly magnetized pulsar (∼10¹⁵ G).

2.6. NGC 1313 X-2

Sathyaprakash et al. (2019) discovered X-ray pulsations in the ULX source NGC 1313 X-2 with a pulse period of ∼1.5 s. Grisé et al. (2008) estimated the mass of the donor star to be ≲12.0 M_⊙ assuming it to be part of a metal-poor star cluster with an age of 2.0 Myr. Pakull & Grisé (2008) estimated the ULX lifetime using observations of the large bubble nebulae around the source, which is ∼1.0 Myr. Based on studies showing that very few NS ULXs would last as long as 1.0 Myr (e.g. Wiktorowicz et al. 2017), Sathyaprakash et al. (2019) suggested that NGC 1313 X-2 is in its last stages of mass transfer.

2.7. Swift J0243.6+6124

Discovered by Kennea et al. (2017) as a transient source, Swift J0243.6+6124 was known to have an accreting NS with a spin period of 9.86 s. Doroshenko et al. (2018) estimated an orbital period of 28.3 days and semi-major axis of $asini={140}_{-3}^{+3}$ lt-s. Zhang et al. (2019) studied Insight-HXMT data for this source and estimated its magnetic field as ∼10¹³ G and an X-ray luminosity of > 10³⁹ erg s⁻¹, confirming the source to be the first Galactic ULX with an NS accretor. Furthermore, these latter authors calculated a high spin-up rate of −2.2 × 10⁻⁸ s s⁻¹. Tao et al. (2019) reported on the spectral behaviour of this source and suggested that if it were located in an external galaxy it would have a similar appearance to the other pulsating ULXs observed.

2.8. M51 ULX-7

Rodríguez Castillo et al. (2020) discovered 2.8 s X-ray pulses in observations of ULX-7 in galaxy M51, therefore finding another pulsating ULX following ULX-8 in the same galaxy. The secular spin-up rate was measured as −10⁻⁹ s s⁻¹ and a variable X-ray luminosity in the range of 10³⁹–10⁴⁰ erg s⁻¹. The authors used the projected semi-major axis of 28.0 lt-s and an assumed accretor mass of 1.4 M_⊙ to infer a donor of mass > 8.0–13.0 M_⊙. Vasilopoulos et al. (2020) studied the X-ray light curves of the source and estimated a magnetic field of the NS (rotating near spin equilibrium) of 2.0–7.0 × 10¹³ G. They also assumed that the NS is freely precessing and estimated the magnetic field to be 3.0–4.0 × 10¹³ G, agreeing with their previous estimate.

3. Numerical tools and calculations

In this section we discuss the numerical code used for the binary evolution calculations along with the adopted model parameters and the code modifications that we introduced.

3.1. Numerical stellar evolution code and progenitor binary

To simulate the evolution of the binaries we use MESA (version 10108; MESASDK version 20180127) which is a stellar structure and binary evolution code developed by Paxton et al. (2011, 2013, 2015, 2018, 2019).

All LMXBs/IMXBs are calculated as systems with an initially zero-age main-sequence (ZAMS) donor and a point-mass NS accretor. It is assumed that the NS was formed during a previous evolutionary stage which we do not study, and the binary survived a possible NS natal kick from the supernova explosion that formed the NS. We compute a grid of models spanning 0.92–8.0 M_⊙ in initial donor mass, and 0.5–100 days in initial orbital period. For reference, initial parameters refer to the orbital parameters at the onset of RLO. The calculation of the mass-transfer rate during RLO is done implicitly using the scheme proposed by Kolb & Ritter (1990).

We use 1.3 M_⊙ and 2.0 M_⊙ for the mass of the accretor; 1.3 M_⊙ is close to the post-supernova Chandrasekhar mass limit for the formation of a NS and 2.0 M_⊙ is on the high-mass end of the NS mass distribution (Lattimer & Prakash 2010; Antoniadis et al. 2013; Rezzolla et al. 2018). The NS companion to the pulsar J0453+1559 has a mass of $1.{174}_{0.004}^{0.004}{M}_{\odot }$ (Martinez et al. 2015)¹ and therefore NS masses below 1.3 M_⊙ have indeed been observed. However, we take 1.3 M_⊙ as the standard NS mass.

For the radiative efficiency of the accretion onto the NS with initial mass $M_{\text{acc}}^i$ (i.e. the release of gravitational energy of the infalling material in the form of radiation; in units of rest-mass energy) we use,

$$\begin{array}{c}\hfill \eta =\frac{G{M}_{\mathrm{acc}}^i}{c^2{R}_{\mathrm{acc}}},\end{array}$$(2)

where c is the speed of light and R_acc is the NS radius which we take as 11.0 km. For the Eddington limit of an accretor with initial mass $M_{\text{acc}}^i$,

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=\frac{4\pi G{M}_{\mathrm{acc}}^{i}c}{\kappa },\end{array}$$(3)

where, κ is the opacity. Using Eq. (2) and simplifying Eq. (3) (for accretion of pure ionised hydrogen), we get

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Edd}}=1.5\times {10}^{-8}(\frac{M_{\mathrm{acc}}^i}{1.3{\mathrm{M}}_{\odot }}){\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}.\end{array}$$(4)

We fix the Eddington limit and the radiative efficiency to initial values as they would not change significantly from the amount of mass that the NSs accrete in our grids. If the mass-transfer rate goes beyond this value, the radiation pressure will prevent any excess material from being accreted. Mass from the donor is transferred conservatively to the accretor and a fraction of this transferred mass is lost from the vicinity of the accretor as an isotropic fast wind or jet with the specific angular momentum of the accretor (see Tauris & van den Heuvel 2006, for a detailed explanation). The efficiency of accretion (ϵ) by the NS is ϵ = 1 − (α + β + δ), where α is the fractional mass lost directly from the donor, β is the fractional mass lost from the vicinity of the accretor, and δ is the fractional mass lost from a circumbinary toroid. We take the values, α = δ = 0, and β = max {0.7, 1 − Ṁ_Edd/Ṁ_donor}.

We consider orbital angular momentum (J_orb) losses via gravitational wave radiation, spin-orbit coupling due to tidal effects (Sect. 3.3), and mass lost from the system. We also include effects due to magnetic braking following the prescription by Rappaport et al. (1983) for donor masses that develop an outer convective envelope at any point. We take the eccentricity to be negligible, as tidal forces would circularise the orbit of a semi-detached binary with a giant star on a relatively short timescale of 10⁴ years (Verbunt & Phinney 1995). This is orders of magnitude shorter than the main sequence lifetime of intermediate-mass stars which are of the order of 10⁸–10⁹ years (main sequence lifetime of low-mass stars would be even longer). Furthermore, tidal forces aim to synchronise the stars with the orbit. We assume the orbit is synchronised by the time RLO begins since the tidal forces would cause the stars to synchronise with the orbit on a relatively short timescale (Huang & Zeng 2000).

We assume solar metallicity, that is Z_⊙ = 0.02, and that any layer in the donor interior is stable against convection if the Ledoux criteria for convection is fulfilled (Ledoux 1947). At the edges of convective zones, we account for overshooting because the convective material slightly enters non-convective zones due to inertia. To describe overshooting we follow the exponential overshooting efficiencies used in the MIST models for low- and intermediate-mass stars, which are f_ov, core = 0.0160 in the core calibrated from properties of the open cluster M67 and f_ov, en = 0.0174 in the envelope calibrated from solar properties (Dotter 2016; Choi et al. 2016). For stellar winds, we use the cool red giant branch wind scheme described by Reimers (1975) with a scaling factor of 0.1. In cases where we get a stripped helium (He) star, we use the prescription for Wolf-Rayet stars by Nugis & Lamers (2000), included in MESA under the Dutch hot wind scheme.

The mass-transfer calculations are carried out until one of the following conditions are met: (i) the donor forms a white dwarf (WD), (ii) the age of the donor star exceeds the Hubble time, (iii) the radius of the donor star extends so far beyond its Roche lobe that L₂ overflow is initiated and the system becomes dynamically unstable (see Sect. 3.4), or (iv) the number of computational steps exceeds a limit of 300 000 (this value was chosen based on previous grid runs). We include condition (iv) for those systems where MESA runs into converging problems and cannot find a solution. This happens for only two types of binaries in our numerical calculations. In the first case, the mass transfer cannot properly remove the last bit of the envelope from the donor. In the second case, MESA is not able to solve for the donor radius which extends quite far beyond the Roche lobe but not enough to trigger the condition of L₂ overflow. Both these cases occur at the end of the mass-transfer phase, and therefore we accept that the binary was stable until the end of its evolution.

3.2. Super-Eddington accretion onto the NS

For non-conservative super-Eddington mass transfer the amount of mass that is accreted is less than that transferred per unit time (0.3× mass transferred, as per our assumptions) until the Eddington limit prevents a further increase in accretion rate. Cygnus X-2 is an example of a system observed to have survived super-Eddington mass transfer while presumably having accreted comparatively less. This source suggested that high mass transfer onto a NS (or BH) could avoid a CE phase with the NS ejecting most of the mass that was transferred at super-Eddington rates (King & Ritter 1999; King & Begelman 1999; Podsiadlowski & Rappaport 2000; King et al. 2001).

Shakura & Sunyaev (1973) studied the observational characteristics of accretion discs in sub-Eddington and super-Eddington regimes of mass transfer in the case of black-hole binaries. In this picture, as the mass transfer approaches the Eddington limit, the structure of the disc changes from a slim disc to a disc with an inner geometrically thick component and an outer thin component. As matter that is transferred to the accretion disc at super-Eddington rates moves radially inwards (transporting angular momentum outwards), strong outflows begin at a certain radius which remove a fraction of the matter, thereby also taking away excess angular momentum (see Vinokurov et al. 2013, for an application of the super-Eddington disc model). The spherisation radius is defined as the radius at which the accretion luminosity first reaches the Eddington limit and strong outflows begin, and one can approximate it as

$$\begin{array}{c}\hfill R_{\mathrm{sph}}={\dot{M}}_{\mathrm{donor}}\frac{G{M}_{\mathrm{acc}}}{L_{\mathrm{Edd}}}·\end{array}$$(5)

Outside R_sph the disc emits X-rays with luminosity L_Edd which depends on the accretion rate following the equation,

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=\eta {\dot{M}}_{\mathrm{Edd}}{c}^2·\end{array}$$(6)

Inside R_sph, the outflowing matter has a velocity which depends on the difference between inward gravity and outward radiation pressure. The velocity of the outflow increases inward which in turn decreases the mass-accretion rate at each radius. This keeps the disc locally Eddington limited. The mass-accretion rate within R_sph at each point in the disc can be described by the following equation:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}^{\mathrm{local}}(R)={\dot{M}}_{\mathrm{donor}}\frac{R}{R_{\mathrm{sph}}}·\end{array}$$(7)

Since the radiation pressure is balanced by the gravitational pressure at each point, the accretion disc inside the spherisation radius has a thickness of the order of the distance from the accretor and emits X-rays with luminosity L_Edd × ln ṁ (where ṁ ≡ Ṁ_donor/Ṁ_Edd; Shakura & Sunyaev 1973). For ṁ > 1, the total radiated luminosity from the accretor can exceed the Eddington limit by

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=L_{\mathrm{Edd}}(1+ln\dot{m}).\end{array}$$(8)

For highly super-Eddington mass-transfer rates another effect could come into play because of the geometrically thick accretion disc. A narrow funnel forms along the rotation axis of the accretor from where radiation can escape as a collimated jet. The observed isotropic-equivalent accretion luminosity as described by King et al. (2001) and King (2009) is then as follows:

$$\begin{array}{c}\hfill L_{\mathrm{acc}}^{\mathrm{iso}}=\frac{L_{\mathrm{Edd}}}{b}(1+ln\dot{m}),\end{array}$$(9)

where b is the beaming factor describing the amount of collimation to the outgoing radiation. The approximated value of b is

$$\begin{array}{c}\hfill b=\{\begin{array}{cc}\frac{73}{{\dot{m}}^2},\hfill & \text{if}\dot{m}>8.5,\\ 1,\hfill & \text{otherwise}.\end{array}\end{array}$$(10)

Because the beaming factor is an approximation, we apply an upper limit (∼10⁴² erg s⁻¹) in calculating the isotropic-equivalent accretion luminosities so that we do not get unphysically high values.

3.3. Spin-orbit coupling and synchronisation timescales

As mentioned before, tidal forces synchronise the stars with the orbit on a timescale which is relatively short. Therefore, whenever there is a change in either the spin angular momentum or the orbital angular momentum, tides will work to synchronise the system again. This action of tides on the orbit may affect the stability and evolution of the binary (Tauris et al. 2001).

When mass is lost from the donor star during mass transfer it also removes spin angular momentum from the star which is supplied to the orbit. If the spin angular momentum removed from the donor and returned to the orbit is non-negligible, then it has a widening effect on the orbit. This is followed by mass loss from the accretor’s vicinity, with the mass lost carrying away the specific angular momentum of the accretor. This competing effect tends to shrink the orbit. Mass that is leaving the system from the accretor’s vicinity carries with it relatively high specific orbital angular momentum (when the accretor is the less massive binary component), having a shrinking effect on the orbit. The donor then becomes sub-synchronous with the orbit and angular momentum has to be transferred from the orbit to the star in order to spin it up, causing the orbit to shrink even further. The interplay between the orbital shrinking and widening effects can reveal how spin-orbit coupling affects the stability of mass transfer. In the presence of strong winds, there is additional loss of angular momentum from the donor. However, for stars with masses of 0.92–8.0 M_⊙, stellar winds are too weak to cause significant orbital change.

The tidal synchronisation timescale is defined as follows (Zahn 1977; Hut 1981):

$$\begin{array}{c}\hfill \frac{1}{T_{\mathrm{sync}}}=3\frac{K}{T}(\frac{q}{r_{\mathrm{g}}}{)}^2(\frac{R_{\mathrm{donor}}}{a}{)}^6,\end{array}$$(11)

where q is the mass ratio (we define q ≡ M_acc/M_donor) and r_g is the gyration radius of the star. K/T is the spin-orbit coupling parameter and is described in two ways for a star with an outer radiative envelope as (K/T)_rad, and a star with an outer convective envelope as (K/T)_conv. For the former case we use,

$$\begin{array}{cc}& (\frac{K}{T}{)}_{\mathrm{rad}}=(\frac{G{M}_{\mathrm{donor}}}{R_{\mathrm{donor}}^3}{)}^{1/2}{(1+q)}^{5/6}{E}_2(\frac{R_{\mathrm{donor}}}{a}{)}^{5/2},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill & E_2={10}^{-0.42}(\frac{R_{\mathrm{conv}}}{R_{\mathrm{donor}}}{)}^{7.5},\hfill \end{array}$$(13)

where R_conv is the radius of the convective core. The value for E₂ is computed from fitting formulae for H-rich stars derived by Qin et al. (2018). For (K/T)_conv, the definition is taken from Hurley et al. (2002) to be

$$\begin{array}{c}\hfill (\frac{K}{T}{)}_{\mathrm{conv}}=\frac{2}{21}\frac{f_{\mathrm{conv}}}{{\tau }_{\mathrm{conv}}}\frac{M_{\mathrm{env}}}{M_{\mathrm{donor}}}{\mathrm{yr}}^{-1},\end{array}$$(14)

where f_conv = min {1.0, (P_tidal/(2τ_conv))²} is a numerical factor and P_tidal is the tidal pumping time-scale defined by |1/P_orb − 1/P_spin|. Here, τ_conv ≈ (MR²/L)^1/3 is the eddy turnover timescale in units of years (Rasio et al. 1996) and M_env is the convective envelope mass. More details about our applied orbital angular momentum evolution, with spin-orbit coupling including tides, are discussed in Sect. 5.3 and Appendix A.

3.4. Angular momentum loss from the second Lagrangian point

Lagrangian points are equilibrium points in space where the gravitational and centrifugal forces in the system cancel each other out. L₁, L₂, and L₃ are unstable equilibrium points from where a test particle, upon small displacement, would move further away. Figure 1 shows these unstable Lagrangian points for a system with a mass ratio of q = 0.26.

Fig. 1. Equipotential lines of the Roche potential for a binary consisting of stars with mass ratio q = 0.26 and binary separation a. The equipotential lines passing through Lagrangian points L1, L2, and L3 are shown.

In most cases of X-ray binaries, analysis has been done for mass transfer via L₁. L₁ lies in between the two stars (hence also known as the inner Lagrangian point) and the equipotential surface passing L₁ is known as the Roche lobe. When a star fills its Roche lobe, any material that crosses L₁ from one star will fall towards the other star. This transfer of matter either decreases the radius (for radiative envelopes) of the donor or increases it (for convective envelopes). In some cases of extreme binary mass ratio, the RLO might not be enough to provide efficient mass-transfer rate and the donor might extend far beyond its Roche lobe to reach the equipotential surface passing through L₂. However, contrary to when mass passes through L₁, the mass that crosses L₂ takes away a large amount of angular momentum from the binary. Once the outer layers of the donor reach L₂ (or the donor obtains a volume equivalent to that of the equipotential lobe passing through L₂; see below) it is expected that the binary orbit will shrink rapidly. We consider this the onset of dynamical instability. For illustration, the L₂ potential surface is shown in Fig. 1. It is the peanut-shaped surface enclosing both binary mass components and passing through the point L₂.

Eggleton (1983) calculated the stellar radius needed in order to initiate mass transfer from the inner Lagrangian point of a binary by calculating the radius of a sphere that will have the same volume as the Roche lobe. This radius for the donor star is referred to as the Roche-lobe radius (R_donor, L1). We take a similar approach to quantify overflow from L₂. In cases where mass transfer via the L₁ point is not sufficient to keep the donor star confined within its Roche lobe we assume that the expanding star needs to fill the entire volume enclosed by the L₂ equipotential surface before the onset of dynamical instability. R_L2 is the volume equivalent radius for the equipotential surface passing through the L₂ point.

Another possibility of mass loss from L₂ point occurs when the radius of the donor star reaches the point L₂ before the donor volume overfills the L₂ potential surface. This case applies only when q ≥ 1 (i.e. M_donor ≤ M_acc) as the point L₂ is much closer to the donor. We refer to the distance between the centre of the donor and L₂ as D_L2. This case is illustrated in Fig. 2 which shows the equipotential surfaces with mass ratio, q = 1.85. In our simulations we assume that binaries experience stable mass transfer when the donor star does not cross any of the two limits discussed above at any point during each evolution (i.e. for stable RLO: R_donor <  min {R_L2,  D_L2} for the entire binary evolution).

Fig. 2. Equipotential lines of the Roche potential for a binary consisting of stars with mass ratio q = 1.85 and binary separation a to illustrate the swapping of the positions of L2 and L3 when Mdonor ≤ Macc (compared to Fig. 1).

To calculate the L₂ volume via numerical integration, we begin by finding the Lagrange points (L₁, L₂ and L₃) for a particular mass ratio. Along the axis joining the centres of the M_donor and M_acc, the entire volume is assumed to be a summation of thin discs. The volume of each disc slice is calculated going from one boundary end to the other using the boundaries of the L₂ equipotential surface, and the subsequent disc volumes are added together to cover the entire volume. In these calculations we consider two mass ratio regimes, q <  1 and q ≥ 1.

We use both R_L2 and D_L2 normalised to R_donor, L1 in order to remove the dependence on the binary orbital separation. Once these radii are calculated we find a fit of R_L2/R_donor, L1 and D_L2/R_donor, L1 on the mass ratio (q = M_acc/M_donor) of the binary. For q <  1, we find that the R_L2 and D_L2 follow monotonically increasing trends toward q = 1 (shown in Fig. 3 in the left panel), which is fitted by the following functions,

$$\begin{array}{c}\hfill \frac{R_{{\mathrm{L}}_2}(q<1)}{R_{\mathrm{donor},{\mathrm{L}}_1}}=0.784q^{1.05}{e}^{-0.188q}+1.004,\end{array}$$(15)

$$\begin{array}{c}\hfill \frac{D_{{\mathrm{L}}_2}(q<1)}{R_{\mathrm{donor},{\mathrm{L}}_1}}=3.334q^{0.514}{e}^{-0.052q}+1.308.\end{array}$$(16)

Fig. 3. Our results form the volume integration of L2 equipotential surfaces with respect to different mass ratios q. All distance values are presented in units of RL1. Here, RL2 is the radius of a sphere with volume equal to that of the L2 equipotential surface (solid orange line). DL2 is the distance of the L2 point from Mdonor (solid blue line). The dashed black line shows the orbital separation in units of RL1. Systems where the donor’s radius exceeds any of these limiting radii are considered as undergoing dynamical instability.

Here, R_donor, L1 is also calculated from the volume of the L₁ equipotential surface using the method described above and the result is consistent with the calculations from Eggleton (1983).

For the second mass ratio regime, we calculate R_L2(q ≥ 1) using

$$\begin{array}{c}\hfill R_{{\mathrm{L}}_2}(q)=R_{{\mathrm{L}}_2}(\frac{1}{q})\forall q>0.\end{array}$$(17)

We find D_L2(q ≥ 1) using the distance between L₂ and M_acc from the case q <  1. As seen in Fig. 3 on the right, there is a sudden jump in the values of D_L2: when crossing q = 1, it goes from around 4.45 for q = 0.997 to around 1.82 for q = 1.003. We fitted functions to the calculated R_L2 and D_L2 values as follows,

$$\begin{array}{c}\hfill \frac{R_{{\mathrm{L}}_2}(q\ge 1)}{R_{\mathrm{donor},{\mathrm{L}}_1}}=0.290q^{0.829}{e}^{-0.016q}+1.362,\end{array}$$(18)

$$\begin{array}{c}\hfill \frac{D_{{\mathrm{L}}_2}(q\ge 1)}{R_{\mathrm{donor},{\mathrm{L}}_1}}=-0.040q^{0.866}{e}^{-0.040q}+1.883.\end{array}$$(19)

The relative errors between our calculations and the corresponding fits are less than 1%.

Volume equivalent L₂ radii calculations have also been done by Marchant et al. (2016) for q <  1 but with an entirely different approach as these latter authors considered the case of overcontact binaries. They calculated the L₂ volume assuming that both the stars expand and fill their respective L₂ sub-volume. In order to test our numerical volume-integrating scheme against that of Marchant et al. (2016), we split the L₂ volume at the L₁ point and compared the resulting calculations to their R_L2 radii, finding good agreement. However, we should stress again that the approach by Marchant et al. (2016) is only applicable to overcontact binaries. Our work was followed by Ge et al. (2020) who, in a different context, derived fits to the volume equivalent L₂ radius using a similar approach.

4. Results

We explore the evolution of LMXBs/IMXBs with different initial conditions, taking into account the physics described in the earlier sections. In both our grids (for NS masses of 1.3 M_⊙ and 2.0 M_⊙ respectively) the binaries that interacted via mass transfer can undergo either stable or unstable mass transfer, excluding the systems with $P_{\text{orb}}^i\lesssim 0.50$ days where the donor star already overflows its Roche lobe at ZAMS which we do not further evolve. The superscript i stands for initial values which corresponds to the orbital parameters at the onset of RLO. We flag systems as ‘stable’ when either the donor has detached from its Roche lobe at the end of mass transfer or the hydrogen in the outer layer has been almost completely removed (remaining hydrogen in the outer layer < 0.005 M_⊙). Donors in stable binaries formed a WD at the end of the mass-transfer sequence resulting in a neutron star–white dwarf (NS–WD) binary. The ‘unstable’ binaries underwent the onset of L₂ overflow.

In Fig. 4 we present this allowed initial parameter space for LMXBs/IMXBs to undergo stable mass transfer for grids containing 1.3 M_⊙ (orange squares) and 2.0 M_⊙ NS accretors (blue squares) along with the unstable sequences (grey squares). The dark red squares towards the left are the systems that never initiated RLO. The general shape of the stable region resembles the work done by Tauris et al. (2000) where they explored the allowed parameter space to form binary millisecond pulsars while avoiding a CE phase. The different dashed white lines separate the grids based on the type of mass-transfer phase that the system underwent. Case A is when the donor is on the main sequence at the onset of mass loss, that is, it is burning hydrogen in the core. Case B is when the donor has exhausted H in its core and H-shell burning phase (post-main sequence). The threshold between the two cases A and B (lower dashed white line in figure) depends more on the initial orbital period than the donor mass; the limit for case A being in the range 2.0–2.8 days for both grids. The middle dashed white line separates two subsets of case B RLO: early case B and late case B. This threshold depends almost linearly on both the initial orbital period and the initial donor mass. At RLO onset, if the donor has a radiative envelope it is termed early case B and if the donor has developed a convective envelope it is termed late case B RLO. The upper dashed white curve encloses systems that undergo case C RLO, which means the donor has exhausted He in its core at the onset of RLO.

Fig. 4. Allowed initial parameter space for LMXBs/IMXBs to undergo stable mass transfer with 1.3 M⊙ accretors (orange squares) and 2.0 M⊙ accretors (blue squares). Grey squares correspond to systems that encountered dynamical instability. The dark red squares towards the left edge of the figure correspond to systems that never initiated RLO. The lower dashed white line separates systems undergoing case A mass transfer from those undergoing case B mass transfer (same boundary in both grids). The middle dotted white line separates early case B (where the donor has yet to form a deep convective envelope) from late case B RLO (where donor has formed a deep convective envelope at onset of RLO). The upper dashed white curve encloses systems that undergo case C mass transfer (same boundary in both grids). Green stars correspond to those systems that undergo a second mass-transfer phase from a stripped Helium-giant star, i.e. case BB RLO. The three red crosses in case BB correspond to systems that terminated due to numerical issues.

The stability region increases for higher accretor mass to include higher donor masses and orbital periods. A similar effect of increase in the parameter space with increasing the accretor mass was obtained by Shao & Li (2012) for the formation of recycled pulsars from IMXBs and LMXBs. The overall shape of the parameter space for both the grids depends a lot on the structure of the donor envelope and the response of the donor to mass loss. A radiative envelope would shrink on mass loss and contribute to the stability of the binary while a convective envelope would expand rapidly on mass loss and make the system increasingly unstable. King & Begelman (1999) showed that mass transfer on a thermal timescale would avoid the CE phase in a binary as long as the envelope is mostly radiative. Case A and early case B have radiative envelopes, with binaries of initial mass ratio greater than ∼0.28, undergoing stable mass-transfer depending on the initial orbital period. In contrast, late case B and case C have convective envelopes at RLO with the stability region being defined by a fixed critical mass ratio which we find to be at around ∼0.5. Istrate et al. (2014) studied the binary evolution in close LMXBs and also found an increase in the depth of the convective envelope with increasing $P_{\text{orb}}^i$ (Fig. 4 in their paper).

In Fig. 4, green stars correspond to systems that go through case BB RLO, which is when a binary, after having lost its hydrogen envelope in a case B mass-transfer phase, detaches and evolves as a stripped helium star and initiates a second RLO phase during the helium-shell burning stage. This subset of case B occurs for only a small part of the parameter space (and only for the grid with a 2.0 M_⊙ accretor). This is because low-mass helium stars (< 0.8 M_⊙) do not expand by any significant amount (e.g. Heusgen 2016; Kruckow et al. 2018) and thus only donor stars ≳5.8 M_⊙ leave behind stripped helium stars massive enough to eventually lead to case BB RLO. The three red crosses in case BB correspond to systems that terminated due to numerical issues.

During the RLO phase, all LMXBs/IMXBs are expected to be bright X-ray binaries, and in many cases the mass-transfer rate from the donor star to the vicinity of the NS can exceed the Eddington limit (L_Edd) significantly. As described in Sect. 3.2, we only allow accretion onto the NS up to the Eddington limit, but we do consider the transition from a thin accretion disc to an inner thick disc model when the mass-transfer rate supplied from the donor star exceeds the Eddington limit (Shakura & Sunyaev 1973), as well as the geometric beaming model proposed by King et al. (2001).

Figure 5 shows an example of a stable system with $M_{\text{donor}}^i=2.77M_{\odot }$, $M_{\text{acc}}^i=1.30M_{\odot }$, $P_{\text{orb}}^i=5.38$ days, where we demonstrate the estimated observed X-ray luminosity under different assumptions. The solid green curve is the accretion luminosity corresponding to the amount of mass accreted by the NS (i.e. capped at exactly the Eddington limit), assuming a radiative efficiency following Eq. (2). The solid orange curve is the accretion luminosity corresponding to a super-Eddington accretion disc following Eq. (8), where although the Eddington limit is locally satisfied at every point in the disc, the integrated luminosity of the disc can exceed the Eddington limit by a small factor. Finally, the solid blue curve is the isotropic-equivalent accretion luminosity corresponding to beamed super-Eddington emission following Eq. (9) where the estimated observed luminosity can exceed the Eddington limit by up to a few orders of magnitude. For comparison, we mark with horizontal lines the luminosity of ULXs corresponding to L_Edd (solid black line), 10 L_Edd (dashed black line), 100 L_Edd (dotted black line), and 1000 L_Edd (dot-dashed black line). For simplicity, we use Eq. (4) which is fixed for an initial accretor mass, because the amount of mass accreted does not change the accretion luminosity significantly. We note that since the prescription used for beamed emission (Sect. 3.2) is an approximation, we fix an upper limit for the calculated isotropic-equivalent accretion luminosities at 10⁴² erg s⁻¹ in order to avoid unphysically high values. In any case, the time that our binaries spend at such high luminosities combined with the inferred very small beaming factors make binaries on that phase effectively non-detectable. In the remainder of the paper, we use the beamed model (i.e. equivalent to the blue curve based on Eqs. (9) and (10)) as our estimate of the observed X-ray luminosity of the binaries.

Fig. 5. Estimated observed X-ray luminosity for a binary with $M_{\text{donor}}^i=2.77M_{\odot }$, $M_{\text{acc}}^i=1.30M_{\odot }$, and $P_{\text{orb}}^i=5.38$ days under different assumptions. The solid green curve is the accretion luminosity corresponding to the mass-accretion rate onto the NS, assuming a fixed radiative efficiency (η = 0.1). The solid orange curve is the accretion luminosity corresponding to a super-Eddington accretion disc following the Shakura & Sunyaev (1973) disc model given by Eq. (8). The solid blue curve is the isotropic-equivalent accretion luminosity corresponding to beamed super-Eddington emission following the King et al. (2001) geometric beaming model (limited at 1042 erg s−1) given by Eqs. (9) and (10). Four reference luminosity values, corresponding to LEdd (solid black line), 10 LEdd (dashed black line), 100 LEdd (dotted black line), and 1000 LEdd (dot-dashed black line), are also plotted. The initial properties of this binary are highlighted with a magenta star in Fig. 6.

Figure 6 shows the time-averaged isotropic-equivalent accretion luminosities with respect to the initial parameters in both the grids. The magenta star is the binary shown in Fig. 5. For accreting NSs of either 1.3 M_⊙ or 2.0 M_⊙, we find time-averaged isotropic-equivalent X-ray luminosities of $⟨L_{\mathrm{acc}}^{\mathrm{iso}}⟩\simeq {10}^{36}-{10}^{41}{\mathrm{erg}\mathrm{s}}^{-1}$ (in some cases even up to 10⁴² erg s⁻¹). These luminosities were calculated following Eq. (9) and averaged over the entire RLO lifetime for systems that reached instantaneous luminosity above 10 L_Edd. The stable systems are enclosed within the solid black boundary (for this and all subsequent figures). Comparing Fig. 4 with Fig. 6, even the systems that undergo dynamical instability (L₂ overflow) are included. This is because, when using detailed binary evolution calculation, we are able to resolve the onset of the dynamical instability which is not instantaneous and is most often preceded by a short but intense phase of mass transfer. In fact, the systems that reach the highest luminosities correspond to the unstable systems on the higher donor-mass end in both panels of Fig. 6. Therefore, binaries that would suffer dynamical instabilities and coalesce (likely producing Thorne-Żytkow objects hypothesised by Thorne & Zytkow 1977) could also be observed as NS ULXs earlier in their evolution. However, the lifetimes of these binaries as ULXs is very short, and so this would act against their detectability.

Fig. 6. Left: time-averaged isotropic-equivalent accretion luminosities ($⟨L_{\mathrm{acc}}^{\mathrm{iso}}⟩$) in LMXBs/IMXBs with a 1.3 M⊙ NS accretor that reached an instantaneous luminosity above 10 LEdd. These luminosities were calculated based on Eqs. (9) and (10) for each binary and at each time-step, and averaged over the entire RLO lifetime. Grey colour denotes the sequences that never achieved an instantaneous luminosity above 10 LEdd. The systems that reach the highest accretion luminosities correspond to unstable systems, on the higher donor mass ends in both panels, before the onset of dynamical instability. The stable systems (enclosed by the solid black boundary) have a relatively low accretion luminosity on average. The magenta star corresponds to the binary shown in Fig. 5. Right: same as left but for LMXBs/IMXBs with 2.0 M⊙ NSs. The three red crosses correspond to systems that terminated due to numerical issues.

We define three X-ray luminosity ranges > 10 L_Edd, > 100 L_Edd, and > 1000 L_Edd, and calculate how long each system spends in each luminosity range. We refer to this time duration as ULX lifetime. The defined luminosity ranges are based on the observed pulsating ULX luminosities, which are in the range of 10–1000 L_Edd (Sect. 2). We compare the isotropic-equivalent accretion luminosities with these luminosity ranges. Figure 7 presents these results for LMXBs/IMXBs with 1.3 M_⊙ and 2.0 M_⊙ accretors. In systems with 1.3 M_⊙ NSs, the longest ULX lifetime is 1.6 × 10⁶ years, corresponding to $M_{\text{donor}}^i=2.3M_{\odot }$ and $P_{\text{orb}}^i=2.16$ days, for observed luminosities of > 10 L_Edd. In systems with 2.0 M_⊙ NSs, the longest ULX lifetime is 1.1 × 10⁶ years, corresponding to $M_{\text{donor}}^i=2.77M_{\odot }$ and $P_{\text{orb}}^i=2.16$ days, again for observed luminosities of > 10 L_Edd. The upper limits to ULX lifetimes are comparable to the ULX age estimate of ∼1 Myr for NGC 1313 X-2 by Pakull & Grisé (2008) (Sect. 2.6). Looking at similar initial donor masses and initial orbital periods in both sets of LMXBs/IMXBs (going from top to bottom row in Fig. 7), higher accretor mass corresponds to a much higher ULX lifetime as the stability area increases. As an example, for $M_{\text{donor}}^i=5.0M_{\odot }$ and $P_{\text{orb}}^i=1.0$ days, the ULX time increases from 3.0 × 10⁴ years to 3.0 × 10⁵ years (going from lower to higher NS accretor mass).

Fig. 7. Top row: ULX lifetime for LMXB/IMXB systems with a NS accretor mass of 1.3 M⊙. Left, middle, and right panels: time that systems spent with $L_{\mathrm{acc}}^{\mathrm{iso}}$ above 10, 100, and 1000 LEdd, respectively. White arrows (middle panels) enclose the potential properties of M82 X-2 at the onset of RLO (see Sects. 2.1 and 5.1). Bottom row: same as the top row but for LMXB/IMXB systems with a NS accretor mass of 1.3 M⊙.

In Fig. 7 (top row), from the leftmost panel to the rightmost (from > 10 L_Edd to > 1000 L_Edd) the systems which initially have long ULX lifetimes ($M_{\text{donor}}^i~2.3M_{\odot }$ and $P_{\text{orb}}^i~2.0$ days) either no longer appear on the plot or have a smaller ULX lifetime. Their isotropic-equivalent luminosities barely reach the higher cutoff values. Similarly, in the bottom row of Fig. 7, most LMXBs/IMXBs with long ULX lifetimes ($M_{\text{donor}}^i~2.77M_{\odot }$ and $P_{\text{orb}}^i~2.0$ days) for > 10 L_Edd, do not reach luminosities > 1000 L_Edd. This implies that binaries that are in the ULX phase for the longest time do not always achieve the highest luminosities. This effect is not seen in the unstable systems, which maintain their short ULX lifetime across the different ULX criteria. There is a part of the stable parameter space where the ULX lifetime decreases significantly, the ULX lifetime going from about 10⁵ years for systems outside this region to 10¹–10² years (the boundary corresponding to $M_{\mathrm{donor}}^i<3.0M_{\mathit{odot}}$ and $P_{\mathrm{orb}}^i<2.0$). Their accretion luminosity is sub-Eddington for almost the entire mass-transfer phase. After a mass-transfer episode, these systems are left with a He core and thin H envelope which expands briefly during H-shell burning, which is enough to increase their luminosity to exceed the Eddington limit albeit for a very short period of time.

Most NS ULXs observed so far have been in the higher luminosity range (100–1000 L_Edd) potentially due to selection effects as more luminous ULXs have a higher chance of being observed. On the other hand, the short lifespans of NS ULXs compared to the age of the universe (∼1.4 × 10¹⁰ yr) implies that it is unlikely to observe these systems in large numbers, which is consistent with the existing small NS ULX sample. The extra physical argument that the emission is highly beamed, especially for very luminous ULXs, further decreases the chances of detecting a large number of NS ULXs.

As alluded to above, the probability of observing any of these LMXB/IMXB sources as ULXs depends on the beaming factor, the ULX lifetime, and the probability of the particular system being formed. We evaluate the first two factors in the following equation describing the likelihood of a ULX observation, which is the beaming factor integrated over the three defined ULX lifetimes:

$$\begin{array}{c}\hfill \mathcal{L}={\displaystyle {\int }_{\mathrm{ULX}}b(t)\mathrm{d}t.}\end{array}$$(20)

Calculating the likelihood using this equation for both the grids and for all three ULX lifetime criteria, we obtain a measure of the relative chance of observing each system in Fig. 8 for the 1.3 M_⊙ NS (top row) and the 2.0 M_⊙ NS (bottom row) grids. All values below 10⁻⁶ are shown in black as they correspond to insignificant likelihood. The highest likelihood is for the stable systems with lower luminosities (10 L_Edd), where the mass-transfer rate never increases to extreme values and thus collimation is very small to almost negligible. Figure 8 shows that the chance of a system being observed as a ULX decreases at higher luminosities. This is due to the fact that for higher luminosities the mass-accretion rates (and the mass-transfer rates) are extremely high which causes the emission from the system to be highly beamed. Also, the probability depends on the structure of the donor at the onset of RLO. Case A systems have lower mass-transfer rates and more isotropic emission and thus slightly higher probabilities of being observed than case B. As mentioned before, there is a part in the stable parameter space where the ULX lifetime drops by many orders of magnitude (Fig. 7). Looking at the same systems in Fig. 8, the likelihood to observe them is negligible. These systems are super-Eddington for a very brief moment in time which is to the detriment of their detectability. Looking at Figs. 7 and 8, it is evident that even though LMXBs are included in the initial parameter space, IMXBs (with donor masses ≳2.0 M_⊙) are better candidates for NS ULXs.

Fig. 8. Top row: relative likelihood (ℒ) to observe a system as a ULX, as described by Eq. (20), for LMXBs/IMXBs with a 1.3 M⊙ NS. The panels are arranged as in Fig. 7. Going from the leftmost to the rightmost panel (> 10 to > 1000 LEdd), ℒ is higher for stable binaries with lower luminosities as their emission is not as strongly beamed. There is in addition some effect by the transition from case A to case B RLO. All values below 10−6 are shown in black as they correspond to insignificant likelihood. Bottom row: same as the top row but for LMXBs/IMXBs with a 2.0 M⊙ NS.

One caveat is that we assume that all systems have an equal probability of formation whereas in reality many systems might be formed at a higher rate than others. To account for these effects we would need to do population synthesis studies which would include exploring the formation probability of close NS + main sequence binaries and the distribution of their binary parameters. However, this is outside the scope of this study. We intend to study the formation rate of these systems in a future work.

Furthermore, one should always be careful when directly comparing our results to the observed population as our models give predictions about the whole population of ULXs with NS accretors, and not only the pulsating ones. The conditions for a NS to produce coherent pulses while also reaching super-Eddington luminosities are still unclear. If for example a strong magnetic field is required then perhaps the pulses are only observable at the beginning of the mass-transfer phase before any significant amount of material is accreted onto the NS, which could bury the magnetic field.

In cases where the mass-transfer sequence is stable throughout the binary evolution, we expect a NS–WD binary to be formed. Figure 9 shows the type of white dwarf formed at the end and the final white dwarf masses for the 1.3 M_⊙ NS (left panel) and the 2.0 M_⊙ NS (right panel) grids. The superscript f stands for final values. Overall the white dwarfs resulting from the 1.3 M_⊙ NS grid are in the mass range of 0.23–0.71 M_⊙. For systems with 2.0 M_⊙ NS accretors, the white dwarfs have an overall mass range of 0.25–0.95 M_⊙. We used final total mass fractions of carbon to distinguish between different WD types. We define WDs with > 95% carbon mass fraction as CO white dwarfs, 0.01–95% as hybrid white dwarfs, and < 0.01% as He white dwarfs. The initially higher donor masses and longer orbital periods result in a degenerate CO core with negligible helium on the surface. For some donor masses and orbital periods, the donor forms a degenerate CO core with a relatively large helium-rich envelope leading to hybrid WDs. In systems with low donor masses and short orbital periods, the donors end up as helium WD systems. The reason for the final fate of the latter class of systems is that due to deep envelope stripping early in the evolution of the donor, combined in some cases with the low initial mass of the donor, the final stripped helium cores are not massive enough to ignite helium.

Fig. 9. Left: resulting distribution of final WD masses for LMXBs/IMXBs with 1.3 M⊙ NSs. The He WDs are the least massive, followed by hybrid WDs, and the CO WDs are the most massive. The outermost black boundary encloses the stable systems and the inner boundaries differentiate between the different WD types. Right: same as left but for LMXBs/IMXBs with 2.0 M⊙ NSs.

Figure 10 shows the mass accreted by the NSs, for systems with 1.3 M_⊙ (left) and 2.0 M_⊙ NS (right) accretors. We should reiterate here our assumption of non-conservative mass transfer with β = max[0.7, 1 − Ṁ_Edd/Ṁ_donor] and that the accretion by the NS is Eddington limited. The 1.3 M_⊙ NSs end up accreting 0.004–0.415 M_⊙ of mass by the end of the mass-transfer phase, the systems with values around $M_{\text{donor}}^i=2.0M_{\odot }$ and $P_{\text{orb}}^i=0.87$ day accreting the most. The amount of mass accreted is not enough to collapse the NS. The 2.0 M_⊙ NSs accrete mass in the range of 0.002–0.585 M_⊙, with the systems with values around $M_{\text{donor}}^i=3.0M_{\odot }$ and $P_{\text{orb}}^i=3.1$ days accreting the highest amount. We compare our results to the maximum NS mass known. The most massive precisely measured NS is PSR J0348+0432 with a mass of 2.01 ± 0.04 M_⊙ (Antoniadis et al. 2013). However, recently, a candidate with a higher NS mass of $2.{17}_{-0.10}^{+0.11}{M}_{\odot }$ was announced (PSR J0740+6620; Cromartie et al. 2020). Although the error bar of the latter source is relatively large, we take 2.17 M_⊙ as our assumed upper limit. This value is also supported by constraints on GW170817 based on combined gravitational wave and electromagnetic observations (Margalit & Metzger 2017).

Fig. 10. Left: mass accreted by the NS ($M_{\mathrm{acc}}^f-M_{\mathrm{acc}}^I$) for LMXBs/IMXBs with a 1.3 M⊙ NS and β = max[0.7, 1 − ṀEdd/Ṁdonor]. Binaries producing the lowest WD masses show the highest accretion on the NS. Right: same as left but for LMXBs/IMXBs with a 2.0 M⊙ NS. The systems below the white boundary are those where the accretor has accreted enough mass to exceed the maximum mass limit for a NS, here assumed to be 2.17 M⊙.

Looking at the final NS masses in Fig. 10 (right panel), we see a certain population enclosed by a white boundary. In these systems the accretor has accreted enough material to overcome the neutron degeneracy pressure that is supporting the star against gravitational collapse (assuming an upper NS mass limit of 2.17 M_⊙). Therefore, according to the maximum NS mass limit assumed, these NSs will collapse to form BHs.

Net accretion on the NS can be broadly described by,

$$\begin{array}{c}\hfill \mathrm{\Delta }{M}_{\mathrm{acc}}=⟨{\dot{M}}_{\mathrm{acc}}⟩\times \mathrm{\Delta }{t}_{\mathrm{X}},\end{array}$$(21)

where Δt_X is the lifetime as an X-ray binary (or the overall mass-transfer phase) and ⟨Ṁ_acc⟩ is the average accretion rate onto the NS. A high amount of accreted mass is achieved with a combination of both a large accretion rate and a long time duration over which it occurs. For LMXBs/IMXBs that reach up to Eddington mass-transfer rates for only a small part of the mass-transfer phase, to zeroth order, ΔM_acc ≃ 0.3 × ΔM_donor. In contrast, for the brightest ULXs, where most of the mass-transfer happens at a highly super-Eddington rate, ⟨Ṁ_acc⟩ ≃ Ṁ_Edd and ΔM_acc ≪ ΔM_donor. Even if the accretion onto the NS was allowed to reach a few times the Eddington limit, the net amount of accreted material by the NS in the brightest ULXs would still be small.

The final orbits in most systems that underwent stable mass transfer are wide because during the mass-transfer phase, after the orbit shrinks initially, the donor evolves to become less massive than the accretor causing the system to widen significantly. These widened orbits can be seen in Fig. 11 which shows that the final orbital periods can be as large as 5–7 times the initial orbital periods. Some of these binaries may be observable by Gaia. Andrews et al. (2019) postulate that Gaia can detect and measure the properties of hidden wide binaries with WD and NS components using astrometric observations. The comparison of NS–WD observations to final orbital parameters from our calculations might help in understanding which binaries have undergone a super-Eddington mass-transfer phase in the past and could have been observed as ULXs, as well as constrain binary evolution physics such as the accretion efficiency.

Fig. 11. Left: final orbital periods $(P_{\mathrm{orb}}^f)$ for LMXBs/IMXBs with a 1.3 M⊙ NS. The longest final orbital periods result from initial periods ≳10.0 days and donors with masses less than about 2.75 M⊙. The shortest final orbital periods are for systems that formed the heaviest WDs, originating from donors of ≳3.5 M⊙. Right: same as left but for LMXBs/IMXBs with 2.0 M⊙ NSs.

5. Discussion

5.1. Comparison to observations of M82 X-2

Since M82 X-2 has fairly well constrained parameters (see Sect. 2.1 and Table 1), we can compare it to our results and see if our results help to explain the observations. We look at the middle panels in both rows in Figs. 7 and 8. Since the currently observed donor is more massive than the accretor, the orbit will shrink as mass is lost following the orbital angular momentum balance equation as described in Eq. (22) for a non-conservative mass-transfer phase. Therefore, we use the observed parameters as lower limits of the initial binary configuration.

We only consider the luminosity range of 100 L_Edd because it is comparable to the observed luminosity of about 1.8 × 10⁴⁰ erg s⁻¹. Bachetti et al. (2014) assumed an accretor of 1.4 M_⊙ and estimated the donor to be ≳5.2 M_⊙. For a 1.3 M_⊙ accretor mass, the donor mass is ≳5.1 M_⊙, using the same binary mass function. The initial parameter estimates are enclosed by the white arrows in Figs. 7 and 8 (middle panels in both rows). The potential initial systems have ULX lifetimes as long as 0.7 × 10⁴ years and a relative peak likelihood of 0.7 × 10⁻⁴. However, if we assume a higher accretor mass of 2.0 M_⊙, the donor mass is ≳5.83 M_⊙. The longest ULX lifetime in this case for the potential initial systems is about 1.1 × 10⁵ years with a peak likelihood of 0.7 × 10⁻². Taking these numbers at face value, it is clear that an initially heavy NS is preferred in order to explain the currently observed properties of M82 X-2. This cannot be excluded as the 1.4 M_⊙ NS used by Bachetti et al. (2014) was an assumption. We should note, however, that for both NS masses, the peak of the relative likelihood does not lie very close to the observed limits we have for the current properties of M82 X-2. This is not necessarily problematic, as in order to calculate the actual probability density distribution of what the properties of NS ULXs ought to be, based on our model, we need to multiply the relative likelihood shown in Fig. 8 with the “prior” probability of forming a NS binary with these initial conditions. This convolution might significantly shift the peak of the resulting probability distribution. We leave this calculation for future work. Fragos et al. (2015), who followed such an approach, estimated the most probable initial donor mass of any NS ULX to be 3.0–8.0 M_⊙, and the initial orbital period to be 1.0–3.0 days. This parameter space lies within our results.

5.2. Comparison with high-mass X-ray binary ULXs

With our work we aim to explore the possibility of an LMXB or IMXB origin for pulsating ULXs. However, there are at least three known NS ULXs which are HMXBs, namely NGC 7793 P13, M51 ULX-7, and more recently NGC 300 ULX1.

NGC 7793 P13 has a luminosity of about 100 L_Edd. With the presence of an observed high-mass donor (> 18 M_⊙, Sect. 2.2 and Table 1), this system cannot be explained by our LMXB/IMXB models. In our NS ULXs, extreme mass ratios produce unstable binaries which reach the ULX luminosity observed by the source for a timescale of ∼10²–10⁴ years before initiating L₂ overflow. However, this system is a HMXB with a donor of mass 18.0 M_⊙ <  M_donor <  23.0 M_⊙ and is outside the range of our simulations.

M51 ULX-7 is another HMXB in a ULX phase and containing an accreting NS. Rodríguez Castillo et al. (2020) suggest that the binary consists of an OB giant star (M_donor ≲ 8.0–13.0 M_⊙) and a NS with a magnetic field of 10¹²–10¹³ G. The most recent NS ULX observed to have a super-giant donor is NGC 300 ULX1, with a donor star mass of M_donor ≳ 8.0 M_⊙ (Heida et al. 2019). These two systems are also not satisfactorily explained by our calculations, since we find such systems have a ULX lifetime of ≲10⁴ years (for donors of 8.0 M_⊙).

Quast et al. (2019) explored the stability of super-giant X-ray binaries with a very high donor mass compared to accretor mass (donors up to 20 times more massive than the accretors), and found that the ULX phase could be long lasting (∼0.4 × 10⁶ years) for such an extreme mass ratio if the super-giant star has a H/He gradient in the layers beneath its surface. These latter authors demonstrated that such systems can evolve on a nuclear timescale with a BH or NS accretor, even for binaries where the donor mass is up to 20 times the accretor mass. In their binary evolution models, the donor stars rapidly decrease their thermal equilibrium radius and can therefore cope with the inevitably strong orbital contraction imposed by such a high mass ratio. These binaries could be post-CE systems, where the super-giant donor star has lost most of its hydrogen-rich envelope and the remaining one is enriched in helium. Recent 1D hydrodynamical simulations of a CE phase between a super-giant donor and a NS predict the formation of such binary configurations (Fragos et al. 2019).

5.3. Effect of spin-orbit coupling on orbital evolution

As we already described in Sect. 3.3, our binary evolution calculations take into account the effect of spin-orbit coupling by modelling the tidal interactions, internal rotation, and mass transfer and/or loss. Here we want to quantify the effect of spin–orbit coupling on the orbital evolution and hence on the stability of mass transfer. We perform tests to calculate the orbital evolution in some limiting cases and compare our numerical calculations with others available in the literature, as well as analytical solutions.

For simplicity, we assume a fully non-conservative mass transfer (β = 1.0) in the following calculations, i.e. no mass is accreted by the NS. Also, we assume that the binary remains always in synchronous rotation with its components so that the spin period of the donor is the same as the orbital period. This assumption corresponds to an infinitely efficient spin–orbit coupling. Under these assumptions, we obtain the analytical solution for the orbital evolution as follows (see Appendix A for the detailed derivation):

$$\begin{array}{cc}\hfill \frac{\dot{a}}{a}& =[\frac{M_{\mathrm{donor}}{M}_{\mathrm{acc}}{a}^2}{M_{\mathrm{donor}}{M}_{\mathrm{acc}}{a}^2-3I_{\mathrm{donor}}M}]\times [{\dot{M}}_{\mathrm{donor}}(\frac{2M_{\mathrm{donor}}}{M_{\mathrm{acc}}M}-\hfill \\ \hfill & \frac{2}{M_{\mathrm{donor}}}+\frac{1}{M})-(\frac{2{\dot{I}}_{\mathrm{donor}}M+I_{\mathrm{donor}}{\dot{M}}_{\mathrm{donor}}}{M_{\mathrm{donor}}{M}_{\mathrm{acc}}{a}^2})],\hfill \end{array}$$(22)

where I_donor is the moment of inertia of the donor. When the system loses mass from the vicinity of the accretor, according to Eq. (22) the orbit shrinks. This decrease occurs as long as M_donor ≳ 1.28 M_acc (Tauris et al. 2001). The orbit expands if mass is being lost from the approximate vicinity of the more massive star.

We run models using MESA, including angular momentum losses from spin–orbit coupling and mass lost from the system, while ignoring gravitational wave radiation and magnetic braking. The results of this test are shown in Fig. 12 (blue and orange curves), and they are compared to the analytical solution of Eq. (22) (solid and dashed black curves), where I_donor and İ_donor are taken from numerical stellar structure calculations. For two extreme cases of spin–orbit coupling, we first consider the case where the donor star is an extended body (I_donor ≠ 0) and there is efficient coupling of donor spin and orbit and thus an efficient transfer of angular momentum between spin and orbit (orange curve in the figure). This keeps the spin of the donor synchronised with the orbit throughout the binary evolution and the binary is “tidally locked”. The second case is where both binary components are approximated as point masses and there is no exchange of angular momentum between the donor (blue curve in the figure) and the orbit (I_donor = 0). The two cases described above are extremes of spin–orbit coupling; any realistic solution would be in-between them. The MESA calculations are consistent with an analytic solution and the coupling does have a slight effect on the orbital evolution.

Fig. 12. Binary separation evolutionary sequences for a system with initially Mdonor = 4 M⊙, Macc = 1.3 M⊙ and 8 days orbital period (β = 1.0). Different cases showing spin–orbit and no spin–orbit coupling are shown in comparison with the analytical equation calculations (described by Eq. (22)).

We also compare two MESA runs with the results of Tauris et al. (2001, their Fig. 2). In this comparison we assume that β = max[0, 1 − Ṁ_Edd/Ṁ_donor], and therefore that mass is transferred conservatively until the Eddington limit is reached. Furthermore, gravitational wave radiation and magnetic braking are included. This comparison is shown in Fig. 13. We again consider two cases, one with a point mass approximation (blue and dashed black curves) and the other one with efficient spin-orbit coupling (orange and solid black curves). While the evolutionary sequences for the point mass approximation (blue curve) match those from Tauris et al. (2001) (dashed black curve), the sequences corresponding to an efficient spin–orbit coupling show a discrepancy. The effect of angular momentum transport between the donor star and the orbit seems to be much larger in Tauris et al. (2001) (solid black curve) compared to our calculated MESA tracks and analytic solutions (orange curve). In an attempt to understand this discrepancy, we also considered the case of highly inefficient synchronisation during RLO. In this scenario, during the high mass-loss phase, the donor loses mass (and angular momentum) which causes its spin to decrease. Since the synchronisation is not efficient, the orbit does not regulate the spin of the star as quickly as angular momentum is being lost from the star. Pratt & Strittmatter (1976) explore this case as a means to have a wider orbit and we also present this run in Fig. 13 as the solid red line. This model has highly efficient spin–orbit coupling until the point the donor almost overflows its Roche lobe and begins mass transfer after which the spin–orbit coupling is highly inefficient. However, this reasoning does not seem to explain the described discrepancy either. Unfortunately, it is impossible to further investigate this apparent discrepancy, as the original calculations and code setup from Tauris et al. (2001) were not available.

Fig. 13. Binary separation evolutionary sequences for a system with initially Mdonor = 4 M⊙, Macc = 1.3 M⊙ and 8 days orbital period (β = max {0, 1 − ṀEdd/Ṁdonor}). Different cases showing efficient spin–orbit and point mass approximation are shown in comparison with similar tracks from Tauris et al. (2001) (Fig. 2 in their paper).

In conclusion, due to the fact that the widening effect of spin–orbit coupling is weak in our calculations, we infer that spin–orbit coupling can affect the mass-transfer stability but is a second-order effect.

6. Conclusions

There is much intrigue surrounding the nature of ULXs because they challenge our understanding of accretion onto compact objects at high rates and of the evolution of binaries containing compact objects. After the discovery of the first pulsating ULX (M82 X-2), NS X-ray binary systems were known to constitute a (possibly major) sub-population of ULXs.

NS IMXBs, with donor masses of 2.0 − 10.0 M_⊙, have been proposed in the past as possible ULXs (Tauris et al. 2000; Fragos et al. 2015). Their extreme mass ratio leads to intense mass transfer, which can become highly super-Eddington. Due to changes in the structure of the accretion disc as the mass-transfer rate exceeds the Eddington limit, face-on observers can infer super-Eddington isotropic-equivalent accretion luminosities, even if the accretion onto the NS is actually limited to the Eddington rate. This is because the increased radiation pressure causes the inner parts of the disc to puff up, creating a funnel, which in turn can cause beaming of the outgoing radiation. These effects can increase the observed isotropic-equivalent luminosities of some accreting NSs to a few orders of magnitude above the Eddington limit. However, the higher the accretion luminosity the lower the chances of actually observing the source since the emission is highly beamed and depends on the angle of viewing.

In this work, we attempt to explain super-Eddington luminosities observed in NS ULXs by calculating extensive grids of binary evolution models of LMXBs/IMXBs with donor stars in the mass range of 0.92–8.0 M_⊙ using the MESA code. We specifically investigate the stability and rate of mass transfer, taking into account rotation and tidal effects, as well as the observable properties and the final outcomes of these binaries. The results we obtain can be summarised as follows:

• It has been demonstrated previously that the stability of mass transfer in a binary system primarily depends on the structure of the envelope of the donor star at the onset of RLO. From our numerical modelling, including stellar rotation and tides, we find that the most extreme mass ratios (q = M_acc/M_donor) that are able to evolve through long-term stable mass transfer are 0.5 for donors with convective envelopes and 0.28 for donors with radiative envelopes.

• For more extreme mass ratios, mass transfer via the first Lagrangian point (L₁) cannot keep the donor star contained within its Roche lobe, as the Roche-lobe radius is shrinking faster than the radius of the donor star. In these cases, the donor star expands well outside its Roche lobe, and its volume can even exceed the volume enclosed by the equipotential surface that passes from the second Lagrangian point (L₂). At this stage, mass loss from the L₂ point will occur, which causes significant loss of orbital angular momentum, thereby accelerating the shrinkage of the orbit and leading to the onset of a CE event. We have calculated L₂ lobe volumes and derived fitting formulae (accurate to errors less than 1% for all values of q) to calculate the volume-equivalent radii of donor stars at the onset of mass loss via L₂. In our calculations, we self-consistently identify the onset of the L₂ overflow and consider this to be the beginning of an irreversible dynamical instability.

• We demonstrate that LMXBs/IMXBs can produce NS-ULXs with typical time-averaged isotropic-equivalent X-ray luminosities between 10³⁹ and 10⁴¹ erg s⁻¹ on a timescale of up to ∼1.0 Myr for the lower luminosities. We present plots of their likelihood of detection, taking into account beaming effects. Based on these results, IMXBs are better candidates for NS ULXs than LMXBs.

• We find that the beamed super-Eddington accretion disc model provides the very high isotropic-equivalent X-ray luminosities that have been observed for NS ULXs. The model also explains the low detectability of these sources as highly beamed emission and short ULX lifetimes are common properties of the sources.

• LMXBs/IMXBs cannot explain the origin of all NS ULXs given that some of their donor stars are massive (i.e. therefore HMXBs). However, we confirm that IMXBs, specifically, are strong candidates for a large fraction of NS ULXs. For the most well-studied NS ULX, M82 X-2, our LMXB/IMXB mass-transfer sequences with 2.0 M_⊙ NSs better explain the observations compared to systems with 1.3 M_⊙ NSs.

• We obtain three populations of WDs at the end of the stable mass-transfer phase: He WDs, CO WDs, and hybrid WDs which have a CO core and a significant amount of He in the envelope. Comparing observations of NS–WD binaries could help to infer whether or not a binary underwent a ULX phase in the past.

• The tidal coupling between the spin and the orbit has a second-order effect on the stability of the binary. Our models contradict some previous results in the literature.

Acknowledgments

The authors thank the anonymous referee for their constructive comments that helped improve the manuscript. The authors are also grateful to Pablo Marchant, Vasilopoulos Georgios, and Andrew King for their useful discussions. This work was supported by the Swiss National Science Foundation Professorship grant (project number PP00P2 176868). This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie RISE action, grant agreement no. 691164 (ASTROSTAT). The authors acknowledge the International Space Science Institute (ISSI) for supporting and funding the international teams program “Ultraluminous X-ray Sources: From the Local Group to the Very First Galaxies”. D.M. thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; her participation in the program has benefited this work. T.M.T. acknowledges an AIAS-COFUND Senior Fellowship funded by the European Union’s Horizon 2020 Research and Innovation Programme (grant no. 754513) and the Aarhus University Research Foundation. E.Z. acknowledges support from the Federal Commission for Scholarships for Foreign Students for the Swiss Government Excellence Scholarship (ESKAS No. 2019.0091) for the academic year 2019-2020.

References

Appendix A: Analytical solution of orbital evolution

The orbital angular momentum (J_orb) of a binary with component masses as M_donor and M_acc can be described as,

$$\begin{array}{c}\hfill J_{\mathrm{orb}}=M_{\mathrm{donor}}{M}_{\mathrm{acc}}\sqrt{\frac{Ga(1-e^2)}{M}},\end{array}$$(A.1)

where M = M_donor + M_acc, a is the binary separation, and e is the eccentricity of the orbit. For a circular orbit (e = 0), a logarithmic differentiation of Eq. (A.1) gives the rate of change of the orbital separation (ȧ),

$$\begin{array}{c}\hfill \frac{\dot{a}}{a}=2\frac{{\dot{J}}_{\mathrm{orb}}}{J_{\mathrm{orb}}}-2\frac{{\dot{M}}_{\mathrm{donor}}}{M_{\mathrm{donor}}}-2\frac{{\dot{M}}_{\mathrm{acc}}}{M_{\mathrm{acc}}}+\frac{{\dot{M}}_{\mathrm{donor}}+{\dot{M}}_{\mathrm{acc}}}{M}.\end{array}$$(A.2)

We consider a binary system which is in synchronous rotation with its components as tides are efficiently transferring angular momentum from donor to the orbit and vice versa. That is, it is tidally locked. Therefore, the orbital angular velocity (ω_orb) is equal to the spin angular velocity (ω_donor), which we take as ω_orb = ω_donor = ω.

The mass lost from the binary carries away angular momentum from the orbit. For a such a system containing a mass losing donor star (M_donor) and an accreting neutron star (or any other compact object with mass, M_acc), the changes in the orbital angular momentum are described by

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{orb}}=\omega a_2^2{\dot{M}}_{\mathrm{donor}}+{\dot{J}}_{\mathrm{ls}},\end{array}$$(A.3)

where a₂ = aM_donor/M is the distance from where the mass transferred is lost from the system, that is from around M_acc. The term ${\dot{J}}_{\mathrm{ls}}$ stands for the contribution to the rate of change of the total angular momentum (${\dot{J}}_{\mathrm{orb}}$) from the donor spin (${\dot{J}}_{\mathrm{spin}}$) and vice versa, due to the coupling between the spin and the orbit. Changes in the spin angular momentum of the donor are defined by

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{spin}}=-{\dot{J}}_{\mathrm{ls}}.\end{array}$$(A.4)

The donor is undergoing solid-body rotation and so its angular momentum can be described as J_spin = I_donorω, and therefore ${\dot{J}}_{\mathrm{spin}}$ can also be expressed in general as

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{spin}}={\dot{I}}_{\mathrm{donor}}\omega +I_{\mathrm{donor}}\dot{\omega },\end{array}$$(A.5)

where I_donor and ω are the moment of inertia and angular velocity of the donor. Using Eqs. (A.3) and (A.5), the rate of change in total orbital angular momentum comes out to be

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{orb}}=\omega a_2^2{\dot{M}}_{\mathrm{donor}}-{\dot{I}}_{\mathrm{donor}}\omega -I_{\mathrm{donor}}\dot{\omega }.\end{array}$$(A.6)

Substituting Eq. (A.6) into Eq. (A.2) and using the assumption of fully non-conservative mass transfer, i.e. Ṁ_acc = 0,

$$\begin{array}{c}\hfill \frac{\dot{a}}{a}=2\frac{\omega a_2^2{\dot{M}}_{\mathrm{donor}}-{\dot{I}}_{\mathrm{donor}}\omega -I_{\mathrm{donor}}\dot{\omega }}{J_{\mathrm{orb}}}-2\frac{{\dot{M}}_{\mathrm{donor}}}{M_{\mathrm{donor}}}+\frac{{\dot{M}}_{\mathrm{donor}}}{M}.\end{array}$$(A.7)

Using the definitions of $J_{\mathrm{orb}}=M_{\mathrm{donor}}{M}_{\mathrm{acc}}\sqrt{Ga/M}$, $\omega =\sqrt{GM/a^3}$, the logarithmic differentiation of ω, $\dot{\omega }=\omega ({\dot{M}}_{\mathrm{donor}}/M-3\dot{a}/a)/2$, and simplifying further we get the analytical solution for the orbital separation,

$$\begin{array}{cc}\hfill \frac{\dot{a}}{a}& =[\frac{M_{\mathrm{donor}}{M}_{\mathrm{acc}}{a}^2}{M_{\mathrm{donor}}{M}_{\mathrm{acc}}{a}^2-3I_{\mathrm{donor}}M}]\times [{\dot{M}}_{\mathrm{donor}}(\frac{2M_{\mathrm{donor}}}{M_{\mathrm{acc}}M}-\frac{2}{M_{\mathrm{donor}}}+\frac{1}{M})\hfill \\ \hfill & -(\frac{2{\dot{I}}_{\mathrm{donor}}M+I_{\mathrm{donor}}{\dot{M}}_{\mathrm{donor}}}{M_{\mathrm{donor}}{M}_{\mathrm{acc}}{a}^2})].\hfill \end{array}$$(A.8)

This equation describes the evolution of orbital separation in a binary where there is efficient exchange of angular momentum between the donor and the orbit and the donor is undergoing solid-body rotation.

Point mass approximation

In case of no exchange of angular momentum between spin and orbit, the binary components are treated as point masses, therefore I_donor = 0. Using this assumption, the orbital evolution equation simplifies as follows,

$$\begin{array}{c}\hfill \frac{\dot{a}}{a}={\dot{M}}_{\mathrm{donor}}(\frac{2M_{\mathrm{donor}}}{M_{\mathrm{acc}}M}-\frac{2}{M_{\mathrm{donor}}}+\frac{1}{M}).\end{array}$$(A.9)

If we look at Eq. (16.18) in (Tauris & van den Heuvel 2006) which is,

$$\begin{array}{c}\hfill \frac{{\dot{J}}_{\mathrm{orb}}}{J_{\mathrm{orb}}}=[\frac{\alpha +\beta {(M_{\mathrm{donor}}/M_{\mathrm{acc}})}^2+\delta \gamma (1+{(M_{\mathrm{donor}}/M_{\mathrm{acc}})}^2)}{1+M_{\mathrm{donor}}/M_{\mathrm{acc}}}]\frac{{\dot{M}}_{\mathrm{donor}}}{M_{\mathrm{donor}}},\end{array}$$(A.10)

where α is the fractional mass lost directly from the donor, β is the fractional mass lost from the vicinity of the accretor, δ is the fractional mass lost from a circumbinary toroid, and the radius of the circumbinary toroid is defined by γ²a.

Using the assumptions Ṁ_acc = 0 ⇒ β = 1, and α = δ = 0, we substitute Eq. (A.10) in Eq. (A.2) and simplify to get the following,

$$\begin{array}{c}\hfill \frac{\dot{a}}{a}={\dot{M}}_{\mathrm{donor}}(\frac{2M_{\mathrm{donor}}}{M_{\mathrm{acc}}M}-\frac{2}{M_{\mathrm{donor}}}+\frac{1}{M}).\end{array}$$(A.11)

This equation is consistent with our calculations where we obtain Eq. (A.9).

All Tables

All Figures

	Fig. 1. Equipotential lines of the Roche potential for a binary consisting of stars with mass ratio q = 0.26 and binary separation a. The equipotential lines passing through Lagrangian points L1, L2, and L3 are shown.
In the text

	Fig. 2. Equipotential lines of the Roche potential for a binary consisting of stars with mass ratio q = 1.85 and binary separation a to illustrate the swapping of the positions of L2 and L3 when Mdonor ≤ Macc (compared to Fig. 1).
In the text

Fig. 3. Our results form the volume integration of L2 equipotential surfaces with respect to different mass ratios q. All distance values are presented in units of RL1. Here, RL2 is the radius of a sphere with volume equal to that of the L2 equipotential surface (solid orange line). DL2 is the distance of the L2 point from Mdonor (solid blue line). The dashed black line shows the orbital separation in units of RL1. Systems where the donor’s radius exceeds any of these limiting radii are considered as undergoing dynamical instability.

In the text

Fig. 4. Allowed initial parameter space for LMXBs/IMXBs to undergo stable mass transfer with 1.3 M⊙ accretors (orange squares) and 2.0 M⊙ accretors (blue squares). Grey squares correspond to systems that encountered dynamical instability. The dark red squares towards the left edge of the figure correspond to systems that never initiated RLO. The lower dashed white line separates systems undergoing case A mass transfer from those undergoing case B mass transfer (same boundary in both grids). The middle dotted white line separates early case B (where the donor has yet to form a deep convective envelope) from late case B RLO (where donor has formed a deep convective envelope at onset of RLO). The upper dashed white curve encloses systems that undergo case C mass transfer (same boundary in both grids). Green stars correspond to those systems that undergo a second mass-transfer phase from a stripped Helium-giant star, i.e. case BB RLO. The three red crosses in case BB correspond to systems that terminated due to numerical issues.

In the text

Fig. 5. Estimated observed X-ray luminosity for a binary with $M_{\text{donor}}^i=2.77M_{\odot }$, $M_{\text{acc}}^i=1.30M_{\odot }$, and $P_{\text{orb}}^i=5.38$ days under different assumptions. The solid green curve is the accretion luminosity corresponding to the mass-accretion rate onto the NS, assuming a fixed radiative efficiency (η = 0.1). The solid orange curve is the accretion luminosity corresponding to a super-Eddington accretion disc following the Shakura & Sunyaev (1973) disc model given by Eq. (8). The solid blue curve is the isotropic-equivalent accretion luminosity corresponding to beamed super-Eddington emission following the King et al. (2001) geometric beaming model (limited at 1042 erg s−1) given by Eqs. (9) and (10). Four reference luminosity values, corresponding to LEdd (solid black line), 10 LEdd (dashed black line), 100 LEdd (dotted black line), and 1000 LEdd (dot-dashed black line), are also plotted. The initial properties of this binary are highlighted with a magenta star in Fig. 6.

In the text

Fig. 6. Left: time-averaged isotropic-equivalent accretion luminosities ($⟨L_{\mathrm{acc}}^{\mathrm{iso}}⟩$) in LMXBs/IMXBs with a 1.3 M⊙ NS accretor that reached an instantaneous luminosity above 10 LEdd. These luminosities were calculated based on Eqs. (9) and (10) for each binary and at each time-step, and averaged over the entire RLO lifetime. Grey colour denotes the sequences that never achieved an instantaneous luminosity above 10 LEdd. The systems that reach the highest accretion luminosities correspond to unstable systems, on the higher donor mass ends in both panels, before the onset of dynamical instability. The stable systems (enclosed by the solid black boundary) have a relatively low accretion luminosity on average. The magenta star corresponds to the binary shown in Fig. 5. Right: same as left but for LMXBs/IMXBs with 2.0 M⊙ NSs. The three red crosses correspond to systems that terminated due to numerical issues.

In the text

Fig. 7. Top row: ULX lifetime for LMXB/IMXB systems with a NS accretor mass of 1.3 M⊙. Left, middle, and right panels: time that systems spent with $L_{\mathrm{acc}}^{\mathrm{iso}}$ above 10, 100, and 1000 LEdd, respectively. White arrows (middle panels) enclose the potential properties of M82 X-2 at the onset of RLO (see Sects. 2.1 and 5.1). Bottom row: same as the top row but for LMXB/IMXB systems with a NS accretor mass of 1.3 M⊙.

In the text

Fig. 8. Top row: relative likelihood (ℒ) to observe a system as a ULX, as described by Eq. (20), for LMXBs/IMXBs with a 1.3 M⊙ NS. The panels are arranged as in Fig. 7. Going from the leftmost to the rightmost panel (> 10 to > 1000 LEdd), ℒ is higher for stable binaries with lower luminosities as their emission is not as strongly beamed. There is in addition some effect by the transition from case A to case B RLO. All values below 10−6 are shown in black as they correspond to insignificant likelihood. Bottom row: same as the top row but for LMXBs/IMXBs with a 2.0 M⊙ NS.

In the text

	Fig. 9. Left: resulting distribution of final WD masses for LMXBs/IMXBs with 1.3 M⊙ NSs. The He WDs are the least massive, followed by hybrid WDs, and the CO WDs are the most massive. The outermost black boundary encloses the stable systems and the inner boundaries differentiate between the different WD types. Right: same as left but for LMXBs/IMXBs with 2.0 M⊙ NSs.
In the text

Fig. 10. Left: mass accreted by the NS ($M_{\mathrm{acc}}^f-M_{\mathrm{acc}}^I$) for LMXBs/IMXBs with a 1.3 M⊙ NS and β = max[0.7, 1 − ṀEdd/Ṁdonor]. Binaries producing the lowest WD masses show the highest accretion on the NS. Right: same as left but for LMXBs/IMXBs with a 2.0 M⊙ NS. The systems below the white boundary are those where the accretor has accreted enough mass to exceed the maximum mass limit for a NS, here assumed to be 2.17 M⊙.

In the text

	Fig. 11. Left: final orbital periods $(P_{\mathrm{orb}}^f)$ for LMXBs/IMXBs with a 1.3 M⊙ NS. The longest final orbital periods result from initial periods ≳10.0 days and donors with masses less than about 2.75 M⊙. The shortest final orbital periods are for systems that formed the heaviest WDs, originating from donors of ≳3.5 M⊙. Right: same as left but for LMXBs/IMXBs with 2.0 M⊙ NSs.
In the text

	Fig. 12. Binary separation evolutionary sequences for a system with initially Mdonor = 4 M⊙, Macc = 1.3 M⊙ and 8 days orbital period (β = 1.0). Different cases showing spin–orbit and no spin–orbit coupling are shown in comparison with the analytical equation calculations (described by Eq. (22)).
In the text

	Fig. 13. Binary separation evolutionary sequences for a system with initially Mdonor = 4 M⊙, Macc = 1.3 M⊙ and 8 days orbital period (β = max {0, 1 − ṀEdd/Ṁdonor}). Different cases showing efficient spin–orbit and point mass approximation are shown in comparison with similar tracks from Tauris et al. (2001) (Fig. 2 in their paper).
In the text