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

X-ray binaries (XRBs) are thought to regulate cosmic thermal and ionisation histories during the Epoch of Reionisation and Cosmic Dawn ($z\sim 5-30$). Theoretical predictions of the X-ray emission from XRBs are important for modeling such early cosmic evolution. Nevertheless, the contribution from Be-XRBs, powered by accretion of compact objects from decretion disks 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 disks. 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\in [10^{-4},0.03]$. We find that Be-XRBs can contribute a significant fraction ($\sim 60\%$) of the total X-ray budget from HMXBs observed in nearby galaxies for $Z\sim 0.0003-0.02$. A similar fraction of observed ultra-luminous ($\gtrsim 10^{39}\ \rm erg\ s^{-1}$) X-ray sources can also be explained by Be-XRBs. Moreover, the predicted metallicty dependence in our fiducial model is consistent with observations, showing a factor of $\sim 8$ increase in X-ray luminosity per unit star formation rate from $Z=0.02$ to $Z=0.0003$.

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 et al. 2023), which dominate the cosmic XRB luminosity density at  ≳ 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 et al. 2013;Xu et al. 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  ∼ 2) galaxies (e.g.Antoniou et al. 2010Antoniou et al. , 2019;;Basu-Zych et al. 2013;Prestwich et al. 2013;Douna et al. 2015;Antoniou & Zezas 2016;Brorby et al. 2016;Lehmer et al. 2016Lehmer et al. , 2019Lehmer et al. , 2021Lehmer et al. , 2022;;Aird et al. 2017;Fornasini et al. 2019Fornasini et al. , 2020) ) 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 Xray 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% 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 (≳ 6 M ⊙ ), main sequence (MS) star (reviewed by e.g.Reig 2011;Rivinius et al. 2013;Rivinius 2019).
Here the massive star is typically of spectral type B (and O) with an rotation velocity above ∼ 70 percent of the equatorial Keplerian limit and shows Balmer emission lines, which can be well explained by a viscous decretion disk 2 (VDD) around the star.This VDD is formed by materials ejected from the star due to redistribution of angular momentum caused by fast rotation 3 (Rivinius et al. 2013).In most observed Be-XRBs, neutron stars (NSs) are identified as the compact companion, although there are two candidate binaries, MWC 656 and AS 386, that contain a black hole (BH) and a Be star 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 et al. 2013;Franchini & Martin 2021).
Be-XRBs make up the largest class of HMXBs identified in observations (Fornasini et al. 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 2 The disks are very light structures compared with the stars, with typical masses  d ∼ 10 −11 − 10 −8 M ⊙ (Granada et al. 2013;Klement et al. 2017;Rivinius 2019).Formation of VDDs can be regarded as simply the means of losing angular momentum for a massive star getting close to the critical limit of rotation (Rivinius 2019). 3The detailed mechanisms for the formation of VDDs around O/B stars are still in debate.Possible mechanisms include mechanical mass loss at critical rotation (e.g.Granada et al. 2013;Hastings et al. 2020;Zhao & Fuller 2020), pulsations (e.g.Cranmer 2009;Rogers et al. 2013;Lee et al. 2014) and small-scale magnetic fields (e.g.Ressler 2021).Nevertheless, in all these scenarios, rapid rotation is required (Rivinius et al. 2013).
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 et al. 2014;Misra et al. 2023).However, previous BPS studies of Be-XRBs (e.g.Zhang et al. 2004;Belczynski & Ziolkowski 2009;Linden et al. 2009;Shao & Li 2014, 2020;Zuo et al. 2014;Vinciguerra et al. 2020;Xing & Li 2021;Misra et al. 2023) 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 et al. 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 ionisation 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 et al. 2013;Panoglou et al. 2016;Cyr et al. 2017;Brown et al. 2018Brown et al. , 2019;;Suffak et al. 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 metallicity4 range  ∈ [10 −4 , 0.03] where observational constraints for HMXBs are available.Our model can also be applied to more metal-poor regimes (e.g. ≲ 10 −6 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 (Sec.5.1), mass and orbital parameter distributions (Sec.5.2), and X-ray outputs of Be-XRBs (Sec.5.3), focusing on how they evolve with metallicity.Finally, we summarize our main findings in Section 6, and discuss their caveats and our outlook to future work in Section 7. The key physical quantities used in this paper are summarized in Table 1.total stellar mass underlying the Be-XRB population N X number of Be-XRBs in the outburst phase per unit SFR L  specific X-ray luminosity per unit SFR L X X-ray luminosity per unit SFR for a certain band

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(Izzard et al. , 2006(Izzard et al. , 2009;;Izzard et al. 2017;Izzard & Halabi 2018;Izzard & Jermyn 2023;Mirouh 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 (ZAMS) for 15 Gyr in the (absolute) metallicity range  ∈ [10 −4 , 0.03] with randomly sampled initial binary properties through the python interface binary_c-python (Hendriks & Izzard 2023) of binary_c.
The Be-XRB model will be explained in detail in the next Section 3.An X-ray spectral model based on observations (Sec.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 main sequence (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 et al. 2002) with the default setup of binary_c 5 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 (Sec.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 Sec.2.2 according to the default setup of binary_c detailed in Izzard et al. (2017).It will be shown in Appendix A that our choices of BSE parameters, combined with standard initial conditions of binary stars (Sec.2.3), can reproduce the population of observed Be-XRBs in the SMC at the metallicity  SMC = 0.0035 (Davies et al. 2015).For simplicity, we assume that the BSE parameters and initial conditions do not evolve with metallicity for  ∈ [10 −4 , 0.03].
Although the BSE models used here keep track of the rotation of each star, the effects of rotation on stellar evolution are not considered, which can be significant (particularly for initially fastrotating stars) and complex, covering various aspects (e.g.mass loss, timescales of evolution phases, stellar structure, nucleosynthesis and remnant masses), especially at low , as shown in detailed stellar-evolution simulations (e.g.Ekström et al. 2012;Georgy et al. 2013;Choi et al. 2017;Groh et al. 2019;Murphy et al. 2021).Such effects may also be important in the modeling of Be-XRBs, as fastrotating 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 et al. 2013).

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 6 .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  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  acc,max ∼ ( g,acc / acc ) −1 ∼  acc  acc /(  acc ) ∼  acc / KH,acc , which is limited by the thermal (Kelvin-Helmholtz) timescale  KH ∼   2 acc /( acc  acc ) of the accretor, given  g,acc ∼   acc / acc the specific energy carried by the in-falling matter that needs to be radiated away, the luminosity  acc , mass  acc and radius  acc of the accretor.binary_c adopts 5 binary_c has been updated recently with a new treatment of pairinstability SNe (Farmer et al. 2019), an improved stellar wind prescription (Schneider et al. 2018;Sander & Vink 2020) and stellar evolution of zerometallicity stars based on mesa data (Paxton et al. 2018(Paxton et al. , 2019)).It is also used to study XRBs from zero-metallicity stars (Sartorio et al. 2023), but only considering XRBs powered by RLO and stellar winds (without Be-XRBs). 6Otherwise mass and angular momentum loss during mass transfer shrink binary orbits too efficiently leading to over-prediction of low-period systems.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  KH during mass transfer by expansion and increase of luminosity (Paczyński & Sienkiewicz 1972;Hurley et al. 2002;Vinciguerra et al. 2020).Similarly, we also increase the upper limit on  from the dynamical timescale of the accretor as well as the thermal and dynamical timescales 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 ( ≳ 0.1) whose progenitor primary stars only undergo weak mass loss (i.e. = 1 for  donor ≪  acc,max , independent of  thermal ), as discussed in Sec. 5.

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.Here we briefly describe the default choices for these parameters adopted in our work for binary_c.
The stability of mass transfer is evaluated with the critical mass ratio  crit =  accretor / donor at the onset of mass transfer.We adopt  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 et al. (2002, see their sec. 2.6).Mass transfer beyond what is allowed by the mass-transfer efficiency (Eq. 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 orbital angular momentum of the donor.
The masses of compact object remnants are determined by the CO core masses of progenitors using the original BSE models in Hurley et al. (2000) and Hurley et al. (2002).We apply natal kicks to Type II and Ib/c SNe that follow a Maxwellian distribution with a 1D dispersion of  kick = 190 km s −1 (Hansen & Phinney 1997), while electron-capture SNe have no natal kicks.The latter are typically happen to highly stripped stars in progenitor binaries of Be-XRBs (see Sec. 5.2), which are expected to have weak natal kicks (≲ 30 km s −1 , see the discussion in sec.3.1 of Vinciguerra et al. 2020).Therefore, we use zero natal kicks for simplicity.

Binary sample and initial conditions
To construct the input catalog of binary stars, we sample  B = 3 × 10 5 binaries randomly from widely used distributions of mass and orbital parameters.To be specific, the primary stellar mass  1 is drawn from the Kroupa (2001) initial mass function (IMF) in the mass range of [5, 100] M ⊙ and the mass ratio  ≡  2 / 1 is generated from a uniform distribution in the range  ∈ Here we only consider binaries with  1 ≥ 5 M ⊙ because only massive primary stars can form the compact objects considered in our Be-XRB model7 (see Sec. 3.1).In the way, the binary stars in our catalog only make up a small fraction of the whole underlying stellar population.Following Misra et al. (2023), we assume that the whole stellar population is made of 70% binary stars (Sana et al. 2012) and 30% 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  1 ∈ [0.01, 100] M ⊙ , and the mass ratio distribution for the entire binary star population is uniform in  ∈ [0.01 M ⊙ / 1 , 1].Under these assumptions8 , we estimate that the total mass of stars in our binary sample accounts for  sample = 0.2 of the total mass of the whole stellar population, using the method in appendix A of Bavera et al. (2020).
Here  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 For orbital parameters, by default we follow Izzard et al. (2017) to draw the initial semi-major axis  from a log-flat distribution for  ∈ [3, 10 4 ] R ⊙ and the initial eccentricity  from a thermal distribution for  ∈ [0, 1).We assume no correlations between ,  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  orb (in the unit of day) from a hybrid distribution   orb ≡ / log  orb based on observations of low-mass (Kroupa 1995) and massive (Sana et al. 2012) where m ≡  1 / O , given  O = 16 M ⊙ as the minimum mass of O stars, and (Kroupa 1995;Sana et al. 2012) It turns out that 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  rot,0 , which can be characterized by the parameter  0 ≡  rot,0 / Kep,0 given the initial Keplerian velocity  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 Sec. 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  0 : In the slowly-rotating (SR) model,  0 is a function of initial stellar mass and metallicity (which determines the initial stellar radius given the initial stellar mass) as in Hurley et al. (2000, see their sec.7.2) and Hurley et al. (2002).In the fast-rotating (FR) model, we set  0 = 0.9 for all stars to obtain an upper limit on the formation efficiency as primary and secondary stars in binaries, do not strictly follow the Kroupa (2001) IMF.Nevertheless, for stars above 1 M ⊙ that are relevant for Be-XRBs, the mass distribution is very close to the Kroupa (2001) IMF with small (≲ 20%) deviations above ∼ 60 M ⊙ and a minute Wasserstein distance of 0.04.Therefore, we expect this imperfect sampling of the Kroupa (2001) IMF to have negligible effects on our results.well as 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 fast-rotating (e.g.Ekström et al. 2008;Bastian et al. 2017;Schootemeĳer et al. 2022).

Identification of Be-XRBs
Inspired by previous BPS studies on Be-XRBs (e.g.Zhang et al. 2004;Belczynski & Ziolkowski 2009;Linden et al. 2009;Shao & Li 2014, 2020;Zuo et al. 2014;Vinciguerra et al. 2020;Xing & Li 2021;Misra et al. 2023), we identify a binary as in the Be-XRB phase with the following criteria: (i) The binary is made of a massive MS (O/B) donor star with  ★ > 6 M ⊙ , and a compact object (NS/BH) companion with a mass  X > 1.29 M ⊙ .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 M ⊙ as a conservative estimate of the minimum mass of Be stars in Be-XRBs (Hohle et al. 2010), which is larger than the minimum mass of single Be stars (∼ 3 M ⊙ , 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).
(ii) There is a VDD around the donor star, which we assume to be present when the following conditions are satisfied: (a) The donor star is fast-rotating with  ≡  rot / Kep > 0.7 to eject mass that can potentially settle into a VDD, where  rot is the rotation velocity and  Kep is the Keplerian velocity at the equator of the donor star.Here we adopt the minimum rotation rate  crit = 0.7 for decrection disk formation suggested by Rivinius et al. (2013) based on observations.
(iii) The VDD overfills the Roche lobe of the O/B star at periastron to allow accretion by the compact object from the VDD:  out >  L1 (, ,  ★ ), where (Eggleton 1983) is the Roche lobe radius given the semi-major axis , eccentricity  and mass ratio  ★ ≡  ★ / X , and  out is the effective disk boundary.
In previous studies,  out is either fixed to a typical value (e.g. 100 R ⊙ in Misra et al. 2023), or set to the average tidal truncation radius (Zhang et al. 2004;Xing & Li 2021) given  X ≡  X / ★ , and  trunc = 3 assuming a typical viscosity parameter = 0.63 (Rímulo et al. 2018) (Okazaki et al. 2002;Panoglou et al. 2016), such that accretion is still possible, although weaker, when  trunc <  L1 (, ,  ★ ).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 disk become subsonic (Krtička et al. 2011): according to the equatorial radius  ★ ≡  eq of the O/B star, the sound speed   = √︁ 2 B  d / p in the ionised isothermal VDD of a temperature  d = 0.6 eff (Carciofi & Bjorkman 2006) given the stellar effective temperature  eff , where  B is the Boltzmann constant and  p is proton mass.In this case, the tidal truncation effect is considered in the calculation of peak accretion rate (see Sec. 3.2.1).Given this optimistic definition of VDD boundary, we are able to reproduce the nearly circular ( = 0.06 ± 0.06) Be-XRB, CPD-29 2176, from progenitor binaries of stars in the mass range  1,2 ∼ 10 − 12 M ⊙ with weak SN natal kicks, similar to the progenitors identified in the BPASS models (Richardson et al. 2023, see their table 2).
(iv) The O/B star itself does not fill the Roche lobe at the periastron:  ★ <  L1 (, ,  ★ ).Otherwise, the system will be classified as a RLO XRB (Reig 2011).

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  bol =   acc  2 during outbursts where  is the radiative efficiency and  acc is the peak accretion rate, and (2) the duty cycle  duty , i.e. the effective fraction of time the binary spends in X-ray outbursts during which the average luminosity is  bol .Previous studies (e.g.Zuo et al. 2014;Misra et al. 2023) usually adopt empirical scaling laws or typical values for  bol and  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 10 with a physically motivated X-ray outburst model.To be specific, our model adopts simulation results calibrated to observational data to calculate  bol (Sec.3.2.1),and also considers the classification of X-ray outbursts which, combined with observational constraints, is used to estimate  duty (Sec.3.2.2).

Peak accretion rate & luminosity
We start with the steady-state peak accretion rate (assumed to be identical to the gas capture rate) predicted by the hydrodynamic relative scatter with respect to Eq. 6 following a Gaussian distribution of  = 10% to capture the variations of  (and  trunc ) from system to system. 10 This is an essential step in our Be-XRB modeling since the stellar and orbital properties of Be-XRBs can vary with metallicity and contribute to the metallcicty evolution of the total X-ray output.simulations in Brown et al. (2019), which satisfies a simple powerlaw  acc,sim ∝ [(1 − )] −2  ej as shown in Fig. 1, where  ej ∝ Σ 0 / (Carciofi & Bjorkman 2008) is the steady-state mass ejection rate given the base surface density Σ 0 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  X = 1.4 M ⊙ and a Be star of  ★ = 18 M ⊙ and  ★ = 7 R ⊙ with constant mass ejection rates, covering  ∈ [0, 0.6] and  orb ∼ 40 − 400 days, which will be extrapolated to broader ranges of  and  in our model.For such a Be star, we estimate the stellar luminosity as  ★ ∼ 3.2 × 10 4 L ⊙ and the disk temperature as  d ∼ 2 × 10 4 K, from which we derive the disk surface density as Σ 0,ref ∼ 0.015 g cm −2 for  ej = 10 −10 M ⊙ yr −1 (given the volume density  0 ∼ 5 × 10 −13 g cm −3 of the disk at the stellar surface for  = 0.63), which sets the normalization of Eq. 8.
In reality, mass ejection can be highly variable, even leading to disk dissipation/formation at timescales 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 Eq. 8 should be regarded as an order-of-magnitude estimate that captures the increasing trend with decreasing pericenter distance (1 − ).Finally, we multiply the peak accretion rate from Eq. 8 by a factor of ( trunc / L1 ) 8 for systems with  trunc <  L1 to capture the steepening of disk density profile beyond  trunc (Okazaki et al. 2002).
Next, we associate the VDD base density Σ 0 with the donor star mass  ★ by fitting observational data (Vieira et al. 2017;Rímulo et al. 2018).The obtained empirical scaling laws capture the increasing trend of Σ 0 with  ★ (Arcos et al. 2017; Klement et al. 37 of ≃ 0.52 dex errors (shaded region around the solid line).The squares show the median base disk surface densities in 8 mass bins measured from 54 Be stars in the SMC (Rímulo et al. 2018, R18), where the errorbars show the bin size and 25-75% percentiles of Σ 0 in each bin.These results can also be described with a power-law log 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 disks, such that Eq. 10 should be regarded as an upper limit.Besides, to consider the large scatter in Σ 0 at similar  ★ , for each Be-XRB, we draw a random number  from a Gaussian distribution of a standard deviation  = 0.52 (0.17) dex, and set Σ 0 = 10  Σ0,MW (SMC) ( ★ ), given the prediction of the best-fit model Σ0,MW (SMC) ( ★ ) for the MW (SMC).In addition to the donor mass dependence, we also consider the metallicity  dependence of Σ 0 with two cases.In the conservative (CS) case, we always use the MW model independent of , motivated by the finding that the X-ray luminosity per luminous HMXB is insensitive to metallicity for  ∼ 0.0004 − 0.03 in nearby galaxies (Douna et al. 2015, see their fig.5), while in the optimistic (OP) case, we assume that Σ 0 increases with decreasing metallicity with a linear relation between log(Σ 0 ) and  from solar to SMC metallicities, i.e.  ∼ 0.0035 − 0.0142, and adopt the MW model for  > Z ⊙ = 0.0142 (Asplund et al. 2009) and the SMC model for  <  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  in the range of [10 − 1000] R ⊙ , a uniform distribution of eccentricity for  ∈ [0, 0.6] and a log-flat distribution of From this sample we select a mock population of Be-XRBs with  orb ∼ 10 − 300 days and  X > 10 34 erg s −1 to be compared with observations.These conditions are chosen to mimic the statistics of most (∼ 90%) well observed Be-XRBs in the MW (Raguzova & Popov 2005;Cheng et al. 2014;Brown et al. 2018).The calibration target is the relation between the (outburst) X-ray luminosity  X and orbital periods  orb , derived by Dai et al. (2006, D06) based on 36 observed Be-XRBs from Raguzova & Popov (2005): The best-fit model indicates that the typical X-ray luminosity follows , which is similar to the  dependence in  acc,sim ∝  −2 .In light of this, we assume that  X is proportional to the bolometric luminosity predicted by simulations with a calibration parameter  X :  X =  X   acc,sim  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 realitic, we should have  X ≃ 1.However, we find that the empirical  X - orb scaling relation (Eq.11) 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  corr ∈ [ X , 1] (fixing  X = 0.25) for the peak accretion rate  acc .We also include a power-law term of  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  corr = 1 corresponds to the optimistic case in which the discrepancy is only caused by observational effects, while  corr = 0.25 is the opposite end where the overestimation in simulations needs to be corrected the most.By default, we adopt  corr =  X /  BC,0 = 0.5, assuming a typical bolometric correction (BC) factor  BC,0 =  X / 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% of the bolometric luminosity (see Sec. 4 below).Besides, we do not cape  acc at the Eddington limit considering that Be-XRBs are promising candidates of ultraluminous X-ray sources (ULXs, with  X ≳ 10 39 erg s −1 , reviewed by e.g.Kaaret et al. 2017;Fabrika et al. 2021;King et al. 2023) and a few Be-XRBs with outburst luminosities above the Eddington  (Vinciguerra et al. 2020) to estimate  orb given the spin period  s .The dashed line represents the best fit to observations from Dai et al. (2006) with the shaded region denoting the 1 uncertainties in the fitting parameters (see Eq. 11).The solid line shows the best-fit power-law model for our mock population.It turns out that the mock population reproduces well observational data with a slightly shallower best fit and a similar scatter given  X = 0.25.luminosity  Edd =   Edd  2 ∼ 2 × 10 38 erg s −1 for typical NSs with  = 0.2 and  X = 1.4 M ⊙ 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 3 ) mass transfer and potentially become ULXs (e.g.Marchant et al. 2017;Wiktorowicz et al. 2017Wiktorowicz et al. , 2019Wiktorowicz et al. , 2021;;Shao et al. 2019;Shao & Li 2020;Abdusalam et al. 2020;Kuranov et al. 2020;Misra et al. 2020).It will be discussed below that ULXs are important in our Be-XRB populations.

Classification of X-ray outbursts & duty cycle
Now we can calculate the outburst strength by Eqs. 12 and 13.We further classify the outbursts into two categories following the convention in observations to estimate the duty cycle (Reig 2011;Rivinius et al. 2013): (i) Type I outbursts are regular, (quasi-)periodic and short-lived (∼ 0.1 − 0.3  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  X ≲ 10 37 erg s −1 .The duty cycle is typically  duty,I ∼ 0.1 − 0.3 (Reig 2011;Sidoli & Paizis 2018).
(ii) Type II outbursts are major enhancements of X-ray flux, by a factor of 10 3 − 10 4 , 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 disk 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:  duty,II ∼ 10 −3 −0.1 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  s,crit = 40 s (Cheng et al. 2014;Haberl & Sturm 2016;Xu & Li 2019).To explain this correlation, it is proposed by Okazaki et al. (2013) with considerations of accretion timescale 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 disk and is spun up efficiently to have spin periods  s ≲ 40 s.In Type I outbursts, accretion is in the form of advection dominated accretion flow (ADAF) resulting in low spins with  s ≳ 40 s.It is further shown by Cheng et al. (2014) that disk 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 et al. (2014) for the power-law+Gaussian VDD model (see sec.2.2 of Martin et al. 2011), the tidal warping criterion is satisfied when the disk truncation radius at periastron  trunc,p =  trunc (1 − ) is larger than the tidal warping radius (eq.30 in Martin et al. 2011): where is the (base) disk viscosity at the stellar surface, Rd = (1 −  2 ) 1/2 is the average separation, Wood et al. 1997) is the (base) scale height of the disk at the stellar surface (Klement et al. 2017), and the power-law index  = 2/(11 − 2) ∼ 0.5 is given by  ∼ 3.5 which is the slope 11 of the disk (mid-plane) density profile (see eqs. 12, 16 and 29 in Martin et al. 2011).Substituting the formula of  tid (Eq.15) to the iniquity  trunc (1 − ) >  tid , we have in which  trunc is given by Eq. 6.Here we use  = 0.63 to evaluate the critical seperation  crit ∝  −/(3−1) ∼  −1 in Eq. 16, motivated by the finding in Cheng et al. (2014) that the observed populations of Be-XRBs with low ( s > 40 s) and high ( 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 = 2.66 considered in Martin et al. (2011).It will be shown below that our choice of  is justified by comparing 11 For simplicity, we adopt  = 3.5 throughout this work assuming that the part of the disk that interacts with the compact object can be well approximated with the steady-state solution (with constant  ej ).In fact, at the inner disk structure can vary significantly (with  ∼ 2 − 5) in response to the variations of  ej , magnetorotational instability and/or the presence of a companion object (Carciofi & Bjorkman 2008;Haubois et al. 2012 the mock population of Be-XRBs with observations (Fig. 4).The discrepancy here may be caused by the fact that the VDD flares (reaching / ≳ 0.1 at  ≳ 10 ★ ∼ 100 R ⊙ ) while the value in Martin et al. (2011) is derived for thin (/ ≪ 1), flat disks (Ogilvie 1999;Lodato & Price 2010).The disk flaring may reduce the viscosity for vertical shear and enhance vertical diffusion, making the disk more vulnerable for tidal warping (with  tid smaller by a factor of ∼ 2).
(2) The accretion rate criterion can be written as Here we adopt the typical radiative efficiencies  = 0.2 for NSs and  = 0.1 for BHs to calculate the Eddington rate  Edd (Eq.14).
We expect the transition Eddington ratio to be in the range  crit ∼ 0.05 − 0.2, where the upper limit is adopted in Okazaki et al. (2013) to explain the outburst strengths of Type I and II in simulations, while the lower limit is consistent with the theoretical thin disk 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.2  corr , because with this choice the  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  s > 40 s and  s < 40 s (see e. Given the above classification, we combine the optimistic duty cycles fduty in observations for the two types of outbursts, fduty,I = 0.3 and fduty,II = 0.1, with a physical limit  duty,max =  ej /  acc to estimate the (average) duty cycle as  duty = min( fduty ,  duty,max ) where is the mass ejection rate for  = 0.63 based on the results from ).In the mock population, 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  acc is comparable to f −1 duty  ej (such that  duty ≲ fduty ), 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  duty .Therefore, we also regard the Be-XRBs classified as Type I according to Eqs. 16 and 17 with  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, disk warping and/or disruption.Recent observations find that the outburst behaviors 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 fduty , and we have verified by numerical experiments that our results are insensitive to outburst classification.The reason is that the majority (∼ 60 − 70%) of X-ray emission comes from systems with  duty ≃  duty,max ≲ fduty in all Be-XRB populations considered in this work (see Sec. 5).

X-RAY SPECTRAL MODEL
Once  bol is known, we only need to determine the spectral shape to obtain the full (intrinsic) spectral energy distribution (SED) of Xray 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 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 2 Here we consider the photon energy range  ∼ 10 −4 − 10 4 keV assuming that this contains most of the power.We take the recent measurements of BC factors of the energy bands  ∼ 0.5 − 2, 2 − 10 and 12 − 15 keV by Anastasopoulou et al. (2022, see their table 4, A1, A6, A9 and A10) for both NS and BH HMXBs 15 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 spectral energy distributions (SEDs) that serve as the targets of spectral fitting.Since soft X-ray photons in the  ∼ 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 up-scattered component with the input spectrum  ,in via two parameters: the fraction  scatter of photons in the input (photon number) spectrum  in () ≡  in / = (ℎ) −1  ,in that are up scattered, and the 14 Theoretical calculation of the X-ray spectra from accreting compact objects (see e.g.Yang et al. 2017;Chashkina et al. 2019;Qiao & Liu 2020;Sokolova-Lapa et al. 2021;Pradhan et al. 2021;Mushtukov & Portegies Zwart 2023) is beyond the scope of our phenomenological model for Be-XRBs in BPS.We therefore adopt simple observation-based models (i.e.black body or thin disk + power law) to capture the general trends. 15The NS HMXB sample in (Anastasopoulou et al. 2022) is purely made of Be-XRBs.power-law index Γ in the spectrum of up-scattered photons: where  (,  0 ) = (Γ − 1)(/ 0 ) −Γ / 0 is the Green's function of inverse Compton scattering.The final specific luminosity can then be derived from the output photon number spectrum with  out () ≡  in / = (ℎ) −1  ,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 channeled by magnetic fields.The top row of Fig. 6 shows the resulting best-fit spectral models in the three regimes with  B  ∼ 0.4 − 0.9 keV,  scatter ∼ 1 and Γ ∼ 2.5 − 2.7.
While for BHs, we use the thin disk (TD) model 16 (Pringle 1981;Takhistov et al. 2022) where  max = 0.488  ,  o ≃  i ( i / o ) 3/4 and  max =  B  max /ℎ, given  i as the temperature at the inner edge ( i = 6  X / 2 ) 16 In principle, the thin disk solution is only valid at high accretion rates (e.g. ≳ 0.07, Pringle 1981;Takhistov et al. 2022), which are expected to cover most cases.Besides, we find that the contribution of BHs to the total X-ray output from Be-XRBs is no more than a few percent in all cases explored.Therefore, we do not consider the ADAF solution for BHs with lower accretion rates.
of the TD, and  o as the temperature at the outer disk boundary  o ∼  L1 (, ,  X ).Since very close binaries ( orb ≲ 7 days) and RLOF are forbidden for Be-XRBs (Panoglou et al. 2016(Panoglou et al. , 2018;;Rivinius 2019) by tidal forces that can slow down the rotation of donor stars, we have  o ≳ 2 R ⊙ and  B  o ≲ 10 −3 keV in most cases.Therefore  o is unimportant for X-rays that we are concerned with ( > 0.1 keV) and the input TD spectrum is controlled by a single parameter  i during the fitting process.The best-fit models for  X = 4 M ⊙ and  o = 10 R ⊙ are shown in the bottom row of Fig. 6 with  B  i ∼ 0.7−1.5 keV,  scatter ∼ 0.3−1 and Γ ∼ 2.1−2.8.

RESULTS
Combining the two choices for initial rotation parameter  0 : SR or FR (Sec.2.3), and the two models for VDD density Σ 0 : CS or OP (Sec.3.2), we have 4 models to explore, which are summarized in Table 2.For each model, we consider ten cases at  = 10 −4 , 3 × 10 −4 10 −3 , 0.0035, 0.005, 0.007, 0.01, 0.0142, 0.02 and 0.03, where  = 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  ≥ 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 properties17 , 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  2002), while the fast-rotating (FR) model uses a constant high value.The third column shows the choices of the (base) surface density Σ 0 of decretion disk (see Sec. 3.2), where the metallicity-dependent optimistic (OP) model uses a linear interpolation for log(Σ 0 ) between the MW fit for at  = 0.0142 (Eq.9) and the SMC fit at  = 0.0035 (Eq.10), while the conservative model adopts the MW fit at all metallicities.
Model  0 Σ 0 SR_CS  0 (  ★ ,  ) Σ 0 (  ★ ) (Eq. 9) FR_CS 0.9 Σ 0 (  ★ ) (Eq. 9) SR_OP  0 (  ★ ,  ) Σ 0 (  ★ ,  ) (Eqs. 9 and 10) FR_OP 0.9 Σ 0 (  ★ ,  ) (Eqs. 9 and 10) of X-ray outburst properties during the Be-XRB phase18 , which is a simplification of the reality that VDDs can be highly variable structures with disk dissipation/formation at timescales of a few years (Reig 2011;Rivinius et al. 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  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  corr are included in Appendix B

Formation efficiency
We first look into the formation efficiency of active Be-XRBs, N X , as a function of  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 timescale ( SF ≳ 100 Myr).Given  Be-XRBs predicted by a BPS run for a single-age stellar population of a total mass  tot , the formation efficiency can be written as where   is the duration of the Be-XRB phase for binary .
The results for all the 4 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  > 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 et al. 2004;Xing & Li 2021) that tidal truncation of the VDD is sharp, leading to a smaller disk boundary at  trunc (Eq.6) than the one adopted in our model ( crit , see Eq. 7).In both cases, N X generally increases with decreasing , but the evolution is not fully monotonic for all Be-XRBs, which has a small peak at  = 0.005 and a small dip at  = 0.01.The general trend is driven by the stronger stellar winds at higher  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 Sec. 5.2).We also find that nearly-circular ( ≤ 0.1) systems make up a significant fraction (∼ 40 − 80%) of active Be-XRBs at  ≲ 0.02.Naïvely, this seems in tension with the rareness (≲ 10%) of such Be-XRBs in observations (Cheng et al. 2014;Sidoli & Paizis 2018).However, observations are very likely incomplete at the low eccentricity end because the measurement of  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  bol ∼ 10 33 − 10 37 erg s −1 due to the suppressed accretion rate from disk truncation, and therefore, difficult to detect.In fact, they produce much less X-rays compared with their more eccentric counterparts such that ignoring them has little (up to a few percent) impact on the overall X-ray output.
Last but not the least, N 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  but remains below ∼ 40% and becomes even comparable to the uncertainties 19 in N X at  ≲ 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 et al. 2020Hastings et al. , 2021;;Wang et al. 2023), which is supported by observations that find a large fraction of classical Be stars with disk truncation (i.e.SED turndown; Klement et al. 2019), the lack of close Be binaries with MS companions (Bodensteiner et al. 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 et al. 2022).On the other hand, the difference between the FR and SR models is generally larger when we focus on Be-XRBs with  > 0.1 at  ≳ 10 −3 .The reason is that mass transfer is weaker in these systems leading to less efficient spin-up of the secondary star, as discussed below (Sec.5.2).Besides, the difference is smaller at lower , 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 ( ∼ 0.5 − 8 keV) X-ray luminosities  [0.5−8] keV > 10 38 erg s −1 and  [0.5−8] keV > 10 39 erg s −1 during outbursts, respectively, for which constraints from observations of HMXBs in nearby galaxies are available down to  ∼ 0.0003 (Douna et al. 2015;Lehmer et al. 2021).The results for the SR_CS model with  corr = 0.5 are shown in Fig. 8.We find that the numbers of Be-XRBs with  [0.5−8] keV > 10 39 erg s −1 predicted by our BPS runs are lower than those inferred from observations 20 that include all types of HMXBs (Lehmer et al. 2021) by about 40% at  ∼ 0.0003−0.02,which indicates significant (∼ 60%) contribution of Be-XRBs to ultra-luminous HMXBs.Besides, the predicted evolution with  for Be-XRBs follows a very similar trend as that seen in observations.The observed numbers from Lehmer et al. (2021) can be fully explained by our Be-XRB population at  ∼ 0.0003 − 0.02, if given a higher correction factor  corr = 0.8 than that assumed by default (  corr = 0.5), as shown in Appendix B. For instance, we have N X = 0.47 M −1 ⊙ yr at  = Z ⊙ given  corr = 0.8, which agrees 20 When comparing our Be-XRB populations with observations, we ignore the effects of anisotropic emission which can be important for ULXs with NS accetors (Wiktorowicz et al. 2019;Khan et al. 2022).To the zeroth order, such effects can be absorbed into the duty cycle  duty .
well with the observed value N X = 0.45 +0.06 −0.09 M −1 ⊙ yr (Kovlakas et al. 2020).However, at  = 0.03, the predicted number from Be-XRBs drops greatly with respect to the  = 0.02 case, becoming much lower than that in observations by a factor of a few.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  = 0.03.If this feature is true, it means that Be-XRBs play much less important roles at  > 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  [0.5−8] keV > 10 38 erg s −1 is similar but less clear considering the large uncertainties in observations (Douna et al. 2015).Although not shown here for conciseness, we find similar trends for Be-XRBs with  [0.5−8] keV > 10 37 erg s −1 .In this case, we have N X ∼ 10 M −1 ⊙ yr at  ∼ Z ⊙ , consistent with observations of HMXBs in nearby galaxies (Gilfanov et al. 2022).Moreover, for Be-XRBs with  [0.5−8] keV ≳ 10 35 erg s −1 , we predict N X ∼ 60−130 M −1 ⊙ yr at  ∼ 0.0035−0.0142,again consistent with the observed rate ∼ 135 M −1 ⊙ yr for all types of HMXBs in nearby star-forming galaxies (Mineo et al. 2012;Gilfanov et al. 2022;Lazzarini et al. 2023).

Masses and orbital parameters
For conciseness, in this section we only show the statistics of Be-XRBs for the SR_CS model at three representative metallicities,  = 10 −4 , 0.0035 (SMC) and 0.0142 (MW), to illustrate the general trends.Each Be-XRB  has a weight   / tot , so that the numbers in the following plots of marginalized distributions correspond to the expected numbers of binaries in the Be-XRB phase per unit SFR.
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 ( 1 ) is the progenitor of the compact object ( X ), and the (initial) secondary star ( 2 ) is the progenitor of the O/Be star ( ★ ).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 .The reason is that more massive primary stars tend to have more massive remnants (i.e. higher  X ) and also require more massive secondary stars to have stable mass transfer.Both the typical primary and secondary masses increase with  which, combined with the bottom-heavy IMF, explains the decreasing formation efficiency of Be-XRBs at higher  (Fig. 7).This trend is caused by two effects: The formation of NSs and BHs requires more massive primary stars at higher  due to stronger stellar winds, and the secondary stars that are more compact at lower  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 .This is caused by an effect that involves less massive stars given stronger winds at higher : 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  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  from  X,max ≃ 4.3 M ⊙ at  = 10 −4 to  X,max ≃ 1.6 M ⊙ at  = 0.03.There are no BHs (with  X > 2.2 M ⊙ ) at Here the (initial) primary star ( 1 ) is the progenitor of the compact object ( X ), and the (initial) secondary star ( 2 ) is the progenitor of the O/Be star ( ★ ).In the scatter plots of joint distribution, the area of each data point is proportional to the weight   / tot , which is also used to derive the marginalized distributions, where the solid, dashed and dotted lines mark the 10%, 50% and 90% percentiles.The dash-dotted lines in the scatter plots of initial-remnant mass relation for primary stars (bottom) approximately divide the Be-XRB population into two groups with lower and higher  X in which the primary stars experience different levels of mass loss (see the main text in Sec.5.2).The data point areas are normalized within each scatter plot separately and should not be compared across different panels.so that  ★ generally decreases with .The fraction of Be-XRBs with low-mass ( ★ < 10 M ⊙ ) O/Be stars (dashed curve in Fig. 10) decreases from ∼ 60% at  = 10 −4 to ∼ 1.5% at  = 0.001, and then increases quasi-monotonically with , reaching ∼ 87% at  = 0.03.This is consistent with the evolution of median O/Be star mass with  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  X − log  1 space with lower and higher  X (approximately divided by the dash-dotted lines in the bottom panels of Fig. 9) that are more distinct at higher .These two groups are also related to the complex features of Be-XRB mass distributions in the log  X − log  ★ space.The lower- X group is made of binaries with relatively low-mass ( 1 ≲ 25 M ⊙ , i.e. log( 1 [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 .The higher- 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- X group.We expect the primary stars in the lower- X group to completely lose their hydrogen envelope 21 and even undergo mass transfer during the helium HG phase.In this way, most of them form low-mass NSs ( X ∼ 1.3 M ⊙ ) 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- 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 21 Recently, a new sample of low-mass helium stars has been discovered in Magellanic Clouds, which are expected to originate from stars of initial masses ∼ 8 − 25 M ⊙ that are stripped by binary interactions (Drout et al. 2023;Gotberg et al. 2023).These helium stars are similar to the NS progenitors of Be-XRBs in the low- X group generated by our BPS runs.
becomes bi-polar at  ≳ 0.001, which is more obvious for O/Be stars from the higher- 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 through discussion on this aspect to future work.
Fig. 11 shows the orbital parameter distribution of Be-XRBs.Similarly to the results in Sec.5.1, nearly-circular ( ≲ 0.1) binaries with  ∼ 100 − 10 3 R ⊙ make up a significant fraction of Be-XRBs, reaching ∼ 60% at  = 10 −4 .This is a natural consequence of binary interactions that tend to circularize binary orbits and our optimistic definition of the VDD boundary (Eq.7) that allows faint objects with little accretion from beyond the tidal truncation radius to be counted as Be-XRBs (Sec.3.1).These systems mostly belong to the aforementioned lower- X group and contain lowmass ( X ∼ 1.3 M ⊙ ) 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  > 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 ∼ 30 Myr for  =  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 (Eq.6) such that nearly-circular systems are excluded.That is to say, wide ( ≳ 300 R ⊙ ) binaries with longer  need to have higher  to interact with the VDD at the pericenter.On the other hand, there is an upper limit of  that increases with  in close binaries ( ≲ 300 R ⊙ ) to avoid RLOF of the O/B star.Finally, we find that the distribution of  is broader at higher , while the median is almost constant at  ∼ 300 R ⊙ .The increase of the fraction of wide Be-XRBs at higher  can be explained by the stronger winds that widen binary orbits more significantly.The increasing relative importance of very close binaries ( ≲ 100 R ⊙ ) at higher  corresponds to the decreasing importance of the lower- X group that produces a stronger peak around  ∼ 500 R ⊙ at lower , and the trend will disappear if we exclude nearly-circular binaries.

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  = 0): where Be-XRB  lives from  ,ini to  ,fin with a (specific) luminosity  , during outbursts for a certain energy (band) 22 and a duty cycle  duty, , and Θ is the Heaviside step function.Fig. 12 shows the results for the 0.5 − 8 keV band from the SR_CS model with  corr = 0.5 22 The specific luminosity is defined as   ≡ /, while for a given Evolution of X-ray luminosity in the 0.5 − 8 keV band per unit stellar mass from Be-XRBs formed by an instantaneous starburst in the SR_CS model with  corr = 0.5, for  = 10 −4 (solid), 0.0035 (dashed), 0.0142 (dash-dotted) and 0.03 (dotted).Here  = 0 corresponds to the moment of starburst.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  ≲ 0.01, while the peak (as well as the onset) of X-ray emission shifts to later stages at higher , up to ∼ 20 Myr for  = 0.03.This is caused by the suppression of Be-XRBs from massive binaries by stellar winds that is more significant at higher  (see Sec. 5.2) and longer lifetimes of (initially less massive) primary stars in Be-XRB progenitors with higher .The X-ray luminosity drops by at least 2 orders of magnitude within 100 Myr post the peak.The most metalpoor model with  = 10 −4 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 lumi-nosity per unit SFR: where   =  ,fin −  ,ini is the duration of the Be-XRB phase for binary .Fig. 13 shows the full (intrinsic) SEDs in terms of L  for  = 10 −4 and 0.0142, in the SR_CS model with  corr = 0.5, where we also plot the contributions of different types of outbursts (Type I and II, see Sec. 3.2), accretion regimes (LH, HS and SE, see Sec. 4) and compact objects (NS and BH).In general, for the photon energy range  ∼ 0.5 − 8 keV that typically contains ∼ 60% 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%.This means that the spectral shape is almost independent of  in this energy range and the majority of X-ray emission comes from luminous (≳ 10 38 erg s −1 ) systems.In fact, SE systems contribute ∼ 55 − 68% of the total bolometric luminosity per unit SFR for  ∈ [10 −4 , 0.03] with mildly increasing importance at lower .
Similarly, Type II outbursts make up ∼ 56 − 71% 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)%.
For more energetic photons ( ≳ 8 keV), the contribution from NSs in the HS regime ( ∼ 0.05 − 2) becomes important due to their harder spectra, especially at higher , while for  ≲ 0.1 keV, BHs (mostly in the SE regime during Type II outbursts) dominate the spectrum when they exist in Be-XRBs at  ≲ 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  ≳ 20 and  ∼ 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 comparison with literature results (see below).Fig. 14 shows the X-ray luminosity per unit SFR as a function of  in these 4 bands for the SR_CS model with  corr = 0.5.The results for all the 4 models in Table 2 are summarized in Table 3.In the  = 10 −4 case where BHs are present, we show the NS and BH contributions with the dash-dot-dotted and densely-dashed curves.The shaded regions show different levels of absorption as defined in Fig. 6.

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 ) the X-ray luminosity (per unit SFR) from Be-XRBs in our SR_CS model with  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  ∼ 0.001 − 0.02, where the evolution with  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 disks 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  ≲ 0.001 and  ≳ 0.02, our results for Be-XRBs show stronger evolution with .In Fragos et al. (2013b), the Xray luminosity from other types of HMXBs is almost constant for  ≲ 0.001, while in our case, the X-ray luminosity from Be-XRBs increases by a factor of ∼ 2 from  = 0.001 to  = 10 −4 , indicating Be-XRB (0.1 2 keV) Be-XRB (0.5 2 keV) Be-XRB (2 10 keV) Be-XRB (0.5 8 keV) F13 (0.5 2 keV) F13 (2 10 keV) F13 (0.5 8 keV) L21 (0.5 8 keV) Figure 14.X-ray luminosity per unit SFR of Be-XRBs in the SR_CS model with  corr = 0.5, in the energy bands  ∼ 0.1 − 2 (solid), 0.5 − 2 (dashed), 2 − 10 (dash-dotted) and 0.5 − 8 keV (dotted).For comparison, we plot the best-fit models from BPS results for other types of HMXBs in Fragos et al. (2013b, F13) with the thick curves following the same line styles for the latter three bands.Here for the 0.5 − 2 keV and 2 − 10 keV bands, we use the fitting formulae and parameters in their eq.3 and table 2, while the results for the 0.5 − 8 keV band are taken from Lehmer et al. (2021,  that Be-XRBs likely play more important roles at lower metallicities.The rapid drop of X-ray luminosity at  > 0.02 in our case is caused by the stronger stellar winds23 at higher  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.
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  [0.5−2] keV / [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  ≳ 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).

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 our Be-XRBs can explain ∼ 60% 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  ∼ 0.02 to  ≲ 0.0003.This is consistent with the results in Sec.5.1 for number counts of ultra-luminous HMXBs with  [0.5−8] keV > 10 39 erg s −1 (Fig. 8).It will be shown in Appendix B that the observed number of (ultra-luminous) HMXBs as well as overall X-ray luminosity from Lehmer et al. (2021) can be fully reproduced by our SR_CS model of Be-XRBs at  ∼ 0.0003−0.02,if given a higher correction factor  corr = 0.8 than that assumed by default (  corr = 0.5).Considering the uncertainties in our results embodied by the correction factor  corr ∼ 0.25 − 1, we conclude that Be-XRBs can contribute at least 30% and up to 100% of the X-ray emission of observed HMXBs at  ∼ 0.0003 − 0.02, and meanwhile make up a similar fraction of the observed ultra-luminous HMXB population.In fact, we find that systems with  [0.5−8] keV > 10 39 erg s −1 only account for ∼ 15 − 50% of the total X-ray output in our Be-XRB populations at  ∼ 0.0003 − 0.02 from the SR_CS model with  corr ∼ 0.25 − 0.8.It is therefore non-trivial that we obtain similar contributions Be-XRBs in the number count of ultra-luminous HMXBs and total X-ray luminosity with respect to observations across 2 orders of magnitude in metallicity.
Finally, we calculate the X-ray luminosity per luminous Be-XRB with  [0.5−8] keV > 10 38 erg s −1 in the 0.  in the range of [0.25, 1] (not shown for conciseness).Our results are consistent with observations of HMXBs in nearby galaxies for  ∼ 0.0004 − 0.03 (Douna et al. 2015, see their fig.5), which further highlights the striking similarity between the luminous X-ray sources in our Be-XRB populations and those in observed HMXB populations.

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 et al. 2013).Using the standard BSE models (Hurley et al. 2002) in the BPS code bi-nary_c, we evolve large populations of randomly sampled massive binaries in the (absolute) metallicity range  ∈ [10 −4 , 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.
(i) Be-XRBs can contribute a significant fraction (∼ 60% in our fiducial case, up to 100% and at least ∼ 30% considering uncertainties) of the total X-ray output from HMXBs observed in nearby galaxies (Lehmer et al. 2021) with a similar metallicity dependence for  ∼ 0.0003 − 0.02, assuming that observations are complete (see Sec. 5.3.2 and Appendix B).In our fiducial model with no  dependence of VDD densities nor enhanced initial stellar rotation, the X-ray luminosity per unit SFR decreases by a factor of ∼ 8 from  = 0.0003 to  = 0.02, which is mainly driven by the stronger stellar winds at higher  that reduce the number of binary stars able to become luminous Be-XRBs.
(ii) Our results for Be-XRBs are also similar to the -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  ∼ 0.001 − 0.02 (Sec.5.3.1).However, at  ≲ 0.001, the X-ray luminosity from Be-XRBs keeps increasing towards lower  (e.g. by a factor of ∼ 2 from  = 0.001 to  = 10 −4 ), while there is almost no evolution in that from other types of HMXBs predicted by Fragos et al. (2013b).This indicates that Be-XRBs likely play more important roles at lower metallicities, such that they can have strong impact on the thermal history and thus 21-cm signal from the IGM at Cosmic Dawn.
(iii) Luminous (≳ 10 38 erg s −1 ) systems with non-periodic, major outbursts of Type II (see Sec. 3.2.2 for a detailed definition) and super-Eddington accretion (mostly on eccentric orbits with  ≳ 0.4), including ULXs (≳ 10 39 erg s −1 ), are important X-ray powerhouses in our Be-XRB populations across the metallicity range considered.In our fiducial model, ∼ 55 − 75% of the total X-ray output of Be-XRBs comes from luminous sources with  [0.5−8] keV > 10 38 erg s −1 , while ULXs with  [0.5−8] keV > 10 39 erg s −1 contribute ∼ 24 − 41% of the total output.Interestingly, we find that our Be-XRBs can account for similar fractions (∼ 30 − 100%) of the number count of luminous and ultra-luminous HMXBs as their contribution to the total X-ray luminosity of HMXBs in observations (Douna et al. 2015;Kovlakas et al. 2020;Lehmer et al. 2021).Such luminous systems have been overlooked in previous BPS studies of Be-XRBs (e.g. Zuo et al. 2014;Misra et al. 2023), which use a simple empirical scaling law between outburst X-ray luminosity and orbital period (Dai et al. 2006, see our Eq.11 and Fig. 3) that imposes a (likely artificial) cutoff at ∼ 10 38 erg s −1 , while our im- Figure 16.X-ray luminosity function (per unit SFR) of Be-XRBs in the 0.5-8 keV band predicted by our fiducial model (SR_CS with  corr = 0.5) for  = 10 −4 (solid), 0.0035 (dashed) and 0.0142 (dash-dotted).Here each Be-XRB is weighted by  duty,   / tot .The power-law fit ( N X / X ∝  −1.6 X ) to the average X-ray luminosity function of observed HMXBs in star-forming galaxies (Mineo et al. 2012, see their fig.5) is plotted with the dotted line for comparison.The X-ray luminosity functions of our Be-XRBs are generally consistent with the observed power-law shape for  X ∼ 10 36 − 10 39 erg −1 , while the number fractions of systems in both the faint and luminous ends are lower for our Be-XRBs compared with the whole population of HMXBs in observations.proved 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 erg s −1 , as shown in Fig. 16.
(iv) 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 Sec. 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 ∼ 30%.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 et al. 2020;Dorigo Jones et al. 2020;Hastings et al. 2020Hastings et al. , 2021;;Dallas et al. 2022;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.

OUTLOOK
To quantify more accurately the roles played by Be-XRBs, particularly in the ULX population (e.g.Kaaret et al. 2017;Kovlakas et al. 2020;Fabrika et al. 2021;Walton et al. 2022;King et al. 2023;Sal-vaggio et al. 2023;Tranin et al. 2023), 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).
(i) 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 (Eq.18) and viscosity.However, in reality VDDs of O/Be stars, especially under the influence of companions, can be highly dynamical structures (Carciofi & Bjorkman 2008;Haubois et al. 2012;Krtička et al. 2015;Panoglou et al. 2016;Vieira et al. 2017;Rímulo et al. 2018), with variable mass ejection and disk dissipation/formation episodes of a few years (Reig 2011).Such dynamical evolution of disk structures are particularly important for the most luminous Be-XRBs with non-periodic, strong Type II outbursts (see Sec. 3.2.2).In fact, according to our steady-state model, the majority (∼ 60 − 70%) of X-ray emission from Be-XRBs are produced by the systems in which the compact object is able to accrete almost all (≳ 90%) of the materials ejected from the O/B star, so that significant variations of disk 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 et al. 2023).
(ii) 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, 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 et al. 2014;Shao & Li 2014, 2020;Vinciguerra et al. 2020;Xing & Li 2021;Misra et al. 2023).
(iii) 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 observations24 and simulations (e.g.Gunawardhana et al. 2011;Marks et al. 2012;Duchêne & Kraus 2013;Moe & Di Stefano 2017;Moe et al. 2019;Jeřábková et al. 2018;Lacchin et al. 2020;Chon et al. 2021Chon et al. , 2022;;Tanvir & Krumholz 2023;Rusakov et al. 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.
(iv) 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 et al. 2017;Groh et al. 2019;Murphy et al. 2021) that fast rotation can impact the mass loss, timescales of evolution phases, stellar structure, nucleosynthesis and remnant masses in a complex and metallicity-dependent manner, particularly for initially fast-rotating stars, which can also change the properties of Be-XRBs that by definition contain fast-rotating O/B stars.
(v) 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 disk and ADAF) that may also correlate with the stellar and orbital parameters as well as types of outbursts (Martin et al. 2011;Okazaki et al. 2013;Cheng et al. 2014;Haberl & Sturm 2016;Xu & Li 2019;Franchini & Martin 2021).
Indeed, robust modeling of XRBs with BPS is a challenging task 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 fast-rotating massive stars (for candidate scenarios, see e.g.Granada et al. 2013;Hastings et al. 2020;Zhao & Fuller 2020;Cranmer 2009;Rogers et al. 2013;Lee et al. 2014;Ressler 2021).Nevertheless, such theoretical efforts are worthwhile considering the important roles played by XRBs in the Epoch of Reionisation and Cosmic Dawn (e.g. Fialkov et al. 2014;Fialkov & Barkana 2014;Pacucci et al. 2014), as well as the wealth of high- observational data to come in the next decades for the 21cm 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 et al. 2021) and SKA (Koopmans et al. 2015) which will place constraints on high- 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 et al. 2023;Rossland et al. 2023).We will also extend the metallicity range down to the extremely metal-poor ( ≲ 10 −6 ) regime of the first stars 25 , which are expected to be promising progenitors of Be-XRBs (Sartorio et al. 2023) given their compact and fast-rotating 26 nature that favors 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 syn- 25 The properties of the first star binaries are still highly uncertain in the lack of direct observations.Recent radiative hydrodynamic simulations (Sugimura et al. 2020;Park et al. 2023a,b;Sugimura et al. 2023) have found that close binaries are likely rare for the first stars (due to outward migrations of protostars and their circumstellar disks by accretion of gas with high angular momentum).This implies that the formation efficiency and X-ray outputs of XRBs from the first stars can be much lower than those of XRBs from metal-rich stars (considered in this paper) that are dominated by close binaries (Liu et al. 2021).Nevertheless, current simulations are still limited (by resolution and treatments of magnetic fields), whose results are yet to be confirmed by observations. 26According to hydrodynamic simulations of primordial star formation, Pop III stars tend to be born as fast ( 0 ∼ 0.5 − 1) rotators (Stacy et al. 2011(Stacy et al. , 2013;;Hirano & Bromm 2018), which is also supported by stellar archaeology observations (Chiappini et al. 2006;Chiappini 2013;Maeder et al. 2015;Choplin et al. 2019;K et al. 2023).
ergy between 21-cm and gravitational wave observations to better constrain binary stellar evolution in high- galaxies 27 .over-prediction of long-period ( orb ≳ 200 days) systems may be explained by incompleteness in observations but the discrepancy at  orb ∼ 10 − 100 days 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  BeXRB = 98 +18 −15 , allowing more space for the potential incompleteness of observations.With optimistic star formation rates of the SMC, the predicted number of Be-XRBs is larger than the observed number in almost every bin of  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 ( orb ∼ 10 days) systems is over-predicted 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 (Eqs.2-4) rather than the default log-flat initial separation distribution, for the SR_CS model we obtain  BeXRB = 65 +12 −10 , lower than the number 74 +13 −11 in the default case, but still consistent with observations.The predicted  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.

Figure 1 .
Figure 1.Relation between the peak accretion rate and pericenter distance based on the simulation results in Brown et al. (2019, see their fig.7)for eccentricity  = 0 (solid), 0.2 (thin dashed), 0.4 (dash-dotted) and 0.6 (dotted), in the unit of  ej,−10 M ⊙ yr −1 given  ej,−10 ≡  ej /(10 −10 M ⊙ yr −1 ).The relation can be well described by a power-law  acc,sim ∝ [ (1 − ) ] −2 as shown by the thick dashed line.These simulations consider a NS with  X = 1.4 M ⊙ around a Be star of  ★ = 18 M ⊙ and  ★ = 7 R ⊙ with  = 0.63 for  orb ∼ 40 − 400 days and derive the median peak accretion rate from 5 orbits after the system settles to steady state.

Figure 2 .
Figure 2. Relation between base disk surface density and Be star mass.The dots denote the 80 Be stars observed in the MW fromVieira et al. (2017,  V17), which can be fit with a power-law relation log(Σ 0 [g cm −2 ] ) ≃ 1.44 log(  ★ [M ⊙ ] ) − 2.37 of ≃ 0.52 dex errors (shaded region around the solid line).The squares show the median base disk surface densities in 8 mass bins measured from 54 Be stars in the SMC(Rímulo et al.  2018, R18), where the errorbars show the bin size and 25-75% percentiles of Σ 0 in each bin.These results can also be described with a power-law log(Σ 0 [g cm −2 ] ) ≃ 1.03 log(  ★ [M ⊙ ] ) − 0.99 of ≃ 0.17 dex errors (shaded region around the dashed line).
Figure 2. Relation between base disk surface density and Be star mass.The dots denote the 80 Be stars observed in the MW fromVieira et al. (2017,  V17), which can be fit with a power-law relation log(Σ 0 [g cm −2 ] ) ≃ 1.44 log(  ★ [M ⊙ ] ) − 2.37 of ≃ 0.52 dex errors (shaded region around the solid line).The squares show the median base disk surface densities in 8 mass bins measured from 54 Be stars in the SMC(Rímulo et al.  2018, R18), where the errorbars show the bin size and 25-75% percentiles of Σ 0 in each bin.These results can also be described with a power-law log(Σ 0 [g cm −2 ] ) ≃ 1.03 log(  ★ [M ⊙ ] ) − 0.99 of ≃ 0.17 dex errors (shaded region around the dashed line).

Figure 3 .
Figure3.Relation between outburst X-ray luminosity and orbital period for the mock population of Be-XRBs with  orb ∼ 10 − 300 days,  ∼ 0 − 0.6 and  X > 10 34 erg s −1 , given the calibration paremeter  X = 0.25.Type I/II Be-XRBs are shown in blue/orange.These mock Be-XRBs are identified from 10000 randomly generated NS-Be star binaries with a logflat distribution of separations for  ∈ [10−1000] R ⊙ , a uniform distribution of eccentricities for  ∈ [0, 0.6] and a log-flat distribution of Be star masses for  ★ ∈ [6, 20] M ⊙ given a fixed NS mass  X = 1.4 M ⊙ .We show the observed Be-XRBs compiled byRaguzova & Popov (2005, RP05, excluding V615 Cas and 2206+543) with filled squares and unfilled triangles.The filled squares have  orb measured by observations, while for the unfilled triangles without  orb measurements, we use the empirical scaling law log(  orb [day] ) = 0.4329 log(  s [s] ) + 1.043(Vinciguerra et al. 2020) to estimate  orb given the spin period  s .The dashed line represents the best fit to observations fromDai et al. (2006) with the shaded region denoting the 1 uncertainties in the fitting parameters (see Eq. 11).The solid line shows the best-fit power-law model for our mock population.It turns out that the mock population reproduces well observational data with a slightly shallower best fit and a similar scatter given  X = 0.25.
g. fig.3inCheng et al. 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 byCheng et al. (2014), it does not reproduce the quasi-bi-modal feature of the observed  X distribution.

Figure 6 .
Figure6.Best-fit spectral models for NSs (top) and BHs (bottom) in the low-hard (LH,  < 0.05, left), high-soft (HS,  ∼ 0.05 − 2, middle) and super-Eddington (SE,  > 2, right) regimes, in terms of   /(ℎ bol ).Here we assume  X = 4 M ⊙ and  o = 10 R ⊙ for the BH spectra.In each case the input spectrum is shown with the solid curve, the final spectrum with the dashed curve and the up-scattered component with the dash-dotted curve, respectively.The data points with error bars show the reference SEDs inferred from the BC factor measurements inAnastasopoulou et al. (2022).The dark ( ≲ 0.1 keV), intermediate dark ( ≲ 0.5 keV) and light ( ∼ 0.5 − 2 keV) shaded regions denote the parts of the spectrum that are expected to be absorbed by the ISM in typical high- HMXB-hosting minihaloes(Sartorio et al. 2023), atomic cooling haloes(Das et al. 2017), and the IGM at  ∼ 15(Pritchard & Furlanetto 2007).

Figure 7 .
Figure7.Number of Be-XRBs (in the outburst phase) per unit SFR as a function of metallicity, for the SR_CS (solid), FR_CS (dashed), SR_OP (dash-dotted) and FR_OP (dotted) models.The results for all Be-XRBs are shown with the thin curves marked by '★', while the unmarked thick curves show the results for Be-XRBs with  > 0.1.The thin vertical lines label the metallicities of the MW, LMC and SMC (from right to left).

Figure 8 .
Figure 8. Number of (ultra-)luminous Be-XRBs (in outbursts) per unit SFR as a function of metallicity in the SR_CS model with  corr = 0.5.The BPS results for luminous ( [0.5−8] keV > 10 38 erg s −1 ) and ultra-luminous ( [0.5−8] keV > 10 39 erg s −1 ) sources are shown by the dashed and solid curves, respectively.For comparison, we show the observational data from Douna et al. (2015, D15, see their fig.4) for  [0.5−8] keV > 10 38 erg s −1 with the crosses.The uncertainties in these data are typically large (∼ 0.5 dex), as implied by their scatter around similar metallicities.The number counts for ultra-luminous sources ( [0.5−8] keV > 10 39 erg s −1 ) from the observed HMXB sample in Lehmer et al. (2021, L21, see their table3and fig.6) are denoted by the data points with (1) error bars, and the shaded region denotes the 16-84% confidence range obtained with mock populations of XRBs sampled from their best-fit model for -dependent X-ray luminosity functions.The thin vertical lines label the metallicities of the MW, LMC and SMC (from right to left).

Figure 9 .
Figure 9. Mass distributions of Be-XRBs (top) and their progenitors (middle), and the initial-remnant mass relation for primary stars of Be-XRB progenitors (bottom) in the SR_CS model for  = 10 −4 (left),  = 0.0035 (middle) and  = 0.0142 (right).Here the (initial) primary star ( 1 ) is the progenitor of the compact object ( X ), and the (initial) secondary star ( 2 ) is the progenitor of the O/Be star ( ★ ).In the scatter plots of joint distribution, the area of each data point is proportional to the weight   / tot , which is also used to derive the marginalized distributions, where the solid, dashed and dotted lines mark the 10%, 50% and 90% percentiles.The dash-dotted lines in the scatter plots of initial-remnant mass relation for primary stars (bottom) approximately divide the Be-XRB population into two groups with lower and higher  X in which the primary stars experience different levels of mass loss (see the main text in Sec.5.2).The data point areas are normalized within each scatter plot separately and should not be compared across different panels.

𝑍 ≳ 0 Figure 10 .
Figure 10.Median O/Be star mass (solid curve for the left axis) in Be-XRBs and fraction (dashed curve for the right axis) of Be-XRBs with low-mass (< 10 M ⊙ ) O/Be stars in the SR_CS model.The thin vertical lines label the metallicities of the MW, LMC and SMC (from right to left).

Figure 11 .
Figure 11.Distribution of Be-XRB orbital parameters in the SR_CS model for  = 10 −4 (left),  = 0.0035 (middle) and  = 0.0142 (right).In each scatter plot of joint distribution, the area of each data point is proportional to the weight   / tot , and the solid, dashed and dotted contours enclose 90, 50 and 10% of systems from a Gaussian smoothed density field produced by the weighted data.This weight is also considered in the marginalized distributions, where the solid, dashed and dotted lines mark the 10%, 50% and 90% percentiles.The data point areas are normalized within each scatter plot separately and should not be compared across different panels.

Figure 13 .
Figure 13.SED of Be-XRBs per unit SFR from the SR_CS model with  corr = 0.5 at  = 10 −4 (top) and 0.0142 (bottom).The total spectrum is shown with the thick solid curve.Contributions from Type I and II outbursts are denoted by the solid and dashed curves.Components of the LH, HS and SE states are shown with the dash-dotted, dotted and long-dashed curves.In the  = 10 −4 case where BHs are present, we show the NS and BH contributions with the dash-dot-dotted and densely-dashed curves.The shaded regions show different levels of absorption as defined in Fig.6.
Figure14.X-ray luminosity per unit SFR of Be-XRBs in the SR_CS model with  corr = 0.5, in the energy bands  ∼ 0.1 − 2 (solid), 0.5 − 2 (dashed), 2 − 10 (dash-dotted) and 0.5 − 8 keV (dotted).For comparison, we plot the best-fit models from BPS results for other types of HMXBs inFragos et al.  (2013b, F13)  with the thick curves following the same line styles for the latter three bands.Here for the 0.5 − 2 keV and 2 − 10 keV bands, we use the fitting formulae and parameters in their eq.3 and table 2, while the results for the 0.5 − 8 keV band are taken from Lehmer et al. (2021, see their fig.4).The observational results for the 0.5 − 8 keV band from Lehmer et al. (2021, L21, see their fig.4 and table 3) are shown with the shaded region.The thin vertical lines label the metallicities of the MW, LMC and SMC (from right to left).

Figure 15 .
Figure 15.X-ray luminosity per luminous Be-XRBs with  [0.5−8] keV > 10 38 erg s −1 , for the SR_CS (solid), FR_CS (dashed), SR_OP (dash-dotted) and FR_OP (dotted) models with  corr = 0.5.The observational data from nearby galaxies and the corresponding power-law fit from Douna et al. (2015, D15, see their fig.5) are shown with the crosses and dashed line surrounded by the shaded region.The thin vertical lines label the metallicities of the MW, LMC and SMC (from right to left).

Table 1 .
Key physical quantities. rot / Kep with the initial value denoted by  0 Edd Eddington limit of accretion rate (Eq.14)  ≡  acc /  Edd , Eddington ratio during outbursts  duty effective fraction of time the Be-XRB spends in outbursts  X outburst X-ray luminosity for a certain band  X calibration parameter for the observed  X −  orb relation  corr correction factor for outburst luminosity (Sec.3.2)  lifetime of the Be-XRB  tot Richardson et al. 2023)hough such systems have been found in observations (e.g.CPD-29 2176,Richardson et al. 2023).The reason is that tidal truncation at  trunc is not an absolute cutoff.What happens is that materials accumulate within  trunc and the disk density profile becomes much steeper beyond  trunc than within  trunc 9In the calculation of  trunc for each Be-XRB, we further introduce a Distribution of outburst X-ray luminosity for the mock population of Be-XRBs (dotted contour, see Fig.3) in comparison with that of the 71 observed Be-XRBs (long-dashed contour) complied byCheng et al. (2014,  C14, see their fig.3).The observed distribution is re-normalized to have the same total number of Be-XRBs.For the mock population, contributions from Type I and II outbursts are plotted with the right and left shaded histograms in blue and orange, respectively, while for the observed population, the sub groups of X-ray pulsars with spin periods  s above and below 40 seconds are shown with the solid and dashed contours.The  X distribution and fraction of Type I (II) events from our mock population are generally consistent with those of the observed Be-XRBs with ; Krtička et al. 2015; Panoglou et al. 2016; Vieira et al. 2017; Rímulo et al. 2018).s > (<) 40 s.
(Reig 2011;al.2012;Rímuloetal. 20Xu & Li 2019es the simple requirement that the compact object does not accrete more than what is ejected from the O/B star.When  duty,max < fduty , 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 disk 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).Itmeans that mass replenishment of VDDs is much more efficient 13 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 disk timescales 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  duty,I ∼ 0.1 − 0.3 and  duty,II ∼ 0.01 − 0.1, which is generally consistent with observations(Reig 2011; Sidoli  & Paizis 2018;Xu & Li 2019

Table 2 .
Summary of models.The second column shows the choices of the initial rotation parameter  0 ≡  rot,0 / Kep,0 given the initial rotation velocity  rot,0 and Keplerian velocity  Kep,0 at the stellar equator (see Sec. 2.3), where the slowly-rotating (SR) model  0 (  ★ ,  ) with mass and metallicity dependence is based on Hurley et al. (2000, see their sec.7.2) and Hurley et al. (