Population synthesis of Be X-ray binaries: metallicity dependence of total X-ray outputs

ABSTRACT

X-ray binaries (XRBs) are thought to regulate cosmic thermal and ionization histories during the Epoch of Reionization and Cosmic Dawn (z ∼ 5–30). Theoretical predictions of the X-ray emission from XRBs are important for modelling such early cosmic evolution. Nevertheless, the contribution from Be-XRBs, powered by accretion of compact objects from decretion discs around rapidly rotating O/B stars, has not been investigated systematically. Be-XRBs are the largest class of high-mass XRBs (HMXBs) identified in local observations and are expected to play even more important roles in metal-poor environments at high redshifts. In light of this, we build a physically motivated model for Be-XRBs based on recent hydrodynamic simulations and observations of decretion discs. Our model is able to reproduce the observed population of Be-XRBs in the Small Magellanic Cloud with appropriate initial conditions and binary stellar evolution parameters. We derive the X-ray output from Be-XRBs as a function of metallicity in the (absolute) metallicity range Z ∈ [10⁻⁴, 0.03] with a large suite of binary population synthesis (BPS) simulations. The simulated Be-XRBs can explain a non-negligible fraction (⁠|$\gtrsim 30{{\ \rm per\ cent}}$|⁠) of the total X-ray output from HMXBs observed in nearby galaxies for Z ∼ 0.0003–0.02. The X-ray luminosity per unit star formation rate from Be-XRBs in our fiducial model increases by a factor of ∼8 from Z = 0.02 to Z = 0.0003, which is similar to the trend seen in observations of all types of HMXBs. We conclude that Be-XRBs are potentially important X-ray sources that deserve greater attention in BPS of XRBs.

1 INTRODUCTION

During the Epoch of Reionization and Cosmic Dawn (z ∼ 5–30), X-ray binaries (XRBs) are expected to be the dominant sources of X-rays that regulate the thermal and ionization evolution and small-scale structure of the interstellar medium (IGM; e.g. Fragos et al. 2013b; Fialkov, Barkana & Visbal 2014; Fialkov & Barkana 2014; Pacucci et al. 2014; Madau & Fragos 2017; Eide et al. 2018), as well as early star formation (e.g. Jeon et al. 2012; Artale, Tissera & Pellizza 2015; Hummel et al. 2015; Ricotti 2016; Park, Ricotti & Sugimura 2021a, b, 2023b).¹ They can leave unique signatures in the 21-cm signal from neutral hydrogen, which is one of the most promising probes of early structure/galaxy/star formation and cosmology (e.g. Fialkov et al. 2017; Bowman et al. 2018; Ewall-Wice et al. 2018; Ma et al. 2018; Madau 2018; Mirocha & Furlanetto 2019; Schauer, Liu & Bromm 2019; Chatterjee et al. 2020; Qin et al. 2020; Gessey-Jones et al. 2022; Kamran et al. 2022; Kaur et al. 2022; Kovlakas et al. 2022; Magg et al. 2022; Muñoz et al. 2022; Acharya, Cyr & Chluba 2023; Bevins et al. 2023; Hassan et al. 2023; Lewis et al. 2023; Ma et al. 2023; Mondal & Barkana 2023; Shao et al. 2023; Ventura et al. 2023; Yang, Li & Li 2023), and in the Lyman-α forest from long-lasting relics of reionization (Montero-Camacho, Zhang & Mao 2023). To fully unleash the power of the 21-cm probe and break potential degeneracy between astrophysics, dark matter physics, and cosmology (e.g. Barkana 2018; Liu, Schauer & Bromm 2019; Yang 2021; Yang, Li & Li 2023; Ghara, Mellema & Zaroubi 2022; Acharya, Cyr & Chluba 2023; Mondal, Barkana & Fialkov 2023; Shao et al. 2023), it is necessary to model the X-ray emission from XRBs accurately.

The metallicity dependence of X-ray outputs from high-mass XRBs (HMXBs; reviewed by, e.g. Walter et al. 2015; Kretschmar et al. 2019; Fornasini, Antoniou & Dubus 2023), which dominate the cosmic XRB luminosity density at z ≳ 3 (Fragos et al. 2013a), is particularly important in the early Universe when the metal content of XRB host galaxies evolves rapidly (e.g. Wise et al. 2012; Johnson, Dalla Vecchia & Khochfar 2013; Xu, Wise & Norman 2013; Pallottini et al. 2014; Liu & Bromm 2020; Ucci et al. 2023), and the X-rays from active galactic nuclei are subdominant (see e.g. Fragos et al. 2013b). In fact, X-ray observations of nearby and distant (up to z ∼ 2) galaxies (e.g. Antoniou et al. 2010, 2019; Basu-Zych et al. 2013; Prestwich et al. 2013; Douna et al. 2015; Antoniou & Zezas 2016; Brorby et al. 2016; Lehmer et al. 2016, 2019, 2021, 2022; Aird, Coil & Georgakakis 2017; Fornasini et al. 2019; Fornasini, Civano & Suh 2020; Riccio et al. 2023) and theoretical predictions by binary population synthesis (BPS) of XRBs (e.g. Linden et al. 2010; Fragos et al. 2013b; Sartorio et al. 2023) all suggest a strong metallicity dependence of the X-ray output from HMXBs. This drives the redshift evolution of the scaling relation between X-ray luminosity and star formation rate (SFR) and can have significant impact on the 21-cm signal (e.g. Kaur et al. 2022). For instance, Fragos et al. (2013b) provides fitting formulae for the X-ray luminosity of HMXBs per unit SFR as a function of metallicity from BPS simulations (Fragos et al. 2013a). In their case, the X-ray luminosity increases by a factor of ∼6 from solar metallicity to ∼1 per cent solar, consistent with the trend seen in observations. Sartorio et al. (2023) derives the X-ray outputs of XRBs from metal-free stars, i.e. the so-called Population III (Pop III). They find that, in optimistic cases, the X-ray emission of Pop III XRBs can be significantly stronger (up to a factor of 40) compared with that of XRBs from metal-enriched stars predicted by Fragos et al. (2013b).

However, the aforementioned BPS studies only consider the XRBs powered by Roche lobe overflow (RLO) and (spherical) stellar winds but ignore an important type of HMXBs, Be-XRBs, likely due to their transient nature. A Be-XRB is made of a compact object and a rapidly rotating, massive (⁠|$\gtrsim 6\ \rm M_{\odot }$|⁠), main sequence (MS) star (reviewed by, e.g. Reig 2011; Rivinius, Carciofi & Martayan 2013; Rivinius 2019). Here, the massive star is typically of spectral type B (and O) with an rotation velocity above ∼70 per cent of the equatorial Keplerian limit and shows Balmer emission lines, which can be well explained by a viscous decretion disc² (VDD) around the star. This VDD is formed by materials ejected from the star due to redistribution of angular momentum caused by fast rotation³ (Rivinius, Carciofi & Martayan 2013). In most observed Be-XRBs, neutron stars (NSs) are identified as the compact companion, although there are two binaries, MWC 656⁴ and AS 386, that contain a Be star and a black hole (BH) candidate but show very faint X-ray emission (Casares et al. 2014; Munar-Adrover et al. 2014; Grudzinska et al. 2015; Khokhlov et al. 2018; Zamanov et al. 2022), and one Be-XRB that contains a white dwarf (Swift J011511.0-725611, Kennea et al. 2021). The X-ray output of Be-XRBs is dominated by X-ray outbursts produced by strong accretion of the compact object from the VDD, which typically occurs close to periastron (Okazaki 2001; Okazaki & Negueruela 2001) and/or when the compact object crosses a tidally warped (eccentric) VDD (Okazaki, Hayasaki & Moritani 2013; Franchini & Martin 2021).

Be-XRBs make up the largest class of HMXBs identified in observations (Fornasini, Antoniou & Dubus 2023), especially in metal-poor environments. In the Milky Way (MW), 74 Be-XRBs have been found among the total 152 known HMXBs (Fortin et al. 2023). In the Large Magellanic Cloud (LMC) at approximately half solar metallicity, there are 33 Be-XRBs among the 40 confirmed HMXBs (Antoniou & Zezas 2016), while in the Small Magellanic Cloud (SMC) at about one quarter solar metallicity, 69 X-ray pulsars are identified as Be-XRBs among the 121 HMXB candidates (Coe & Kirk 2015; Haberl & Sturm 2016). The HMXB population of M33 is also dominated by Be-XRBs (Lazzarini et al. 2023). Besides, theoretical models find that Be-XRBs can be an important component of the X-ray luminosity function of HMXBs in the MW (Zuo, Li & Gu 2014; Misra et al. 2023b). However, previous BPS studies of Be-XRBs (e.g. Zhang, Li & Wang 2004; Belczynski & Ziolkowski 2009; Linden et al. 2009; Shao & Li 2014, 2020; Zuo, Li & Gu 2014; Vinciguerra et al. 2020; Xing & Li 2021; Misra et al. 2023b) focus on the cases of solar and SMC metallicities in which they either do not model the X-ray emission or use rough estimates and empirical scaling laws to characterize the VDDs and X-ray outbursts of Be-XRBs (e.g. Dai, Liu & Li 2006; Coe & Kirk 2015; Klement et al. 2017). The overall X-ray outputs from Be-XRB populations as a function of metallicity has not been investigated quantitatively, while a strong metallicity dependence is expected from the reduced mass loss at low metallicities (Linden et al. 2010; Fragos et al. 2013a; Sartorio et al. 2023). Therefore, it is crucial to include Be-XRBs in BPS models to evaluate the metallicity dependence of X-ray outputs from the entire population of HMXBs.

In light of the potential importance of Be-XRBs for the cosmic thermal and ionization history, we build a physically motivated Be-XRB model to predict the X-ray output from Be-XRBs as a function of metallicity with BPS. Inspired by the recent advancements in hydrodynamic simulations of VDDs in Be-XRBs (Okazaki, Hayasaki & Moritani 2013; Panoglou et al. 2016; Cyr et al. 2017; Brown et al. 2018, 2019; Suffak, Jones & Carciofi 2022), our model fully captures for the first time the dependence of X-ray outburst properties (i.e. strength and duty cycle) on stellar and orbital parameters of Be-XRBs by combining simulation results (Brown et al. 2019) with VDD properties inferred from observations of Be stars (Vieira et al. 2017; Rímulo et al. 2018). In this paper, we focus on the absolute metallicity⁵ range Z ∈ [10⁻⁴, 0.03] where observational constraints for HMXBs are available. Our model can also be applied to more metal-poor regimes (e.g. Z ≲ 10⁻⁶ for Pop III stars) that are likely more important at Cosmic Dawn (Sartorio et al. 2023).

The paper is structured as follows. In Section 2, we provide an overview of our method and discuss the setup of BPS parameters, binary sample, and initial conditions. In Section 3, we explain our physically motivated model for the identification and characterization of Be-XRBs. In Section 4, we build a empirical model for the X-ray spectra of Be-XRBs. In Section 5, we present our predictions on the formation efficiency (Section 5.1), mass and orbital parameter distributions (Section 5.2), and X-ray outputs of Be-XRBs (Section 5.3), focusing on how they evolve with metallicity. Finally, we summarize our main findings in Section 7 and discuss their caveats and our outlook to future work in Section 6. The key physical quantities used in this paper are summarized in Table 1.

2 BINARY POPULATION SYNTHESIS

We add a new module for the identification and characterization of Be-XRBs to the BPS code binary_c (Izzard et al. 2004, 2006, 2009; Izzard et al. 2017; Izzard & Halabi 2018; Izzard & Jermyn 2023; Mirouh et al. 2023; Hendriks et al. 2023; Hendriks & Izzard 2023b; Yates et al. 2023), which simulates the evolution of stars in each binary and the binary orbit governed by binary interactions (e.g. mass transfer and tidal effects) and stellar evolution processes such as winds and supernovae (SNe). We evolve large populations of binaries from zero-age main-sequence for 15 Gyr in the (absolute) metallicity range Z ∈ [10⁻⁴, 0.03] with randomly sampled initial binary properties through the python interface binary_c-python (Hendriks & Izzard 2023a) of binary_c. The Be-XRB model is explained in detail in the next Section 3. An X-ray spectral model based on observations (Section 4) is applied to the Be-XRB populations generated by binary_c in post-processing to calculate their X-ray outputs.

As shown in previous BPS studies (e.g. Vinciguerra et al. 2020; Xing & Li 2021), the main channel of Be-XRB formation is expected to be stable mass transfer during the MS and Hertzsprung gap (HG) phases, where the initial secondary star grows by accretion from the initial primary star and is meanwhile spun up to become an O/Be star. Thereafter, if the initial primary star collapses into a compact object when the O/Be star is still on MS, and the system remains bound, we can obtain a Be-XRB. Therefore, formation of Be-XRBs is sensitive to binary stellar evolution parameters governing the stability and efficiency of mass transfer, angular momentum loss, as well as natal kicks of SNe (see e.g. Shao & Li 2014; Vinciguerra et al. 2020; Xing & Li 2021). The initial binary properties may also play an important role.

In this work, we use the standard BSE models (Hurley, Tout & Pols 2002) with the default setup of binary_c⁶ with an updated stellar wind model from Schneider et al. (2018) and Sander & Vink (2020) as well as a special treatment of mass-transfer efficiency (Section 2.1). In addition to the mass-transfer efficiency, we briefly explain our choices of select BSE parameters that are important for Be-XRBs in Section 2.2 according to the default setup of binary_c detailed in Izzard et al. (2017). It is shown in Appendix  A that our choices of BSE parameters, combined with standard initial conditions of binary stars (Section 2.3), can reproduce the population of observed Be-XRBs in the SMC at the metallicity Z_SMC = 0.0035 (Davies et al. 2015). For simplicity, we assume that the BSE parameters and initial conditions do not evolve with metallicity for Z ∈ [10⁻⁴, 0.03].

The BSE models used here keep track of the spin evolution of each star regulated by mass loss, accretion, and tidal interactions, as detailed in Hurley, Tout & Pols (2002). In particular, during stable mass transfer via RLO, the accretor gains angular momentum from the accreted material, which is assumed to come from the inner edge of an accretion disc with the specific angular momentum of the circular orbit on the surface of the accretor. However, the effects of rotation on stellar evolution are not considered, which can be significant (particularly for initially FR stars) and complex, covering various aspects (e.g. mass loss, time-scales of evolution phases, stellar structure, nucleosynthesis, and remnant masses), especially at low Z, as shown in detailed stellar-evolution simulations (e.g. Ekström et al. 2012; Georgy et al. 2013; Choi, Conroy & Byler 2017; Groh et al. 2019; Murphy et al. 2021). Such effects may also be important in the modelling of Be-XRBs, as FR O/B stars are involved by definition. However, it is beyond the scope of this work to take into account these effects because the detailed mechanisms that connect the formation and properties of VDDs (i.e. the so-called ‘Be phenomenon’) with stellar evolution processes are still unresolved (Rivinius, Carciofi & Martayan 2013).

For simplicity, we ignore the mass growth of compact objects via accretion in Be-XRBs, so that our Be-XRB module does not affect binary stellar evolution. This approximation is justified by the fact that VDDs are very light structures (Rivinius 2019) that cannot supply much mass to compact objects. We have verified that among all compact objects in the Be-XRBs simulated in this work, the mass accreted from the VDD is less then a few per cent of the initial mass in most (⁠|$\gtrsim 99{{\ \rm per\ cent}}$|⁠) cases and remains below |$50{{\ \rm per\ cent}}$| of the initial mass in the most extreme systems with massive VDDs.

2.1 Mass-transfer efficiency

It is found by Vinciguerra et al. (2020) using the compas code (Riley et al. 2022) that efficient accretion during stable mass transfer is required to reproduce the observed orbital period distribution of Be-XRBs in the SMC.⁷ Since compas is also based on the BSE models, we expect that an enhancement of mass-transfer efficiency with respect to the default prescription of binary_c is necessary in our case. Therefore, we set the mass-transfer efficiency parameter, i.e. the ratio of the accreted mass-to-mass lost by the donor, as

Here, |$\dot{M}_{\rm donor}$| is the mass loss rate of the donor, β_thermal is the maximal mass-transfer efficiency defined with respect to the maximal steady-state mass acceptance rate |$\dot{M}_{\rm acc,\max }\sim (\epsilon _{\rm g,acc}/L_{\rm acc})^{-1}\sim L_{\rm acc}R_{\rm acc}/(GM_{\rm acc})\sim M_{\rm acc}/t_{\rm KH,acc}$|⁠, which is limited by the thermal (Kelvin–Helmholtz) time-scale |$t_{\rm KH}\sim GM_{\rm acc}^{2}/(R_{\rm acc}L_{\rm acc})$| of the accretor, given ϵ_{g, acc} ∼ GM_acc/R_acc the specific energy carried by the in-falling matter that needs to be radiated away, the luminosity L_acc, mass M_acc, and radius R_acc of the accretor. binary_c adopts a conservative choice β_thermal = 1 by default, while here we use β_thermal = 30, which is approximately the value required to match observations inferred by Vinciguerra et al. (2020). This rather large value of β_thermal captures the variation of t_KH during mass transfer by expansion and increase of luminosity (Paczyński & Sienkiewicz 1972; Hurley, Tout & Pols 2002; Vinciguerra et al. 2020). Similarly, we also increase the upper limit on β from the dynamical time-scale of the accretor as well as the thermal and dynamical time-scales of the donor by a factor of 30. Although such efficient mass transfer is required to reproduce the observed Be-XRBs in the SMC (see Appendix  A), we find by numerical experiments that the total X-ray output from Be-XRBs is insensitive to β_thermal. The reason is that the total X-ray output is dominated by luminous Be-XRBs mainly on eccentric orbits (e ≳ 0.1) whose progenitor primary stars only undergo weak mass loss (i.e. β = 1 for |$\dot{M}_{\rm donor}\ll \dot{M}_{\rm acc,\max }$|⁠, independent of β_thermal), as discussed in Section 5.

2.2 Other key BSE parameters

In addition to mass-transfer efficiency, the properties of Be-XRBs are also expected to depend on the prescriptions for mass transfer stability, angular momentum loss, remnant masses, and SN natal kicks (Vinciguerra et al. 2020). We plan to explore their effects in the future (see Section 6). Here, we briefly describe the default choices for these parameters adopted in our work for binary_c.

Mass transfer is stable when M_accretor/M_donor > q_crit at the onset of mass transfer given the critical mass ratio q_crit. We adopt q_crit = 5/8, 1/3, and 1/4 in the hydrogen MS, helium MS, and HG phases for both hydrogen and helium burning of the donor, respectively. In the giant phase, we use the prescription in Hurley, Tout & Pols (2002, see their section 2.6). Mass transfer beyond what is allowed by the mass-transfer efficiency (equation 1) is lost from the system. We adopt the fast (also called Jeans) model (Huang 1963) to calculate the angular momentum loss in this process, i.e. the specific angular momentum carried by the lost mass is equal to the specific angular momentum of the donor. We have verified by numerical experiments using the isotropic re-emission model (Soberman, Phinney & van den Heuvel 1997), in which the lost material carries the specific angular momentum of the accretor, that the prescription of angular momentum loss has minor effects on our results (with |$\lesssim 20{{\ \rm per\ cent}}$| and |$\lesssim 50{{\ \rm per\ cent}}$| changes in the total X-ray output for Z ≲ 0.01 and Z ∼ 0.01–0.03, respectively). The reason is that under the high mass-transfer efficiency described in Section 2.1, the mass/angular momentum loss during stable mass transfer has little impact on the orbital parameters (and luminosities) of Be-XRBs, which are more sensitive to stellar winds and SN natal kicks (see below and Section 5.3.1).

The masses of compact object remnants are determined by the CO core masses of progenitors using the original BSE models in Hurley, Pols & Tout (2000) and Hurley, Tout & Pols (2002). We apply natal kicks to Type II and Ib/c SNe that follow a Maxwellian distribution with a dispersion of |$\sigma _{\rm kick}=190\ \rm km\ s^{-1}$| (Hansen & Phinney 1997), while electron-capture SNe have no natal kicks. The latter typically happen to highly stripped stars in progenitor binaries of Be-XRBs (see Section 5.2), which are expected to have weak natal kicks (⁠|$\lesssim 30\ \rm km\ s^{-1}$|⁠, see the discussion in section 3.1 of Vinciguerra et al. 2020). Therefore, we use zero natal kicks for simplicity.

2.3 Binary sample and initial conditions

To construct the input catalogue of binary stars, we sample N_B = 3 × 10⁵ binaries randomly from widely used distributions of mass and orbital parameters. To be specific, the primary stellar mass M₁ is drawn from the Kroupa (2001) initial mass function (IMF) in the mass range of |$[5,100]\ \rm M_{\odot }$| and the mass ratio q ≡ M₂/M₁ is generated from a uniform distribution in the range q ∈ [0.1 M_⊙/M₁, 1]. Here, we only consider binaries with |$M_{1}\ge 5\ \rm M_{\odot }$| because only massive primary stars can form the compact objects considered in our Be-XRB model⁸ (see Section 3.1). In the way, the binary stars in our catalogue only make up a small fraction of the whole underlying stellar population. Following Misra et al. (2023b), we assume that the whole stellar population is made of |$70{{\ \rm per\ cent}}$| binary stars (Sana et al. 2012) and |$30{{\ \rm per\ cent}}$| single stars. For the whole stellar population, the single stars and the primary stars in binaries also follow the Kroupa (2001) IMF but in the range |$M_{1}\in [0.01,100]\ \rm M_{\odot }$|⁠, and the mass ratio distribution for the entire binary star population is uniform in q ∈ [0.01 M_⊙/M₁, 1]. Under these assumptions,⁹ we estimate that the total mass of stars in our binary sample accounts for f_sample = 0.2 of the total mass of the whole stellar population, using the method in appendix A of Bavera et al. (2020). Here, f_sample serves as a normalization factor for the calculation of X-ray outputs per unit stellar mass or SFR. The mass of the whole stellar population corresponding to our binary sample is |$M_{\rm tot}=f_{\rm sample}^{-1}\sum _{i}^{N_{\rm B}}(M_{1,i}+M_{2,i})\sim 3\times 10^{7}\ \rm M_{\odot }$|⁠.

For orbital parameters, by default, we follow Izzard et al. (2017) to draw the initial semimajor axis a from a log-flat distribution for |$a\in [3,10^{4}]\ \rm R_{\odot }$| and the initial eccentricity e from a thermal distribution for e ∈ [0, 1). We assume no correlations between a, e and masses of stars, while evidence of such correlations has been found in observations (e.g. Moe & Di Stefano 2017). We also consider an alternative model in which we draw the initial orbital period P_orb (in the unit of day) from a hybrid distribution |$f_{P_{\rm orb}}\equiv {dN}/{d\log P_{\rm orb}}$| based on observations of low-mass (Kroupa 1995) and massive (Sana et al. 2012) stars, motivated by the ideas in Izzard et al. (2017, see their appendix B1) and Sartorio et al. (2023, see their section 3.2.2):

where |$\tilde{m}\equiv M_{1}/M_{\rm O}$|⁠, given |$M_{\rm O}=16\ \rm M_{\odot }$| as the minimum mass of O stars, and (Kroupa 1995; Sana et al. 2012)

It turns out that the results in this case are very similar to those of the default model (see Appendix  B). Therefore, we only show the results of the default model in our main text. We defer a more detailed investigation of initial binary parameters to future work.

Finally, another initial condition parameter that can be important for Be-XRBs is the initial stellar rotation velocity v_{rot, 0}, which can be characterized by the parameter W₀ ≡ v_{rot, 0}/v_{Kep, 0} given the initial Keplerian velocity v_{Kep, 0} at the stellar equator. The reason is that rapid rotation is required to make O/Be stars, which is also used in our model to identify Be-XRBs (see Section 3.1). The chance of forming O/Be stars is expected to be higher for stars with faster initial rotation. Here, we consider two models for W₀. In our slowly rotating (SR) model, we adopt the fit formula for v_{rot, 0} from Hurley, Pols & Tout (2000, see their sec. 7.2) based on the MS data in Lang (1992):

given the initial stellar mass M. In this case, W₀ is a function of M and metallicity (which determines the initial stellar radius and v_{Kep, 0} given M). In our fast-rotating (FR) model, we set W₀ = 0.9 for all stars to obtain an upper limit on the formation efficiency and X-ray output of Be-XRBs. The FR model can also be regarded as the asymptotic situation when we decrease metallicity, since more metal-poor stars are more likely to be FR (e.g. Ekström et al. 2008; Bastian et al. 2017; Schootemeijer et al. 2022).

3 BE-XRB MODEL

3.1 Identification of Be-XRBs

Inspired by previous BPS studies on Be-XRBs (e.g. Zhang, Li & Wang 2004; Belczynski & Ziolkowski 2009; Linden et al. 2009; Shao & Li 2014, 2020; Zuo, Li & Gu 2014; Vinciguerra et al. 2020; Xing & Li 2021; Misra et al. 2023b), we identify a binary as in the Be-XRB phase with the following criteria:

• The binary is made of a massive MS (O/B) donor star with |$M_{\star }\gt 6\ \rm M_{\odot }$|⁠, and a compact object (NS/BH) companion with a mass |$M_{\rm X}\gt 1.29\ \rm M_{\odot }$|⁠. We ignore Be-XRBs with white dwarfs for simplicity considering their faintness and rareness (Kennea et al. 2021). We set the mass threshold for donor stars at |$6~\rm M_{\odot }$| as a conservative estimate of the minimum mass of Be stars in Be-XRBs (Hohle, Neuhäuser & Schutz 2010), which is larger than the minimum mass of single Be stars (⁠|$\sim 3\ \rm M_{\odot }$|⁠, Vieira et al. 2017). This choice is supported by the fact that only early spectral types of Be stars (e.g. no later than B5 in the Coe & Kirk 2015 SMC catalogue) that are expected to be massive have been found in observations of Be-XRBs (Antoniou et al. 2009; Reig 2011; Maravelias et al. 2014; Shao & Li 2014).

• There is a VDD around the donor star, which we assume to be present when the following conditions are satisfied:

  • The donor star is FR

    with W ≡ v_rot/v_Kep > 0.7 to eject mass that can potentially settle into a VDD, where v_rot is the rotation velocity and v_Kep is the Keplerian velocity at the equator of the donor star. Here, we adopt the minimum rotation rate W_crit = 0.7 for decrection disc formation suggested by Rivinius, Carciofi & Martayan (2013) based on observations.

  • The orbital period is not too small, i.e. P_orb > 7 d, otherwise the VDD cannot form due to tidal forces from the companion (Panoglou et al. 2016, 2018; Rivinius 2019). All of the ∼170 Be-XRBs detected so far have P_orb > 10 d (Raguzova & Popov 2005; Coe & Kirk 2015; Antoniou & Zezas 2016) except for one object in the SMC, [MA93] 798, with P_orb ∼ 0.7–2.7 d (Schmidtke, Cowley & Udalski 2013).

• The VDD overfills the Roche lobe of the O/B star at periastron to allow accretion by the compact object from the VDD: R_out > R_L1(a, e, q_⋆), where (Eggleton 1983)

  $$\begin{eqnarray} R_{\rm L1}(a,e,q_{\star })= \frac{0.49q_{\star }^{2/3}(1-e)a}{0.6q_{\star }^{2/3}+\ln (1+q_{\star }^{1/3})} \end{eqnarray}$$

  (6)

  is the Roche lobe radius given the semimajor axis a, eccentricity e and mass ratio q_⋆ ≡ M_⋆/M_X, and R_out is the effective disc boundary.

  In previous studies, R_out is either fixed to a typical value (e.g. |$100\ \rm R_{\odot }$| in Misra et al. 2023b), or set to the average tidal truncation radius (Zhang, Li & Wang 2004; Xing & Li 2021)

  $$\begin{eqnarray} R_{\rm trunc}=f_{\rm trunc}a\ ,\quad f_{\rm trunc}=N_{\rm trunc}^{-2/3}(1+q_{\rm X})^{-1/3}\ , \end{eqnarray}$$

  (7)

  given q_X ≡ M_X/M_⋆, and N_trunc = 3 assuming a typical viscosity parameter =0.63 (Rímulo et al. 2018).¹⁰ The latter definition excludes Be-XRBs with low eccentricities e ≲ 0.1 (see e.g. fig. 6 of Xing & Li 2021), even though such systems have been found in observations (e.g. CPD-29 2176, Richardson et al. 2023). The reason is that tidal truncation at R_trunc is not an absolute cutoff. What happens is that materials accumulate within R_trunc, and the disc density profile becomes much steeper beyond R_trunc than within R_trunc (Okazaki et al. 2002; Panoglou et al. 2016), such that accretion is still possible, although weaker, when R_trunc < R_L1(a, e, q_⋆). Such weak accretion beyond the truncation radius can also explain the presence of persistent low-luminosity Be-XRBs in observations (e.g. Sguera et al. 2023).

  In light of this, we consider an optimistic and physically motivated definition of VDD boundary as the radius beyond which gas flows in the disc become subsonic (Krtička, Owocki & Meynet 2011):

  $$\begin{eqnarray} R_{\rm crit}=0.3(v_{\rm Kep}/c_{s})^{2}R_{\star }\ \end{eqnarray}$$

  (8)

  according to the equatorial radius R_⋆ ≡ R_eq of the O/B star, the sound speed |$c_{s}=\sqrt{2k_{\rm B}T_{\rm d}/m_{\rm p}}$| in the ionised isothermal VDD of a temperature T_d = 0.6T_eff (Carciofi & Bjorkman 2006) given the stellar effective temperature T_eff, where k_B is the Boltzmann constant and m_p is proton mass. In this case, the tidal truncation effect is considered in the calculation of peak accretion rates (see Section 3.2.1). Given this optimistic definition of VDD boundary, we are able to reproduce the nearly circular (e = 0.06 ± 0.06) Be-XRB, CPD-29 2176, from progenitor binaries of stars in the mass range |$M_{1,2}\sim 10-12\ \rm M_{\odot }$| with weak SN natal kicks, similar to the progenitors identified in the BPASS models (see their table 2 Richardson et al. 2023).

• The O/B star itself does not fill the Roche lobe at the periastron: R_⋆ < R_L1(a, e, q_⋆). Otherwise, the system will be classified as an RLO XRB (Reig 2011).

3.2 X-ray outbursts of Be-XRBs

To the first order (ignoring the contributions from quiescent phases), X-ray emission of a Be-XRB can be described by (1) the bolometric luminosity of accretion flows around the compact object |$L_{\rm bol}=\epsilon \dot{M}_{\rm acc}c^{2}$| during outbursts where ϵ is the radiative efficiency and |$\dot{M}_{\rm acc}$| is the peak accretion rate, and (2) the duty cycle f_duty, i.e. the effective fraction of time the binary spends in X-ray outbursts during which the average luminosity is L_bol. Previous studies (e.g. Zuo, Li & Gu 2014; Misra et al. 2023b) usually adopt empirical scaling laws or typical values for L_bol and f_duty, which do not fully take into account the dependence of X-ray emission on stellar and orbital properties (e.g. eccentricity) of Be-XRBs. Here, we fully capture such dependence¹¹ with a physically motivated X-ray outburst model. To be specific, our model adopts simulation results calibrated to observational data to calculate L_bol (Section 3.2.1) and also considers the classification of X-ray outbursts which, combined with observational constraints, is used to estimate f_duty (Section 3.2.2).

3.2.1 Peak accretion rate and luminosity

We start with the steady-state peak accretion rate (assumed to be identical to the gas capture rate) predicted by the hydrodynamic simulations in Brown et al. (2019), which satisfies a simple power-law |$\dot{M}_{\rm acc,sim}\propto [(1-e)a]^{-2}\dot{M}_{\rm ej}$| as shown in Fig. 1, where |$\dot{M}_{\rm ej}\propto \Sigma _{0}/\alpha$| (Carciofi & Bjorkman 2008) is the steady-state mass ejection rate given the base surface density Σ₀ and viscosity parameter α of the VDD. For simplicity, we fix α = 0.63 throughout our calculation based on the measurements by Rímulo et al. (2018), so that the simulation results can be well described by

This relation is obtained from a series of simulations for a Be-XRB made of a NS with |$M_{\rm X}=1.4\ \rm M_{\odot }$| and a Be star of |$M_{\star }=18\ \rm M_{\odot }$| and |$R_{\star }=7\ \rm R_{\odot }$| with constant mass ejection rates, covering e ∈ [0, 0.6] and P_orb ∼ 40–400 d, which will be extrapolated to broader ranges of e and a in our model. For such a Be star, we estimate the stellar luminosity as |$L_{\star }\sim 3.2\times 10^{4}\ \rm L_{\odot }$| and the disc temperature as |$T_{\rm d}\sim 2\times 10^{4}\ \rm K$|⁠, from which we derive the disc surface density as |$\Sigma _{0,\rm ref}\sim 0.015\ \rm g\ cm^{-2}$| for |$\dot{M}_{\rm ej}=10^{-10}\ \rm M_{\odot }\ yr^{-1}$| (given the volume density |$\rho _{0}\sim 5\times 10^{-13}\ \rm g\ cm^{-3}$| of the disc at the stellar surface for α = 0.63), which sets the normalization of equation (9).

In reality, mass ejection can be highly variable, even leading to disc dissipation/formation at time-scales of a few years (Reig 2011), and the viscosity parameter α can vary from system to system (Vieira et al. 2017; Rímulo et al. 2018), so that VDDs are more complex in reality than simulated by Brown et al. (2019) at steady state. Therefore, the accretion rate predicted by equation (9) should be regarded as an order-of-magnitude estimate that captures the increasing trend with decreasing pericentre distance (1 − e)a. Finally, we multiply the peak accretion rate from equation (9) by a factor of (R_trunc/R_L1)⁸ for systems with R_trunc < R_L1 to capture the steepening of disc density profile beyond R_trunc (Okazaki et al. 2002).

Next, we associate the VDD base density Σ₀ with the donor star mass M_⋆ by fitting observational data (Vieira et al. 2017; Rímulo et al. 2018). The obtained empirical scaling laws capture the increasing trend of Σ₀ with M_⋆ (Arcos et al. 2017; Klement et al. 2017; Vieira et al. 2017; Rímulo et al. 2018), as shown in Fig. 2. To be specific, we have

with ≃ 0.52 dex scatter by fitting the data of 80 Be stars observed in the MW (Vieira et al. 2017), and

of ≃ 0.17 dex scatter for 54 Be stars observed in the SMC (Rímulo et al. 2018). The observations by Rímulo et al. (2018) are likely biased towards dense discs, such that equation (11) should be regarded as an upper limit. Besides, to consider the large scatter in Σ₀ at similar M_⋆, for each Be-XRB, we draw a random number χ from a Gaussian distribution of a standard deviation σ = 0.52 (0.17) dex, and set |$\Sigma _{0}=10^{\chi }\tilde{\Sigma }_{0,\rm MW\ (SMC)}(M_{\star })$|⁠, given the prediction of the best-fit model |$\tilde{\Sigma }_{0,\rm MW\ (SMC)}(M_{\star })$| for the MW (SMC). In addition to the donor mass dependence, we also consider the metallicity Z dependence of Σ₀ with two cases. In the conservative (CS) case, we always use the MW model independent of Z, motivated by the finding that the X-ray luminosity per luminous HMXB is insensitive to metallicity for Z ∼ 0.0004–0.03 in nearby galaxies (see their fig. 5 of Douna et al. 2015), while in the optimistic (OP) case, we assume that Σ₀ increases with decreasing metallicity with linear relation between log (Σ₀) and Z from solar to SMC metallicities, i.e. Z ∼ 0.0035–0.0142, and adopt the MW model for |$Z\gt \rm Z_{\odot }=0.0142$| (Asplund et al. 2009) and the SMC model for Z < Z_SMC = 0.0035 (Davies et al. 2015).

To model the peak accretion rate and luminosity more precisely, we calibrate our model with observations of Be-XRBs in the MW. To do so, we apply the above formalism (at solar metallicity) to a randomly generated sample of 10000 NS-Be star binaries with a log-flat distribution of a in the range of |$[10\!-\!1000]\ \rm R_{\odot }$|⁠, a uniform distribution of eccentricity for e ∈ [0, 0.6] and a log-flat distribution of M_⋆ for |$M_{\star }\in [6,20]\ \rm M_{\odot }$|⁠, given fixed |$M_{\rm X}=1.4\ \rm M_{\odot }$|⁠. From this sample we select a mock population of Be-XRBs with |$P_{\rm orb}\sim 10-300\ \rm days$| and |$L_{\rm X}\gt 10^{34}\ \rm erg\ s^{-1}$| to be compared with observations. These conditions are chosen to mimic the statistics of most (∼90 per cent) well-observed Be-XRBs in the MW (Raguzova & Popov 2005; Cheng, Shao & Li 2014; Brown et al. 2018). The calibration target is the relation between the (outburst) X-ray luminosity L_X and orbital periods P_orb, derived by Dai, Liu & Li (2006, D06) based on 36 observed Be-XRBs from Raguzova & Popov (2005):

The best-fitting model indicates that the typical X-ray luminosity follows |$L_{\rm X}\propto P_{\rm orb}^{-3/2} \propto a^{-9/4}$|⁠, which is similar to the a dependence in |$\dot{M}_{\rm acc,sim}\propto a^{-2}$|⁠. In light of this, we assume that L_X is proportional to the bolometric luminosity predicted by simulations with a calibration parameter ψ_X: |$L_{\rm X}= \psi _{\rm X}\epsilon \dot{M}_{\rm acc,sim}c^{2}$|⁠, given ϵ = 0.2 the typical radiative efficiency for NSs. This calibration factor captures the difference between X-ray luminosity and bolometric luminosity in observations as well as the difference between the peak accretion rates predicted by simulations and those in reality.

If the observed X-ray luminosity completely dominates the bolometric luminosity and the simulations are realistic, we should have ψ_X ≃ 1. However, we find that the empirical L_X-P_orb scaling relation (equation 12) can be reproduced by the mock population with ψ_X = 0.25 ≪ 1, as shown in Fig. 3. There are two possible reasons for this low ψ_X value: (1) The accretion rates predicted by the simulations in Brown et al. (2019) are overestimated and/or the mass ejection rates of Be stars in Be-XRBs are lower than those of Be stars at large in the MW (Vieira et al. 2017). (2) The X-ray luminosity derived from observations in fact only represents a (small) fraction of the bolometric luminosity. To capture these uncertainties, we introduce a correction factor f_corr ∈ [ψ_X, 1] (fixing ψ_X = 0.25) for the peak accretion rate |$\dot{M}_{\rm acc}$|⁠. We also include a power-law term of M_X with index ξ to account for the mass dependence, so that

where we adopt ξ = 2 for Bondi-like accretion. Then the bolometric luminosity during outbursts of a Be-XRB can be written as

Here, f_corr = 1 corresponds to the optimistic case in which the discrepancy is only caused by observational effects, while f_corr = 0.25 is the opposite end where the overestimation in simulations needs to be corrected the most. By default, we adopt f_corr = ψ_X/f_{BC, 0} = 0.5, assuming a typical bolometric correction (BC) factor f_{BC, 0} = L_X/L_bol = 0.5. The motivation is that most observations of Be-XRBs in Raguzova & Popov (2005) come from the 2–10 keV band, which typically counts for ∼50 per cent of the bolometric luminosity (see Section 4 below).

Besides, we do not cape |$\dot{M}_{\rm acc}$| at the Eddington rate

considering that Be-XRBs are promising candidates for ultra-luminous X-ray sources (ULXs, with |$L_{\rm X}\gtrsim 10^{39}\ \rm erg\ s^{-1}$|⁠, reviewed by, e.g. Kaaret, Feng & Roberts 2017; Fabrika et al. 2021; King, Lasota & Middleton 2023) and a few Be-XRBs with outburst luminosities above the Eddington luminosity |$L_{\rm Edd}=\epsilon \dot{M}_{\rm Edd}c^{2}\sim 2\times 10^{38}\ \rm erg\ s^{-1}$| for typical NSs with ϵ = 0.2 and |$M_{\rm X}=1.4\ \rm M_{\odot }$| have been observed (see table 1 of Karino 2022 and Fig. 3). Moreover, BPS studies show that XRBs (with both BH and NS accretors) can undergo episodes of highly super-Eddington (up to η ∼ 10³) mass transfer and potentially become ULXs (e.g. Marchant et al. 2017; Wiktorowicz et al. 2017, 2019, 2021; Shao, Li & Dai 2019; Shao & Li 2020; Abdusalam et al. 2020; Kuranov, Postnov & Yungelson 2020; Misra et al. 2020). It is discussed below that ULXs are important in our Be-XRB populations.

We still use |$L_{\rm bol}=\epsilon \dot{M}_{\rm acc}c^{2}$| when |$\dot{M}_{\rm acc}\gt \dot{M}_{\rm Edd}$|⁠, ignoring any possible suppression of L_bol in the super-Eddington regime by, e.g. radiation-driven winds from the accretion disc¹² (Shakura & Sunyaev 1973), so that our results should be regarded as optimistic estimates. We also ignore the beaming effects of accretion disc geometry that can boost the apparent luminosities (and reduce the observed duty cycles or number counts) of ULXs with |$\eta \equiv \dot{M}_{\rm acc}/\dot{M}_{\rm Edd}\gtrsim 8.5$| (King et al. 2001; King 2009; Lasota & King 2023). Our simple approach is motivated by the lack of features around the Eddington limits of NSs and BHs in the observed luminosity function of HMXBs, which implies that super-Eddington systems are most likely ‘normal’ XRBs similar to their sub-Eddington counterparts (Gilfanov et al. 2022).

3.2.2 Classification of X-ray outbursts and duty cycle

Now we can calculate the outburst strength by equations (13) and (14). We further classify the outbursts into two categories following the convention in observations to estimate the duty cycle (Reig 2011; Rivinius, Carciofi & Martayan 2013):

• Type I outbursts are regular, (quasi-)periodic and short-lived (∼0.1–0.3 P_orb) increases of X-ray flux by a factor of ∼10–100, peaking at or close to the pericentric passage of the compact object with X-ray luminosities |$L_{\rm X}\lesssim 10^{37}\ \rm erg\ s^{-1}$|⁠. The duty cycle is typically f_{duty, I} ∼ 0.1–0.3 (Reig 2011; Sidoli & Paizis 2018).

• Type II outbursts are major enhancements of X-ray flux, by a factor of 10³–10⁴, even reaching the Eddington limit. They do not have preferred orbital phases and last longer than Type I outbursts (up to a few orbital periods). During a Type II outburst, a radiatively efficient thin accretion disc is expected to form around the compact object, and the VDD structure can be significantly disrupted. The duty cycle is usually lower than the Type I case: f_{duty, II} ∼ 10⁻³ − 0.1 (Sidoli & Paizis 2018; Xu & Li 2019).

These two types of outbursts generally correlate with the two peaks in the observed bi-model spin period distribution of NSs in Be-XRBs, which can be divided by a critical spin period |$P_{\rm s,crit}=40\ \rm s$| (Cheng, Shao & Li 2014; Haberl & Sturm 2016; Xu & Li 2019). To explain this correlation, it is proposed by Okazaki, Hayasaki & Moritani (2013) with considerations of accretion time-scale and spin-up efficiency that the two types of outbursts experience different modes of accretion: During a Type II outburst, the NS accretes at a high rate via a radiatively efficient thin accretion disc and is spun up efficiently to have spin periods |$P_{\rm s}\lesssim 40\ \rm s$|⁠. In Type I outbursts, accretion is in the form of advection-dominated accretion flow (ADAF) resulting in low spins with |$P_{\rm s}\gtrsim 40\ \rm s$|⁠. It is further shown by Cheng, Shao & Li (2014) that disc warping plays an important role in the spin evolution of NSs, such that Type II outbursts tend to occur when NSs interact with tidally warped VDDs. Motivated by these results (and generalizing them to BHs), we assume that a Be-XRB will have Type II outbursts when two criteria for (1) tidal warping and (2) accretion rate are satisfied, as defined below.

1. Following the analysis in Cheng, Shao & Li (2014) for the power-law + Gaussian VDD model (see section 2.2 of Martin et al. 2011), the tidal warping criterion is satisfied when the disc truncation radius at periastron R_{trunc, p} = f_trunc(1 − e)a is larger than the tidal warping radius (equation 30 in Martin et al. 2011):

   $$\begin{eqnarray} R_{\rm tid}&=&\left[\frac{2\nu _{\star }(G M_{\star })^{1/2}\bar{R}_{\rm b}^{3}}{3GM_{\rm X}R_{\star }^{n-2}}\right]^{\gamma } \\ &=& a^{3\gamma }(1-e^{2})^{3\gamma /2}\left[\frac{2\alpha H_{\star }^{2}M_{\star }}{3M_{\rm X}R_{\star }^{n-1/2}}\right]^{\gamma }\ , \end{eqnarray}$$

   (16)
   where |$\nu _{\star }=\alpha H_{\star }^{2}R_{\star }^{3/2}(GM_{\star })^{1/2}$| is the (base) disc viscosity at the stellar surface, |$\bar{R}_{\rm d}=a(1-e^{2})^{1/2}$| is the average separation, |$H_{\star }=\sqrt{2}c_{s}v_{\rm Kep}^{-1}R_{\star }$| (≃ 0.04R_⋆, Wood, Bjorkman & Bjorkman 1997) is the (base) scale height of the disc at the stellar surface (Klement et al. 2017), and the power-law index γ = 2/(11 − 2n) ∼ 0.5 is given by n ∼ 3.5 which is the slope¹³ of the disc (mid-plane) density profile (see eqs. 12, 16, and 29 in Martin et al. 2011). Substituting the formula of R_tid (equation 16) to the iniquity f_trunca(1 − e) > R_tid, we have

   $$\begin{eqnarray} a\lt a_{\rm crit}\equiv \left[\frac{f_{\rm trunc}(1-e)}{(1-e^{2})^{\frac{3\gamma }{2}}}\left(\frac{3}{2\alpha }\frac{M_{\rm X}}{M_{\star }}\frac{R_{\star }^{n-1/2}}{ H_{\star }^{2}}\right)^{\gamma }\right]^{\frac{1}{{(3\gamma -1)}}}\ , \end{eqnarray}$$

   (17)

   in which f_trunc is given by equation (7). Here, we use α = 0.63 to evaluate the critical seperation a_crit ∝ α^{−γ/(3γ − 1)} ∼ α⁻¹ in equation (17), motivated by the finding in Cheng, Shao & Li (2014) that the observed populations of Be-XRBs with low (P_s > 40 s) and high (P_s < 40 s) spins, roughly corresponding to Type I and II outbursts, can be well divided by the tidal warping criterion with α ∼ 0.5–1. However, the adopted value of α here is much lower than the viscosity parameter for vertical shear α₂ = 2.66 considered in Martin et al. (2011). It is shown below that our choice of α is justified by comparing the mock population of Be-XRBs with observations (Fig. 4). The discrepancy here may be caused by the fact that the VDD flares (reaching H/R ≳ 0.1 at |$R\gtrsim 10 R_{\star }\sim 100\ \rm R_{\odot }$|⁠) while the value in Martin et al. (2011) is derived for thin (H/R ≪ 1), flat discs (Ogilvie 1999; Lodato & Price 2010). The disc flaring may reduce the viscosity for vertical shear and enhance vertical diffusion, making the disc more vulnerable for tidal warping (with R_tid smaller by a factor of ∼2).

2. The accretion rate criterion can be written as

   $$\begin{eqnarray} \eta \equiv \dot{M}_{\rm acc}/\dot{M}_{\rm Edd}\gt \eta _{\rm crit}\ . \end{eqnarray}$$

   (18)

   Here we adopt the typical radiative efficiencies ϵ = 0.2 for NSs and ϵ = 0.1 for BHs to calculate the Eddington rate |$\dot{M}_{\rm Edd}$| (equation 15). We expect the transition Eddington ratio to be in the range η_crit ∼ 0.05–0.2, where the upper limit is adopted in Okazaki, Hayasaki & Moritani (2013) to explain the outburst strengths of Type I and II in simulations, while the lower limit is consistent with the theoretical thin-disc formation criterion η > 0.07α adopted in Takhistov et al. (2022) based on Pringle (1981), given the viscosity parameter α = 0.63 in our case. We set η_crit = 0.2f_corr, because with this choice the L_X distributions of Type I and II outbursts from the mock population are generally consistent with those of observed Be-XRBs corresponding to the two peaks of spin period distribution at P_s > 40 s and P_s < 40 s (see e.g. fig. 3 of Cheng, Shao & Li 2014), as shown in Fig. 4. Since our mock population of Be-XRBs is not meant to fully capture the statistics of observed Be-XRBs complied by Cheng, Shao & Li (2014), it does not reproduce the quasi-bi-modal feature of the observed L_X distribution.

Given the above classification, we combine the optimistic duty cycles |$\hat{f}_{\rm duty}$| in observations of the two types of outbursts, |$\hat{f}_{\rm duty,I}=0.3$| and |$\hat{f}_{\rm duty,II}=0.1$|⁠, with a physical limit |$f_{\rm duty,\max }=\dot{M}_{\rm ej}/\dot{M}_{\rm acc}$| to estimate the (average) duty cycle as |$f_{\rm duty}=\min (\hat{f}_{\rm duty},f_{\rm duty,max})$| where

is the mass ejection rate for α = 0.63 based on the results from Brown et al. (2019).¹⁴ The limit |$f_{\rm duty,\max }$| captures the simple requirement that the compact object does not accrete more than what is ejected from the O/B star. When |$f_{\rm duty,\max }\lt \hat{f}_{\rm duty}$|⁠, we assume that the compact object is able to accrete all materials ejected by the O/B star during X-ray outbursts, despite the fact that for classical O/Be stars only a small fraction (∼0.01) of the ejected materials is expected to settle into the disc according to the standard steady-state VDD model (e.g. Haubois et al. 2012; Rímulo et al. 2018), while the majority will fall back to the star. This optimistic assumption is required to explain the observed high duty cycles (up to 0.3, Reig 2011; Sidoli & Paizis 2018). It means that mass replenishment of VDDs is much more efficient¹⁵ for O/Be stars in Be-XRBs than predicted by the standard VDD model (for O/Be stars in isolation), which is supported by the shorter disc time-scales of Be stars in Be-XRBs compared with single Be stars in observations (Reig 2011).

With this model, we find that most Be-XRBs in the mock population do satisfy f_{duty, I} ∼ 0.1–0.3 and f_{duty, II} ∼ 0.01–0.1, which is generally consistent with observations (Reig 2011; Sidoli & Paizis 2018; Xu & Li 2019). In the mock population, the evolution of binary and stellar parameters during the Be-XRB phase is ignored, whereas in the BPS runs such evolution can change the outburst type of a Be-XRB. Be-XRBs with both types of outbursts also exist in observations. We classify the systems that experience both types of outbursts as Type I/II. Besides, when |$\dot{M}_{\rm acc}$| is comparable to |$\hat{f}_{\rm duty}^{-1}\dot{M}_{\rm ej}$| (such that |$f_{\rm duty}\lesssim \hat{f}_{\rm duty}$|⁠), significant disruption of the VDD by the compact object is expected to happen, such that the system will show Type II features especially with low f_duty. Therefore, we also regard the Be-XRBs classified as Type I according to equations (17) and (18) with f_duty < 0.1 as Type I/II in post-processing. Finally, we count Type I/II systems into the general Type II category, which refers to all Be-XRBs that once undergo major outbursts with high accretion rates, disc warping and/or disruption. Recent observations find that the outburst behaviours of Be-XRBs are likely more diverse and complex than the conventional two types (Sidoli & Paizis 2018). Nevertheless, our consideration of the outburst type only affects the final X-ray output indirectly by the optimistic duty cycle |$\hat{f}_{\rm duty}$|⁠, and we have verified by numerical experiments that our results are insensitive to outburst classification. The reason is that the majority (∼60–70 per cent) of X-ray emission comes from systems with |$f_{\rm duty}\simeq f_{\rm duty,max}\lesssim \hat{f}_{\rm duty}$| in all Be-XRB populations considered here.

4 X-RAY SPECTRAL MODEL

Once L_bol is known, we only need to determine the spectral shape to obtain the full (intrinsic) spectral energy distribution (SED) of X-ray outbursts. For simplicity, we consider three regimes of accretion rates: low-hard (LH, η < 0.05), high-soft (HS, η ∼ 0.05–2) and super-Eddington (SE, η > 2), for both NSs and BHs. Motivated by the ideas in Fragos et al. (2013a), we find the typical spectral shape in each regime by fitting simple spectral models to the BC factors measured in observations (e.g. McClintock & Remillard 2006; Wu et al. 2010; Anastasopoulou et al. 2022) for select energy bands.¹⁶ Here we consider the photon energy range E ∼ 10⁻⁴−10⁴ keV assuming that this contains most of the power.

We take the recent measurements of BC factors of the energy bands E ∼ 0.5–2, 2–10, and 12–15 keV by Anastasopoulou et al. (2022, see their tables 4, A1, A6, A9, and A10) for both NS and BH HMXBs¹⁷ as references. We estimate the typical value (and uncertainty) of the BC factor for each band in each regime, as shown in Fig. 5. These estimates are used to construct reference SEDs that serve as the targets of spectral fitting. Since soft X-ray photons in the E ∼ 0.5 − 2 keV band are mostly responsible for X-ray heating of the IGM (Das et al. 2017), we increase the weight of this band by a factor of 4 in the fitting process to better reproduce the corresponding BC factor.

For simplicity, we assume that for both NSs and BHs the final spectrum is made of two components: an input spectrum and a spectrum of inverse Compton scattered photons (by corona electrons) with a power-law tail in the high-energy end. Following Sartorio et al. (2023), we use the SIMPL-1 Comptonizon model in Steiner et al. (2009, see their eqs. 1 and 4) to connect the upscattered component with the input spectrum |$L_{\nu ,\rm in}$| via two parameters: the fraction f_scatter of photons in the input (photon number) spectrum |$n_{\rm in}(E)\equiv dN_{\rm in}/dE=(hE)^{-1}L_{\nu ,\rm in}$| that are upscattered, and the power-law index Γ in the spectrum of upscattered photons:

where G(E, E₀) = (Γ − 1)(E/E₀)^−Γ/E₀ is the Green’s function of inverse Compton scattering. The final specific luminosity can then be derived from the output photon number spectrum with |$n_{\rm out}(E)\equiv dN_{\rm in}/dE=(hE)^{-1}L_{\nu ,\rm out}$|⁠.

For the input spectrum, we adopt the black body (BB) spectrum for NSs assuming that the majority of X-ray radiation is produced at hot spots on the NS surface from inflows channelled by magnetic fields. The top row of Fig. 6 shows the resulting best-fit spectral models in the three regimes with k_BT ∼ 0.4−0.9 keV, f_scatter ∼ 1 and Γ ∼ 2.5−2.7.

While for BHs, we use the thin disc (TD) model¹⁸ (Pringle 1981; Takhistov et al. 2022)

where T_max = 0.488T_i, T_o ≃ T_i(R_i/R_o)^3/4 and ν_max = k_BT_max/h, given T_i as the temperature at the inner edge (R_i = 6GM_X/c²) of the TD, and T_o as the temperature at the outer disc boundary R_o ∼ R_L1(a, e, q_X). Since very close binaries (P_orb ≲ 7 d) and RLO are forbidden for Be-XRBs (Panoglou et al. 2016, 2018; Rivinius 2019) by tidal forces that can slow down the rotation of donor stars, we have |$R_{\rm o}\gtrsim 2\ \rm R_{\odot }$| and k_BT_o ≲ 10⁻³ keV in most cases. Therefore R_o is unimportant for X-rays that we are concerned with (E > 0.1 keV) and the input TD spectrum is controlled by a single parameter T_i during the fitting process. The best-fit models for |$M_{\rm X}=4\ \rm M_{\odot }$| and |$R_{\rm o}=10\ \rm R_{\odot }$| are shown in the bottom row of Fig. 6 with k_BT_i ∼ 0.7–1.5 keV, f_scatter ∼ 0.3–1 and Γ ∼ 2.1–2.8.

5 RESULTS

Combining the two choices for initial rotation parameter W₀: SR or FR (Section 2.3), and the two models for VDD density Σ₀: CS or OP (Section 3.2), we have four models to explore, which are summarized in Table 2. For each model, we consider ten cases at Z = 10⁻⁴, 3 × 10⁻⁴ 10⁻³, 0.0035, 0.005, 0.007, 0.01, 0.0142, 0.02, and 0.03, where Z = 0.0035, 0.07, and 0.0142 correspond to the situations of the SMC, LMC and MW, respectively. Since the CS and OP models are identical at Z ≥ 0.0142, we only need to run 34 BPS simulations in total. We record the time-averaged values of the Be-XRB properties including properties of donors and accretors, orbital parameters as well as X-ray outburst properties,¹⁹ which are then used to calculate the total X-ray output. Considering the short lifetimes of Be-XRBs, for simplicity, when calculating the total X-ray emission of a Be-XRB population, we assume no variation of X-ray outburst properties during the Be-XRB phase,²⁰ which is a simplification of the reality that VDDs can be highly variable structures with disc dissipation/formation at time-scales of a few years (Reig 2011; Rivinius, Carciofi & Martayan 2013). Since our purpose is to evaluate the overall X-ray output from a population of Be-XRBs, higher-order effects are expected to be unimportant. In this section we mainly show the results from the SR_CS model with f_corr = 0.5, defined as the fiducial case, because it achieves the best agreement with observations and the key trends in the metallicity dependence of Be-XRB properties are similar in the other cases. Select results for the other 3 models in Table 2 and different values of f_corr are included in Appendix  B.

5.1 Formation efficiency

We first look into the formation efficiency of active Be-XRBs, |$\mathcal {N}_{\rm X}$|⁠, as a function of Z which, considering the short lifetimes (a few to a few tens Myr) of Be-XRBs, is defined as the number of Be-XRBs in the outburst phase per unit SFR for a long enough star formation time-scale (τ_SF ≳ 100 Myr). Given N Be-XRBs predicted by a BPS run for a single-age stellar population of a total mass M_tot, the formation efficiency can be written as

where τ_i is the duration of the Be-XRB phase for binary i.

The results for all four models in Table 2 are shown in Fig. 7, where we count both the number of all Be-XRBs (thin curves with markers) and that of Be-XRBs with e > 0.1 (thick curves without markers). The latter case is meant to capture the situation assumed in previous BPS studies of Be-XRBs (e.g. Zhang, Li & Wang 2004; Xing & Li 2021) that tidal truncation of the VDD is sharp, leading to a smaller disc boundary at R_trunc (equation 7) than the one adopted in our model (R_crit, see equation 8). In both cases, |$\mathcal {N}_{\rm X}$| generally increases with decreasing Z, but the evolution is not fully monotonic for all Be-XRBs, which has a small peak at Z = 0.005 and a small dip at Z = 0.01. The general trend is driven by the stronger stellar winds at higher Z that increasingly widen binary orbits and reduce the number of stars that can become NSs and BHs. The non-monotonic features are likely caused by the complex interplay between stellar winds and mass transfer rate (see Section 5.2). We also find that nearly-circular (e ≤ 0.1) systems make up a significant fraction (⁠|$\sim 40-80{{\ \rm per\ cent}}$|⁠) of active Be-XRBs at Z ≲ 0.02. Naïvely, this seems in tension with the rareness (⁠|$\lesssim 10{{\ \rm per\ cent}}$|⁠) of such Be-XRBs in observations (Cheng, Shao & Li 2014; Sidoli & Paizis 2018). However, observations are very likely incomplete at the low eccentricity end because the measurement of e is difficult for nearly-circular binaries, and only a small fraction of observed Be-XRBs have eccentricity measurements. Besides, these binaries are typically faint with |$L_{\rm bol}\sim 10^{33}-10^{37}\ \rm erg\ s^{-1}$| due to the suppressed accretion rate from disc truncation and therefore difficult to detect. In fact, they produce much fewer X-rays compared with their more eccentric counterparts such that ignoring them has little (up to a few per cent) impact on the overall X-ray output.

Last but not least, |$\mathcal {N}_{\rm X}$| is higher in the FR models with higher initial rotation rates than in the SR models as expected. The difference between the FR and SR models is larger at higher Z but remains below |$\sim 40{{\ \rm per\ cent}}$| and becomes even comparable to the uncertainties²¹ in |$\mathcal {N}_{\rm X}$| at Z ≲ 0.0035. The overall small difference indicates that in most binaries that can potentially become Be-XRBs the (initial) secondary star will be spun up to become an O/Be star (via stable mass transfer during the MS and HG phases) regardless of its initial rotation rate. This is consistent with the scenario that all or most (young) O/Be stars are produced by mass and angular momentum transfer from companion stars (Shao & Li 2014; Hastings, Wang & Langer 2020; Hastings et al. 2021; Dodd et al. 2023; Wang et al. 2023), which is supported by observations that find a large fraction of classical Be stars with disc truncation (i.e. SED turndown; Klement et al. 2019), the lack of close Be binaries with MS companions (Bodensteiner, Shenar & Sana 2020), and a higher run-away/field frequency of O/Be stars (or fast rotators) compared with normal O/B stars (Dorigo Jones et al. 2020; Dallas, Oey & Castro 2022). On the other hand, the difference between the FR and SR models is generally larger when we focus on Be-XRBs with e > 0.1 at Z ≳ 10⁻³. The reason is that mass transfer is weaker in these systems leading to less efficient spin-up of the secondary star, as discussed below (Section 5.2). Besides, the difference is smaller at lower Z, where the secondary stars are more compact and are more easily spun up to become O/Be stars.

In addition to the general Be-XRB population, we also derive the number of luminous and ultra-luminous sources with broad-band (E ∼ 0.5–8 keV) X-ray luminosities |$L_{[0.5-8]~\rm keV}\gt 10^{38}\ \rm erg\ s^{-1}$| and |$L_{[0.5-8]~\rm keV}\gt 10^{39}\ \rm erg\ s^{-1}$| during outbursts, respectively, for which constraints from observations of HMXBs in nearby galaxies are available down to Z ∼ 0.0003 (Mapelli et al. 2010; Douna et al. 2015; Kovlakas et al. 2020; Lehmer et al. 2021). The results in our fiducial case are shown in Fig. 8. For ultra-luminous sources with |$L_{[0.5-8]~\rm keV}\gt 10^{39}\ \rm erg\ s^{-1}$|⁠, we find that the numbers of Be-XRBs predicted by our BPS runs are lower than those inferred from observations²² of all types of HMXBs (Lehmer et al. 2021) by about 40 per cent at Z ∼ 0.0003–0.02, which means that the simulated Be-XRBs can explain ∼60 per cent of all ultra-luminous HMXBs, assuming that observations are complete. For instance, we have |$\mathcal {N}_{\rm X}=0.29\ \rm M_{\odot }^{-1}\ yr$| at |$Z=\rm Z_{\odot }$| for Be-XRBs, while the observed value for all types of HMXBs is |$\mathcal {N}_{\rm X}=0.45_{-0.09}^{+0.06}\ \rm M_{\odot }^{-1}\ yr$| (Kovlakas et al. 2020). Interestingly, the predicted evolution with Z of Be-XRBs follows a similar trend as that seen in observations of all types of HMXBs at Z ∼ 0.0003 − 0.02. However, at Z ≳ 0.02 the decrease of |$\mathcal {N}_{\rm X}$| with Z is much stronger for our simulated Be-XRBs compared with the observed trend for all types of HMXBs, such that the predicted number for Be-XRBs becomes lower than the observed number for all types of HMXBs by a factor of a few at Z = 0.03. This is likely caused by strong stellar winds and poor statistics of the small number (∼300) of Be-XRBs identified from the BPS run for Z = 0.03. If this feature is true, it means that Be-XRBs play much less important roles at Z > 0.02, where wind-fed XRBs make up the majority of ultra-luminous sources (e.g. Wiktorowicz et al. 2021). The situation for sources with |$L_{[0.5-8]~\rm keV}\gt 10^{38}\ \rm erg\ s^{-1}$| (Douna et al. 2015) is similar. Although not shown here for conciseness, we find similar trends for Be-XRBs with |$L_{[0.5-8]~\rm keV}\gt 10^{37}\ \rm erg\ s^{-1}$|⁠. In this case, we have |$\mathcal {N}_{\rm X}\sim 10\ \rm M_{\odot }^{-1}\ yr$| at |$Z\sim \rm Z_{\odot }$|⁠, consistent with observations of HMXBs in nearby galaxies (Gilfanov et al. 2022). Moreover, for Be-XRBs with |$L_{[0.5-8]~\rm keV}\gtrsim 10^{35}\ \rm erg\ s^{-1}$|⁠, we predict |$\mathcal {N}_{\rm X}\sim 60-130\ \rm M_{\odot }^{-1}\ yr$| at Z ∼ 0.0035 − 0.0142, again below the observed rate |$\sim 135\ \rm M_{\odot }^{-1}\ yr$| for all types of HMXBs in nearby star-forming galaxies (Mineo, Gilfanov & Sunyaev 2012; Gilfanov et al. 2022; Lazzarini et al. 2023).

5.2 Masses and orbital parameters

For conciseness, in this section we only show the detailed statistics of Be-XRBs for the SR_CS model at three representative metallicities, Z = 10⁻⁴, 0.0035 (SMC), and 0.0142 (MW), to illustrate the general trends. Each Be-XRB i has a weight τ_i/M_tot, so that the number for each bin in the plots of distributions corresponds to the expected number of binaries in the Be-XRB phase per unit SFR (⁠|$\rm M_\odot \ yr^{-1}$|⁠).

Fig. 9 shows the mass distributions of Be-XRBs and their progenitors as well as the initial-remnant mass relation for primary stars of Be-XRB progenitors. Here, the (initial) primary star (M₁) is the progenitor of the compact object (M_X), and the (initial) secondary star (M₂) is the progenitor of the O/Be star (M_⋆). We find a positive correlation between the O/Be star mass and compact object mass, as well as between the progenitor masses, and that this correlation is stronger at lower Z. The reason is that more massive primary stars tend to have more massive remnants (i.e. higher M_X) and also require more massive secondary stars to have stable mass transfer. Both the typical primary and secondary masses increase with Z which, combined with the bottom-heavy IMF, explains the decreasing formation efficiency of Be-XRBs at higher Z (Fig. 7). This trend is caused by two effects: The formation of NSs and BHs requires more massive primary stars at higher Z due to stronger stellar winds, and the secondary stars that are more compact at lower Z are more easily spun up to become O/Be stars. These two effects further complement each other by the stability of mass transfer. Nevertheless, the high-mass tail of the progenitor mass distribution shrinks with increasing Z. This is caused by an effect that involves less massive stars given stronger winds at higher Z: Strong stellar winds from the most massive stars drive significant expansion of binary orbits so that mass transfer (which is usually required to make O/Be stars in the SR models) is suppressed, and the chance of forming close enough binaries to allow accretion from VDDs is also reduced. In fact, the trend of shrinking high-mass tail with increasing Z is much weaker in the FR models where the secondary stars rotate rapidly from the beginning and do not need to be spun up by mass transfer to become O/Be stars.

Because of stellar winds, the maximum compact object mass decreases with Z from |$M_{\rm X,\max }\simeq 4.3\ \rm M_{\odot }$| at Z = 10⁻⁴ to |$M_{\rm X,\max }\simeq 1.6\ \rm M_{\odot }$| at Z = 0.03. There are no BHs (with |$M_{\rm X}\gt 2.2\ \rm M_{\odot }$|⁠) at Z ≳ 0.0035, and the fraction of BH systems is |$\sim 1-3{{\ \rm per\ cent}}$| at lower Z. In general, our Be-XRBs are dominated by low-mass NSs with |$M_{\rm X}\sim 1.3\ \rm M_{\odot }$|⁠. The mass distribution of O/Be stars shows complex Z dependence. To demonstrate the general trends, we show the median O/Be star mass in Be-XRBs (solid curve) as well as the fraction of Be-XRBs hosting low-mass (⁠|$M_{\star }\lt 10\ \rm M_{\odot }$|⁠) O/Be stars (dashed curve) as a function of Z in Fig. 10. The median O/Be star mass increases from |$\sim 9\ \rm M_{\odot }$| at Z = 10⁻⁴ to |$\sim 16\ \rm M_{\odot }$| at Z = 0.001 and shows little evolution for Z ∼ 0.001–0.01 before dropping rapidly for Z > 0.01 down to |$\sim 7\ \rm M_{\odot }$| at Z = 0.03. The increasing trend at low Z can be explained by the increase of secondary mass with Z and higher mass transfer rates from more massive primary stars with larger radii (as |$\dot{M}_{\rm acc,\max }\propto R_{\rm acc}$|⁠) at higher Z. For Z ≳ 0.01, the formation of massive Be-XRBs from massive progenitors is increasingly suppressed by orbital expansion with stronger stellar winds at higher Z, which also reduce the masses of O/Be stars, so that M_⋆ generally decreases with Z. The fraction of Be-XRBs with low-mass (⁠|$M_{\star }\lt 10\ \rm M_{\odot }$|⁠) O/Be stars (dashed curve in Fig. 10) decreases from ∼60 per cent at Z = 10⁻⁴ to ∼1.5 per cent at Z = 0.001, and then increases quasi-monotonically with Z, reaching ∼87 per cent at Z = 0.03. This is consistent with the evolution of median O/Be star mass with Z and can also be explained by the above arguments about secondary mass, mass transfer rate and stellar winds.

We notice in Fig. 9 that there are two groups of Be-XRBs in the log M_X − log M₁ space with lower and higher M_X (approximately divided by the dash-dotted lines in the bottom panels of Fig. 9) that are more distinct at higher Z. These two groups are also related to the complex features of Be-XRB mass distributions in the log M_X − log M_⋆ space. The lower-M_X group is made of binaries with relatively low-mass (⁠|$M_{1}\lesssim 25\ \rm M_{\odot }$|⁠, i.e. log (M₁ [M_⊙]) ≲ 1.4) primary stars that always transfer significant mass to the secondary star. These binaries also cluster around a primary mass that increases with Z. The higher-M_X group contains both low-mass and massive primary stars. The resulting O/Be star masses cover a broader range than those from the lower-M_X group. We expect the primary stars in the lower-M_X group to completely lose their hydrogen envelope²³ and even undergo mass transfer during the helium HG phase. In this way, most of them form low-mass NSs (⁠|$M_{\rm X}\sim 1.3\ \rm M_{\odot }$|⁠) via electron-capture SNe with no natal kicks. We find that such systems account for about half of the Be-XRBs currently in the SMC (with each simulated Be-XRB re-weighted according to the star formation history of the SMC, see Appendix  A), which is consistent with the results of Vinciguerra et al. (2020). In contrast, the primary stars in the higher-M_X group keep a fraction of their hydrogen envelopes before collapse either due to higher initial masses, shorter lifetimes and/or less mass transfer in wider orbits. There is evidence of this scenario from observations of partially-stripped star + Be star binaries such as HR 6819 (Frost et al. 2022) and SMCSGS-FS 69 (Ramachandran et al. 2023). In general, both groups produce O/Be stars of a broad mass range, and the mass distribution of O/Be stars becomes bi-polar at Z ≳ 0.001, which is more obvious for O/Be stars from the higher-M_X group. This likely results from the complex dependence of the mass loss/accretion rates on primary/secondary masses and orbital parameters in BSE models, which also correlate with the natal kicks and remnant masses of primary stars, as hinted by observations (Zhao et al. 2023). We defer a thorough discussion on this aspect to future work.

Fig. 11 shows the orbital parameter distribution of Be-XRBs. Similarly to the results in Section 5.1, nearly-circular (e ≲ 0.1) binaries with |$a\sim 100-10^{3}\ \rm R_{\odot }$| make up a significant fraction of Be-XRBs, reaching |$\sim 60{{\ \rm per\ cent}}$| at Z = 10⁻⁴. This is a natural consequence of binary interactions that tend to circularize binary orbits and our optimistic definition of the VDD boundary (equation 8) that allows faint objects with little accretion from beyond the tidal truncation radius to be counted as Be-XRBs (Section 3.1). These systems mostly belong to the aforementioned lower-M_X group and contain low-mass (⁠|$M_{\rm X}\sim 1.3\ \rm M_{\odot }$|⁠) NSs born in electron-capture SNe of helium stars with no natal kicks. Since strong mass transfer happens in their evolution histories, the O/Be stars are initially less massive and live longer than in the case of e > 0.1. We also find that it is necessary to take into account such nearly-circular systems in order to explain the large population of Be-XRBs currently observed in the SMC, because the number of nearly-circular Be-XRBs as a function of time after a starburst has a strong peak at |$\sim 30\ \rm Myr$| for Z = Z_SMC = 0.0035, and the SMC experienced a starburst just ∼20–40 Myr ago (Rubele et al. 2015). This result is consistent with the finding in Linden et al. (2009) that Be-XRBs currently in the SMC primarily form through electron-capture SNe with low natal kicks.

If we ignore the nearly circular binaries, our results are consistent with those in Xing & Li (2021, see their fig. 6) who define the VDD boundary with the tidal truncation radius (equation 7) such that nearly-circular systems are excluded. That is to say, wide (⁠|$a\gtrsim 300\ \rm R_{\odot }$|⁠) binaries with longer a need to have higher e to interact with the VDD at the pericenter. On the other hand, there is an upper limit of e that increases with a in close binaries (⁠|$a\lesssim 300\ \rm R_{\odot }$|⁠) to avoid RLO of the O/B star. Finally, we find that the distribution of a is broader at higher Z, while the median is almost constant at |$a\sim 300\ \rm R_{\odot }$|⁠. The increase of the fraction of wide Be-XRBs at higher Z can be explained by the stronger winds that widen binary orbits more significantly. The increasing relative importance of very close binaries (⁠|$a\lesssim 100\ \rm R_{\odot }$|⁠) at higher Z corresponds to the decreasing importance of the lower-M_X group that produces a stronger peak around |$a\sim 500\ \rm R_{\odot }$| at lower Z, and the trend will disappear if we exclude nearly-circular binaries.

5.3 X-ray outputs

To characterize the X-ray outputs from the simulated Be-XRB populations, we start with the time evolution of total X-ray luminosity per unit stellar mass from an instantaneous starburst (at t = 0):

where Be-XRB i lives from |$t_{i,\rm ini}$| to |$t_{i,\rm fin}$| with a (specific) luminosity L_{ν, i} during outbursts for a certain energy (band)²⁴ and a duty cycle f_{duty, i}, and Θ is the Heaviside step function. Fig. 12 shows the results for the 0.5 − 8 keV band from the SR_CS model with f_corr = 0.5 at 4 representative metallicities. There is a delay of ∼3 − 10 Myr between the starburst and the onset of X-ray emission from Be-XRBs, which reflects the evolutionary time for (the most) massive primary stars in Be-XRB progenitors to become compact objects. The total X-ray luminosity peaks at a few Myr after the starburst for Z ≲ 0.01, while the peak (as well as the onset) of X-ray emission shifts to later stages at higher Z, up to ∼20 Myr for Z = 0.03. This is caused by the suppression of Be-XRBs from massive binaries by stellar winds that is more significant at higher Z (see Section 5.2) and longer lifetimes of (initially less massive) primary stars in Be-XRB progenitors with higher Z. The X-ray luminosity drops by at least 2 orders of magnitude within 100 Myr post the peak. The most metal-poor model with Z = 10⁻⁴ has a much slower drop compared with the other cases due to a higher fraction of Be-XRBs with low-mass Be stars from low-mass progenitors (see Figs 9 and 10).

Now, given the short-lived nature of Be-XRBs, their overall X-ray output can be well characterized by the (specific) X-ray luminosity per unit SFR:

where |$\tau _{i}=t_{i,\rm fin}- t_{i,\rm ini}$| is the duration of the Be-XRB phase for binary i. Fig. 13 shows the full (intrinsic) SEDs in terms of |$\nu \mathcal {L}_{\nu }$| for Z = 10⁻⁴ and 0.0142, in the SR_CS model with f_corr = 0.5, where we also plot the contributions of different types of outbursts (Type I and II, see Section 3.2), accretion regimes (LH, HS, and SE, see Section 4) and compact objects (NS and BH). In general, for the photon energy range E ∼ 0.5–8 keV that typically contains |$\sim 60{{\ \rm per\ cent}}$| of the integrated bolometric luminosity from a Be-XRB population, the SED is always dominated by NSs in the SE regime with η > 2 mostly via Type II outbursts, and the contribution of BHs remains below 10 per cent. This means that the spectral shape is almost independent of Z in this energy range and the majority of X-ray emission comes from luminous (⁠|$\gtrsim 10^{38}\ \rm erg\ s^{-1}$|⁠) systems. In fact, SE systems contribute ∼55–68 per cent of the total bolometric luminosity per unit SFR for Z ∈ [10⁻⁴, 0.03] with mildly increasing importance at lower Z. Similarly, Type II outbursts make up ∼56–71 per cent of the total X-ray output due to their luminous nature, even though the number fraction of all (active) Type II systems is (much) lower, ∼13–44 (2–16) per cent. For more energetic photons (E ≳ 8 keV), the contribution from NSs in the HS regime (η ∼ 0.05–2) becomes important due to their harder spectra, especially at higher Z, while for E ≲ 0.1 keV, BHs (mostly in the SE regime during Type II outbursts) dominate the spectrum when they exist in Be-XRBs at Z ≲ 0.001.

Next, to better evaluate the metallicity dependence, we calculate the X-ray luminosity per unit SFR in select energy bands for each simulated Be-XRB population. We start with two energy bands with very-soft (0.1–2 keV) and soft (0.5–2 keV) X-rays that are expected to escape the host haloes and interact with the IGM at z ≳ 20 and z ∼ 10, respectively (Das et al. 2017; Sartorio et al. 2023). We further consider a hard band (2–10 keV) and a broad band (0.5–8 keV) to make comparisons with literature results (see below). Fig. 14 shows the X-ray luminosity per unit SFR as a function of Z in these four bands for the SR_CS model with f_corr = 0.5. The results for all four models in Table 2 are summarized in Table 3.

5.3.1 Be-XRBs VS other types of HMXBs

We first compare our predictions for Be-XRBs with the BPS results for other types of HMXBs powered by RLO and (spherical) stellar winds (Fragos et al. 2013b). As shown in Fig. 14, we find that (at the same Z) the X-ray luminosity (per unit SFR) from Be-XRBs in our SR_CS model with f_corr = 0.5 is comparable (lower by up to a factor of 3) to that from other types of HMXBs predicted by Fragos et al. (2013b) for Z ∼ 0.001–0.02, where the evolution with Z is also similar. This indicates that Be-XRBs can be as important as other types of HMXBs. Moreover, the striking similarity between the metallicity dependence of X-ray outputs in the two populations of HMXBs powered by distinct mechanisms (decretion discs of O/Be stars VS RLO and stellar winds) can be understood by the fact that their key properties (number counts and accretion rates of compact objects) are determined by the same binary stellar evolution processes (e.g. stellar winds, mass transfer, and SN natal kicks) in this metallicity range.

However, at Z ≲ 0.001 and Z ≳ 0.02, our results for Be-XRBs show stronger evolution with Z. In Fragos et al. (2013b), the X-ray luminosity from other types of HMXBs is almost constant for Z ≲ 0.001, while in our case, the X-ray luminosity from Be-XRBs increases by a factor of ∼2 from Z = 0.001 to Z = 10⁻⁴. The rapid drop of X-ray luminosity at Z > 0.02 in our case is caused by the stronger stellar winds at higher Z that significantly reduce the number of (luminous) Be-XRBs (see Figs 7 and 8), such that the predicted X-ray output has large errors (∼0.5 dex) due to the poor statistics of luminous Be-XRBs. In fact, the X-ray output of the Z = 0.03 model can be enhanced by a factor of ∼10 when using the wind model in Hurley, Tout & Pols (2002) that predicts weaker winds at high Z than considered in our default wind prescription based on Schneider et al. (2018) and Sander & Vink (2020). However, the effect of varying wind models on the total X-ray output is much weaker at lower metallicities (within factors of ∼3 and ∼20 per cent for Z ≲ 0.02 and Z ≲ 0.005, respectively).

The difference between our results and those in Fragos et al. (2013b) only varies slightly with the energy band considered, with the broad band (0.5–8 keV) having the smallest difference and the hard band (2–10 keV) the largest. This indicates that the overall SEDs of Be-XRB populations in our case are similar to those in Fragos et al. (2013b) for other types of HMXBs. For instance, we have L_{[0.5-2] keV}/L_[2-10] keV ≃ 0.41, and this ratio is ≃0.35 in Fragos et al. (2013b). The small difference is caused by the different X-ray spectral models adopted in our work (Fig. 13) and by Fragos et al. (2013b, see their fig. 1). The latter predict slightly harder spectra with stronger X-ray emission at E ≳ 2 keV. The reason is that their spectral model (Fragos et al. 2013a) only considers the LH and HS regimes and is calibrated to the BC factors derived from the early samples of Galactic XRBs in McClintock & Remillard (2006) and Wu et al. (2010), while our model further includes the SE regime (which is dominant for Be-XRBs) and adopts the recent observational data for BC factors in Anastasopoulou et al. (2022).

5.3.2 Comparison with HMXBs observed in nearby galaxies

Next, we compare the overall X-ray outputs from our Be-XRB populations with observational results for HMXB populations in nearby galaxies (Douna et al. 2015; Lehmer et al. 2021). Assuming that observations are complete, we find that the simulated Be-XRBs in our fiducial model (SR_CS with f_corr = 0.5) can explain |$\sim 60{{\ \rm per\ cent}}$| of the observed X-ray luminosity of (all types of) HMXBs in nearby galaxies (Lehmer et al. 2021), with a similar metallicity dependence (see the shaded region in Fig. 14). The X-ray luminosity increases by about one order of magnitude from Z ∼ 0.02 to Z ≲ 0.0003. This is consistent with the results in Section 5.1 for number counts of ultra-luminous HMXBs with |$L_{[0.5-8]~\rm keV}\gt 10^{39}\ \rm erg\ s^{-1}$| (Fig. 8). It is shown in Appendix  B that both the observed number of (ultra-)luminous HMXBs and overall X-ray luminosity per unit SFR from Lehmer et al. (2021, for all types of HMXBs) can be fully reproduced by our SR_CS model of Be-XRBs at Z ∼ 0.0003–0.02, if given a higher correction factor f_corr = 0.8 than that assumed by default (f_corr = 0.5). This is non-trivial because systems with |$L_{[0.5-8]~\rm keV}\gt 10^{39}\ \rm erg\ s^{-1}$| only account for ∼15–50 per cent of the total X-ray output in our Be-XRB populations at Z ∼ 0.0003–0.02 from the SR_CS model with f_corr ∼ 0.25–0.8. Moreover, this indicates that f_corr ≳ 0.8 is ruled out in our SR_CS model, considering that other types of HMXBs also contribute to the observed values.

Finally, we calculate the X-ray luminosity per luminous Be-XRB with |$L_{[0.5-8]~\rm keV}\gt 10^{38}\ \rm erg\ s^{-1}$| in the 0.5–8 keV band |$\langle L_{\rm [0.5-8]\ keV}^{38}\rangle /\langle N_{\rm X}^{38}\rangle$|⁠, as shown in Fig. 15 for all four models in Table 2 assuming f_corr = 0.5. In general, we have |$\langle L_{\rm [0.5-8]\ keV}^{38}\rangle /\langle N_{\rm X}^{38}\rangle \sim 0.5-1\times 10^{39}\ \rm erg\ s^{-1}$|⁠, insensitive to both metallicity and model parameters even with variations of f_corr in the range of [0.25,1] (not shown for conciseness). Our results are consistent with observations of HMXBs in nearby galaxies for Z ∼ 0.0004–0.03 (see fig. 5 of Douna et al. 2015), which further implies that the luminous X-ray sources in our Be-XRB populations are similar to those in observed HMXB populations.

6 CAVEATS

To quantify more accurately the roles played by Be-XRBs, particularly in the ULX population (e.g. Kaaret, Feng & Roberts 2017; Kovlakas et al. 2020; Fabrika et al. 2021; Walton et al. 2022; King, Lasota & Middleton 2023; Salvaggio et al. 2023; Tranin, Webb & Godet 2023; Misra et al. 2023a), and their imprints in the early Universe, more work needs to be done to overcome the following caveats of our study (with descending order of importance).

• Although being more complete compared with previous studies, our model of the X-ray outbursts in Be-XRBs still relies on a phenomenological model based on simulations of steady-state VDDs with constant mass ejection rates (equation 19) and viscosity. However, in reality VDDs of O/Be stars, especially under the influence of companions, can be highly dynamic structures (Carciofi & Bjorkman 2008; Haubois et al. 2012; Krtička, Kurfürst & Krtičková 2015; Panoglou et al. 2016; Vieira et al. 2017; Rímulo et al. 2018), with variable mass ejection and disc dissipation/formation episodes of a few years (Reig 2011). Such dynamical evolution of disc structures are particularly important for the most luminous Be-XRBs with non-periodic, strong Type II outbursts (see Section 3.2.2). In fact, according to our steady-state model, the majority (⁠|$\sim 60-70{{\ \rm per\ cent}}$|⁠) of X-ray emission from Be-XRBs are produced by the systems in which the compact object is able to accrete almost all (≳90 per cent) of the materials ejected from the O/B star, so that significant variations of disc structures are expected to occur. Additional uncertainties of our model may reside in the interpretation of simulation results and comparison with observations, as we do not model the detailed structure of accretion flows and geometry of X-ray emission, which are complex and still in debate, especially for ULXs (Wiktorowicz et al. 2019; Fabrika et al. 2021; Mushtukov & Tsygankov 2022; King, Lasota & Middleton 2023; Lasota & King 2023). For instance, taking into account the suppression of accretion by radiation-driven winds in the super-Eddington regime (Shakura & Sunyaev 1973) reduces the X-ray outputs and number counts of ULXs from our Be-XRB populations by up to ∼60 per cent and a factor of ∼10, respectively.

• We do not fully explore the parameter space of (binary) stellar evolution. Although the mass-transfer efficiency in our BPS runs is calibrated to reproduce the number and range of orbital periods of Be-XRBs observed in the SMC, the agreement with observations is imperfect, and the other parameters (fixed in our case) governing the mass loss rate from stellar winds, stability of mass transfer, angular momentum loss, spin-up of the accretor, as well as natal kicks of SNe may affect Be-XRB properties and final X-ray outputs greatly (see, e.g. Linden et al. 2009; Zuo, Li & Gu 2014; Shao & Li 2014, 2020; Vinciguerra et al. 2020; Xing & Li 2021; Misra et al. 2023b; Willcox et al. 2023). For instance, recent parallax and proper-motion measurements of young, isolated radio pulsars and Galactic Be-XRBs have placed new constraints on SN natal kicks. These support a kick-velocity distribution that is the sum of two Maxwellian distributions with σ_kick ∼ 50 and |$340\ \rm km\ s^{-1}$|⁠, where the low-σ_kick component accounts for 20 ± 10 per cent of all systems (Verbunt, Igoshev & Cator 2017; Igoshev et al. 2021). Our model of natal kicks is qualitatively consistent with these results given the separate treatments of electron-capture SNe (with σ_kick = 0) and high-σ_kick cases (Type II and Ib/c with σ_kick = 190), although our σ_kick are smaller. We find by numerical experiments that compared with our default case, using the best-fit model of kick-velocity distribution in Igoshev et al. (2021) with generally stronger kicks will reduce the X-ray output from Be-XRBs by up to a factor of ∼3 for Z ≲ 0.02. However, with the kick-velocity distribution from Igoshev et al. (2021), Be-XRBs on nearly circular (e ≲ 0.1) orbits become extremely rare (≲ 4 per cent), and the predicted number of Be-XRBs in the SMC is lower by a factor of ∼2 than that observed (see Appendix  A). We defer a more detailed investigation into the effects of SN natal kicks on the X-ray output from Be-XRBs to future work. The prescription of stellar winds is also important, especially for metal-rich stellar populations, causing up to one order of magnitude variations in the X-ray output at Z = 0.03 (see Section 5.3.1).

• In the construction of our binary populations, we do not consider the correlations between initial binary properties (i.e. binary fraction, distributions of masses, and orbital parameters) nor their metallicity dependence, while there is evidence/hints of such correlations and evolution with metallicity from observations²⁵ and simulations (e.g. Gunawardhana et al. 2011; Marks et al. 2012; Duchêne & Kraus 2013; Moe & Di Stefano 2017; Moe, Kratter & Badenes 2019; Jeřábková et al. 2018; Lacchin et al. 2020; Chon, Omukai & Schneider 2021; Chon et al. 2022; Tanvir & Krumholz 2023; Rusakov, Steinhardt & Sneppen 2023). Taking these effects into account may modify the predicted evolution of X-ray luminosity with metallicity and change our finding that the total X-ray output from Be-XRBs is insensitive to initial conditions.

• The effects of rotation on stellar evolution are not considered in our BPS runs, although it is shown by detailed stellar-evolution simulations (e.g. Ekström et al. 2012; Georgy et al. 2013; Choi, Conroy & Byler 2017; Groh et al. 2019; Murphy et al. 2021) that fast rotation can impact the mass loss, time-scales of evolution phases, stellar structure, nucleosynthesis and remnant masses in a complex and metallicity-dependent manner, particularly for initially FR stars, which can also change the properties of Be-XRBs that by definition contain FR O/B stars.

• We use simple models calibrated to observations to derive the X-ray spectra of Be-XRBs, which only consider three regimes of accretion rates (low-hard, high-soft, and super-Eddington) for two types of compact objects (NSs and BHs), while in reality, the X-ray spectra of Be-XRBs can have more complex dependence on compact object properties (e.g. masses, spins and magnetic fields) and accretion modes (e.g. thin disc and ADAF) that may also correlate with the stellar and orbital parameters as well as types of outbursts (Martin et al. 2011; Okazaki, Hayasaki & Moritani 2013; Cheng, Shao & Li 2014; Haberl & Sturm 2016; Xu & Li 2019; Franchini & Martin 2021).

7 SUMMARY AND CONCLUSIONS

We improve population synthesis of Be-XRBs with a physically motivated model that combines recent hydrodynamic simulations of VDDs in Be-XRBs (Brown et al. 2019) and VDD properties inferred from observations of classical Be stars (Vieira et al. 2017; Rímulo et al. 2018). Our model for the first time fully takes into account the dependence of X-ray outburst properties (i.e. strength and duty cycle) on stellar and orbital parameters of Be-XRBs, and also considers the classification of X-ray outbursts into the two conventional types in observations (Reig 2011; Rivinius, Carciofi & Martayan 2013). Using the standard BSE models (Hurley, Tout & Pols 2002) in the BPS code binary_c, we evolve large populations of randomly sampled massive binaries in the (absolute) metallicity range Z ∈ [10⁻⁴, 0.03] with initial conditions and BSE parameters calibrated to reproduce the population of observed Be-XRBs in the SMC (see Appendix  A). Finally, we apply an X-ray spectral model based on recent observations of HMXBs (Anastasopoulou et al. 2022) to our simulated Be-XRBs to evaluate the X-ray luminosity per unit SFR from Be-XRBs as a function of metallicity. The effects of free model parameters for initial conditions and VDD densities are also explored. Our main findings are summarized as follows.

• Be-XRBs can probably explain a non-negligible fraction of

  the total X-ray output from HMXBs observed in nearby galaxies (Lehmer et al. 2021). Assuming that current observations are complete, the X-ray luminosities per unit SFR from our simulated Be-XRB populations are at least |$\sim 30{{\ \rm per\ cent}}$| of the observed values for all types of HMXBs at Z ∼ 0.0003–0.02. The minimum fraction |$\sim 30{{\ \rm per\ cent}}$| is reached with the minimum correction factor f_corr = 0.25 allowed in our Be-XRB model calibrated to the observed relation between outburst X-ray luminosity and orbital period (Section 3.2.1), and is subject to uncertainties from our imperfect modelling of VDDs, the limited parameter space of BPS explored, and the fact that our BPS simulations only consider Be-XRBs rather than all types of HMXBs (see Section 6). In our fiducial model with no Z dependence of VDD densities nor enhanced initial stellar rotation, the X-ray luminosity per unit SFR decreases by a factor of ∼8 from Z = 0.0003 to Z = 0.02, which is mainly driven by the stronger stellar winds at higher Z that reduce the number of binary stars able to become luminous Be-XRBs. Interestingly, this trend is similar to that seen in observations of all types of HMXBs.

• Our results for Be-XRBs are also similar to the Z-dependent X-ray outputs from other types of HMXBs powered by RLO and (spherical) stellar winds in previous BPS studies (e.g. Fragos et al. 2013b) for Z ∼ 0.001 − 0.02 (Section 5.3.1). However, at Z ≲ 0.001, the X-ray luminosity from Be-XRBs keeps increasing towards lower Z (e.g. by a factor of ∼2 from Z = 0.001 to Z = 10⁻⁴), while there is almost no evolution in that from other types of HMXBs predicted by Fragos et al. (2013b). This implies that Be-XRBs might play more important roles at lower metallicities, such that they can have strong effects on the thermal history and thus 21-cm signal from the IGM at Cosmic Dawn. More comprehensive BPS studies that consider both Be-XRBs and other types of HMXBs self-consistently are required to verify if this is indeed the case.

• Luminous (⁠|$\gtrsim 10^{38}\ \rm erg\ s^{-1}$|⁠) systems with non-periodic, major outbursts of Type II (see Section 3.2.2 for a detailed definition) and super-Eddington accretion (mostly on eccentric orbits with e ≳ 0.4), including ULXs (⁠|$\gtrsim 10^{39}\ \rm erg\ s^{-1}$|⁠), are important X-ray powerhouses in our Be-XRB populations across the metallicity range considered. In our fiducial model, |$\sim 55-75{{\ \rm per\ cent}}$| of the total X-ray output of Be-XRBs comes from luminous sources with |$L_{[0.5-8]\ \rm keV}\gt 10^{38}\ \rm erg\ s^{-1}$|⁠, while ULXs with |$L_{[0.5-8]\ \rm keV}\gt 10^{39}\ \rm erg\ s^{-1}$| contribute |$\sim 24-41{{\ \rm per\ cent}}$| of the total output. Similarly to the case of total X-ray output, the simulated Be-XRBs can also account for non-negligible fractions of the number counts of luminous and ultra-luminous HMXBs in observations (Mapelli et al. 2010; Douna et al. 2015; Kovlakas et al. 2020; Lehmer et al. 2021), although more work needs to be done to fully evaluate the contributions of Be-XRBs. Such luminous systems have been overlooked in previous BPS studies of Be-XRBs (e.g. Zuo, Li & Gu 2014; Misra et al. 2023b), which use a simple empirical scaling law between outburst X-ray luminosity and orbital period (see our equation 12 and fig. 3 of Dai, Liu & Li 2006) that imposes a (likely artificial) cutoff at |$\sim 10^{38}\ \rm erg\ s^{-1}$|⁠, while our improved Be-XRB model further considers the dependence of outburst accretion rate/luminosity on VDD density (and its scatter) as well as orbital eccentricity, so that the luminosity function of our Be-XRBs has a luminous tail up to |$10^{41}\ \rm erg\ s^{-1}$|⁠, as shown in Fig. 16.

• The major uncertainty in the total X-ray output from Be-XRBs arises from the uncertainties in accretion rates and VDD properties rather than those in initial conditions when fixing the IMF (see Section 2.3). The former can change the final X-ray output by a factor of a few, while the difference made by the latter is within ∼50 per cent (see Section 5.1 and Appendix  B). The reason is that the X-ray luminosity of an accreting compact object is proportional to the accretion rate, which is then proportional to the VDD density (or mass ejection rate of the O/B star) to the first order, and there are still large uncertainties in these quantities from idealized simulations and limited observations, especially at low metallicity. On the other hand, a key feature of our BPS simulations is that in most binaries that can potentially become Be-XRBs the secondary star will be spun up to become an O/Be star (via stable mass transfer) regardless of its initial rotation rate, which is consistent with the scenario that most (young) O/Be stars are produced by binary interactions as predicted/hinted by theoretical models and observations (Shao & Li 2014; Klement et al. 2019; Bodensteiner, Shenar & Sana 2020; Dorigo Jones et al. 2020; Hastings, Wang & Langer 2020; Hastings et al. 2021; Dallas, Oey & Castro 2022; Dodd et al. 2023; Wang et al. 2023).

In general, our BPS results highlight the possibility that Be-XRBs constitute an important component in the population of HMXBs, making significant contributions to the overall X-ray output and number count of ULXs.

Indeed, robust modelling of XRBs with BPS is a challenging task (see Section 6) because of large uncertainties in binary stellar evolution, particularly for Be-XRBs due to their transient nature and our lack of understanding on the detailed mechanisms that drive the formation and evolution of VDDs around FR massive stars (for candidate scenarios, see, e.g. Granada et al. 2013; Hastings, Wang & Langer 2020; Zhao & Fuller 2020; Cranmer 2009; Rogers et al. 2013; Lee, Neiner & Mathis 2014; Ressler 2021). Nevertheless, such theoretical efforts are worthwhile considering the important roles played by XRBs in the Epoch of Reionization and Cosmic Dawn (e.g. Fialkov, Barkana & Visbal 2014; Fialkov & Barkana 2014; Pacucci et al. 2014), as well as the wealth of high-z observational data to come in the next decades for the 21-cm signal by radio telescopes such as HERA (DeBoer et al. 2017), SARAS (Singh et al. 2022), REACH (de Lera Acedo et al. 2022), NenuFAR (Mertens, Semelin & Koopmans 2021), MIST (Monsalve et al. 2023), and SKA (Koopmans et al. 2015) which will place constraints on high-z XRBs and the underlying star/galaxy/structure formation processes.

The present paper is a small step that points out the importance of Be-XRBs. In future work, we plan to model Be-XRBs (with possible improvements of the aforementioned caveats) and other types of XRBs in one BPS framework to predict their impact on the 21-cm signal as well as cosmic X-ray background and derive constraints from existing observational data (e.g. Lehmer et al. 2012; Moretti et al. 2012; Cappelluti et al. 2017; Fialkov et al. 2017; Bowman et al. 2018; Bevins et al. 2023; Lazare, Sarkar & Kovetz 2023; Rossland et al. 2023). We will also extend the metallicity range down to the extremely metal-poor (Z ≲ 10⁻⁶) regime of Pop III stars,²⁶ which are expected to be promising progenitors of Be-XRBs (Sartorio et al. 2023) given their compact and fast-rotating²⁷ nature that favours stable mass transfer (Inayoshi et al. 2017) and formation of O/Be stars. Another interesting topic to investigate is the connection between (Be-)XRBs and compact object mergers (see e.g. Marchant et al. 2017; Mondal et al. 2020; Fishbach & Kalogera 2022; Kotko & Belczynski 2023; Liotine et al. 2023) that may foster the synergy between 21-cm and gravitational wave observations to better constrain binary stellar evolution in high-z galaxies.²⁸

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for the insightful comments and David D. Hendriks for his help with binary_c-python. BL and AF are supported by the Royal Society University Research Fellowship. NS gratefully acknowledges the support of the Research Foundation - Flanders (FWO Vlaanderen) grant 1290123N. RGI acknowledges funding by the Science and Technology Facilities Council (STFC) consolidated grants ST/L003910/1 and ST/R000603/1, and thanks the Brexit Research and Interchange on Differentiated Governance in Europe (BRIDGE) network. This work excessively used the public packages numpy (van der Walt, Colbert & Varoquaux 2011), matplotlib (Hunter 2007), and scipy (Virtanen et al. 2020). The authors wish to express their gratitude to the developers of these packages and to those who maintain them.

DATA AVAILABILITY

The data of Be-XRBs and python scripts used to manage BPS runs and conduct post-processing will be shared on reasonable request to the corresponding authors. The BPS code binary_c (including the Be-XRB module) is available at https://gitlab.com/binary_c/binary_c/-/tree/development-BeXRB. The corresponding python interface binary_c-python is available at https://gitlab.com/binary_c/binary_c-python/-/tree/feature/source_file_sampling.

Footnotes

References

APPENDIX A: BE-XRBS IN THE SMC

As mentioned in Section 2, the BSE parameters and initial conditions of our BPS runs are chosen to reproduce the observed population of Be-XRBs in the SMC (Coe & Kirk 2015) at the metallicity Z_SMC = 0.0035 (Davies et al. 2015). The agreements with observations are evaluated in two aspects: (1) The total number N_BeXRB of binaries currently in the Be-XRB phase with orbital periods P_orb < 10³ d²⁹ is larger than 69, given the 69 Be-XRBs already found in the SMC (Coe & Kirk 2015). (2) The orbital period distribution of these binaries should cover P_orb ∼ 10–10³ d with a peak around 100 days (Coe & Kirk 2015; Vinciguerra et al. 2020).

The BPS results (i.e. N Be-XRBs from a single-age stellar population of a total mass M_tot) are convolved with the star formation history of the SMC (see their fig. 16 Rubele et al. 2015) to calculate N_BeXRB as

where SFR(t) is the star formation history as a function of look-back time t, and |$t_{i,\rm ini\ (fin)}$| marks the beginning (end) of the Be-XRB phase of binary i. Similarly, the orbital period distribution is derived by applying the weight w_i in equation (A1) to each Be-XRB.

The results for the SR_CS model (see Table 2) with the default initial binary statistics are shown in Fig. A1, which are almost the same as those for the SR_OP model because our Be-XRB routine (Section 3) does not affect binary stellar evolution. The number count requirement is satisfied with |$N_{\rm BeXRB}= 74_{-11}^{+13}$|⁠, very close to the observed number 69, implying that the observed sample of Be-XRBs in the SMC is nearly complete. Here, the uncertainties in N_BeXRB arise from the errors in SFR(t). The observed orbital period distribution is well reproduced by the SR_CS model at P_orb < 10³ days with a Wasserstein distance of 0.11. Our results are also consistent with those from the BPS study by Vinciguerra et al. (2020, see their fig. 3 and table 2).³⁰ Nevertheless, the predicted distribution is broader and peaks at ∼300 d compared with the observed distribution with a narrower peak around 100 d. The overprediction of long-period (P_orb ≳ 200 d) systems may be explained by incompleteness in observations but the discrepancy at P_orb ∼ 10–100 d is difficult to reconcile, which indicates that our binary stellar evolution and Be-XRB models are still imperfect.

For comparison, we show the results for the FR_CS model in Fig. A2 (with the default initial conditions), which are also almost the same as those for the FR_OP model. Here, we have |$N_{\rm BeXRB}= 98_{-15}^{+18}$|⁠, allowing more space for the potential incompleteness of observations. With optimistic SFRs of the SMC, the predicted number of Be-XRBs is larger than the observed number in almost every bin of P_orb. However, the discrepancy in the shape of the orbital period distribution is larger with a Wasserstein distance of 0.16. In particular, the number of low-period (P_orb ∼ 10 d) systems is overpredicted by up to a factor of ∼10.

Finally, we find that the orbital period distribution of Be-XRBs in the SMC is sensitive to the initial binary properties. When using the hybrid initial orbital period distribution (equations 2–4) rather than the default log-flat initial separation distribution, for the SR_CS model we obtain |$N_{\rm BeXRB}= 65_{-10}^{+12}$|⁠, lower than the number |$74_{-11}^{+13}$| in the default case, but still consistent with observations. The predicted P_orb distribution now also peaks around 100 days as in observations, achieving a slightly smaller Wasserstein distance of 0.1, as shown in Fig. A3. The results for the FR models are similar to those under the default initial conditions (Fig. A2) and not shown.

APPENDIX B: ADDITIONAL RESULTS FROM ALTERNATIVE MODELS

To illustrate the weak dependence of our results on the distribution of initial binary orbital parameters, Fig. B1 shows the relative difference between the X-ray outputs in the 0.5–8 keV band from the two initial condition models considered in our study, i.e. the default model with a log-flat distribution of initial separations and the alternative model with an observation-based hybrid initial orbital period distribution (Section 2.3). In general, the relative difference is small, i.e. within |$20{{\ \rm per\ cent}}$| for Z ≲ 0.02 and up to |$\sim 50{{\ \rm per\ cent}}$| for the Z = 0.03 case (with poor statistics). Therefore, in the rest of this appendix, we only show the results from the default initial condition model.

To understand the impact of model parameters on the total X-ray output from Be-XRBs, in Fig. B2 we compare the X-ray luminosity per unit SFR in the 0.5–8 keV band for all four models in Table 2 with f_corr = 0.5. We find that models with higher initial rotation rates (FR) and VDD densities (OP) produce slightly (⁠|$\lesssim 50{{\ \rm per\ cent}}$|⁠), and significantly (up to a factor of ∼4) higher X-ray luminosities than in the fiducial case SR_CS, respectively. In the OP models with f_corr = 0.5, the luminosity from Be-XRBs is higher than that observed for all types of HMXBs (Lehmer et al. 2021) by up to a factor of ∼2. This means that the OP models can only be consistent with observations with f_corr ≲ 0.25. Indeed, it is shown in B3 that the X-ray luminosity and number of (ultra-)luminous HMXBs per unit SFR from Be-XRBs in the SR_OP model reach the observed values for all types of HMXBs given f_corr = 0.25. Since f_corr = 0.25 is already the lower limit allowed by the calibration of accretion rate/luminosity for Be-XRBs (Section 3.2.1), and other types of HMXBs also contribute to the observed luminosity/number, we conclude that the OP models are strongly disfavoured. Similarly, we find that f_corr ≳ 0.8 is ruled out for the CS models, as illustrated in Fig. B4 where the observed X-ray luminosity and number of (ultra-)luminous HMXBs per unit SFR for all types of HMXBs from Lehmer et al. (2021) are reproduced solely by Be-XRBs in the SR_CS model with f_corr = 0.8.