© The Authors 2024

1. Introduction

The existence of primordial black holes (PBHs), a hypothetical type of black hole (BH) that formed soon after the Big Bang, could shed light on the mystery surrounding the nature of dark matter (DM), which has persisted for almost a century (Zwicky 1933). DM appears to be the most abundant mass component in the Universe, dominating the dynamics of collapsed objects and providing the underlying structure for the entire visible cosmos. To date, it is widely assumed that DM may consist of so far unknown particles that predominantly interact via gravity and possibly through weak interactions with visible matter (e.g., Feng 2010). However, despite numerous efforts to uncover weakly interacting massive particles as potential DM candidates, no particle has been found to match the predicted properties, such as the allowed range in interaction cross section and energy parameter space (Baudis 2012; Boveia & Doglioni 2018).

The detection of gravitational waves (GWs) from merging black hole binaries by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and the Virgo interferometer (VIRGO) (Abbott et al. 2016) has rekindled the idea, initially proposed by Hawking (1971), that PBHs could serve as DM candidates. It was shown that the merger rates of PBHs could be compatible with this observation without violating the requirement that PBHs are at most as abundant as the DM in the Universe (Bird et al. 2016; Clesse & García-Bellido 2017; Sasaki et al. 2018). Notably, the best-fit power law + peak model determined a Gaussian peak at 33 M_⊙ for black hole masses detected by LIGO (Abbott et al. 2021).

The widely studied hypothesis suggests that PBHs originate from high-density perturbations generated by inflationary dynamics at small scales (Green & Kavanagh 2021; Carr et al. 2021). Depending on the specific formation mechanism, PBHs can emerge with different masses. The nature of the fluctuation enhancement defines the shape of their mass distribution function (Sureda et al. 2021). Peaks in the power spectrum yield nearly monochromatic distributions (Magaña et al. 2023). For instance, scenarios such as chaotic new inflation may generate relatively narrow peaks (Saito et al. 2008). Conversely, models featuring inflection points in the potential plateau (García-Bellido & Ruiz Morales 2017) or hybrid inflation (Clesse & García-Bellido 2015) lead to extended mass functions spanning a wide range of PBH masses.

The potential discovery of PBHs carries significant implications for our understanding of the Universe, even though they only constitute a fraction of the total DM. Notably, PBHs may influence the formation of large-scale structures through Poisson fluctuations (Meszaros 1975; Padilla et al. 2021) and serve as seeds for supermassive black holes (SMBHs) with masses ranging from 10⁵ M_⊙ to 10¹⁰ M_⊙ (Carr & Rees 1984). SMBHs have been identified at the centers of galaxies at redshifts as high as z ∼ 7, corresponding to a Universe age younger than 0.8 billion years (Carr & Silk 2018). This challenges explanations that solely rely on accreting stellar remnant black holes (Volonteri 2010).

This scenario gains plausibility when the seed for these SMBHs is a sufficiently massive PBH that grows through accretion. Moreover, PBHs can potentially serve as seeds for the formation of massive galaxies, providing a potential explanation for the recent detection of unusually massive galaxy candidates at z > 8 by the James Webb Space Telescope (JWST) (Liu & Bromm 2022; Labbé et al. 2023; Glazebrook et al. 2024).

Additionally, PBHs could offer a solution to the cusp-core problem observed in low-mass galaxies without relying on baryonic processes alone (Boldrini et al. 2020). Furthermore, they provide a valuable avenue for gaining insights into the physics governing small scales and early epochs (Cole & Byrnes 2018; Sato-Polito et al. 2019; Kalaja et al. 2019; Gow et al. 2021). The varied masses and densities of PBHs could lead to a range of potentially observable effects.

One particularly intriguing case supporting the cosmological significance of PBHs is SDSS J1004+4112. This system features strong microlensing in a quadruple system with a massive cluster lens, comprising a centrally dominant galaxy surrounded by several smaller galaxies. In this scenario, the quasar images exhibit a remarkable separation of 16 arcsec, which is well beyond the reach of detectable starlight. The observed strong microlensing effect in this system challenges conventional explanations and prompts consideration of a cosmological distribution of solar-mass compact bodies, argued to be PBHs (Hawkins 2020).

As PBHs have not yet been detected, there are constraints on their maximum number density and, consequently, on the fraction of DM they comprise (Carr 2019; Sureda et al. 2021). Specifically, PBHs with masses exceeding M_PBH ∼ 5 × 10⁻¹⁹ M_⊙ could accrete baryonic matter. Additionally, their lifetimes might surpass the age of the Universe (Page & Hawking 1976; MacGibbon 1987; Araya et al. 2021).

Recently, Liu et al. (2022) investigated the impact of PBHs on the formation of the first stars in the Universe. The authors developed a subgrid model integrated into a customized version of GIZMO (Hopkins 2015) and conducted zoom-in simulations of small (∼50 kpc) overdense regions where first star formation is expected to occur at z ∼ 20 − 30. These simulations used high resolutions m_DM ∼ 2 M_⊙ and m_gas ∼ 0.4 M_⊙ for DM and gas particles, respectively. Focusing their study on PBHs with masses around 30 M_⊙, they concluded that the conventional understanding of first star formation in molecular-cooling minihalos remains largely unaltered when PBHs constitute f_PBH ≈ 10⁻⁴ − 0.1 of the DM. However, they observed effects on the clumpiness and density profile of the gas distribution due to the presence of PBHs. These effects arise from the potential of PBHs to heat the gas and/or decelerate its collapse by accretion feedback and substructures around PBHs, ultimately resulting in a shallower gas density profile.

Moreover, Ziparo et al. (2022) explored cosmic radiation backgrounds from PBHs, presenting a semianalytical model that predicts X-ray and radio emissions by these sources. Their model considers the distribution of PBHs within DM halos and in the intergalactic medium (IGM). Notably, emission from PBHs in the IGM is not as significant as the overall PBH emission. It contributes between 1% and 40% to the total. The majority of cosmic X-ray and radio background emission comes from PBHs within DM mini-halos with masses below 10⁶ M_⊙ during early epochs (z > 6). While PBHs (comprising approximately 0.3% of the DM) could explain the observed cosmic X-ray background excess, they cannot explain the cosmic radio background (see also Zhang et al. 2024).

Furthermore, Kashlinsky (2021) employed a hydrodynamical formalism to analyze the interaction between DM and baryons in the early Universe. He found that supersonic advection, particularly in scenarios involving LIGO-type BHs as DM, leads to early collapse, potentially explaining the SMBHs in high-redshift (z > 7) quasars.

Our primary objective is to investigate whether PBHs are plausible constituents of DM and the feedback effect they might produce in the interstellar medium (ISM) and IGM. We accomplish this by examining their potential mass ranges and evaluating their contribution to the overall content of DM. Specifically, we explore monochromatic PBH mass scenarios of 1, 33, and 100 M_⊙, where PBHs might play a significant role in the DM composition, consistent with findings in previous research (Villanueva-Domingo et al. 2021). Furthermore, we explore various levels of PBH contributions to DM, using the parameter f_PBH = 10⁻⁴, 10⁻³, 10⁻², 10⁻¹ and 1 to represent the mass fraction of DM that is composed of PBHs. To assess their impact on the ISM and IGM, we consider PBH accretion feedback, a process involving the injection of energy that alters the properties of the surrounding gas. To achieve this goal, we introduce a new semianalytical model integrated into a hydrodynamical simulation. For this purpose, we use galaxies from the ChemodynamIcal propertiEs of gaLaxies and the cOsmic web (CIELO; Tissera et al. in prep.). These simulations were used previously to study the impact of a local group environment on infalling disk satellites (Rodríguez et al. 2022) and the impact of baryons on the shape of dark matter halos (Cataldi et al. 2023).

To enhance the versatility and applicability of our model while also addressing resolution limitations, we introduce a novel treatment inspired by the work of Liu et al. (2022). This approach involves adapting and integrating their high-resolution halo gas-density profile into our simulation, thereby enabling the study of the accretion regimes of PBHs within the central densest regions of the halos. This method not only enables a comprehensive evaluation of the impact of PBHs on the ISM, but also opens the possibility of its application to larger simulated volumes.

This paper is structured as follows. Section 2 outlines the theoretical framework for PBH accretion physics. Section 3 details the implementation of the semianalytical model in the zoom-in hydrodynamical simulation, and it specifies our methods for gas-particle heating and cooling. Section 4 presents our findings and discusses their implications, including analyses of variations within our model. Finally, Section 5 summarizes our results and conclusions.

2. Primordial black hole accretion physics

In this section, we describe the theoretical model we employed to compute PBH accretion and emission.

2.1. Bondi-Hoyle-Lyttleton accretion

We assumed that in the IGM and in halos, PBHs accrete via the Bondi-Hoyle-Lyttleton model with an accretion rate given by

$$\begin{array}{c}\hfill \dot{M}=4\pi r_{\mathrm{B}}^2\stackrel{\sim }{v}\rho =\frac{4\pi G^2{M}^{2}n\mu m_{\mathrm{p}}}{{\stackrel{\sim }{v}}^3},\end{array}$$(1)

where $r_{\mathrm{B}}=GM/{\stackrel{\sim }{v}}^2$ is the Bondi radius, G = 6.6743 × 10⁻⁸cm³/g/s² is the gravitational constant, M is the PBH mass, n is the gas number density, μ is the mean molecular weight, m_p = 1.6726219 × 10⁻²⁴g is the proton mass, and $\stackrel{\sim }{v}\equiv {(v^2+c_{\mathrm{s}}^2)}^{1/2}$ represents the characteristic velocity. Here, v stands for the relative velocity between the gas and the PBH, which in the simulation can be taken as the relative velocity between the gas and the DM, and c_s stands for the sound speed in the medium, defined for an ideal gas as

$$\begin{array}{c}\hfill c_{\mathrm{s}}=\sqrt{\frac{\gamma k_{\mathrm{B}}T}{\mu m_{\mathrm{p}}}},\end{array}$$(2)

where γ = 5/3 is the adiabatic index, k_B = 1.380649 × 10⁻¹⁶ erg/K is the Boltzmann constant, and T is the absolute gas temperature.

2.2. Condition for accretion disk formation

To determine whether an accretion disk forms around a black hole (BH) with mass M, we analyzed whether the outer edge of the BH accretion disk, denoted as r_out, resides within the radius of the innermost stable circular orbit (ISCO). No efficient disk formation is expected when the outer edge is situated within the ISCO radius. The radius of the outer disk edge is given by Agol & Kamionkowski (2002)

$$\begin{array}{c}\hfill r_{\mathrm{out}}\simeq 5.4\times {10}^9{r}_{\mathrm{s}}{\left(\frac{M}{100{\mathrm{M}}_{\odot }}\right)}^{\frac{2}{3}}{\left(\frac{\stackrel{\sim }{v}}{10\mathrm{km}/\mathrm{s}}\right)}^{-\frac{10}{3}},\end{array}$$(3)

where r_s = 2GM/c² represents the Schwarzschild radius of the BH, and c = 2.9979246 × 10¹⁰cm/s is the speed of light in vacuum.

For a Schwarzschild BH, that is, a BH without rotation that is spherically symmetric and uncharged, the ISCO radius is defined as

$$\begin{array}{c}\hfill r_{\text{ISCO}}=3r_{\mathrm{s}}.\end{array}$$(4)

This radius is the minimum distance at which an object can maintain a stable circular orbit around the BH, resisting the overpowering gravitational pull.

2.3. Accretion regimes

According to the accretion rate, the flow of material undergoing accretion could lead to the formation of a thin disk or of an advection-dominated accretion flow (ADAF), characterized by an uncertain outer accretion disk radius. Following Takhistov et al. (2022), we assumed that when the dimensionless accretion parameter ṁ = Ṁ/Ṁ_Edd exceeds a threshold of 0.07α, where α = 0.1 represents the viscosity parameter, a standard thin α disk forms. Here, Ṁ_Edd is the Eddington accretion rate, calculated assuming Eddington luminosity and a radiative efficiency of 0.1, characteristic of thin disks (Kato et al. 2008).

Furthermore, for a more comprehensive analysis, following Takhistov et al. (2022), we divided the ADAF into three distinct subregimes. The first subregime is the luminous hot accretion flow (LHAF), where the accretion parameter falls within the range Ṁ_ADAF = 0.1α² < ṁ < Ṁ_LHAF = 0.7α. The second regime corresponds to the standard ADAF, applicable when Ṁ_ADAF = 10⁻³α² < ṁ < Ṁ_ADAF. Last, the electron ADAF (eADAF) regime applies when ṁ < Ṁ_eADAF. In the following sections, we study the details of each accretion regime.

2.3.1. Thin disk

The standard thin disk formed by the accretion onto PBHs is optically thick and efficiently emits blackbody radiation (Shakura & Sunyaev 1973), enabling a fully analytical description. The variation in disk temperature with radius outside the inner disk region is described by (Pringle 1981)

$$\begin{array}{c}\hfill T(r)=T_{\mathrm{i}}{\left(\frac{r_{\mathrm{i}}}{r}\right)}^{3/4}{[1-{\left(\frac{r_{\mathrm{i}}}{r}\right)}^{\frac{1}{2}}]}^{\frac{1}{4}},\end{array}$$(5)

where r_i is the inner disk radius, which is taken to be the ISCO radius, and

$$\begin{array}{c}\hfill T_{\mathrm{i}}=\left(\frac{3GM\dot{M}}{8\pi r_{\mathrm{i}}^3{\sigma }_{\mathrm{B}}}\right)=53.3\mathrm{eV}{\left(\frac{n}{1{\mathrm{cm}}^{-3}}\right)}^{1/4}{\left(\frac{\stackrel{\sim }{v}}{10\mathrm{km}/\mathrm{s}}\right)}^{-3/4},\end{array}$$(6)

where σ_B = 5.6703744 × 10⁻⁵g/s³/K⁴ is the Stefan-Boltzmann constant.

Beyond the inner disk edge, the disk reaches a maximum temperature of T_max = 0.488T_i at radius r = 1.36r_i.

The thin-disk emission spectrum results from a combined contribution of the blackbody spectra from each radius. Using the scaling relations of Pringle (1981) and requiring continuity, the resulting spectrum can be approximately described as

$$\begin{array}{c}\hfill \begin{array}{cc}& \nu <T_{\mathrm{o}}:L_{\nu }=c_{\alpha }{\left(\frac{T_{max}}{T_{\mathrm{o}}}\right)}^{\frac{5}{3}}{\left(\frac{\nu }{T_{max}}\right)}^2\hfill \\ \hfill & T_{\mathrm{o}}<\nu <T_{max}:L_{\nu }=c_{\alpha }{\left(\frac{\nu }{T_{max}}\right)}^{\frac{1}{3}}\hfill \\ \hfill & T_{max}<\nu :L_{\nu }=c_{\alpha }{\left(\frac{\nu }{T_{max}}\right)}^2{e}^{1-\frac{\nu }{T_{max}}},\hfill \end{array}\end{array}$$(7)

where ν is the photon energy, T_o = T(r_out) is the temperature at the outer disk radius, and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill c_{\alpha }=& 1.27\times {10}^{29}\mathrm{erg}{\mathrm{eV}}^{-1}{\mathrm{s}}^{-1}{\left(\frac{M}{{\mathrm{M}}_{\odot }}\right)}^2\hfill \\ \hfill & {\left(\frac{n}{1{\mathrm{cm}}^{-3}}\right)}^{3/4}{\left(\frac{\stackrel{\sim }{v}}{10\mathrm{km}/\mathrm{s}}\right)}^{-9/4}.\hfill \end{array}\end{array}$$(8)

c_α is normalized such that the emission has the maximum possible efficiency for a Schwarzschild BH (∫L_νdν = 0.57 Ṁc²) (Takhistov et al. 2022).

2.3.2. ADAF

When an ADAF disk forms, the heat generated by viscosity is not efficiently radiated away, and a significant amount of energy is advected into the BH event horizon by matter heat capture and the gas influx. In this scenario, the disk exhibits a complex multicomponent emission spectrum, with the dominant components stemming from electron cooling, synchrotron radiation, inverse-Compton (IC) scattering, and bremsstrahlung processes.

To characterize the flow, we adopted the following parameter values: the fraction of the viscously dissipated energy that directly heats electrons, δ = 0.3; the ratio of gas pressure to total pressure, β = 10/11; the minimum flow radius, set equal to the ISCO radius, r_min = 3r_s; and the viscosity parameter, α = 0.1. These parameter choices are consistent with recent numerical simulations and observations (Yuan & Narayan 2014), and they were also employed in the model presented by Takhistov et al. (2022).

The synchrotron emission is self-absorbed and peaks at a photon energy of

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\nu }_p=& 1.83\times {10}^{-2}\mathrm{eV}{\left(\frac{\alpha }{0.1}\right)}^{-\frac{1}{2}}{\left(\frac{1-\beta }{1/11}\right)}^{\frac{1}{2}}\hfill \\ \hfill & {\theta }_{\mathrm{e}}^2{\left(\frac{r_{min}}{3r_{\mathrm{s}}}\right)}^{-\frac{5}{4}}{\left(\frac{M}{{\mathrm{M}}_{\odot }}\right)}^{-\frac{1}{2}}{\left(\frac{\dot{m}}{{10}^{-8}}\right)}^{\frac{3}{4}},\hfill \end{array}\end{array}$$(9)

where θ_e represents the temperature in units of the electron mass m_e,

$$\begin{array}{c}\hfill {\theta }_{\mathrm{e}}=\frac{k{T}_{\mathrm{e}}}{m_{\mathrm{e}}{c}^2}=\frac{T_{\mathrm{e}}}{5.93\times {10}^9\mathrm{K}}.\end{array}$$(10)

The peak luminosity is given by

$$\begin{array}{c}\hfill L_{{\nu }_{\mathrm{p}}}=5.06\times {10}^{38}\frac{\mathrm{erg}}{\mathrm{s}·\mathrm{eV}}{\alpha }^{-1}(1-\beta )\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\dot{m}}^{\frac{3}{2}}{\theta }_{\mathrm{e}}^5{\left(\frac{r_{\text{min}}}{r_{\mathrm{s}}}\right)}^{-\frac{1}{2}}.\end{array}$$(11)

The synchrotron photons undergo IC scattering with the surrounding electron plasma. The resulting IC spectrum is given by

$$\begin{array}{c}\hfill L_{\nu ,\mathrm{IC}}=L_{{\nu }_{\mathrm{p}}}{\left(\frac{\nu }{{\nu }_{\mathrm{p}}}\right)}^{-{\alpha }_{\mathrm{c}}},\end{array}$$(12)

where

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{c}}=-\frac{ln{\tau }_{\mathrm{e}\mathrm{s}}}{lnA}\end{array}$$(13)

with ${\tau }_{\mathrm{e}\mathrm{s}}=12.4\dot{m}{\alpha }^{-1}{\left(\frac{r_{\text{min}}}{r_{\mathrm{s}}}\right)}^{-1/2}$ the electron scattering optical depth, and $A=1+4{\theta }_{\text{e}}+16{\theta }_{\text{e}}^2$.

Finally, the emission from bremsstrahlung of the thermal spectrum is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L_{\nu ,\text{brems}}=& 1.83\times {10}^{17}{\left(\frac{\alpha }{0.1}\right)}^{-2}{\left(\frac{c_1}{0.5}\right)}^{-2}ln\left(\frac{r_{max}}{r_{min}}\right)F\left({\theta }_{\mathrm{e}}\right)\hfill \\ \hfill & {\left(\frac{T_{\mathrm{e}}}{5\times {10}^9\mathrm{K}}\right)}^{-1}{e}^{-(h\nu /k{T}_{\mathrm{e}})}\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\dot{m}}^2\mathrm{ergs}{\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1},\hfill \end{array}\end{array}$$(14)

where c₁ = 0.5 and F(θ_e) is given by Mahadevan (1997)

$$\begin{array}{c}\hfill F({\theta }_e)=\{\begin{array}{cc}\begin{array}{cc}& 4{\left(\frac{2{\theta }_e}{{\pi }^3}\right)}^{\frac{1}{2}}(1+1.78{\theta }_e^{1.34})\hfill \\ \hfill & +1.73{\theta }_e^{\frac{3}{2}}(1+1.1{\theta }_e+{\theta }_e^2-1.25{\theta }_e^{\frac{5}{2}}),\hfill \end{array}\hfill & {\theta }_e\le 1,\hfill \\ \begin{array}{cc}& \left(\frac{9{\theta }_e}{2\pi }\right)(ln[1.12{\theta }_e+0.48]+1.5)\hfill \\ \hfill & +2.30{\theta }_e(ln\left[1.12{\theta }_e\right]+1.28),\hfill \end{array}\hfill & {\theta }_e\ge 1.\hfill \end{array}\end{array}$$(15)

Each ADAF subregime has different temperature dependences because θ_e is determined by balancing the heating and radiation processes (Mahadevan 1997). Direct viscous electron heating is dominant at low accretion rates, while ion-electron collisional heating replaces it at higher accretion rates. The method for calculating θ_e in each subregime can be found in Section 3.1.3 of Takhistov et al. (2022).

Figure 1 illustrates the boundary lines that define the different accretion regimes using a fixed PBH mass of 33 M_⊙. The interplay between the gas density and the characteristic velocity of the surrounding medium determines to which regime the system belongs. The figure illustrates a clear correlation, indicating that lower characteristic velocities and higher gas densities are associated with higher accretion rates, delineating the boundaries between different accretion regimes.

Fig. 1. Continuous color map of the parameter ṁ across a grid of characteristic velocities and hydrogen number densities considering a fixed PBH mass of 33 M⊙. The dashed black lines indicate the thresholds for MeADAF, MADAF, and MLHAF as a function of the viscosity parameter α = 0.1. The region above MLHAF corresponds to the thin-disk formation regime, while in the region between MLHAF and MADAF, LHAF occurs. The area between MADAF and MeADAF corresponds to the standard ADAF. The eADAF regime exists below MeADAF.

Figure 2 displays examples of the resultant spectra (calculated using the formalism described above) generated by PBHs within varying gas density environments and having different masses, all at a fixed characteristic velocity in the medium of 10 km/s. The curves showing a knee are indicative of thin-disk emission. With higher densities at a certain PBH mass, the power produced by PBHs increases.

Fig. 2. Spectra generated by accreting PBHs of 1, 33, and 100 M⊙ (in yellow, red, and purple, respectively) in different gas density environments (1 cm−3 as dotted lines, 102 cm−3 as dash-dotted lines, and 104 cm−3 as solid lines). This example assumes a constant characteristic velocity $\stackrel{\sim }{v}=10$ km/s and a mean molecular weight μ = 1.22.

3. Implementation

In this section, we provide an overview of our semianalytical model implementation. We start by introducing the hydrodynamical simulation we used and then describe the essential subgrid physics. Given our focus on implementing PBHs on galactic scales, where the dynamical and temporal timescales vary widely, we required simplifications and adaptations to mimic the expected effects at the scales we can resolve.

3.1. CIELO simulations

We worked with the high-resolution run of the CIELO suite of zoom-in hydrodynamical simulations. This project aimed to study the assembly and evolution of galaxies in different environments (Tissera in prep.; Rodríguez et al. 2022; Cataldi et al. 2023).

These simulations adopted a ΛCDM cosmology with Ω₀ = 0.317, Ω_Λ = 0.6825, Ω_B = 0.049, and h = 0.6711 (Planck Collaboration I 2014).

The version of GADGET-3 we used to run these simulations includes a multiphase model for the gas component, metal-dependent cooling (Sutherland & Dopita 1993), star formation, and supernova (SN) feedback as described in Scannapieco et al. (2005) and Scannapieco et al. (2006). These multiphase and SN-feedback models were used to successfully reproduce the star formation activity of galaxies during quiescent and starburst phases and can drive violent mass-loaded galactic winds with a strength reflecting the depth of the potential well (Scannapieco et al. 2005, 2006). This physically motivated thermal SN-feedback scheme does not include any ad hoc mass-scale-dependent parameters. As a consequence, it is particularly well suited for the study of galaxy formation in a cosmological context.

For this work, we used the highest-resolution run of the Pehuen halos¹ of the CIELO suite, which corresponds to 29 zoom-in initial conditions generated with the MUSIC code (Hahn & Abel 2011). Pehuen halos were taken from a DM-only run of a cosmological periodic cubic box with a side length of L = 50 Mpc h⁻¹. In the highest-resolution run, the DM mass resolution was 1.36 × 10⁵ M_⊙ h⁻¹, and the initial gas particle mass was 2.1 × 10⁴ M_⊙ h⁻¹.

Virialized halos were identified by applying the friends-of-friends (FoF) algorithm (Davis et al. 1985), and the substructures within each dark matter halo were selected by using the SUBFIND code (Springel et al. 2001; Dolag et al. 2009). Figure 3 illustrates the projected distribution of hydrogen number density and temperature in the simulation at z ∼ 23. The location of the halos considered for this analysis is highlighted. These halos have more than 15 and 20 gas an dark matter particles, respectively.

Fig. 3. Projected distribution of the hydrogen number density (left panel) and temperature (right panel) of gas in the Pehuen halos at z ∼ 23. Each halo considered in this paper is distinguished by a labeled circle with a radius equal to ten times the virial radius of the respective halo.

Considering the nature of this work, it is relevant to mention that CIELO does not have a molecular cooling implementation, which is a relevant factor at the redshift under consideration. However, we applied this cooling semianalytically in post-processing to correctly follow the changes in the gas temperature due to PBHs at very high redshift. All physical processes were coded so that they could be switched off or on, as required.

The following sections provide a comprehensive description of the gas heating and cooling model (Sections 3.2 and 3.3, respectively), along with the application of these models to the simulation (Section 3.4). In this implementation, we considered each DM particle to consist of a specific constant fraction of PBHs (f_PBH), and we assumed that the phase-space distribution of PBHs at the scales resolved by CIELO simulations is identical to that of standard cold DM. To select the gas particles affected by PBH feedback, for a DM particle i, we used the smoothed-particle hydrodynamics (SPH) smoothing length of the nearest gas particle, denoted as h_i. Importantly, for the purposes of this paper, we chose to demonstrate the model outcomes at z ≈ 23, a time preceding the formation of the first stars in the simulation. This is because we do not expect PBH effects to be significant in comparison to SN feedback. Nonetheless, the versatility of our model extends to any other simulation and redshift.

3.2. Gas heating

We calculated the resulting heating power from the PBH emission as follows:

$$\begin{array}{c}\hfill {\mathcal{H}}_{\mathrm{disc}}(M,n,v)={\displaystyle {\int }_{E_{min}}^{E_{max}}}{L}_{\nu }(M,n,v)f_{\mathrm{h}}(1-e^{-\tau })\mathrm{d}\nu ,\end{array}$$(16)

where L_ν(ν) is the specific luminosity at frequency ν, calculated using the formalism presented in Section 2.3. Here, f_h is the fraction of energy deposited as heat. Following Takhistov et al. (2022), we adopted f_h = 1/3, which is consistent with the estimates of Shull & van Steenberg (1985), Ricotti et al. (2002), Furlanetto & Stoever (2010), and Kim (2021). Additionally, τ is the optical depth for the gas component. The integral was computed over the energy range from E_min = 13.6 eV to E_max = 5T_i, with T_i as the characteristic temperature of the inner disk. For energies below E_min, hydrogen gas behaves as optically thin to continuous emission, leading to inefficient ionizing radiation absorption. Given the exponential decline of thin-disk emission spectra at higher energies, contributions beyond E_max become negligible.

For 13.6 eV < E < 30 eV, we employed the photoionization cross section given by (Bethe & Salpeter 1957; Band et al. 1990),

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \sigma (E)={\sigma }_0{y}^{-\frac{3}{2}}{(1+y^{\frac{1}{2}})}^{-4},\end{array}\end{array}$$(17)

where y = E/E₀, E₀ = 1/2E_i, and σ₀ = 605.73 Mb = 6.06 × 10⁻¹⁶ cm². Above 30 eV, we used attenuation length data from Figure (32.16) of Olive (2014).

The optical depth for a gas system of size l and with number density n is then given by

$$\begin{array}{c}\hfill \tau (n,E)=\sigma (E)nl.\end{array}$$(18)

3.3. Gas cooling

Since PBH feedback heats the surrounding gas, it is essential to correctly account for temperature loss due to cooling processes to determine the effective heating rate. In this section, we outline the cooling mechanisms we considered in our model.

3.3.1. Atomic cooling

We calculated the atomic cooling considering the processes described in Cen (1992), that is, collisional excitation of neutral hydrogen (H⁰) and singly ionized helium (He⁺) by electrons, collisional ionization of H⁰, He⁰, and He⁺, standard recombination of H⁺, He⁺, and He⁺⁺, dielectronic recombination of He⁺, and free-free emission (bremsstrahlung) for all ions.

We adopted the cooling rates obtained by following the approach described in Katz et al. (1996), that is, considering collisional equilibrium abundances at any given temperature: The collisional excitation contributions from H and He dominate the cooling until the free–free transitions become important. Table 1 of Katz et al. (1996) contains the set of radiative cooling rates adopted in this work (the formulas were taken from Black 1981, with the modifications introduced by Cen 1992).

The atomic cooling functions we implemented are valid for gas temperatures higher than 5 × 10³ K. At temperatures outside this range, the formation of negative ions (H⁻) and molecules (H₂, ${\text{H}}_2^{+}$, ${\text{H}}_3^{+}$, and HeH⁺) can affect the thermal properties of the gas (Black 1981).

3.3.2. Molecular cooling

For gas temperatures lower than 5 × 10³ K, we considered molecular cooling by molecular hydrogen (H₂) and hydrogen deuteride (HD). We implemented the H₂ cooling following the model of Galli & Palla (1998). In this model, the final functional form of the total cooling in units of ergcm⁻³ s⁻¹ is given by

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{{\mathrm{H}}_2}=\frac{n_{{\mathrm{H}}_2}{\mathrm{\Lambda }}_{{\mathrm{H}}_2,\mathrm{LTE}}}{1+{\mathrm{\Lambda }}_{{\mathrm{H}}_2,\mathrm{LTE}}/{\mathrm{\Lambda }}_{{\mathrm{H}}_2,\mathrm{n}\to 0}},\end{array}$$(19)

where Λ_H2, LTE= H_R + H_V is the high-density limit, expressed as the sum of the rotational and the vibrational cooling at high densities, given by (Hollenbach & McKee 1979)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_R=& (9.5\times {10}^{-22}{T}_3^{3.76})/(1+0.12T_3^{2.1})\hfill \\ \hfill & \times exp[-{(0.13/T_3)}^3]+3\times {10}^{-24}exp[-0.51/T_3]\hfill \\ \hfill H_V=& 6.7\times {10}^{-19}exp(-5.86/T_3)\hfill \\ \hfill & +1.6\times {10}^{18}exp(-11.7/T_3)\hfill \end{array},\end{array}$$(20)

with T₃ = T/10³, and Λ_{H2, n → 0} is the low-density limit, given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill log\left({\mathrm{\Lambda }}_{{\mathrm{H}}_2,\mathrm{n}\to 0}\right)=& [-103+97.59log(T)-48.05log{(T)}^2\hfill \\ \hfill & \times 10.8log{(T)}^3-0.9032log{(T)}^4]n_{\mathrm{H}}.\hfill \end{array}\end{array}$$(21)

The final H₂ cooling is valid in the range 13 K < T < 10⁵ K.

For HD cooling, we used the modified version of Lipovka et al. (2005) given by Coppola et al. (2011), which is valid for gas temperatures below 2 × 10⁴ K. Lipovka et al. (2005) provided a two-parameter fit of the collisional cooling function that depends on temperature and H fractional abundance. Coppola et al. (2011) evaluated a radiative cooling function that takes into account both dipole and quadrupole electric transitions, for which they used the calculations by Abgrall et al. (1982).

For completeness, we note that in primordial gas, another cooling agent is provided by LiH, but this is found to be always subdominant (Liu & Bromm 2018) and was thus not considered here.

3.3.3. Inverse-Compton cooling

In addition to the radiative cooling processes, we included IC cooling of the microwave background given by Ikeuchi & Ostriker (1986),

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{C}}=5.41\times {10}^{-36}{n}_{\mathrm{e}}(T-T_{\mathrm{CMB},\mathit{z}}){(1+z)}^4{\text{erg s}}^{-1}{\mathrm{cm}}^{-3},\end{array}$$(22)

where T_CMB, z = 2.73 × (1 + z) is the CMB temperature at redshift z.

Finally, we took into account the volume expansion (contraction) of gas particles in each model iteration, resulting in adiabatic cooling (heating).

Given the heating and cooling theoretical framework, we describe below how we implemented them in our model.

3.4. Subgrid model

We worked with per-particle data extracted from the CIELO simulation, encompassing both gas and DM particles. Under the assumption that each DM particle contains a specific constant fraction of PBHs, we calculated the gas heating by PBH accretion over a period of 5 Myr starting at z = 22.85. This short period was chosen because it allowed us to assume that the particles did not undergo significant changes in their positions and were not strongly affected by gravitational effects among themselves. The starting redshift was chosen at the last available snapshot in the CIELO simulation before the onset of star formation. This setup allowed us to explore the possible impact of PBH feedback on gas properties, and hence, on the regulation of the star formation activity. Additionally, this also ensured that there were virialized halos containing at least ten gas particles. Given the substantial sensitivity of cooling rates to gas temperature, we executed gas heating and cooling computations through 135 discrete time steps. The selection of this time step count was determined through a convergence test that determined the minimum number of steps required for convergence while ensuring computational efficiency.

We recalculated the mass accreted by the PBHs at each time step by considering Eq. (1). We find that the mass accreted by PBHs is negligible in all our calculations, ensuring that the masses of the PBHs remains unchanged throughout the implementation.

Below, we provide a detailed breakdown of the iterative procedure steps.

• We computed the density and characteristic velocity of the gas surrounding the PBHs. To accomplish this, we encompassed all gas particles within a sphere centered on a DM particle, with a radius equivalent to 2h_i. The attributes of these gas particles were weighted using the spline kernel such that for a gas particle j (Monaghan & Lattanzio 1985),

  $$\begin{array}{c}\hfill W_j(r,h_j)=\frac{8}{\pi h_j^3}\{\begin{array}{cc}1-6{\left(\frac{r}{h_j}\right)}^2+6{\left(\frac{r}{h_j}\right)}^3,\hfill & 0⩽\frac{r}{h_j}⩽\frac{1}{2}\hfill \\ 2{(1-\frac{r}{h_j})}^3,\hfill & \frac{1}{2}<\frac{r}{h_j}⩽1\hfill \\ 0,\hfill & \frac{r}{h_j}>1\hfill \end{array},\end{array}$$(23)

  where r represents the distance between DM particle i and gas particle j, and h_j is the smoothing length of gas particle j. All properties were estimated by using the kernel to weight the contributions of the neighboring gas particles. These particles are called neighbors of a given DM particle.

• We injected energy contributed by the PBHs associated with a given DM particle into all the gas neighbors of a given DM particle using the weights of Eq. (23).

• After all the heating contributions resulting from accretion onto PBHs were incorporated into the thermal energy of each individual gas particle, we calculated the temperature increase and updated n_H0, n_H+, n_He0, n_He+, n_He++, and n_e number densities using Eqs. (30)–(38) of Katz et al. (1996), but considering n_H = n_H0 + n_H+ + n_H2, with n_H fixed at the value given by the CIELO simulation to conserve the numbers of hydrogen nuclei (see Appendix A for more details). For simplicity, n_H2 was calculated by using the approximation of Tegmark et al. (1997), such that the fraction of H₂ with respect to the total hydrogen number density was given by

  $$\begin{array}{c}\hfill f_{{\mathrm{H}}_2}=3.5\times {10}^{-4}{\left(\frac{T}{{10}^3\mathrm{K}}\right)}^{1.52}.\end{array}$$(24)

  This approximation is valid for gas temperatures T ≲ 10³K, the temperature range in which cooling by H₂ works.

• Then, we cooled the gas depending on its temperature, as described in Sections 3.3.1, 3.3.2, and 3.3.3, and recalculated the per-particle gas volume due to adiabatic expansion or contraction.

3.5. Gas density treatment

Figure 3 shows that CIELO has resolution limitations that prevent us from considering the high densities expected at scales below 10 pc. To overcome this limitation in gas resolution, we adopted an innovative approach. We used the density profile at the onset of star formation, at z ≈ 30, within a halo of 2 × 10⁵ M_⊙ from the case A zoom-in region for LCDM from Liu et al. (2022), which offers better resolution (hereafter referred to as the reference density profile). First, we rescaled the reference density profile so that it aligned with the virial mass of our halos. We achieved this by multiplying the reference profile by the ratio of the enclosed mass within our CIELO halos at their respective virial radii to the enclosed mass given by the reference halo at its virial radius. Although our implementation started at z ∼ 23, the mass of the reference density profile typically was within an order of magnitude of the masses of the halos present in CIELO. This effectively integrated the reference density profile into our simulations, compensating for resolution limitations below about 10 pc. Subsequently, by discerning the percentage of gas with a particular density within the halos, we determined the proportion of accretion inside the halo that occurred through each mechanism (eADAF, ADAF, LHAF, and thin disk), as these mechanisms are contingent upon the density of accreted gas. Figure 4 illustrates an example of the rescaled density profile for a CIELO halo (upper panel) and the enclosed mass profile (lower panel). With this information, we calculated in this example the accretion rate of a 33 M_⊙ PBH using a mechanism determined by the gas density at different radial locations within the halo.

Fig. 4. Gas density profile and enclosed mass profile for a 33 M⊙ PBH. The upper panel displays the gas density profile taken from Liu et al. (2022) (dashed black line) alongside the rescaled density profile for one of our CIELO halos (red line). In the lower panel, the enclosed mass profile over the virial mass for this halo is shown. In both panels, the profiles are color-coded based on the gas-accretion rate of a 33 M⊙ PBH. For this particular example, we set $\stackrel{\sim }{v}=10$ km/s and μ = 1.22, but in our implementation, these values varied depending on per-particle information from CIELO. Each panel is segmented into colored zones that align with different radii and their respective gas densities, so that the corresponding accretion mechanisms (eADAF, ADAF, LHAF, and thin disk) are active within the segments. The profiles are color-coded based on the accretion rate of the PBH.

We employed a subgrid model to characterize the density distribution of gas particles, focusing on the positions of neighboring gas surrounding DM particles. This method aids in identifying regions in which it is suitable to apply the halo density profile. To accomplish this, we employed a criterion that identifies halo-dominated regions when 70% of the neighboring gas particles are gravitationally bound to the halo. Following this, we used the rescaled reference density profile to the regions identified as halo-dominated.

This novel method successfully takes the gas density distributions within dark matter halos into account beyond the spatial scales that can be resolved in a cosmological context.

In Fig. 5, we present a comprehensive view of PBH accretion and heating across various masses (1 M_⊙, 33 M_⊙, and 100 M_⊙ from left to right) within the density versus velocity plane. The contour lines delineate specific accretion regimes (eADAF, ADAF, and LHAF) that are aligned with the scheme in Fig. 1. The figure excludes the 10⁻¹² M_⊙ PBH case because the dominance of the eADAF regime across the entire velocity and density range for this mass leads to significantly lower accretion rates and heating values than for the other PBH masses explored in this study. All depicted results are representative of conditions at z ∼ 23, giving out the initial gas conditions in our implementation, hereafter referred to as the large-scale implementation.

Fig. 5. Accretion and heating rates of PBHs with varying masses in the density vs. velocity plane. The color-coded panels illustrate the behavior of PBHs with varying masses (1 M⊙, 33 M⊙, and 100 M⊙ columns from left to right) in the density vs. velocity plane, depicting the accretion rate (upper row) and the resulting heating rate (bottom row). The dashed white lines delineate the upper limits of specific accretion regimes (eADAF, ADAF, and LHAF) as shown in Fig. 1. The red rectangles enclose the density and characteristic velocity ranges around each PBH within the halos, employing a rescaled version of the gas density profiles from Liu et al. (2022) and the CIELO particle velocities. The blue rectangles encapsulate weighted values around PBHs in the IGM. The dashed yellow line marks the peak density of a gas particle in the CIELO simulation. These results are representative of gas conditions at z ∼ 23.

4. Results and discussion

4.1. Semianalytical model outcomes

This section presents our semianalytical model outcomes applied to the CIELO simulation, based on which, we explored the PBH impact on gas properties within and outside halos. The analysis took place at z ∼ 23 and the model was implemented over 5 Myr. We considered PBH masses of 1, 33, and 100 M_⊙ as potential DM contributors. We also tested varied PBH DM contributions ranging from f_PBH = 10⁻⁴ to 1. Our semianalytical model assessed changes in gas properties (in particular, n_H, n_H2, and T) for each PBH mass and f_PBH combination.

As indicated in Eq. (1), the accretion rate depends not only on the gas density, but also on the mean molecular weight and on the characteristic velocity in the medium. Figure 6 presents a binned spatial distribution of DM particles surrounding a halo with a virial mass of 4.6 × 10⁶ M_⊙ and a virial radius of 0.2 kpc. The color of each bin represents the average value of the relative velocity between the DM particle and the surrounding gas, the sound speed in the medium, and the mean molecular weight (from left to right). These properties distribution leads to different PBH accretion regimes and, consequently, different emission intensities. The characteristic velocity in the medium is primarily determined by the relative velocity between the gas and the DM particles.

Fig. 6. Key gas properties around dark matter particles in projected spatial coordinates. The upper panels illustrate the spatially projected distribution of dark matter particles in the (x, y) coordinates, and the lower panels depict the distribution in the (x, z) coordinates. The colors indicate the mean values of the weighted (from left to right) relative velocity between the DM particle and the neighboring gas particles, the sound speed in the medium, and the mean molecular weight. The particle distribution is centered on a halo with Mvir = 4.6 × 106 M⊙, Tvir = 8768 K, and rvir = 0.2 kpc, highlighted by cyan circles in each panel.

In Fig. 7, the spatial distribution of gas particles is depicted, colored according to properties obtained through the implementation of our model. The left panel illustrates the ratio of the temperature when applying only our cooling model without PBH heating, and the virial temperature of the halo. The middle and right panel present results obtained by applying our reference model (M_PBH, f_PBH) = (33 M_⊙, 1). In the middle panel, we show the ratio of the final temperature with this model and the virial temperature of the halo, and in the right panel, we display the final neutral hydrogen number density.

Fig. 7. As in Fig. 6, but for the spatial distribution of gas particles and the mean values (from left to right) of the ratio of the temperature obtained by applying only the cooling model (fPBH = 0) and the virial temperature of the halo, the ratio of the temperature considering the entire dark matter component as composed of 33 M⊙ PBHs (fPBH = 1) and the virial temperature of the halo (Tvir = 8768 K), and the final neutral gas density of this last model.

The differences with and without PBHs are noticeable. While the gas can cool down in the cooling-only model (i.e., without the PBH model and only the cooling models), in the reference model, the gas experiences significant heating in regions with higher gas densities and lower characteristic velocities. Inside the halo, the temperature increases by roughly one order of magnitude compared to the virial temperature, and the neutral hydrogen gas density decreases considerably, rendering these conditions unsuitable for star formation.

We extended this analysis to six halos, covering a range of virial masses from 1.2 × 10⁶ to 4.7 × 10⁶ M_⊙ and virial temperatures from 3800 K to 8800 K. We applied our implementation to all DM particles around approximately 14 × r_vir, which was chosen to encompass gas particles well within the ISM and IGM.

In the bottom panels of Figs. 8 and 9, the medians of the final gas temperature and neutral hydrogen number density are displayed, respectively. These medians are presented as functions of the distance from the center of the halos. Initially, individual profiles were calculated for each halo, and subsequently, the median values were derived among the six halos. Each column corresponds to the results obtained considering different PBH masses and f_PBH values. The middle and upper panels of Figs. 8 and 9 depict the ratios of the final temperature and neutral hydrogen number density values with and without PBHs, respectively. The upper panel of Fig. 8 shows the final gas temperature of the PBHs model over the virial temperature of the halo. For reference, the case with the cooling-only model implementation (f_PBH = 0) is also included (yellow stars). The error bars represent the standard deviation within each distance bin.

Fig. 8. Median temperature ratios and gas temperature as a function of distance from the halo center. From top to bottom, each row shows the median ratio of the final (i.e., after implementing our model) and the halo virial temperature, the median ratio of the final gas temperature and the gas temperature when only our cooling model is applied, and the median gas temperature. All as a function of distance from the halo center, expressed as a fraction of the virial radius. Each column from left to right corresponds to considering PBHs with masses of 1, 33, and 100 M⊙, with different fractions of DM composed of these PBHs, as indicated in each panel. For reference, we include the case of fPBH = 0, where we applied our cooling model alone. The error bars represent the standard deviation in each distance bin. Individual profiles were first calculated around each halo separately, and then the median values of the temperature ratios were derived from these individual profiles of the six halos. The virial temperatures range from 3800 to 8800 K.

Fig. 9. Same as the middle and lower panels of Fig. 8, but for the neutral hydrogen number density.

However, for higher PBH masses, an increase in PBH contributions to the DM budget, denoted by a higher value of f_PBH, leads to more pronounced alterations in gas properties. Specifically, in the cases of PBH masses of 1 M_⊙, 33 M_⊙, and 100 M_⊙, along with fractions exceeding or equal to 10⁻², 10⁻³, and 10⁻³ respectively, the temperature rises to over 10⁴ K, and the neutral hydrogen density is significantly affected toward the halo center. These changes in gas properties become less prominent at greater distances from the halo due to the properties of the medium (density and characteristic velocity), resulting in lower PBH accretion rates and, consequently, reduced energy injection into the medium. For this reason, our analysis was focused toward the halo center.

For a PBH mass of 1 M_⊙, fractions of 10⁻², 10⁻¹, and 1 yield temperatures of (1.7, 2.0, and 2.0)T_vir, respectively, with these values representing the median within the halos. In terms of the temperature of the model without PBHs and employing our cooling model (T_fPBH = 0) alone, these fractions result in temperatures toward the halo center of (324.8, 436.6, and 449.9)T_fPBH = 0. Correspondingly, the neutral hydrogen density values as a fraction of the values considering our cooling implementation alone (n_{H0, fPBH = 0}) are (2.3 × 10⁻³, 9.1 × 10⁻⁴, and 5.6 × 10⁻⁴)n_{H0, fPBH = 0}.

For a 33 M_⊙ PBH monochromatic mass function, varying fractions of 10⁻³, 10⁻², 10⁻¹, and 1 lead to temperatures proportional to (1.6, 1.9, 4.2, and 22.0)T_vir, respectively. In terms of T_fPBH = 0, these fractions result in temperatures within regions near the halo center of (298.9, 407.2, 729.1, and 4949.4)T_fPBH = 0. The corresponding neutral hydrogen density fractions are (2.7 × 10⁻³, 1.2 × 10⁻³, 1.0 × 10⁻⁵, and 2.0 × 10⁻⁹)n_{H0, fPBH = 0}.

Similarly, for a 100 M_⊙ PBH monochromatic mass function, fractions of 10⁻³, 10⁻², 10⁻¹, and 1 correspond to temperatures of (1.8, 2.2, 6.1, and 29.0)T_vir, respectively. In terms of T_fPBH = 0, these fractions result in temperatures within regions near the halo center of (353.9, 460.0, 984.1, and 5847.6)T_fPBH = 0. The corresponding neutral hydrogen density fractions are (1.9 × 10⁻³, 8.4 × 10⁻⁴, 1.2 × 10⁻⁶, and 1.4 × 10⁻⁹)n_{H0, fPBH = 0}.

Our results indicate that the changes in temperature in the region closest to the halo center cause significant variations in the density of neutral hydrogen, which is crucial for the formation of molecular hydrogen and, consequently, for star formation. Due to the drastic alterations in gas properties, our findings suggest that PBHs with masses of 1 M_⊙ and fractions greater than or equal to approximately 10⁻² are ruled out. The same applies to PBHs with a mass of 33 M_⊙ and 100 M_⊙ and fractions greater than approximately 10⁻³. This is a more restricted limit than was found by Liu et al. (2022) for 33 M_⊙ PBHs of 10⁻¹.

It is noteworthy that the accretion rate depends on the mass of the PBH. Therefore, there are low PBH masses that do not produce any impact on the gas properties, and hence our method provides no constraints in these cases. We conducted a test using M_PBH = 10⁻¹² M_⊙, a potential mass window for PBH DM. In this scenario, even if the entire DM component is composed of PBHs, the expected heating is negligible and is found not to affect the final gas temperature.

In a near future, we plan to carry out a new simulation that incorporates our model in a self-consistent way with the evolution of the gas component and subsequent stellar formation, which could provide insight into the effects of PBHs on the regulation of the transformation of gas into stars at the very early stages of galaxy formation.

4.2. Variations in the implementation

In addition to the results presented above, we explored alternative approaches in our implementation. This section summarizes our findings.

4.2.1. Energy injection

In addition to the standard SPH method for energy injection, we also tested a spherically symmetric radiative transfer approximation. This involved determining the mean free path as a function of photon energy from the accretion disk. Subsequently, we sorted all gas particles based on their distance to the PBH, and using an interpolated function for the photon energy as a function of the mean free path, we injected energy into each gas particle accordingly. While the computational cost of this method was significantly higher, the main trends observed remained consistent with our previous results.

4.2.2. Halo threshold

As described in Section 3.5, we used a particle selection criterion to identify halo-dominated regions. Regarding the application of the halo density profile, we assumed that the density profile is considered valid when 70% of particles around each DM particle reside within the halo. We also conducted tests using more conservative thresholds, such as 80%. However, due to the resolution limitations of our simulations, imposing a stricter limit often caused the profile to not being applied to the majority of halos. Consequently, this resulted in negligible temperature changes within the halos, even for a significant fraction of massive PBHs.

Notably, opting for lower limits frequently resulted in a substantial expansion of regions classified as ‘halo-dominated’, encompassing even areas far from the halo where the validity of the halo density profile is not expected. Specifically, for a 60% threshold for our halo sample, the count of DM particles encircled by a ‘halo-dominated’ region increased by as much as 50% to 400%, selecting 1.5 to 5 times the actual number of particles in the halo. This sensitivity emphasizes the importance of carefully setting this parameter to avoid overestimating the high-density regions.

4.2.3. Halo-scale analytical implementation

In this section, we study an alternative iteration of our implementation that was specifically crafted to solve the resolution challenge related to gas density in the cooling equations. It is important to note that although we employed the halo density profile from Liu et al. (2022) to calculate the heating rate of PBHs within halos, we refrained from adjusting the individual gas density of CIELO particles according to this profile. This decision was influenced by the mass resolution constraints inherent to these CIELO particles and also for the sake of simplicity. Consequently, the cooling rate in the higher-density regions within the halos might be underestimated.

To explore this limitation, we developed a fully analytical model at a halo scale. We used temperature and density data from Liu et al. (2022), as depicted in the upper panel of Fig. 10. To spatially distribute PBHs within the halo, we integrated a Navarro, Frenk, and White (NFW; Navarro et al. 1997) density profile. The lower panel of Fig. 10 illustrates the DM profile along with the distribution of PBHs at various radii based on their respective masses.

Fig. 10. Temperature, numerical densities, and DM density profiles as functions of distance from the halo center. In the upper panel, the continuous line illustrates the temperature (left y-axis), and the dashed lines represent numerical densities (right y-axis), both as functions of distance from the halo center in units of the virial radius, as obtained from Liu et al. (2022). In the lower panel, the continuous line denotes the DM density profile (NFW) on the left y-axis, and on the right y-axis, the number of PBHs enclosed within each radius for a specific monochromatic mass is shown, as labeled. The halo under consideration is at z = 30.3 and has Mvir = 2 × 105 M⊙, with rvir = 60 pc.

The method we adopted involves the creation of spherical shells within the halo, identifying regions with specific densities and temperatures. Instead of employing the SPH method, we employed the approach outlined in Section 4.2.1. The corresponding heating rates due to PBH accretion were calculated by using the same method as applied in our large-scale implementation and subsequently injected into each gas shell.

It is important to note that this model differs from the previous one in that it does not consider the energy contribution of gas heating due to accretion onto PBHs in the vicinity of the halo. This distinction may lead to an underestimation of heating. Additionally, a characteristic velocity of 10 km/s was fixed for the analysis, different from the previous implementation, where the characteristic velocity of CIELO particles was employed.

However, our primary focus was on cooling. Within this framework, the cooling rates were calculated by considering the high gas densities provided by the Liu et al. (2022) density profile.

The results illustrating the final gas temperature for a halo at z = 30.3, with a virial mass of M_vir = 2 × 10⁵ M_⊙ and a virial radius of r_vir = 60 pc are shown in Fig. 11. In regions of high density, close to the halo center, the temperature increase is generally more pronounced, peaking at densities near 10³ cm⁻³. While high densities are considered in the cooling equations, it is observed that cooling does not surpass heating in the cases of interest. For instance, for M_PBH = 1 M_⊙ and f_PBH = 10⁻⁴, the median temperature calculated using this implementation is 100 K, while the median temperature within the halo region depicted in Fig. 8 is 250 K.

Fig. 11. Final gas temperature inside the halo characterized in Fig. 10 as a function of the initial neutral hydrogen number density of each gas shell. The final temperature is computed considering different PBHs masses and various fractions of DM composed of these PBHs, as labeled. The initial temperature is included with black circles.

We also investigated the sensitivity of the results to variations in the virial halo masses, considering a range within the scale of CIELO halos. The observed trends in the final gas temperature remained consistent with our main results.

5. Summary and conclusions

We developed a semianalytical model that works on top of the CIELO hydrodynamical zoom-in simulation to investigate the influence of PBHs on gas properties within and outside halos at z ∼ 23. We comprehensively examined the evolution of gas properties, including temperature, molecular hydrogen density, and neutral hydrogen density, in the presence of PBHs around z ∼ 23. In the simulated region, we identified six halos with virial masses of about 10⁶ M_⊙.

We explored diverse combinations of PBH masses (1, 33, and 100 M_⊙) and fractions of DM composed of these PBHs (10⁻⁴, 10⁻³, 10⁻², 10⁻¹, and 1). Our results revealed that PBHs with masses of 1 M_⊙ and fractions greater than or equal to approximately 10⁻² induce significant changes in gas properties. The same applies to PBHs with a mass of 33 M_⊙ and 100 M_⊙ and fractions greater than approximately 10⁻³. These alterations could potentially delay star-forming activities to lower redshifts than expected or might even inhibit star formation, particularly within the halos.

Notably, for 33 M_⊙ PBHs with fractions exceeding 10⁻², the final temperature increases by more than ten times compared to the initial temperature, and the neutral hydrogen density decreases to less than 10⁻³ times its initial value toward the halo center. Due to the temperature increase, the cooling time in this case is reduced by about ten times. Similar trends emerge with 100 M_⊙ PBHs; a fraction exceeding approximately 10⁻³ leads to comparable variations in the gas properties to the latter PBH mass, however.

We further explored variations in our implementation, including a spherically symmetric radiative transfer calculation for energy injection. While this method exhibited consistent trends with our previous results, it was computationally more demanding. Additionally, we tested different particle selection criteria to consider a halo density profile. We found that setting the threshold at 70% yielded more applicable profiles for the thin-disk formation in halos with massive PBHs given our simulation resolution. On the other hand, more conservative thresholds, such as 80%, often resulted in negligible changes, underlining the sensitivity of the method to this parameter.

Moreover, we developed a fully analytical model based on an NFW profile and the gas density profile from Liu et al. (2022), allowing us to test the density treatment adopted in our semianalytical model.

In conclusion, our model implementation demonstrated significant variations in gas properties due to the presence of PBHs. Notably, PBHs with masses of 1 M_⊙ and fractions greater than 10⁻² markedly influence gas properties. This significant impact is also observed for PBHs with masses of 33 M_⊙ and 100 M_⊙ and fractions exceeding roughly 10⁻³. One important consideration is the potential influence of the Lyman-Werner radiation produced during PBH accretion, which could dissociate H₂ and HD molecules. This process may result in less efficient cooling, and thus, in potentially higher final gas temperatures and more restrictive allowed fractions. This process will be analyzed in a future work.

A potential caveat regarding the temperature increase observed in our results is that the dynamical reaction of gas to PBH heating was not considered. When simulating heating or cooling and hydrodynamics on the fly, the feedback from PBHs can be self-regulated, so that strong accretion in dense regions will be rare. In this case, the heating by PBHs mostly operates at relatively low densities. When cooling wins at such low densities, the gas will collapse and form compact dense clumps. Some of these clumps may not contain PBHs and could still form stars, a phenomenon referred to as survivorship bias in Liu et al. (2022). In our case, we implicitly assumed that PBHs always exist in dense regions that can potentially form stars, which is not necessarily true in reality because of the discrete nature of PBHs and the self-regulation of BH accretion feedback, as shown in Liu et al. (2022). Moving forward, we aim to address this limitation by conducting self-consistent hydrodynamical simulations that incorporate the PBH accretion feedback model. This simulation will aim to validate our approach and offer deeper insights into the influence of PBHs on star formation processes.

Additionally, in the near future, our focus will be to use this PBH model to calculate the extensive radiation background generated by PBHs across X-ray and radio frequencies, employing semianalytical techniques. Furthermore, we plan to move away from assuming isotropic accretion by PBHs and instead consider the angular momentum of the accreted gas.

Despite the uncertainties in our idealized modeling, it is evident that the feedback exerted by PBHs on the primordial gas in the high-redshift Universe may lead to strong constraints on the properties of a possible PBH dark matter component. This thermal accretion-driven coupling with the gas fundamentally distinguishes PBH dark matter from the more standard particle-based alternatives, giving us an additional empirical probe to ascertain the nature of dark matter.

Acknowledgments

CC is supported by ANID-PFCHA/Doctorado Nacional/2020 – 21202137. PBT acknowledges partial funding by Fondecyt-ANID 1200703/2020, 1240465/2024, and ANID Basal Project FB210003. NDP was supported by a RAICES, a RAICES-Federal and PICT-2021-I-A-00700 grants from the Ministerio de Ciencia, Tecnología e Innovación, Argentina. BL is supported by the Royal Society University Research Fellowship. SP acknowledges partial support from MinCyT through BID PICT 2020 00582 and from CONICET through PIP 2022 0214 RDT thanks the Ministerio de Ciencia e Innovación (Spain) for financial support under Project grant PID2021-122603NB-C21. This project has received funding from the European Union Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement No 734374- LACEGAL. We acknowledge the use of the Ladgerda Cluster (Fondecyt 1200703/2020). The numerical simulation is part of the CIELO project, was performed at the Barcelona Supercomputer Center.

References

Appendix A: Chemistry in the model

For a specified density, temperature, and ionizing background spectrum, we used the slightly modified equations from Katz et al. (1996), accounting for the presence of H₂, to determine the values of n_H0, n_H+, n_He0, n_He+, n_He++, and n_e,

$$\begin{array}{c}\hfill n_{{\mathrm{H}}_0}=\frac{n_{\mathrm{H}}{\alpha }_{{\mathrm{H}}_{+}}}{{\alpha }_{{\mathrm{H}}_{+}}+{\mathrm{\Gamma }}_{e{\mathrm{H}}_0}+\frac{{\mathrm{\Gamma }}_{\gamma {\mathrm{H}}_0}}{n_e}}\end{array}$$(A.1)

$$\begin{array}{c}\hfill n_{{\mathrm{H}}_{+}}=n_{\mathrm{H}}-n_{{\mathrm{H}}_0}-n_{{\mathrm{H}}_2}\end{array}$$(A.2)

$$\begin{array}{c}\hfill n_{\mathrm{He}+}=\frac{y{n}_{\mathrm{H}}}{1+\frac{{\alpha }_{{\mathrm{He}}_{+}}+{\alpha }_d}{{\mathrm{\Gamma }}_{e{\mathrm{He}}_0}+\frac{{\mathrm{\Gamma }}_{\gamma {\mathrm{He}}_0}}{n_e}}+\frac{{\mathrm{\Gamma }}_{e{\mathrm{He}}_{+}}+\frac{{\mathrm{\Gamma }}_{\gamma {\mathrm{He}}_{+}}}{n_e}}{{\alpha }_{{\mathrm{He}}_{++}}}}\end{array}$$(A.3)

$$\begin{array}{c}\hfill n_{{\mathrm{He}}_0}=\frac{n_{{\mathrm{He}}_{+}}({\alpha }_{{\mathrm{He}}_{+}}+{\alpha }_d)}{{\mathrm{\Gamma }}_{e{\mathrm{He}}_0}+\frac{{\mathrm{\Gamma }}_{\gamma {\mathrm{He}}_0}}{n_e}}\end{array}$$(A.4)

$$\begin{array}{c}\hfill n_{{\mathrm{He}}_{++}}=\frac{n_{{\mathrm{He}}_{+}}({\mathrm{\Gamma }}_{e{\mathrm{He}}_{+}}+\frac{{\mathrm{\Gamma }}_{\gamma {\mathrm{He}}_{+}}}{n_e})}{{\alpha }_{{\mathrm{He}}_{+}}}\end{array}$$(A.5)

$$\begin{array}{c}\hfill n_e=n_{{\mathrm{H}}_{+}}+n_{{\mathrm{He}}_{+}}+2n_{{\mathrm{He}}_{++}},\end{array}$$(A.6)

where y = (n_He0 + n_He+ + n_He++)/n_H, and we assumed Γ_γHe0 = Γ_γHe+ = 0.

The collisional ionization rates (Γ_eH0, etc.) and recombination rates (α_H+, etc.) are listed in in Table A.1.

All Tables

All Figures

Fig. 1. Continuous color map of the parameter ṁ across a grid of characteristic velocities and hydrogen number densities considering a fixed PBH mass of 33 M⊙. The dashed black lines indicate the thresholds for MeADAF, MADAF, and MLHAF as a function of the viscosity parameter α = 0.1. The region above MLHAF corresponds to the thin-disk formation regime, while in the region between MLHAF and MADAF, LHAF occurs. The area between MADAF and MeADAF corresponds to the standard ADAF. The eADAF regime exists below MeADAF.

In the text

	Fig. 2. Spectra generated by accreting PBHs of 1, 33, and 100 M⊙ (in yellow, red, and purple, respectively) in different gas density environments (1 cm−3 as dotted lines, 102 cm−3 as dash-dotted lines, and 104 cm−3 as solid lines). This example assumes a constant characteristic velocity $\stackrel{\sim }{v}=10$ km/s and a mean molecular weight μ = 1.22.
In the text

	Fig. 3. Projected distribution of the hydrogen number density (left panel) and temperature (right panel) of gas in the Pehuen halos at z ∼ 23. Each halo considered in this paper is distinguished by a labeled circle with a radius equal to ten times the virial radius of the respective halo.
In the text

Fig. 4. Gas density profile and enclosed mass profile for a 33 M⊙ PBH. The upper panel displays the gas density profile taken from Liu et al. (2022) (dashed black line) alongside the rescaled density profile for one of our CIELO halos (red line). In the lower panel, the enclosed mass profile over the virial mass for this halo is shown. In both panels, the profiles are color-coded based on the gas-accretion rate of a 33 M⊙ PBH. For this particular example, we set $\stackrel{\sim }{v}=10$ km/s and μ = 1.22, but in our implementation, these values varied depending on per-particle information from CIELO. Each panel is segmented into colored zones that align with different radii and their respective gas densities, so that the corresponding accretion mechanisms (eADAF, ADAF, LHAF, and thin disk) are active within the segments. The profiles are color-coded based on the accretion rate of the PBH.

In the text

Fig. 5. Accretion and heating rates of PBHs with varying masses in the density vs. velocity plane. The color-coded panels illustrate the behavior of PBHs with varying masses (1 M⊙, 33 M⊙, and 100 M⊙ columns from left to right) in the density vs. velocity plane, depicting the accretion rate (upper row) and the resulting heating rate (bottom row). The dashed white lines delineate the upper limits of specific accretion regimes (eADAF, ADAF, and LHAF) as shown in Fig. 1. The red rectangles enclose the density and characteristic velocity ranges around each PBH within the halos, employing a rescaled version of the gas density profiles from Liu et al. (2022) and the CIELO particle velocities. The blue rectangles encapsulate weighted values around PBHs in the IGM. The dashed yellow line marks the peak density of a gas particle in the CIELO simulation. These results are representative of gas conditions at z ∼ 23.

In the text

Fig. 6. Key gas properties around dark matter particles in projected spatial coordinates. The upper panels illustrate the spatially projected distribution of dark matter particles in the (x, y) coordinates, and the lower panels depict the distribution in the (x, z) coordinates. The colors indicate the mean values of the weighted (from left to right) relative velocity between the DM particle and the neighboring gas particles, the sound speed in the medium, and the mean molecular weight. The particle distribution is centered on a halo with Mvir = 4.6 × 106 M⊙, Tvir = 8768 K, and rvir = 0.2 kpc, highlighted by cyan circles in each panel.

In the text

Fig. 7. As in Fig. 6, but for the spatial distribution of gas particles and the mean values (from left to right) of the ratio of the temperature obtained by applying only the cooling model (fPBH = 0) and the virial temperature of the halo, the ratio of the temperature considering the entire dark matter component as composed of 33 M⊙ PBHs (fPBH = 1) and the virial temperature of the halo (Tvir = 8768 K), and the final neutral gas density of this last model.

In the text

Fig. 8. Median temperature ratios and gas temperature as a function of distance from the halo center. From top to bottom, each row shows the median ratio of the final (i.e., after implementing our model) and the halo virial temperature, the median ratio of the final gas temperature and the gas temperature when only our cooling model is applied, and the median gas temperature. All as a function of distance from the halo center, expressed as a fraction of the virial radius. Each column from left to right corresponds to considering PBHs with masses of 1, 33, and 100 M⊙, with different fractions of DM composed of these PBHs, as indicated in each panel. For reference, we include the case of fPBH = 0, where we applied our cooling model alone. The error bars represent the standard deviation in each distance bin. Individual profiles were first calculated around each halo separately, and then the median values of the temperature ratios were derived from these individual profiles of the six halos. The virial temperatures range from 3800 to 8800 K.

In the text

	Fig. 9. Same as the middle and lower panels of Fig. 8, but for the neutral hydrogen number density.
In the text

Fig. 10. Temperature, numerical densities, and DM density profiles as functions of distance from the halo center. In the upper panel, the continuous line illustrates the temperature (left y-axis), and the dashed lines represent numerical densities (right y-axis), both as functions of distance from the halo center in units of the virial radius, as obtained from Liu et al. (2022). In the lower panel, the continuous line denotes the DM density profile (NFW) on the left y-axis, and on the right y-axis, the number of PBHs enclosed within each radius for a specific monochromatic mass is shown, as labeled. The halo under consideration is at z = 30.3 and has Mvir = 2 × 105 M⊙, with rvir = 60 pc.

In the text

	Fig. 11. Final gas temperature inside the halo characterized in Fig. 10 as a function of the initial neutral hydrogen number density of each gas shell. The final temperature is computed considering different PBHs masses and various fractions of DM composed of these PBHs, as labeled. The initial temperature is included with black circles.
In the text