Primordial black hole formation with full numerical relativity

We study the formation of black holes from subhorizon and superhorizon perturbations in a matter dominated universe with 3+1D numerical relativity simulations. We find that there are two primary mechanisms of formation depending on the initial perturbation's mass and geometry -- via $\textit{direct collapse}$ of the initial overdensity and via $\textit{post-collapse accretion}$ of the ambient dark matter. In particular, for the latter case, the initial perturbation does not have to satisfy the hoop conjecture for a black hole to form. In both cases, the duration of the formation the process is around a Hubble time, and the initial mass of the black hole is $M_\mathrm{BH} \sim 10^{-2} H^{-1} M_\mathrm{Pl}^2$. Post formation, we find that the PBH undergoes rapid mass growth beyond the self-similar limit $M_\mathrm{BH}\propto H^{-1}$, at least initially. We argue that this implies that most of the final mass of the PBH is accreted from its ambient surroundings post formation.


INTRODUCTION
Primordial black holes (PBHs) form in the early stages of the universe, and their idea was first conceived in the late sixties and early seventies [1][2][3]. It is notable that it was the potential existence of small black holes from primordial origin that led Hawking to theorize black hole evaporation [4]. It was realised shortly after that PBHs could constitute a significant part of cold dark matter [5], and interest in PBHs has spiked in the recent past as a result. Evaporating PBHs have been suggested as explanations for galactic and extra-galactic γ-ray backgrounds, short γ-ray bursts and anti-matter in cosmic rays [6][7][8][9][10][11][12] and PBHs could provide seeds for the formation of supermassive black holes and large-scale structure [13,14]. Moreover, PBHs could be responsible for certain lensing events [15,16], with recent analysis suggesting that the population of black holes (BHs) detected by the LIGO/Virgo/KAGRA (LVK) observatories [17] may be primordial [18,19]. Additionally, work is underway to use next generation gravitational wave experiments to detect PBH formation and mergers [20,21]. Results obtained by the NANOGrav Collaboration [22] have been associated to PBHs, as well [23][24][25][26].
In the standard picture, these fluctuations collapse post inflation, while the universe is dominated by radiation energy. The nonzero radiation pressure resists collapse, meaning that the inhomogeneities must be fairly large for PBHs to form.
PBH formation in matter dominated epochs has also been extensively studied analytically and semianalytically. In various non-standard universe histories, inflation is followed by a period of matter domination [106][107][108][109]. PBH formation in such an early epoch of matter domination was considered early on [110]. More recently, a threshold amplitude for the collapse of a scalar field overdensity was found [111], the effects of nonsphericity [112] and inhomogeneity [113] on the collapse were investigated, the resulting spin of the PBHs was studied [114], the duration of an early epoch of matter domination was constrained by considering the PBH abundance [115] and constraints on the amplitude and spectral index of the collapsing scalar field were obtained [116].
In this work, we use full 3+1D numerical relativity simulations to investigate the collapse of subhorizon and superhorizon non-linear perturbations in an expanding universe that is dominated by matter. We model the expanding background and the collapsing perturbation using an oscillating massive scalar field and massless scalar field respectively. The massless scalar field's initial energy is thusly contained purely in its gradients. We will show that there are two broad mechanisms of black hole formation -via direct collapse for the case where the FIG. 1. Direct collapse and accretion driven mechanisms: The figure summarizes the two distinct processes of PBH formation studied in this paper. Top panel shows the direct collapse mechanism, where a high density initial superhorizon perturbation collapses and forms a black hole as soon as the perturbation reaches the centre. Bottom panel depicts the accretion driven collapse mechanism, where a lower density initial perturbation collapses and acts as a seed to trigger the accretion of the background dark matter, which subsequently collapses to form a black hole. Both start from the same initial radius R0, but with different initial amplitudes ∆ξ. In the leftmost figures, we show the initial size of the Hubble horizon (white solid line), which will grow as time evolves. In the other figures, the Hubble horizon has grown larger than the box size. Colourbars are shown in the top right, with lighter (darker) colours signifying higher (lower) energy densities, and scales fixed per mechanism. Video comparisons of these mechanisms can be found here [117] and here [118].
overdensity is sufficiently large that it will form a black hole, and via post-collapse accretion for the case where the overdensity is smaller. In both cases, the process is rapid and its duration is around a single Hubble time, forming PBHs with initial masses of M BH H ∼ 10 −2 M 2 Pl . We illustrate these mechanisms in Fig. 1.
Our choice of fundamental scalar fields as dark matter, instead of a pressureless cold fluid (as suggested by [119]), is prompted by our focus on an early time (i.e. pre-BBN) matter dominated phase instead of the present late time matter dominated phase. Such early era matter domination is often driven by a non-thermal fundamental scalar or moduli dynamics [120] instead of the more familiar cold pressureless fluid such as thermal WIMP dark matter. Furthermore, early matter phases will eventually transition into a radiation domination epoch, such that the standard Hot Big Bang cosmological evolution can proceed. Such a phase transition from matter domination into radiation domination can then be achieved through the decay of the scalar field into either standard model particles or intermediaries. This paper is organised as follows. In section I we explain the numerical setup we use to evolve a scalar perturbation in a dark matter dominated background. In section II we introduce the two aforementioned formation mechanisms and their characteristics, and we study the properties of the black holes that are formed post collapse. In section III we comment on the post formation growth of the PBHs, and we conclude in section IV.

I. EARLY MATTER DOMINATION EPOCH WITH SCALAR FIELDS
We will use a metric with − + ++ signature, in Planck units = c = 1 such that G = M −2 Pl . The action we will consider is involving a massive scalar field φ with mass m that models the ambient dark matter, and a massless scalar field ξ that sources the initial perturbation. They are both minimally coupled to gravity but not otherwise coupled to one another, i.e.
Since the field ξ has no potential, it will only influence dynamics via its gradients. Furthermore, it will dilute much more rapidly than dark matter, and hence it does not affect the long term dynamics of the system once its initial job of sourcing a perturbation is done 1 . When the gradients in ξ are negligible, the spacetime dynamics are dominated by the behaviour of the background scalar field φ. When φ is additionally homogeneous on a given spatial hyperslice, the metric of the spacetime is well described by the Friedman-Lemaître-Robertson-Walker (FLRW) line element where dΩ 2 2 = dθ 2 +sin 2 θdφ 2 . The scale factor a(t) evolves according to the Friedmann equation H 2 = 8πρ/3M 2 Pl , where H(t) ≡ȧ/a is the Hubble parameter 2 . The equation of motion for φ reduces to the Klein-Gordon equation where the Hubble parameter is and the corresponding pressure is given by If the oscillation of φ is sufficiently undamped, which is the case if 2m 3H 3 , the friction term in Eqn. (5) can be neglected. The dynamics of φ are then approximately given by a simple harmonic oscillator φ(t) = φ 0 cos mt , whose pressure is 1 In principle, we could use a single massive scalar φ. However, in practice, we find that large perturbations of the massive scalar would introduce a large infusion of potential energy into the dynamics of the background resulting in non-matter dominated evolution, at least initially. 2 Dotted variables are derivatives with respect to cosmic time t. 3 This condition is obtained by interpreting Eqn. (5) as the equation of motion for a damped oscillator of the form m dx 2 dt 2 + c dx dt + kx = 0. The condition for undamped oscillation then becomes As long as the oscillation period T is sufficiently smaller than one Hubble time, this averages to zero over one Hubble time, i.e. p DM = 0, resulting in a dark matter dominated expansion, which can be interpreted as a model for pressureless dust [121] at large scales.
Meanwhile, the massless scalar field ξ provides the energy density perturbation that will trigger BH formation. In this paper, we exclusively consider initially static spherically symmetric perturbations and we leave the generalisation to fewer degrees of symmetry for future work. We choose the initial configuration of ξ to be space dependent as where ∆ξ, σ and R 0 are the amplitude, width and the initial size of the perturbation respectively. We comment further on this perturbation shape in appendix A 2. The mass of the initial perturbation scales roughly as R 2 0 . We emphasise that this perturbation is non-linear, despite its moniker. Nevertheless, its massless nature means that it will propagate very close to the speed of light.
The initial staticity of the perturbation reflects the limited dynamics of a superhorizon perturbation that is frozen out. However, frozen out does not imply motionless. Nevertheless, the ξ field is massless and accelerates rapidly to velocities around c, regardless of its initial velocity. Thus, combined with the fact that a static configuration simplifies the constraints the initial data must satisfy (which is discussed more detail in appendix A 2), we deem this an acceptable simplification.
Given the initial static configuration, we expect to see the perturbation split into an infalling mode, which drives the PBH formation, and an outgoing mode, which rapidly disperses. Our simulation box has periodic boundary conditions -we ensure that the simulation domain is sufficiently large that the outgoing mode does not reach the boundary before PBH formation takes place. The background scalar field φ starts from rest, so thaṫ φ = 0 and the initial Hubble parameter in the absence of inhomogeneities is . Since the configuration of ξ breaks the homogeneity of the initial spatial hyperslice, to set up the correct initial conditions for the metric, we will solve the Hamiltonian constraint. We choose a conformally flat ansatz for the 3-metric γ ij , Then, the Hamiltonian constraint reduces to an equation for the conformal factor ψ where Here the local expansion K is the trace of the extrinsic curvature, K = TrK ij . Eqn. (11) then becomes We choose an initially expanding spacetime with K = −3H 0 , so that the periodic integrability condition is satisfied 4 [122][123][124]. The eventual Hamiltonian constraint only depends on the radial coordinate r due to the spherical symmetry of the setup, and we solve for the conformal factor ψ numerically

II. PRIMORDIAL BLACK HOLE FORMATION
Our main scale of reference will be the initial size of the unperturbed Hubble horizon H 0 , which is fixed for all simulations by choosing the initial value of the scalar field φ to be φ 0 = 7.8×10 −3 M Pl , with m ≈ 10 2 H 0 . In the following, we will vary the initial size of the perturbation from subhorizon to superhorizon, R 0 H 0 ∈ [0.575, 1.6]. We will also vary the perturbation amplitude within the range ∆ξM −1 Pl ∈ [0.075, 0.12], whilst keeping the initial width fixed to σ 0 = 0.15H −1 0 , such that the ratio between the maximum gradient energy density to dark matter energy density is ρ ξ /ρ DM ∼ 1. As the perturbation mass is small compared to the Hubble mass in all scenarios we consider, the background expansion is not significantly influenced by the perturbation's presence.
Our simulations reach just into the superhorizon regime and at the moment, we are limited from exploring this regime further numerically by the computational costs of such simulations. Nonetheless, we will argue in section II B that our simulations already capture some superhorizon dynamics. Moreover, we argue that it may not be necessary to probe the superhorizon regime much further, as our model is effectively pressureless. The corresponding speed of sound is therefore small, and physics already propagate very slowly on scales just past the horizon size.
We find that, for both subhorizon and superhorizon perturbations, PBH formation occurs via two possible mechanisms -a direct collapse mechanism whereby the PBH is formed by the initial perturbation of ξ itself, and a post-collapse accretion mechanism whereby the initial perturbation of ξ sources a gravitational potential that then accretes the background dark matter φ until a PBH 4 The initial energy density of the system is completely dominated by the scalar potential of the homogeneous dark matter field φ, which allows us to neglect the perturbation field ξ. Then, the main contribution to the initial energy density is given by the (homogeneous) value of the potential and thus K 2 = 24πV (φ 0 ). forms. What determines the type of PBH formation process depends (unsurprisingly) on both the geometry and mass of the initial perturbation shell, as well as the expansion rate of the background cosmology. We will discuss these two mechanisms below 5 .

A. Direct collapse
In the direct collapse scenario, the perturbation collapses towards its geometric centre (to which we will henceforth simply refer as the centre) and forms a black hole directly on its own, without significant accretion of the background DM density. We will now estimate the time a shell takes to undergo direct collapse. The perturbation field ξ is massless, and hence if we ignore the backreaction of the shell on the background geometry, it propagates along null-like geodesics 6 . In an FLRW background, the scale factor a is given by the null element dt 2 = a 2 (t)dr 2 . Solving this kinematic equation, the co-moving radius of the shell is then where we set a 0 ≡ 1 to be the initial scale factor at the initial time. The value of the scale factor at the moment the shell collapses to the centre a * is then the solution to the equation r shell (a * ) = 0, namely which is roughly a Hubble time. Notice that a * is independent of the initial mass and depends only on R 0 . We show in Fig. 2 that this analytic estimate is in good agreement with our numerical results.
To determine whether or not a given initial perturbation shell will undergo direct collapse into a black hole, we consider the width of the shell at the time when the shell reaches the centre σ * = σ(a * ). Ignoring backreaction again, since the field ξ is massless, the width of the shell as it collapses towards the centre scales as the expansion rate, i.e. σ(a) = σ 0 a .
Thus the width of the shell when it reaches the centre is simply σ(a * ) = σ 0 a * . At this moment, applying the hoop conjecture 7 suggests that if the condition is satisfied, where M infall is half 8 the initial mass of the shell obtained by integrating the gradient energy of the 6 If the perturbation field ξ instead has mass m ξ ∼ m φ , then the collapse time is roughly the free fall time τ ff ∼ M Pl /m ξ ξ which is also roughly one Hubble time. 7 The presence of an expanding background modifies the hoop conjecture somewhat in general [125], but we checked that the effects are negligible in our analytic estimates. 8 The infalling mass is half the initial mass, since the other half will radiate outwards to infinity, so the shell's initial (vanishing) momentum is conserved.

FIG. 3.
Evolution of the local expansion K and energy densities ρDM and ρ ξ at the centre of the collapse r = 0, as a function of the Hubble parameter H(t) at infinity -recall that K > 0 corresponds to locally collapsing spacetime. A representative subhorizon (superhorizon) is shown in thin light blue (thick dark blue) in the accretion collapse case. The top, middle and bottom panels show the evolution of the expansion, the background energy density and gradient energy density respectively. Initially, the background energy density decays as ρDM ∼ a(t) −3 . When the perturbation reaches the centre (dotted vertical lines) and disperses, gravitational effects decouple the system and stop the local expansion, acting as a seed for the accretion of the background matter ρDM. The accretion of the background matter continues until and after a black hole forms (dashed vertical lines).
profile ξ(r) roughly given by Eqn. (9) in flat space then a black hole will form. This result is again in good agreement with our numerical results, as shown in Fig.  2.
The fact that such simple estimates agree with our numerical results suggests that the backreaction of the perturbation on the background dynamics is not very important, at least at the level of determining when and how a black hole will form, even if the shell density is large and locally ρ ξ > ρ DM . This is backed up by our numerical simulations, where we see that the ρ DM profile is not strongly affected by the presence of ρ ξ , at least initially, as can be seen in a video of the numerical evolution of the energy densities here [118].

B. Accretion collapse
On the other hand, if σ(a * ) > 2GM infall , a black hole does not form directly. In this case, the energy density of the perturbation ρ ξ disperses after reaching the centre and becomes locally sub-dominant to the background energy density ρ DM . Nevertheless, the presence of ξ generates a gravitational potential well in the centre, which seeds accretion of the background DM and eventually causes a collapse into a black hole. We illustrate this process in Fig. 3.
In this phase, the initially homogeneous and expanding background spacetime is made to locally collapse by the perturbation, with the expansion K locally changing sign from negative (expanding) to positive (contracting), decoupling the region near the centre from the rest of the expanding background. The local dark matter begins to accrete at an extremely high rate δρ DM /ρ DM ∝ a 34 , as shown in Fig. 4. Once sufficient DM mass has accumulated, a PBH forms. This process takes around an e-fold to complete. This rapid accretion rate is much higher than that predicted from linear theory, which is δρ DM /ρ DM ∝ a, indicating that the process is highly nonlinear.
From our simulations, we note two salient points. Firstly, if we consider shells that undergo accretion collapse, for fixed initial amplitude ∆ξ, the smaller the initial R 0 (and therefore the smaller the mass) of the initial perturbation, the more massive the initial mass of the PBH. This somewhat counter-intuitive result is due to FIG. 5. Evolution of the dark matter energy density at the centre of the collapse for a set of initial radii R0 and amplitudes ∆ξ. For same radii perturbations, accretion begins at the same time. However, the accretion rate is larger for larger amplitudes, which results in the formation of a black hole at an earlier time.
the fact that the PBH forms via accreting DM, thus a more massive seed will generate a steeper potential well, and hence the Schwarzschild radius is reached earlier and at a smaller value for the PBH mass. To confirm this, we checked that keeping R 0 fixed but increasing ∆ξ also yields a less massive initial PBH -this is true for both subhorizon and superhorizon cases, as can be seen in Fig.  2. In Fig. 5 we plot the dark matter energy density for two different values for R 0 and ∆ξ. We confirm that larger amplitudes (and thus more massive seeds) result in a faster accretion rate.
Secondly, as R 0 approaches H −1 0 , the expansion rate of the universe begin to exert a competing effect. For shells with larger R 0 , it takes longer for the shell to reach the centre, and thus a smaller ρ ξ and less steep potential when accretion begins. This leads to an increase in the initial mass of the PBH following our argument aboveresulting in the "bump" in the initial mass of the black holes (e.g. the black dots in Fig. 6).

III. PBH GROWTH AND FINAL MASS
In the cases of both direct and accretion collapse, the initial mass of the PBH formed is small compared to the Hubble horizon, M BH H 0 ∼ 10 −2 M 2 Pl -see Fig. 6. Once the initial PBH has formed, the PBH accretes DM from its surroundings in the growth phase at a rate that depends on the steepness of the potential and the density of the surrounding DM "scalar cloud" [126][127][128][129]. In general and regardless of the details of the parameters, we find that the initial accretion rate is much higher than the linear theory prediction of δρ DM /ρ DM ∝ a, as mentioned above. This growth rate is roughly constant, at least initially, and its contribution to the mass of the PBH will rapidly dwarf that of its initial mass. In a matter dominated universe, naive Newtonian collapse suggests that the maximum mass of the black hole is bounded by M BH H ∼ αM 2 Pl [1], where α 1 is some constant which depends on the exact details of the accretion. This suggests a self-similar growth at some equilibrium point. In references [130,131], it was demonstrated numerically that while the initial growth can be rapid, it will not achieve self-similar growth as accretion is not efficient once the black hole decouples from the background spacetime. However, these works used a stiff massless scalar field as ambient matter instead of a massive scalar field, which more accurately models the ambient DM. From Fig. 6, we find that M ∼ H −β where β 1. As M approaches the Hubble horizon, we expect β ≤ 1 although unfortunately, we were unable to track the growth of PBH beyond a few factors of their initial mass, as the numerical cost becomes prohibitive.
As long as the universe is dominated by DM, the black hole will continually accrete and grow without end. This would be the case if the PBH is formed in the present late time DM dominated epoch -however such late time PBH has already been ruled out [27,28]. As we mentioned in the introduction, we consider instead an early phase of DM domination before transitioning into the era of radiation domination prior to the onset of Big Bang Nucleosynthesis (BBN), i.e. before the temperature of the universe is around 1 MeV. This provides a natural cut-off for the growth of the PBH.
Nevertheless, if we assume that the rapid growth we observe continues until M BH H ∼ M 2 Pl , and that the BH grows self-similarly after, it is implied that the final mass of the PBH is independent of when it forms. This means the final mass of the PBH is given by where T is the temperature of the universe at the onset of radiation domination. Taking T BBN = 1 MeV as the natural cut-off for the growth of the PBH, the most massive black holes that can be formed via this accretion mechanism are M BH ≈ 10 38 g ≈ 10 5 M [132,133].
On the other hand, if the PBH growth asymptotes to a slower rate than the self-similar rate, or achieve selfsimilarity before M BH ∼ H −1 M 2 Pl , then our simulations suggests that M BH 10 −2 H −1 M 2 Pl , where H is the Hubble parameter when the PBH forms. This means that PBH formed around T ∼ 5 MeV, M BH 40M could form the basis of the population of massive BH that are being detected today by the LVK observatories.

IV. SUMMARY AND DISCUSSION
In this paper, we demonstrate that superhorizon nonlinear perturbations can collapse and form PBHs in a Pl = 0.075. The growth rate of the black hole mass is larger for larger shells, because they source a larger gravitational potential. Black dots correspond to the initial black hole masses at formation, identified using an apparent horizon finder. . matter dominated universe, using full numerical relativity. We show that, depending on the mass of the initial perturbation shell, this happens via either the direct collapse or the accretion collapse mechanisms. We provide an analytic criterion Eqn. (18) using the hoop conjecture to determine which mechanism is relevant for a given setting, and compute the timescale of collapse using the same prescription. Despite the O(1) non-linearity, we find that the dynamics of collapse can be modeled as a simple superhorizon mass shell collapsing in an expanding background. This suggests that semi-analytic estimates of PBH formation in a matter dominated era are broadly accurate.
On the other hand, details matter. We showed that even in the cases where the perturbation is insufficient on its own to form a PBH in a direct collapse, non-linear accretion rates are far higher than what standard linear theory predicts, causing a rapid collapse into a PBH via accretion of ambient DM. In both the direct collapse and accretion collapse formation cases, the initial mass of the PBH is roughly M BH ∼ 10 −2 H −1 M 2 Pl , but formation is followed by an extremely rapid growth M ∝ H −β where β 1. Presumably, this growth will asymptote to either the self-similar rate β = 1 or the decoupled rate β < 1 [130,131]. Interestingly, even if the self-similar rate is not achieved, the fact that most of the mass of the PBH is gained through post-formation accretion suggests that there might be a mechanism to generate PBHs with non-trivial spin. Such non-trivial spin might for exam-ple be generated by the collapse of a non-spherically symmetric shell, even if the shell is initially spinless. In that case, the PBH might not form in the centre of the initial mass distribution and thus form with spin, whilst outgoing radiation carries away angular momentum of opposite sign, such that angular momentum is still globally conserved, as suggested by [112,114]. We will explore this possibilty in a future publication.  7. Initial profiles of the massless ξ field, as a function of radial distance away from the center. Depicted are initially subhorizon, horizon and superhorizon sized perturbations, with same amplitude and thus same energy densities.

ACKNOWLEDGMENTS
Appendix A: Numerical methodology

Evolution equations
This work was written based on simulations run using GRChombo [134], with the CCZ4 formulation of the Einstein equations [135]. This formulation relies on the 4dimensional spacetime being foliated into 3-dimensional non-intersecting hyperslices, whose intrinsic curvature is described by a spatial metric γ ij , whilst their embedding in the 4-dimensional spacetime is encoded in the extrinsic curvature K ij . The line element is decomposed as where the lapse α and the shift vector β i are userspecified gauge functions. The spatial metric γ ij is often written as a product of a conformal factor ψ and a background (or conformal ) metricγ ij , so that the determinant of the conformal metric equals one, Time evolution proceeds along the vector t a = αn a + β a , where n a is the unit normal vector to the hyperslice that is being evolved. The gauge functions α and β i are specified on the initial hyperslice, and then evolved using evolution equations suitable for long-term stable numerical simulations. The choice in this work is il γ lk,j +γ lj,k −γ jk,l .

Initial data
The matter content of this work is comprised by a massive φ field that dominates the background dynamics, and an inhomogeneous massless ξ field which provides the local overdensity through its gradients. We choose potentials and spherically symmetric initial field configurations Examples of initial field profiles for the massless ξ field are given in Fig. 7. The hyperbolic tangents allow us to localize the gradients in the ξ field, and thereby its energy density, as the field is without potential. In our simulations, mH −1 0 = 62.6. We also choose a conformally flat metricγ ij = δ ij , so the energy density on the initial hyperslice is given by which corresponds to a shell-like overdensity in a dark matter environment. This setup is spherically symmetric and we reduce the computational cost of evolution by simulating one eighth of the system using symmetric boundary conditions. A schematic depiction of the initial setup is given in Fig. 8.
The equation of motion for a massive homogeneous field φ is given by the Klein-Gordon equation where H ≡ȧ/a is the Hubble function defined with the scale factor of the universe a(t).
For the universe to stay matter dominated for an extended period of time, we require the solution of Eqn. (A8) to be an undamped oscillation. This constrains the initial value of φ to and we use φ 0 = 7.8 × 10 −3 M Pl for all our simulations. Momentum constraints are trivially satisfied, and we solve the Hamiltonian constraint (14) numerically to find the conformal factor ψ.