© The Authors 2025

1. Introduction

Ultra-luminous X-ray sources (ULXs) are extragalactic point-like sources that exhibit extraordinarily high luminosities in the X-ray (L_X > 10³⁹ erg/s), exceeding the Eddington luminosity for black holes (BHs) of mass ≈10 M_⊙. The first ULXs were discovered using the Einstein Observatory, designed to obtain images of distant galaxies (for a review, see Fabbiano 1989). Early X-ray surveys revealed a considerable quantity of these sources. Fabbiano (1989) identified 16 sources with X-ray luminosities L_X > 10³⁹ erg/s. To date, over 1800 ULXs have been identified, increasing our knowledge and understanding of ULX sources (Kaaret et al. 2017; Walton et al. 2022; King et al. 2023).

Since their discovery, the debate over what powers these systems has been a continuous topic of both theory and observation. Theoretically, there are three possibilities as to how these systems reach such high luminosities: (1) accretion by intermediate mass black holes (IMBH), 10² to 10⁵ M_⊙, (Colbert & Mushotzky 1999): this increases the mass of the BH, and the system accretes at sub-Eddington rate; (2) beamed emission from a stellar-mass BH or neutron star (NS) (Georganopoulos et al. 2002): this avoids super-Eddington accretion and keeps the mass of the compact object below 100 M_⊙; and (3) super-Eddington accretion (e.g. Begelman 2001; Finke & Böttcher 2007; Sądowski & Narayan 2015): accretion onto stellar mass BHs or NSs emitting with a super-Eddington luminosity. We note that a system can simultaneously be beamed and supercritical. Observations have shown evidence for IMBHs (Lasota et al. 2011; Mezcua et al. 2013), stellar mass BHs (Middleton et al. 2013), and pulsations detected in a few systems indicate the presence of an NS (Bachetti et al. 2014; Fürst et al. 2016) as the primary object in these systems. It has been suggested that ULXs could be the progenitors of compact object mergers detected by LIGO (Inoue et al. 2016; Finke & Razzaque 2017; Mondal et al. 2020).

The debate on ULXs centres around the mass of the compact objects and their accretion rate. While X-ray data have been instrumental in our understanding of ULXs, they have limitations. They often do not provide information about the companion of the compact object. Identifying the companion can offer valuable insights into the masses of these binary systems. An example of the significance of this approach can be found in the confirmation of the existence of Cygnus X-1, the first known Galactic BH (Murdin & Webster 1971). ULXs also emit radiation in other wavelengths, including the ultraviolet (UV) and optical bands (Gladstone et al. 2013).

Previous observational studies have provided valuable insights into the connection between X-ray and UV radiation in ULXs, suggesting a common origin or interconnected emission mechanisms. Sonbas et al. (2019) investigated the relationship between the X-ray and UV wavebands in several ULXs with known optical counterparts by comparing the ratio of flux densities and by utilising the spectral optical–X-ray index, α_ox. Their analysis reveals a significant anti-correlation between α_ox and the UV luminosity at 2500 Å (L_2500 Å). However, due to the complexity of ULX systems and the limitations of observational data, it is challenging to fully explore and understand the correlation between X-ray and UV radiation based only on observations. Therefore, simulated studies that follow the evolution of binary systems and model the emission processes to provide detailed insights into the spectral behaviour are essential. Simulated studies enable us to disentangle the complex interplay of physical mechanisms in astrophysical sources, allowing us to explore a wide parameter space and investigate various emission scenarios. By employing a population synthesis code, we simulated a range of ULX scenarios and systematically studied the spectral correlation between X-ray and UV radiation.

In this work, we used the population synthesis code COSMIC (Breivik et al. 2020) to simulate the evolution of binary systems and study the relationship between the UV and X-ray emission of ULXs. We explored the effects of metallicity on the synthetic population of ULXs and the α_ox relation for BH-ULXs and NS-ULXs and compared them with observational data. The layout of the paper is as follows: in Sect. 2, we discuss our COSMIC simulation set-up; the UV- and X-ray-emission models are discussed in Sect. 3; and in Sect. 4 we present and discuss our results. We summarise and conclude our study in Sect. 5.

2. Synthetic binary population

To study the correlation between X-ray and UV radiation in ULXs, we employed the COSMIC (Breivik et al. 2020) binary population synthesis code. COSMIC combines theoretical models of stellar evolution, binary interactions, and compact object mergers to simulate the formation and evolution of binary systems. We modeled the evolution of binaries starting from the zero-age main sequence (ZAMS) until the formation of the X-ray binary systems.

In our simulations, we drew the primary (most massive) ZAMS mass, M₁, from a Kroupa initial mass function (IMF; Kroupa et al. 1993). This IMF is given by dN/dM₁ ∝ M₁^−α, with α = 1.3 for 0.08 M_⊙ < M₁ ≤ 0.5 M_⊙; α = 2.2 for 0.5 M_⊙ < M₁ ≤ 1 M_⊙; and α = 2.7 for M₁ > 1.0 M_⊙. The ZAMS mass of the secondary star, M₂ (ranging from 0.5 M_⊙ to 150 M_⊙), was drawn from a flat distribution of the binary mass ratio q = M₂/M₁ within the range of [0.1, 1.0]. The orbital period (P) was drawn from a distribution of f(log₁₀P/d)≈(log₁₀P/d)^−0.55, with log₁₀P/d ranging from 0.15 to 5.5, and the eccentricity (e) was drawn from a distribution that followed f(e)≈e^−0.42 within the interval of [0.0, 0.9] (Sana et al. 2013).

In our study, we simulated 2 × 10⁶ binary systems with varying metallicities, ranging from 0.5% of the solar metallicity (Z_⊙) to 1.5 Z_⊙. The precise value of Z_⊙ remains uncertain (Vagnozzi et al. 2017). We adopted Z_⊙ = 0.02 as the benchmark value for consistency with previous studies (e.g. Mondal et al. 2020). Additionally, we set the binary fraction to be 100% for a primary mass greater than 10 M_⊙ (Sana et al. 2013). The configurations in our standard COSMIC model parameters are as follows.

1. Common envelope (CE) phase: CE efficiency of α_CE = 1, we follow the pessimistic scenario by letting stars without a core-envelope boundary to automatically lead to merger during the CE phase (Belczynski et al. 2007).

2. Supernova model: to estimate the mass of the final compact object after the supernova explosion, we followed the rapid supernova model (Fryer et al. 2012).

3. Accretion: studies using simulations have demonstrated that mass accretion onto a compact object can exceed the Eddington limit (McKinney et al. 2014). It is now considered possible for super-Eddington accretion to occur in X-ray binaries with both BHs and NSs as compact objects (Van Den Eijnden et al. 2019; Vasilopoulos et al. 2020). Therefore, we allowed super-Eddington accretion to reach values of up to 1000 times the Eddington limit.

4. Critical mass ratios (q_c): q_c, is the mass ratio between the donor and accretor at the beginning of Roche-lobe overflow (RLOF). It is used to determine whether a system will experience a CE phase or continue with stable mass transfer (MT) based on the mass ratio of the components at the beginning of RLOF. Mass transfer through RLOF is only allowed for a thermal timescale if the mass ratio at the beginning of RLOF is below the q_c. If the mass ratio ≥q_c, the mass transfer occurs on a dynamic timescale, resulting in the initiation of the CE phase. To determine q_c, we followed the treatment in Sect. 5.1 of Belczynski et al. (2008).

We sampled a population of main-sequence binary stars that evolved from the ZAMS until the formation of X-ray binary systems. Our ULX population is selected based on the value of X-ray luminosity, L_X > 10³⁹ erg/s. In the next section, we discuss how our X-ray and UV luminosities are calculated.

3. Emission from the accretion disc and donor star

3.1. Total X-ray emission

The X-ray emission from ULXs primarily originates from the accretion processes onto the compact objects: either stellar-mass BHs or NSs. In these systems, a massive donor star transfers mass to the compact object through RLOF or wind accretion. The high luminosities observed in ULXs, exceeding the Eddington limit for typical stellar-mass BHs, suggest the presence of highly efficient accretion processes. Accretion models play a crucial role in our understanding of the X-ray emission from ULXs and require the careful consideration of their physical parameters and radiative processes.

In a compact binary system, mass is transferred from one star to another. This process leads to the creation of an accretion disc. We calculate an isotropic bolometric X-ray luminosity (L_X, iso) from accretion as (Shakura & Sunyaev 1973; Poutanen et al. 2007)

$$\begin{array}{c}\hfill L_{\mathrm{X},\mathrm{iso}}=\{\begin{array}{cc}{L}_{\mathrm{Edd}}[1+ln\dot{m}]\hfill & \mathrm{if}\dot{m}>1\hfill \\ \eta \dot{M}{c}^2\hfill & \mathrm{if}\dot{m}\le 1.\hfill \end{array}\end{array}$$(1)

Here, Ṁ is the mass-accretion rate, c is the speed of light, ṁ is the ratio Ṁ/Ṁ_Edd between the mass accretion rate and the Eddington accretion rate Ṁ_Edd = 12L_Edd/c², L_Edd is the Eddington luminosity of a star or compact object with mass M given by

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=1.26\times {10}^{38}\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\mathrm{erg}\mathrm{s}}^{-1},\end{array}$$(2)

η is the accretion efficiency defined as

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

for NSs, and

$$\begin{array}{c}\hfill {\eta }_{\mathrm{BH}}=1-E(R_{\mathrm{ISO}})\end{array}$$(4)

for BHs, where G is Newton’s gravitational constant. The function E(R) is given by

$$\begin{array}{c}\hfill E(R)=\frac{R^2-2\frac{\mathit{GM}}{c^2}R+a\frac{\mathit{GM}}{c^2}{\left(\frac{\mathit{GM}}{c^2}R\right)}^{1/2}}{R{[R^2-3\frac{\mathit{GM}}{c^2}R+2a\frac{\mathit{GM}}{c^2}{\left(\frac{\mathit{GM}}{c^2}R\right)}^{1/2}]}^{1/2}},\end{array}$$(5)

where a is the BH spin and R_ISO is the radius of the innermost stable orbit,

$$\begin{array}{cc}& R_{\mathrm{ISO}}=\frac{R_{\mathrm{Sch}}}{2}(3+Z_2-{\left[(3-Z_1)(3+Z_1+2Z_2)\right]}^{1/2}),\hfill \end{array}$$(6)

$$\begin{array}{cc}& Z_1=1+{(1-\frac{4a^2}{R_{\mathrm{Sch}}^2})}^{1/3}[{(1+\frac{2a}{R_{\mathrm{Sch}}})}^{1/3}+{(1-\frac{2a}{R_{\mathrm{Sch}}})}^{1/3}],\hfill \end{array}$$(7)

$$\begin{array}{cc}& Z_2={(\frac{12a^2}{R_{\mathrm{Sch}}^2}+Z_1^2)}^{1/2},\hfill \end{array}$$(8)

where R_Sch = 2GM/c² is the Schwarzschild radius (Brown et al. 2000).

When the accretion rate is high, the luminosity can be focused into narrow cones, resulting in greater observed luminosity compared to the isotropic emission. This phenomenon is referred to as ‘beaming’. In order to explore the potential impact of the beaming factor, we followed King (2009) to calculate the beaming factor (b), which is given by

$$\begin{array}{c}\hfill b=\{\begin{array}{cc}\frac{73}{{\dot{m}}^2}\hfill & \mathrm{if}\dot{m}>8.7\hfill \\ 1\hfill & \mathrm{if}\dot{m}<8.7.\hfill \end{array}\end{array}$$(9)

The beamed X-ray luminosity is then given by

$$\begin{array}{c}\hfill L_{\mathrm{X},\mathrm{b}}=\frac{L_{\mathrm{X},\mathrm{iso}}}{b},\end{array}$$(10)

where L_X, iso is the X-ray luminosity in Eq. (1).

3.2. ULX spectrum

In addition to their intense X-ray emission, ULXs emit in the optical and UV. The presence of UV emission in conjunction with X-ray radiation opens a window to understanding the broader radiative processes at play in these enigmatic systems. This dual emission may hold clues to the intricate interconnections between the various components of the accretion process, the accretion disc, companion stars, and the surrounding environment. To study ULXs in both the UV and X-ray bands, it is essential to consider contributions from both the accretion disc and the stellar companion. The X-ray emission primarily originates from the accretion disc, whereas the UV emission has two components: one from the disc and another from the stellar companion. In this section, we describe the models utilised in our calculations.

3.2.1. Accretion disc

Ultra-luminous-X-ray-source spectral modelling has been a focus in several studies using models such as DISKBB (Mitsuda et al. 1984) and DISKIR (Gierliński et al. 2008, 2009). These models often imply the presence of high-mass black holes in the ULX systems, which are not the focus of this study. Vinokurov et al. (2013) proposed a physically motivated model by examining the spectral energy distribution (SED) of supercritical accretion discs (SCADs) and applied this model to fit the observed SEDs of five ULXs.

According to the SCAD model, the accretion discs become thick and develop a strong wind, forming a wind funnel above the disc at distances shorter than the spherisation radius, with a temperature dependence given by T ∝ r^−1/2. The disc’s luminosity in this region is approximated by L_bol ≈ L_Edd ln (ṁ). Outside the spherisation radius, the disc remains thin, and its total luminosity equals the Eddington luminosity, L_Edd. The thin disc heats the wind from below, and from the inner side of the funnel the wind is heated by the supercritical disc (see Fig. 2 of Vinokurov et al. 2013).

When applied to ULX spectra, the SCAD model suggests stellar-mass BHs as the compact object, whereas the DISKIR model often yields much higher BH masses, suggesting IMBHs. The SCAD model accounts for the high-energy curvature observed in ULX spectra and the presence of outflows without invoking IMBHs. It explains the observed luminosities and spectral features using stellar-mass BHs undergoing supercritical accretion, which is consistent with the lack of a large population of IMBHs observationally. Since all our simulated ULXs have either stellar-mass BHs or NSs as the compact objects, we used the SCAD model to describe the supercritical accretion discs in our simulated ULX population. The temperature of the disc in this model is given by

$$\begin{array}{c}\hfill T(r)=T_{\mathrm{in}}\{\begin{array}{cc}{(rsin{\theta }_f)}^{-1/2},\hfill & 1\le r\le r_{\mathrm{sp}}\hfill \\ [\frac{f_{\mathrm{out}}}{sin{\theta }_f}(1+ln\dot{m}){]}^{1/4}{r}^{-1/2},\hfill & r_{\mathrm{sp}}\le r\le r_{\mathrm{ph}}\hfill \end{array},\end{array}$$(11)

where T_in⁴ = L_Eddsin(θ_f)/(4πσR_in²), R_in = 3R_Sch = 6GM/c², M is the compact object mass, θ_f is the funnel angle, σ is the Stefan-Boltzmann constant, f_out is the fraction of bolometric flux thermalised in the disc, and r = R/R_in is the radius of the accretion disc. Here, r_sp = Ṁ/Ṁ_Edd is the spherisation radius, and $r_{\mathrm{ph}}=\dot{M}\kappa /(4\pi v_{w}cos{\theta }_f{R}_{\mathrm{in}})$ is the photospheric radius, where κ is the opacity (κ ≈ κ_es = 0.34 cm²/g for fully ionised gas) and v_w = 1000 km s⁻¹ is the wind velocity. With our COSMIC data, this leaves only f_out and θ_f as free parameters in this model. The effect of these two free parameters is discussed in the results section. The luminosity, at a given frequency ν and radius r within a thickness of dr of the disc, can be expressed as

$$\begin{array}{c}\hfill \mathrm{d}{L}_{\nu }(r)=2\pi R_{\mathrm{in}}^2·r{B}_{\nu }(T(r))·\mathrm{d}r,\end{array}$$(12)

where

$$\begin{array}{c}\hfill B_{\nu }(T)=\frac{2h{\nu }^3}{c^2}·\frac{1}{e^{\frac{h\nu }{\mathit{kT}}}-1}\end{array}$$(13)

is the Planck function describing black-body radiation at a given temperature, T, and a given frequency, ν. Here, h is the Planck constant and k is the Boltzmann constant. Therefore, the total luminosity as a function of frequency, considering the temperature gradient, is

$$\begin{array}{c}\hfill L_{\nu }=2\pi R_{\mathrm{in}}^2{\displaystyle {\int }_1^{r_{\mathrm{ph}}}}{B}_{\nu }(T(r))r\mathrm{d}r.\end{array}$$(14)

3.2.2. Stellar companion

Another component of the ULX spectrum is the companion star, which mainly contributes to the UV band. Observations have indicated that the most likely donor stars are massive and hot OB-type stars, which can contribute significantly to the UV emission (Liu et al. 2004; Roberts et al. 2008; Gladstone et al. 2013). Other studies have found evolved red super-giant companions (Heida et al. 2015, 2019). This emission can also be attributed to the high-energy radiation absorbed and re-emitted by the donor star’s surface (Middleton et al. 2015). The UV luminosity from a companion star of radius (R_⋆) and effective temperature (T_⋆) is obtained by

$$\begin{array}{c}\hfill L_{\mathrm{UV},\star }=\pi R_{\star }^2\nu B_{\nu }(T_{\star }),\end{array}$$(15)

where we use λ = c/ν = 2500 Å. The UV contribution from the disc (L_UV, disk) is obtained from Eq. (14) at λ = 2500 Å. Then, the total UV emission from the ULX is given by

$$\begin{array}{c}\hfill L_{\mathrm{UV}}=L_{\mathrm{UV},\star }+L_{\mathrm{UV},\mathrm{disk}}.\end{array}$$(16)

4. Results and discussion

Our study aims to simulate a population of ULXs, investigate their properties, and compare them to the observed population. In this section we show our results and compare them to the observational data where available.

4.1. Number of ULXs

Various studies have explored the relationship between the number of ULXs, star formation rate (SFR), and metallicity. Using the population synthesis code StarTrack, Wiktorowicz et al. (2017) simulated a robust correlation between the number of ULXs and the SFR history. Meanwhile, Mondal et al. (2020), also using StarTrack, found an increase in the number of ULXs with decreasing metallicity in their simulations, specifically in the context of stars born at the same time (i.e. a single burst). In an X-ray survey of a sample of 66 galaxies, Mapelli et al. (2011) found that there is no significant direct correlation between the number of ULXs (N_ULX) and metallicity. However, when normalising the number of ULXs by the SFR, (N_ULX/SFR), a marginally significant anti-correlation with metallicity is observed. This indicates that lower metallicity environments produce more ULXs per unit of star formation compared to higher metallicity environments. These studies show that metallicity has some influence on the efficiency of ULX formation, but this effect is secondary to the influence of SFR. Mapelli et al. (2011) also found the power-law index of $\alpha =0.{55}_{-0.19}^{+0.21}$ for N_ULX/SFR, where the number of ULX is normalised by the SFR and $\alpha =0.{16}_{-0.28}^{+0.28}$ for N_ULX as a function of metallicity.

In Fig. 1, we show N_ULX as a function of metallicity for our simulations. We note that we define N_ULX from our synthetic population at the end of the ULX phase. The time evolution of ULX population properties will be the focus of future studies. The total number of ULXs, as well as the sub-populations of BH-ULXs and NS-ULXs, show a declining trend with increasing metallicity. The slope for BH-ULXs is steeper (α = 0.21 ± 0.02) compared to the total ULX population (α = 0.16 ± 0.01) and NS-ULXs (α = 0.09 ± 0.01). At higher metallicities, the populations of NS-ULXs and BH-ULXs become comparable. The steeper decline for BH-ULXs indicates that the formation of BH-ULXs is more sensitive to metallicity compared to NS-ULXs. The shallower slope for NS-ULXs suggests that while metallicity does affect their formation, it is not as critical as for BH-ULXs. This aligns with the theoretical understanding that massive metal-poor stars are more likely to form massive BHs through direct collapse (Heger et al. 2003). Our power-law index of α = 0.16 ± 0.01 for the combined population of BH- and NS-ULXs is consistent with the slope from Mapelli et al. (2011).

Fig. 1. Number of ULXs formed in our simulation (NULX) as a function of metallicity. The y-axis indicates the number of ULXs derived from the simulations, with emphasis placed on the trend’s shape rather than the absolute values. The solid lines represent fits to the simulated data, with the corresponding slopes mentioned in the legend.

The parameters f_out and θ_f in Eq. (11) play a crucial role in shaping the UV and X-ray emission characteristics of ULXs. A higher fraction of the bolometric flux thermalised in the disc (i.e. larger f_out) results in stronger UV emission. Conversely, for any observer viewing the system down the funnel formed by the disc combination, smaller θ_f angles will result in enhanced X-ray emission by directing more of the X-ray emission into the increasingly narrow funnel. For our simulated populations, we set f_out to the average value derived from the systems analysed in Vinokurov et al. (2013), adopting f_out = 0.03. A detailed discussion of the impact of varying f_out on our results can be found in Sect. 4.3. This value of f_out ensures that a small but significant fraction of the bolometric flux is thermalised in the accretion disc, contributing to UV emission. Given the high accretion rates typical of ULXs, this thermalisation fraction effectively reproduces the observed emission without overestimating the flux.

We studied our population with θ_f set to 45° and treat this as a non-beamed population (see Sect. 4.2). It should be noted that θ_f = 45° was selected as a standard, intermediate value within the allowed range of 0° −90° for exploratory purposes and does not imply any inherent physical property of the system. In Sect. 4.3.1, we calculate θ_f from the accretion rate using the beaming model by King (2009) and compare the two scenarios.

4.2. Results with constant θ_f

4.2.1. X-ray luminosity function

One approach to studying the ULX population involves studying their X-ray luminosity function. Figure 2 shows the cumulative XLF of the simulated ULX population at 40% of Z_⊙ in the 0.2−12 keV band. We fit this population with a function,

$$\begin{array}{c}\hfill N(>L_{\mathrm{X}})=C{L}_{\mathrm{X}}^{-\alpha }exp[-{(L_{\mathrm{X}}/L_s)}^{\beta }],\end{array}$$(17)

Fig. 2. Cumulative XLF of the simulated ULX population at 40% Z⊙ in the 0.2−12 keV band.

and obtained power-law index α = 1.99 ± 0.02, β = 3.92 ± 3.70, and L_s = (43.7 ± 9.94)×10³⁹ erg s⁻¹ at this metallicity. However, the power-law index changes with metallicity. The XLF for a high-metallicity population is steeper than a low-metallicity population; it ranges from α = 1.77 ± 0.01 at 0.05 Z_⊙ to 2.15 ± 0.38 at 1.5 Z_⊙, showing a relative excess of high-luminosity ULXs at lower metallicities (see Table 1). Similarly, Lehmer et al. (2021) observed an increasing number of luminous sources with declining metallicity, indicating that low-metallicity environments have more high-luminosity high-mass X-ray binary (HMXB) sources.

4.2.2. Optical–X-ray spectral index (α_ox)

The X-ray-to-UV ratio is described as the interband spectral index (Tananbaum et al. 1979),

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{ox}}=\frac{{log}_{10}(L_{\nu ,\mathrm{X}}/L_{\nu ,\mathrm{UV}})}{{log}_{10}({\nu }_{\mathrm{X}}/{\nu }_{\mathrm{UV}})},\end{array}$$(18)

where both L_ν, UV and L_ν, X are expressed as monochromatic luminosities (in units of erg/s/Hz) defined at ν_UV = 1.2 × 10¹⁵ Hz (2500 Å) and ν_X = 4.8 × 10¹⁷ Hz (2 keV), respectively, following the work of Sonbas et al. (2019), who found an anti-correlation between α_ox and L_ν,UV in observational data. In the context of the SCAD model, discussed in Sect. 3.2.1, such a relationship is expected, where the UV luminosity is dominated by emission from the funnel and depends on the accretion rate (ṁ), while the X-ray luminosity is dominated by emission from the inner accretion disc and does not depend on ṁ (see Fig. 3 in Vinokurov et al. 2013). Both the components depend on the mass of the primary in the same way. Systems with greater ṁ have greater UV luminosity but the same X-ray luminosity (for a given primary mass M). So, as ṁ increases, L_ν,UV also increases, and the ratio $L_{\nu ,X}/L_{\nu ,\mathrm{UV}}\propto {\alpha }_{\mathrm{ox}}$ decreases. This leads to the anti-correlation between α_ox and L_ν,UV. The diversity of this relation, however, depends on various factors such as the metallicity of the systems, the accretion disc’s wind velocity, the fraction of flux thermalisation in the disc, and so on. We present those results here.

Figure 3 presents the assumed SEDs for three BH-ULXs and three NS-ULXs generated from our simulations (see Fig. 3 of Vinokurov et al. 2013). The spectra are calculated at a distance of 3 Mpc using Eq. (14) without the contribution from the stellar companion. The figure also includes the key parameters for each system: M and ṁ. In Fig. A.1¹, we show the effect of varying ṁ for a constant M, where it is evident that the 2 keV X-ray flux is independent of ṁ, a crucial point made above for the α_ox and L_ν,UV anti-correlation. Other studies have shown that the ULX wind velocity can be as high as 0.1c − 0.3c (Pinto & Walton 2023; Pinto et al. 2020; Kosec et al. 2018). Increasing the wind velocity from 1000 km s⁻¹ to 60 000 km s⁻¹ (0.2c) leads to a reduction in flux at lower frequencies and the disappearance of the UV shoulder. However, the X-ray emission at higher frequencies is unaffected by the change in wind velocity (see Fig. A.2). This results in more positive α_ox values and a steeper α_ox − L_ν,UV relation (see Fig. A.4).

Fig. 3. Spectra of three BH-ULX and three NS-ULX systems from our simulation, showing the mass of the compact object (M) and the ratio ṁ = Ṁ/ṀEdd. Flux is calculated for all systems located at 3 Mpc, and for all curves fout = 0.03 and θf = 45°.

Figure 4 shows the index α_ox, defined in Eq. (18), as a function of UV luminosity, combining the accretion disc and the companion star, at different metallicities (2.5%, 20%, and 40% of solar metallicity from top to bottom) of our simulated ULX population. We fit the BH- and NS-ULX simulated populations together and separately with the function

Fig. 4. αox index as function of Lν,UV (2500 Å) luminosity for our simulated ULX population (BH or NS primaries). The best fit to αox − Lν,UV simulated data in our population (black solid line) and observed ULXs data in Sonbas et al. (2019) (black points) are also shown. The top, middle, and bottom panels correspond to metallicities of 0.025 Z⊙, 0.20 Z⊙, and 0.40 Z⊙, respectively. We also show the αox − Lν,UV fits for the BH-ULX and NS-ULX systems separately.

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{OX}}=m{log}_{10}{L}_{\nu ,\mathrm{U}V}+C,\end{array}$$(19)

where the slope m and intercept C are free parameters. The slopes from these fits are shown in Table 2.

The results presented in Table 2 reveal that the slope m of the α_ox − L_ν,UV relation varies with metallicity and across different ULX populations. For BH-ULXs, the slope becomes less negative as metallicity increases, indicating that higher metallicities reduce the sensitivity of the X-ray to UV spectral index to UV luminosity. In contrast, NS-ULXs exhibit a relatively stable slope around −0.33 across all metallicities. This consistency implies that the SED of NS-ULXs is less affected by changes in metallicity compared to BH-ULXs. The combined population shows a similar trend to BH-ULXs, with the slope decreasing slightly with increasing metallicity. Combined slopes are closer to the BH-ULX slopes, indicating that BH-ULXs may play a slightly more dominant role in determining the overall behaviour. This is in line with the simulated BH-ULXs outnumbering the NS-ULXs in our simulated populations, as evident from Fig. 1.

Also shown in Fig. 4 are the observational data from Sonbas et al. (2019). We note that the α_ox values for the observed ULXs are described better by the BH-ULX systems with f_out = 0.03, except, possibly, for two cases that can be explained by NS-ULX systems. The metallicities of the host galaxies of the ULXs listed in Sonbas et al. (2019) are between 0.10 Z_⊙ and 0.36 Z_⊙, with an average metallicity of 0.22 Z_⊙. This is approximately the value corresponding to the middle panel of Fig. 4.

Figure 4 shows that most of the NS-ULX systems in our simulation are in the lower range of UV luminosities compared to the data points from Sonbas et al. (2019) with f_out = 0.03. The range of L_ν, UV is affected by the choice of f_out. Increasing the f_out value leads to high L_ν, UV (see Sect. 4.3). This suggests a potential observational bias, where the observed sample might preferentially include ULXs with higher UV fluxes. This bias could occur because systems with lower UV fluxes are harder to detect and may be under-represented in the observational sample, leading to a steeper observed slope. Furthermore, our synthetic population includes systems with higher UV luminosities not present in the observational data, indicating that our simulations predict a broader range of UV luminosity for ULX systems.

4.2.3. Disc versus stellar-companion contribution to UV emission

We calculated UV emission both from the accretion disc and the stellar companion star. We find that a higher percentage of ULX systems are dominated by the disc component in the UV, which is consistent with observational studies (Tao et al. 2012). The percentage of disc-dominated ULXs increases with metallicity, starting from 77% at the lowest metallicity simulated (0.005 Z_⊙) to 93% at the highest metallicity simulated (1.5 Z_⊙). The percentages of ULXs with disc-dominated UV emission in Fig. 4 are 85%, 82%, and 86% from the top to the bottom panel, respectively. We also find that the majority of the NS-ULXs are dominated by the disc in the UV, leaving the slope of the companion-dominated UV systems being close to that of the BH-ULXs. We note that the companion stars that do dominate the total UV emission are primarily OB-type stars with very high surface temperatures (from 10 000 K to 50 000 K).

Figure 5 shows the α_ox − L_ν,UV relation at 0.2 Z_⊙. We note here that this is specifically for the binaries with L_UV dominated by the companion. The donor-dominated BH-ULX systems display a stronger dependence on metallicity, while NS-ULX systems remain relatively unaffected by changes in metallicity (see Table 3).

Fig. 5. αox for donor-star-dominated UV emission for ULX systems as a function of UV (2500 Å) luminosity for our simulated ULX population (BH or NS primaries). The best fit to αox − Lν,UV simulated data in our population (black solid) at 0.2 Z⊙ and that for observed ULXs in Sonbas et al. (2019) (black points) are also shown.

4.3. Impact of free parameters f_out and θ_f on α_ox

In this model, the key free parameters, f_out and θ_f, play crucial roles in shaping the reprocessing of disc radiation and the observed spectra. The parameter f_out was determined by the average value in systems explored by Vinokurov et al. (2013). In reality, every system will have a different value of f_out. Unfortunately, we do not have a way to determine the value of f_out for each system; however, we can study the effect of varying this parameter in our population properties.

The UV emission is highly sensitive to the value of f_out; higher values result in greater UV contribution (see the SED in Fig. A.3). Increasing f_out makes the values of α_ox more negative and the slope of α_ox − L_ν, UV relation becomes closer to that of Sonbas et al. (2019) for the combined population (see Fig. A.5). With f_out = 0.2, the majority of the observed (Sonbas et al. 2019) ULXs can be explained by both NS-ULXs and BH-ULXs from our simulations.

From Fig. 3, it can be observed that the frequency corresponding to 2500 Å does not lie within the flat part of the UV and optical SED for those systems. The wind temperature scales as $f_{\mathrm{out}}^{1/4}$, and decreasing f_out shifts the flat portion of the spectrum towards higher frequencies (or, equivalently, shorter wavelengths). Hence, to place 2500 Å within the flat part of the spectrum, f_out must be decreased further.

The parameter θ_f primarily influences the X-ray emission, making α_ox values more positive with decreasing θ_f, though it has a negligible impact on the slope of the α_ox − L_ν,UV relation for a fixed f_out. Decreasing the funnel angle from 45° to 5° represents a tenfold increase in luminosity. With θ_f = 45°, the most luminous BH-ULX and NS-ULX at the lowest metallicity have luminosities of L_X = 1.4 × 10⁴¹ erg s⁻¹ and L_X = 5.4 × 10³⁹ erg s⁻¹, respectively, while the most luminous systems have L_X = 4.2 × 10⁴⁰ erg s⁻¹ and L_X = 5.5 × 10³⁹ erg s⁻¹ at the highest metallicity simulated. To fully explore the effects of varying θ_f, we calculated the funnel angle θ_f based on the accretion rate using the model proposed by King (2009) in Sect. 4.3.1.

4.3.1. Results with beaming and comparative analysis

In the initial phase of our study, we assumed a fixed funnel angle of θ_f = 45°, allowing us to explore the fundamental properties of ULXs and their X-ray luminosity under this assumption. However, this assumption does not fully capture the dynamic nature of accretion processes, particularly at high accretion rates where the emission is expected to be highly beamed (King et al. 2001). According to King et al. (2001) and later refined by King (2009), the funnel angle θ_f varies with the accretion rate ṁ according to Eq. (9). When the accretion rate is high, the luminosity is focused into narrower cones, significantly enhancing the observed luminosity. The beaming factor b is defined as

$$\begin{array}{c}\hfill b=\frac{\mathrm{\Omega }}{4\pi }=1-cos{\theta }_f,\end{array}$$(20)

where Ω is the solid angle of emission. For a conical beam with an opening angle of θ_f, the solid angle is given by

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{total}}=4\pi (1-cos{\theta }_f).\end{array}$$(21)

We now calculate the funnel angle θ_f based on the accretion rate using the model proposed by King (2009). This equation allows us to compute θ_f directly from the beaming factor as

$$\begin{array}{c}\hfill {\theta }_f={cos}^{-1}(1-b).\end{array}$$(22)

For b = 1, θ_f = 90°, and r_ph tends towards infinity. To avoid this, we set θ_f, max = 87°, preventing r_ph from becoming infinitely large. With a maximum value of ṁ capped at 1000, the minimum achievable value of b is 7 × 10⁻⁵. However, our simulations show a minimum b value of 10⁻⁴. Previous studies have suggested a saturation in the geometric beaming due to a high accretion rate, indicating that beyond a certain mass transfer rate (Ṁ), further increases in Ṁ do not significantly impact the scale height of the accretion disc (Lasota et al. 2016; Wiktorowicz et al. 2017). We note that the SCAD model does not put a limit on the geometric thickness of the disc and hence the funnel angle θ_f.

The distribution of b values from our simulations is illustrated in Fig. A.6. We find that the most luminous NS-ULX achieves L_X, b = 5.3 × 10⁴² erg s⁻¹, while the most luminous BH-ULX reaches L_X, b = 7.5 × 10⁴³ erg s⁻¹ at the lowest metallicity simulated, assuming an observer is always looking down the funnel.

Figure 6 shows the cumulative XLF of the simulated ULX population at 40% Z_⊙ in the 0.2−12 keV band, incorporating the effects of beaming. The comparison between beamed and non-beamed scenarios reveals that beaming has a significant impact on the luminosity distribution of ULXs. The most luminous sources become even more prominent due to the narrowing of the funnel angle at higher accretion rates, which focuses the emission into a smaller solid angle and therefore, increases the observed luminosity if the source is viewed down the wind funnel. The power-law index α in the beamed scenario at 40% of Z_⊙ is 1.50, compared to 1.99 without beaming. This indicates that beaming increases the number of high-luminosity ULXs, shifting the XLF towards higher luminosities and flattening the slope. This flattening is a direct result of the increased apparent luminosity due to beaming in systems with high accretion rates.

Fig. 6. Cumulative XLF of simulated ULX population at 40% of Z⊙ in the 0.2−12 keV band.

Walton et al. (2011) constructed XLFs for observed ULXs in spiral and elliptical galaxies and found α = 1.85 ± 0.11 and a steeper index of α = 2.5 ± 0.4 for the two galaxy types, respectively. Salvaggio et al. (2023) studied the XLF of ULXs found in the Cartwheel galaxy and found a slope of $\alpha =1.{76}_{-0.43}^{+0.51}$. Our simulated power-law indices (beamed and non-beamed) fall within the uncertainties of the indices reported for the observed data by these authors.

Table 4 presents the α, L_s, and β in Eq. (17) for XLF of ULXs at different metallicities, including the effects of beaming. These parameters highlight the significant impact of beaming on the observed properties of ULX populations. The slope α, which characterises the distribution of ULX luminosities, shows a notable difference between the beamed and non-beamed cases. With a constant θ_f = 45° (Table 1), the slope α is consistently steeper across all metallicities – ranging from 1.77 to 2.15 – indicating a rapid decline in the number of high-luminosity ULXs. In contrast, with beaming (Table 4), the slope remains nearly constant across all metallicities.

The dependence of N(L_X) on the beaming angle θ_f flattens the slope of the XLF and reduces its sensitivity to metallicity. Hence, beaming masks the underlying variations in mass distributions across metallicities that would otherwise steepen the XLF, as higher metallicities tend to produce lower mass compact objects in the case of black holes.

With beaming, L_s values are generally higher. For instance, at 0.5% of Z_⊙, L_s = 218 × 10³⁹ erg s⁻¹, which is much higher than L_s = 59.5 × 10³⁹ erg s⁻¹ with θ_f = 45° and at 0.5% of Z_⊙.

The introduction of beaming also impacts the optical–X-ray spectral index, α_ox. Figure 7 presents α_ox as a function of L_ν,UV (2500 Å luminosity) for our simulated ULX population, including beaming. The slopes (m) of the α_ox versus L_ν,UV relation differ significantly between the beamed (Table 5) and non-beamed ones (Table 2) at high metallicity, highlighting the substantial impact of beaming on the relative X-ray and UV emissions of ULXs.

Fig. 7. αox index as function of Lν,UV (2500 Å) luminosity for our simulated ULX population (BH or NS primaries) with beaming included in the calculation. The best fit to αox − Lν,UV simulated data in our population (black solid line) and observed ULX data in Sonbas et al. (2019) (black points) are also shown. The top, middle, and bottom panels correspond to metallicities of 0.25 Z⊙, 0.2 Z⊙, and 0.4 Z⊙, respectively.

Including beaming leads to high-beamed systems resulting in more positive α_ox values compared to the observational data and, hence, more data points lying above those of Sonbas et al. (2019). In general, introducing beaming reduces the magnitude of the slope, m, particularly at higher metallicities. This suggests that at higher metallicities, beaming results in a larger number of systems with focused emission, leading to a flattening of the α_ox − L_ν,UV relation. As a result, metallicity plays a diminished role in determining the observed correlation, with beaming effects becoming the dominant factor influencing the properties of these ULX systems. Furthermore, we notice that the inclusion of beaming leads to a large scatter in the α_ox − L_ν,UV plane and a linear fit may not be appropriate as can be inferred from the near-flat slope of the fits.

4.3.2. Neutron star ULXs

Pulsations found in some ULX systems have provided compelling evidence that at least a fraction of ULXs are powered by NSs (Bachetti et al. 2014; Fürst et al. 2016). In Fig. 8, we present the XLF of our simulated NS-ULX population at Z = 40% Z_⊙. Our simulations reveal that the XLF of NS-ULXs has a slope similar to that of the total ULX population, though with slightly higher L_s and β; values in Fig. 8 are given by 381 ± 29 × 10³⁹ erg s⁻¹ and 1.3 ± 0.1, respectively.

Fig. 8. Cumulative XLF of simulated NS-ULX population at 40% of Z⊙ in the 0.2−12 keV band including beaming.

We find that the average value of θ_f ≈ 22° for all NS-ULXs at Z = 40% Z_⊙, which corresponds to 1/b ≈ 13. Therefore, a beaming factor of at least 10 is required to achieve extreme luminosities (L_X > 10⁴⁰ erg/s) observed in some NS-ULXs. However, the high pulsed fractions observed in some pulsating ULXs (Bachetti et al. 2022) give a strong argument against strong beaming in these systems (Mushtukov et al. 2021). Furthermore, Bachetti et al. (2022) showed that the extreme luminosity of M82 X-2 can be explained by super-Eddington accretion and high mass-transfer rates onto a highly magnetised NS.

The omission of magnetic fields from our model presents a limitation. Including the magnetic field can significantly affect the accretion geometry and emission properties of the NS-ULX systems (e.g. Middleton et al. 2019; Mushtukov et al. 2021; Bachetti et al. 2022). Incorporating these effects is, however, beyond the scope of this work.

4.3.3. Observable ULX population

Beamed emission from a ULX is observable if the observer’s line of sight is within θ_f. To study what fraction of our ULX population is observable, we assigned a random viewing angle (β) to all ULXs using the distribution P(β) = sin(β). The cumulative distribution function (CDF) of P(β) is then

$$\begin{array}{c}\hfill C(\beta )={\displaystyle {\int }_0^{\beta }}P(\beta \prime )\mathrm{d}\beta \prime ={\displaystyle {\int }_0^{\beta }}\frac{1}{2}sin\beta \prime \mathrm{d}\beta \prime =\frac{1}{2}(1-cos\beta ).\end{array}$$(23)

Here, C(β) ranges from 0 (at β = 0) to 1 (at β = π) and β = 0 in a face-on case, and β = π/2 is an edge-on case. Therefore,

$$\begin{array}{c}\hfill \beta =arccos(1-2C).\end{array}$$(24)

We used Eq. (24) to calculate β, where C is randomly sampled from a uniform distribution in [0, 1]. A ULX is observable if θ_f > β or θ_f > π − β.

Table 6 shows the percentage of ULX population observable after viewing angle correction. A relatively small fraction (≲10%) of the NS ULXs are observable as compared to the BH ULXs. Notably, the NS-ULX fraction remains largely independent of metallicity, while the BH ULXs exhibit a much higher fraction, ranging from 43% at 0.5% Z_⊙ to 78% at 150% Z_⊙. This difference in the observability of the NS ULXs and BH ULXs can be understood from their beaming distributions shown in Fig. A.7.

Although the fraction of observable ULXs remains below 44% across all metallicities, their properties are not substantially different from those of the entire simulated population. Figure 9 displays the luminosity function of the observable ULXs at 40% Z_⊙, while Fig. 11 presents their α_ox − L_ν,UV relation. The observable population shows a slightly steeper α_ox − L_ν,UV slope. It is important to note, however, that our analysis of the observable population does not account for any plausible selection effects of the detector.

Fig. 9. Cumulative XLF of simulated observable ULXs at 40% of Z⊙ in the 0.2−12 keV band including beaming.

To investigate the effect of selection criteria on the XLF, we computed the selection probability as a function of luminosity. To calculate this probability, we divided the luminosity range into logarithmically spaced bins and computed the total number of systems and the number of observable systems in each bin. The selection probability is then defined as the ratio of the observable to the total number of systems.

Figure 10 shows the selection probability for the population in Fig. 6 and the observable population in Fig. 9. The results reveal that the probability remains approximately flat between 10³⁹ and 2 × 10⁴⁰ erg s⁻¹, and this is where most of the ULX systems are located.

Fig. 10. Selection probability, ratio of the observable to the total number of all ULX systems (blue), BH-ULX (green), and NS-ULX (red), as a function of X-ray luminosity at 40% Z⊙.

Fig. 11. αox versus Lν,UV (2500 Å) for our simulated observable ULX population (BH or NS primaries) with beaming included in the calculation. The best fit to the simulated data (black solid line) and observed ULXs in Sonbas et al. (2019) (black points) are also shown. The top, middle, and bottom panels correspond to metallicities of 0.25 Z⊙, 0.2 Z⊙, and 0.4 Z⊙, respectively.

The invariance of the XLF slope is primarily driven by the dominance of ULX systems within the flat-probability range in Fig. 10. In this region, the probability of making an observation is high due to low beaming, which results in a wider opening angle. At high luminosities (L > 2 × 10⁴⁰ erg s⁻¹), the NS-ULX emission becomes increasingly beamed into narrow cones, significantly reducing the likelihood of observation and resulting in a drop in the total probability. In contrast, the BH-ULX selection probability remains constant below luminosities of L ∼ 10⁴¹ erg s⁻¹, where beaming is unimportant.

The selection process does not introduce significant luminosity-dependent biases in this range, allowing the intrinsic XLF properties to be preserved in the observed sample. It is also important to remind the reader that the viewing angle is randomly assigned from a sin(β) distribution. Hence, the luminosity does not depend on the angle β.

5. Summary and conclusions

We summarise our results below.

1. In this study, we used the COSMIC population synthesis code to simulate the evolution of binary systems and investigate the relationship between UV and X-ray emission during the ULX phase.

2. We adopted the SCAD model by Vinokurov et al. (2013) to calculate the SEDs of our simulated ULX population from the accretion disc. While X-rays are emitted solely by the accretion disc, the UV emission is a combination of heated wind from the disc and the surface emission from the companion star.

3. Our findings indicate that the number of ULXs decreases with increasing metallicity, which is consistent with observational data. We found that BH-ULXs decrease at a faster rate than NS-ULXs with increasing metallicity and obtained a power-law index of α = 0.16 ± 0.01 for the combined population of BH- and NS-ULXs, which is consistent with values reported in the literature. The decreasing number of ULXs with increasing metallicity has significant implications for our understanding of ULX populations in different galactic environments. In metal-poor galaxies, we expect a higher number of ULXs, which can serve as important probes of early galaxy evolution.

4. Luminosity Function: the XLF slopes of our simulated ULX population vary with metallicity; that is, with lower metallicity populations exhibiting a relative excess of high-luminosity ULXs. Without beaming, the slope ranges from α = 1.77 ± 0.01 at 0.5% solar metallicity to α = 2.15 ± 0.38 at 1.5 times solar metallicity. The inclusion of beaming flattens the XLF slope across all metallicities, reflecting the boosted luminosity due to the beaming effect. Our simulated slopes are consistent with observational data from Walton et al. (2011) and Salvaggio et al. (2023), who reported similar values for ULXs in different types of galaxies.

5. The optical–X-ray spectral index α_ox showed a significant anti-correlation with UV luminosity across different metallicities. This relation is found using both the disc and the stellar companion contributions to the UV emission. Our results suggest that UV emission in ULXs is predominantly disc-dominated, with the percentage of disc-dominated ULXs increasing as metallicity rises. The slope of the α_ox − L_ν,UV relation for BH-ULXs changes with metallicity; it becomes steeper at higher metallicities, while the slope for NS-ULXs remains relatively constant. With beaming included, the slopes of the α_ox − L_ν,UV relation become significantly flatter and remain relatively constant with metallicity.

6. For a constant θ_f, our results suggest that the majority of the ULXs with unique optical counterparts selected by Sonbas et al. (2019) are likely powered by black holes when f_out = 0.03; however, for high f_out > 0.1 more NS-ULX systems are found to be compatible with the observational data. Including beaming in our calculation resulted in high α_ox values compared to Sonbas et al. (2019) data points.

7. Our results support the hypothesis that stellar-mass BHs and NSs, undergoing supercritical accretion, can account for the observed properties of ULXs without necessitating the presence of IMBHs.

8. While the X-ray emission remains unaffected, higher wind velocity reduces UV emission, causing the α_ox values to become more positive.

9. Due to their high beaming, NS-ULX systems produce narrow emission cones, which significantly reduce their likelihood of being observed.

We note that COSMIC has a few limitations in simulating the evolution of binary populations. For example, the relationship between remnant mass and initial mass or carbon-oxygen core mass is not always straightforward. Other theoretical models have shown that some massive stars directly implode into black holes without an explosion (Patton & Sukhbold 2020), which COSMIC does not fully account for. Unlike advanced stellar evolution codes such as Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011), which provides detailed simulations of the physical processes governing the life cycle of stars, COSMIC is a rapid population-synthesis code that relies on simplified models for stellar and binary evolution, introducing systematic biases in the resulting binary populations (Marchant et al. 2021; Gallegos-Garcia et al. 2021). These limitations can be addressed by population synthesis codes such as POSYDON (Fragos et al. 2023), which incorporate a detailed MESA code. However, current versions of these codes are constrained to a limited range of metallicity and are not yet suitable for studies similar to ours, where metallicity is varied over a wide range.

We also note a few limitations of our radiation calculations. The SCAD model has been developed for BH-ULX systems. Its application to NS-ULX systems in our study is limited. One significant issue is that we did not account for the effects of magnetic fields in NS-ULX systems. Magnetic fields can play a critical role in the dynamics of the accretion disc, potentially leading to disc truncation and the occurrence of pulsations, both of which are characteristic features of NS-ULXs (Middleton et al. 2019). Calculations by Mushtukov et al. (2021) suggest that strong magnetic fields (B > 10¹³ G) in systems such as M82 X-2 facilitate super-Eddington accretion and reduce reliance on beaming to achieve high luminosities. Similarly, Bachetti et al. (2022) argued that moderate beaming combined with accretion funnel effects and magnetic-field interactions better explain the observed properties of pulsating ULXs. Additionally, the irradiation of the donor star by the accretion disc, which can increase UV emission (Motch et al. 2014), was not considered, potentially affecting the accuracy of the α_ox versus L_ν,UV relation. Our future studies will aim to improve on some of those limitations.

Data availability

Appendix A, which includes Figs. A.1 through A.7, is available on Zenodo: https://doi.org/10.5281/zenodo.14715088

Acknowledgments

We are grateful to the referee for helpful comments that have improved this paper. We would like to thank Kalvir Dhuga, Eda Sonbas and Andrzej Zdziarski for helpful discussions. L.N. is grateful to Katelyn Breivik for advice on the COSMIC code. L.N. and S.R. were partially supported by a BRICS STI grant and by a NITheCS grant from the National Research Foundation, South Africa. J.D.F. was supported by NASA through contract S-15633Y and the Office of Naval Research.

References

All Tables

All Figures

	Fig. 1. Number of ULXs formed in our simulation (NULX) as a function of metallicity. The y-axis indicates the number of ULXs derived from the simulations, with emphasis placed on the trend’s shape rather than the absolute values. The solid lines represent fits to the simulated data, with the corresponding slopes mentioned in the legend.
In the text

	Fig. 2. Cumulative XLF of the simulated ULX population at 40% Z⊙ in the 0.2−12 keV band.
In the text

	Fig. 3. Spectra of three BH-ULX and three NS-ULX systems from our simulation, showing the mass of the compact object (M) and the ratio ṁ = Ṁ/ṀEdd. Flux is calculated for all systems located at 3 Mpc, and for all curves fout = 0.03 and θf = 45°.
In the text

Fig. 4. αox index as function of Lν,UV (2500 Å) luminosity for our simulated ULX population (BH or NS primaries). The best fit to αox − Lν,UV simulated data in our population (black solid line) and observed ULXs data in Sonbas et al. (2019) (black points) are also shown. The top, middle, and bottom panels correspond to metallicities of 0.025 Z⊙, 0.20 Z⊙, and 0.40 Z⊙, respectively. We also show the αox − Lν,UV fits for the BH-ULX and NS-ULX systems separately.

In the text

	Fig. 5. αox for donor-star-dominated UV emission for ULX systems as a function of UV (2500 Å) luminosity for our simulated ULX population (BH or NS primaries). The best fit to αox − Lν,UV simulated data in our population (black solid) at 0.2 Z⊙ and that for observed ULXs in Sonbas et al. (2019) (black points) are also shown.
In the text

	Fig. 6. Cumulative XLF of simulated ULX population at 40% of Z⊙ in the 0.2−12 keV band.
In the text

Fig. 7. αox index as function of Lν,UV (2500 Å) luminosity for our simulated ULX population (BH or NS primaries) with beaming included in the calculation. The best fit to αox − Lν,UV simulated data in our population (black solid line) and observed ULX data in Sonbas et al. (2019) (black points) are also shown. The top, middle, and bottom panels correspond to metallicities of 0.25 Z⊙, 0.2 Z⊙, and 0.4 Z⊙, respectively.

In the text

	Fig. 8. Cumulative XLF of simulated NS-ULX population at 40% of Z⊙ in the 0.2−12 keV band including beaming.
In the text

	Fig. 9. Cumulative XLF of simulated observable ULXs at 40% of Z⊙ in the 0.2−12 keV band including beaming.
In the text

	Fig. 10. Selection probability, ratio of the observable to the total number of all ULX systems (blue), BH-ULX (green), and NS-ULX (red), as a function of X-ray luminosity at 40% Z⊙.
In the text

	Fig. 11. αox versus Lν,UV (2500 Å) for our simulated observable ULX population (BH or NS primaries) with beaming included in the calculation. The best fit to the simulated data (black solid line) and observed ULXs in Sonbas et al. (2019) (black points) are also shown. The top, middle, and bottom panels correspond to metallicities of 0.25 Z⊙, 0.2 Z⊙, and 0.4 Z⊙, respectively.
In the text