© The Authors 2026

1 Introduction

The mass distribution of observed black holes features a population of stellar-mass black holes at the low-mass end and a population of supermassive black holes (SMBHs) at the high-mass end. In the framework of hierarchical structure formation, supermassive black holes are expected to form through repeated merging and accretion of matter from a seed population of lighter (stellar-mass) black holes. In this picture, objects in the intermediate-mass range with 10²–10⁵ M_⊙ are naturally expected to exist, but they remain to be conclusively identified. However, there are several observations that are possibly explained by the existence of intermediate-mass black holes (IMBHs). The following main points of motivation are often discussed.

1. Observations of larger-than-expected velocity dispersions in select globular clusters hint at central IMBHs “stirring up” the stars and introducing additional contributions to the velocity dispersion (Frank & Rees 1976). An analysis of the stellar kinematics of the globular cluster ω Centauri by Noyola et al. (2008) reveals an increase in the mass-to-light ratio toward the center of the cluster and suggests the existence of a central black hole (BH) with a mass of 4 × 10⁴ M_⊙. A study by Häberle et al. (2024) has identified seven fast-moving stars in the central 3 arcseconds of ω Centauri, with velocities exceeding the escape velocity of the cluster, indicative of a central BH with a lower bound mass of 8200 M_⊙.

2. Another point of motivation comes from the observation of ultra-luminous X-ray sources (ULXs). These are sources that emit X-rays with luminosities of L_X ~ 10³⁹–10⁴⁰erg/s, exceeding the luminosities of accreting stellar BHs at the Eddington limit while being below the luminosities generated by typical active galactic nuclei (AGN) by a margin of several orders of magnitude. Common explanations for ULXs are given by X-ray binaries (King et al. 2001; Lasota & King 2023), although accreting IMBHs provide an alternate explanation for select sources (Miller et al. 2003).

3. Possibly the strongest point of motivation for the existence of IMBHs comes from the rapid growth of SMBHs and the great masses they managed to acquire in the early universe. Bogdán et al. (2024) reported the observation of the high redshift quasar UHZ1 with an estimated mass of ~10⁷–10⁸ M_⊙ using the Chandra X-ray observatory. The object was later spectroscopically confirmed to have a redshift of z = 10.1 using the James Webb Space Telescope (JWST) (Goulding et al. 2023; Natarajan et al. 2024). With having only a few 100 Myr to grow, IMBHs with ~10⁴–10⁵ M_⊙ are thought to act as seed BHs in the early universe that can reach supermassive scales by mere accretion of matter.

Irrespective of mass, BHs are expected to modify the dark matter (DM) distribution and cause the adiabatic growth of a DM overdensity, called a “spike”. If DM is self-annihilating, the resulting gamma-ray luminosity of the spike will be enhanced compared to an unmodified DM mass density profile (Zhao & Silk 2005; Bertone et al. 2005), increasing the reach for indirect searches with gamma-ray telescopes as demonstrated by HESS Collaboration (2008). Indications for the existence of these types of overdensities are drawn from the rapid orbital decay rates observed in two low-mass X-ray binary systems (LMXBs), which are suggestive of the companion stars experiencing dynamical friction in a DM spike environment (Chan & Lee 2023). With IMBHs dressed with a DM spike providing another way to probe DM, there is interest in studying the IMBH population itself (Aschersleben et al. 2024).

This work aims to estimate the number of IMBHs in Milky Way-like galaxies that might be dressed with DM spikes. Specifically, unmerged IMBHs are of interest, because the spike might be disturbed or disrupted during a merger with another BH. With the merging behavior of stellar binary BHs (BBHs) studied observationally using gravitational wave detectors, they provide constraints on the BH-forming stellar population. For this reason, the presented work considers the formation of 10²–10³ M_⊙ “light seed” IMBHs following the collapse of high-mass Population III (Pop III) stars. In addition to being IMBH progenitors, Pop III stars are expected to be short-lived and luminous and so might have left detectable footprints in present-day stellar observables, although latest efforts have deemed this possibility unlikely (Sun et al. 2021). It is therefore tested to what extent an IMBH-forming population of Pop III stars contributes to the optical and infrared bands of the diffuse cosmic background radiation called the extragalactic background light (EBL) (Hill et al. 2018).

The methods for simulating Pop III stars and for assessing the merging behavior of their remnant BBHs are discussed in Sect. 2. The primary result on the number of unmerged IMBHs in Milky Way-like galaxies and the implications of the required Pop III stellar population on stellar observables are presented in Sect. 3. In Sect. 4, the caveats and uncertainties of this work are discussed and compared with previous work conducted on Pop III IMBHs.

2 Methods

2.1 Simulation of Population III stellar numbers

In order to estimate the number of unmerged IMBHs from Pop III stars in Milky Way-like galaxies, it is necessary to track the evolution of their progenitors. For this purpose, the A-SLOTH code has been used (Hartwig et al. 2022). It is a highly efficient and parallelized semi-analytical model (SAM) that is calibrated with six local observables of the Milky Way. Sophisticated implementations of feedback mechanisms governing star formation and treatment of metallicity allow for individual sampling and tracking of Pop III and Pop II stars until their supernova phases. In this work, A-SLOTH Version 1.2.1 was used¹. Since in its default configuration, the simulation code aims to reproduce the Milky Way, the set of simulation parameters was largely left unchanged. Following recent numerical simulations (Klessen & Glover 2023), the Pop III IMF was chosen to be log-flat with the following functional form

$$\xi (m)=\frac{\mathrm{d}N}{\mathrm{d}m}=N_0{m}^{-\gamma },$$(1)

where γ = 1 is the power-law index and N₀ is a normalization factor such that

$$1={\int }_{m_{min}}^{m_{max}}m\xi (m)\mathrm{d}m.$$(2)

To allow the formation of high-mass BHs, the stellar mass range was set to 5 M_⊙–500 M_⊙ and the simulation system mass was updated to match the latest measurement of the Milky Way mass, which is taken to be the upper limit of ~2.5 × 10¹¹ M_⊙ determined by Jiao et al. (2023). An addition has been made to the source code that writes formation time, mass, number, and relative probability of occurrence, according to the logarithmically binned IMF, into an output file, which allows tracking of the whole Pop III stellar population and further post-processing. In total, 100 simulations have been run with the same parameter set such that fluctuations in star formation emerging from differences in initial conditions could be mitigated. The star formation rates (SFRs) from the individual simulations and their average are shown in Fig. 1. The parameterization of the average, which will from now on be called the A-SLOTH SFR, is given by

$${\psi }_A(z)={\psi }_2(z)+{\psi }_4(z),$$(3)

with

$${\log }_{10}\left(\frac{{\psi }_2(z)}{{\mathrm{M}}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}{\mathrm{M}\mathrm{p}\mathrm{c}}^{-3}}\right)=a+b{z}^{\mathrm{\prime }}+c{z}^{\mathrm{\prime }2},$$(4)

$${\log }_{10}\left(\frac{{\psi }_4(z)}{{\mathrm{M}}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}{\mathrm{M}\mathrm{p}\mathrm{c}}^{-3}}\right)=a+b{z}^{\mathrm{\prime }}+c{z}^{\mathrm{\prime }2}+d{z}^{\mathrm{\prime }3}+e{z}^{\mathrm{\prime }4},$$(5)

where z′ = z + z₀. The parameterization fit parameters are given in Table 1. The parameterization is purposefully cut off at both high and low redshifts to account for the start of Pop III star formation and its termination. While the starting redshift results from the analytical prescriptions implemented in the A-SLOTH code, the termination redshift is inferred from the absence of Pop III stars in our cosmic neighborhood. This indicates that Pop III star formation has effectively stopped. For this work, the formation of metal-free stars is assumed to occur until the redshift of the closest observed candidate object that indicates the existence of Pop III stars. This is the Lynx Arc, a strongly gravitationally lensed structure at redshift z = 3.357 (Fosbury et al. 2003)².

For comparison, the A-SLOTH SFR is viewed along with two other Pop III SFRs in Fig. 2: the SFR presented by Jaacks et al. (2019) results from a hydrodynamical simulation of the early universe that was halted at redshift z ~ 7. It suggests that robust ongoing Pop III star formation is possible at z < 7 in the presence of Pop II stars. The parametrized Jaacks SFR follows

$$\begin{array}{rl}& {\log }_{10}\left(\frac{{\psi }_J(z)}{{\mathrm{M}}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}{\mathrm{M}\mathrm{p}\mathrm{c}}^{-3}}\right)\\ & =-3.7626+1.5430\times {10}^{-1}z-9.7130\times {10}^{-3}{z}^2.\end{array}$$(6)

The other Pop III SFR shown in Fig. 2 is a simplified parametrization of the Skinner & Wise (2020) Pop III SFR that was adopted by Tanikawa et al. (2022) in their study of reproducing the Gravitational Wave Transient Catalog 2 (GWTC-2) BBH merger rate. It shall be referred to as the SWT SFR. Shown in Fig. 2 is also the best fit model of total star formation by Madau & Fragos (2017) constructed from multi-wavelength measurements in the redshift range z = 0–10 that will be called the MF17 SFR.

Fig. 1 SFRs from 100 A-SLOTH simulations (blue and green lines). The averaged SFR is shown in white with the overlaid parametrization in red.

Fig. 2 Comparison of Pop III A-SLOTH, Jaacks, and SWT SFRs with MF17. A-SLOTH Pop III star formation is terminated at z = 3.357.

2.2 Remnants of Population III stars

The evolution of stars and the type of remnant they leave behind is highly dependent on their mass. Zero-metallicity progenitor stars ≲20 M_⊙ are generally incapable of forming a BH remnant (Woosley et al. 2002). For stars above ≳20 M_⊙ BH formation is possible either by fallback or direct collapse (Fryer et al. 2012), where in the latter case the full stellar mass is assumed to be retained, as metal-free stars do not suffer significant mass loss from intense winds (Belczynski et al. 2016). Stars with helium cores above ~30–50 M_⊙, experience a phenomenon known as pair instability (Woosley 2017), where photons produced in the core of the star are sufficiently energetic to produce electronpositron pairs via γγ ⟶ e⁺e⁻. The subsequent loss of radiation pressure causes an explosive ignition of fusion to counteract the stellar collapse leading to a recurring pulsational pair-instability supernova (PPISN) or, if the first pulse is sufficiently energetic to disrupt the star, a pair-instability supernova (PISN). A PPISN leaves behind a BH remnant with a mass given by the maximal stable helium-core mass, whereas a PISN does not leave behind a remnant at all.

To create cores with ~30–50 M_⊙, the zero age main sequence (ZAMS) masses range from ~70 M_⊙ (Woosley 2017) to ~95 M_⊙ (Tanikawa et al. 2021). In this work, the start of the pair instability is marked by a helium core of 45 M_⊙, equivalent to the ZAMS mass of 95 M_⊙. This model was assumed by Tanikawa et al. (2022) to simulate the Pop III BBH merger rate that is used to determine the number of unmerged IMBHs. PISN occurs for stars in the range 130–260 M_⊙. Beyond ZAMS masses of ~260 M_⊙ photodisintegration of the collapsing core caused by pair instability prevents explosive fusion such that a direct collapse to a BH is possible, again devouring the whole star (Fryer et al. 2001; Renzo et al. 2020). Furthermore, the term “high-mass Pop III” is reserved for masses beyond 260 M_⊙, which are capable of forming IMBH remnants. The different evolutionary paths are indicated in Figs. 3 and 4, which show the remnant mass function and the sampling of Pop III stars by the A-SLOTH code, respectively.

Fig. 3 Remnant mass function of metal-free stars following Fryer et al. (2012). The mass range above 260 M⊙ is not shown, though it simply follows the dotted line (BH retains full progenitor mass). For a comparison of the stellar model with inefficient convective overshoot (this plot) with a model with efficient convective overshoot, see Fig. 2 in Tanikawa et al. (2021).

Fig. 4 Averaged mass distribution of Pop III stars from 100 A-SLOTH simulations. The sampling reflects the log-flat IMF. Coloring indicates different stellar evolutionary paths depending on the ZAMS mass.

2.3 Population III binary black holes

The number of unmerged IMBHs follows from the assessment of BBH mergers that occur in a population of Pop III remnants. The nature and occurrence of BBH mergers have been studied using gravitational wave observations (LIGO Scientific Collaboration, Virgo Collaboration and KAGRA Collaboration 2023), which, however, are limited to lower-mass systems. Instead, this work refers to simulations by Tanikawa et al. (2022). The authors reproduce the BBH merger rate from GWTC-2, from which the merging behavior is extrapolated to higher masses. Tanikawa et al. (2022) account for a great variety of stellar progenitor types including Pop III stars, which allows the viewing of mergers of Pop III BBHs in isolation. They find that the differential merger rate with respect to (w.r.t.) the primary mass above m₁ = 45 M_⊙ is dominated by Pop III remnants. With the simulated merger rate being in agreement with GWTC-2 at high primary masses, this puts a very strict constraint on the amount of allowed Pop III star formation. An SFR that is too large would overproduce the measured BBH merger rate at high primary masses, and therefore one is limited to the SWT SFR used by Tanikawa et al. (2022) or an SFR with equivalent cumulative star formation. The number of Pop III stars simulated using the A-SLOTH code follows from the A-SLOTH SFR that well exceeds SWT in terms of cumulative star formation.

Since star formation is assumed to take place with a constant log-flat IMF throughout cosmic time for both the A-SLOTH and SWT SFR, the A-SLOTH stellar population is simply down-scaled to match the population size that would be expected for the SWT SFR. The required scaling factor is determined by the ratio of the SFR integrals as given by

$$S_{SWT/A}=\frac{{\int }_{z_{min}}^{z_{max}}{\psi }_{SWT}(z)\mathrm{d}z}{{\int }_{z_{min}}^{z_{max}}{\psi }_A(z)\mathrm{d}z}.$$(7)

2.3.1 Effective volume occupied by stellar population

For the BBH merger rate, simulated by Tanikawa et al. (2022), to be applied to the population of remnant BHs from Pop III stars, a certain spatial understanding of the distribution of the initial stellar population is required. Due to A-SLOTH being unable to track the positions of the various objects throughout the simulation, spatial information on the stars in the Pop III catalog and their remnant BHs is unavailable. Thus, one must rely on the SFR to relate the population of stars to a volume in which merger events are expected to occur. Following Madau & Dickinson (2014) the SFR can be translated into a stellar remnant mass density or “dead” mass density (DMD) ρ_†. The conversion is given by

$${\rho }_{†}(z)=D{\int }_z^{z_{max}}\psi \left(z^{\mathrm{\prime }}\right)\left|\frac{\mathrm{d}{t}^{\mathrm{\prime }}}{\mathrm{d}{z}^{\mathrm{\prime }}}\right|\mathrm{d}{z}^{\mathrm{\prime }}.$$(8)

The “dead fraction” D can be understood as the fraction of mass of a stellar population that ends up as stellar remnants and is given by the following equation

$$D={\int }_{m_{min}}^{m_{max}}w(m)\xi (m)\mathrm{d}m,$$(9)

where w(m) is the remnant mass function.

Making the restriction for masses in the calculation of Eq. (9) to include only high-mass Pop III stars (>260 M_⊙), Eq. (8) yields the mass density of IMBHs. This IMBH DMD is further split to obtain individual DMDs for each A-SLOTH high-mass bin i by weighting the IMBH DMD with the ratio of the sampling probability of each bin p_i and the probability of sampling a high-mass star p_HM as given by

$${\rho }_{†,i}(z)=\frac{p_i}{p_{HM}}{\rho }_{†}(z).$$(10)

The complete IMBH DMD and the individual high-mass bin DMDs are illustrated in Fig. 5. The division of each individual DMD by its bin mass then converts it to the number density for each bin. The derived number densities at present day (z = 0) can finally be used to relate the population of stars, and subsequently their BH remnants, to an “effective volume” V_eff which they would occupy under the assumption of a homogeneous spatial distribution of the objects

$$V_{eff,i}=s_{SWT/A}\frac{N_i{m}_i}{{\rho }_{†,i}(z=0)},$$(11)

where N_i is the number of objects in the high-mass bin i and m_i is the bin mass. The factor s_SWT/A is used to scale the population size corresponding to the A-SLOTH SFR to match the SWT SFR. In view of the number densities, V_{eff, i} is the same for all bins.

Fig. 5 Total dead mass density of all IMBH remnants (green) and dead mass densities for the individual IMBH capable bins (blue lines).

2.3.2 Binary black hole combinatorics

The BBH merger rate can be applied to the effective volume V_eff to determine the number of mergers, given the extent of the initial stellar population. A comparison of the number of mergers with the initial population size reveals how many IMBHs remain unmerged. For this, the Pop III BBH merger rate R_T22 from Fig. 1 of Tanikawa et al. (2022) is integrated over redshift, as given by

$$R_{merge}={\int }_{z=0}^{z_{max}}{R}_{T22}(z^{\mathrm{\prime }})\left|\frac{\mathrm{d}{t}^{\mathrm{\prime }}}{\mathrm{d}{z}^{\mathrm{\prime }}}\right|\mathrm{d}{z}^{\mathrm{\prime }}.$$(12)

The number of mergers per unit volume amounts to R_merge = 17.72 Mpc⁻³. This merger number density contains events with all the various pairings of primary and secondary BH masses that are possible. In order to make a statement on the merger number density for events specifically containing IMBHs, an overview of all different BBH mass combinations and their relative likelihoods of occurrence is necessary.

Tanikawa et al. (2022) operate in the scenario where BBH systems arise as a consequence of evolved Pop III stellar binary systems. It follows that the frequency of occurrences is given by the initial stellar pairings where the BH masses in question are determined via the remnant mass function. Naturally, the formation of a Pop III binary with a particular pairing of masses results from the complicated physical processes that govern the fragmentation of their common protostellar disk. Since studying this demands extensive and resource intensive N-body simulations as performed by Stacy & Bromm (2013) and Stacy et al. (2016), a simplified approach is chosen where the occurrence of a mass pairing is solely given by the occurrence of the single stars themselves. One might picture an urn filled with stars of different masses from which two stars are picked at random to comprise a binary system. For our purposes, this urn is to be identified with the population of stars as tracked by A-SLOTH. Since all n_bin mass bins are sampled virtually equally, one ends up with $n_{bin}^2$ combinations, all sharing the same relative occurrence. Fig. 6 illustrates the combinatorics for a simplified case with one bin for every type of stellar evolution (1. BH incapable, 2. BH capable, 3. PPISN, 4. PISN, 5. IMBH capable). Note that only a subset of the 64 mass bins in A-SLOTH are considered for this calculation. This is for the following two reasons:

1. Certain initial stellar binaries do not lead to mergers. For one, Pop III stars with masses < 20 M_⊙ are too light to form BHs in the first place. Hence, all bins with m_i < 20 M_⊙ are excluded from the calculation (all binary combinations containing 1-type stars in Fig. 6 are excluded). The effects of pair-instability demand another augmentation in the number of possible combinations. Stars ending their lives in a PISN do not leave any remnants behind. As such, all combinations, where at least one of the binary partners has a mass 130 M_⊙ ≤ m_i ≤ 260 M_⊙, do not contribute (all binary combinations containing 4-type stars in Fig. 6 are excluded). Finally, systems in which both stars undergo PPISN and leave behind a 45 M_⊙ BH also do not contribute (combination (3, 3) in Fig. 6). This is due to the immense mass loss that stars can experience during a PPISN. Mass loss can significantly shrink the stars and thus eliminate the possibility of stellar interaction, which under normal circumstances would accelerate the rate of inspiral or disrupt the binary altogether. All n_PPISN bins with $n_{PPISN}^2$ pairings of PPISN-capable stars are subtracted from the list of merging combinations.

2. Tanikawa et al. (2022) has defined the upper mass limit of Pop III stars as 150 M_⊙, while this work includes stars with masses up to 500 M_⊙. Considering the effects of pair instability, as discussed in Sect. 2.2, IMBHs cannot be formed in Tanikawa et al. (2022) at all. This implies that the mergers per unit volume of around 17.72 Mpc⁻³ do not account for the mergers of IMBHs. Rather, it can be understood as a reduced number of mergers corresponding to objects <150 M_⊙ (combinations (2, 2) and (2, 3) in Fig. 6).

Statements about the mergers of IMBHs can be made by exploiting the equal relative occurrence of all initial stellar mass combinations. If the merger number density for a single pairing is determined, it can be multiplied by the number of combinations containing IMBH progenitors to yield the number of mergers missing in Tanikawa et al. (2022). Following Fig. 6 and the augmentations to the contributing BBH combinations to the mergers in Tanikawa et al. (2022), the number of relevant initial stellar bins equates to n_bin = 26, with n_PPISN = 4 being the number of bins that yield 45 M_⊙ remnants. The resulting number of combinations contributing to the cumulative mergers per unit volume is

$$n_{comb}=n_{bin}^2-n_{PPISN}^2=660.$$(13)

The individual merger number density for any one of the 660 combinations R₁ is naturally

$$R_1=\frac{R_{merge}}{n_{comb}}\approx 0.03{\mathrm{M}\mathrm{p}\mathrm{c}}^{-3}.$$(14)

One can compute the number of objects per bin “lost” in the mergers throughout cosmic time by multiplying R₁ by the object “loss multiplicity”, that is, the aggregated number of objects from a bin lost in the various combinations due to mergers. Generalizing the simplified scheme shown in Fig 6 and extending the number of bins to n_bin,HM = 35 to include the high-mass stars, in other words, the IMBHs, which were previously missing in n_bin, the loss multiplicity N_loss is given by

$$N_{loss}=2(n_{bin,HM}-1)+2=2n_{bin,HM}.$$(15)

The first term of the sum represents the number of objects lost in combinations in which the binary partners have different masses. The factor 2 is included to account for the two combinatorial pairings that yield the same physical binary (e.g., pairings (5, 3) and (3, 5), while combinatorially unique with equal relative occurrence, both yield a physical system (5, 3) with the combined relative occurrence). The second term of the sum accounts for the pairing of objects with identical masses (the case (5, 5)). In this instance, a merger event will reduce the number of objects in the bin by 2.

Combining Eqs. (11), (14), and (15), the number of merged objects per bin can be obtained and is given by

$$N_{merged,i}=R_1{N}_{loss}{V}_{\text{eff},i}.$$(16)

A sum of Eq. (16) over all high-mass bins returns the number of IMBHs that have been lost to mergers throughout cosmic time

$$N_{merged}=\sum _{i,m_i\ge 260}{N}_{merged,i}.$$(17)

Naturally, the difference between the initial number of IMBHs N_HM and N_merged gives the number of IMBHs that have remained unmerged until the present day

$$N_U=N_{HM}-N_{merged}.$$(18)

Fig. 6 Combinations of binary stellar systems shown for five mass bins. Combinations inside dotted outline simulated in Tanikawa et al. (2022). Combinations inside thin solid outline lead to BBH mergers in Tanikawa et al. (2022). Combinations inside thick solid outline lead to mergers of IMBHs (this work).

3 Results

3.1 Number of unmerged IMBHs

The A-SLOTH simulation yields an average number of high-mass Pop III stars ⟨N⟩ = 2242.85 with a statistical uncertainty of σ = 103.55 per Milky Way-like galaxy arising from 100 simulation runs. In view of the stellar evolutionary paths of high-mass Pop III stars, as discussed in Sect. 2.2, that number is immediately to be identified as the population size of their remnant IMBHs. Compared to the A-SLOTH SFR, SWT yields about a tenth of cumulative star formation following Eq. (7) such that the IMBH population reduces to N^SWT = 222.62 ± 10.28. This shall be regarded as the total number of IMBHs formed over cosmic time before mergers are taken into account.

The number of unmerged IMBHs, as inferred from simulations of the merger rate in Tanikawa et al. (2022) using the combinatorial assessment of BBH mergers in an effective volume, equates to $N_U^{SWT}$ = 131.87 ± 6.09. The uncertainty again denotes the statistical 1σ spread about the mean value. Fig. 7 illustrates the distribution of the number of unmerged IMBHs across the 100 A-SLOTH simulations.

The systematic uncertainties on the number of unmerged IMBHs are most likely greater than the variability resulting from 100 simulation runs. As the population size is sensitive to the underlying SFR, it is the primary source of systematic uncertainty in this work. Observations with gravitational wave detectors provide the most stringent constraints on the SFR yet. By reproducing high-mass BBH mergers using Pop III remnants, the result from this work is to be understood as an upper bound. Any SFR that yields a population size below the one presented here cannot be excluded.

3.2 Pop III contribution to extragalactic background light

In view of an IMBH-forming population of Pop III stars constrained by observations of BBH mergers, the extent to which their light contributes to the EBL was investigated. Their contribution is determined following the equation for the specific emissivity ε_ν presented by Kneiske et al. (2002)

Fig. 7 Distribution of the number of unmerged IMBHs across the 100 A-SLOTH simulations for the SWT SFR. Light blue lines indicate the standard deviation for each bin.

$${ϵ}_{\nu }(z)={\int }_z^{z_{max}}{L}_{\nu }(t(z)-t\left(z^{\mathrm{\prime }}\right))\psi \left(z^{\mathrm{\prime }}\right)\left|\frac{\mathrm{d}{t}^{\mathrm{\prime }}}{\mathrm{d}{z}^{\mathrm{\prime }}}\right|\mathrm{d}{z}^{\mathrm{\prime }}$$(19)

and the resulting spectral energy distribution of the EBL

$$P_{\nu }(z)=\nu I_{\nu }(z)=\nu \frac{c}{4\pi }{\int }_z^{z_{max}}{ϵ}_{{\nu }^{\mathrm{\prime }}}\left(z^{\mathrm{\prime }}\right)\left|\frac{\mathrm{d}{t}^{\mathrm{\prime }}}{\mathrm{d}{z}^{\mathrm{\prime }}}\right|\mathrm{d}{z}^{\mathrm{\prime }}.$$(20)

The quantities L_ν(t(z)) and ψ(z) are the Pop III luminosity and the SFR, respectively. The cosmological parameters enter via

$$\left|\frac{dt}{dz}\right|=\frac{1}{H_0(1+z)E(z)},$$(21)

with

$$E(z)=\sqrt{{\mathrm{\Omega }}_{\mathrm{\Lambda }}+{\mathrm{\Omega }}_k(1+z{)}^2+{\mathrm{\Omega }}_m(1+z{)}^3+{\mathrm{\Omega }}_r(1+z{)}^4}$$(22)

where a flat ΛCDM cosmology with parameters from the Planck Collaboration XIII (2016) was assumed. The intrinsic Pop III spectra and the reprocessed spectra from the ISM were modeled using the fitting functions presented in Fernandez & Komatsu (2006). The resulting Pop III contribution to the EBL is shown in Fig. 8 for the A-SLOTH, Jaacks, and SWT SFRs. To assess the Pop III contribution to the EBL for the full star formation period, the Jaacks SFR has been extended for z < 7. It has been extrapolated following Eq. (6) to meet the redshift of the Lynx Arc at z = 3.357, which is assumed to be the time of termination of Pop III star formation. Comparing the solid lines with the EBL measurements by JWST (Windhorst et al. 2022) shows that the Pop III contribution to the EBL is, at best, at the 0.1% level. In the case of the SWT SFR it is safe to assume that the contribution is negligible. Although this would lead to overproduction of BBH mergers, scaling the Pop III SFRs such that they match the furthermost available data of the measured MF17 SFR at redshift z = 10 raises the Pop III EBL contributions to the 1% level. The scalings are 20 for A-SLOTH, 4 for Jaacks, and 250 for SWT. While the calculation of the EBL from only Milky Way-like galaxies, as would be the case for the A-SLOTH SFR, would certainly introduce a bias, the above statement likely holds true for a generalized universe, since the numerical simulations of the Jaacks and SWT SFRs did not constrain the size and shape of the formed galaxies. It can be concluded that Pop III stars are, in effect, hidden in both the amount of star formation as well as the EBL. The number of unmerged IMBHs from Pop III stars determined in this work is therefore not in conflict with measured stellar observables. This result is also in agreement with Finke et al. (2022), who have presented a model of the optical and infrared bands of the EBL. While no explicit comments were made on the impact of Pop III stars, the contributions of light from different cosmic times to the present-day EBL suggest that the most significant impact comes from light emitted in the local universe. In fact, the light from redshift z ≥ 3 does not seem to play a role at all. A similar conclusion was reached by Sun et al. (2021) regarding the mean near-infrared radiation background (NIRB) when employing a Pop III SFR model similar to the scaled Jaacks SFR in this work. Their modeling of the angular EBL fluctuations provides promising prospects for discerning the Pop III contribution, as the angular fluctuations show strong wavelength dependence for high multipoles due to prominent Lyman-α emission, which further increases with growing stellar mass. In the scenario suggested here, the expected EBL contribution is predominantly limited by the constraints on the SFR obtained through gravitational wave measurements of the merger rate. A detection of a signature of Pop III stars in the EBL appears to be unlikely given these constraints.

Fig. 8 EBLs for configurations 1-A-0, 1-J-0, and 1-SWT-0 both for unscaled SFRs (solid lines) and scaled SFRs (dashed lines) with factors 20, 4, and 250 for A-SLOTH, Jaacks, and SWT respectively. For comparison, the most recent model of EBL lower limits by Finke et al. (2022) and JWST galaxy count and direct measurements (Windhorst et al. 2022).

4 Discussion

4.1 Caveats

Given the physical complexity of BBH mergers, as they emerge from Pop III stars, the method presented in this work provides a simplified approach to estimate the fraction of BBH merged systems. An effort was made to account for initial pairings that would not result in BBH mergers; however, the assessment of individual binary systems was not possible. Rather, combinations were categorically included or excluded using robust arguments based on our current understanding of BBH formation. With all of this in mind, the uncertainty of the number of unmerged IMBHs with this method is bound to be considerably larger than the statistical 1σ spread among the 100 Milky Way-like galaxy simulations.

With these shortcomings in mind, the question arises as to how reasonable the calculated number of $N_U^{SWT}$ = 131.87 ± 6.09 unmerged IMBHs in a Milky Way-like galaxy is. By constructing binary progenitor Pop III systems with a binary fraction of f_B = 0.5, the same one used by Tanikawa et al. (2022) in their simulations, and making the naive assumption that all binary remnant systems are guaranteed to merge irrespective of stellar evolution and orbital dynamics, it is possible to determine a lower limit for the number of unmerged IMBHs. This number corresponds to the population size of solitary high-mass Pop III stars, which, following a binary fraction of f_B = 0.5, is 1/3 of all high-mass objects. Therefore, the lower limit is $N_{U,lower}^{SWT}$ = 74.21 ± 3.43. The number of unmerged IMBHs calculated with the method presented in this work thus neatly falls between the total number of high-mass objects and the lower limit of unmerged objects.

4.2 Other IMBH formation scenarios

Unmerged IMBHs from Pop III stars present only a subset of objects that might be able to acquire a DM spike. The presented method, for example, does not account for giant stars that form following a common envelope stage of a binary stellar system. These do not create BBH merger signatures as the remnant of a giant star is just a solitary IMBH, even though the initial system was a stellar binary. This IMBH might then go on to acquire a DM spike. On the other hand, a second-generation IMBH that is produced in a merger might also dress with a DM spike if the merger happens sufficiently early to allow the DM spike to grow. Though likely, these mergers would not be picked up by gravitational wave detectors as they would happen too far away. In this scenario, dressed IMBHs can also be produced from the merger of stellar BHs. Additionally, primordial black holes may contribute to the distribution of IMBHs (Blinnikov et al. 2016; Postnov & Chekh 2024).

4.3 Comparison with other works

Compared with previous works conducted on determining the number of unmerged IMBHs from Pop III stars, this project yields comparatively few objects. Bertone et al. (2005) have used the method of populating simulated halos at redshift z = 18 with 100 M_⊙ BHs if the smoothed primordial density field featured a 3σ peak, resulting in a number of 1027 ± 84 unmerged IMBHs. While stellar observables such as the SFR and the EBL would allow for a greater population of Pop III IMBHs, advances in gravitational wave astronomy put very strict constraints on the amount of allowed star formation and, subsequently, the IMBH population size.

Complementary to the presented work is the assessment of “heavy seed” IMBH populations that form from the direct collapse of giant pristine gas clouds. A recent examination of the matter was conducted by Aschersleben et al. (2024), where the results of the large-scale hydrodynamical simulation EAGLE (Schaye et al. 2014) were analyzed to create a mock catalog of unmerged direct collapse IMBHs in Milky Way-like galaxies. The seed BHs are assumed to form with a mass of 10⁵ M_⊙ in host galaxies that exceed a mass of 10¹⁰ M_⊙/h. This study reveals that, on average, a Milky Way-like galaxy and its respective satellites host ${15}_{-6}^{+9}$ heavy seed unmerged IMBHs, an order of magnitude less than the light seed IMBHs that form as Pop III remnants. A comparison with the 101 ± 22 heavy seeds determined by Bertone et al. (2005), using a similar method as discussed above, once again shows that previous analyzes had overestimated the number of IMBHs and that the latest experimental observations and advances in numerical simulations present a more reserved estimate. Even though heavy seeds appear to be 10 times less common than the light seed Pop III remnants, heavy seeds are expected to have far bigger DM spikes, potentially being sufficiently luminous to be spotted by current and future gamma-ray observatories, making them ideal for “smoking gun” detections of DM. While present-day gravitational wave detectors have limited capabilities in making observations of IMBH mergers, the next-generation ground- and space-based detectors will open up the intermediate mass range such that more stringent limits on both light seed as well as heavy seed populations and their formation mechanisms will be set (Fragione & Loeb 2023).

Acknowledgements

We thank Tilman Hartwig and Simon Glover for their assistance in operating the A-SLOTH simulation code. We also thank Justin Finke for providing information on methods and results regarding the contribution of early-type stars to the EBL. Many thanks are expressed towards Sara Porras Bedmar, whose help in performing EBL calculations is greatly appreciated. Special thanks are given to Ataru Tanikawa for extensive discussions on stellar evolutionary models and the merging behavior of stellar remnants which make up the foundation of this work. DH acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306.

References

All Tables

All Figures

	Fig. 1 SFRs from 100 A-SLOTH simulations (blue and green lines). The averaged SFR is shown in white with the overlaid parametrization in red.
In the text

	Fig. 2 Comparison of Pop III A-SLOTH, Jaacks, and SWT SFRs with MF17. A-SLOTH Pop III star formation is terminated at z = 3.357.
In the text

	Fig. 3 Remnant mass function of metal-free stars following Fryer et al. (2012). The mass range above 260 M⊙ is not shown, though it simply follows the dotted line (BH retains full progenitor mass). For a comparison of the stellar model with inefficient convective overshoot (this plot) with a model with efficient convective overshoot, see Fig. 2 in Tanikawa et al. (2021).
In the text

	Fig. 4 Averaged mass distribution of Pop III stars from 100 A-SLOTH simulations. The sampling reflects the log-flat IMF. Coloring indicates different stellar evolutionary paths depending on the ZAMS mass.
In the text

	Fig. 5 Total dead mass density of all IMBH remnants (green) and dead mass densities for the individual IMBH capable bins (blue lines).
In the text

	Fig. 6 Combinations of binary stellar systems shown for five mass bins. Combinations inside dotted outline simulated in Tanikawa et al. (2022). Combinations inside thin solid outline lead to BBH mergers in Tanikawa et al. (2022). Combinations inside thick solid outline lead to mergers of IMBHs (this work).
In the text

	Fig. 7 Distribution of the number of unmerged IMBHs across the 100 A-SLOTH simulations for the SWT SFR. Light blue lines indicate the standard deviation for each bin.
In the text

	Fig. 8 EBLs for configurations 1-A-0, 1-J-0, and 1-SWT-0 both for unscaled SFRs (solid lines) and scaled SFRs (dashed lines) with factors 20, 4, and 250 for A-SLOTH, Jaacks, and SWT respectively. For comparison, the most recent model of EBL lower limits by Finke et al. (2022) and JWST galaxy count and direct measurements (Windhorst et al. 2022).
In the text