© The Authors 2023

1. Introduction

X-ray binaries (XRBs) are systems with bright X-ray emission arising from the accretion of matter from a non-degenerate star onto a compact object (CO), such as a neutron star (NS) or a black hole (BH). Observations of XRBs with X-ray telescopes, including Chandra, NuSTAR, and XMM-Newton, have significantly contributed to our understanding of accreting CO sources (e.g., Sarazin et al. 2000; Fabbiano 2006; Bozzo et al. 2012; Douna et al. 2015; Kaaret et al. 2017; Vulic et al. 2018; Lazzarini et al. 2019). Furthermore, the accumulation of large samples of extragalactic XRBs has allowed for statistical studies of their populations. X-ray luminosity functions (XLFs) are the primary tool for studying the statistical properties of observed XRBs. There is observational evidence that the properties of XRB populations (through the study of their XLFs) depend heavily on the host galaxy properties. For instance, the number and integrated luminosity of high-mass X-ray binaries (HMXBs; XRBs with typical donor star masses ≳8 M_⊙) scale with the star-formation rate (SFR) of the host galaxy due to their short lifetimes of a few million years (Grimm et al. 2003; Ranalli et al. 2003; Lehmer et al. 2010; Antoniou et al. 2010; Mineo et al. 2012; Antoniou & Zezas 2016). Low-mass X-ray binaries (LMXBs; with typical donor star masses ∼1 M_⊙), on the other hand, tend to have longer lifetimes (∼10⁹ yr) and thus do not correlate with the current local SFR. Instead, the integrated luminosity of LMXBs depends on the entire past star-formation history (SFH) of the host galaxy and scales with its total stellar mass (integrated SFR over time; Ghosh & White 2001; Ptak et al. 2001; Grimm et al. 2002; Gilfanov 2004; Lehmer et al. 2010; Boroson et al. 2011).

Observations of XRB populations in star-forming galaxies show that the number of XRBs, and their integrated X-ray luminosity, anticorrelate with the metallicity of the host galaxy (Kaaret et al. 2011; Basu-Zych et al. 2013, 2016; Brorby et al. 2016; Lehmer et al. 2019, 2021; Kovlakas et al. 2020). Metallicity affects the evolution of a binary through several processes. Lower metal content in stellar atmospheres leads to a decreased line-driven wind-mass loss (Vink et al. 2001). This results, on average, in relatively heavy COs (exceeding 30 M_⊙) at the end of the binary’s life (Mapelli et al. 2009; Zampieri & Roberts 2009). Furthermore, weaker stellar winds result in decreased orbital expansion due to orbital angular-momentum losses. Both of these factors lead to a more efficient accretion onto the COs and hence, on average, to more luminous XRB populations (Fragos et al. 2013a,b; Basu-Zych et al. 2016). In addition, low metallicity stars tend to expand less and later in their evolution than those with higher metallicities, reducing the chances of encountering a dynamically unstable mass-transfer episode, or if they do so, increasing their chances of surviving the common-envelope (CE) phases (Linden et al. 2010; Basu-Zych et al. 2016). Both effects tend to increase the formation efficiency of XRBs with decreasing metallicity. Dray (2006) used Monte Carlo simulations and found the number of HMXBs increased by a factor of 3 in the Small Magellanic Cloud (SMC), with Z/Z_⊙ = 0.2, compared to solar (however, age effects are also in play; Antoniou et al. 2019). Similarly, Linden et al. (2010) used the StarTrack code (Belczynski et al. 2008) and found an increase of a factor of 3.5 when going to Z/Z_⊙ = 0.2 for young bright HMXBs. The increased number of XRBs in low-metallicity environments, coupled with the metallicity evolution of the Universe (e.g., Fragos et al. 2013a,b), points to the possibility of XRBs having a significant effect on the ionization of the interstellar medium in the local Universe as well as the intergalactic medium in the early Universe (e.g., Madau & Fragos 2017; Schaerer et al. 2019; Kaur et al. 2022; Kovlakas et al. 2022).

Generally, the XLFs of XRB populations are fitted with either a single or a broken power law depending on the relative contribution of HMXBs and LMXBs in the XRB population. Studying observations of XRBs in the Milky Way and Magellanic Clouds, Grimm et al. (2003) proposed first that the HMXB XLF was in the form of a single, smooth power law. They found no specific features in the XLF correlating with the Eddington limits of NSs and BHs but did find evidence of a high luminosity cutoff at 10⁴⁰ erg s⁻¹. Belczynski et al. (2004) constructed the first XRB population synthesis models that closely reproduced the XLF of NGC 1569 (a starburst galaxy) when using a combination of an old (metal-poor) and young (metal-rich) stellar population. The findings of Grimm et al. (2003) were later confirmed by Mineo et al. (2012), who studied HMXB populations in 29 nearby star-forming galaxies. They found that the constructed XLF follows a power law with a slope of 1.6 in the X-ray luminosity range log(L_X/erg s⁻¹)∼35 to 40. Zuo et al. (2014) used an updated version of the Binary Stellar Evolution code (BSE; Hurley et al. 2000, 2002; Liu & Li 2007; Zuo et al. 2008) to study XLFs of HMXBs in star-forming galaxies. They found that for a wide range of model parameters, the XRB population is dominated by wind-fed XRBs with BH accretors, which produces an XLF that closely resembles the observed one.

On the contrary, the XLFs of LMXB-dominated populations tend to fit better with broken power laws. Kim & Fabbiano (2004) used a broken power law to fit the XLF of extragalactic LMXB populations observed with Chandra and found that the best fit involved a break at ∼5 × 10³⁸ erg s⁻¹, around the Eddington luminosity of a CO with mass 3.0 M_⊙, and power law slopes of ∼1.8 and ∼2.8 below and above the break, respectively. Fragos et al. (2008) constructed the theoretical XLFs of LMXBs and compared them with the elliptical galaxies NGC 3379 and NGC 4278 using the population synthesis code StarTrack (Belczynski et al. 2008). They found that the XLF is dominated by binaries with NS accretors and giant (main sequence) donors above (below) luminosities of ∼10³⁷ erg s⁻¹. Furthermore, they concluded that although the XLFs depend on the adopted combination of the population synthesis parameters, the treatment of the transient sources is a crucial factor for the determination of the XLF – for instance, modeling their variability instead of assuming a constant duty cycle for all sources, which would affect the number of LMXBs identified (Fragos et al. 2009). Studying the XRBs in the bulge and ring of the ring galaxy NGC 1291 (which is dominated by an old stellar population with little recent star formation), Luo et al. (2012) reached similar conclusions. In a related study, Tzanavaris et al. (2013) compared synthetic with observed XLFs from a sample of 12 nearby, late-type galaxies from the Spitzer Infrared Nearby Galaxy Survey (Kennicutt et al. 2003). They performed a broad model parameter space study, varying parameters such as the CE efficiency, stellar wind strengths, and supernova (SN) kick velocity distributions, and explored the relative contribution to the XLF from different XRB subpopulations as a function of the SFH of the parent stellar population.

The evolution with the age of the XLF shape and integrated luminosity of an LMXB population has been a topic of both theoretical and observational studies. Theoretical studies of LMXB populations predict that their XLF flattens in younger stellar populations, while the integrated X-ray luminosity can increase by up to one order of magnitude (e.g., Fragos et al. 2008, 2013a,b). Initial observational results contradicted the theoretical predictions, with Zhang et al. (2012) finding a higher number of LMXBs per unit of stellar mass in older elliptical galaxies compared to younger ones. This discrepancy arose from the inclusion of globular cluster LMXBs in the Zhang et al. (2012) study and the underlying correlation of a specific globular cluster frequency with stellar mass in the observed galaxy sample, which biased the results. Other observational investigations (Kim et al. 2009; Lehmer et al. 2014, 2020) found that when removing the globular cluster LMXB subpopulations, more luminous LMXBs indeed dominate younger populations, confirming the theoretical predictions. Lehmer et al. (2017) studied far-UV to far-IR data of XRB sources in the spiral galaxy M 51 and constructed an age-dependent XLF. They found that both single and broken power laws were good fits for the data. Generally, XLFs (for X-ray luminosities ≳10³⁷ erg s⁻¹) of older populations are found to be steeper than for younger populations, which complements the fact that older populations do not reach X-ray luminosities as high as those reached by younger populations that have higher numbers of HMXBs (Kim & Fabbiano 2010; Lehmer et al. 2014; Gilbertson et al. 2022). Lehmer et al. (2019) used a sample of 38 local galaxies (within a radius of 30 Mpc) from Chandra Advanced CCD Imaging Spectrometer (ACIS) imaging data, covering a broad range of SFRs and stellar masses to constrain the XLFs of HMXBs and LMXBs. They fitted both single and broken power law models to their data. For HMXB XLFs, they found the shape to be more complex than a single power law. The XLF slope rapidly declined beyond luminosities of 10³⁸ erg s⁻¹, which is consistent with results from Tzanavaris et al. (2013), Zuo et al. (2014), and Artale et al. (2019).

The aforementioned studies imply that generalizing the observed XLF of an arbitrary XRB population is a difficult task and involves many implicit assumptions about the underlying binary population. Many of these studies focused on either LMXBs or HMXBs, as they have distinct observational and physical characteristics (Tauris & van den Heuvel 2006). High-mass XRBs are some of the brightest X-ray sources in the local Universe. In our study, we focus on young XRB populations that are dominated by HMXBs and the effects of various physical assumptions on their XLFs. There are generally three broad types of HMXBs: (i) binaries with a super-giant donor undergoing wind-fed accretion, (ii) binaries with a super-giant donor star undergoing Roche-lobe overflow (RLO), and (iii) binaries with a fast spinning Be-star donor. In the case of wind-fed accretion, a fraction of the mass leaving the donor star in the form of strong stellar winds is captured by the gravitational pull of the CO (Bondi & Hoyle 1944). For XRBs undergoing RLO, the donor overfills its Roche lobe and transfers mass to the CO through the inner Lagrangian point. The transferred material carries enough angular momentum, which prevents it from being accreted immediately, forming an accretion disk around the CO. Be XRBs are characterized by rapidly spinning O/B donors with a disk around their equator, called the decretion disk (Balona & Ozuyar 2020). If the binary is close enough that the CO interacts with the decretion disk, the CO may accrete material from the donor without the donor filling its Roche lobe or the presence of any strong stellar wind.

Naturally, the observed properties of XRB populations will be affected by the formation environment and subsequent evolution of these sources. Some studies have explored the effects of binary physics on XRB observables. Brandt & Podsiadlowski (1995) studied the post-SN kick velocities for XRB progenitors and found a tight correlation between the post-SN orbital periods and eccentricities of LMXBs, with LMXBs that have periods longer than a few days getting large eccentricities. They also found that a significant fraction of HMXBs would have precessing accretion disks, which would have implications on XRB observations. Linden et al. (2011) and Zuo & Li (2014) investigated the effect of the CE phase on XRB populations and find that the constrained quantities for the CE efficiency depend on various assumptions, such as SN kicks, metallicity, or bolometric correction factors, that introduce uncertainties.

Despite the uncertainties, it is critical that theoretical models are able to reproduce observed XLFs if we aim to use them to infer the properties of the underlying XRB populations. X-ray luminosity functions contain signatures of the prior evolution of the XRBs that can be leveraged to constrain uncertain stellar and binary physics. The different processes involved in binary evolution directly impact the properties of additional sources that are of interest to the scientific community, such as gravitational wave sources and electromagnetic transients (e.g., gamma-ray bursts and stripped SNe), as all of them will have one or more XRB phases in their evolution (Zezas et al. 2019).

In this work we present, for the first time, population synthesis models of HMXBs that employ detailed stellar structure and binary evolution calculations throughout the entire evolution of the binaries. We carry out a parameter study of physical assumptions that affect the observed HMXB population. The different physical assumptions that we investigate iterate through various prescriptions, dictating aspects of evolution, such as SN kicks, the CE phase, and wind-fed accretion. In Sect. 2 we discuss the population synthesis code used throughout this work (POSYDON), the physics employed for the simulations, and the different model parameters we investigate. In Sect. 3 we construct XLFs from the simulated populations and investigate the main effects of the different model parameters. Section 4 discusses the comparison of different models to the observed XLF from Lehmer et al. (2019), the combined effects of certain parameters on the simulated XLF, and some possible reasons for any disparity between simulations and observations. Finally, we present our concluding remarks in Sect. 5.

2. Method and physical assumptions

For our study, we used the POpulation SYnthesis with Detailed binary-evolution simulatiONs code (POSYDON; see Fragos et al. 2023, for a detailed description of the code), a new binary population simulation framework that combines the flexibility of parametric binary population synthesis codes with detailed stellar structure and binary evolution models. Crucial to our study, POSYDON enables a physically accurate and self-consistent determination of the stability and mass-transfer rate evolution of mass-transfer phases thanks to extensive precomputed grids of binary-star models (using the Modules for Experiments in Stellar Astrophysics code or MESA code; Paxton et al. 2011, 2013, 2015, 2018, 2019). Furthermore, the available information on the internal structure of binary components enables an improved treatment of evolutionary phases, such as the CE and the core collapse. Of special interest here is also the capability of POSYDON to track the angular momentum content of both stellar components throughout the evolution of a binary, which allows us to identify potential Be XRBs, as we describe further below.

2.1. Initial binary properties

The target of this study is to construct synthetic XLFs of XRB populations dominated by HMXBs and compare them to observed HMXB XLFs from Lehmer et al. (2019). High-mass XRBs are associated with young (≲100 Myr) stellar populations, it could be argued that there would be observational effects from the wide range of ages of the consisting populations that would interfere with the comparison. While the SFH is expected to play a significant role in the XRB content of a galaxy, XLFs computed by combining observations from multiple galaxies are averaging out the effects from the varying ages in the populations. The observed XRB sample is composed of 38 galaxies, the majority of which have masses in the same order of the mass as the Milky Way (log(M/M_⊙)≈10.80; Licquia & Newman 2015). These galaxies are expected to have nearly uniform SFHs, and consequently, we expect the cumulative SFH of the sample to be nearly constant. Now, the star formation indicators that were used to measure the SFR of the galaxies in the sample, have some implicit assumptions about the SFH of the respective populations, which was a constant SFR for a duration of 100 Myr. To compare our synthetic XLFs to the observed one in the most self-consistent way, we assumed the same SFH. Thus, using POSYDON, we simulate synthetic binary populations corresponding to continuous star formation for 100 Myr. All the populations are run at solar metallicity (0.0142; Asplund et al. 2009).

The primary star is the initially more massive binary component at the zero-age main sequence (MS) and is sampled from the initial mass function from Kroupa (2001) within the range 7.0 − 120.0 M_⊙. The secondary’s mass is calculated from the mass ratio, which is sampled from a flat distribution (Sana et al. 2013). We assumed independent distributions of the mass ratio and period; however, recent studies show that there might be a dependence of mass ratios on the orbital periods (Moe & Di Stefano 2017). The initial orbital period is drawn from a power law in log-space (Sana et al. 2013) with the same limits as the POSYDON grids (Fragos et al. 2023) and the initial orbits are taken as circular since the current POSYDON version only follows the mass-transfer phases for circular binaries. The binary fraction assumed is 0.7 (Sana et al. 2012). We normalize the generated population by f_corr(≈5.89) to compensate for the fact that our initial parameters do not cover the entire existing masses in the Universe (see Appendix A in Bavera et al. 2020 for the detailed calculation of f_corr).

2.2. X-ray luminosity calculation

We selected binaries that consist of a CO and a non-degenerate star and identified as XRBs those that have a calculated bolometric X-ray luminosity greater than 10³⁵ erg s⁻¹. For XRBs undergoing RLO or wind-fed mass accretion, the bolometric X-ray luminosity () for sub-Eddington mass-transfer rates is calculated as follows,

(1)

where Ṁ_acc is the mass-accretion rate, η is the radiative efficiency of accretion, c is the speed of light, and ṁ ≡ Ṁ_tr/Ṁ_Edd is the Eddington ratio. The mass-transfer rate, Ṁ_tr, from the donor to the vicinity of the accretor is assumed to be equal to the accretion rate for ṁ ≤ 1.0. The mass-transfer rate, for XRBs undergoing accretion from stellar wind (not including Be XRBs that are treated separately), is calculated by following the Bondi & Hoyle (1944) accretion model, while for XRBs undergoing RLO, it is computed self-consistently in our grids of detailed binary evolution simulations. The radiative efficiency, η, is the fraction of the rest mass energy of accreted matter that is radiated away and is defined as

(2)

where M_acc is the CO mass, R_acc is the CO radius (for NS accretors) or the spin-dependent innermost stable circular orbit (from Podsiadlowski et al. 2003, for BH accretors), and G is the gravitational constant.

We use the accretion disk model by Shakura & Sunyaev (1973) to define the luminosity from super-Eddington mass-transfer rates. This model describes the X-ray luminosity for an accretion disk receiving material at super-Eddington rates, with strong outflows keeping the disk locally Eddington limited. For very high mass-transfer rates, Ṁ_tr > 8.5 × Ṁ_Edd; Ṁ_Edd being the Eddington limit, King et al. (2001) and King (2009) introduced a beaming effect, due to the presence of a geometrically thick accretion disk, in the form of a beaming factor to account for the collimation of outgoing emission. Therefore, X-ray luminosity from a super-Eddington accretion disk (that includes geometrical beaming) is described as follows,

(3)

where b is the beaming factor describing the amount of collimation to the outgoing radiation due to geometrical beaming from the accretion disk. The approximate value of b following King (2009) is defined as

(4)

The bolometric luminosity calculated by Eq. (3) would be perceived as super-Eddington if assumed to be isotropic. Since the model described above is only approximate, very high mass-transfer rates may provide unphysical strong beaming (b ≪ 10⁻³). To avoid this, we set a lower limit for b at 3.2 × 10⁻³, which roughly corresponds to a ṁ ≈ 150.0 (see also Lasota et al. 2016; Wiktorowicz et al. 2017).

Using similar criteria as Belczynski & Ziolkowski (2009) and Zuo et al. (2014), we identified Be XRBs by selecting binaries with fast spinning (surface velocity ≳70% of the critical surface velocity) MS donors, orbital periods in the range 10 − 300 days, and in a detached state (not undergoing RLO). We only included binaries with donor masses ≳6.0 M_⊙, which is a conservative estimate for the minimum mass of B stars (Hohle et al. 2010). An additional criterion that we used to identify Be XRBs is the donor decretion disk radius exceeding the Roche-lobe radius of the Be star, at which point the CO would start to accrete matter from the decretion disk. The Roche-lobe radius of the Be stars is calculated at the periastron of their eccentric orbits. Klement et al. (2017) estimated the outer edges of the decretion disk at around stellar radii of 30 to 150 R_*, following which we use 100 R_* as an approximate decretion disk radius (R_* being the donor radius). The X-ray luminosities of Be XRBs are calculated following Dai et al. (2006, Eq. (11)), which describes the correlation between peak observed luminosities for Be XRBs and the orbital periods (observational data from Raguzova & Popov 2005),

(5)

As Be XRBs are transient in nature, we included a duty cycle of 10% in our calculations (Sidoli & Paizis 2019, however see also Reig 2011 where duty cycles up to 30% are observed).

In order to compare the calculated bolometric luminosities () to X-ray luminosities that are in the Chandra band 0.5 to 8 keV (L_X), they are normalized using a bolometric normalization of L_{X, [0.5 − 8 KeV]}/L_bolometric ≈ 0.5, which is a zeroth order approximation of the Chandra bolometric correction estimates from Fragos et al. (2013a) and Anastasopoulou et al. (2022). The following subsections, Sects. 2.3–2.7, describe all the considered parameters (summarized in Table 1).

2.3. Remnant mass prescription

At the end of its evolution, a massive star creates an iron core that collapses resulting in the formation of a NS or a BH, depending on the properties of the pre-collapse star. The grids of detailed single- and binary-star models in POSYDON follow the stars’ evolution until their core-carbon exhaustion, beyond which we employ prescriptions that dictate the core-collapse outcomes. We used the following two prescriptions to determine the “explodability” and final remnant properties of the pre-SN stars.

First, the prescription from Patton & Sukhbold (2020) determines the final CO by treating the carbon-oxygen core of an evolved star separately from the rest of the star. Patton & Sukhbold (2020) calculated the explodability of stripped carbon-oxygen stars, for a range of masses and compositions, and derived a prescription for the outcome of the core collapse based on the final carbon-oxygen core mass and the composition of the carbon-oxygen core at carbon ignition. The masses of the newly formed COs are taken as the pre-SN He-core masses of the stars (accounting for any fall back, further described in the next section), assuming loss of any envelope during previous evolution.

Second, Fryer et al. (2012) predicted the core-collapse outcome based on the pre-SN carbon-oxygen core mass. Following their prescription, the newly formed CO masses are taken as the pre-SN masses of the stars, after taking account of any fall back. The authors provided “rapid” and “delayed” mechanisms for core collapse. The rapid mechanism predicts the presence of a mass gap in remnant masses in the range 2.0 − 5.0 M_⊙. The delayed mechanism produces a continuous range of NS and BH masses. Due to observations of COs that fall within the aforementioned BH mass gap (Thompson et al. 2019; Abbott et al. 2021), we focus on the delayed mechanism in our study.

2.4. Asymmetric natal kick normalization

The collapse of the iron core is accompanied by a SN (unless the core implodes without any explosion taking place). If the SN explosion is asymmetric, due to either asymmetric mass loss (Janka & Mueller 1994; Burrows & Hayes 1996; Janka 2017) or neutrino loss (Bisnovatyi-Kogan 1993; Janka 2017), the newly formed CO experiences a natal kick. We assume the magnitude of the kicks imparted are drawn from a Maxwellian distribution with the velocity dispersion σ = 265.0 km s⁻¹, for core-collapse SNe (from observed velocities of radio pulsars; Hobbs et al. 2005), or σ = 20.0 km s⁻¹, for electron-capture SNe (Giacobbo & Mapelli 2019). However, the true distribution of natal kicks is still uncertain (Verbunt et al. 2017; Igoshev et al. 2021; Willcox et al. 2021; Kapil et al. 2023).

Furthermore, we normalized the velocities based on the properties of the CO formed to calculate the final kick velocity, v_kick,

(6)

where n is the normalization and v_mag is the velocity magnitude drawn from the Maxwellian distribution. The distribution of the kick directions is isotropic. We study three prescriptions for natal asymmetric kicks to calculate the kick normalization (n): (i) The natal kicks are normalized by the mass of the created BH using the following normalization,

(7)

where 1.4 M_⊙ represents the mass of a typical NS and M_BH is the BH mass. With increasing BH masses, the imparted kicks are weaker. The NS kicks are not rescaled. (ii) The natal kicks are normalized by the fall-back fraction using the following normalization,

(8)

where f_fb is the fall-back fraction, which is the fraction of the stellar material that is not driven away by the SN explosion and falls back on the newly formed CO. Following the Patton & Sukhbold (2020) mechanism, f_fb = 1.0 for all BHs. Following the Fryer et al. (2012) mechanism, f_fb = 1.0 for stellar core masses (of the pre-collapse star) greater than 11 M_⊙ and for masses lower than 11.0 M_⊙, f_fb < 1.0. Hence, whenever f_fb = 1.0, there is no baryonic mass loss considered, while still accounting for neutrino losses (which is capped at 0.5 M_⊙Fragos et al. 2023). For NSs, when the SN mechanism follows the prescription from Fryer et al. (2012), the fall-back fraction for NSs is f_fb < 1.0 (however, for NSs from electron-capture SN, f_fb = 0). When the SN mechanism follows the prescription from Patton & Sukhbold (2020), it is assumed that there is no fall-back on the NS (f_fb = 0.0; Fragos et al. 2023). (iii) The natal kicks are not normalized for BHs (f_fb = 0, therefore, n = 1). Hence, the BHs receive strong kicks. Neutron star kicks are normalized using the fall-back fraction.

2.5. Circularization of the orbit at the onset of Roche-lobe overflow

We check for potential initiation of the RLO in a binary at the periastron of the orbit. At the onset of RLO, we assume that the orbit is circularized as tidal forces circularize the orbit at short timescales (about 10⁴ years for a semidetached binary with a giant donor; Zahn & Bouchet 1989; Verbunt & Phinney 1995). Following the circularization, POSYDON evolves the mass-transfer phase utilizing the precomputed model grids of circular binaries. The separation of the circularized binary will determine the stability and duration of the mass-transfer phase, and thus its observability as an XRB, thereby affecting the XLF. We investigated two options for the circularized binary separation. First, the final separation of the circular binary was taken to be the “periastron” of the eccentric orbit at the onset of RLO. The orbit will lose both energy and angular momentum as it circularizes and the new circularized binary separation (a_circ) will be

(9)

where a is the semimajor axis and e is the eccentricity of the eccentric orbit.

Second, the binary conserves its angular momentum and loses only energy as it circularizes, and the new circularized orbital separation is larger than the periastron distance of the originally eccentric orbit. For an eccentric binary with a semimajor axis a and an eccentricity e, the new circular binary separation a_circ) is calculated imposing the circular angular momentum is equal to the eccentric angular momentum:

(10)

In practice, both of these approaches are approximations of a potentially much more complicated process, encompassing the range of outcomes of the detailed treatment, and thus will give us some insight into the potential signatures that this process imprints on the XLF. In order to effectively check the validity of these assumptions, we would need a better treatment of eccentric orbits in mass-transferring binaries, which we do not currently have. A comprehensive theoretical framework of the secular evolution of eccentric, mass-transferring binaries has been recently developed (Sepinsky et al. 2007, 2009, 2010; Dosopoulou & Kalogera 2016a,b) but has not been applied yet in the context of detailed binary evolution models.

2.6. Common-envelope evolution

During a mass-transfer phase, if a dynamical instability is encountered, a CE phase is initiated, where the companion of the expanding star is engulfed by its envelope. The star that is spiraling inside the envelope loses orbital energy as frictional heat is deposited into the envelope, which increases the internal temperature of the envelope and it expands. When enough energy is deposited into the envelope during the inspiral, the envelope is expelled (refer to Ivanova et al. 2013, for a review of the CE). We parameterize CE using the so-called α_CE − λ_CE prescription (Webbink 1984; Livio & Soker 1988) to predict the outcome of the CE phase. This prescription calculates the change in the orbital energies before and after the CE () and compares it to the binding energy (E_bind) of the envelope. It follows as

(11)

(12)

where α_CE (known as the CE efficiency) is the fraction of orbital energy that is injected into the envelope, λ_CE is a factor that accounts for the structure of the envelope of the mass-losing star (de Kool 1990) known as the binding energy parameter, M₁ and M₂ are the masses of the two stars (with the subscripts “env” and “c” denoting the envelope and core mass of the respective star). If the predicted final orbital separation after the envelope ejection a_f is such that either the stellar core of the donor or the accretor is overfilling their respective Roche lobes, we consider the envelope ejection to be unsuccessful assuming the two bodies would merge. If neither is filling their Roche lobe, we consider the envelope ejection to be successful resulting in a detached, circular, and close binary.

For a high CE efficiency parameter (α_CE ∼ 1.0) since most of the available orbital energy can be used to unbind the CE, the envelope can be ejected with the binary shrinking less than for lower efficiency. Ultimately, an efficient CE leads to relatively wide post-CE orbits and allows lower mass companions to eject the envelope of the CO progenitor star.

We investigate the following two values for α_CE. First, to investigate the case when the energy conversion is not perfectly efficient, we study a population with α_CE = 0.3 to see the effect of a low efficiency on the CE phase, as some studies show α_CE in the range 0.1−0.5 (Zorotovic et al. 2010; Davis et al. 2010; Nandez et al. 2015). Second, to investigate the effect of full availability of orbital energy to eject the envelope, we used α_CE = 1.0.

A value of α_CE > 1.0 might point to extra sources of energy to eject the envelope due to uncertainties in the CE energy budget. Equivalently, it might also imply that the binding energy of the envelope is lower than that estimated from single-star models, as the stellar structure will differ for stars undergoing CE compared to undisturbed single stars (Fragos et al. 2019; Klencki et al. 2021; Marchant et al. 2021; Renzo et al. 2023). In our study, we adopt the second interpretation and explore different estimates of the binding energy further below. Additionally, for CE donors in the Hertzsprung gap, there is no well-formed, steep boundary between the core and envelope. In this case, the inspiral never stops and the CE results in a binary merger (Belczynski et al. 2007). This is the “pessimistic” approach when estimating the outcome of the CE (Taam & Sandquist 2000; Ivanova & Taam 2004).

The binding energy (E_bind) of a stellar envelope is often characterized by the binding energy parameter (λ_CE), which in turn is dependent on how the He-rich core of the donor is defined. To begin with, λ_CE is defined as follows,

(13)

where M and M_env are the stellar mass and envelope mass of the donor, and R is the stellar radius. The envelope binding energy, and hence also the λ_CE, are dependent on where we set the core-envelope boundary in the mass profile of the donor star and its value changes as the star evolves (de Kool 1990). In the context of the CE phase, the core-envelope boundary separates the part of the star that will be ejected, the envelope, from the remaining stripped core, which will contract and allow the binary to detach. The position of the boundary is often parameterized as the first point in mass where the fraction of H drops below an adopted threshold. The core-envelope boundary is not a well-constrained value, depending highly on the evolutionary stage of the donor and investigations show that the choice of the boundary will have different results on the CE outcome (Dewi & Tauris 2000; Tauris & Dewi 2001; Ivanova 2011; Kruckow et al. 2016). While traditionally this boundary is set where the hydrogen (H) abundance (X_H) drops at least below 10%, recent studies of COs inspiraling within the envelope of a supergiant star suggested that the core-envelope boundary might be located further out (Fragos et al. 2019; Quast et al. 2019; Klencki et al. 2021; Marchant et al. 2021). The binding energy of the envelope will be higher (or lower) if we set the boundary at a deeper (or outer) mass coordinate. Therefore, by changing the boundary we implicitly change the λ_CE, and the values for the boundary that we study are at H-abundance of 0.01 and 0.3. The use of stellar profiles from detailed grids in POSYDON aids in a physically motivated estimation of λ_CE.

2.7. Criteria for observable wind-fed accretion

For wind-fed HMXBs, X-ray emission is produced when the accretor captures and accretes a fraction of the donor’s stellar wind. One can estimate the amount of accreted material using the Bondi-Hoyle approximation, which assumes that the accretor moves through the locally isotropic wind of the mass-losing companion (Bondi & Hoyle 1944). In this mechanism, the amount of accreted material is proportional to the rate at which the companion star is losing mass in stellar winds, inversely proportional to the orbital separation squared, and also inversely proportional to the fourth power of the wind velocity. This makes it more efficient for binaries with giant-star companions, which tend to have sufficiently intense and slow stellar winds.

In contrast to NS accretors, which have a hard surface, material accreting onto a BH has to first form an accretion disk before falling into the BH’s horizon, in order to efficiently dissipate its gravitational energy in the form of electromagnetic radiation and become a bright X-ray source. Hirai & Mandel (2021) investigated the BH HMXB Cygnus X-1, which does not follow the simplified picture presented by the Bondi-Hoyle model, as it has a focused accretion stream (Miller et al. 2005; Hanke et al. 2009). The accretion stream is produced due to the high Roche-lobe filling factor of the companion and is aided by the rotation of the donor, which reduces wind terminal velocity and promotes the capturing of the wind leaving the equator.

When this stream has high enough angular momentum, material forms an accretion disk before falling into the BH, producing X-rays in the process. Hirai & Mandel (2021) estimated that for an HMXB to have observable X-ray emission from wind-fed accretion, the donor should fill its Roche lobe by at least 80 to 90%. We identify observable wind-fed XRBs as XRBs (having MS donors) where the donor fills its Roche lobe by at least 80%. In our population study, we explore two cases. First, we used the Hirai & Mandel (2021) criteria and identified wind-fed XRBs where the donor fills its Roche lobe by at least 80%. Second, we did not use this criterion for defining the creation of an accretion stream.

3. Results

Each of the parameters explored (Table 1) probes different aspects of binary evolution and might have an effect on the resulting population of simulated HMXBs, and thus the synthetic XLF. We present all the 96 constructed synthetic XLFs (cumulative and differential) in Appendix A. To investigate these effects, we compared the XLFs from the simulated populations to a “default” XLF that follows the parameters shown in Table 2 and corresponds to model 15 in Appendix A. In constructing the synthetic XLFs we need to take into account any directional X-ray emission. For instance, Eq. (4) describes the geometrical beaming of the X-ray emission in the case of super-Eddington mass transfer, which implies that some XRBs would be less visible due to non-preferential line-of-sight angles with respect to the rotational axis of the accretion disk. Therefore, in constructing the XLF we down-weight the beamed sources using the beaming factor described in Eq. (4). We down-weighted the Be-XRB sources in accordance with a duty cycle of 10%. Additionally, we approximated all calculated bolometric X-ray luminosities to X-ray luminosities in the Chandra band in the range 0.5−8 keV.

The resulting synthetic XLF of the default population (shown as N(> L_{X, [0.5 − 8 Kev]})/SFR) split to show different types of mass transfer, COs, and donor types, is shown in Fig. 1. The figure also compares the synthetic XLF to the observed HMXB XLF from Lehmer et al. (2019). The assumptions we made in the simulations regarding the SFH and metallicity are similar to the measured properties of the host galaxies of the XRBs that were used to construct the observed XLF (see the discussion in Sect. 2.1). Around the luminosity of 10³⁸ erg s⁻¹, we see an overabundance of XRBs in the synthetic XLF compared to observations, by a factor of up to 10. However, the synthetic XLF shows a break at an approximate luminosity of 10³⁸ erg s⁻¹ where the shape of the XLF changes, a feature also present in the observed XLF. For X-ray luminosities below the break, Be XRBs dominate the population (dotted gray line) and for luminosities above the break, it is dominated by RLO BH XRBs with H-rich donors (shown as the solid brown line). Among the Be XRBs, ∼96% of them have NS accretors, while only ∼4% appear to be BH Be XRBs (that is, a ratio of NS Be XRBs to BH Be XRBs of 24). This model prediction is consistent with other theoretical models that found the ratio of NS Be XRBs to BH Be XRBs in the range 10−50 (Belczynski & Ziolkowski 2009; Shao & Li 2014), depending on the specifications of the physical model assumed.

Fig. 1. XLF of the default population (solid blue line) compared to the observed XLF from Lehmer et al. (2019, pink shaded region). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.

There is evidence that Be XRBs are formed from binary interactions that occurred previously in the evolution, with the B-type stars spinning up due to stable transfer of matter from the NS-progenitors (Shao & Li 2014; Vinciguerra et al. 2020). The number of NS Be XRBs is most likely affected by the slope of the initial-mass function; a high number of lower-mass stars produced would lead to a higher number of lower-mass COs (Belczynski & Ziolkowski 2009). Another effect that could dictate the Be-XRBs demographic is the CE phase, as binaries forming BHs have more massive primaries that initiate CE earlier (often when they are in the Hertzsprung gap) than binaries with lower-mass primaries, which merge due to the application of the pessimistic approach during the CE phase (see Belczynski & Ziolkowski 2009; Dominik et al. 2012). However, in our models, when we compare against a population run with the optimistic approach (where we allow all donor types to successfully eject the CE), we see no significant differences in XRB demographics. The total default XLF has a steeper slope above the break, which is qualitatively similar to observations, primarily dominated by the population of RLO BH XRBs with H-rich donors. As an additional check, we run a population with an initial binary fraction of 0.5 (that would affect the normalization) and see no significant difference.

Looking at the combined XLF in Fig. 1 for X-ray luminosities about 10³⁸ to 10³⁹ erg s⁻¹, the number of XRBs exceeds the observed sample by factors up to ∼10. We refer to this excess as the “XLF bump” in the rest of this study. Looking at the subpopulations, the bump seems to be primarily coming from the RLO BH XRBs with H-rich donors. At lower luminosities (below 10³⁸ erg s⁻¹) XRBs are dominated mainly by Be XRBs, which are highly dependent on our assumptions for Be XRBs, with some contribution from wind BH XRBs with H-rich donors. There is also a slight difference at the higher luminosity end of the XLF with respect to the observations, a result of the scarce observational data at these luminosities. The observed XRBs reach up to 3 × 10⁴⁰ erg s⁻¹ while out simulations produce XRBs with luminosities up to 3 × 10⁴² erg s⁻¹. Since at a given luminosity, the cumulative XLF carries the uncertainties of the following luminosities, we refer to Fig. 2, which shows the differential form of the XLF. We see that the populations of Be XRBs, RLO BH XRBs, and wind-fed BH XRBs (both with H-rich donors) dominate the XLF bump. Since we do not model part of the binary populations (initial zero-age MS binary consists of mainly massive stars), luminosities below ≲10³⁶ erg s⁻¹ are under-filled.

Fig. 2. Differential form of the synthetic XLFs, showing the default population (solid blue line) compared to the observed XLF from Lehmer et al. (2019, pink shaded region). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.

Additionally, we can look at the distribution of the XRB luminosities with respect to the binary separation in Fig. 3, distinguishing between massive donors (≳8.0 M_⊙) and less massive donors (< 8.0 M_⊙). We see that the different parts of the parameter space are occupied by different types of XRBs, which reflects the dependence of their luminosities on their respective formation channels. Generally, the XRBs undergoing RLO with H-rich donors have close orbits centered around 10 days and bright luminosities (≳5 × 10³⁷ erg s⁻¹), while wind-fed XRBs with H-rich donors have wider orbits going as high as 10⁵ R_⊙ and lower luminosities (≲5 × 10³⁸ erg s⁻¹). XRBs with He-rich donors are present at close orbits (≲10 R_⊙) because He stars are more compact compared to H-rich stars, which enables the binaries to evolve to close orbits.

Fig. 3. Distribution of the calculated X-ray luminosities of the simulated XRBs in the default population with respect to the orbital separation. The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with gray circles. Large and small symbols denote massive (with ≳8.0 M⊙) and less massive (with < 8.0 M⊙) XRBs, respectively.

In Fig. 3, the less massive donors in the population of RLO BH XRBs are at the lower end of the distribution, going as bright as 10³⁹ erg s⁻¹, as the massive donors can easily undergo high mass-transfer rates without exceeding the super-Eddington limit (or exceeding it by a relatively small amount); this makes their RLO phase more stable at higher luminosities. In the population with wind-fed BH XRBs, the less massive donors occupy the entire luminosity range. However, they are mostly present around large separations of around 10³ R_⊙, with the donors being mostly in the post-MS, giant phase, during which they can drive strong winds. At closer orbits (≲10³ R_⊙), the systems undergoing wind-fed accretion are mostly massive MS stars that have strong stellar winds. The estimated X-ray luminosities for Be XRBs follow Eq. (5), which expresses the observed correlation between the peak luminosity and orbital period, hence the tight correlation in Fig. 3.

In the following sections, we explore how different assumptions (mentioned in Sect. 2 and summarized in Table 1) about physical processes in the formation and evolution of XRBs may leave an imprint on the XLF of the whole population. We also explore the best-fitting combination of parameters that reduces the XLF bump.

3.1. Effect of varying the remnant mass prescription

We explored two prescriptions that describe the masses of SN remnants, from Patton & Sukhbold (2020) and Fryer et al. (2012; see our Sect. 2.3). We compare the synthetic XLFs of XRB populations simulated with two remnant formation prescriptions in Fig. 4, one following Patton & Sukhbold (2020) (the default population) and the other using the delayed prescription from Fryer et al. (2012; which corresponds to model 7 in Appendix A). We find no significant differences in the resulting XLFs between the two mechanisms. In order to see how these prescriptions affect the resulting populations we look into the mass distribution of the COs produced.

Fig. 4. XLFs of the default population with the parameters in Table 2 showing the Patton & Sukhbold (2020) SN mechanism (blue line) and the population with the Fryer et al. (2012, delayed) SN mechanism (orange line).

Figure 5 shows the normalized distributions of CO masses (by type) at the end of the SN phase, for binaries that remain bound post-SN, for both prescriptions. The CO mass range covered by the two prescriptions is similar. However, Patton & Sukhbold (2020) produce a gap in the CO mass distribution between NSs and BHs and a narrow NS mass distribution, while the delayed prescription produces a continuous CO mass distribution. In both cases, the NS mass distribution peaks around 1.4 M_⊙ (similar to Özel et al. 2012; Özel & Freire 2016). Our results are also similar to the more detailed comparative study by Patton et al. (2022), and they indicate that the HMXB XLF does not provide sufficient information to constrain the mass spectrum of COs.

Fig. 5. Distribution of newly formed CO masses (MCO) for binaries that remain bound post-SN, for the default population with the parameters in Table 2, including the Patton & Sukhbold (2020) SN mechanism (blue line) and the population with the Fryer et al. (2012, delayed) SN mechanism (orange line). The populations are further split by the type of CO, NS (solid line), and BH (dashed line). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the default population for NS formation (NS: Patton & Sukhbold 2020, solid blue line) to 1.

3.2. Effect of varying the natal kick normalization

As detailed in Sect. 2.4, we investigated three different prescriptions for the SN kick normalizations (see Table 1). The kick normalization will influence the future evolution of the first-born CO and change the fraction of binaries that disrupt or have significantly altered orbits. These effects will reflect in the synthetic XLFs. In the default BH-XRB population with mass-weighted kicks, the strength of the kicks decreases with increasing BH mass, implying that low-mass BHs get considerable kicks. For BH kicks normalized by the fall-back fraction (corresponding to model 31 in Appendix A), the BHs either receive no kicks or weak kicks (depending on the SN mechanism and the progenitor carbon-oxygen mass). Out of the three prescriptions described in Sect. 2.4, the strongest BH kicks are for the case with no additional normalization (corresponding to model 47 in Appendix A), where BHs receive kicks drawn directly from a Maxwellian distribution with σ = 265 km s⁻¹, the same as observationally constrained for NSs (Hobbs et al. 2005).

Figure 6 shows the comparison of these three normalizations with respect to our default model. With increasing kick velocities, the number of XRBs below an X-ray luminosity of 3 × 10³⁸ erg s⁻¹ decreases. At the same time, an excess of XRBs emerges at X-ray luminosities in the range 3 × 10³⁸ to 10³⁹ erg s⁻¹ for populations with BH mass normalized kicks (that is, the default model) and with kicks that are not normalized. Above 10³⁹ erg s⁻¹, the XLF slope for the population with no normalization on the kicks drops below the other two normalization prescriptions. Figure 7 shows again the XLFs of the three models, now further split by the type of mass transfer the XRBs currently undergo (wind or RLO). The first striking difference is the lower number of the wind-fed XRBs for the population with no kick normalization shown by the green-dashed line (reduced by a factor of 20 compared to the default population). Another difference is the lower number of RLO XRBs for the population with fall-back normalized kicks. The reasons for these differences can be understood on the basis of the XRB orbital period distribution right after the SN that forms the BH.

Fig. 6. XLFs of the default population with the parameters in Table 2 showing the population with BH mass normalized SN kicks (blue line), the population with fall-back normalized SN kicks (orange line), and the population with no normalization on the SN kicks (green line).

One complication here is that the default remnant mass prescription, which follows Patton & Sukhbold (2020), assigns f_fb = 1.0 for all BHs, resulting in no BH kicks. For the prescription from Fryer et al. (2012), low-mass BHs receive weak kicks. Figure 8 shows the XLFs of the three BH kick models, further split by the type of mass transfer. Comparing to Fig. 7 we see mostly similarities, except the higher number of RLO XRBs for the population with no kick normalization (by a factor of ∼2). However, the luminosity regime in which this effect takes place (below 10³⁹ erg s⁻¹) is dominated by wind XRBs, and hence this increase is not reflected in the simulated XLF. Generally, we focus more on the default remnant prescription when discussing BH kick normalizations.

Fig. 7. XLFs of the default population with the parameters in Table 2 showing the population with BH mass normalized SN kicks (blue lines), the population with fall-back normalized SN kicks (orange lines), and the population with no normalization on the SN kicks (green lines). The populations follow the remnant mass prescription from Patton & Sukhbold (2020). The populations are further split by the type of mass transfer that occurred, RLO (solid line) or wind-fed accretion (dashed line). For the sake of clarity, the Be-XRB population is not shown as it is not affected by the kick normalization and obscures the effects of the BH kick normalization.

Fig. 8. XLFs of the default population with the parameters in Table 2 showing the population with BH mass normalized SN kicks (blue lines), the population with fall-back normalized SN kicks (orange lines), and the population with no normalization on the SN kicks (green lines). The populations follow the remnant mass prescription from Fryer et al. (2012, delayed). The populations are further split by the type of mass transfer that occurred, RLO (solid line) or wind-fed accretion (dashed line). For the sake of clarity, we did not include Be XRBs in this figure as they are not affected by the kick normalization and obscure the effects of the BH kick normalization.

Since the kick prescriptions we use affect only the BH XRBs, we only focus on BH binaries in the following discussion. Figure 9 shows the distribution of the orbital periods for the three kick normalization prescriptions, right after the SN of the primary star for the progenitor binaries of BH XRBs. For BH mass-weighted kicks, low-mass BHs will receive a weak kick, inversely proportional to their mass, which will disrupt some of the binaries at large orbits with orbital velocities comparable to the kicks (≳100 days in Fig. 9). Binaries at periods ≳100 days are surviving with fall-back weighted kicks (refer to Sect. 2.4) as they have full fall-back and receive no kicks. The effect of strong kicks disrupting binaries at larger orbits is even more pronounced in BH XRBs with no kick normalization where most binaries above periods of 100 days are disrupted (in addition to an overall lower number of surviving binaries). The disruption of wider orbits leads to a lower number of wind-fed BH XRBs. This argument also explains the higher number of wind-fed XRBs for fall-back weighted kicks compared to those for mass-weighted kicks, and the lower number of XRBs with no kick normalization compared to the other two populations at high-tail ends of the XLF.

Fig. 9. Distributions of orbital periods (Porb) at the end of the SN for the progenitors of BH XRBs, for the default population with the parameters in Table 2, including the population with BH mass normalized SN kicks (solid blue), the population with fall-back normalized SN kicks (hatched orange), and the population with no normalization on the SN kicks (hatched green). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the default population (BH mass normalized kicks shown in solid blue) to 1.

Looking at Fig. 9 for orbital periods of less than 100 days, mass-weighted normalization and fall-back normalization (with no BH kicks, effectively) allow more BH binaries to survive the SN kick than un-normalized kicks. This occurs despite the mass-weighted and un-normalized kicks imparting stronger kicks on many BH binaries than the fall-back kicks, which are much weaker (for Fryer et al. 2012), or no kicks (for Patton & Sukhbold 2020). The reason for this is the eccentricity imparted on the binaries surviving the SN kicks. Figure 10 shows the distribution of the eccentricities of the BH XRBs that survive SNe. There, we see that the surviving binaries in mass-weighted kicks and un-normalized kicks have acquired increased eccentricities. When the secondary stars evolve in these binaries and fill their Roche lobe in a highly eccentric orbit, our default assumption of instantaneous circularization at the periastron distance leads to close circular RLO XRBs (refer to Eq. (9)). This causes the overall higher number of close-orbit RLO XRBs in the populations with mass-weighted and un-normalized kicks compared to the fall-back XLF (Figs. 7 and 9). This higher number of RLO XRBs, in combination with the suppression of wind-fed binaries, produces the bump in the XLF in the 3 × 10³⁸ to 10³⁹ erg s⁻¹ X-ray luminosity range (see Fig. 6) for the XRBs with BH mass weighted kicks and kicks with no normalization. Additionally, there is a small population of SN-surviving BH-binaries at periods ≲0.3 days (for all three normalization approaches) that arise from a double CE phase in the evolution prior to the BH formation, for binaries that started with an initial mass ratio ∼1.

Fig. 10. Distribution of the eccentricity (e) versus log10 of the orbital period (Porb) for the same populations as in Fig. 9.

One thing to be noted is that in the population with fall-back normalized BH kicks and the Patton & Sukhbold (2020) SN mechanism, since there is no baryonic mass loss during the formation of the BH, the change in eccentricity is provided by the neutrino loss (which is limited up to 0.5 M_⊙). Therefore, in the absence of kicks by mass loss, all eccentric systems would be assumed to result from neutrino losses. However, an additional effect might come into play as the resultant BH mass is taken as the mass of the pre-SN He-core, and whatever H-envelope is left on the star at this stage is assumed to be lost. If the star loses considerable mass in the form of its envelope, an eccentricity is introduced into the orbit. The difference between the pre- and post-SN masses of BH-resulting stars goes up to 8.0 M_⊙ for 13% of the BH-forming binaries (with a peak around 2.0 M_⊙), clearly more than that accounted for by neutrino mass losses. Hence, there donors retain their H-envelope during their pre-SN evolution, which was then lost during the SN event.

Therefore, in Fig. 10, we see some binaries with increased eccentricities (going up to 0.4) around the orbital period of 10⁴ days, in the population with fall-back normalized kicks. These systems end up as wide wind-fed XRBs in Fig. 7 and their eccentricity comes from the combined effect of H-envelope ejection and neutrino loss. The BH masses in these systems span a range the 4.00−14.00 M_⊙, and their pre-SN progenitors were late giants with radii ∼1000 R_⊙. The reason these systems are centered around 10⁴ days is that they form the tail end of the population with eccentric post-SN orbits. These systems did not undergo an RLO phase previously as they had wide orbital periods at zero-age MS ≲5 × 10³ days, which further widened with stellar winds and were on the order of 10⁴ days when the first SN occurred. Orbits wider than 10⁴ days at SN were disrupted due to the kick from the H-envelope loss, and orbits narrower than these had interacting binaries that resulted in loss of the H-envelope (the pre-SN stars in these binaries are in their giant phase with radii in the order of 1000 R_⊙). In general, Our findings highlight the importance of our assumptions for the circularization of eccentric orbit when RLO occurs. We discuss this further in Sect. 3.3.

3.3. Effect of varying the circularization of the orbit at the onset of Roche-lobe overflow

In Sect. 2.5, we mentioned that we can vary the assumptions for the circularization of the orbits at the onset of RLO. In this section we look at the effects of changes in orbit circularization on the XLF. Figure 11 shows the synthetic XLFs for models following the different assumptions, and an additional model that is explained further below. A circular orbit with the same angular momentum as the eccentric orbit at the start of RLO would have a larger separation than an orbit with the separation equal to the periastron distance (as is in the default population), as is evident on comparing Eqs. (9) and (10). The population where we conserve the orbital angular momentum during the circularization at the onset of RLO corresponds to model 16 in Appendix A. Intuitively, wider orbits should reduce the number of binaries undergoing RLO early on in the evolutionary life of the donor (in the long-lasting MS phase) and increase the binaries undergoing RLO with giant donors. Since the latter phase is a shorter-lived phase, the number of observed RLO XRBs would be reduced.

Fig. 11. XLFs of the default population with the parameters in Table 2 showing circularization at RLO onset at the periastron of the eccentric orbit (blue line), the population with the circularization at an orbit conserved angular momentum (orange line), and the population with circularization at the binary separation of the eccentric orbit (green line).

To investigate an even more extreme assumption of conserving the orbital energy during the circularization process, we run a single population where binaries that overfill their Roche lobe in eccentric orbits, result in circularized orbits with the same separations as the pre-RLO eccentric orbits (shown as the green line in Fig. 11). The number of XRBs around luminosities 10³⁸ to 10³⁹ erg s⁻¹ reduces by a factor of up to 2 with the wide orbits, a difference that, when considering the total XLF, is not pronounced. This approach is similar to the assumption made in some population synthesis codes, for example BSE (Hurley et al. 2000)¹. The rest of the parameters for this model are the same as the default population. We do not use this circularization option in the rest of the iterations of POSYDON models.

Figure 11 shows a slight reduction in the bump of the XLF around ∼10³⁸ erg s⁻¹ for the models where we assume that binaries overfilling their RLO in eccentric orbits instantaneously circularize, conserving either angular momentum or orbital energy. To further investigate how the treatment of the circularization affects the populations, in Fig. 12 we divide the XLFs of the three models by the type of mass transfer undergone by the XRBs (namely, wind-fed XRBs, RLO XRBs, or Be XRBs). There is a noticeable effect of increasing the circularized orbital separation in the RLO systems. Systems with wider orbits lead to fewer RLO XRBs.

Fig. 12. XLFs of the default population with the parameters in Table 2 showing the population with pre-RLO circularization at the periastron (blue lines), at an orbit with conserved angular momentum (orange lines), and at the binary separation of the pre-RLO eccentric orbit (green lines). The populations are further split by the type of mass transfer that occurred, RLO (solid line), and wind-fed accretion (dashed line). Be XRBs are shown with a dotted gray line, showing no effect as they are not circularized. AM stands for angular momentum.

We can see the distributions of the circularized orbital periods of the resulting XRBs at RLO onset, for all three approaches in Fig. 13. The presented systems were identified as initiating RLO at the periastron, and then the separation of the circularized orbit was adjusted following the respective model used. The wider circularized orbits suppress the number of XRBs below periods of ∼30 days. This results in the decrease of RLO XRBs seen in the synthetic XLFs (Fig. 12). On the other hand, looking at periods ≳100 days, circularization at the binary separation leads to a larger number of XRBs compared to the other two populations, as expected. About 90% of the wide orbits (with periods ≳100 days) lead to wind-fed XRBs, for all three circularization options.

Fig. 13. Distributions of orbital periods (Porb) of circularized orbits at the onset of RLO for the resulting XRBs. The presented systems were identified as initiating RLO at the periastron, and then the separation of the circularized orbit was adjusted following the respective model used. The figure shows the default population with the parameters in Table 2 including circularization at RLO onset at the periastron of the eccentric orbit (solid blue), the population with the circularization at an orbit conserved angular momentum (hatched orange), and the population with circularization at the binary separation of the eccentric orbit (hatched green). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the default population (with orbits circularized at the periastron shown in solid blue) to 1.

3.4. Common-envelope evolution

3.4.1. Effect of varying the common-envelope efficiency

We investigate the effect of different assumptions for the α_CE − λ_CE prescription on the synthetic XLF. Firstly, we look into two different values of α_CE. Figure 14 shows the synthetic XLFs of populations with α_CE = 1.0 (default population) and α_CE = 0.3 (corresponding to model 11 in Appendix A). Overall, there is no significant difference between the two populations when looking at the total XLF. This is similar to the results by Linden et al. (2011) and Zuo & Li (2014), who found that α_CE does not significantly alter the HMXB numbers. However, in their study, a high value of α_CE is needed to fit the XLF. Figure 15 shows the XLF of the population with α_CE = 0.3, split into the different types of XRBs. Compared to Fig. 1 with α_CE = 1.0, there is a decrease in XRBs with helium-rich (He-rich) donors, as those binaries have undergone a CE in their prior evolution and the α_CE governs this phase. For instance, the population of wind-fed BH XRBs with He-rich donors disappears with a lower α_CE (shown as the dashed purple line in Fig. 1) and the population of RLO NS XRBs with He-rich donors decreases by a factor of 3 (shown as the solid yellow line). The lower the α_CE, the more the accretor inspirals in the envelope in order to eject it. Alternatively, higher α_CE corresponds to a larger post-CE orbital separation and a higher probability of surviving the CE phase. However, the variation of α_CE has little effect on the total synthetic XLF since, based on our models, BH XRBs with hydrogen-rich (H-rich) donors and Be XRBs, which dominate the XLF, do not go through a CE phase during their formation.

Fig. 14. XLF of the default population with the parameters in Table 2 showing the population with CE efficiency (αCE) of 1.0 (blue line) and the population with αCE of 0.3 (orange line).

Fig. 15. XLF of the population with αCE = 0.3 (solid orange line). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.

3.4.2. Effect of varying the definition of the core-envelope boundary in CE evolution

The second CE parameter that we look into is the definition of the core-envelope boundary, via which we indirectly investigate the binding energy parameter λ_CE. Figure 16 compares the synthetic XLFs of the populations with the core-envelope boundary definition during the CE at 0.30 (default population) and at X_H = 0.01 (corresponding to model 13 in Appendix A; refer to Sect. 2.6 for more details). Similarly to varying the CE efficiency, we do not see a distinct difference between the two XLFs. We can investigate the populations in detail by looking at Fig. 17 where we show the XLF of the model with core-envelope boundary definition at X_H = 0.01, split by the different XRB subpopulations. XRBs with He-rich donors, which have undergone a CE in their prior evolution, increase in number when compared to the default model shown in Fig. 1.

Fig. 16. XLFs of the default population with the parameters in Table 2 showing the population with core-envelope definition at XH = 0.3 (blue line) and the population with core-envelope definition at XH = 0.01 (orange line).

Fig. 17. XLF of the population with core-envelope definition at XH = 0.01 (solid orange line). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.

Comparing the post-CE populations in Figs. 1 and 17, we get similar fractions of surviving NS and BH binaries for both definitions of the core-envelope boundary. The fraction of binaries with NSs and BHs that enter a CE and successfully eject it is ∼5% and ∼18%, respectively. However, some XRB subpopulations with He-rich donors show significant variations. Specifically, there is an increase in the number of RLO NS XRBs with He-rich donors (increase by a factor of 30) and in the number of wind-fed BH XRBs with He-rich donors (increase by a factor of 6), going from the core-envelope boundary definition during the CE at X_H = 0.3 (Fig. 1) to X_H = 0.01 (Fig. 17). Even then, this increase is not reflected in the overall synthetic XLFs as it is dominated by systems not undergoing CE. Additionally, the fractional increase in the two subpopulations is not the same, with RLO NS XRBs with He-donors increasing by a larger factor than wind-fed BH XRBs with He-rich donors.

The orbital energy released into the CE by the inspiraling CO depends on the mass of the CO (∝M_CO), with BHs providing more energy than NSs. Hence, more BHs survive the CE phase than NSs (cf. ∼18% versus ∼5%). The XRB population that results from the further evolution of binaries that survive the CE, is determined by the distribution of the orbital periods of the surviving binaries. Figure 18 shows the normalized distributions of the orbital periods post-CE for the models with the core-envelope boundary definitions at X_H = 0.30 and at 0.01, for different types of CO (NSs or BHs). We see that the distribution shifts toward lower periods when going from X_H = 0.3 to X_H = 0.01, with a greater difference for BH-He-star binaries than NS-He-star binaries (which maintains the peak of its distribution around ∼1 day). The resulting closer orbits affect different types of binaries differently.

Fig. 18. Distribution of the post-CE orbital periods (Porb) for the default population with the parameters in Table 2 including the core-envelope definition at XH = 0.30 (blue histogram) and the population with the core-envelope definition at XH = 0.01 (orange histogram). The populations are further split by the type of CO in the surviving binary, showing NSs (step) and BHs (hatched). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the NS-He post-CE binary population to 1.

For NS binaries, closer orbits (for the model with core-envelope boundary definition during the CE at X_H = 0.01) result in a higher number of RLO systems, as seen in Fig. 17. The reason is that the subsequent donors for these systems typically have masses ≲4 M_⊙ and will expand overfilling their Roche lobes earlier in their evolutionary lifetime. This leads to more donors in their He-MS phase than in their He-giant phase. Mass-transfer phases with He-MS donors have a longer duration than with He-giant donors and therefore, there would be a significantly higher chance of observing them. Within our Milky Way, so far, there are no observations of RLO NS XRBs with He-rich donors. If they are present they would be expected to lie along the Galactic plane where younger populations dominate. However, the extinction along the Galactic plane would reduce any chances of observations. Since our models already predict a small number of these sources within our galaxy, even in the most optimistic case (≲10), the lack of observed systems cannot put a robust constraint on our models.

For BH binaries, closer orbits (for the model with core-envelope boundary definition at X_H = 0.01) result in more efficient Bondi-Hoyle wind-fed accretion and hence a higher number of bright wind-fed BH XRBs with He-star donors. However, there are no RLO BH XRBs with He-rich donors as seen in Figs. 1 and 17. The reason is that the stripped He-rich stars in these post-CE binaries have masses ≳4 M_⊙, and do not expand much (to fill their Roche lobe) during their subsequent evolution (refer to Fig. 8 in Fragos et al. 2023; effect also seen by Habets 1986). Hence, no RLO BH XRBs with He-rich donors are produced in our synthetic populations for either core-envelope boundary definitions.

3.5. Effect of the criteria for observable wind-fed accretion

Figure 19 presents the effect of limiting the observable X-ray in wind-fed HMXBs, following the study by Hirai & Mandel (2021, see Sect. 2.7 for more information). We limit the subpopulation of wind-fed BH HMXBs with the criterion that filling of the donor Roche lobe by at least 80% is needed to have observable X-ray luminosity. This criterion is applied only for systems with BH accretors undergoing wind-fed accretion, as the hard surface of a NS lifts the requirement of the accretion disk formation in order to produce X-ray emission. Since wind-fed BH XRBs dominate the XLF below luminosities ∼10³⁸ erg s⁻¹, the Hirai & Mandel (2021) criteria only affect the number of BH XRBs at the lower luminosity end of the XLF. Looking at the combined XLFs in Fig. 19, XRBs with X-ray luminosities below ∼10³⁸ erg s⁻¹ are indeed reduced in number.

Fig. 19. XLFs of the default population with the parameters in Table 2 showing no criterion for wind-fed accretion (blue line) and the population with limited observability of wind-fed XRBs following Hirai & Mandel (2021, orange line).

The XLF of the model with the criterion by Hirai & Mandel (2021) is shown in Fig. 20 (and corresponding to model 63 in Appendix A), divided into the types of XRBs (showing the types of mass transfer, donor, and accretor). Compared to the default model (Fig. 1) we can see a clear decrease in the number of wind-fed BH XRBs. Wind-fed BH XRBs with H-rich donors decreased by a factor of 5, while wind-fed BH XRBs with He-rich donors decreased by a factor of 10. The majority of the population of Be XRBs is not affected by the limited observability as the condition is applicable to only BH XRBs.

Fig. 20. XLF of population with limited observability of wind-fed XRBs following Hirai & Mandel (2021). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS) and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.

4. Discussion

In addition to the physical assumptions investigated above, there are certain assumptions that have been implicitly made when carrying out the population synthesis study, the effect of which was not explored in detail. When the binary consists of two non-degenerate stars at the beginning of the evolution, any mass transfer that occurs is considered to be mostly nonconservative, as the spin-up of the accretor close to critical rotation prevents further accretion (Packet 1981; Petrovic et al. 2005; Ritchie et al. 2012; Renzo & Götberg 2021; Sen et al. 2022). The excess material is lost from the vicinity of the accretor, carrying away its angular momentum in the form of a fast isotropic wind. A more conservative mass-transfer phase would lead to more massive accretors with increased mass-transfer stability and consequently, more massive companions in the resulting XRBs. Additionally, the envelopes of these massive accretors would be easier to eject during a CE phase that might follow, increasing their chances of survival (Renzo et al. 2023). After either of the two stars has undergone a SN event and formed a CO, there might be a second RLO phase for which the mass-accretion rate is limited by the CO Eddington limit. The excess material, again, leaves the system taking away the angular momentum of the accretor. With a change in the critical mass-accretion rate, the stability of the mass transfer and the resulting population would be affected.

We assume stellar wind loss even during RLO, which would affect the orbit, particularly for massive donor stars with lower-mass accretors where the widening effect from strong winds would counter the contraction from RLO and stabilize the mass transfer (refer to Fragos et al. 2023, for a detailed description of the wind-loss prescriptions used). Additionally, there could be increased RLO stability due to irradiation effects from the accretion luminosity of the CO driving extra mass loss from the donor surface (Ruderman et al. 1989; Tavani & London 1993), which is not considered in the present study. Finally, the assumption of the fast isotropic wind generated when the mass-transfer rate exceeds the Eddington limit, which takes away the accretor angular momentum without interacting with the rest of the binary is only an approximation that would break down in cases of close binaries where the wind velocity is comparable to the orbital velocity and the outgoing wind will interact with the orbit, thereby taking away some orbital angular momentum (Mohamed & Podsiadlowski 2011; Schrøder et al. 2021).

As seen in Sect. 3, many of the various physical parameters that were explored in this study do not leave a distinct imprint on the combined synthetic XLF. However, these parameters may significantly affect some XRB subpopulations that are subdominant overall. With the advent of multiwavelength observational data sets, which will not only put constraints on the X-ray luminosity distributions of the overall XRB populations but also on the properties of their donors (Antoniou & Zezas 2016; Lazzarini et al. 2018, 2021; Antoniou et al. 2019) and accretors (Zezas et al. 2014; Wik et al. 2014; Maccarone et al. 2016; Yukita et al. 2016; Lazzarini et al. 2018; Vulic et al. 2018; Yao et al. 2021), we will be able to look into the XLFs of subpopulations and set more stringent constraints in population synthesis models and the associated stellar and binary evolution physics.

Future studies comparing the HMXB synthetic populations (their donor and accretor characteristics) to observations would help us better understand the physics involved. Detailed characterization studies have been done for HMXBs in various metallicity and SFR environments, showing differences in various cases. For instance, in the Large Magellanic Cloud (with an average metallicity of Z_solar/2), HMXBs are dominant in regions with recent star formation bursts (Antoniou & Zezas 2016), while in the SMC (with an average metallicity of Z_solar/4) HMXBs are found more in regions where star formation bursts occurred 25 to 60 Myr ago (Antoniou et al. 2010). We have mentioned various other uncertainties involved in binary evolution, such as mass-transfer efficiency or angular momentum losses, that have not been included in this study as their accurate behavior is uncertain, and varying them to investigate their effect is computationally expensive when detailed binary sequences are involved, as it would require recomputing entire the binary track library.

In our default population (refer to Fig. 1) there is an overabundance of XRBs compared to observations. In this section we discuss some of the effects of combinations of the different physical parameters that we have included in our study and possible causes for this discrepancy (for instance, transient behavior and assumptions of orbital eccentricity). We also identify four models that best describe the observed XLF.

4.1. Combinations of physical parameters

There is a certain level of degeneracy involved in the effects of the different parameters, on the XRB populations. So far we explored how the variation in a single parameter affects the XLF, but sometimes the combined effect of two or more parameters can have a more profound effect on the resulting synthetic XLF. Therefore, we discuss a few instances where the combined effects of two parameters have a significant change in the XLF and explore how some of the combinations reduce the so-called XLF bump. The combinations we discuss below do not correspond to the set of parameters that best reproduce the observed XLF, instead they showcase some examples of combined effects of the parameters on the synthetic XLF. In all the XLFs, the overabundance of XRBs below the luminosity of 10³⁸ erg s⁻¹ is caused by the population of Be XRBs and is highly dependent on the empirical relation used to calculate the luminosities, along with the assumptions regarding the identification of transient Be XRBs in the simulated populations.

For instance, from the three prescriptions of BH kick normalizations, fall-back kicks have the lowest number of RLO XRBs (for the default remnant mass prescription of Patton & Sukhbold 2020) at the point where the XLF moves from being dominated by wind-fed XRBs to RLO XRBs, shown by the change of slope at 10³⁹ erg s⁻¹ for the solid orange line in Fig. 7. On the other hand, fall-back kicks result in a larger number of wind-fed XRBs for luminosities of less than ∼10³⁸ erg s⁻¹ than the other two kick prescriptions. Following Hirai & Mandel (2021), the observability of some of these wind-fed XRBs may actually be limited (see Sect. 2.7). Figure 21 shows the result of these two effects combined, also comparing the resultant XLF to the default model (that has the parameter values shown in Table 2), and to the observations from Lehmer et al. (2019, pink shaded region). Wind-fed XRBs for X-ray luminosities below 10³⁸ erg s⁻¹ are reduced in number when the observable wind-fed accretion criterion is included. The combined effect of these two assumptions leads to a smoother XLF and a less pronounced bump around ∼10³⁸ erg s⁻¹ (up to a factor of ∼6 compared to observations).

Fig. 21. XLFs of populations showing the default population with the parameters in Table 2 (blue line) and the population with a combination of the fall-back SN kicks and the limited observability of wind BH XRBs following Hirai & Mandel (2021, orange line). The two curves are compared to the observed XLF from Lehmer et al. (2019).

Similarly, the combined effect of BH kicks with no normalization (see Sect. 2.4) and eccentric orbits that circularize at the onset of RLO at final orbital separations wider than the periastron distance (e.g., by assuming that the angular momentum of the orbit at the onset of RLO is conserved; see Sect. 3.3) also results in a drastic change in the synthetic XLF. Figure 22 shows the XLF for this case, comparing it to the default model and the observed XLF from (Lehmer et al. 2019, pink shaded region). The normalization of the synthetic XLF with the combined effects compares well to that of the observed XLF, particularly at 10³⁸ erg s⁻¹ where the slope changes for both the XLFs. However, the slope at luminosities greater than 10³⁹ erg s⁻¹ is steeper than the observed XLF.

Fig. 22. XLFs of populations showing the default population with the parameters in Table 2 (blue line) and the population when we do not normalize the BH kicks and conserve the angular momentum of the orbit at the onset of RLO (orange line). The two curves are compared to the observed XLF from Lehmer et al. (2019).

4.2. High-mass X-ray binary transients

Generally, known transient HMXBs are either Be XRBs or super-fast X-ray transients. Peak luminosities of Be XRBs go as high as ∼10³⁸ erg s⁻¹ (Raguzova & Popov 2005; Cheng et al. 2014), which is the low-luminosity end of the observed XLF from Lehmer et al. (2019). In all the XLFs we have studied thus far, we see a consistent population of Be XRBs, as shown in Figs. 1, 17, 15, and 20. Be XRBs are not affected by most of the prescriptions that we are investigating. For instance, BH kicks and wind-fed BH-XRB observability are merely affecting the Be-XRB population, 96% of which are wide NS XRBs in our simulations.

We should stress once again, however, that our treatment of Be-XRB X-ray luminosities is only a rough estimate based on empirical data of the peak luminosities of observed sources in the galaxy, which might overestimate their brightness (luminosities are calculated using the fitting formula from Dai et al. 2006), while still approximately matching the total number of Be XRBs to observations (as seen in Fig. 1). Thus, their contribution to the synthetic XLFs presented throughout the paper should only be considered as an order of magnitude estimate. Super-fast X-ray transients are transients with wind-fed OB-type supergiant donors (peak luminosities ∼10³⁶ erg s⁻¹; Negueruela et al. 2006). However, we do not consider super-fast X-ray transients in our calculations. The mechanisms that cause these systems to be transients are not fully understood and due to their low luminosities (which is at the sensitivity limit of most extragalactic surveys outside the local group), they provide limited insights by comparisons between our synthetic XLFs and the observed XLFs.

4.3. Instantaneous circularization at the onset of Roche-lobe overflow

During the formation of the CO, symmetric mass loss and asymmetric natal kicks often impart a significant eccentricity on the post-SN binary orbit. Tidal interactions between the binary components tend to circularize the orbit (Hut 1981; Zahn & Bouchet 1989; Verbunt & Phinney 1995), especially for the shorter period binaries. However, tides are often not efficient enough in doing so in the evolutionary timescales of the companion stars, which results in the donor stars overfilling their Roche lobes at the periastron of still eccentric orbits. Figure 23 shows the two-dimensional distributions of orbital period and eccentricity of binaries composed of a CO and non-degenerate companion star at the onset of RLO, for the RLO XRBs present in our default model. Also, 70% of the BH XRBs with H-rich donors undergoing RLO (that dominate the XLF bump) have eccentricities greater than 0.2 at the onset of RLO, with the distribution peaking around a relatively low eccentricity of 0.2 to 0.3. Even so, the XLF bump will be affected by how eccentric binaries are treated and since transient XRBs would be less observed compared to persistent XRBs, the appearance of the XLF bump should reduce. Therefore, for binaries with significant eccentricities at the onset of RLO, for example above 0.2, which constitute 68% of the total binaries reaching an RLO XRB phase in our default model, the subsequent treatment of the eccentricity evolution during the mass-transfer phase becomes important.

Fig. 23. Two-dimensional distributions of the eccentricity (e) and log10 of the orbital period (Porb) at the onset of RLO for the resulting XRBs, for the default model with parameters as described in Table 2.

A commonly made assumption in binary population synthesis codes is that the orbit circularizes instantaneously at the onset of RLO (e.g., Belczynski et al. 2008; Fragos et al. 2023). Then different assumptions can be made to estimate the post-circularization orbital separation. We discussed some of these assumptions in Sects. 2.5 and 3.3. We showed that they can have a significant impact on the properties of the resulting XRB population. The general trend we observed is that models assuming close circularized orbits at RLO onset (at the periastron) have a significant boost in the subpopulation of RLO XRBs with orbital periods of ≲10 days, which also tend to be long-lived. In the synthetic XLFs this contributes to the excess or bump at X-ray luminosities in the range 10³⁸ − 10³⁹ erg s⁻¹.

In practice, all of the prescriptions used in our study of the circularization process upon the onset of mass transfer are crude simplifying assumptions, and the variations we observe in the resulting XRB population highlight the need for better treatment of this process. A theoretical framework for the secular evolution of mass-transferring eccentric binaries has been developed in recent years. Additionally, Sepinsky et al. (2007, 2009); Sepinsky et al. (2010) showed that the circularization depends on the relative effect of the RLO phase and the tidal torques, and that mass transfer does not always circularize the orbit on a short enough timescale, in some cases even making the orbit more eccentric. Dosopoulou & Kalogera (2016a,b), Vick & Lai (2020) and Vick et al. (2021) have developed more complete theories for the secular evolution of such eccentric mass-transferring binaries. However, although in principle feasible, the computational complexity has not yet allowed the implementation of these theories within the context of detailed binary evolution models. Our results, in addition to some observed mass-transferring binaries in eccentric orbits (Boffin et al. 2014; Walter et al. 2015), highlight that the assumption of instantaneous circularization at the onset of RLO is a crude one, and further development of our modeling tools in that direction is imperative. A more physically realistic treatment of mass transfer in eccentric orbits will be included in future versions of POSYDON (Fragos et al. 2023).

4.4. Comparison to observations

The primary scientific goal of this work is to study signatures that different physical processes during the formation of XRBs may imprint in the XLF of a population. Throughout the paper, we have used the observational derived XLF of HMXBs by Lehmer et al. (2019) as a qualitative benchmark in order to infer if the aforementioned signatures would be identifiable. Here, we make a first attempt to make a quantitative comparison between the synthetic XLFs calculated from the 96 models run with the different combinations of parameters (see Appendix A) and the observationally derived one.

Using the Kolmogorov–Smirnov test, we measured the probability that samples drawn from our POSYDON populations (with ∼1.2 × 10⁴ XRBs in each population) are drawn from the same underlying population as the samples drawn from the observed XLF (10⁴ samples drawn). Even though this test is not sensitive at the tail ends of the distributions, we are more concerned with comparing the general shape and fit of the distribution. We use the one-sample Kolmogorov–Smirnov test to assess whether a synthetic XLF agrees with the observed XLF from Lehmer et al. (2019). We find disagreement at a level of 0.01 for all XLFs from the simulated populations (P_KS ≪ 0.01 in all cases), which suggests that our synthetic XLFs are not drawn from the same underlying population as the observed XLF.

To identify the “best-fitting” model, we used two approaches. Firstly, we match the slope of the XLFs and disregard the differences in the XLF normalization. We identify models that would define the observed population (requiring a correction for the normalization). We calculate the ratios of the number of XRBs in different bins of the XLF and compare them to ratios in the observed XLF. Secondly, we match the normalization of the XLFs with the observed ones and find the best-fitting models. However, following this approach, either the shape or slope of the selected XLF might not match the observed XLF. To compare the normalizations, we do a chi-squared test using the number of XRBs in each bin of the synthetic XLF and the number of observed XRBs in each bin.

We present the two best-fitting models from each of one of the two approaches in Fig. 24. The specific parameters for the models shown in the figure are described in Table 3. Models 55 and 64 were chosen based on the first approach (i.e., matching the overall shape of the observed XLF), but they overestimate the number of XRBs by up to one order of magnitude. Models 36 and 44 were chosen based on the second approach (i.e., matching the normalization of the observed XLF), and they reproduce the change of slope for the observations at luminosities of 10³⁸ erg s⁻¹. However, their overall shape is not consistent with observations. The two approaches selected two types of models with some key differences between the groups that cause the differences seen in Fig. 24. Models 55 and 64 have BH-mass normalized kicks that result in weaker kicks than in models 36 and 44, which have no kick normalization. The stronger kicks in models 36 and 44 disrupt a part of the close orbit binaries when the SN occurs, resulting in a lower number of XRBs at 10³⁸ erg s⁻¹, where the XLF starts to be dominated by RLO systems.

Fig. 24. XLFs of populations showing the four best-fit models, which were run with the parameters described in Table 3. The synthetic curves are compared to the observed XLF from Lehmer et al. (2019).

5. Conclusions

In this work we have investigated the effects of different physical assumptions in the context of XRB formation in young stellar populations. We ran populations through 96 models with different combinations of parameters (described in Table 1). For this, we used POSYDON, which is a new binary population synthesis code that utilizes detailed stellar structure and binary evolution models throughout the entire evolution of a binary. The populations were run assuming a constant SFR over a time duration of 100 Myr. We created synthetic XLFs to study the simulated XRBs and compared them to observed HMXB XLFs from Lehmer et al. (2019) to obtain insights into the physics involved. The individual effects of the different physics parameters were discussed, and four best-suited models that most closely match the observations were chosen. We summarize our main findings in the following list:

• The HMXB XLFs we generated have a more complex shape than a single power law (for instance, see Fig. 1). This qualitative effect is also seen in the empirical XLF from observations. In our simulations, the shape changes at a luminosity of 10³⁸ erg s⁻¹ (similar to observations) due to the population of RLO BH XRBs with H-rich donors.

• There is an overabundance of XRBs compared to observations (up to a factor of 10) for certain model parameter combinations at intermediate X-ray luminosities (referred to as the XLF bump, seen around luminosities of 10³⁸ to 10³⁹ erg s⁻¹ in Fig. 1), primarily due to RLO BH XRBs with H-rich donors.

• Increasing the strength of SN kicks in BH XRBs does not necessarily disrupt close binaries. Instead, it results in a larger number of highly eccentric binaries that end up as bright RLO XRBs, compared to weaker kicks. However, strong kicks disrupt wide binaries, reducing the number of bright wind-fed BH XRBs (by a factor of ∼20).

• Changing the core-envelope boundary does not have a significant effect on the surviving populations of NS and BH binaries. However, a more inefficient CE phase, due to either lower assumed α_CE values or higher envelope binding energies, would result in closer orbits leading to more XRBs.

• Including a limiting criterion for the observability of wind BH HMXBs, for example by requiring the formation of an accretion disk around the BH as in Hirai & Mandel (2021), does not affect XRB luminosities greater than approximately 10³⁸ erg s⁻¹, and decreases the slope of the XLF at lower luminosities. However, this criterion does not affect most Be XRBs (∼96% of which are NS XRBs) as it applies to only wide BH XRBs.

• Commonly used assumptions such as instantaneous circularization at the periastron with the onset of RLO might not be accurate for binaries that have high eccentricities when the donor fills its Roche lobe (∼68% have eccentricities above 0.2). These eccentric sources might appear as transient HMXBs, which should be taken into account when simulating XLFs.

• No combination of the explored physical parameters (combinations of parameters described in Table 1 and the XLFs shown in Fig. A.1 and described in Table A.1) matches well both the shape and the normalization of the observed XLF simultaneously. However, in this work we did not explore any variations of the distributions of initial binary properties, such as variations in the initial mass function, orbital separation distributions, binary mass ratios, and the binary fraction, which are expected to primarily change the overall normalization of the modeled XLFs, rather than their shape.

Even though our synthetic XLFs do not always agree with observations, this work reveals the importance of large-scale parameter studies. Future comparative work that includes multiwavelength characterization studies of XRB populations across various metallicities and SFHs would help constrain the physics involved in binary evolution.

Acknowledgments

The authors thank the anonymous referee for their constructive comments that helped improve the manuscript. This work was supported by the Swiss National Science Foundation Professorship Grant (PP00P2_176868; PI: Fragos). This work was also supported by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie RISE action, grant agreement No. 873089 (ASTROSTAT-II). D.M. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101002352). K.K. acknowledges support from the Federal Commission for Scholarships for Foreign Students for the Swiss Government Excellence Scholarship (ESKAS No. 2021.0277). E.Z. acknowledges funding support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 772086).

References

Appendix A: XLFs of the simulated population

We show all the cumulative XLFs from the 96 simulations carried out for this study in Fig. A.1, and all the corresponding differential forms of the XLFs in Fig. A.2. Further information about the combination of parameters that were used to run these populations are shared in Table A.1.

Fig. A.1. Synthetic XLFs of populations 1 to 24 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

Fig. A.2. Synthetic XLFs of populations 25 to 48 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

Fig. A.3. Synthetic XLFs of populations 49 to 72 from Tables A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

Fig. A.4. Synthetic XLFs of populations 73 to 96 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

Fig. A.5. Synthetic differential forms of the XLFs for populations 1 to 24 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

Fig. A.6. Synthetic differential forms of the XLFs for populations 25 to 48 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

Fig. A.7. Synthetic differential forms of the XLFs for populations 49 to 72 from Tables A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

Fig. A.8. Synthetic differential forms of the XLFs for populations 73 to 96 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.

All Tables

All Figures

	Fig. 1. XLF of the default population (solid blue line) compared to the observed XLF from Lehmer et al. (2019, pink shaded region). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.
In the text

	Fig. 2. Differential form of the synthetic XLFs, showing the default population (solid blue line) compared to the observed XLF from Lehmer et al. (2019, pink shaded region). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.
In the text

Fig. 3. Distribution of the calculated X-ray luminosities of the simulated XRBs in the default population with respect to the orbital separation. The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with gray circles. Large and small symbols denote massive (with ≳8.0 M⊙) and less massive (with < 8.0 M⊙) XRBs, respectively.

In the text

	Fig. 4. XLFs of the default population with the parameters in Table 2 showing the Patton & Sukhbold (2020) SN mechanism (blue line) and the population with the Fryer et al. (2012, delayed) SN mechanism (orange line).
In the text

Fig. 5. Distribution of newly formed CO masses (MCO) for binaries that remain bound post-SN, for the default population with the parameters in Table 2, including the Patton & Sukhbold (2020) SN mechanism (blue line) and the population with the Fryer et al. (2012, delayed) SN mechanism (orange line). The populations are further split by the type of CO, NS (solid line), and BH (dashed line). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the default population for NS formation (NS: Patton & Sukhbold 2020, solid blue line) to 1.

In the text

	Fig. 6. XLFs of the default population with the parameters in Table 2 showing the population with BH mass normalized SN kicks (blue line), the population with fall-back normalized SN kicks (orange line), and the population with no normalization on the SN kicks (green line).
In the text

Fig. 7. XLFs of the default population with the parameters in Table 2 showing the population with BH mass normalized SN kicks (blue lines), the population with fall-back normalized SN kicks (orange lines), and the population with no normalization on the SN kicks (green lines). The populations follow the remnant mass prescription from Patton & Sukhbold (2020). The populations are further split by the type of mass transfer that occurred, RLO (solid line) or wind-fed accretion (dashed line). For the sake of clarity, the Be-XRB population is not shown as it is not affected by the kick normalization and obscures the effects of the BH kick normalization.

In the text

Fig. 8. XLFs of the default population with the parameters in Table 2 showing the population with BH mass normalized SN kicks (blue lines), the population with fall-back normalized SN kicks (orange lines), and the population with no normalization on the SN kicks (green lines). The populations follow the remnant mass prescription from Fryer et al. (2012, delayed). The populations are further split by the type of mass transfer that occurred, RLO (solid line) or wind-fed accretion (dashed line). For the sake of clarity, we did not include Be XRBs in this figure as they are not affected by the kick normalization and obscure the effects of the BH kick normalization.

In the text

Fig. 9. Distributions of orbital periods (Porb) at the end of the SN for the progenitors of BH XRBs, for the default population with the parameters in Table 2, including the population with BH mass normalized SN kicks (solid blue), the population with fall-back normalized SN kicks (hatched orange), and the population with no normalization on the SN kicks (hatched green). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the default population (BH mass normalized kicks shown in solid blue) to 1.

In the text

	Fig. 10. Distribution of the eccentricity (e) versus log10 of the orbital period (Porb) for the same populations as in Fig. 9.
In the text

	Fig. 11. XLFs of the default population with the parameters in Table 2 showing circularization at RLO onset at the periastron of the eccentric orbit (blue line), the population with the circularization at an orbit conserved angular momentum (orange line), and the population with circularization at the binary separation of the eccentric orbit (green line).
In the text

Fig. 12. XLFs of the default population with the parameters in Table 2 showing the population with pre-RLO circularization at the periastron (blue lines), at an orbit with conserved angular momentum (orange lines), and at the binary separation of the pre-RLO eccentric orbit (green lines). The populations are further split by the type of mass transfer that occurred, RLO (solid line), and wind-fed accretion (dashed line). Be XRBs are shown with a dotted gray line, showing no effect as they are not circularized. AM stands for angular momentum.

In the text

Fig. 13. Distributions of orbital periods (Porb) of circularized orbits at the onset of RLO for the resulting XRBs. The presented systems were identified as initiating RLO at the periastron, and then the separation of the circularized orbit was adjusted following the respective model used. The figure shows the default population with the parameters in Table 2 including circularization at RLO onset at the periastron of the eccentric orbit (solid blue), the population with the circularization at an orbit conserved angular momentum (hatched orange), and the population with circularization at the binary separation of the eccentric orbit (hatched green). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the default population (with orbits circularized at the periastron shown in solid blue) to 1.

In the text

	Fig. 14. XLF of the default population with the parameters in Table 2 showing the population with CE efficiency (αCE) of 1.0 (blue line) and the population with αCE of 0.3 (orange line).
In the text

	Fig. 15. XLF of the population with αCE = 0.3 (solid orange line). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.
In the text

	Fig. 16. XLFs of the default population with the parameters in Table 2 showing the population with core-envelope definition at XH = 0.3 (blue line) and the population with core-envelope definition at XH = 0.01 (orange line).
In the text

	Fig. 17. XLF of the population with core-envelope definition at XH = 0.01 (solid orange line). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS), and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.
In the text

Fig. 18. Distribution of the post-CE orbital periods (Porb) for the default population with the parameters in Table 2 including the core-envelope definition at XH = 0.30 (blue histogram) and the population with the core-envelope definition at XH = 0.01 (orange histogram). The populations are further split by the type of CO in the surviving binary, showing NSs (step) and BHs (hatched). The distributions have been weighted with the same factor that normalizes the maximum counts per bin for the NS-He post-CE binary population to 1.

In the text

	Fig. 19. XLFs of the default population with the parameters in Table 2 showing no criterion for wind-fed accretion (blue line) and the population with limited observability of wind-fed XRBs following Hirai & Mandel (2021, orange line).
In the text

	Fig. 20. XLF of population with limited observability of wind-fed XRBs following Hirai & Mandel (2021). The population is further split by the type of XRB (RLO or wind), type of accretor (BH or NS) and type of donor (H-rich or He-rich). Be XRBs are shown with a dotted gray line.
In the text

	Fig. 21. XLFs of populations showing the default population with the parameters in Table 2 (blue line) and the population with a combination of the fall-back SN kicks and the limited observability of wind BH XRBs following Hirai & Mandel (2021, orange line). The two curves are compared to the observed XLF from Lehmer et al. (2019).
In the text

	Fig. 22. XLFs of populations showing the default population with the parameters in Table 2 (blue line) and the population when we do not normalize the BH kicks and conserve the angular momentum of the orbit at the onset of RLO (orange line). The two curves are compared to the observed XLF from Lehmer et al. (2019).
In the text

	Fig. 23. Two-dimensional distributions of the eccentricity (e) and log10 of the orbital period (Porb) at the onset of RLO for the resulting XRBs, for the default model with parameters as described in Table 2.
In the text

	Fig. 24. XLFs of populations showing the four best-fit models, which were run with the parameters described in Table 3. The synthetic curves are compared to the observed XLF from Lehmer et al. (2019).
In the text

	Fig. A.1. Synthetic XLFs of populations 1 to 24 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text

	Fig. A.2. Synthetic XLFs of populations 25 to 48 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text

	Fig. A.3. Synthetic XLFs of populations 49 to 72 from Tables A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text

	Fig. A.4. Synthetic XLFs of populations 73 to 96 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text

	Fig. A.5. Synthetic differential forms of the XLFs for populations 1 to 24 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text

	Fig. A.6. Synthetic differential forms of the XLFs for populations 25 to 48 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text

	Fig. A.7. Synthetic differential forms of the XLFs for populations 49 to 72 from Tables A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text

	Fig. A.8. Synthetic differential forms of the XLFs for populations 73 to 96 from Table A.1, showing the types of mass transfer that occurred, the donors, and the accretors.
In the text