# The Observed versus Total Population of ULXs

Grzegorz Wiktorowicz, Jean-Pierre Lasota, Matthew Middleton, Krzysztof Belczynski

► To cite this version:

Grzegorz Wiktorowicz, Jean-Pierre Lasota, Matthew Middleton, Krzysztof Belczynski. The Observed versus Total Population of ULXs. The Astrophysical Journal, 2019, 875 (1), pp.53. 10.3847/1538-4357/ab0f27. [hal-01952653](https://hal.science/hal-01952653v1)

HAL Id: hal-01952653

https://hal.science/hal-01952653v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Observed versus Total Population of ULXs

Grzegorz Wiktorowicz¹,², Jean-Pierre Lasota³,⁴, Matthew Middleton⁵, and Krzysztof Belczynski³

¹National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, People's Republic of China; gwiktoro@astruw.edu.pl

²School of Astronomy & Space Science, University of the Chinese Academy of Sciences, Beijing 100012, People's Republic of China

³Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warsaw, Poland

⁴Institut d'Astrophysique de Paris, CNRS et Sorbonne Université, UMR 7095, 98bis Bd Arago, F-75014 Paris, France

⁵Department of Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ, UK

Received 2018 November 21; revised 2019 March 6; accepted 2019 March 11; published 2019 April 16

## Abstract

We have analyzed how anisotropic emission of radiation affects the observed sample of ultraluminous X-ray sources (ULXs) by performing simulations of the evolution of stellar populations, employing recent developments in stellar and binary physics, and by utilizing a geometrical beaming model motivated by theory and observation. While ULXs harboring black hole accretors (BH ULXs) are typically emitting isotropically, the majority of ULXs with neutron star accretors (NS ULXs) are found to be beamed. These findings confirm previous assertions that a significant fraction of ULXs are hidden from view due to a substantial misalignment of the emission beam and the line of sight. We find the total number of NS ULXs in regions with constant star formation (SF), solar metallicity, and ages above ~1 Gyr to be higher than the BH ULXs, although observationally both populations are comparable. For lower metallicities, BH ULX dominate both the total and observed ULX populations. As far as burst SF is concerned, young ULX populations are dominated by BH ULXs, but this changes as the population ages and, post SF, NS ULXs dominate both the observed and total ULX populations. We also compare our simulation output to a previous analytical prediction for the relative ratio of BH to NS ULXs in idealized flux-limited observations and find broad agreement for all but the lowest metallicities. In so doing, we find that in such surveys the observed ULX population should be heavily dominated by black hole systems rather than by systems containing neutron stars.

**Key words:** stars: black holes – stars: neutron – methods: statistical – X-rays: binaries

## 1. Introduction

Geometrical beaming occurs in an accreting system when radiation preferentially escapes along a beam with an opening solid angle <4π steradians. As a result, an observer located in the cone of emission will infer a higher luminosity by assuming isotropic emission, than the real total integrated luminosity of the system. There are observational and theoretical arguments for the presence of beaming in systems with very high accretion rates onto a compact object. For example, despite its very high accretion rate ($\dot{M} \approx 10^{-4} M_{\odot} yr^{-1}$ Fabrika 2004), the binary system SS 433 is faint in the X-rays ($L_X \approx 10^{36} erg s^{-1}$) with recent evidence (Middleton et al. 2018) indicating that most of the radiation escapes at high inclinations to the line of sight (similar to the jets, Begelman et al. 2006; Medvedev & Fabrika 2010).

From a theoretical point of view, basic physics and detailed numerical calculations indicate that at super-Eddington accretion rates, the disk ceases to be geometrically thin around the spherisation radius ($R_{sph}$, e.g., Shakura & Sunyaev 1973; Ohsuga et al. 2009; Ohsuga & Mineshige 2011; Sądowski et al. 2014) where radiation pressure dominates and the Eddington limit is reached locally. The location of $R_{sph}$ is expected to be linearly proportional to the accretion rate, $R_{sph} \sim \dot{m}$, where $\dot{m}$ is the mass accretion rate in Eddington units, $\dot{m}_{Edd} = L_{Edd}/\eta c^2$, with a radiative efficiency, $\eta \approx 0.1$, and an Eddington luminosity ($L_{Edd}$) for an accretor mass ($M_{acc}$) and hydrogen abundance in the accretion flow ($X$):

$$L_{Edd} = 2.6 \times 10^{38} \frac{1}{1+X} \frac{M_{acc}}{M_{\odot}} \left[ \frac{\text{erg}}{\text{s}} \right] \quad (1)$$

Due to the large aspect ratio of the disk and ease with which material is lost in a wind (needed to keep the accretion rate at the Eddington value for smaller radii in the absence of advection), emission from the inner-most regions (where the most energetic photons are formed) is trapped in a conical, optically thick structure.

Ultraluminous X-ray sources (ULXs) are defined as point-like, off-nuclear sources with—isotropically assumed—observed X-ray luminosities above $L_{ULX} = 10^{39} erg s^{-1}$ (for a recent review, see Kaaret et al. 2017) and are particularly important in the context of super-Eddington accretion. One interpretation involves sub-Eddington accretion onto intermediate-mass BHs (Colbert & Mushotzky 1999). However, King et al. (2001) showed that globular clusters on average cannot produce the necessary number of IMBHs to explain all ULXs and, instead, argue that only the presence of beaming in ULXs avoids serious formation difficulties. Observationally, evolution in the X-ray spectra of ULXs, coupled with the short timescale variability would also argue for geometrical beaming in a super-critical flow (Middleton et al. 2015a; Middleton & King 2016).

The population synthesis of ULXs has been performed in several studies. Rappaport et al. (2005) showed that population of ULXs in spiral galaxies can be explained by short-high mass-transfer (MT) phases in BH+MS binaries. However, they neglected pre-supernova evolution in their calculations and did not take into account other types of binaries, e.g., containing NS accretors, or evolved donors. Their study was expanded in Madhusudhan et al. (2008) where they additionally predicted the observational properties of the ULX population. Meanwhile, intermediate-mass BHs were proposed as potential (if
hypothetical) accretors in ULXs for which a violation of the Eddington limit is not necessary (Colbert & Mushotzky 1999; Madhusudhan et al. 2006). More recently, Linden et al. (2010) utilized the StarTrack population synthesis code to show that the bulk of ULXs can be explained as a high-luminosity tail of high-mass X-ray binaries. Since the recent discovery of NS in ULXs (Bachetti et al. 2014), NS accretors were included in studies of ULXs. For these objects, the Eddington limit is apparently surpassed more than several times (e.g., Fragos et al. 2015; Shao & Li 2015; Wiktorowicz et al. 2015). Similarly, the detection of double compact object mergers (Abbott et al. 2016) triggered an investigation of potential connections between double compact objects and ULXs (e.g., Finke & Razzaque 2017; Marchant et al. 2017; Klencki et al. 2018). Massive stars (mainly red super-giants) were detected in optical and infrared bands as potential donors in a few ULXs (e.g., Kaaret et al. 2004; Liu et al. 2004; Heida et al. 2014, 2016), whereas Wiktorowicz et al. (2017) predicted main-sequence donors (typically 5.9–11 $M_\odot$ for BH accretors and 0.9–1.5 $M_\odot$ for NS accretors) for the majority of ULXs.

The discovery of pulsing ULXs (PULX; Bachetti et al. 2014) has called into question the role of beaming versus strong dipole magnetic fields (e.g., Mushtukov et al. 2015b), but observations of cyclotron resonance lines (Brightman et al. 2018; Walton et al. 2018; but see Koliopanos et al. 2019) would broadly support “normal”, pulsar-like dipole field strengths $10^{11} \lesssim B \lesssim 10^{13}$ G. The combination of such field strengths and super-Eddington accretion rates allows for a consistent—if not complete—description of the observed properties of PULXs (King & Lasota 2016, 2019; King et al. 2017; M. Middleton et al. 2019, in preparation). In particular, King & Lasota (2016) found that neutron star ULXs are likely to have higher apparent luminosities than black hole ULXs (BHULXs) for a given MT rate, as their increased beaming outweighs their lower Eddington luminosities. For example, using methods provided in Section 2, for a typical NS and BH mass ($1.4 M_\odot$, and $7 M_\odot$, respectively e.g., Özel et al. 2010) and an MT rate $\dot{M} = 10^{-5} M_\odot$ yr$^{-1}$, the real, total integrated luminosity ($L_X$) is higher for the BH ($\sim 5.3 \times 10^{39}$ erg s$^{-1}$) than for a NS ($\sim 1.4 \times 10^{39}$ erg s$^{-1}$). However, the apparent luminosity for an observer located in the emission cone is higher for the NS ($\sim 4.3 \times 10^{41}$ erg s$^{-1}$) than for the BH ($\sim 2.0 \times 10^{41}$ erg s$^{-1}$). We note that some GRMHD simulations show the opposite result. For example, Abarca et al. (2018) performed a simulation of super-Eddington accretion ($\dot{M} = 200 L_{Edd}/c^2$) onto a non-magnetized non-rotating neutron star ($M_{NS} = 1.4 M_\odot$) and, contrary to the BH case, found no significant beaming. In this model, the luminosity is about $L_X \approx 10^{38}$ erg s$^{-1}$ and pulsations would not be visible, and so it is not directly related to PULXs. Although magnetar-strength magnetic fields may also be responsible for the emergence of super-Eddington levels of radiation (e.g., Basko & Sunyaev 1976; Mushtukov et al. 2015a) and strongly influence the mass flow in the accretion disk (e.g., Parfrey & Tchekhovskoy 2017), in this paper we focus on geometrical beaming only. We are encouraged to do so by the fact that all known magnetars are single (Olausen & Kaspi 2014) and although they are supposed to be formed in binaries, their birth leads to the binary orbit disruption (see, e.g., Clark et al. 2014). The rare survivors are unlikely to be present in PULXs (Popov 2016).

For a given funnel opening angle $\theta$, the probability of an observer being located in its cone of emission is given by

$$P_{\text{obs}}(\theta) = 1 - \cos \theta/2. \qquad (2)$$

Therefore, the stronger the beaming (lower $\theta$), the lower the fraction of observed systems in the total population. On the other hand, beamed sources may be visible from much larger distances due to higher apparent luminosities (see Middleton & King 2017, Section 4).

In this paper, we analyze the impact of beaming on the relation between the observed and total sample of ULXs. Our calculations are based on results presented in Wiktorowicz et al. (2017) where beaming was already included, but not analyzed in detail. Our main motivations are the recent discoveries of NSs in ULXs (Bachetti et al. 2014; Fürst et al. 2016; Israel et al. 2017; Carpano et al. 2018) and observational hints that many non-pulsing ULXs may host NSs (e.g., Koliopanos et al. 2017; Mushtukov et al. 2017; Pintore et al. 2017; Walton et al. 2018) as predicted by previous works, e.g., King et al. (2001) and King & Lasota (2016).

## 2. Methods

In Wiktorowicz et al. (2017), ULX populations in different environments were analyzed; however, that work focused only on observed ULXs, i.e., those that are predicted to be visible from the Earth. The total population of ULXs (including so called “misaligned” or “hidden” sources such as SS433: Middleton et al. (2018), or MQ1 in M83: Soria et al. (2014) was not analyzed.

We utilized the Startrack population synthesis code (Belczynski et al. 2002, 2008) with further updates (see Wiktorowicz et al. 2017, and references therein). For the initial primary masses, we used the Kroupa initial mass function with $P(M_{ZAMS}) \propto M_{ZAMS}^{\Gamma}$ (Kroupa & Weidner 2003) across a range $5-150 M_\odot$. The power-law index, $\Gamma = -1.3$ for stars with $M_{ZAMS} \le 0.5 M_\odot$, $\Gamma = -2.2 M_\odot$ for stars with $0.5 < M_{ZAMS} \le 1$, and $\Gamma = -2.3$ for stars heavier than $1 M_\odot$ on zero-age main sequence (ZAMS). The distribution of mass ratios ($q = M_2/M_1$, where $M_{1/2}$ is the primary/secondary mass) was assumed to be uniform between $q_{min} = 0.08 M_\odot/M_1$ and 1. The initial distribution of orbital periods ($P$) and eccentricities ($e$) are $P(\log P) \sim (\log P)^{-0.55}$ and $P(e) \sim e^{-0.42}$ (Sana et al. 2012), which is the main difference in comparison to Wiktorowicz et al. (2017), but strong differences in the resulting binary populations are not expected (de Mink & Belczynski 2015; Klencki et al. 2018).

Every binary formed is evolved over 10 Gyr in isolation, i.e., no dynamical interactions with third bodies are taken into account, with special attention paid to interactions such as common envelope (CE; Ivanova et al. 2013) and MT (for details see Belczynski et al. 2008). To estimate the final compact object mass after a supernova, we use the “rapid” supernova formation mechanism (Belczynski et al. 2012; Fryer et al. 2012). For both NSs and BHs, we draw natal kicks from a Maxwellian distribution with $\sigma = 265$ km s$^{-1}$ (Hobbs et al. 2005), but scaled proportionally to the fraction of ejected mass that falls back onto the compact object. The kick velocity applied to a newly formed compact object ($v_{\text{kick,fin}}$) is obtained from $v_{\text{kick,fin}} = v_{\text{kick}}(1 - f_{\text{fb}})$, where $v_{\text{kick}}$ is the kick velocity that was drawn from a Maxwellian distribution with $\sigma = 265$ km s$^{-1}$, and $f_{\text{fb}}$ is the fraction of mass that was ejected in the SNa explosion that is accreted back onto the compact
object. We assumed that BHs forming via direct collapse obtain no natal kick.

We focus exclusively on sources undergoing Roche lobe overflow (RLOF) MT, during which, the X-ray luminosity is assumed to be (Shakura & Sunyaev 1973; Poutanen et al. 2007)

$$L_X = \begin{cases} L_{\text{Edd}}(1 + \ln \dot{m}_{\text{tr}}) & \dot{m}_{\text{tr}} > 1 \\ L_{\text{Edd}}\dot{m}_{\text{tr}} & \dot{m}_{\text{tr}} \le 1' \end{cases} \quad (3)$$

where $\dot{m}_{\text{tr}} = \dot{M}_{\text{tr}}/\dot{M}_{\text{Edd}}$ is the MT rate in Eddington units and $\dot{M}_{\text{tr}}$ is the MT rate. The apparent (spherical) luminosity is

$$L_{\text{app}} = L_X/b, \quad (4)$$

where the beaming factor is defined as $b \stackrel{\text{def}}{=} \Omega/4\pi = P_{\text{obs}}(\theta)$ and $\Omega$ is the combined solid angle of both beams. King (2009) showed that the observed relation of soft X-ray excess ($L_{\text{soft}}$) and disk temperature ($T_{\text{disk}}$) in ULXs, $L_{\text{soft}} \propto T_{\text{disk}}^{-4}$, implies

$$b \sim \frac{73}{\dot{m}_{\text{tr}}^2} \quad \text{for } \dot{m}_{\text{tr}} > 8.5, \quad (5)$$

whereas for $\dot{m} < 8.5$, emergent radiation is essentially unbeamed. In the following, the beaming is included as

$$b = \begin{cases} 1 & \dot{m}_{\text{tr}} \le 8.5, \\ \frac{73}{\dot{m}_{\text{tr}}^2} & 8.5 < \dot{m}_{\text{tr}}, \end{cases} \quad (6)$$

which provides monotonicity and continuity. These formulae have been successfully used to describe various classes of ULXs and hyperluminous X-ray sources (see, e.g., King & Lasota 2014, 2016; Lasota et al. 2015; King et al. 2017).

The prescription Equation (5) for the beaming factor is supported both by observations and theory (King et al. 2001; King 2009). Miller et al. (2013) analyzed a sample of ULXs finding them to broadly adhere to a $L_{\text{soft}} \propto T^4$ relation. While this may be due to freezing the absorption column—which, in at least one object, is observed to change, possibly with precession phase (Middleton et al. 2015b)—it may also indicate a changing fraction of energy lost via winds or very massive black holes very close to the Eddington limit. In Wiktorowicz et al. (2017), we assumed an ad hoc beaming saturation in order to avoid exceedingly small values of $b$ (and correspondingly large values of $L_{\text{app}}$). In the context of the present paper, this is not necessary because, as we have confirmed, extremely beamed sources are not only hard to observe, but also extremely rare and short-lived.

In Wiktorowicz et al. (2017), we have tested different beaming models applicable to population synthesis: a model with no beaming at all, a model with constant beaming ($b = 0.1$) for all super-Eddington sources, and a model based on the relation of the photosphere height to disk radius $H/R = 1.6/1 + r/\dot{m}_{\text{tr}}$ obtained with the GRMHD code KORAL (Lasota et al. 2016). The resulting numbers of ULXs, and ratio of BH to NS ULXs are highly similar (differing by a factor of $\lesssim 2-3$), except in the case of the least physical model with constant beaming, where the differences were more significant.

The results from other detailed GRMHD simulations, which are not applicable to population synthesis due to scarcity of tested configurations, frequently show results that diverge from the prescriptions of King et al. (2001; e.g., Sadowski et al. 2015). In these simulations, a highly collimated outflow never forms and the beam remains wide even for very high MT rates (Lasota et al. 2016). We see a similar problem when analyzing

the simulation of super-Eddington accretion on super-massive BHs (e.g Dai et al. 2018), which are relevant here due to the scale-free behavior of the most important relations (except for the radiation-to-gas-pressure ratio). Nevertheless, these codes give results comparable to those obtained with the KORAL code, so, as stated above, we predict a small effect on our results and especially on our conclusions. The only situation where the results concerning ULX are significantly different from detailed codes and prescription of King (2009) are the most luminous ULXs (so called extreme ULXs; Wiktorowicz et al. 2015), but as we have already said, these extreme ULX are unimportant for our general conclusions.

It should be stressed, however, that the King (2009) scenario uses the Shakura & Sunyaev (1973) model in which an accretion flow fed at a super-Eddington mass-transfer rate consists of an external, geometrically thin disk down to the spherization radius, below which the mass accretion is effectively only Eddington, the excess being ejected (advection of energy is negligible) in a quasi-spherical outflow whose collimating structure might be the cause of the anisotropic luminosity. In the case of neutron star accretors, this picture is strongly supported by the fact that the low-mass X-ray binary Cygnus X-2, has survived being fed $\sim 3 M_{\odot}$ from its (initially more massive) companion star at very super-Eddington rates ($\sim 10^{-5} M_{\odot} \text{yr}^{-1}$), but has evidently gained no more than $\sim \text{few} \times 0.2 M_{\odot}$ (King & Ritter 1999), evidently ejecting all the surplus. On the other hand, none of the GRMHD codes simulate the external thin disk but use a torus outside the inner $\sim 50$ gravitational radii, so their results are not directly comparable with the Shakura & Sunyaev (1973) model, since, e.g., the definitions of the spherization radii are not comparable. To confuse matters further, the luminosity formula $L = L_{\text{Edd}}(1 + \ln \dot{m}_{\text{tr}})$ is the same for the Shakura & Sunyaev (1973) and the advection dominated models (Poutanen et al. 2007; Lasota 2015), but in the first case $\dot{m}(R) = \dot{m}_{\text{tr}}$, whereas in the second $\dot{m}(R) \propto R$.

### 3. Results I: Volume-limited Surveys

The results presented in this section are for a uniform population of initial binaries (same metallicity of all stars: either solar $Z = Z_{\odot} = 0.02$, $10\% Z_{\odot} = 0.002$, or $1\% Z_{\odot} = 0.0002$) and assuming a simple model of star formation (SF): a constant star formation rate (SFR), with duration of 10 Gyr, or burst-like, with a duration of 100 Myr. Although both models form a total stellar mass of $6 \times 10^{10} M_{\odot}$, which corresponds to the stellar mass of the Milky Way (including bulge and disk; Licquia & Newman 2015), they cannot be directly compared to the Milky Way or any other complex stellar system that has various episodes of SF and nonuniform chemical composition. Nevertheless, estimates for more realistic systems may be obtained through use of our presented results.<fnref n="6" />

We focus on a volume-limited case in which we constrain the volume within which we observe ULXs. This then applies to volume-limited catalogs (e.g., Swartz et al. 2011; Walton et al. 2011), or galaxy-focused observations (e.g., Wolter et al. 2018). We note that the same volume of space may contain different amounts of stellar mass, (or, equivalently, the same stellar mass may occupy different volumes). In the following, we use the volume that contains the stellar mass of the Milky

<footnote n="6">⁶ Data files will be soon available at https://universeathome.pl/universe/ulx.php. In case of questions concerning the usage, please contact the first author.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Number of ULXs as a function of time since the beginning of star formation. The three upper panels present the results for constant star formation, whereas the three lower ones are for a star formation burst that lasted 100 Myr. The total stellar mass formed is the same for all plots and equals $6 \times 10^{10} M_{\odot}$, which is approximately the total stellar mass of the Milky Way (Licquia & Newman 2015). Three metallicities were considered (solar $Z = Z_{\odot} = 0.02$, 10% $Z_{\odot} = 0.002$, and 1% $Z_{\odot} = 0.0002$). On each panel, four lines present the total number of ULXs with BH, or NS accretors (gray solid and dashed line, respectively) and observed ULXs (red lines), i.e., ULXs whose beam intercepts the Earth (see Equation (2)). Although for BH ULXs, the difference between visible and total sample is small, for NS ULXs, the total sample is typically 5–15 times larger than the observed one.</figcaption></figure>

Way. For comparison with observations containing more/less stellar mass, the results should be scaled up/down accordingly, as the number of ULXs is directly proportional to stellar mass.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The toy model of the Milky Way constructed as a mix of 50% solar metallicity stars ($Z = 0.02$) and 50% sub-solar metallicity stars ($Z = 0.002$). Although the SFH of the Galaxy is much more complicated, here we present two simple cases: constant SF and burst SF (similarly to Figure 1) and discuss them in the text.</figcaption></figure>

Our results are presented in Figure 1. This is an updated and expanded version of Figure 2 from Wiktorowicz et al. (2017), where only the observed ULXs were presented. Here we also show the total number of ULXs (i.e., the sum of both observable and hidden sources; based on Equation (2)).

In order to demonstrate how our results can be applied for practical purposes, we have constructed a simple (“toy”) model of the Galaxy. We assumed that 50% of its stars have formed with solar metallicity ($Z = 0.02$) and 50% with a lower one ($Z = 0.002$). If we impose a constant SFR (Figure 2; upper panel), the resulting synthetic population (~100 observed ULXs) is strongly inconsistent with the actual one (no observed ULXs). However, measurements of the star formation history (SFH) in the Milky Way suggest that the SFR may have been significantly higher in the past (≥a few Gyr ago, e.g., Silva Aguirre et al. 2018). The SFH is, therefore, more burst-like. For such a case (Figure 2; bottom panel), we predict that only 2 Gyr after the burst, more than 80% of ULXs are hidden from our view, for which SS433 may be representative (Middleton et al. 2018).

Nevertheless, observations also suggest recent, small but significant SF in the Galaxy (~1 $M_{\odot}$ yr$^{-1}$ for ~1 Gyr, e.g., Maciel et al. 2012; The burst that might have happened ~2 Gyr, or earlier bursts, have little effect on the number of ULXs; see Figure 1; lower panels), which, according to the “toy” model, should produce ~20 (or ~3 if the SF occurs only in metal-rich environments; Figure 1, upper-most panel) observable ULXs at the current time. In contrast, in the Milky Way we observe only sources which become ULXs during their outbursts’ peaks when their emission goes above $10^{39}$ erg s$^{-1}$, such as the BeX binary Swift J0243.6+6124 ($L_{X,peak} \approx 5 \times 10^{39}$ erg s$^{-1}$; Tsygankov et al. 2017) or some low-mass-X-ray-binary transients (LMXBTs; Tetarenko et al. 2016). These “transient ULXs” are not a part of our results because BeX binaries do not belong to our sample, and the MT rates of the LMXBTs are well below the Eddington value.

A better agreement with observations is obtained when we do a simple scaling of our results to the observed stellar mass of galaxies within 14.5 Mpc ($M_{tot} \approx 3.5 \times 10^{12} M_{\odot}$) where 107 ULXs were found (Swartz et al. 2011). Assuming that in this volume the recent (~1 Gyr) SFR was small ~1 $M_{\odot}$ yr$^{-1}$ and chemical composition in this volume is similar to those used in
the “toy” model, we obtain a prediction of 175 ULXs, which is a less than a factor of two difference.

Possible sources of discordance include both theory and observation. It was shown in, e.g., Wiktorowicz et al. (2017) that the predicted number of ULXs may vary depending on the accretion model used. However, the ratio of NS to BH ULXs is only slightly affected, which agrees with our general conclusions. Some evolutionary phases that are important for the formation of ULXs (and X-ray binaries in general) like the CE are not well understood, and we do not have good models for them. Nevertheless, in Wiktorowicz et al. (2014), we showed that different CE models, which give significantly different CE outcomes, result in very similar predictions for X-ray binary (including ULX) populations, although populations of progenitors may differ significantly. Therefore, other evolutionary models may improve the fit to the simple model of the Milky Way presented above, but will not change our general conclusions. We also note that estimates of the observational parameters, which are necessary for population synthesis studies, like the total stellar mass and SFH, are not very precise. For example, the recent estimates of the Milky Way's stellar mass vary by a factor of ~2 (compare, e.g., Bovy & Rix 2013; Licquia & Newman 2015; Xiang et al. 2018).

### 3.1. Ratio of Observed to Total Sample

The initial NS mass in ULXs according to our simulations is typically around $1.3 M_{\odot}$. After formation, the NS's mass increases due to accretion, and in the ULX population is typically around $1.4 M_{\odot}$. If the ULX phase occurs early after ZAMS ($t_{age} \lesssim 500$ Myr), donors are typically main-sequence (MS) stars, and HG/RG (HG = Hertzsprung gap, RG = red giant) if it occurs later ($t_{age} \gtrsim 1$ Gyr). Typically, the donor mass is $M_{don} \lesssim 2 M_{\odot}$.

Within our simulation, we find that the majority of NS ULXs are beamed sources (Figures 3 and 4) that result in a low average ratio of the observed to the total number of NS ULXs (1/15–1/5, depending on the metallicity, Figure 1). The beaming results from the fact that an NS typically requires a strongly super-critical accretion rate to appear as a ULX because sub-critically accreting NSs will have Eddington-limited apparent luminosities ($L_X \lesssim L_{Edd,NS} \lesssim 5 \times 10^{38}$ erg s⁻¹), well below the empirically defined ULX luminosity of $L_{ULX} = 10^{39}$ erg s⁻¹. Indeed, a typical NS observed in a binary system ($M_{NS} \approx 1.4 M_{\odot}$) requires $\dot{m} > 10$ to reach $L_{app} > L_{ULX}$ (see Equation (4)), so the required beaming factor is always lower than $b \lesssim 0.7$. This means that the probability of observing a typical NS ULX is always lower than ~70%. Our results indicate that the average probability is, actually, only between ~7% and 20%.

Exceptions to the above do occur; we note that some NS ULXs may emit isotropically, especially when the metallicity of the environment is $Z \gtrsim 10\% Z_{\odot}$. In these cases, an NS may undergo a long phase of MT during which its mass increases to ~1.8 $M_{\odot}$ and the donor loses its hydrogen envelope. In such systems, the donor is typically a low-mass (~0.1 $M_{\odot}$) hybrid WD with a C–O–He rich core and a He rich envelope, in a very close orbit with the NS (orbital period $P < 1$ hr; see Belczynski & Taam 2004). Noting that the Eddington limit is higher for helium rich donors (Equation (1)), which implies that the isotropic emission of such a system may surpass $L_{ULX}$ when $\dot{m} \gtrsim 3$, which is significantly lower than required for an NS with a typical mass of $1.4 M_{\odot}$ accreting from a hydrogen-rich donor. These systems are good candidates for the brightest

ultra-compact X-ray binaries (see King 2011); however, their fraction among NS ULXs is only rarely expected to be higher than 4% (see Figure 3).

Unlike the condition for NSs, as the defining ULX luminosity ($L_{ULX} = 10^{39}$ erg s⁻¹) is the Eddington luminosity for a ~7 $M_{\odot}$ BH—a typical mass of a stellar-mass BH in the Milky Way (e.g., Özel et al. 2010)—a large fraction of BHs can become ULXs without the need for highly super-critical MT rates. When we take into account the lack of beaming in our models up to $\dot{m} < 8.5$, binaries with all stellar-mass BH may obtain ULX luminosities without additional amplification. According to the adopted beaming model (Equation (6)), these objects will emit isotropically, i.e., $P_{obs} = 1$. However, a fraction of the BH ULX population may be beamed due to highly super-critical MT rates and naturally results in high apparent luminosities ($\gtrsim 10^{40}$ erg s⁻¹; see Figure 1). Typically, the beamed fraction is $\lesssim 10\%$ of the total population (e.g., Figure 3, $Z = 0.02$), but it may exceed 50% in extreme cases (e.g., Figure 3: $Z = 0.002$, shortly after SF ceases) and can reach nearly 100% for extremely young populations (Figure 3: $t_{age} < 5$ Myr). BH ULXs most commonly appear in low-metallicity environments ($Z \ll Z_{\odot}$) and very young stellar populations ($t_{age} \lesssim 6$ Myr), and the majority of these systems will still have MT rates below $\dot{m} \approx 8.5$, because they are easier to obtain as a result of binary evolution. These “low-luminosity” BH ULXs (see Middleton et al. 2012; Middleton et al. 2013) outweigh the extreme (beamed) BH ULXs and, as a consequence, the ratio of the observed to the total sample of BH ULXs is typically only ~0.8. The average probability of observing a BH ULX from the Earth is therefore ~80%, far higher than the case for NS ULXs discussed above.

A current, important question in the field is “What is the ratio of observed NS to BH ULXs?”, and our simulations allow us to consider this under our model assumptions. For starburst systems (Figure 1, bottom panels) during SF, BH ULXs dominate the total and observed ULX population, while NS ULXs dominate these populations once SF has ceased. For prolonged/continuous SF (Figure 1, top panels), BH ULXs always dominate the total and observed ULX populations for low metallicities ($Z = 0.002$ and $Z = 0.0002$ models). For high metallicities ($Z = 0.2$; typical of the Milky Way disk), NS ULXs dominate the total population with an equal fraction of observed BH and NS ULXs. We note that, following the discussion above, the observed number of BH ULXs is always (in our simulated test cases) very similar to the total number of BH ULXs while the observed number of NS ULXs is always significantly below the total number of NS ULXs.

### 3.2. Beamed versus Isotropic as a Function of Stellar Population Age

Here we present a more detailed description of how the relative number of ULXs with different levels of beaming in the total population changes when the stellar population ages (depicted in Figure 3).

Young ULXs ($t_{age} \lesssim 10$ Myr) typically harbor BHs and massive (~10 $M_{\odot}$) MS donors that filled their Roche lobes (RL) due to nuclear evolution. However, the higher the metallicity, the stronger the mass loss in the stellar wind (e.g., Vink 2015); therefore, donors in a low-metallicity environment are usually more massive and, as a result, usually provide a higher thermal-timescale MT rate (as this is proportional to $M_{don}$). Additionally, low-metallicity stars do
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Fraction of BH/NS ULXs in the total population divided into four categories: sub-Eddington sources ($b = 1, \dot{m} \le 1$), isotropic sources ($b = 1, \dot{m} < 8.5$), mildly beamed sources ($1 < b < 0.1$), and highly beamed sources ($b \le 0.1$). Two star formation (SF) models are included: SF burst (duration 100 Myr; three upper plots) and constant SF (three lower plots). Three different metallicities are presented: $Z = Z_{\odot} = 0.02$, $Z = 10\% Z_{\odot} = 0.002$, and $Z = 1\% Z_{\odot} = 0.0002$. There are no NS ULXs in the sub-Eddington category.</figcaption></figure>

not expand as significantly during the MS compared to solar metallicity stars (e.g., Pols et al. 1998), so those that managed to fill their RL will tend to be more massive. As a result, for $Z < Z_{\odot}$, BHULXs are mostly significantly beamed ($b \le 0.1$), whereas for $Z = Z_{\odot}$ the emission is mainly isotropic. With time, the fraction of highly beamed BH ULXs quickly drops,

and at an age of ~10 Myr, the total ULX population is dominated by mildly beamed ($1 < b < 0.1$) and isotropic BH ULXs. A noticeable fraction of BH ULXs (up to 20% of the total population for sub-solar metallicities, i.e., $Z = 0.002$ and $Z = 0.0002$ models) harbor a relatively massive BH with $M_{BH} \gtrsim 10 M_{\odot}$, accreting at sub-Eddington rates. When the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The same as Figure 3, but for the *observed* population of ULXs.</figcaption></figure>

population is a few 10 Myr old, highly and mildly beamed NS ULXs emerge and for $Z = Z_{\odot}$ start to dominate the total ULX population at $t_{age} \approx 100$ Myr.

The continued evolution of the ULX population depends on the adopted SF history. In the case of constant SF (three lower plots on Figure 3). The fraction of highly and mildly beamed sources (mostly NS ULXs) grows steadily to becomes nearly constant after $t_{age} \approx 1$ Gyr. The fraction of beamed sources is

then ~80%, 40%, 50% for $Z = 0.02, Z = 0.002, Z = 0.0002$, respectively. The fractions of highly and mildly beamed sources are comparable. Highly beamed ULXs in an old stellar population ($t_{age} > 1$ Gyr) are mostly NS ULXs.

A different situation occurs in post-burst populations when the SF is extinguished. The number of BH ULXs drops quickly as no new BHs are produced and massive companions quickly evolve off the MS. NS ULXs, which are predominantly
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Histograms of beaming factor (b) in the total population of ULXs for three tested metallicities (marked in the upper left-hand corner) and three population ages: 100 Myr (at the end of a SF burst), 1 Gyr (900 Myr after the end of a SF burst), and after 10 Gyr of continuous SF. The rightmost peak (detached for the sake of clarity) represents isotropic emission (b = 1; for both super- and sub-Eddington sources). Other bins represent beamed emission (b < 1).</figcaption></figure>

beamed, quickly become dominant and in a few 100 Myr constitute nearly 100% of all ULXs. Therefore, post-burst ULX populations are predicted to be mostly beamed, NS systems.

Figure 5 presents the distributions of beaming factor in the total population of ULXs for different metallicities, accretor types, and three representative population ages. Two are for burst SF: 100 Myr after the end of SF burst and 900 Myr after the end, respectively. The third case represents the distribution for a population age of 10 Gyr and continuous SF. ULXs emitting isotropically were subtracted from the distribution and are presented in a detached bin. The left-most tail is formed by extremely beamed sources, which, although very luminous, contribute little to volume-limited populations due to a very small observation probability (Equation (2)).

### 3.2.1. The Most Highly Beamed Sources

ULXs characterized by the highest beaming of their radiation are the most luminous systems in the population. Although BH ULXs are predominantly isotropic emitters (except very young, $t_{age} \lesssim 10$ Myr, populations), some of them may exhibit strong beaming ($b < 0.1$). Such systems have typically Hertzsprung gap (HG) donors (in contrast to MS donors, which are the typical companions in BH ULXs in general; Wiktorowicz et al. 2015; Wiktorowicz et al. 2017) with masses between 5 and $10 M_{\odot}$. Typical BH masses of these ULXs are similar to those in isotropic ones ($\sim 6-10 M_{\odot}$) for all metallicities.

## 4. Results II: Flux-limited Surveys

Until now, we have considered the observational properties of ULXs that would be found in a volume-limited survey.

However, flux-limited surveys (e.g., the ROSAT all-sky survey and the forthcoming eROSITA all-sky survey) provide the broadest indication of the population of ULXs. Using our numerical results, we performed a comparison between simulated and analytical estimates (Middleton & King 2017) for the ratio of the number of NS ULXs to BH ULXs in idealized flux-limited observations, assuming that all ULXs directed toward the observer will be detectable and recognizable. Such observations will naturally be biased toward more strongly beamed, more luminous sources, which can be detected out to larger distances. As a result, the ratio of beamed to isotropic sources increases in comparison to volume-limited observations.

We assume that within large volumes ($D \gtrsim 10$ Mpc$^3$), the distribution of stellar mass is homogeneous and that light travel time does not influence the results significantly. Additionally, we use uniform (“toy”) universe models with the same metallicity and SFH at any place for an easier comparison of simulated and analytical estimates. After such a simplification, the limiting distance for detecting a source can be expressed as $D_{lim} = \sqrt{L_{X,app}/4\pi f_{lim}}$, where $L_{X,app}$ is the apparent luminosity of the source and $f_{lim}$ is the limiting observable flux. Consequently, the volume within which the source will be observable is $V \propto L_{X,app}^{3/2}$, where the scaling factor depends on $f_{lim}$, which we assume to be the same for all sources (i.e., the conditions of a flux-limited survey). If we additionally define the mean number density of stars as $n$, the probability of observing a particular source will be (See Middleton & King 2017)

$$P \propto f_{\text{SFH}} \cdot n \cdot b \cdot V \propto f_{\text{SFH}} \cdot n \cdot b \cdot L_{\text{X,app}}^{3/2} \cdot f_{\text{lim}}^{-3/2}, \quad (7)$$
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The upper panel shows the estimated ratio of NS ULXs to BH ULXs for our simulations ($E(N_{NS})/E(N_{BH})$; Equation (12)) and the analytic formula of Middleton & King (2017; $P_{NS}/P_{BH}$; Equation (6)). Solar metallicity ($Z = 0.02$) and constant SF through the last 10 Gyr were assumed. The middle and bottom panels present the number densities of NS/BH ULXs ($n(\text{NS}/\text{BH})$; middle panel), ratio of average masses of NS/BH accretors ($M_{\text{avr}}(\text{NS}/\text{BH})$; middle panel), and average masses of NS/BH separately (bottom panel) obtained from our simulations and used for the calculation of $P_{\text{NS}}/P_{\text{BH}}$.</figcaption></figure>

where $f_{\text{SFH}}$ is the probability that the particular system will presently be in the ULX phase. More precisely,

$$f_{\text{SFH}} = \frac{1}{\int_{0}^{10 \text{ Gyr}} \text{SFH}(t') dt'} \int_{t_{\text{age}}-t}^{t_{\text{age}}-t+dt} \text{SFH}(t') dt', \quad (8)$$

where $t_{\text{age}}$ is the age of the population since SF started, $t$ is the age of a given system (time since ZAMS) during the ULX phase, $dt$ is the length of the ULX phase, and SFH($t'$) is the SFH. For the simplified case, we defined SFH($t'$) as

$$\text{SFH}(t') = 6 \frac{M_{\odot}}{\text{yr}}, \quad (9)$$

for constant SF and as

$$\text{SFH}(t') = \begin{cases} 600 & t' \le 100 \text{ Myr} \frac{M_{\odot}}{\text{yr}}, \\ 0 & t' > 100 \text{ Myr} \frac{M_{\odot}}{\text{yr}}, \end{cases} \quad (10)$$

for burst SF. The estimated number of NS or BH ULXs ($N_{\text{NS/BH}}$) may then be calculated from

$$E(N_{\text{NS/BH}}) = \sum P_{\text{NS/BH}} \propto \sum_{\text{NS/BH ULXs}} f_{\text{SFH}} \cdot b \cdot L_{\text{app}}^{\frac{3}{2}}, \quad (11)$$

where $P_{\text{NS/BH}}$ is the probability of observation of a particular NS or BH ULX (Equation (7)), and the summation is performed over entire ULX lifetime for all NS and BH ULXs. The ratio of NS ULXs to BH ULXs is then

$$\frac{E(N_{\text{NS}})}{E(N_{\text{BH}})} = \frac{\sum_{\text{NS ULXs}} f_{\text{SFH}} \cdot b \cdot L_{\text{app}}^{\frac{3}{2}}}{\sum_{\text{BH ULXs}} f_{\text{SFH}} \cdot b \cdot L_{\text{app}}^{\frac{3}{2}}}. \quad (12)$$

In Figure 6, we show the comparison of our simulated estimate of the relative number of NS ULXs to BH ULXs

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as Figure 6, upper panel, but for different metallicities. For lower metallicities, there are relatively far fewer NS ULXs (both for our results and using the Middleton & King (2017) prescription) than for solar metallicity ($Z = 0.02$).</figcaption></figure>

(Equation (12)), to the simplified analytical formula of Middleton & King (2017, Equation (6)). The number densities and mean masses of accretors enter into their formula, and we calculate these directly from our simulations (which are also provided in Figure 6). Clearly the results are consistent at $t_{\text{age}} < 1$ Gyr and only differ across the entire simulation by a factor $\lesssim 2$; this is a clear validation of the simplified approach by Middleton & King (2017). We note that when the various contributions to the overall population are considered, our results diverge to a greater extent at low metallicities, and at late times, they can differ by up to an order of magnitude for $Z = 0.0002$ (see Figure 7). There are two reasons for this discrepancy: first, Middleton & King (2017) assumed that $\dot{m}_{\text{tr}}$ ($\dot{m}_0$ in their work) is always $\ge 1$ (see their Equation (1)), whereas we consider BH ULXs that can be classified as ULXs with $\dot{m}_{\text{tr}} < 1$ and in the population are mostly unbeamed. Second, in obtaining their Equation (4), it was assumed that MT rates for NS ULXs and BH ULXs are similar whereas this is not always true because BH ULXs can achieve much higher stable MT rates (e.g., from massive stars, $M_{\text{donor}} > 10 M_{\odot}$), which would otherwise lead to a dynamical instability for NS accretors. These assumptions contribute to the discrepancy the most where the metallicity of the environment is the lowest, as masses of BHs are then higher on average, so there are more ULXs emitting isotropically.

## 5. Discussion

Our simulations should allow us to answer the fundamental question, “Is beaming necessary to explain the observed population of ULXs?” Unfortunately, without detailed information on the environment in which observed ULXs reside, such as SFH and metallicity distribution, we are unable to make reliable predictions. However, the overall results from our simulations regarding the populations of NS and BH ULXs and
the role of beaming should be robust to different evolutionary models (Wiktorowicz et al. 2017) and these are important for understanding the observed population in both volume and flux-limited surveys.

There are two key results from our simulations: the fraction of observed to total sample of ULXs and the ratio of NS versus BH ULXs. In general, for star-forming regions, the fraction of observed to the total sample of ULXs is ~0.8 (independently of metallicity) due to a high abundance of BH ULXs that typically emit isotropically, except for the very early ages ($t_{age} < 10$ Myr) when most of the BH ULX are significantly beamed and the fraction of observed to total sample is smaller. Conversely, for a solar metallicity environment after a long SF episode (>$\sim$100 Myr of continuous SF), the observed population is dominated by NS ULXs. For old stellar populations where the SF ceased ~1 Gyr ago, the ratio of observed to total population of ULXs is typically ~0.2 because it consists nearly exclusively of NS ULXs, which are nearly always beamed.

As we show in Figure 1, the relative fraction of NS to BH ULXs changes as a function of both metallicity, SF model (continuous or burst), and age since SF commenced. In the case of continuous SF, BH ULXs dominate the observed population of ULXs (the lowest ratio to NS ULXs is 1:1 for solar metallicities at late times). For burst SF, young ULX populations are dominated by BH ULXs, but this changes as the population ages and, post star formation, NS ULXs dominate both the observed and total population of ULXs.

Our adopted geometrical beaming model predicts extremely strong beaming for some ULXs (Figure 5). The corresponding apparent luminosities are well above those observed for extreme ULXs ($L_{X,max} \approx 1 \times 10^{42}$ erg s⁻¹). To avoid such a situation, in the previous work (Wiktorowicz et al. 2017), we applied a saturation threshold for beaming at $b_{lim} = 3.2 \times 10^{-3}$ (ṁ$_{lim} \approx 150$), which capped the luminosities to $L_X \lesssim 10^{42}$ erg s⁻¹. In this paper, we found that adopting the saturation does not change the results and conclusions significantly.

King et al. (2001) suggested that the thermal-timescale MT is the best process to fuel a ULX. Indeed, part of the ULXs in our results, including the most luminous ones, experience a thermal-timescale MT. In these sources, the donor is not in thermal equilibrium and its characteristics (especially the mass–radius relation) may differ strongly from those calculated for thermally stable stars, which are used in population synthesis codes. Consequently, the predicted observational properties of ULXs (e.g., the duration of the ULX phase, or peak luminosity) may differ in contrast to the results of detailed codes (e.g., Belczynski et al. 2008, find a difference by a factor of ~2 in the duration of the RLOF phase of the $1.3 M_\odot$ NS and $1.6 M_\odot$ MS donor). However, in our results, the observed and total populations are dominated by ULXs undergoing nuclear timescale MT (the exceptions are very old populations; see Figure 8). This is the direct consequence of a much longer duration of the MT phases. In spite of the problems with calculating the realistic MT rates and durations of thermal-timescale mass transfer, we show that it can affect our results only for old stellar populations, whereas ULXs are mostly found in star-forming environments (e.g., Fabbiano et al. 2003; Gao et al. 2003).

Specifically, we found that typically a BH ULX (except the most luminous ones with $L_X \gtrsim 10^{41}$ erg s⁻¹), undergo a nuclear timescale MT. Therefore, ULX populations dominated by BH

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Fraction of ULXs (with both BH and NS accretors) undergoing thermal-timescale mass transfer in the total (tot; continuous lines) and observed (obs; dashed lines) population. Two star formation (SF) models are included: constant SF (upper panel) and SF burst (duration 100 Myr; lower panel). Three different metallicities are presented: $Z = Z_\odot = 0.02$, $Z = 10\% Z_\odot = 0.002$, and $Z = 1\% Z_\odot = 0.0002$.</figcaption></figure>

accretors are, consequently, dominated by nuclear timescale MT (typically, the fraction is $\gtrsim$90% for both total and observed populations). As we predicted, BH-dominated ULX populations are typical for SF bursts and constant SF environments (except the total population in prolonged constant SF in solar metallicity environment; see Figure 1), where most of the ULXs are found (e.g., Swartz et al. 2004, 2009). Also, the majority of ULXs are expected to undergo an MT on a nuclear timescale.

NS ULXs may have a much higher fraction of systems undergoing a thermal-timescale MT that depends on SFH (see Figure 8). For the solar metallicity environment, the fraction of thermal-timescale MT among ULXs reaches ~40% (~5% for observed population), mostly due to a high fraction of NS accretors. Also, in older environments (a few 100 Myr after the cessation of SF), the population is dominated by NS ULX; thus, the fraction of thermal-timescale MT is higher (up to ~50% for total population of up to ~20% for observed population). We note that the fraction reaches 100% for very old stellar populations ($t_{age} \gtrsim 9$ Gyr) in low-metallicity environments ($Z = 0.002$ and $Z = 0.0002$), but it may be a result of very limited statistics.

## 6. Summary

In this paper, we have analyzed how geometrical beaming, i.e., anisotropic emission of radiation, affects the observed sample of ULXs when compared to the total sample of these objects, some part of which is hidden from our view. Our simulated results are based on the previous analysis of ULX populations in different environments published in Wiktorowicz et al. (2017) and may be seen as an extension of the previous work. The utilized beaming model (beaming
factor $b \propto \dot{m}^{-2}$) is based on theoretical and observational grounds (King 2009).

We show that ULXs harboring BH accretors are typically emitting isotropically ($b = 1$) and undergo a nuclear timescale MT, whereas those with NS accretors are predominantly beamed (typically $b = 0.07-0.2$) and in most cases MT occurs on a thermal timescale. Our analysis shows that the beaming is dependent on different stellar environments; very young (burst) populations ($t_{age} < 10$ Myr), dominated by BH ULXs, are significantly beamed while BH ULXs in older stellar populations are usually isotropic emitters. However, the majority of NS ULXs are always beamed, irrespective of stellar environment. In terms of the relative ratio of species, we find that the ratio of NS ULXs to BH ULXs is higher in the total sample than in the observed sample. In the case of continuous SF, BH ULXs typically outnumber the NS ULXs in the observed sample. While BH ULXs also outnumber the NS ULXs in the observed sample for burst SF at early times, post SF, NS ULXs tend to dominate the observed population instead. In the case of the latter, the observed NS ULXs represent only 20% of the total NS ULX population, and many are expected to be obscured from view (in the absence of precession, which may act to bring some into view—see Dauser et al. 2017; Middleton et al. 2018).

Finally, we found that the ratio of the number of NS ULXs to BH ULXs in idealized flux-limited observations is consistent within a factor $\lesssim 2$ of that found by Middleton & King (2017), with divergence at late times and lower metallicities where, in the case of the latter, large MT rates can lead to instabilities for accreting NS systems, and very massive BHs are more common.

We are grateful to the anonymous referee who helped us to improve the paper. J.P.L. and M.M. thank Andrew King for many inspiring discussions. We are thankful to thousands of volunteers, who took part in the *Universe@Home* project<fnref n="7" /> and provided their computers for this research. G.W. is partly supported by the Presidents International Fellowship Initiative (PIFI) of the Chinese Academy of Sciences under grant No. 2018PM0017 and by the Strategic Priority Research Program of the Chinese Academy of Science Multi-waveband Gravitational Wave universe (grant No. XDB23040000). J.P.L. and K.B. acknowledge support by the National Science Centre, Poland grant 2015/19/B/ST9/01099 and JPL by a grant from the French Space Agency CNES. M.M. thanks STFC for support via an Ernest Rutherford Fellowship.

## ORCID iDs

Matthew Middleton https://orcid.org/0000-0002-8183-2970

## References

* Abarca, D., Kluźniak, W., & Sadowski, A. 2018, MNRAS, 479, 3936

* Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, ApJL, 818, L22

* Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Natur, 514, 202

* Basko, M. M., & Sunyaev, R. A. 1976, MNRAS, 175, 395

* Begelman, M. C., King, A. R., & Pringle, J. E. 2006, MNRAS, 370, 399

* Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407

* Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223

* Belczynski, K., & Taam, R. E. 2004, ApJ, 603, 690

* Belczynski, K., Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012, ApJ, 757, 91

* Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115

* Brightman, M., Harrison, F. A., Fürst, F., et al. 2018, NatAs, 2, 312

* Carpano, S., Haberl, F., Maitra, C., & Vasilopoulos, G. 2018, MNRAS, 476, L45

* Clark, J. S., Ritchie, B. W., Najarro, F., Langer, N., & Negueruela, I. 2014, A&A, 565, A90

* Colbert, E. J. M., & Mushotzky, R. F. 1999, ApJ, 519, 89

* Dai, L., McKinney, J. C., Roth, N., Ramirez-Ruiz, E., & Miller, M. C. 2018, ApJL, 859, L20

* Dauser, T., Middleton, M., & Wilms, J. 2017, MNRAS, 466, 2236

* de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58

* Fabbiano, G., King, A. R., Zezas, A., et al. 2003, ApJ, 591, 843

* Fabrika, S. 2004, ASPRv, 12, 1

* Finke, J. D., & Razzaque, S. 2017, MNRAS, 472, 3683

* Fragos, T., Linden, T., Kalogera, V., & Sklias, P. 2015, ApJL, 802, L5

* Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91

* Fürst, F., Walton, D. J., Harrison, F. A., et al. 2016, ApJL, 831, L14

* Gao, Y., Wang, Q. D., Appleton, P. N., & Lucas, R. A. 2003, ApJL, 596, L171

* Heida, M., Jonker, P. G., Torres, M. A. P., et al. 2014, MNRAS, 442, 1054

* Heida, M., Jonker, P. G., Torres, M. A. P., et al. 2016, MNRAS, 459, 771

* Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974

* Israel, G. L., Belfiore, A., Stella, L., et al. 2017, Sci, 355, 817

* Ivanova, N., Justham, S., Chen, X., et al. 2013, A&ARv, 21, 59

* Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303

* Kaaret, P., Ward, M. J., & Zezas, A. 2004, MNRAS, 351, L83

* King, A. 2011, ApJL, 732, L28

* King, A., & Lasota, J.-P. 2014, MNRAS, 444, L30

* King, A., & Lasota, J.-P. 2016, MNRAS, 458, L10

* King, A., Lasota, J.-P., & Kluźniak, W. 2017, MNRAS, 468, L59

* King, A. R. 2009, MNRAS, 393, L41

* King, A. R., Davies, M. B., Ward, M. J., Fabbiano, G., & Elvis, M. 2001, ApJL, 552, L109

* King, A. R., & Lasota, J.-P. 2019, MNRAS, 485, 3588

* King, A. R., & Ritter, H. 1999, MNRAS, 309, 253

* Klencki, J., Moe, M., Gladysz, W., et al. 2018, A&A, 619, A77

* Koliopanos, F., Vasilopoulos, G., Buchner, J., Maitra, C., & Haberl, F. 2019, A&A, 621, A118

* Koliopanos, F., Vasilopoulos, G., Godet, O., et al. 2017, A&A, 608, A47

* Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076

* Lasota, J.-P. 2015, arXiv:1505.02172

* Lasota, J.-P., King, A. R., & Dubus, G. 2015, ApJL, 801, L4

* Lasota, J.-P., Vieira, R. S. S., Sadowski, A., Narayan, R., & Abramowicz, M. A. 2016, A&A, 587, A13

* Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96

* Linden, T., Kalogera, V., Sepinsky, J. F., et al. 2010, ApJ, 725, 1984

* Liu, J.-F., Bregman, J. N., & Seitzer, P. 2004, ApJ, 602, 249

* Maciel, W. J., Rocha-Pinto, H. J., & Costa, R. D. D. 2012, in IAU Symp. 284, The Spectral Energy Distribution of Galaxies—SED 2011, ed. R. J. Tuffs & C. C. Popescu (Cambridge: Cambridge Univ. Press), 379

* Madhusudhan, N., Justham, S., Nelson, L., et al. 2006, ApJ, 640, 918

* Madhusudhan, N., Rappaport, S., Podsiadlowski, P., & Nelson, L. 2008, ApJ, 688, 1235

* Marchant, P., Langer, N., Podsiadlowski, P., et al. 2017, A&A, 604, A55

* Medvedev, A., & Fabrika, S. 2010, MNRAS, 402, 479

* Middleton, M. J., Heil, L., Pintore, F., Walton, D. J., & Roberts, T. P. 2015a, MNRAS, 447, 3243

* Middleton, M. J., & King, A. 2016, MNRAS, 462, L71

* Middleton, M. J., & King, A. 2017, MNRAS, 470, L69

* Middleton, M. J., Miller-Jones, J. C. A., Markoff, S., et al. 2013, Natur, 493, 187

* Middleton, M. J., Sutton, A. D., Roberts, T. P., Jackson, F. E., & Done, C. 2012, MNRAS, 420, 2969

* Middleton, M. J., Walton, D. J., Alston, W., et al. 2018, arXiv:1810.10518

* Middleton, M. J., Walton, D. J., Fabian, A., et al. 2015b, MNRAS, 454, 3134

* Miller, J. M., Walton, D. J., King, A. L., et al. 2013, ApJ, 776, L36

* Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Ingram, A. 2017, MNRAS, 467, 1202

* Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015a, MNRAS, 454, 2539

* Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015b, MNRAS, 447, 1847

* Ohsuga, K., & Mineshige, S. 2011, ApJ, 736, 2

* Ohsuga, K., Mineshige, S., Mori, M., & Kato, Y. 2009, PASJ, 61, L7

* Olausen, S. A., & Kaspi, V. M. 2014, ApJS, 212, 6

* Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918

<footnote n="7">⁷ https://universeathome.pl/universe/</footnote>
Parfrey, K., & Tchekhovskoy, A. 2017, ApJL, 851, L34

Pintore, F., Zampieri, L., Stella, L., et al. 2017, ApJ, 836, 113

Pols, O. R., Schröder, K.-P., Hurley, J. R., Tout, C. A., & Eggleton, P. P. 1998, MNRAS, 298, 525

Popov, S. B. 2016, A&AT, 29, 183

Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187

Rappaport, S. A., Podsiadlowski, P., & Pfahl, E. 2005, MNRAS, 356, 401

Sądowski, A., Narayan, R., McKinney, J. C., & Tchekhovskoy, A. 2014, MNRAS, 439, 503

Sądowski, A., Narayan, R., Tchekhovskoy, A., et al. 2015, MNRAS, 447, 49

Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Sci, 337, 444

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Shao, Y., & Li, X.-D. 2015, ApJ, 802, 131

Silva Aguirre, V., Bojsen-Hansen, M., Slumstrup, D., et al. 2018, MNRAS, 475, 5487

Soria, R., Long, K. S., Blair, W. P., et al. 2014, Sci, 343, 1330

Swartz, D. A., Ghosh, K. K., Tennant, A. F., & Wu, K. 2004, ApJS, 154, 519

Swartz, D. A., Soria, R., Tennant, A. F., & Yukita, M. 2011, ApJ, 741, 49

Swartz, D. A., Tennant, A. F., & Soria, R. 2009, ApJ, 703, 159

Tetarenko, B. E., Sivakoff, G. R., Heinke, C. O., & Gladstone, J. C. 2016, ApJS, 222, 15

Tsygankov, S. S., Doroshenko, V., Lutovinov, A. A., Mushtukov, A. A., & Poutanen, J. 2017, A&A, 605, A39

Vink, J. S. 2015, in Astrophysics and Space Science Library, Vol. 412, Very Massive Stars in the Local Universe (Cham: Springer International), 77

Walton, D. J., Fürst, F., Heida, M., et al. 2018, ApJ, 856, 128

Walton, D. J., Roberts, T. P., Mateos, S., & Heard, V. 2011, MNRAS, 416, 1844

Wiktorowicz, G., Belczynski, K., & Maccarone, T. 2014, Binary Systems, their Evolution and Environments, 37

Wiktorowicz, G., Sobolewska, M., Lasota, J.-P., & Belczynski, K. 2017, ApJ, 846, 17

Wiktorowicz, G., Sobolewska, M., Sądowski, A., & Belczynski, K. 2015, ApJ, 810, 20

Wolter, A., Fruscione, A., & Mapelli, M. 2018, ApJ, 863, 43

Xiang, M., Shi, J., Liu, X., et al. 2018, ApJS, 237, 33