# The ionizing and heating power of ultraluminous X-ray sources under the geometrical beaming model

K. Kovlakas, T. Fragos, D. Schaerer, A. Mesinger

► To cite this version:

K. Kovlakas, T. Fragos, D. Schaerer, A. Mesinger. The ionizing and heating power of ultraluminous X-ray sources under the geometrical beaming model. Astronomy & Astrophysics - A&A, 2022, 665, pp.A28. 10.1051/0004-6361/202244252. [hal-03758781](https://hal.science/hal-03758781v1)

HAL Id: hal-03758781

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

Submitted on 9 Mar 2023

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 CC BY 4.0 - Attribution - International License
# The ionizing and heating power of ultraluminous X-ray sources under the geometrical beaming model*

K. Kovlakas<fnref n="1" />, D. Fragos<fnref n="1" />, D. Schaerer<fnref n="1" />,<fnref n="2" />, and A. Mesinger<fnref n="3" />

Received 11 June 2022 / Accepted 15 July 2022

## ABSTRACT

While there is now a consensus that X-ray binaries (XRBs) are the dominant X-ray sources in the early Universe and play a significant role during the epoch of heating of the intergalactic medium (IGM), recent studies report contradicting results regarding their contribution in the nebular emission of local Universe galaxies. Ultraluminous X-ray sources (ULXs), which dominate the X-ray budget of normal galaxies, may be important interstellar-medium (ISM) ionizing sources. However, their output in the extreme-ultraviolet (EUV) and soft X-ray part of the spectrum remains observationally unconstrained. In this paper, we predict the ionizing and heating power from ULX populations under the geometrical beaming scenario, and three models describing the emission from super-critical accretion disks. We find that our theoretical spectra for ULX populations cannot (can) explain the He II (Ne V) emission observed in some galaxies, with their contribution being less (more) important than the underlying stellar population. Stochastic fluctuations in the number of ULXs may allow for equal contributions in the He II emission, in a fraction of galaxies. We provide average spectra of ULX populations as an input to local, and early-Universe studies. We find that the soft X-ray emission arising from super-critical accretion is significant for the heating of the IGM, and consistent with recent constraints from the 21-cm cosmic signal. Based on the dependence on the adopted compact-object (CO) mass and accretion model, we encourage efforts in modeling ULX spectra via simulations, and their combination with detailed binary population synthesis models.

**Key words.** accretion, accretion disks – X-rays: binaries – X-rays: ISM – X-rays: galaxies – ultraviolet: galaxies

## 1. Introduction

Low-metallicity galaxies have been found to exhibit high nebular He II and Ne V emission (e.g., Garnett et al. 1991; Izotov et al. 2012). The origin of this emission remains unknown, with multiple explanations in the recent literature: radiative shocks (e.g., Plat et al. 2019), binary stars (e.g., Götberg et al. 2017), active galactic nuclei in dwarf galaxies (e.g., Mezcua et al. 2018), X-ray binaries (XRBs; e.g., Power et al. 2013), cluster winds (e.g., Oskinova & Schaerer 2022), among others (see discussion in Olivier et al. 2021). Understanding ionizing sources is crucial for the study of the interstellar medium (ISM; e.g., Nanayakkara et al. 2019), and the heating of the intergalactic medium (IGM; e.g., Mesinger et al. 2013; Madau & Fragos 2017).

Observations show that XRB populations and their integrated luminosity are anticorrelated with the metallicity of the host galaxy (e.g., Douna et al. 2015); this is also supported by theoretical work (e.g., Fragos et al. 2013a,b). This has motivated the study of the contribution of XRBs in the ionization of the ISM in low-metallicity galaxies, with contradicting results (e.g., Schaerer et al. 2019; Saxena et al. 2020; Senchyna et al. 2020; Umeda et al. 2022). In the case of the extremely low-metallicity galaxy IZw 18, exhibiting strong He II emission, the X-ray luminosity is dominated by an ultraluminous X-ray source (ULX; see review of Kaaret et al. 2017), presenting an interesting case

for the study of individual XRBs as ionizing sources. However, recent observational work is inconclusive: the source position and spectrum is not consistent with the He II emission (e.g., Kehrig et al. 2021), unless a geometrical beaming effect is considered (e.g., Rickards Vaught et al. 2021).

One of the most important difficulties in understanding the contribution of X-ray sources, which are typically observed at energy bands 0.2–10 keV, is that their spectra are poorly constrained at the extreme-ultraviolet (EUV) regime (He II: ~54 eV, Ne V: ~97.1 eV; see Lehmer et al. 2022). However, ULXs are strong candidates for four reasons: (i) ULXs exhibit softer spectra than typical XRBs and can potentially produce a high rate of He-ionizing photons (e.g., Simmonds et al. 2021); (ii) they have been shown to dominate the X-ray emission from normal galaxies due to the shallow slope of the high-mass XRB (HMXB) luminosity function (XLF; e.g., Lehmer et al. 2019); (iii) they are abundant in low-metallicity environments (e.g., Mapelli et al. 2010); and (iv) there is evidence of geometrical beaming (e.g., King 2009), which indicates a larger underlying population of ULXs, and hence a stronger contribution in the ionization of the ISM than what is expected by the observed population (e.g., Rickards Vaught et al. 2021).

An interesting case of a powerful ionizing source in the Galaxy is SS 433 (Stephenson & Sanduleak 1977). Its high optical and UV emission (e.g., Waisberg et al. 2019) is believed to be powered by a super-critical accretion disk, despite its low observed X-ray luminosity (e.g., Middleton et al. 2021). This can be explained through models of super-critical disks, which predict that close to the compact object (CO), the disk becomes

<footnote n="1">¹ Departement d'Astronomie, Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland
e-mail: konstantinos.kovlakas@unige.ch</footnote>

<footnote n="2">² CNRS, IRAP, 14 avenue E. Berlin, 31400 Toulouse, France</footnote>

<footnote n="3">³ Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy</footnote>

<footnote>* Full Table 1 is only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/665/A28</footnote>
geometrically thick, as strong optically-thick winds flow from its surface (e.g., Shakura & Sunyaev 1973, henceforth SS73). The highly anisotropic radiation pattern results in face-on observers detecting collimated emission from the inner regions of the disk (beaming), and edge-on observers seeing “softer” emission escaping from the wind’s photosphere extending to large radii with respect to the inner disk radius (e.g., Poutanen et al. 2007, henceforth P07). Therefore, the common interpretation is that SS 433 belongs to the general population of ULXs with observed luminosities $L_{\text{obs}} \ge 10^{39} \text{ erg s}^{-1}$, but, being viewed at a high inclination, its strong X-ray emission is invisible to us (e.g., Begelman et al. 2006). Consequently, a unified picture arises where the main difference between ULXs and soft ULXs is the viewing angle (e.g., Sutton et al. 2013; Pinto et al. 2017).

In this paper, we combine the properties of observed ULX populations and theoretical models describing their intrinsic spectra. We constrain the ionization power of ULXs under the super-Eddington accretion and beaming models.

## 2. Methodology

The contribution of ULX populations in the ionization of the ISM and the heating of the IGM is hampered by three observational challenges. Firstly, the hard EUV and soft X-ray part of the spectrum cannot be detected directly due to absorption. Secondly, under the geometrical-beaming scenario, extragalactic edge-on ULXs are too faint to be detected by X-ray telescopes, and therefore, we may underestimate the ULX content of galaxies. Thirdly, the distributions of the physical properties of ULXs (e.g., accretor and donor masses, structure of the accretion disk, etc.) remain unknown. To overcome these challenges in our analysis, we adopted models for super-critically accreting sources that provide the luminosity, geometrical beaming, and spectral energy distribution (SED), given the CO mass and mass-transfer rate of the systems. Below, we show that by considering different CO masses, the beaming factor and the spectrum only depend on the luminosity, and therefore, the underlying ULX population and its integrated spectra can be anchored on empirical constraints of the HMXB XLF.

### 2.1. Luminosity distribution and beaming in ULXs

Super-critical disks are, by definition, encountered in systems with luminosities exceeding the Eddington limit:

$$L_{\text{Edd}} = 1.26 \times 10^{38} m \text{ erg s}^{-1}, \qquad (1)$$

where $m = M/M_{\odot}$ is the mass of the accreting object in solar units (e.g., Abramowicz et al. 1980). Based on the typical range of black-hole (BH) masses in XRBs ($\le 15 M_{\odot}$; Remillard & McClintock 2006), sources with luminosities exceeding $L_{\text{lim}} = 2 \times 10^{39} \text{ erg s}^{-1}$ are likely to be super-Eddington. Although super-critically accreting COs with near-solar masses may exhibit lower luminosities than the aforementioned limit, when we anchored our results on XLFs, we ignored the part below $L_{\text{lim}}$ to avoid mixing under- and super-Eddington sources. In SS73 it is shown that when the mass-transfer rate $\dot{M}$ exceeds the rate corresponding to the Eddington limit, $\dot{M}_{\text{Edd}}$, the bolometric luminosity is

$$L_{\text{bol}} = L_{\text{Edd}} (1 + \ln \dot{m}), \qquad (2)$$

where $\dot{m} = \dot{M}/\dot{M}_{\text{Edd}}$. This is the result of a fraction of accretion power being spent for driving strong outflows, which keep the accretion on the CO Eddington-limited. Even at extreme mass-transfer rates (e.g., $\dot{m} \sim 1000$), the bolometric luminosity cannot reach $10^{41} \text{ erg s}^{-1}$ without invoking COs of $M > 100 M_{\odot}$.

In the beaming scenario of King et al. (2001), the most extreme ULXs are explained by the fact that their observed X-ray luminosities ($L_{\text{obs}}$) are the isotropic equivalent of moderately super-Eddington, beamed sources viewed face-on:

$$L_{\text{obs}} = b^{-1} L_{\text{bol}}, \qquad (3)$$

where $b$ is the beaming factor. Furthermore, King (2009) hinted at a dependence of the beaming factor on the mass-transfer rate,

$$b = \begin{cases} 1 & \dot{m} \le 8.5 \\ (8.5/\dot{m})^2 & \dot{m} > 8.5' \end{cases}, \qquad (4)$$

which explains the observed anticorrelation between the luminosity and temperature of the soft X-ray emission in ULXs (e.g., Kajava & Poutanen 2009; Feng et al. 2016).

We used the $L_{\text{obs}} > 2 \times 10^{39} \text{ erg s}^{-1}$ part of the HMXB XLF from Lehmer et al. (2019) to model the distribution of ULX luminosities:

$$\frac{dN}{dL_{38}} = \text{SFR} \times K_{\text{HMXB}} \times \begin{cases} L_{38}^{-\gamma} & 20 < L_{38} < L_c \\ 0 & \text{elsewhere} \end{cases}, \qquad (5)$$

where $L_{38} = L_{\text{obs}}/10^{38} \text{ erg s}^{-1}$, SFR is the star-formation rate of the parent stellar population, and the normalization $K_{\text{HMXB}}$, power-law index $\gamma$, and cut-off luminosity $L_c$ are fitted parameters. Equations (1)–(5) interconnect the observable $L_{\text{obs}}$ with the mass-transfer rate $\dot{m}$ via the XLF, assuming a distribution of masses $m$. However, ULXs are rare sources (typically one per galaxy; e.g., Kovlakas et al. 2020), and their distances are prohibiting when constraining their parameters, such as the mass of the accretor. For this reason, we considered three different accretor masses: 1.4 (corresponding to neutron stars<fnref n="1" />), 8, and 20 $M_{\odot}$. Using Eqs. (2)–(4) and the fact that by definition at sub-Eddington luminosities $L_{\text{bol}} = \dot{m}L_{\text{Edd}}$, we inferred $\dot{m}$ from the ratio $f = L_{\text{obs}}/L_{\text{Edd}}$:

$$\dot{m} = \begin{cases} f & f \le 1 \\ e^{f-1} & 1 < f \le 3.14, \\ 8.5 \left( \frac{f}{3.14} \right)^{\frac{4}{9}} & f > 3.14 \end{cases}, \qquad (6)$$

where the case for $f > 1 + \ln 8.5 \approx 3.14$ is an approximation of the inverse of $f = (\dot{m}/8.5)^2 (1 + \ln \dot{m})$ with an accuracy of <2% for $f < 800$ (a sufficiently high value corresponding to $10^{41} \text{ erg s}^{-1}$ for $m = 1$).

### 2.2. The blackbody component of the ULX spectrum

The spectral properties of super-Eddington sources are highly dependent on the geometry of the disk, as well as the interplay between the outflowing gas and the radiation. Specifically, the wind’s photosphere is expected to produce a soft component in the spectrum of ULXs (e.g., King & Pounds 2003).

According to the SS73 model, the local mass-transfer rate increases as we go towards the inner regions of the accretion disk. However, at super-Eddington rates, there is a radius, the spherization radius ($R_{\text{sp}}$), at which the local Eddington limit is reached, initiating outflows. As a result, the mass-transfer rate decreases as we approach the CO further. At the $R_{\text{sp}}$ the accretion disk is optically and geometrically thick, forming a

<footnote n="1">¹ The NS mass was used as a lower limit; we did not account for the presence of the magnetic field and alternative accretion modes (e.g., Brice et al. 2021).</footnote>
nearly-spherical blackbody source. The R<sub>sp</sub> was computed by equalizing the accretion and Eddington luminosities (cf., Frank et al. 2002):

$$R_{\text{sp}} = \frac{GM\dot{M}}{\eta \dot{M}_{\text{Edd}} c^2} = \frac{GM_{\odot}}{\eta c^2} m\dot{m} = 1.5 \times 10^5 m\dot{m} \text{ cm}, \quad (7)$$

where η is the accretion efficiency, for which we adopted the value of 0.1. Using the Stefan-Boltzmann law, we calculated the effective temperature at the R<sub>sp</sub>:

$$T_{\text{eff}} = \left( \frac{L_{\text{Edd}}}{4\pi R_{\text{sp}}^2 \sigma} \right)^{\frac{1}{4}}. \quad (8)$$

However, this model might underestimate the radius of the blackbody. Theoretical work considering various physical processes in super-critical disks (for example wind and advection; e.g., Lipunova 1999; P07) have provided more realistic prescriptions for the temperature and radius of the photosphere from which a soft blackbody component originates in ULXs (see review in Fabrika et al. 2021). In P07 the radius of the photosphere was estimated as follows:

$$R_{\text{ph}} = 3 \frac{e_w}{\zeta \beta} \dot{m}^{\frac{3}{4}} R_{\text{in}}, \quad (9)$$

where R<sub>in</sub> is the inner disk radius, and e<sub>w</sub>, ζ, and β are model parameters (see below). We adopted an inner radius equal to three Schwarzschild radii for BHs (assuming nonspinning BHs):

$$R_{\text{in}} = \frac{6GM_{\odot}}{c^2} m, \quad (10)$$

and 11 km for neutron stars (NSs; Özel et al. 2016). The T<sub>eff</sub> is

$$T_{\text{ph}} = 9.28 \times 10^6 \text{ K} \left( \frac{\zeta \beta}{\varepsilon_w} \right)^{\frac{1}{2}} m^{-\frac{1}{4}} \dot{m}^{-\frac{3}{4}}, \quad (11)$$

where ζ = √(ξ² - 1) is a function of the ratio of the perpendicular velocity of the ejected gas over the orbital velocity (ξ = u<sub>z</sub>/u<sub>k</sub>), and β is the ratio of the mean wind radial velocity over the orbital velocity at the R<sub>sp</sub>. Due to energy constraints (cf. P07), β ≈ 1 ≈ ζ, and therefore we varied only e<sub>w</sub>, which is the fraction of the accretion energy powering the outflow. We adopted two values for e<sub>w</sub>, 0.3 and 0.7, to study the effect of this parameter.

## 2.3. The disk component of the ULX spectrum

For all models, the soft blackbody spectrum and its normalization were computed using Eqs. (7)–(11). The remaining emission was modeled by a multicolor disk blackbody, using the xspec model diskbb from sherpa (Doe et al. 2007), which was parametrized by the inner temperature of the accretion disk, computed using the formula from P07:

$$T_{\text{in}} = 1.6 \times m^{-\frac{1}{4}} \left(1 - 0.2 \dot{m}^{-\frac{1}{3}}\right) \text{keV}. \quad (12)$$

# 3. Results

## 3.1. Model predictions for the ionizing photons from ULXs

For a grid of observed luminosities (2 × 10<sup>39</sup>–10<sup>41</sup> erg s<sup>-1</sup>), and for three different CO masses (1.4, 8, and 20 M<sub>☉</sub>), using Eqs. (1)–(12), we computed the ṁ and b, as well as the photosphere radius and temperature for three cases: the SS73 model

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Solutions of Eqs. (1)–(12) for different CO masses (see legend) with L<sub>obs</sub> in the range 2–100 × 10<sup>39</sup> erg s<sup>-1</sup>. From top to bottom: (i) the mass-transfer rate (ṁ) of a ULX as a function of the isotropic-equivalent luminosity when observed face-on (L<sub>obs</sub>). (ii) The beaming factor (b) as a function of L<sub>obs</sub>; the flattening occurs at ṁ = 8.5. (iii) The radius of the region producing the soft component as a function of ṁ for different models (line styles; see lower legend). (iv) The effective temperature of the blackbody region as a function of the ṁ; the dotted-dashed line shows the inner temperature of the disk.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** He II ionization from ULXs with different CO masses (colors; see top legends) and spectrum models (line styles; see bottom legends), as a function of the face-on luminosity of the source. Top: fraction of the photons emerging from the blackbody (bb) component with respect to the total (sum of the blackbody and multicolor disk). Bottom: rate of He II-ionizing photons from the blackbody and disk component in the range 54.4–300 eV.</figcaption></figure>

and the P07 model with two values for the parameter $e_w$ (see Fig. 1). As expected, the $T_{ph}$ from the P07 model is lower than the $T_{eff}$ from SS73 because of the larger photosphere radius compared to the $R_{sp}$. We find that the total energy of the soft component from P07 is $\sim L_{Edd}/3$.

For three different ionization potentials, $E_{ion}$ (13.6 eV for H I, 54.4 eV for He II, and 96.6 eV for Ne V), we computed the rate of ionizing photons $Q_{ion} = \int_{E_{ion}}^{E_{max}} N_E(E) dE$, where $N_E$ is the number of emitted photons with energy between E and $E + dE$, and $E_{max} = 300$ eV. In Fig. 2 we show the number of He II-ionizing photons as a function of the observed luminosity, as well as the relative contribution of the blackbody component. The number of H I and Ne V-ionizing photons are higher by 0–5%, and lower by 5–20%, respectively, depending on the observed luminosity and the model.

## 3.2. Comparison between ULX and stellar populations

We computed the rate of ionizing photons from underlying stellar populations using BPASS (Eldridge et al. 2017) spectra for different metallicities (see Fig. 3), assuming continuous star-formation for 100 Myr. We adopted the BPASS value for the solar metallicity, $Z_{\odot} = 0.02$.

The ionizing power of a ULX population in a given galaxy depends on the number of ULXs and their individual spectra. Our results on the rate of ionizing photons from ULXs as a func-

tion of the observed X-ray luminosity (Sect. 3.1), in combination with XLFs (e.g., Lehmer et al. 2019), provide a handle on the ionizing power of ULX populations.

To get a representative galaxy sample, such as the one used in the recent ULX demographic study of Kovlakas et al. (2020), we selected all star-forming galaxies within a distance of 40 Mpc, with reliable SFR and metallicity estimates in the Heraklion Extragalactic Catalogue (HECATE; Kovlakas et al. 2021). Since HECATE uses infrared SFR indicators and metallicities based on optical emission-line ratios, which are biased in the case of passive galaxies or in the presence of active galactic nuclei, we only selected the 1061 objects without nuclear activity and SDSS colors $g - r < 0.65$ mag (with a cut in the uncertainties in the photometry: $e_g, e_r < 0.1$ mag; Kyritsis et al., in prep.). We used the “homogenized SFR indicator” which combines five different SFR indicators in HECATE. The SFR in this sample spans from $3.9 \times 10^{-2}$ to $9.9 M_{\odot} yr^{-1}$, which is consistent with the range of SFRs in the calibration of the homogenized SFR indicator in HECATE (see Fig. 5 in Kovlakas et al. 2021), but most importantly it falls in the region where the indicator is linearly correlated with SED-based SFR estimates for galaxies with $g - r < 0.65$ (see Fig. 6 in Kovlakas et al. 2021). The metallicity<fnref n="2" /> in our sample is in the range [0.0025, 0.032].

For each galaxy, we sampled the XLF (see Eq. (5)) above $2 \times 10^{39}$ erg s$^{-1}$. Each XLF was scaled for the SFR of the galaxy and for its metallicity using the fit from Douna et al. (2015). However, we note that both the normalization and the shape of the HMXB XLF might depend on the metallicity (Lehmer et al. 2021).

The XLFs were also corrected for geometrical beaming. Under the beaming scenario, not all ULXs are observed due to the nonfavorable angles towards us. While this has no effect on the integrated observed X-ray luminosity (the number of the sources is decreased by a factor of b, but their observed luminosity increases by b as well), the intrinsic population of ULXs was underestimated. For a given CO mass, the beaming factor is a function of the luminosity (cf. Fig. 1), and therefore the “unbeamed” XLF is the following:

$$ \frac{dN_{corr}(L_{38})}{dL_{38}} = b^{-1}(L_{38}) \frac{dN(L_{38})}{dL_{38}}. \qquad (13) $$

We sampled each galaxy’s XLF 1000 times to quantify the stochasticity characterizing the high-end of the XLF. For each galaxy, we estimated the expected number of ULXs by integrating the XLF above $L_{lim}$, and we used this value as the mean of the Poisson distribution (typically 0.1–10 ULXs depending on the properties of the galaxies) from which we sampled the 1000 numbers of ULXs corresponding to each iteration. Then, for each iteration, we sampled from the XLF to get a list of luminosities. Using the results from Sect. 3.1, we summed the $Q_{ion}$ of the ULXs (based on their observed luminosity; see bottom panel of Fig. 2), in each iteration and galaxy, for all ionization potentials, CO masses, and models. In order to show the trend with the metallicity, we binned the galaxies by metallicity with bin edges at metallicity values that correspond to BPASS results 0.002, 0.003, 0.004, 0.006, 0.008, 0.01, 0.014, 0.02, 0.03, and 0.04, and we computed the average rate of ionizing photons, as well as the 99% percentile.

In Fig. 3 we show the rate of ionizing photons for only three combinations of CO masses and models (everything else falls

<footnote n="2">² We converted the HECATE gas-phase metallicities $12 + \log_{10}(\text{O/H})$ to Z, using $12 + \log_{10}(\text{O/H})_{\odot} = 8.69$ (Asplund et al. 2009) and $Z_{\odot} = 0.02$ (Eldridge et al. 2017).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Rate of He II (**left**) and Ne V (**right**) ionizing photons from ULX populations in galaxies as a function of the gas-phase metallicity (stellar value indicated with vertical magenta lines). The black lines show the prediction from stellar populations using BPASS models. The colored solid lines depict the mean value from ULXs assuming three different combinations of CO masses and spectrum models (see legend), covering the full range of our (nine) estimates. The dashed lines show the upper 99% limit to give a sense of the scatter due to the stochastic nature of ULXs.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Average intrinsic spectrum from ULX populations assuming different CO masses (color) and spectrum models (line style), as well as stellar populations from BPASS (black line for solar metallicity and gray for one-tenth of the solar metallicity). All spectra were normalized by the SFR and the metallicity (Z); since we assumed that the XLF scales with $Z_{\odot}/Z$, there is no variation in the ULX spectra, contrary to the stellar populations. The units in the *y* axis are the same as the provided data (see Table 1). The ionization potentials of H I, He II, and Ne V are indicated with vertical lines. The shaded areas denote the regions in the spectrum that we do not consider as ionizing when calculating $Q_{\text{HI}}$, $Q_{\text{HeII}}$, and $Q_{\text{NeV}}$ (see Figs. 2 and 3).</figcaption></figure>

in between), and only for the He II and Ne V (the *Q* for the H I is smaller than the stellar one by ~4 dex). The solid lines correspond to the mean value in the aforementioned metallicity bins, while the dashed ones, to the 99% percentile, show how the stochastic nature of ULXs may lead to He II emitters. Only a small fraction of galaxies, and under a favorable scenario (massive CO, $e_w = 0.7$), host ULX populations with a He II-ionizing

power comparable to that from the stellar populations. Depending on the CO mass and the adopted model, the Ne V-ionizing power of ULX populations can match, or even exceed that of the stellar populations by 2 dex. The scatter of an order of magnitude is caused by the stochastic nature of ULXs and the shape of the XLF. Specifically, in low-metallicity galaxies, which in general are characterized by low SFRs, the ULX content is small (~0.1
<table><caption>**Table 1.** Average ULX spectra scaled by the SFR and metallicity (Z), corresponding to all nine combinations of CO masses (1.4, 8, and 20 M⊙), and blackbody models (SS73, *e*<sub>w</sub> = 0.3, and *e*<sub>w</sub> = 0.7).</caption><thead><tr><th colspan="2">CO mass (M<sub>☉</sub>)</th><th>1.4</th><th>1.4</th><th>1.4</th><th>8.0</th><th>8.0</th><th>8.0</th><th>20.0</th><th>20.0</th><th>20.0</th></tr><tr><th colspan="2">Model</th><th>SS73</th><th>e<sub>w</sub> = 0.3</th><th>e<sub>w</sub> = 0.7</th><th>SS73</th><th>e<sub>w</sub> = 0.3</th><th>e<sub>w</sub> = 0.7</th><th>SS73</th><th>e<sub>w</sub> = 0.3</th><th>e<sub>w</sub> = 0.7</th></tr><tr><th>E<sub>ph</sub><br>(eV)</th><th>λ<br>(Å)</th><th colspan="9">Decimal logarithm of photon flux in (λ, λ + 1Å) bin, scaled by SFR and metallicity<br>(photon s<sup>-1</sup> Myr<sup>-1</sup> yr × Z/Z<sub>☉</sub>)</th></tr></thead><tbody><tr><td>0.123984</td><td>100 000</td><td>32.11996</td><td>32.20639</td><td>32.21280</td><td>33.28700</td><td>33.39010</td><td>33.39158</td><td>33.84989</td><td>33.98321</td><td>33.98385</td></tr><tr><td>0.123986</td><td>99 999</td><td>32.11997</td><td>32.20640</td><td>32.21281</td><td>33.28701</td><td>33.39010</td><td>33.39159</td><td>33.84989</td><td>33.98321</td><td>33.98386</td></tr><tr><td>0.123987</td><td>99 998</td><td>32.11998</td><td>32.20640</td><td>32.21281</td><td>33.28702</td><td>33.39011</td><td>33.39160</td><td>33.84990</td><td>33.98322</td><td>33.98386</td></tr><tr><td>0.123988</td><td>99 997</td><td>32.12506</td><td>32.21149</td><td>32.21790</td><td>33.29210</td><td>33.39519</td><td>33.39668</td><td>33.85498</td><td>33.98830</td><td>33.98895</td></tr><tr><td>0.123989</td><td>99 995</td><td>32.11999</td><td>32.20642</td><td>32.21283</td><td>33.28703</td><td>33.39012</td><td>33.39161</td><td>33.84991</td><td>33.98323</td><td>33.98387</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** The first and second columns are the photon energy (*E*<sub>ph</sub>) and wavelength (λ), respectively, while the rest are the decimal logarithm of the photon flux in the (λ, λ + 1Å) bin for each combination of parameters. In other words, SFR × (Z/Z<sub>☉</sub>) × Σ<sub>i=1</sub><sup>10<sup>4</sup></sup> 10<sup>4</sup> E<sub>ph,i</sub>, where *i* is the row index and *y*<sub>*i*</sub> is a given model column, gives the total intrinsic luminosity (in eV) in the range [1, 10<sup>4</sup>] Å for a ULX population in an underlying stellar population of a given SFR and metallicity. The full table with 100 000 entries is available at the CDS.

ULXs) and therefore they often do not host a ULX. On the other hand, actively star-forming galaxies, with higher metallicity, are expected to host ≥1 ULXs, but due to the shallow slope of the XLF, their luminosities and ionizing power covers a wide range of values.

### 3.3. Average spectrum of ULX populations

While in the previous section we focused on the ionizing power of stochastic ULX populations hosted in individual galaxies, here we construct a SFR-scaled and metallicity-dependent average spectrum of ULXs in star-forming galaxies that can be used as input to IGM-heating and cosmic X-ray background studies (e.g., Upton Sanderbeck et al. 2018). To do so, we summed the XLF-weighted theoretical spectra for ULXs with $2 \times 10^{39}$ erg s<sup>-1</sup> < *L*<sub>obs</sub> < $10^{41}$ erg s<sup>-1</sup>. We also normalized the spectra for the metallicity since we assumed that the XLF is linearly anticorrelated with the metallicity (e.g., Douna et al. 2015). Figure 4 shows the spectra (see Table 1 for the data) along with the stellar models from BPASS for two different metallicity values.

### 3.4. Average effect on IGM heating

Using the average spectra in Sect. 3.3, we quantified the effect of ULX populations in IGM-heating at the epoch of cosmic heating. We focused on the 0.3–1 keV part of the intrinsic spectrum since at lower energies the photons are absorbed by the ISM of the host galaxies, whereas higher energy photons are expected to penetrate the IGM and not significantly contribute in the heating. In addition, we studied the effect of the metallicity-dependent photo-electric absorption of the ISM by applying the model vphabs from xspec for different values of hydrogen column density (*n*<sub>H</sub>).

For all CO masses, the soft X-ray luminosity as a function of metallicity weakly depends on the spectrum model (SS73, P07; within a factor two). In Fig. 5 we show the results for the two extreme CO masses (1.4 and 20 M☉) and for *e*<sub>w</sub> = 0.3. The normalization is roughly linear to the accretor mass since the models predict harder spectra for lower CO masses (the blackbody component is nearly Eddington-limited). While the exact distribution of CO masses is important to constrain the soft X-ray part of average ULX spectra, we find that the blackbody and disk com-

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

**Fig. 5.** Average soft X-ray luminosity (*L*<sub>0.3–1 keV</sub>) of ULX populations as a function of the host galaxy’s metallicity (solar value indicated with a vertical magenta line), assuming CO masses of 20 M☉ (top, intense lines) and 1.4 M☉ (bottom, pale lines) and spectrum described by the P07 model with *e*<sub>w</sub> = 0.3. Black lines correspond to the intrinsic (unabsorbed) X-ray luminosity, while blue and orange lines account for the absorption in the galaxy’s ISM, with *n*<sub>H</sub> = $10^{21}$ and $10^{22}$ cm<sup>2</sup>, respectively. Dashed, dotted, and solid lines indicate the blackbody component, disk component, and their sum, respectively. (For 20 M☉, the black dashed and dotted lines coincide, indicating equal component contributions.) The gray band denotes the $L_{X<2 keV}$ 68% confidence interval constrained from the cosmic 21-cm signal for galaxies at z ~ 8 (Abdurashidova et al. 2022).

ponents have a comparable contribution, with a factor of 1.1–9.3 underestimation if we neglect the blackbody component.

## 4. Discussion

Assuming that the spectra of ULXs are the sum of a blackbody (using the SS73 and P07 models) and a disk component, we quantified the ionizing power of individual ULXs for different CO masses. We find a wide range in *Q*<sub>HeII</sub> ~ $10^{45–48}$ photon s<sup>-1</sup>, which is considerably lower than the rates from He-emitters such as the IZw 18 galaxy, ~1.33×10<sup>50</sup> photon s<sup>-1</sup> (Kehrig et al. 2015). The ULX in this galaxy has been investigated as the origin of the ionizing radiation (e.g., Schaerer et al. 2019); however, difficulties in constraining the EUV part of the spectrum
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Disk (thick lines) and blackbody (BB; thin lines) components of the intrinsic spectra of a ULX with $L_{\text{obs}} = 10^{40} \text{ erg s}^{-1}$ for two different CO masses, 1.4 $M_{\odot}$ (orange) and 20 $M_{\odot}$ (blue), and two different values of $e_w$ for the P07 model (line styles, see legend). For comparison, in purple, we overplotted the DIS (dotted line) and BMC (dashed line) models from Simmonds et al. (2021), as well as the model for a BH of mass 10 $M_{\odot}$ (continuous line) from Senchyna et al. (2020), scaled for the same source luminosity. We note that in our models, the normalization of the intrinsic spectra depends on the CO mass: the smaller the CO mass is, the smaller the beaming factor is for the same observed, face-on luminosity.</figcaption></figure>

of a ULX (e.g., Simmonds et al. 2021) might result in stark differences in its ionization power. Our theoretical ULX SEDs are not capable of producing strong He II ionization ($Q_{\text{HeII}}/L_X \le 4 \times 10^7 \text{ photon erg}^{-1}$) as invoked in previous studies (e.g., $2 \times 10^{10} \text{ photon erg}^{-1}$; Schaerer et al. 2019).

In Fig. 6, we compare the model spectra with previous studies for a ULX with $L_{\text{obs}} = 10^{40} \text{ erg s}^{-1}$. Specifically, we show the spectrum of the disk and the blackbody component of a ULX with an observed (face-on) luminosity of $10^{40} \text{ erg s}^{-1}$ for two different CO masses, and two different values of $e_w$ in the P07 model. For the same observed, isotropic-equivalent luminosity, different CO masses correspond to different beaming factors. For this reason, in the case of the 1.4 $M_{\odot}$ accretor ($b \approx 0.076$), the bolometric luminosity of the source is lower than in the case of the 20 $M_{\odot}$ accretor ($b \approx 0.81$), leading to a lower normalization for the former. While the shape of the blackbody component is affected by the adopted value for $e_w$, the disk component is the same for a given CO mass (the corresponding lines in Fig. 6 are not repeated for different values of $e_w$.) We also show the models from Simmonds et al. (2021) and Senchyna et al. (2020), which have been scaled at the same luminosity, the normalization of which do not depend on the beaming factor since geometrical collimation is not considered in these studies. We note

that these models have been overplotted to allow for qualitative comparisons against recent studies of the contribution of ULXs and XRBs in nebular emission. Our individual ULX SEDs show that the photon flux at ~54.4 eV varies by one order of magnitude depending on the parameters of the models. However, they are weaker by many orders of magnitude compared to the literature models. We should stress again that no direct observations exist in this range and, in all cases, the ULX SEDs are extrapolated from theoretical or empirical models, calibrated to higher energies (≥300 eV).

When considering the ionizing power of ULX populations, the geometrical beaming model results into two opposing effects. On the one hand, it results in the bolometric luminosity of the sources and, consequently, their ionization power being lower than what is inferred from the observed luminosities. On the other hand, a fraction of the ULXs cannot be observed in the X-rays, but nevertheless ionizes the ISM. While the X-ray budget is preserved, the picture is more complicated in the EUV part of the spectrum: its shape and the beaming factor are correlated through their dependence on the *m*. Consequently, the use of observed spectra and population synthesis techniques (e.g., Fragos et al. 2013b) might underestimate the ionizing power of ULXs. Using the beaming-corrected HMXB XLF to anchor our
analysis, we derived the underlying distribution of *m* and calculated the ionizing power of ULX populations in a realistic sample of galaxies. We find that in only a small fraction (~1%) of galaxies do the ULXs compete with the underlying stellar populations in He II emission, which is in agreement with recent studies (e.g., Senchyna et al. 2020), while they are significant in the case of Ne V (e.g., Simmonds et al. 2021). On the other hand, observational studies may also put constraints on the contribution of the hot gas component in star-forming galaxies, which is found to be comparable to the stellar component in ionizing the He II in low-metallicity galaxies (e.g., Lehmer et al. 2022). From our modeling, the Ne V-ionizing power of ULXs points at the possibility of constraints for the EUV emission from ULXs with the use of emission lines associated with high-ionization potentials. We encourage the study of high-ionization emission line galaxies, their emission processes, and in general systematic studies constraining the ionizing power and spectra of ULXs (e.g., Izotov et al. 2021).

In addition, we computed the average spectrum of the ULX population in a galaxy normalized by its SFR and metallicity, for the different theoretical spectra, and adopted CO masses. These spectra can be used as input for IGM-heating studies and be compared to observational constraints from the cosmic 21-cm signal (e.g., Abdurashidova et al. 2022). As an example, we calculated the heating power of ULX populations via their soft X-ray emission. We find that the contribution from the blackbody component is comparable to the disk component, highlighting the importance of the former in quantifying the effect of ULXs in the early Universe. However, the normalization in this case strongly depends on the CO mass, and only weakly in the spectrum model among the ones investigated here.

In conclusion, the dependence of our estimates on the shape of the spectra, the stochastic nature of ULXs, and the CO mass distribution result in a 2 dex scatter in the ionizing and heating power of ULX populations. Despite the success of analytical models in interpreting key properties of accreting sources, hydrodynamical simulations are necessary for investigating the radiative and mechanical feedback in super-critical accretion disks (e.g., Sądowski & Narayan 2015). Due to the computational cost of these simulations, only a select number of fiducial systems have been investigated. We encourage such efforts to continue, especially in the case of NS accretors where the presence of magnetic fields is important, with some sources showing a lack of beaming (e.g., Binder et al. 2018). Combining accurate spectral models with detailed binary population synthesis techniques (e.g., Fragos et al. 2022) will provide stringent constraints on the ionizing power of ULXs, and their contribution in the heating of the early Universe.

**Acknowledgements.** We would like to thank the anonymous referee for critical comments which improved the paper. This work was supported by the Swiss National Science Foundation Professorship Grant (PP00P2_176868; PI Fragos). K.K. acknowledges support from the Federal Commission for Scholarships for Foreign Students for the Swiss Government Excellence Scholarship (ESKAS No. 2021.0277).

## References

* Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 924, 51
* Abramowicz, M. A., Calvani, M., & Nobili, L. 1980, ApJ, 242, 772

* Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

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

* Binder, B., Levesque, E. M., & Dorn-Wallenstein, T. 2018, ApJ, 863, 141

* Brice, N., Zane, S., Turolla, R., & Wu, K. 2021, MNRAS, 504, 701

* Doe, S., Nguyen, D., Stawarz, C., et al. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, ASP Conf. Ser., 376, 543

* Douna, V. M., Pellizza, L. J., Mirabel, I. F., & Pedrosa, S. E. 2015, A&A, 579, A44

* Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058

* Fabrika, S. N., Atapin, K. E., Vinokurov, A. S., & Sholukhova, O. N. 2021, Astrophys. Bull., 76, 6

* Feng, H., Tao, L., Kaaret, P., & Grisé, F. 2016, ApJ, 831, 117

* Fragos, T., Lehmer, B., Tremmel, M., et al. 2013a, ApJ, 764, 41

* Fragos, T., Lehmer, B. D., Naoz, S., Zezas, A., & Basu-Zych, A. 2013b, ApJ, 776, L31

* Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2022, ApJS, submitted [arXiv:2202.05892]

* Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics, 3rd edn. (Cambridge: Cambridge University Press)

* Garnett, D. R., Kennicutt, R. C.Jr., Chu, Y. H., & Skillman, E. D. 1991, ApJ, 373, 458

* Götberg, Y., de Mink, S. E., & Groh, J. H. 2017, A&A, 608, A11

* Izotov, Y. I., Thuan, T. X., & Privon, G. 2012, MNRAS, 427, 1229

* Izotov, Y. I., Thuan, T. X., & Guseva, N. G. 2021, MNRAS, 508, 2556

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

* Kajava, J. J. E., & Poutanen, J. 2009, MNRAS, 398, 1450

* Kehrig, C., Vílchez, J. M., Pérez-Montero, E., et al. 2015, ApJ, 801, L28

* Kehrig, C., Guerrero, M. A., Vílchez, J. M., & Ramos-Larios, G. 2021, ApJ, 908, L54

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

* King, A. R., & Pounds, K. A. 2003, MNRAS, 345, 657

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

* Kovlakas, K., Zezas, A., Andrews, J. J., et al. 2020, MNRAS, 498, 4790

* Kovlakas, K., Zezas, A., Andrews, J. J., et al. 2021, MNRAS, 506, 1896

* Lehmer, B. D., Eufrasio, R. T., Tzanavaris, P., et al. 2019, ApJS, 243, 3

* Lehmer, B. D., Eufrasio, R. T., Basu-Zych, A., et al. 2021, ApJ, 907, 17

* Lehmer, B. D., Eufrasio, R. T., Basu-Zych, A., et al. 2022, ApJ, 930, 135

* Lipunova, G. V. 1999, Astron. Lett., 25, 508

* Madau, P., & Fragos, T. 2017, ApJ, 840, 39

* Mapelli, M., Ripamonti, E., Zampieri, L., Colpi, M., & Bressan, A. 2010, MNRAS, 408, 234

* Mesinger, A., Ferrara, A., & Spiegel, D. S. 2013, MNRAS, 431, 621

* Mezcua, M., Civano, F., Marchesi, S., et al. 2018, MNRAS, 478, 2576

* Middleton, M. J., Walton, D. J., Alston, W., et al. 2021, MNRAS, 506, 1045

* Nanayakkara, T., Brinchmann, J., Boogaard, L., et al. 2019, A&A, 624, A89

* Olivier, G. M., Berg, D. A., Chisholm, J., et al. 2021, ApJ, submitted [arXiv:2109.06725]

* Oskinova, L. M., & Schaerer, D. 2022, A&A, 661, A67

* Özel, F., Psaltis, D., Güver, T., et al. 2016, ApJ, 820, 28

* Pinto, C., Alston, W., Soria, R., et al. 2017, MNRAS, 468, 2865

* Plat, A., Charlot, S., Bruzual, G., et al. 2019, MNRAS, 490, 978

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

* Power, C., James, G., Combet, C., & Wynn, G. 2013, ApJ, 764, 76

* Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49

* Rickards Vaught, R. J., Sandstrom, K. M., & Hunt, L. K. 2021, ApJ, 911, L17

* Sądowski, A., & Narayan, R. 2015, MNRAS, 453, 3213

* Saxena, A., Pentericci, L., Schaerer, D., et al. 2020, MNRAS, 496, 3796

* Schaerer, D., Fragos, T., & Izotov, Y. I. 2019, A&A, 622, L10

* Senchyna, P., Stark, D. P., Mirocha, J., et al. 2020, MNRAS, 494, 941

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

* Simmonds, C., Schaerer, D., & Verhamme, A. 2021, A&A, 656, A127

* Stephenson, C. B., & Sanduleak, N. 1977, ApJS, 33, 459

* Sutton, A. D., Roberts, T. P., & Middleton, M. J. 2013, MNRAS, 435, 1758

* Umeda, H., Ouchi, M., Nakajima, K., et al. 2022, ApJ, 930, 37

* Upton Sanderbeck, P. R., McQuinn, M., D’Aloisio, A., & Werk, J. K. 2018, ApJ, 869, 159

* Waisberg, I., Dexter, J., Olivier-Petrucci, P., Dubus, G., & Perraut, K. 2019, A&A, 624, A127