# How do Primordial Black Holes change the Halo Mass Function and Structure?

Saiyang Zhang(張賽暘） Department of Physics, University of Texas at Austin,
Austin, TX 78712, USA  Weinberg Institute for Theoretical Physics, Texas
Center for Cosmology and Astroparticle Physics,
University of Texas at Austin, Austin, TX 78712, USA  Volker Bromm Department
of Astronomy, University of Texas at Austin, Austin, TX 78712, USA  Weinberg
Institute for Theoretical Physics, Texas Center for Cosmology and
Astroparticle Physics,
University of Texas at Austin, Austin, TX 78712, USA  Boyuan Liu(劉博遠)
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK

###### Abstract

We examine the effects of massive primordial black holes (PBHs) on cosmic
structure formation, employing both a semi-analytical approach and
cosmological simulations. Our simulations incorporate PBHs with a
monochromatic mass distribution centered around $10^{6}\ \rm M_{\odot}$,
constituting a fraction of $10^{-2}$ to $10^{-4}$ of the dark matter (DM) in
the universe, with the remainder being collision-less particle dark matter
(PDM). Additionally, we conduct a $\Lambda$CDM simulation for comparative
analysis with runs that include PBHs. At smaller scales, halos containing PBHs
exhibit similar density and velocity dispersion profiles to those without
PBHs. Conversely, at larger scales, PBHs can expedite the formation of massive
halos and reside at their centers due to the ‘seed effect’. To analyze the
relative distribution of PBH host halos compared to non-PBH halos, we apply
nearest-neighbor (NN) statistics. Our results suggest that PBH host halos,
through gravitational influence, significantly impact the structure formation
process, compared to the $\Lambda$CDM case, by attracting and engulfing nearby
newly-formed minihalos. Should PBHs constitute a fraction of DM significantly
larger than $\sim$$10^{-3}$, almost all newly-formed halos will be absorbed by
PBH-seeded halos. Consequently, our simulations predict a bimodal feature in
the halo mass function, with most of the massive halos containing at least one
PBH at their core and the rest being less massive non-PBH halos.

Cosmological perturbation theory (341) — Early universe(435) — Dark
matter(353) — Large-scale structure of the universe(902)

††facilities: Lonestar6 (TACC)††software: astropy (Astropy Collaboration et
al., 2013, 2018), Colossus (Diemer, 2018)

##  1 Introduction

Since the James Webb Space Telescope (JWST) deployment in 2022, groundbreaking
observations have begun to illuminate the enigmatic early universe and the
genesis of its first galaxies. Recent JWST observations have revealed an
astonishing population of massive galaxies with stellar masses exceeding
$10^{9}\,\mathrm{M}_{\odot}$ at redshifts $z\gtrsim 10$ (e.g., Finkelstein et
al., 2022; Castellano et al., 2022; Naidu et al., 2022; Yan et al., 2023a;
Castellano et al., 2023; Labbé et al., 2023; Yan et al., 2023b), and
supermassive black holes (SMBHs) with masses $\gtrsim
10^{7}\,\mathrm{M}_{\odot}$ at $z=7-10$ (e.g., Goulding et al., 2023; Bogdán
et al., 2024; Larson et al., 2023; Maiolino et al., 2023; Greene et al., 2024;
Natarajan et al., 2024; Kovacs et al., 2024). While most galaxies are
identified through photometric measurements subject to calibration
uncertainties, a select few have been spectroscopically confirmed (i.e.
Curtis-Lake et al., 2023; Arrabal Haro et al., 2023; Robertson et al., 2023).
The existence of those galaxies, with higher than expected stellar mass at
such high redshift (e.g., GN-z11; Tacchella et al. 2023), challenges our
standard $\Lambda$CDM cosmology (Inayoshi et al., 2022; Boylan-Kolchin, 2023).
To alleviate this tension, a re-calibration of the stellar mass has been
proposed, such as using different stellar initial mass functions and their
resulting spectral energy distributions, when interpreting the observations
(e.g., Finkelstein et al., 2023; Wang et al., 2024).

A related challenge is posed by the discovery and follow-up spectroscopic
confirmation of massive black holes in the pre-reionization universe (e.g.,
UHZ1; Goulding et al. 2023; Bogdán et al. 2024). The standard SMBH formation
channel is via gas accretion onto the light black hole seeds from massive
Population III (Pop III) star remnants with $\sim 100-1000\
\mathrm{M}_{\odot}$ (for reviews, see Smith & Bromm, 2019; Inayoshi et al.,
2020). However, rapid growth of the light seeds to
$10^{7}\,\mathrm{M}_{\odot}$ in less than $\sim 500$ Myr by $z\sim 10$
requires an accretion rate exceeding the Eddington limit (e.g., Bogdán et al.,
2024). Recent observations of high-redshift active galactic nuclei (AGN)
suggest that the majority may be dormant with sub-Eddington accretion rates
(e.g., Jeon et al., 2023). If interspersed with episodic super-Eddington
phases, however, such light-seed models may be able to explain at least some
of the newly discovered high-$z$ SMBHs (Juodžbalis et al., 2024; Maiolino et
al., 2024). An alternative SMBH formation channel proceeds via heavy black
hole seeds from the direct collapse of a primordial gas cloud (Bromm & Loeb,
2003; Begelman et al., 2006; Lodato & Natarajan, 2006). The resulting
accelerated growth (e.g., Jeon et al., 2024) may be required to explain part
of the emerging JWST phenomenology (Bogdán et al., 2024; Kovacs et al., 2024).

To understand the challenges posed by the JWST observations for early
structure formation, we need to consider the properties of the dominant dark
matter (DM) component. Of particular interest are models where a fraction of
dark matter is composed of massive primordial black holes (PBHs), with the
possibility of accelerating the early formation of massive galaxies (e.g., Liu
& Bromm, 2022, 2023; Su et al., 2023; Colazo et al., 2024), or of providing
the initial seeds for high-$z$ SMBHs (e.g., Silk et al., 2024).

PBHs are well motivated DM candidates, theorized to form in the early universe
through various mechanisms (reviewed in Carr & Kühnel, 2020). These include
the collapse of overdensities prior to matter-radiation equality (Zel’dovich &
Novikov, 1967; Hawking, 1971; Escrivà, 2022) and collapse of domain walls
during inflation (Belotsky et al., 2019). While PBHs may not account for all
dark matter, even a sub-dominant mass fraction can significantly influence
cosmic history and imprint observable signals (for a review of PBH
phenomenological constraints, see Carr et al., 2021b). Especially, in
structure formation history, PBHs play an important role through the ‘Poisson’
and ‘seed’ effects (Carr & Silk, 2018; Inman & Ali-Haïmoud, 2019). The
’Poisson’ effect stems from the number density fluctuations of discrete PBHs,
which can accelerate the growth of overdensities and the formation of halos
(Meszaros, 1975; Afshordi et al., 2003; Kashlinsky, 2021; Cappelluti et al.,
2022). On smaller scales, the ‘seed’ effect arises from individual PBHs, where
they act as seeds for SMBH and galaxy formation by accretion of the
surrounding material (Mack et al., 2007; Ricotti, 2007; Ricotti et al., 2008).

Previously, we have studied the influence of stellar-mass PBHs
($10-100\,\mathrm{M}_{\odot}$) on the formation of Pop III stars (Liu et al.,
2022), and the contribution of PBH DM to cosmic radiation backgrounds (Zhang
et al. 2024; see also Ziparo et al. 2022). We have also considered the extreme
case of supermassive ($\sim 10^{9}-10^{11}\,\mathrm{M}_{\odot}$) PBHs on the
assembly of the first galaxies (Liu & Bromm, 2022, 2023). All of these studies
indicate that structure formation is indeed accelerated in PBH cosmologies,
with an earlier formation of DM halos. In this paper, motivated by the recent
JWST observational results, we explore the influence of
$10^{6}\,\mathrm{M}_{\odot}$ PBHs on structure formation, via a combination of
‘Poisson’ and ‘seed’ effects. This PBH mass scale is motivated by a formation
scenario that invokes the softening of the equation of state during
$e^{+}e^{-}$ annihilation, when the universe had cooled to a temperature of
$\sim 1\rm\ MeV$ (Carr et al., 2021a).

We run our simulations till $z\simeq 10$, where our results can be compared
with JWST observations, exploring PBH mass fractions of $10^{-2}$ to
$10^{-4}$, and assuming that the remainder consists of particle dark matter
(PDM). When evaluating the resulting halo mass functions, we interpret any
deviations from the Press-Schechter (PS) formalism (Press & Schechter, 1974)
through a nearest-neighbor (NN) statistical approach (Banerjee & Abel, 2021),
and develop a simple heuristic model to understand the structural changes
induced by different PBH cosmologies.

In Section 2, we outline our analytical framework designed to assess the
influence of PBHs on structure formation and the halo mass function. Following
this, Section 3 details the initial settings and relevant parameters of our
simulations. Based on these simulations, in Section 4, we derive the DM halo
abundances within the different models, and explore the impact of PBHs across
various scales. We conclude in Section 5, where we summarize our findings and
discuss directions for future research.

Throughout the paper, we employ the cosmological parameters from the Planck15
measurement (Planck Collaboration et al., 2016), as follows: $\Omega_{\rm
m}=0.3089$, $h=0.6774$, $\sigma_{8}=0.8159$, $n_{s}=0.9667$.

##  2 PBH Halo Mass Function

For simplicity, we assume a monochromatic mass distribution for PBHs with
$M_{\rm PBH}=10^{6}\,\mathrm{M}_{\odot}$. Observations constrain the abundance
of PBHs as a fraction of total DM mass, denoted by $f_{\rm PBH}$, leading to
maximum values allowed for a monochromatic distribution (Carr et al., 2021b).
The current upper limit on $10^{6}\,\mathrm{M}_{\odot}$ PBHs is $f_{\rm
PBH}\lesssim 10^{-3}$, deriving from the effect of dynamical friction (Carr &
Sakellariadou, 1999) and the formation of large-scale structure (Carr & Silk,
2018). We also note that the most stringent limit arises from $\mu$-distortion
of the cosmic microwave background (CMB) when assuming PBH formation from the
high-density tail of a Gaussian random field (Nakama et al., 2018; Hooper et
al., 2024). However, relaxing the Gaussian fluctuation assumption (Nakama et
al., 2018) and considering formation mechanisms dependent on a different
inflation field (i.e. Kawasaki & Murai, 2019; Carr & Kühnel, 2020) could
alleviate this bound. Therefore, in this paper, we adopt $10^{-3}$ as the
fiducial value, consistent with existing bounds, and vary the PBH fraction in
DM within $f_{\rm PBH}\sim 10^{-2}-10^{-4}$ to study its impact on the halo
mass function (HMF).

Following previous work (Inman & Ali-Haïmoud, 2019; Liu et al., 2022; Liu &
Bromm, 2023; Zhang et al., 2024), we briefly summarize the analytical
formalism for deriving the HMF within PBH cosmologies, and we reproduce the
key equations for the convenience of the reader.

The average (comoving) PBH number density is

$\bar{n}_{\mathrm{PBH}}=\frac{f_{\mathrm{PBH}}}{M_{\rm
PBH}}(\Omega_{\mathrm{m}}-\Omega_{\mathrm{b}})\frac{3H_{0}^{2}}{8\pi G}\mbox{\
,}$

defining $k_{\mathrm{PBH}}=\left(2\pi^{2}\bar{n}_{\mathrm{PBH}}\right)^{1/3}$
as the inter-PBH separation scale, and $\Omega_{\rm b}$ is the energy density
fraction of baryons. PBHs are assumed to be randomly distributed over space,
and the fluctuation in PBH number density leads to the ‘Poisson’ effect,
causing isocurvature fluctuations in the density field (Carr & Silk, 2018).
Previous treatments by Afshordi et al. (2003); Kashlinsky (2016) assumed that
PBHs provide all the DM. Here, we relax this assumption by considering
$f_{\mathrm{PBH}}$ as a free parameter. We express the total power spectrum
$P_{\mathrm{tot}}(k)$, extrapolated to $z=0$, as the combination of adiabatic
$P_{\rm ad}(k)$ and isocurvature $P_{\rm iso}(k)$ modes with a mode mixing
term $P_{\rm corr}(k)$ (Inman & Ali-Haïmoud, 2019; Liu & Bromm, 2023), as

$\displaystyle P_{\mathrm{tot}}(k)$ $\displaystyle=P_{\rm ad}(k)+P_{\rm
iso}(k)+P_{\rm corr}(k)\,,$ (1)

where $P_{\rm ad}(k)$ refers to the standard $\Lambda$CDM power spectrum,
generated with the python package colossus (Diemer, 2018). The isocurvature
term $P_{\rm iso}(k)$ arises from the Fourier transform of the overdensities
attributable to the discrete nature of PBHs:

$P_{\rm iso}(k)=\left[f_{\rm PBH}D_{0}\right]^{2}/\bar{n}_{\rm PBH},$ (2)

and it becomes significant at scales smaller than the matter-radiation
equality scale $k_{\rm eq}\sim 0.01\,\mathrm{Mpc}^{-1}$. Here $D_{0}=D_{\rm
PBH}(a=1)$ is the growth factor of PBH-induced perturbations at $z=0$, well
approximated by the analytical formalism in Inman & Ali-Haïmoud (2019):

$\displaystyle D_{\rm PBH}(a)\approx\left(1+\frac{3\gamma a}{2a_{-}a_{\rm
eq}}\right)^{a_{-}}\ ,\gamma=\frac{\Omega_{m}-\Omega_{b}}{\Omega_{m}}\ ,$ (3)
$\displaystyle a_{-}=\frac{1}{4}\left(\sqrt{1+24\gamma}-1\right)\ \ ,$

with $a_{\rm eq}\sim 1/3400$ representing the scale factor at matter-radiation
equality.

For the correlation between adiabatic and isocurvature modes, given by $P_{\rm
corr}(k)$, we adopt the heuristic formula initially developed by Liu et al.
(2022) and updated in Liu & Bromm (2023), applied to the range $k_{\rm
eq}<k<3k_{\mathrm{PBH}}$:

$\displaystyle P_{\rm
corr}(k)=f_{\mathrm{PBH}}D_{0}\left(k/k_{\mathrm{PBH}}\right)^{3}P_{\rm
ad}(k).$ (4)

Further work will be needed to fully explore this heuristic term.

By convolving the total power spectrum with a top-hat filter, we apply the PS
formalism (Press & Schechter, 1974) to compute the HMF, $dn/dM_{\rm h}$,
incorporating corrections for ellipsoidal collapse (Sheth & Tormen, 1999). In
Figure 1, we show the total power spectrum including the PBH isocuvature and
correlation terms, calculated with Equ. (1)-(4), together with the resulting
halo mass function as a function of $f_{\rm PBH}$. At the lower mass end of
the HMF, our analytical formalism exhibits a more complex pattern, while at
the high mass end, it asymptotically approaches the $\Lambda$CDM limit with
smaller $f_{\rm PBH}$.

Figure 1: Structure formation within PBH cosmologies. Left panel: Power
spectrum of the DM density field for different PBH mass fractions, $f_{\rm
PBH}=10^{-4}$ (pink), $10^{-3}$ (green), 0.01 (orange), calculated following
Equ. (1)-(4), in comparison with the standard $\Lambda$CDM case (blue)
measured by Planck Collaboration et al. (2016). Right panel: The analytical
HMF at $z=10$, employing the same conventions for the line color as before. It
is evident that the excess power on small scales boosts the number density of
low-mass DM halos, when PBHs are present.

##  3 Numerical Methodology

###  3.1 Simulation Box and Initial Parameters

In this work, we focus on the formation of large-scale structure within PBH
cosmologies, limiting our simulations to dark matter (DM) only111Accordingly,
we set $\Omega_{\rm b}=0$ in Equ. (3) to calculate the growth factor. This
causes a small difference in the HMF derived from the PS formalism compared
with the case using the Planck15 value $\Omega_{\rm b}=0.0486$ (Planck
Collaboration et al., 2016)..Consequently, we defer the more complex
simulations involving baryonic star formation and feedback physics to future
studies. Our simulations are conducted using the state-of-the-art simulation
code gizmo (Hopkins & Raives, 2016), incorporating the Tree+PM gravity solver
from gadget-3 (Springel, 2005). In total, we have executed 8 simulations under
various PBH initial configurations, as detailed in Table 1.

Table 1: Summary of key properties and parameters for the eight simulations. $m_{\rm PBH}$ represents the initial PBH mass (assuming all PBHs have the same initial mass), set at $10^{6}\,\mathrm{M}_{\odot}$ throughout. $f_{\rm PBH}$ denotes the mass fraction of PBHs relative to the total DM content. $N_{\rm gen}$ indicates the number of PBHs generated from the initial conditions. NONLIN is a flag indicating whether the initial non-linear perturbations of PDM particles by PBHs are included (✓) or not (✗). Similarly, CORR is a flag indicating whether the growth factor of PBH induced perturbations is corrected to reproduce the growth of halos seeded by individual PBHs from the spherical collapse theory (Mack et al., 2007). The final two parameters, $\epsilon_{\rm DM}$ and $\epsilon_{\rm PBH}$, represent the (co-moving) softening lengths for PDM and PBH particles, respectively. Run | $m_{\rm PBH}\ [\rm M_{\odot}]$ | $f_{\rm PBH}$ | $N_{\rm gen}$ | NONLIN | CORR | $\epsilon_{\rm DM}[\rm kpc/h]$ | $\epsilon_{\rm PBH}[\rm kpc/h]$
---|---|---|---|---|---|---|---
CDM | - | - | - | - | - | 0.05 | -
PBH2 | $10^{6}$ | $10^{-2}$ | 989 | ✓ | ✓ | 0.05 | 0.002
PBH3 | $10^{6}$ | $10^{-3}$ | 97 | ✓ | ✓ | 0.05 | 0.002
PBH4 | $10^{6}$ | $10^{-4}$ | 10 | ✓ | ✓ | 0.05 | 0.002
PBH3_PS | $10^{6}$ | $10^{-3}$ | 97 | ✗ | ✓ | 0.05 | 0.002
PBH3_NS | $10^{6}$ | $10^{-3}$ | 97 | ✓ | ✗ | 0.05 | 0.002
PBH3_PS_NS | $10^{6}$ | $10^{-3}$ | 97 | ✗ | ✗ | 0.05 | 0.002
PBH3_smooth | $10^{6}$ | $10^{-3}$ | 97 | ✓ | ✓ | 0.05 | 0.5

To select our simulation box, we consider a trade-off between resolution and
statistical significance. We aim to generate a substantial number of PBHs
within our box. Concurrently, we seek to resolve the forces between PDM and
PBH particles, necessitating that the mass of PDM particles be significantly
lower than that of PBH particles, thereby requiring high simulation
resolution. Given the constraints of data storage and computational resources,
we limit the total particle number to $N=256^{3}$. For PBHs, we consider mass-
dependent constraints on their abundance (for a general review, see Carr et
al., 2021b). We adopt a fiducial value of $f_{\rm PBH}=10^{-3}$ (PBH3),
centered at $M_{\rm PBH}=10^{6}\ \mathrm{M}_{\odot}$. To generate
approximately $\mathcal{O}(10^{2})$ PBHs within the box, we select a box size
of $\sim 1\ \mathrm{(Mpc/h)}^{3}$ in comoving units. Consequently, the mass of
PDM particles within this box is $\sim 5,100\ \mathrm{M}_{\odot}/\mathrm{h}$,
which also meets the force resolution requirement. We generate the initial
conditions using the MUSIC code (Hahn & Abel, 2011) at an initial redshift of
$z_{\rm ini}=1/a_{\rm ini}-1=300$, where the growth of structure is still
linear, and terminate our simulations at $z=10$, above which the box size is
marginally representative. To facilitate a comprehensive comparison of the
impact of PBHs on structure formation, we also consider the standard
$\Lambda$CDM case (CDM) and other potential PBH fractions with the same PBH
mass: $f_{\rm PBH}=10^{-2}$ (PBH2), and $10^{-4}$ (PBH4).

In our simulations, we set the softening length for PDM particles to
$\epsilon_{\rm DM}\sim 0.01L/N^{1/3}\simeq 0.05\ \mathrm{kpc/h}$, emulating
the collisionless nature of PDM. The softening length for PBHs is chosen to be
$\epsilon_{\rm PBH}\sim 0.002\ \mathrm{kpc/h}$, simulating the dynamics of a
point particle to better resolve the small-scale structure around PBHs.
Additionally, we explore an extreme case for PBHs where their discrete nature
is smoothed out (PBH3_smooth), with a softening length for PBHs of $\sim 0.5\
\mathrm{kpc/h}$, to examine the influence of this parameter on structure
formation.

###  3.2 PBH initial conditions

We establish the initial conditions for simulations including PBHs based on
the analytical framework presented by Inman & Ali-Haïmoud (2019), and the
numerical recipes developed in Liu et al. (2022); Liu & Bromm (2023), taking
into account the density perturbations induced by PBHs across various scales.
For the reader’s convenience, we again briefly summarize the procedure (Liu et
al., 2022; Liu & Bromm, 2023), reproducing key expressions from these studies.

Initially, we add PBH particles into the $\Lambda$CDM cosmology initial
conditions generated by music (Hahn & Abel, 2011), according to the following
algorithm:

  * •

Divide the simulation box into $N_{\rm grid}$ grid cells, each containing
approximately one PBH.

  * •

Iterate over all cells. For the $i$-th cell, randomly draw a number $N_{i}$,
following a Poisson distribution, as the total PBHs within the cell and
randomly place them within the cell, noting their positions. Assign each PBH
the velocity of its nearest neighbor PDM particle.

  * •

Reduce the mass of each PDM particle by a fraction $\simeq f_{\rm PBH}$ for
mass conservation.

Next, we incorporate the isocurvature perturbation term induced by PBHs into
the initial positions and velocities of PDM particles, denoted by
$\left[T_{\rm iso}(a)-1\right]f_{\rm PBH}\delta_{\rm iso}^{0}$. Utilizing the
Zel’dovich approximation (Zel’dovich, 1970), we calculate the comoving
displacement $\vec{\psi}$ and velocity $\Delta\mathbf{v}$ fields, resulting
from the comoving acceleration field $\nabla\phi_{\rm iso}(\mathbf{x})$, as
follows:

$\displaystyle\vec{\psi}(\vec{x})$ $\displaystyle=-\frac{2D_{\rm PBH}(a_{\rm
ini})}{3\Omega_{m}H_{0}^{2}}\nabla\phi_{\rm iso}(\vec{x})\ ,\quad$ (5)
$\displaystyle\Delta\vec{v}(\vec{x})$ $\displaystyle=\frac{a_{\rm
ini}\dot{D}_{\rm PBH}(a_{\rm ini})}{D_{\rm PBH}(a_{\rm
ini})}\vec{\psi}(\vec{x})\,$ $\displaystyle\nabla\phi_{\rm iso}(\vec{x})$
$\displaystyle=4\uppi Gm_{\rm
PBH}\sum_{i}\frac{\vec{x}-\vec{x}_{i}}{|\vec{x}-\vec{x}_{i}|^{3}}\,,$

where $\vec{x}_{i}$ represents the comoving coordinate of the $i$-th PBH at
$a_{\rm ini}$, and the growth factor $D_{\rm PBH}(a_{\rm ini})$ is given by
Equ. (3). The comoving acceleration field $\nabla\phi_{\rm iso}(\mathbf{x})$
accounts for contributions from all PBHs. Applying these equations, we adjust
the positions and velocities of each PDM particle $j$ from its original
coordinate $\vec{x}_{j}$ and velocity $\vec{v}_{j}$ in the $\Lambda$CDM
initial conditions to $\hat{\vec{x}}_{j}=\vec{x}_{j}+\vec{\psi}(\vec{x}_{j})$
and $\hat{\vec{v}}_{j}=\vec{v}_{j}+\Delta\vec{v}(\vec{x}_{j})$, respectively.

Given the ‘seed’ effect (as proposed by Mack et al., 2007), perturbations near
PBHs are nonlinear from the outset. For simplicity, we define a nonlinear
scale for a single PBH via $d_{\rm nl}=[D_{\rm PBH}(a_{\rm ini})/\bar{n}_{\rm
PBH}]^{1/3}$, where $\bar{n}_{\rm PBH}$ is the average PBH number density. To
consistently model both ‘seed’ and ‘Poisson’ effects, we introduce two
numerical treatments. On scales significantly larger than $d_{\rm nl}$, we can
effectively use linear approximation to the initial conditions, by applying a
truncated displacement field as

$\displaystyle\vec{\psi}_{\rm Poisson}=\min(1,d_{\rm
PDM}/|\vec{\psi}|)\vec{\psi},$ (6)

where $d_{\rm PDM}$ is the average distance between PDM particles as given by
the initial resolution of the simulation.

In regions very close to PBHs ($r\lesssim d_{\rm nl}$), we account for
nonlinear structures with stronger perturbations than allowed by truncation,
as indicated by NONLIN in Table 1. Considering a virialized halo as the
outcome of the ‘seed’ effect with an overdensity $\Delta_{\rm vir}=18\pi^{2}$,
we apply a specific truncation to capture the nonlinear structure:

$\displaystyle\vec{\psi}_{\rm seed}$ $\displaystyle=\min(1,f_{\rm
shrink}r/|\vec{\psi}|)\vec{\psi},\quad$ (7) $\displaystyle f_{\rm shrink}$
$\displaystyle=1-\Delta_{\rm vir}^{-1/3}\approx 0.822.$

Here, $r$ is the distance of the PDM particle from its nearest PBH. With this
prescription, we create a uniform sphere of PDM surrounding the PBH with an
overdensity of $\Delta_{\rm vir}$ and a radius $d_{\rm nl}$, to approximate
the nonlinear features. By running our gizmo simulations (see the next
section), we find that the mass of the virialized halo around an isolated PBH
grows linearly with $a$, as $M_{h}\sim(4\pi/3)D_{\rm PBH}(a)m_{\rm PBH}$,
which is generally consistent with the spherical collapse theory (see Mo et
al., 2010), where $M_{\rm h}\sim m_{\rm PBH}a/a_{\rm eq}$ (Mack et al., 2007),
but the normalization is off by a factor of a few. To mitigate the difference
in halo mass between the analytical predictions by Mack et al. (2007) and our
simulations using the growth factor from Equ. (3), we correct the growth
factor to $D_{\rm PBH}(a_{\rm ini})\sim 3a_{\rm ini}/(4\pi a_{\rm eq})$, as
indicated by CORR in Table 1. Finally, we combine the two perturbation regimes
to calculate the displacement field induced by PBHs as
$\vec{\psi}=\vec{\psi}_{\rm seed}+\vec{\psi}_{\rm Poisson}$, capturing both
‘seed’ and ‘Poisson’ effects.

For completeness, we conduct simulations with the same PBH parameters but
without correcting the growth factor (PBH3_NS), without the seed effect
(PBH3_PS), and without both (PBH3_PS_NS)222Specifically, when the NONLIN flag
is turned on, the Poisson term $\vec{\psi}_{\rm Poisson}$ include the
contributions from all PBHs but the nearest one, while the seed term
$\vec{\psi}_{\rm seed}$ comes from the nearest PBH. If NONLIN is turned off,
we only consider the Poisson term contributed by all PBHs.. The latter case
replicates our prescription from previous work (Liu et al., 2022).

##  4 Results and Discussion

###  4.1 Halo Statistics in PBH cosmologies

To derive the HMF from simulation snapshots, we rely on the rockstar halo
finder (Behroozi et al., 2013) to assign particles to halos. However, one
caveat is that when reading the snapshot files, rockstar only considers PDM
particles in calculating the halo mass. Therefore, we need to manually assign
the mass of PBHs to their respective host halos in our calculations. To
exclude less massive halos that contain too few PDM particles and are thus
poorly resolved dynamically, we only consider halos above the molecular
cooling threshold, set by $M_{\rm h}\gtrsim 1.5\times
10^{5}\,\mathrm{M}_{\odot}\left(\frac{1+z}{31}\right)^{-2.1}$ (Trenti &
Stiavelli, 2009). Given that this mass is typically much smaller than the
halos seeded by PBHs in the redshift range of $z\simeq 10-30$, we employ it as
an approximate lower mass limit for PDM-only halos333Previously, we derived
the star-forming minihalo mass threshold within the stellar-mass ($10-100\ \rm
M_{\odot}$) PBH regime by comparing the molecular/radiative cooling vs. the
heating rates by PBH accretion (Liu et al., 2022). Calculating the precise
lower mass limit for the host halos where first star formation is triggered
under the influence of PBH accretion feedback within other PBH mass regimes
will require hydrodynamical simulations to assess cooling and heating, which
is beyond the scope of this study., defined as halos containing no PBHs.

Initially, when analyzing the PBH simulations, we observed that halos are
positioned closer together, with some overlapping at small scales, compared to
the morphology in the CDM (pure PDM) simulation. At the small scales involved,
baryonic effects will be important, in particular the feedback from star
formation (e.g., Brooks et al., 2013; Jaacks et al., 2019). As simulations
including baryonic feedback are beyond the scope of our current study, and the
numerical prescriptions used to produce cosmological initial conditions
including PBHs are imperfect (Liu & Bromm, 2023, see the discussion in their
Sec. 2.2), we cannot accurately account for the statistics of subhalos, and we
defer a more detailed analysis of subhalo structure within PBH cosmologies to
future work. For simplicity, we therefore remove all subhalos from our halo
catalog, retaining only their most massive hosts. Next, we summarize our
algorithm to address this cloud-in-cloud problem and to remove subhalos:

  * 1.

Iterate over each PBH in the PBH catalogue, indexed by $i$. For the $i$-th
PBH, iterate over the halo catalog indexed by $j$. Identify and record the
index of the most massive host halo by calculating the Euclidean distance from
this PBH to the center of the $j$-th halo and comparing it with the virial
radius of the $j$-th halo.

  * 2.

Iterate over the halo catalog indexed by $j$. For each halo that is not in the
local maximum-mass host list, calculate the Euclidean distance between the
center of $j$-th halo and the center of the local maximum host. Add this halo
to the removal list, if this separation distance is smaller than the virial
radius of the local maximum host.

  * 3.

Iterate over PBHs, indexed by $i$, to find the index of their most massive
hosts. Add the PBH mass to these most massive hosts.

  * 4.

Calculate the pairwise separations among the remaining halos that are not
local maximum hosts. For each pair, the less massive halo is removed if the
separation is smaller than the virial radius of the more massive halo.

On the other hand, for the CDM run, we only apply step 4 to remove all
subhalos from the halo catalogue, to facilitate a meaningful comparison with
the PBH cases. This method allows us to systematically exclude all subhalos
within each host halo. Clearly, developing a halo finder capable of properly
capturing the subtleties introduced by PBHs in cosmic structure formation is
an important goal for future work.

###  4.2 Large-Scale Structure

Figure 2: Distribution of matter and PBHs. We project the simulation box onto
a 2D plane to depict the cosmic web at $z=10$ (in comoving units). Left panel:
Cosmic web in the baseline CDM simulation. Right panel: Distribution for the
PBH3 run, with $f_{\rm PBH}=10^{-3}$. In this visualization, black dots
represent the positions of PBHs, while red circles indicate the positions of
halos, with their radii corresponding to the comoving virial radii of the
halos. The color bar denotes the number of PDM particles stacked within each
pixel. In the right panel, all subhalos within their more massive hosts have
been removed to better showcase the spatial distribution of halos and PBHs.
Figure 3: Comparing the simulated halo mass function (dashed lines) from
different initial conditions after removing all subhalos at $z=10$. To provide
a reference point, we also reproduce the HMF derived from the analytical
formalism in Section 2 (solid lines) at the same redshift, noting that the
same line color indicates identical initial conditions and PBH parameters.
Left panel: Effect of varying the fraction of PBHs in the DM, specifically
considering $f_{\rm PBH}=10^{-2}$ (PBH2), $10^{-3}$ (PBH3), and $10^{-4}$
(PBH4). Right panel: Effect of varying numerical implementation and
parameters, for a fixed $f_{\rm PBH}=10^{-3}$. Specifically, we consider the
absence of the correction of the growth factor (PBH_NS), non-linear
perturbations from the nearest PBH (PBH_PS), and of both (PBH_PS_NS).
Additionally, the softening parameter for PBH particles is increased in the
PBH3_smooth case. In both panels, the horizontal dashed line represents the
number density resolution limit set by the comoving volume of our box as
$1/V_{\rm com}$.

After removing all subhalos, we generate projection plots at the last snapshot
of the CDM and PBH3 simulations at $z=10$, depicted in Figure 2. These plots
illustrate the structural differences when PBHs constitute $f_{\rm
PBH}=10^{-3}$ of the total DM, compared to a standard $\Lambda$CDM simulation.
On a larger scale, more massive halos have formed by this redshift in the PBH3
simulation, while significantly fewer have formed in the CDM case.
Additionally, the PBH3 run demonstrates clear segregation in halo mass, with
almost all massive halos containing PBHs and being substantially larger than
their pure PDM halo counterparts. Furthermore, the abundance of less massive
halos is suppressed compared to the CDM run. Structurally, the average spacing
between halos is larger in the PBH3 case. We also examined these effects in
the PBH2 and PBH4 runs to understand how PBH abundance impacts structure
formation. Surprisingly, almost all halos in the PBH2 run contain at least one
PBH, and pure PDM halos have almost vanished. The PBH4 simulation, however,
more closely resembles the CDM one, displaying a mix of clumpy pure PDM
($\sim$$10^{6}\,\mathrm{M}_{\odot}$) minihalos and a few PBH-hosting halos.
For further details on how different PBH mass fractions $f_{\rm PBH}$
influence structure formation, we refer the reader to Figure 9 in Appendix A.

Subsequently, we count the halos within each mass bin and plot the HMF after
subhalo removal at $z=10$, as shown in Figure 3. For ease of comparison, we
also reproduce the analytical HMFs, as derived in Section 2. The left panel of
Figure 3 illustrates the effects of varying the PBH mass abundance parameter
$f_{\rm PBH}$, and the right panel displays the effects of different initial
conditions or numerical implementations for the same $f_{\rm PBH}$. Here, all
simulations are represented with dashed lines. The sharp cutoff in our
simulated results can be attributed to the halo selection mass limit at the
lower end and the box size limit at the higher end. From the left panel, we
find that with higher $f_{\rm PBH}$, the abundance of halos with mass $\ll
10^{8}\mathrm{M}_{\odot}$ is greatly suppressed, in contrast to the analytical
results. Comparing the simulation results in the right-hand panel, the
differences are minimal, indicating that variations in the choice of numerical
parameters do not greatly change the HMF. However, in the PBH3_NS
implementation, we observe that the mass of the most massive halos is
approximately twice as large as in the PBH3 run. This effect is likely due to
the larger structure growth factor $D_{\rm PBH}$ without the correction
motivated by the results in Mack et al. (2007, see Sec. 3.2), resulting in
heavier PBH-seeded halos. At the lower mass end ($\ll
10^{8}\mathrm{M}_{\odot}$), the halo abundance in the PBH3_NS run is
suppressed by a factor of 2 as compared to the PBH3 case. Both panels suggest
that more massive halos suppress the abundance of less massive halos, a topic
we will elaborate on in more detail in Section 4.4.

###  4.3 Small-Scale Structure

Figure 4: Structure of massive halos in the CDM and PBH3 runs at $z=10$. From
the PBH3 run, we specifically select the ten most massive halos that do not
contain PBHs (red, mid panels), alongside with the ten most massive PBH
hosting halos (blue, right panels). For comparison, we also select 10 most
massive halos from the CDM run without any influences from PBHs (red, left
panels). The opaque curves represent the most massive halos among the selected
ones, whereas the others are represented by quasi-transparent curves. Top
panels: Halo density profile. Here, we normalize $\rho(r)$ to the average
cosmic density $\rho_{\rm crit}(z)$ at $z=10$, and represent the physical
distance to a halo’s center, $r$, normalized by the halo’s virial radius
($r_{\rm vir}$). Bottom panels: Velocity dispersion profile, given by
$\sigma_{v}$ vs. $r/r_{\rm vir}$.

Among the halo catalogue for the PBH3 run, we find that almost all of the more
massive halos ($\gtrsim 10^{8}\,\mathrm{M}_{\odot}$) host PBHs, with only one
being a pure PDM halo. To fairly compare PBH hosting halos and pure PDM ones,
we select the ten most massive pure PDM halos with mass range of $\sim 3\times
10^{6}-2\times 10^{8}~{}\mathrm{M}_{\odot}$ and compare them with the ten most
massive PBH host halos ranging from $\sim 3-9\times
10^{8}~{}\mathrm{M}_{\odot}$. In this way, the best resolved PDM halos are
selected for comparison, as seen from Figure 4. For comparison, we also select
the ten most massive PDM halos with mass range of $\sim 4\times 10^{7}-1\times
10^{8}~{}\mathrm{M}_{\odot}$ from the CDM run as the simulation without PBH.
Our results are shown in Figure 4. We observe that a few selected PDM halos
exhibit a constant density near the center, forming a core-like profile. These
core-flattened halos are interacting with the nearest more massive PBH host
halos, resulting in perturbed density profiles.

The remaining PDM and PBH host halos typically scale as $\rho\propto r^{-1}$
for the central profile. At larger distances from the halo center, both PBH
host and PDM halos exhibit a steeper profile, with $\rho\propto r^{-2.2}$,
largely consistent with the Navarro-Frenk-White (NFW) profile. A similar
scaling is also observed when sampling PDM halos in the CDM run for
comparison. We also examine the velocity dispersion versus normalized radius
in the right panel of Figure 4. We find that the overall velocity dispersion
is higher for PBH hosts, as the selected cases are more massive, and we note
that PDM halos show a more flattened radial profile.

Given that this work involves a DM-only simulation, whereas baryonic feedback
is known to dominate on small scales (e.g., Brooks et al., 2013; Wetzel et
al., 2016), we defer further studies involving baryonic physics to future work
in order to better simulate smaller-scale structure and stellar feedback
processes within PBH cosmologies.

###  4.4 Nearest-Neighbor Statistics

Figure 5: Nearest-neighbor (NN) statistics for standard $\Lambda$CDM case,
represented by our CDM simulation. We show the relationship between halo mass
and the mass of its NN halo (top panels), as well as the physical distance to
the NN (bottom panels), for all halos in the simulation box (represented by
blue circles), over $z\simeq 25-10$. Additionally, we indicate the virial
radius for different halo masses (green lines). For snapshots larger than a
redshift of $\sim 25$, the scarcity of halos precludes the generation of such
plots. Figure 6: NN analysis for the PBH3 simulation, where the PBH density
fraction is $f_{\rm PBH}=10^{-3}$. The same conventions are employed as in
Fig. 5. Here, four distinct marker colors are used to differentiate between
halos and their NNs based on their association with PBHs: No-PBH halo with No-
PBH halo NN (blue), No-PBH halo with a PBH halo NN (pink), PBH halo with No-
PBH halo NN (purple), and PBH halo with a PBH NN halo (red). The distinction
between the halos that contain PBHs and the pure PDM ones is clearly visible.

Next, we explore a possible way to explain key features in the HMF using
nearest-neighbor (NN) statistics, an attractive choice due to its minimal
parameter dependence, as previously explored by Banerjee & Abel (2021) in the
context of large-scale sky surveys. We plot the halo mass against the mass of,
and distance to, its nearest neighbor halo in Figures 5 and 6, representing
the CDM and PBH3 cases, respectively, to examine the clustering patterns and
their evolution through cosmic time.

In the CDM run (Figure 5), we observe a hierarchical merger pattern, with most
of the minihalos with $\sim 10^{6}\,\mathrm{M}_{\odot}$ forming near each
other and merging into larger halos with $\gtrsim 10^{8}\,\mathrm{M}_{\odot}$
as they evolve with redshift. In contrast, in the PBH3 run, we categorize
halos into four distinct NN pairings: PDM halo with PDM halo, PDM halo with
PBH host halo, PBH host halo with PDM halo and PBH host halo with PBH host
halo. The top panels of Figure 6 show a clustering pattern that indicates the
distinction between PDM halos and PBH host halos. The time evolution of the NN
distances show a more complex pattern, as the number of PDM halos paired with
PBH host halos and vice versa decreases with decreasing redshift.

Further details on NN statistics for additional simulations, exploring our
suite of parameters and initial setups, can be found in Appendix A. NN plots
for the PBH_PS and PBH2 runs (Figures 10 and 11 in Appendix A) exhibit a
consistent evolutionary pattern. However, in the PBH4 run, PBH host halo
pairings with other PBH host halos vanish as redshift decreases, indicating
that all PBH halos are surrounded by PDM halos and disruption of PDM halos by
PBH host halos is less pronounced. For PBH abundances of $f_{\rm PBH}\gg
10^{-3}$ almost all halos host PBHs, and PDM halos formed at earlier redshifts
nearly vanish. Conversely, for $f_{\rm PBH}\ll 10^{-3}$, the clustering
pattern (Figure 12) is similar to that of the CDM run (Figure 5). Moreover,
removing the ‘seed’ effect in the initial condition does not significantly
affect the clustering pattern, as shown in Figure 10.

Figure 7: Halo assembly in the early universe. In the $\Lambda$CDM case,
minihalos are expected to cluster and undergo hierarchical merging to form
larger mass halos (top panel). However, PBHs significantly disrupt this
hierarchical merger process for minihalos formed close to PBH-seeded halos
(bottom panel). With increasing PBH abundance, almost no minihalo can
withstand this ingestion, and nearly all are likely to be absorbed by the
nearest PBH-seeded halo. This dynamic interaction underscores the pivotal role
of PBHs in altering the standard view of cosmic structure formation through
their impact on hierarchical merging within the PBH-CDM setting. Figure 8:
Schematic representation of halo dynamics within a PBH-CDM universe, modeling
their trajectories relative to their nearest PBH-seeded halo. The vertical
axis of ‘Initial Position’ refers to the starting physical distance $r_{0}$ of
a PDM halo in proximity to its closest PBH-seeded halo. We show the regime
where the trajectories of PDM halos are (not) affected by PBHs in yellow
(light orange) according to Equ. (9). We further indicate the zones under PBH
tidal influence (pink), and boundaries of isolated PBH-seeded halos (blue). To
illustrate the dependence of PDM halo dynamics on the PBH mass fraction
$f_{\rm PBH}$, we plot exemplar trajectories of PDM halos starting from
$r_{0}=0.2d_{\rm PBH}$ for $f_{\rm PBH}=10^{-2}$ (red), $f_{\rm PBH}=10^{-3}$
(orange), and $f_{\rm PBH}=10^{-4}$ (purple), indicating that more PDM halos
will be disrupted/engulfed by PBH-seeded halos when $f_{\rm PBH}$ is higher.

To further explore the phenomenology in our NN statistics, we consider an
idealized dynamical model for the motion of PDM halos around nearby PBH-seeded
halos, to explain the bimodal feature in the halo mass function. We here
specifically consider the competition between the self-gravity of a given
region and the Hubble expansion in the presence of PBHs (schematically
illustrated in Fig. 7). For a PDM halo that forms at redshift $z_{0}$ at a
physical distance $r_{0}$ from its NN PBH host halo, we assume that $r_{0}$ is
sufficiently large for the two halos to be effectively treated as point
masses.

Further assuming a vanishing velocity at formation, the acceleration of a PDM
halo at a proper distance $r$ from its NN PBH host halo, and at cosmic time
$t$, can be effectively expressed as444For a similar argument, see Equ.(5) in
Ali-Haïmoud et al. (2017), in the context of binary mergers.

$\frac{d^{2}r}{dt^{2}}=r(H^{2}(t)+\dot{H}(t))-\frac{G\bar{M}_{\rm
h,PBH}(t)}{r^{2}}\mbox{\ .}$ (8)

For simplicity, $\bar{M}_{\rm h,PBH}(t)$, the average mass of the PBH-seeded
NN halo is approximated by $\bar{M}_{\rm h,PBH}\sim m_{\rm PBH}a/a_{\rm eq}$,
based on the spherical collapse theory involving an isolated PBH (Mack et al.,
2007). Numerically solving Equ. (8), we identify a critical initial physical
distance, beyond which the trajectory of a newly-formed halo is hardly
affected by the presence of PBHs, implying that its merger and accretion
history is not impacted as well, approximately given by $r_{\rm crit}\simeq
30\,{\rm kpc}\,[(1+z_{0})/11]^{-1}$.

We first compare this distance to the Hill radius (Shu, 1982), approximated as
$R_{\rm Hill}=2.4R_{\rm vir}(\bar{M}_{\rm h,PBH})\simeq 5.5\,{\rm
kpc}\,[(1+z_{0})/11)]^{-4/3}$, which represents the maximum radius for tidal
influence from the PBH host halos. Since the tidal influence sphere is much
smaller than the region (defined by $r_{\rm crit}$) where a PBH-seeded halo
can disrupt the merger and accretion histories of its neighbors, the latter
effect dominates at larger scales in the formation and growth of halos leading
to a bimodal feature in the HMF. Therefore, the fate of a PDM halo can be told
by comparing $r_{0}$ with $r_{\rm crit}$:

$r_{0}\begin{cases}\gg r_{\rm crit}:&\text{PBH non-influenced region}\\\ \ll
r_{\rm crit}:&\text{PBH influenced region,}\end{cases}$ (9)

i.e., $r_{\rm crit}$ divides the merger histories of PDM halos into those that
will be significantly altered (influenced) by the presence of a nearby PBH,
and those that will not.

Therefore, this critical radius for PBH influenced regions sets a bifurcation
limit for the inter-PBH separation and thus the mass fraction of PBHs within
DM. If, on average, PBHs are too closely spaced, this will result in overlaps
of the perturbation regions. In this case, trajectories of almost all PDM
halos will be impacted. Equating the critical radius $r_{\rm crit}$ to half
the inter-PBH separation (proper) distance $d_{\rm
PBH}\simeq(1/\bar{n}_{\mathrm{PBH}})^{1/3}(1+z)^{-1}\simeq 27\,{\rm
kpc}\,[(1+z_{0})/11]^{-1}[f_{\rm PBH}/10^{-3}]^{-1/3}$, we calculate a rough
bifurcation threshold for the PBH abundance (independent of $z_{0}$) as
$f_{\rm PBH}\simeq 10^{-4}$. If the PBH abundance is much larger than this
limit, almost no PDM halos will survive if formed close to a PBH, and halos
that have already formed will be engulfed by the NN PBH host halo. Conversely,
if $f_{\rm PBH}\ll 10^{-4}$, the existing halos will be much less disrupted,
which facilitates the formation of more massive PDM halos through mergers with
other nearby PDM halos. We illustrate our model in Figure 8, showing different
evolutionary possibilities.

We again emphasize that this is an extremely simplified dynamical model,
serving to illustrate the key governing physics in the complexity of our
simulations, where PBHs interact with standard PDM.

##  5 Summary and Conclusions

In this paper, we investigate how the presence of million-solar-mass PBHs
impact the assembly history of DM halos in the early universe. Building on our
earlier work to implement PBH initial conditions (Liu et al., 2022; Liu &
Bromm, 2023), we carry out DM-only cosmological simulations to study the
impact of such $10^{6}\,\mathrm{M}_{\odot}$ PBHs on the HMF and early
structure formation. Our primary findings indicate that PBHs can seed the
formation of massive halos. When PBHs constitute a fraction $f_{\rm PBH}\gg
10^{-4}$ of DM, the PBH-seeded halos efficiently disrupt and swallow their
pure PDM-halo counterparts, leading to a bimodal feature in halo mass
distribution. Additionally, if $f_{\rm PBH}\gg 10^{-3}$, nearly all halos will
contain at least one PBH, with PDM halos being engulfed by PBH-seeded halos
shortly after their formation. When varying the initial conditions in terms of
perturbation strength while maintaining the same PBH mass abundance, we find
that these changes have non-trivial effects, but do not affect the bimodal
feature in the HMF, compared to simulations that vary $f_{\rm PBH}$.

In modeling this halo formation scenario, we employ NN statistics to observe
the clustering patterns of neighboring halo pairs and their evolution across
redshift. To facilitate the physical interpretation, we develop a halo free-
fall model to elucidate the evolution of structure and clustering patterns
observed in the NN analysis. We identify a critical distance below which PDM
halo formation is impacted, and any already formed PDM halo located much
closer to its NN PBH-seeded halo than this distance will be engulfed.
Furthermore, the critical distance can be compared with the inter-PBH
separation distance $d_{\rm PBH}$, which gives a critical value for the PBH
DM-fraction of $f_{\rm PBH}\simeq 10^{-4}$, to ascertain whether we are
entering a universe dominated by PBH-seeded halos.

Regarding structure formation in this PBH-CDM setting, we use the Press-
Schechter (PS) formalism (Press & Schechter, 1974; Sheth & Tormen, 1999) as an
approximate avenue towards deriving the HMF analytically, assuming a Gaussian
distribution of overdensity peaks. However, the overdensity from PBHs, due to
their discreteness, takes the form of a delta function in real space, and the
PBHs are Poisson distributed at small scales. These features evade the
Gaussian assumption, so PBHs and their perturbations contribute non-trivially
to the mass variance when taking the volume average over the window function
(Kashlinsky, 2016; Desjacques & Riotto, 2018; De Luca et al., 2020). The
standard PS formalism may thus not adequately describe structure formation or
the HMF, due to the deviation from Gaussian statistics in PBH cosmologies. One
consequence is that the PS formalism does not sufficiently explain the
predicted bimodal feature in the HMF. Future work, by carefully revisiting the
evaluation of the filtered-mass variance $\sigma(M)$ and the nonlinear growth
of structure around PBHs, needs to address these shortcomings.

We recognize several limitations in our study that point to future work on
massive PBHs. Firstly, our simulations are DM-only, which are suitable for
studying large-scale structures, but they cannot accurately model the
structure on the scale of individual halos. Moreover, without incorporating
baryonic physics (e.g., baryon-DM streaming motion, see Kashlinsky, 2021;
Atrio-Barandela, 2022; Liu et al., 2022) in our simulations, it is challenging
to model the growth of PBHs and their feedback effects on galaxies (see, e.g.,
Pandey & Mangalam, 2018; Lu et al., 2021; Takhistov et al., 2022; De Luca &
Bellomo, 2023; Casanueva-Villarreal et al., 2024). Additionally, with the
current simulation volume, we do not produce any DM halos with masses $\gg
10^{9}\,\mathrm{M}_{\odot}$ at $z\sim 10$. To compare our model with JWST
observations of high-redshift ($z\simeq 10$) galaxies (e.g., Finkelstein et
al., 2022; Castellano et al., 2023), we need to develop a mock galaxy
catalogue over a much larger cosmic volume, with linear size of $L\gtrsim 10\
h^{-1}\rm Mpc$, capable of capturing the formation of more massive DM halos
within the simulation box. With a larger simulation volume, we could also
extend our study of the structure formation to lower redshifts ($z\simeq 1$),
and use the NN statistics to test against sky surveys like SDSS-V (Kollmeier
et al., 2019) to better constrain the PBH abundance.

We would like to acknowledge the fruitful discussion with Nicola Bellomo,
Steve Finkelstein and Florian Kühnel. BL is supported by the Royal Society
University Research Fellowship. The authors acknowledge the Texas Advanced
Computing Center (TACC) for providing HPC resources under allocation AST23026.

##  Appendix A Cosmic Web and Nearest Neighbor Statistics for other Runs

Figure 9: Structure formation for different PBH abundances. Similar to Fig. 2,
we project our simulation box onto a 2D plane to depict the cosmic web at
$z=10$ in comoving units. Left panel: Cosmic Web from the PBH2 simulation run,
with $f_{\rm PBH}=10^{-2}$. Right panel: Cosmic Web from the PBH4 simulation
run, with $f_{\rm PBH}=10^{-4}$. In both cases, all subhalos within their more
massive hosts have been removed to better showcase the spatial distribution of
halos and PBHs. Figure 10: NN statistics for the PBH_PS case, with a PBH mass
fraction of $f_{\rm PBH}=10^{-3}$. Different from Fig. 6, this simulation does
not consider non-linear perturbations around PBHs in the initial conditions.
Figure 11: NN statistics for large PBH-abundance case. Similar to Fig. 10, we
analyze the PBH2 run, with a PBH mass fraction of $f_{\rm PBH}=10^{-2}$. We
employ the same convention for the colors as before to differentiate halo
groups. Figure 12: NN statistics for small PBH abundance. Similar to Fig. 10,
we plot the same quantities for the PBH4 case, with a PBH mass fraction of
$f_{\rm PBH}=10^{-4}$. We again use the same color convention to differentiate
groups of halos. Note the convergence to the standard $\Lambda$CDM case, as
expected for increasingly small PBH abundances.

## References

  * Afshordi et al. (2003) Afshordi, N., McDonald, P., & Spergel, D. N. 2003, ApJ, 594, L71
  * Ali-Haïmoud et al. (2017) Ali-Haïmoud, Y., Kovetz, E. D., & Kamionkowski, M. 2017, Phys. Rev. D, 96, 123523
  * Arrabal Haro et al. (2023) Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, ApJ, 951, L22
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Atrio-Barandela (2022) Atrio-Barandela, F. 2022, ApJ, 939, 69
  * Banerjee & Abel (2021) Banerjee, A., & Abel, T. 2021, MNRAS, 500, 5479
  * Begelman et al. (2006) Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109
  * Belotsky et al. (2019) Belotsky, K. M., Dokuchaev, V. I., Eroshenko, Y. N., et al. 2019, Eur. Phys. J. C, 79, 246
  * Bogdán et al. (2024) Bogdán, Á., Goulding, A. D., Natarajan, P., et al. 2024, Nature Astronomy, 8, 126
  * Boylan-Kolchin (2023) Boylan-Kolchin, M. 2023, Nature Astronomy, 7, 731
  * Bromm & Loeb (2003) Bromm, V., & Loeb, A. 2003, ApJ, 596, 34
  * Brooks et al. (2013) Brooks, A. M., Kuhlen, M., Zolotov, A., & Hooper, D. 2013, ApJ, 765, 22
  * Cappelluti et al. (2022) Cappelluti, N., Hasinger, G., & Natarajan, P. 2022, ApJ, 926, 205
  * Carr et al. (2021a) Carr, B., Clesse, S., García-Bellido, J., & Kühnel, F. 2021a, Physics of the Dark Universe, 31, 100755
  * Carr et al. (2021b) Carr, B., Kohri, K., Sendouda, Y., & Yokoyama, J. 2021b, Reports on Progress in Physics, 84, 116902
  * Carr & Kühnel (2020) Carr, B., & Kühnel, F. 2020, Annual Review of Nuclear and Particle Science, 70, 355
  * Carr & Silk (2018) Carr, B., & Silk, J. 2018, MNRAS, 478, 3756
  * Carr & Sakellariadou (1999) Carr, B. J., & Sakellariadou, M. 1999, ApJ, 516, 195
  * Casanueva-Villarreal et al. (2024) Casanueva-Villarreal, C., Tissera, P. B., Padilla, N., et al. 2024, arXiv e-prints, arXiv:2405.02206
  * Castellano et al. (2022) Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJ, 938, L15
  * Castellano et al. (2023) —. 2023, ApJ, 948, L14
  * Colazo et al. (2024) Colazo, P. E., Stasyszyn, F., & Padilla, N. 2024, arXiv e-prints, arXiv:2404.13110
  * Curtis-Lake et al. (2023) Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, 7, 622
  * De Luca & Bellomo (2023) De Luca, V., & Bellomo, N. 2023, arXiv e-prints, arXiv:2312.14097
  * De Luca et al. (2020) De Luca, V., Desjacques, V., Franciolini, G., & Riotto, A. 2020, J. Cosmology Astropart. Phys, 2020, 028
  * Desjacques & Riotto (2018) Desjacques, V., & Riotto, A. 2018, Phys. Rev. D, 98, 123533
  * Diemer (2018) Diemer, B. 2018, ApJS, 239, 35
  * Escrivà (2022) Escrivà, A. 2022, Universe, 8, 66
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M. B., Haro, P. A., et al. 2022, ApJ, 940, L55
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13
  * Goulding et al. (2023) Goulding, A. D., Greene, J. E., Setton, D. J., et al. 2023, ApJ, 955, L24
  * Greene et al. (2024) Greene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39
  * Hahn & Abel (2011) Hahn, O., & Abel, T. 2011, MNRAS, 415, 2101
  * Hawking (1971) Hawking, S. 1971, Monthly Notices of the Royal Astronomical Society, 152, 75
  * Hooper et al. (2024) Hooper, D., Ireland, A., Krnjaic, G., & Stebbins, A. 2024, J. Cosmology Astropart. Phys, 2024, 021
  * Hopkins & Raives (2016) Hopkins, P. F., & Raives, M. J. 2016, MNRAS, 455, 51
  * Inayoshi et al. (2022) Inayoshi, K., Harikane, Y., Inoue, A. K., Li, W., & Ho, L. C. 2022, ApJ, 938, L10
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27
  * Inman & Ali-Haïmoud (2019) Inman, D., & Ali-Haïmoud, Y. 2019, Phys. Rev. D, 100, 083528
  * Jaacks et al. (2019) Jaacks, J., Finkelstein, S. L., & Bromm, V. 2019, MNRAS, 488, 2202
  * Jeon et al. (2024) Jeon, J., Bromm, V., Liu, B., & Finkelstein, S. L. 2024, arXiv e-prints, arXiv:2402.18773
  * Jeon et al. (2023) Jeon, J., Liu, B., Bromm, V., & Finkelstein, S. L. 2023, MNRAS, 524, 176
  * Juodžbalis et al. (2024) Juodžbalis, I., Maiolino, R., Baker, W. M., et al. 2024, arXiv e-prints, arXiv:2403.03872
  * Kashlinsky (2016) Kashlinsky, A. 2016, ApJ, 823, L25
  * Kashlinsky (2021) —. 2021, Phys. Rev. Lett., 126, 011101
  * Kawasaki & Murai (2019) Kawasaki, M., & Murai, K. 2019, Phys. Rev. D, 100, 103521
  * Kollmeier et al. (2019) Kollmeier, J., Anderson, S. F., Blanc, G. A., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 274
  * Kovacs et al. (2024) Kovacs, O. E., Bogdan, A., Natarajan, P., et al. 2024, arXiv e-prints, arXiv:2403.14745
  * Labbé et al. (2023) Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, 616, 266
  * Larson et al. (2023) Larson, R. L., Finkelstein, S. L., Kocevski, D. D., et al. 2023, ApJ, 953, L29
  * Liu & Bromm (2022) Liu, B., & Bromm, V. 2022, ApJ, 937, L30
  * Liu & Bromm (2023) —. 2023, arXiv e-prints, arXiv:2312.04085
  * Liu et al. (2022) Liu, B., Zhang, S., & Bromm, V. 2022, MNRAS, 514, 2376
  * Lodato & Natarajan (2006) Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813
  * Lu et al. (2021) Lu, P., Takhistov, V., Gelmini, G. B., et al. 2021, ApJ, 908, L23
  * Mack et al. (2007) Mack, K. J., Ostriker, J. P., & Ricotti, M. 2007, ApJ, 665, 1277
  * Maiolino et al. (2023) Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2308.01230
  * Maiolino et al. (2024) Maiolino, R., Scholtz, J., Witstok, J., et al. 2024, Nature, 627, 59
  * Meszaros (1975) Meszaros, P. 1975, A&A, 38, 5
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution
  * Naidu et al. (2022) Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14
  * Nakama et al. (2018) Nakama, T., Carr, B., & Silk, J. 2018, Phys. Rev. D, 97, 043525
  * Natarajan et al. (2024) Natarajan, P., Pacucci, F., Ricarte, A., et al. 2024, ApJ, 960, L1
  * Pandey & Mangalam (2018) Pandey, K. L., & Mangalam, A. 2018, Journal of Astrophysics and Astronomy, 39, 9
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425
  * Ricotti (2007) Ricotti, M. 2007, ApJ, 662, 53
  * Ricotti et al. (2008) Ricotti, M., Ostriker, J. P., & Mack, K. J. 2008, ApJ, 680, 829
  * Robertson et al. (2023) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, Nature Astronomy, 7, 611
  * Sheth & Tormen (1999) Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119
  * Shu (1982) Shu, F. H. 1982, The Physical Universe
  * Silk et al. (2024) Silk, J., Begelman, M. C., Norman, C., Nusser, A., & Wyse, R. F. G. 2024, ApJ, 961, L39
  * Smith & Bromm (2019) Smith, A., & Bromm, V. 2019, Contemporary Physics, 60, 111
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  * Su et al. (2023) Su, B.-Y., Li, N., & Feng, L. 2023, arXiv e-prints, arXiv:2306.05364
  * Tacchella et al. (2023) Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023, ApJ, 952, 74
  * Takhistov et al. (2022) Takhistov, V., Lu, P., Gelmini, G. B., et al. 2022, J. Cosmology Astropart. Phys, 2022, 017
  * Trenti & Stiavelli (2009) Trenti, M., & Stiavelli, M. 2009, ApJ, 694, 879
  * Wang et al. (2024) Wang, B., Leja, J., Atek, H., et al. 2024, ApJ, 963, 74
  * Wetzel et al. (2016) Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJ, 827, L23
  * Yan et al. (2023a) Yan, H., Ma, Z., Ling, C., Cheng, C., & Huang, J.-S. 2023a, ApJ, 942, L9
  * Yan et al. (2023b) Yan, H., Sun, B., Ma, Z., & Ling, C. 2023b, arXiv e-prints, arXiv:2311.15121
  * Zel’dovich (1970) Zel’dovich, Y. B. 1970, A&A, 5, 84
  * Zel’dovich & Novikov (1967) Zel’dovich, Y. B., & Novikov, I. D. 1967, Soviet Ast., 10, 602
  * Zhang et al. (2024) Zhang, S., Liu, B., & Bromm, V. 2024, MNRAS, 528, 180
  * Ziparo et al. (2022) Ziparo, F., Gallerani, S., Ferrara, A., & Vito, F. 2022, MNRAS, 517, 1086


# Design Principles and Physical Properties of Two-Dimensional
Heterostructured Borides

Alejandro Lopez-Bezanilla<EMAIL_ADDRESS>Theoretical Division, Los
Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

Principles of design to create dynamically stable transition metal,
lanthanide, and actinide based low-dimensional borides are presented. A charge
transfer analysis of donor metal atoms to electron deficient honeycombed B
lattices allows to predict complex covalent heterostructures hosting Dirac
states. The applicable guidelines are supported with the analysis of phonon
spectra computed with first-principles calculations to demonstrate the
physical stability of nanometer-thick heterostructures. Similar or dissimilar
layered borides can be stacked on top of each other in a layer-by-layer
fashion creating an interface that can be fundamentally different from the
individual layers, opening a rich playground to explore novel physical
properties and new materials. Functionalities such as multiple Dirac states,
highly dispersive electronic bands, and decoupled acoustic-optical phonon are
studied. The combination of appealing electronic properties and physical
realization make of predicted layered borides promising materials to integrate
a new generation of two-dimensional materials.

###### pacs:

Valid PACS appear here

##  I Introduction

Progress on synthesis routes, ranging from vapor deposition to solution-
processing Lewis  _et al._ (2017), will allow for expanding the range of
devices that build on the advantages of two-dimensional (2D) materials, which
will bring enormous benefits for technology industryFerrari  _et al._ (2015).
Special efforts are devoted to developing exfoliation and transfer techniques
which allow for piling-up several types of layered materials forming
heterostructures kept together by van der Waals (vdW) forces. With strong
covalent bonding between the atoms within each layer, layered compounds such
as NbSe2, MoS2, WTe2, and TaS2, to name a few, share structural similarity but
exhibit very different electronic properties, ranging from semiconducting
Mattheiss (1973) to metallic Ayari  _et al._ (2007) behaviour depending on
composition, geometry, and thickness. The resulting thick slabs own novel or
enhanced properties with electronic structures often distinctive from their
layered constituentsBritnell  _et al._ (2013); Hamm and Hess (2013). The
broadness in functionality observed on both the monolayers and the
heterostructures originates from the quantum confinement of electrons to the
2D geometry, and also from the modulation of the chemical composition and
electronic structure at the atomic scale.

Strictly within the flatland, it remains to be explored the role of strong
charge transfer (covalent bonding) between stacked layers, as opposed to
transient distortions in the distribution of charge density between layers as
in vdW heterostructures. Solid candidates to form all-covalent 2D materials
are borides, whose parent materials exist in bulk form. Borides are widely
used as the world’s hardest metals, and to store actinide elements in stable
compoundsLupinetti  _et al._ (2002). Layered borides can be prepared at
ambient pressureJames and Jasuja (2017) and, due to their high electrical
conductivity, can be easily cleaved and shaped with electric discharge
devices, an extended tool in metal manufacturing. Conventional methods for the
syntheses of metal borides are based on the high-temperature
carbothermal/borothermal reduction of metal oxides with boron-carbon sources
such as boron oxide/carbon or boron carbideGreenwood  _et al._ (1966). One
might assume intrinsic experimental difficulties to isolate atomic boride
monolayers since they cannot slide past each other due to the strong
interlayer covalent bonding, but experimental achievements in the exfoliation
of metal borides into aqueous dispersion of few-layer-thick structures allows
for envisioning purely covalent layered structures James and Jasuja (2017).

Figure 1: Schematic representation of the multilayered boride systems. Several
honeycombed layers of boron atoms (B) sandwich triangular layers of $d$\- and
$f$-valence orbital based metal elements (Mn).

In this paper, we study metal-atom based borides in 2D to analyze their
electronic, thermal and dynamical properties, expanding the scope of borides
beyond storage and metallurgical purposes. We branch away from the abundant
computational approaches that have propelled the field of materials discovery
by seeking optimal configurations relying on little or no empirical
knowledgeOganov  _et al._ (2019). Instead, we propose an insightful analysis
driven by an intuition-guided materials informatics scheme to search for
structure-property relationships on complex sets of borides. As a result, a
broad range of dynamically stable 2D materials hosting multiple cone-shaped
bands (Dirac states) are identified. In particular, a principle of design to
build $d$\- and $f$-orbital-based element borides by piling up layers of
existing borides with heterogeneous compositions is procured. A comparison of
band and phonon diagrams as well as heat capacity is provided for a series of
isoelectronic compounds. The physical viability of the predicted 2D compounds
is demonstrated within the harmonic approximation by phonon spectrum analysis.

##  II Computational Methods

Correlated effects of $f$-orbitals are specially difficult to capture with
current first-principles methods, and multiple schemes with different models
of screening can be used to provide a description of the band structure of
materials based on lanthanide and actinide elements. The approach of choice
for this study was self-consistent density functional theory (DFT) with
relativistic calculations including spin-orbit coupling (SOC) for all
materials. Even though including SOC in a DFT framework may pose some accuracy
issues Nordström  _et al._ (2000); Söderlind  _et al._ (2019), band structure
calculations exhibit expected results such as small band gap openings and
mixing of electronic states that remove degeneracies. DFT-based calculations
were performed within the Perdew-Burke-Ernzerhof (PBE) generalized gradient
approximation (GGA) functional for the exchange correlation, and the
projector-augmented-wave method as implemented in VASPKresse and Hafner
(1993); Kresse and Furthmüller (1996); Kresse and Joubert (1999). The
electronic wave functions were computed with plane waves up to a kinetic-
energy cutoff of 500 eV. The integration in the 2D k-space was performed using
a 56$\times$56$\times$1 $\Gamma$-centered k-point grid. Atomic coordinates and
lattice constants were fully relaxed until the residual forces were smaller
than 1 meV/Å. Including an on-site Coulomb repulsion term on the U atom leads
to artificial magnetic moments that are absent in experimental observations of
for example paramagnetic bulk UB2 Ohkochi  _et al._ (2008) and, therefore, it
was excluded in this study. The force-constant method and the PHONOPY package
Togo and Tanaka (2015) were employed for computing phonon spectra and heat
capacities.

##  III Design principles

Boron atom is the first and lightest element of group III characterized by an
atomic number $Z$=5 and an $s^{2}p^{1}$ electronic configuration of its three
valence electrons. Boron valence orbitals can combine to form $sp^{2}$ hybrid
orbitals, yielding strong and highly directional $\sigma$ bonds which,
together with the short covalent radius and large coordination numbers, allows
boron to form a large amount of stable clustered speciesBoustani  _et al._
(1999). The nature of the bonding in metal borides has been discussed from the
point of view of gain or loss of electrons by the B atoms relative to the
metal atoms. There is a consensus that electron transfer must be in the same
direction for all types of B-rich borides Greenwood  _et al._ (1966),
following a simple rule based on the relative electronegativities of B and the
d- or f-orbital atoms.

The so-called ’Aufbau principle’Boustani (1997) provides construction rules
for the design of different highly stable boron structures with quasiplanar
geometries that can be viewed as finite-size clusters of planar surfaces. The
lack of the additional stability provided by the hyper-conjugated network of
$p_{z}$ orbitals typical of C based materials such as graphene or carbon
nanotubes is the result of B atom electron deficiency (one valence electron
less than C atom) which prevents boron from exhibiting stand-alone honeycombed
2D structures. Although periodic flat distributions including triangular
motifs and vacanciesFeng  _et al._ (2016) are possible, realization of 2D
extended flat structures composed of B atoms arranged in honeycombed lattices
is always subjected to additional electron supply from another source such as
a metal substrateLi  _et al._ (2018).

The general design rules introduced by Tang and Ismail-BeigiTang and Ismail-
Beigi (2007, 2009) clarified that an optimal configuration is obtained if
electrons fill all the in-plane $\sigma$ bonding states leaving in-plane
antibonding states empty with partial occupancy of the $\pi$ solid-state
orbitals. Self-doping was demonstrated as a possible route to stabilizing 2D
and tubular all-B structures. Detailed explanations were providedTang and
Ismail-Beigi (2007, 2009) on the method to obtain 2D metal borides based on a
charge transfer model in which a metal atom donates electrons to the boron
sublattice, such as layered MgB2.

Following this metal-boron stabilizing interaction, we demonstrated the
viability of 2D ZrB2Lopez-Bezanilla (2018a), CrB4Lopez-Bezanilla (2018b), and
UB4Lopez-Bezanilla (2020a) and the physical properties derived from the
interplay between $p$ and $d$ orbitals. Continuing along these lines, here we
use and extend the charge transfer model of Refs.Tang and Ismail-Beigi (2007,
2009); Liu and Tománek (2019) as a general design principle for boron
heterostructures with dynamical stability and relevant physical properties
based on metal-atom doping. Specifically, we demonstrate with first-principles
calculations that $d$-orbital (transition metal) and $f$-orbital (lanthanide
and actinide) atoms yield the exact amount of charge to several honeycombed
boron lattices so that 2D structures composed of several metal and boron atoms
layers are physical viable. This results is consistent with nuclear
magnetic‐resonance experiment of the groups IV, V, and VI transition‐metal
bulk diboridesSilver and Kushida (1963) that suggested that the stability of
$\sigma$ and $\pi$ bonds in the boron layer is possible by transfer of
electrons in metal-boron donor-acceptor bonds.

In a single B honeycombed lattice with a metal atom (M) sitting on top of each
hexagon, transfer of a whole electron per B atom would give a formulation
M2+(B-)2. B layers can be regarded as isoelectronic with graphene upon charge
transfer, which makes Ti, Zr, Hf based monolayers isostructural transition
metal borides stable compoundsLopez-Bezanilla (2018a). The resulting material
features electronic properties similar to those of graphene such as dispersive
Dirac states in the vicinity of the Fermi energy level, with different
equivalency and slight mixing of states due to SOC effects.

Figure 2:  Phonon spectra of studied layered borides. The absence of negative
frequencies guarantees the dynamical stability of the nanostructures. The same
y-range for the single-metal-atom type borides allows to appreciate the effect
of an increasing atomic mass in reducing the phonon dispersion and compacting
the total frequency range. In most cases, one or two indirect band gaps are
observed between acoustic and optical branches (i.e., UTi6) or optical
branches (i.e., NbTaB6). The presence of phonon gaps has a direct consequence
on the evolution of the heat capacity with temperature, as shown in Figure 4.

Deposition of a second hexagonal B layer on top of the triangular lattice of M
atoms yields a structure where each M atom is covalently attached to four B
atoms simultaneously, two in each layer. The formulation M4+(B-)4 entails the
utilization of d4 of f4 atoms as stabilizing species of both B lattices.
Candidate M atoms are Cr, Mo, W in the d-orbital series, and Nd and U atoms in
the f-orbital series.

Following this metal-boron stabilizing interaction, it is natural to consider
complex structures where a M4+(B-)4 structure serves as a substrate of
additional single layer metal borides. Indeed, over that closed shell compound
a M2+(B-)2 layer could be deposed, conserving the heterostructure the
electronic balance of the parent compounds. According to the representation of
Figure1, the new compound would be represented by M1M2B6. Furthermore, adding
another layer to the opposite side of the M1B4 sub-layer would still preserve
the charge balance of the new M1M2M3B8 heterostructure. The theoretical
prediction of the viability of these heterostructures does not entail a
trivial experimental synthesis but provides evidence of the potential richness
of borides when combined to form multi-element 2D compounds. Their dynamical
stability is analyzed in the next section to show that despite the disparity
of the B and metal-atom masses and radii, the stability is guaranteed as long
as the electronic balance is preserved.

##  IV Dynamical stability

Metal borides as well as boron hydrides are some of the few systems in which
extensive clustering and catenation of boron atoms occursGreenwood  _et al._
(1966). In our case, the dynamical stability of the 2D transitional metal
borides proposed above is practically straightforward if we consider that some
of the single metal-atom compounds analyzed exist in bulk form. For instance,
vertical AA stacking of single hexagonal B layers intercalated with Zr or U
yields bulk ZrB$2$Zimmermann  _et al._ (2008) and UB$2$Kardoulaki  _et al._
(2020) respectively, two ceramic materials with a hexagonal covalent structure
owning good thermal and electrical conductivities.

Figure 3: Separate contribution of each atom to the phonon frequencies,
represented as projected density of phonon states, of TiUTiB8, Nb2B6, and
VTiVB8. Heaviest metal atoms contribute mainly to lowest frequency phonons
whereas B atoms dominate the higher frequency modes. y-axis is in arbitrary
units.

Dynamical stability is analyzed within the linear displacement regime at zero
temperature for predicting vibration properties. Inspection of phonon diagrams
in Figure 2 allow us to verify that the compounds are stable as the absence of
negative vibration frequencies and of anomalous softening modes reveal. Since
the analysis of the Zr, Cr, and U single layer borides were presented
elsewhere by the author of this paperLopez-Bezanilla (2018a, b, 2020b), we
will focus hereafter on the new proposed borides. Phonon diagrams of
heterostructures are indicated by the notation M1M2B6 and M1M2M3B8, where in
some cases M1=M2 (e.g. V2B6) or M2=M3 (e.g. ZrMoZrB8). The Brillouin zone (BZ)
of interest with symmetry P6/mm is the area delimited by the $\Gamma$, $M$,
and $K$ points in the hexagonal reciprocal lattice.

In all cases, restoring forces oppose deformation when atoms are displaced
from their equilibrium positions, indicating that small deformations induced
on the material are inefficient in changing their shape. Phonon modes span a
frequency range of less than 30 THz, depending the total extension on the mass
of the metal atom: the heavier the metal atom(s), the narrower the spectrum
extension. Bi- and tri-metallic atom layered structures exhibit phonon diagram
with 24 and 33 phonon branches respectively: three acoustic (LA, TA and ZA)
and 21 or 30 optical modes. The LA and TA modes correspond to longitudinal and
transverse phonon oscillations of the metallic atoms in the 2D layer plane.
Mode ZA corresponds to vibrations in the direction normal to the direction of
oscillation of the former modes. LA and TA modes progress from a linear
dispersion in the zone center ($\Gamma$) to a non-dispersive evolution in the
zone edge ($M-K$). Contrary to commonly observed in 2D nanostructures Molina-
Sánchez and Wirtz (2011), layered borides exhibit optical branches of larger
dispersion and broader distribution of phonon velocities than acoustical
branches. For instance, Ta2B6 exhibits 10 THz dispersive optical modes,
whereas the acoustic modes span $\sim$6 THz. In the single-type metal-atom
layers, the heavier the metal atom is, the less dispersive the acoustic
phonons are. The $d^{3}$ and $f^{3}$ atom series in the first row of Figure 2
clearly displays the progressive narrowing of the frequency range as one
increases the number of protons in the isomorphic structures.

The large difference of masses between B and M atoms determines the
predominant contribution of each atom to acoustic and optical phonon branches.
Figure 3 shows the separate contribution of each atom to the vibration modes
at all frequencies of three borides. The heavier metal atoms contribute mainly
to the lowest frequencies, an in Nb2B6, were all phonons branches below 8 THz
correspond almost exclusively to vibration modes of the Nb atoms. In UTi2B8
the U atom participates only in the acoustic modes. In turn, Ti atoms are the
main participants in the next set of six optical branches. Although B atoms
can vibrate at all frequencies, their energetic contribution to the vibration
spectrum is mainly above 10 THz. For two metal atoms of similar masses, the
difference is less pronounced but significant, as in VTiVB8, where both V and
Ti peaks in the projected density of states differ by 2.3 THz. Visual
inspection of the eigenvectorshen  generated by PHONOPY reveals that at lower
frequencies metal atoms vibrate with larger amplitude than B atoms. Flat
branches right above the acoustics correspond to vertical oscillations of M
atoms within the honeycomb layers, which explains the low propagation
character of those branches. At higher frequencies B atoms vibrate with a
remarkably large amplitude while M atoms are barely displaced.

The large difference of the atomic masses is the origin of an indirect gap
between optical and acoustic modes in some structures, similarly to the gap
observed in bulk BAsLindsay  _et al._ (2012, 2013); Tian  _et al._ (2018) and
Li2TeMukhopadhyay  _et al._ (2016). This gap and that gap found in the rest of
the structures between bunched optical branches enables selective phonon
frequency transmission, namely no kind of atomic vibration propagation is
allowed in the band gap. In particular UTiB6 displays a tiny gap (which
vanishes when adding a second layer of TiB2) that would the acoustic modes
avoid crossing with optical modes. This may potentially lead to interesting
phenomena in thermal transport, since acoustic-optical phonon backscattering
is not allowed by the gapLindsay  _et al._ (2012).

Figure 4: Specific heat capacities $c$ of the multilayered B compounds
normalized to the mass in each unit cell are compared with the one of
graphene. Inset: The derivative of $c$ with respect to the temperature shows
that $c$ increases differently at several T ranges as a consequence of the
various gaps observed in the phonon spectra.

Figure 4 shows the dependency of the specific heat capacity $c$, in
$JK^{-1}mol^{-1}$ and normalized to the atomic weight $m$ of each cell, with
temperature of all the borides. Specific heat capacity is defined as the ratio
of the supplied heat to the resulting temperature change for a unit mass of
the substance and, therefore it measures the ability of the materials to
absorb heat. Dimensional restriction, small group velocity, and large mass of
the metallic atoms lead to high specific heat capacities. Lighter metal based
borides exhibit enhanced $c$ at any temperature with respect to the heavier-
metal compounds, exceeding in all cases that of graphene below 600 K, also
shown for comparison in Figure 4. Each bundled curves shows how much of an
effect metal atoms can have on the heat capacity of the borides, as the main
difference at temperatures above $\sim$ 100 K is in the composition. At
$\sim$100 K the normalized $c$ of all borides are very similar and much larger
than that of graphene, whereas at $\sim$ 200 K the $c/m$ of lighter compounds
starts diverging. At large temperature, the specific heat capacities serves to
differentiate the mass composition of the layered borides.

At very low temperatures, the heat capacity of some borides undergo two slope
changes. For a clearer representation, in the inset of Figure 4 the derivative
of $c/m$ with respect to the temperature is plotted. With two local maxima,
first at $\sim$33 K and second at $\sim$130 K, this anomalous change is due to
the frequency gap between phonon branches that can also be observed for
compounds such as Pr2B6 and UTiB6, and is absent in the graphene layer.

##  V Electronic properties

Figure 5: Electronic band diagram of V2B6, Nb2B6, Ta2B6, VTaB6, MoZr2B8. Two
Dirac cones in the vicinity of the Fermi level develop a gap as the spin-orbit
coupling effects become larger with metallic atoms size. Very dissimilar
(VTaB6) and similar atomic radius (MoZr2B8) lead to a very different
distribution of electronic bands despite the charge transfer model is the same
for all the borides. Color-resolved electronic state projection (lower panels)
of Ta2B6 show the participation of each atomic orbital to the material’s
electronic states. The color intensity of the lines is proportional to the
contribution. Horizontal dashed lines point out the Fermi level.

Low-energy fermionic excitations behaving as massless Dirac particles is a
subject of intense study for multiple technological applications ranging from
quantum information technologies to novel electronic devices.
GrapheneNovoselov  _et al._ (2004) is one the prototypical materials hosting
electronic bands with linear dispersion in which charge carriers travel at the
speed of light due to their zero massNovoselov  _et al._ (2005). The
electronic properties of the heterostuctured layered borides studied in this
paper resemble to monolayer graphene in the shape of their electronic bands
and the presence of nodes in the excitation spectrum as a result of groups of
two cone-shaped bands touching each other at one single point in the BZ. A
relevant difference of the borides with respect to graphene is the number of
such cones and the location in the BZ. Whereas graphene exhibits two Dirac
points at K and K’ points, borides outnumber with some electronic bands out of
high-symmetry points in the hexagonal BZ. The field of 2D Dirac materials
owning cone-shaped band structure and high carrier mobilities is expanded here
with the description of multiple borides that with strong or weak SOC exhibit
similar features to that of graphene.

Figure 5 shows the relativistic band structure of a selection of the borides.
Focus is on the lighter metal atom borides with no $f$-orbitals, which makes
DFT calculations less controversial. Indeed, multiple schemes with different
models of screening can be used to provide a first-principles description of
the band structure of correlated $f$-orbitals rich borides, they all leading
to various interpretations of results that as of today cannot be compared with
experimental evidence on thin metal borides. Relativistic effects due to the
large size of metallic elements are accounted in the calculations including
SOC effects.

First three panels of Figure 5 show the distribution of electronic states of
the group-V elements with $d^{3}$ electrons. According to our charge transfer
model, a total of six $d$-electrons are yielded to the six empty $p_{z}$
orbitals of the B unit-cell atoms in the three honeycombed lattices. As
opposed to graphene, several bands cross the Fermi level. Similarly to the
all-C monolayer a Dirac point is observed at K point, although slightly below
the Fermi level. In addition, a second set of cone-shaped bands meet at a
single point at Fermi level in the $\Gamma-M$ line. This picture is clearly
visible in V2B6. As the mass of the metal atom increases, a small gap between
Dirac cones develops, becoming a few meV large in Nb2B6, and clearly
noticeable in the band diagram of Ta2B6. The combination of two metal atoms
with dissimilar size and a relative lack of symmetry in the material leads to
a splitting of electronic bands and removal of state degeneracy, as shown in
VTaB6. However, a heterostructure composed of two similar-size atoms in a more
symmetric geometry such as MoZr2B8, partially preserve the dispersion of the
bands found in the one-metal-atom type borides.

It is worth determining the origin of the Dirac cones, and Ta2B6 is taken as a
case of study. Lower panels of Figure 5 display a color-scheme representation
of the electronic state diagrams projected onto each band, which allows for
observing the independent contribution of each M-$d$ and B-$p_{z}$ atomic
orbital. At K point, the main contribution to the upper cone-shaped bands is
from Ta-$d_{yz}$ orbitals, whereas Ta-$d_{xz}$ orbitals are involved in the
lower cone. Two different Ta orbitals are responsible of the two cones in the
middle of the high-symmetry $\Gamma-M$ line, namely Ta-$d_{xy}$ and
Ta-$d_{z^{2}}$. Note that whereas the cones at K point can be projected on two
cone-shaped contributions, the second two cones are projected in a state that
develops from near the zone center in the conduction bands down to the
vicinity of K (Ta-$d_{xy}$), and a state that extend from the valence band up
to high energies, with two gaps (Ta-$d_{z^{2}}$) in the energy range
displayed. From the B atom, the $p_{z}$ orbital also contribute to the
delocalization of charge projecting its contribution to the four cones,
pointing out to strong hybridization of this orbital with all Ta-atom
orbitals.

##  VI Conclusions

A general strategy to mix and match layered boride materials with different
properties has been presented. Covalent heterostructures with a variety of
dynamical and electronic properties have been proposed benefiting of the
ability of 2D layer of similar and dissimilar composition to attached each
other. Theoretical evidence has been provided on the following: i) a series of
borides constructed following a charge transfer model from metal to B atoms
exhibit dynamical stability as 2D compounds, ii) an itinerant character of the
charge carriers is observed due to the B honeycombed layers are isolectronic
to graphene, hosting electronic states with a delocalized character, iii) the
heat capacity of the thin layers is the largest observed in 2D materials. iv)
borides feature an infrequent acoustic-optical and optical-optical gap as a
result of the large difference in the mass between the metal and B atoms.

The large variety of compounds in which layered borides can take form leads to
a rich diversity of physical phenomena usually not found in 2D materials, such
as multiple Dirac states and bundling of acoustic and optical phonon branches.
Although the most accurate predictions must come upon contrast with
experimental observation, first-principles calculations support our charge
transfer model, providing a theoretical ground to predict new compounds. The
design principles of this study contribute to the creation of a theoretical
framework to predict and explain the physics of borides in 2D, where $d$ and
$f$-orbitals have a prominent role in the formation of Dirac states. The
predicted outstanding electronic properties of heterostructured borides
provide an extra motivation for borides integration with the nanometer-scale
silicon complementary metal-oxide-semiconductor (CMOS) technology as well as
beyond-CMOS devices and circuits. Zhang  _et al._ (2019)

##  VII Acknowledgments

Los Alamos National Laboratory is managed by Triad National Security, LLC, for
the National Nuclear Security Administration of the U.S. Department of Energy
under Contract No. 89233218CNA000001. I acknowledge the computing resources
provided on Bebop, the high-performance computing clusters operated by the
Laboratory Computing Resource Center at Argonne National Laboratory.

## References

  * Lewis  _et al._ (2017) E. A. Lewis, J. R. Brent, B. Derby, S. J. Haigh, and D. J. Lewis, Chem. Commun. 53, 1445 (2017).
  * Ferrari  _et al._ (2015) A. C. Ferrari, F. Bonaccorso, V. Fal’ko, K. S. Novoselov, S. Roche, P. Bøggild, S. Borini, F. H. L. Koppens, V. Palermo, N. Pugno, J. A. Garrido, R. Sordan, A. Bianco, L. Ballerini, M. Prato, E. Lidorikis, J. Kivioja, C. Marinelli, T. Ryhänen, A. Morpurgo, J. N. Coleman, V. Nicolosi, L. Colombo, A. Fert, M. Garcia-Hernandez, A. Bachtold, G. F. Schneider, F. Guinea, C. Dekker, M. Barbone, Z. Sun, C. Galiotis, A. N. Grigorenko, G. Konstantatos, A. Kis, M. Katsnelson, L. Vandersypen, A. Loiseau, V. Morandi, D. Neumaier, E. Treossi, V. Pellegrini, M. Polini, A. Tredicucci, G. M. Williams, B. Hee Hong, J.-H. Ahn, J. Min Kim, H. Zirath, B. J. van Wees, H. van der Zant, L. Occhipinti, A. Di Matteo, I. A. Kinloch, T. Seyller, E. Quesnel, X. Feng, K. Teo, N. Rupesinghe, P. Hakonen, S. R. T. Neil, Q. Tannock, T. Löfwander, and J. Kinaret, Nanoscale 7, 4598 (2015).
  * Mattheiss (1973) L. F. Mattheiss, Phys. Rev. B 8, 3719 (1973).
  * Ayari  _et al._ (2007) A. Ayari, E. Cobas, O. Ogundadegbe, and M. S. Fuhrer, Journal of Applied Physics 101, 014507 (2007), https://doi.org/10.1063/1.2407388 .
  * Britnell  _et al._ (2013) L. Britnell, R. M. Ribeiro, A. Eckmann, R. Jalil, B. D. Belle, A. Mishchenko, Y.-J. Kim, R. V. Gorbachev, T. Georgiou, S. V. Morozov, A. N. Grigorenko, A. K. Geim, C. Casiraghi, A. H. C. Neto, and K. S. Novoselov, Science 340, 1311 (2013), http://science.sciencemag.org/content/340/6138/1311.full.pdf .
  * Hamm and Hess (2013) J. M. Hamm and O. Hess, Science 340, 1298 (2013), http://science.sciencemag.org/content/340/6138/1298.full.pdf .
  * Lupinetti  _et al._ (2002) A. J. Lupinetti, J. L. Fife, E. Garcia, P. K. Dorhout, and K. D. Abney, Inorganic Chemistry 41, 2316 (2002), pMID: 11978092, https://doi.org/10.1021/ic015607a .
  * James and Jasuja (2017) A. L. James and K. Jasuja, RSC Adv. 7, 1905 (2017).
  * Greenwood  _et al._ (1966) N. N. Greenwood, R. V. Parish, and P. Thornton, Q. Rev. Chem. Soc. 20, 441 (1966).
  * Oganov  _et al._ (2019) A. R. Oganov, C. J. Pickard, Q. Zhu, and R. J. Needs, NATURE REVIEWS MATERIALS 4, 331 (2019).
  * Nordström  _et al._ (2000) L. Nordström, J. M. Wills, P. H. Andersson, P. Söderlind, and O. Eriksson, Phys. Rev. B 63, 035103 (2000).
  * Söderlind  _et al._ (2019) P. Söderlind, A. Landa, and B. Sadigh, Advances in Physics 68, 1 (2019), https://doi.org/10.1080/00018732.2019.1599554 .
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 48, 13115 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Ohkochi  _et al._ (2008) T. Ohkochi, S.-i. Fujimori, H. Yamagami, T. Okane, Y. Saitoh, A. Fujimori, Y. Haga, E. Yamamoto, and Y. Ōnuki, Phys. Rev. B 78, 165110 (2008).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Boustani  _et al._ (1999) I. Boustani, A. Rubio, and J. A. Alonso, Chemical Physics Letters 311, 21 (1999).
  * Boustani (1997) I. Boustani, Phys. Rev. B 55, 16426 (1997).
  * Feng  _et al._ (2016) B. Feng, J. Zhang, Q. Zhong, W. Li, S. Li, H. Li, P. Cheng, S. Meng, L. Chen, and K. Wu, Nat Chem 8, 563 (2016).
  * Li  _et al._ (2018) W. Li, L. Kong, C. Chen, J. Gou, S. Sheng, W. Zhang, H. Li, L. Chen, P. Cheng, and K. Wu, Science Bulletin 63, 282 (2018).
  * Tang and Ismail-Beigi (2007) H. Tang and S. Ismail-Beigi, Phys. Rev. Lett. 99, 115501 (2007).
  * Tang and Ismail-Beigi (2009) H. Tang and S. Ismail-Beigi, Phys. Rev. B 80, 134113 (2009).
  * Lopez-Bezanilla (2018a) A. Lopez-Bezanilla, Phys. Rev. Materials 2, 011002 (2018a).
  * Lopez-Bezanilla (2018b) A. Lopez-Bezanilla, 2D Materials 5, 035041 (2018b).
  * Lopez-Bezanilla (2020a) A. Lopez-Bezanilla, Journal of Physics: Materials 3, 024002 (2020a).
  * Liu and Tománek (2019) D. Liu and D. Tománek, Nano Letters 19, 1359 (2019).
  * Silver and Kushida (1963) A. H. Silver and T. Kushida, The Journal of Chemical Physics 38, 865 (1963), https://doi.org/10.1063/1.1733774 .
  * Zimmermann  _et al._ (2008) J. Zimmermann, G. Hilmas, W. Fahrenholtz, R. Dinwiddie, W. Porter, and H. Wang, Journal of the American Ceramic Society 91, 1405 (2008), cited By 206.
  * Kardoulaki  _et al._ (2020) E. Kardoulaki, J. White, D. Byler, D. Frazer, A. Shivprasad, T. Saleh, B. Gong, T. Yao, J. Lian, and K. McClellan, Journal of Alloys and Compounds 818, 153216 (2020).
  * Lopez-Bezanilla (2020b) A. Lopez-Bezanilla, Journal of Physics: Materials 3, 024002 (2020b).
  * Molina-Sánchez and Wirtz (2011) A. Molina-Sánchez and L. Wirtz, Phys. Rev. B 84, 155413 (2011).
  * (33) Http://henriquemiranda.github.io/phononwebsite/index.html.
  * Lindsay  _et al._ (2012) L. Lindsay, D. A. Broido, and T. L. Reinecke, Phys. Rev. Lett. 109, 095901 (2012).
  * Lindsay  _et al._ (2013) L. Lindsay, D. A. Broido, and T. L. Reinecke, Phys. Rev. Lett. 111, 025901 (2013).
  * Tian  _et al._ (2018) F. Tian, B. Song, X. Chen, N. K. Ravichandran, Y. Lv, K. Chen, S. Sullivan, J. Kim, Y. Zhou, T.-H. Liu, M. Goni, Z. Ding, J. Sun, G. A. G. Udalamatta Gamage, H. Sun, H. Ziyaee, S. Huyan, L. Deng, J. Zhou, A. J. Schmidt, S. Chen, C.-W. Chu, P. Y. Huang, D. Broido, L. Shi, G. Chen, and Z. Ren, Science 361, 582 (2018), https://science.sciencemag.org/content/361/6402/582.full.pdf .
  * Mukhopadhyay  _et al._ (2016) S. Mukhopadhyay, L. Lindsay, and D. S. Parker, Phys. Rev. B 93, 224301 (2016).
  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004), http://science.sciencemag.org/content/306/5696/666.full.pdf .
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * Zhang  _et al._ (2019) C. Zhang, T. He, S. K. Matta, T. Liao, L. Kou, Z. Chen, and A. Du, The Journal of Physical Chemistry Letters 10, 2567 (2019), https://doi.org/10.1021/acs.jpclett.9b00762 .