Before Reheating
Perhaps all of the dark matter is made of black holes with the mass of an asteroid, 10¹⁹–10²⁰ grams. Stars do not make such black holes. They could only have been born at the very beginning, from slightly denser clumps of matter. We are testing whether this is realistic if, for a while, the early Universe was filled with cold “dust” rather than hot radiation.
The whole calculation comes down to a single number: how strong a clump has to be to collapse into a black hole before the dust era ends. The lower this threshold, the more black holes are born.
A Universe made of dust
After inflation, space may have been filled with a heavy scalar field. It oscillates around its minimum, and on average these oscillations behave like a pressureless gas of very heavy particles. Cosmologists call this dust, and the epoch the dust era, or the early matter era.
For black holes this is a generous environment. Radiation pressure, which scatters clumps in the usual hot Universe, is absent here. Any overdensity keeps growing until the expansion lets go of it, and then it starts to fall in on itself.
But dust has its own catch: the particles do not collide. They do not stick together at the centre; they fly straight through it and out the other side, like a pendulum. Whether a black hole is born depends on whether enough matter piles up at the centre for light to stop escaping from there.
Dust versus a particle swarm
For the fall of ideal dust there is an exact solution: Lemaître–Tolman–Bondi, or LTB. In it the shells of the clump, like the layers of an onion, fall towards the centre and never pass through one another. The innermost shell is the first to collapse to a point, and this happens at a strictly defined moment, t_C(0).
The real field cannot do that. We describe it as a swarm of collisionless particles. These are the Einstein–Vlasov equations, the same ones used to describe the stars in a galaxy, but in full general relativity. The particles pass through the centre, fly apart and come back, and a multi-stream core builds up around the centre. A horizon appears only when this core becomes compact enough.
What the code does
The spherical clump is sliced into tens of thousands of concentric shells, each of which is a stream of particles. The code simultaneously solves Einstein’s equations for the geometry and the equations of motion of the particles in that geometry. The calculation runs while the Universe expands a hundredfold or more. At every step the code looks for a trapped surface, a sphere from inside which even light can no longer get out. Its appearance is the birth of a black hole.
One of the tricks is the choice of time. The code measures time by the clocks of a distant observer in the expanding Universe. These clocks can be “laid” over space in such a way that they run through the future horizon and do not run into the singularity.
The cast
The main knob of the calculation. It is the peak curvature of space at the centre: ζ = μ·e−k²r²/6, the profile of Yoo and colleagues (2026). At horizon entry the density at the centre is about 2.4·μ above the mean, so for μ = 0.04 it is 10% higher. After that this excess grows with the expansion. Weak clumps are incomparably more common, so every step down in μ matters.
While a clump is larger than the horizon, its parts cannot exchange signals and do not “know” that together they are denser than their surroundings. t_H is the moment when the horizon grows to the size of the clump and gravity starts pulling it together as a whole. All times here are in units of t_H. In the dust era the size of the Universe grows as t2/3, so by 1000 t_H it is 100 times larger.
The mass inside the horizon at time t_H, essentially the whole clump. The black hole is born small, about 0.01 M_H, and then grows by swallowing infalling matter.
When the centre of the clump would collapse to a point if the matter were ideal dust. For μ = 0.05 this is 61 t_H, for μ = 0.02 it is 221 t_H, for μ = 0.01 it is 607 t_H. A handy yardstick: real horizons in a particle swarm appear at 1.1–1.3 t_C(0).
The field does not oscillate forever. It decays into hot radiation (reheating), pressure stops the collapse, and whatever has not made it in time is never born. We test D = 100, 300 and 1000 t_H.
The main result: the smallest amplitude at which a horizon manages to appear before the deadline D.
Random particle velocities, as a fraction of the speed of light. At σ = 0 everything falls straight towards the centre. Dispersion acts like thermal pressure: particles miss the centre, and the collapse slows down. The sources are the wave nature of the field (σ ≈ √6/q) and small-scale ripples that “unfreeze” during contraction.
How many times the oscillation frequency of the field exceeds the expansion rate. For asteroid-mass black holes q ≥ 3·10⁹: the field makes billions of oscillations per expansion time. So replacing the field with particles is accurate to about 10⁻⁵, and the wave dispersion is negligible, about 10⁻⁹.
ν is the peak height in units of the rms fluctuation: ν = 4 means a rare peak, four times higher than a typical one. Δ is the width of the perturbation spectrum: a narrow one (Δ = 0.1) gives an almost Gaussian clump, a broad one (Δ = 0.5) a sharp peak with gently sloping wings.
Twice the mass inside a sphere divided by its radius (in units with G = c = 1). For an ordinary star this is a few millionths. When 2m/R reaches 1, the sphere becomes a horizon.
The clocks of a distant observer (coordinate time t) and the clocks of the infalling shells themselves (proper time τ). Gravity slows time inside the clump, and the shells’ clocks lag by 5–8%. We give both values in the results.
Key assumptions
- A1
Sphere
The clump is perfectly round. Real clumps are flattened. A Newtonian 3D simulation shows that typical ellipticity can raise the threshold for the 100 t_H deadline from 0.039 to 0.06–0.08. Exactly where inside this bracket will be decided by a 3D calculation in GR (section “The 3D stage”). So the spherical threshold is a lower bound.
- A2
Cold start
The baseline case has zero velocity dispersion. Dispersion is studied in separate series of runs (in detail in the section “Velocity dispersion” below). If it is present from the very start, the threshold shifts by no more than 15%. If, however, the dispersion switches on right after turnaround, the picture changes strongly (blue squares in the chart below).
- A3
Instant reheating
We assume that collapse is impossible after the deadline D. In reality radiation slows the infall gradually, so the real boundary is blurred.
- A4
Convergence
Every result is re-checked on a finer grid and with more shells. This turned out to be critical: at small μ thousands of shells swirl in the core, and with too few of them the answer shifted by 8%. The answer stops changing at about 1000 shells in the core. And a coarse grid used to draw a false threshold “floor” near μ ≈ 0.014. Once the grid was rebuilt to follow the expansion, this floor disappeared.
Results so far
Each point below is a separate run lasting many hours. The crossing of the orange line with a deadline gives the threshold.
| μ | tC(0), t_H | horizon, t_H | in units of tC(0) | with early dispersion |
|---|---|---|---|---|
| 0.05 | 61 | 68 | 1.12 | 76 t_H |
| 0.04 | 83 | 96 | 1.15 | 125 t_H |
| 0.03 | 124 | 150 | 1.21 | 260 t_H |
| 0.025 | 161 | — | — | none by 450 t_H |
| 0.02 | 221 | 284 | 1.28 | — |
| 0.017 | 280 | 364 | 1.30 | — |
| 0.010 | 607 | 793 | 1.31 | — |
The final cold threshold map: μ_th(100 t_H) ≈ 0.039, μ_th(300 t_H) ≈ 0.019, μ_th(1000 t_H) ≈ 0.0086. The key point: a smooth, cold Gaussian clump has no floor, at least down to μ = 0.010. The longer the dust era lasts, the weaker the clumps that manage to collapse.
This is the result for a smooth, cold clump. Small-scale ripples that stir up the velocities right after turnaround break this picture. Already with moderate ripples (k/k₀ = 10) an effective floor appears near μ ≈ 0.027: at μ = 0.025 there is no horizon even by 450 t_H. With coarser ripples (k/k₀ = 3) there is no spherical black hole for any μ ≤ 0.05. Details in the section “Velocity dispersion” below.
Latest run
The forecast came true: the horizon appeared at the lower edge of the expected window, 207 t_H before the deadline. The delay factor has almost stopped growing: 1.28 → 1.30 → 1.31 for μ = 0.02 → 0.017 → 0.010. With it, the threshold for the 1000 t_H deadline comes out at about 0.0086. So with a long dust era, clumps 4–5 times weaker than for a 100 t_H deadline are enough.
Velocity dispersion: help or hindrance
The cold calculation assumes that all particles fall straight towards the centre. In reality they have small random velocities, and the most important one is sideways. A particle with a sideways velocity misses the centre and passes at some distance from it. This distance is called the pericentre. If the pericentres are larger than the future horizon, the core stays puffy and the horizon is born later.
But dispersion does not only work against collapse. A small kick given at horizon entry changes the binding energy of the inner shells by more than the binding energy itself. Their binding energy is tiny, so the kick easily outweighs it. As a result, some shells become more tightly bound, fall earlier and assemble the core faster. That is why weak dispersion, paradoxically, speeds up the birth of the black hole. On a converged grid (32 shells per cell) this speed-up amounts to 6–13% of t_C(0).
Where the dispersion comes from
Waves of the field itself
The field is a wave, and its “particles” carry momentum, like a quantum particle. At horizon entry this gives σ ≈ √6/q. For real asteroid-mass black holes q ≥ 3·10⁹, and the dispersion is negligible: about 10⁻⁹.
Tested for q from 25 to 2450: the threshold shifts by no more than 15%, and upwards only for q ≲ 30. This axis is closed. A 3D check at μ = 0.3 confirmed it: in the horizon growth above 0.2 M_H a dispersion of σ = 0.027 is not visible.
Small halos, after Harada
Small-scale ripples on top of the clump collapse early into tiny halos. Their internal velocities are released when the density of the clump catches up with the density of the halos. When this happens depends on the scale of the ripples.
Right at the caustic even 0.35c changes nothing. If, however, the release starts early, from 0.61 t_C(0), the horizon is late by 0.27 t_C(0).
Ripples during contraction, after Ebrahimian
Small-scale waves become nonlinear during the contraction itself and stir up the velocities right after each shell turns around. The coarser the ripples (the smaller k/k₀), the stronger the kick.
A kick of 0.03c cancels the horizon for every μ tested. 0.01c delays it by 0.1–1 t_C(0), 0.003c by no more than 0.15.
Dispersion from the start: converged values
| μ | cold | σ = 0.001 | σ = 0.01 | σ = 0.03 | σ = 0.1* |
|---|---|---|---|---|---|
| 0.04 | 1.15 | 1.06 | 1.08 | 1.14 | 1.20 |
| 0.03 | 1.21 | 1.05 | 1.07 | — | 1.15 |
Horizon time in units of t_C(0), 32 shells per cell in the core. *σ = 0.1 was computed with 8 shells per cell, where the cold value is 1.12–1.13, so it is only roughly comparable with the other columns. On the same grid at μ = 0.05–0.07 it gives 1.24–1.33. For σ ≥ 0.03 the newborn black hole is heavier: 0.02–0.04 M_H instead of 0.01.
Small-scale ripples: full scan
| source and kick | μ = 0.05 | μ = 0.04 | μ = 0.03 |
|---|---|---|---|
| cold, no dispersion | 1.12 | 1.12 | 1.13 |
| Ebrahimian, k/k₀ = 3 (≈ 0.03c) | none by 150 t_H | none by 300 t_H | none by 400 t_H |
| Ebrahimian, k/k₀ = 10 (≈ 0.01c) | 1.24 | 1.51 | 2.09 |
| Ebrahimian, k/k₀ = 30 (≈ 0.003c) | — | 1.21 | 1.35 |
| k/k₀ = 10, released at 0.6 R_max | — | 1.42 | 1.79 |
| k/k₀ = 10, released at 0.3 R_max | — | 1.33 | — |
| Harada, k/k̃ = 5 (0.17c from 0.61 t_C(0)) | — | 1.39 | — |
| Harada, k/k̃ = 10 (0.17–0.35c at the caustic) | — | 1.09–1.13 | 1.12 |
Horizon time in units of t_C(0), 4 shells per cell. k/k₀ is how many times smaller the ripples are than the clump. Their amplitude is ζ = 0.05 throughout.
Thresholds including ripples released right after turnaround:
| ripples | μ_th(100 t_H) | μ_th(300 t_H) |
|---|---|---|
| none (cold clump) | 0.039 | 0.019 |
| k/k₀ = 30 | ≈ 0.040 | ≈ 0.025 |
| k/k₀ = 10 | ≈ 0.045 | ≈ 0.028 |
| k/k₀ = 3 | no spherical black hole for any μ ≤ 0.05 within 300 t_H | |
Takeaway. Dispersion from the waves of the field itself is almost harmless: it switches on early and has time to cool down with the expansion. The danger is small-scale ripples, if they stir up the velocities right after turnaround, while the particles are still slow. The threshold “floor” is set by the strength of the kick together with the moment it switches on. The same kick released later does less harm, and right at the caustic it barely matters. Nonlinear ripples 3–10 times smaller than the clump, with ζ ≈ 0.05, are exactly what a realistic spectrum produces. So the cold conclusion “there is no floor” is fragile. This agrees with the Newtonian estimate of Ebrahimian and colleagues and with the 3D bound of Yoo and colleagues. The decisive test is a 3D calculation with a real small-scale spectrum. In it the ripples are not modelled as a kick but computed directly.
Clump shape
Everything above referred to a single profile, a Gaussian one, as used by Yoo and colleagues. But the shape of a clump depends on what the spectrum of primordial perturbations was like. A narrow spectrum gives an almost Gaussian peak with a shallow density dip around it. A broad one gives a sharp peak with long, gently sloping wings. The limiting case is a flat core with a sharp edge. The narrow- and broad-spectrum profiles are the mean shape of a rare peak of height ν = 4 (four times the rms fluctuation). The flat core is taken as a limiting control case.
What decides the outcome is how closely in step the shells fall. If they arrive together, the core gains mass at once. If their arrival is spread out in time, the centre is fed slowly, and the particles have time to fly apart.
| profile | μ = 0.05 | μ = 0.03 | μ = 0.02 | black hole mass |
|---|---|---|---|---|
| Gaussian (Yoo), for comparison | 1.12 | 1.21 | 1.28 | ≈ 0.001–0.01 M_H |
| narrow spectrum (Δ = 0.1) | 1.09 | 1.10 | 1.16 | ≈ 0.001 M_H |
| broad spectrum (Δ = 0.5) | 1.48 | 1.80 | none by 1.8 | ≈ 0.001 M_H |
| flat core | 1.05* | 0.99* | 0.97* | 0.77–0.82 M_H at once |
First-horizon time in units of t_C(0) (*for the flat core, t_C(r_m)), cold case, 32 shells per cell. Each profile has its own t_H: for example, at μ = 0.03, t_C(0) = 90 t_H for the narrow spectrum and 85 t_H for the broad one. A dispersion of σ = 0.03 barely changes the narrow spectrum (1.11) and the flat core (0.98), but delays the broad one by another 12% (1.66 at μ = 0.05).
A small seed
The horizon is born at 1.1–1.3 t_C(0) around a small core of a few thousandths of M_H. Then the black hole grows by swallowing infalling matter.
Late or nothing
The centre is fed too slowly. The horizon is 50–80% late, and at μ = 0.02 there is none even by 1.8 t_C(0): the compactness is stuck at 0.08.
A big black hole at once
The whole edge falls at once and immediately collapses into a black hole with almost the entire horizon mass. This is the classic collapse of a uniform ball (top-hat).
What if the clump is flattened?
Real peaks are not round. A flattened clump collapses along its axes one at a time: first into a pancake, then into a filament, then into a lump. We first estimated this with the Bond–Myers homogeneous ellipsoid model. According to it, the last axis collapses only 5–9% later than a sphere. Then we checked with a Newtonian 3D simulation: 7 million particles, a 192³ grid, one peak with the ellipticity typical of ν = 4 peaks. It showed something quite different.
A peak whose density falls off towards the edges collapses much more anisotropically than a homogeneous ellipsoid. At e = 0.13 the pancake forms at 0.77 of the sphere’s time (the ellipsoid gives 0.91), the filament at 0.94–0.97, and the last axis collapses only at 2.7 times the sphere’s time (the ellipsoid gives 1.05). At e = 0.19 the figures are 0.62, 0.89 and 1.7 respectively. The central lump still assembles, only later.
The Newtonian simulation cannot see the horizon itself: its resolution only reaches a compactness of 2M/R ≈ 0.15. The horizon will be born somewhere between the filament stage (if the filament is compact enough) and the collapse of the last axis. Hence the bracket for the 100 t_H deadline: the threshold lies somewhere from 0.039 to 0.06–0.08. The 3D bound of Yoo and colleagues, 0.045–0.05, lies inside it.
Small-scale ripples in the same 3D simulation (ζ = 0.05 on scales 2–10 times smaller than the peak) broke the region into fragments in one of three realisations, and no compact centre formed. In the other two the centre assembled. This is the 3D analogue of the result with a kick after turnaround, only now from a real spectrum. Three realisations are not yet statistics.
Takeaway. Clump shape is a first-order effect, comparable to small-scale ripples, and any threshold must be quoted together with the shape of the spectrum. Ellipticity may also turn out to be important: the homogeneous ellipsoid model badly underestimates it. Where the true threshold lies inside the 0.039…0.08 bracket can only be decided by a 3D calculation in full GR. That is the next stage.
The 3D stage
The spherical part of the programme is almost complete: velocity dispersion, small-scale ripples and profile shape have been computed. Two questions remain open, and both are three-dimensional: how much ellipticity and real small-scale ripples raise the threshold. The Newtonian 3D simulation gave a bracket but cannot say exactly when the horizon is born. For that, full general relativity in three dimensions is needed.
Why this is hard
The first horizon is tiny. By the time it appears, the Universe has grown about 1000-fold, and its radius is about 0.3/k, a ten-thousandth of the computational domain. A uniform grid, such as the 80³ grid of Yoo and colleagues, will see the horizon only once it has grown to 0.05–0.1 M_H, that is, noticeably later than its actual birth. What is needed is an adaptive grid that refines itself where matter gathers. A horizon finder and a cosmological background are also needed.
Why GRChombo
| code | GR | matter | adaptive grid | role |
|---|---|---|---|---|
| GRChombo | full | scalar fields, no particles | yes, mature | main code |
| CosmoGRaPH | full | fluid, experimental particles | in papers, not in the release | cross-check |
| GRAMSES | approximate near horizons | particles | yes | request the code |
| COSMOS (Yoo) | full | particles only in a private version | none | via a request to the authors (V0) |
GRChombo is an open-source code (BSD-3 licence) with full GR, an adaptive grid and a horizon finder. It has already been used for primordial black holes in a matter era: de Jong, Aurrekoetxea and Lim (2022–2023) modelled the background with a massive scalar field. It had no particles, and neither does its AMReX port, GRTeclyn. We wrote the particle module ourselves (stage 2 below).
Stage 1: scalar field — a bounce instead of a black hole done
GRChombo is built and passes all 17 of its tests, including the horizon finder. The calculation starts from the exact dust solution at the moment the centre turns around (0.5 t_C(0)). The matter is a scalar field with the same density, so the constraint equations are satisfied exactly. Two obstacles had to be worked around:
- !
Standard clocks freeze
For a slow collapse from a dilute state, the standard time gauge (1+log) stops the clocks in the core at about 0.17 t_C(0), and the collapse never happens on the time slices. We introduced a soft gauge tied to the profile: far away it coincides with cosmic time, and in the core it slows down as moderately as our 1D clocks do.
- !
A bug in the GRChombo example
In its cosmology example, the mean curvature is assigned to an object that never reaches the evolution equations.
Why the field bounces and the particles do not. A field with q = 45 is not a particle swarm but a wave. A self-gravitating wave has a maximum mass that its “quantum” pressure can still hold back from collapse: the Kaup mass, 0.633/m. Here it is 0.028 M_H. But the first horizon in the 1D map is born on a core of only 0.001–0.016 M_H. The wave pressure bounces such a core back before enough mass reaches it. For the field to behave like dust, the Kaup mass must be smaller than the mass of the first horizon. For a cold core at μ = 0.1 that means m·t_H ≈ 10⁴, which is prohibitively expensive.
For real asteroid-mass black holes q ≥ 3·10⁹, and the Kaup mass is negligible. The bounce is a property of the cheap numerical stand-in, not of the physics. The lesson for the programme: a scalar field will not do as a substitute for particles in 3D. The planned field run at μ = 0.3 was cancelled, and the work moved straight to particles.
Stage 2: particles in GRTeclyn — the 3D horizon matched the sphere done
The particle module for GRTeclyn had to be written from scratch. The particles move along geodesics in the 3D metric and feed their energy, momentum and pressure back into it. It was tested step by step:
A homogeneous dust Universe
Expansion, curvature and density match the exact solution to 0.2–0.3%; mass and particle number are conserved exactly.
A round clump up to the caustic
At μ = 0.3, density, curvature and geometry match the exact LTB solution to 0.2–0.5% wherever the profile is resolved, right up to 0.96 t_C(0).
Particles on nested grids
At the boundary between coarse and fine grids the particles at first produced a density skew of ±12%. After the particle deposition was reworked it is below 0.1%, and the central density peak is resolved at 587 times the mean instead of 126 on a single grid.
Then came a run through horizon formation: μ = 0.3, round clump, 17.8 million particles, 7 grid levels. The finest cell is 64 times smaller than the coarse one and 512 times smaller than r_m.
First attempt: a horizon, but 4% late. The horizon was found and grew to 0.35 M_H. To compare it with the 1D map, a common clock is needed. The clocks of the distant Universe in 3D and the clocks of our spherical slices in 1D run differently, and the difference reaches 0.35–0.5 t_C(0). So each particle was given its own clock, and the comparison used the time lived by the matter at the horizon itself. On these clocks the 3D horizon lagged the 1D one by 0.04–0.055 t_C(0) at any mass, that is, by about 4%. A run with a grid twice as coarse gave the same shift. So resolution was not the cause.
The cause: a bug in the particle equation of motion. The force produced by the uneven flow of time (the gradient of the lapse α) had been multiplied by α one extra time. All earlier tests ran at α = 1, where the bug is invisible. In the working gauge α is 0.3–0.8 in the core, so this force was 20–70% weaker than it should be, and the collapse was held back. The clue came from the central clock: the lag began exactly when α at the centre dropped below 0.85.
After the fix: agreement. Up to the caustic the shells in 3D reproduce the exact LTB solution to better than 0.5% in radius and mass. The first horizon is seen at 1.23 t_C(0) on the distant clocks, with a mass of 0.065 M_H.
A gauge in which the interior does not blow up. In the soft gauge that suits the field, the calculation broke down at 0.27 M_H: the slice kept pushing into a “puncture” two cells wide. With the standard gauge (coefficient 2 instead of 0.3) the clocks in the core freeze by themselves, and the horizon can be followed to the end, 2.2 t_C(0), up to a mass of 0.58–0.70 M_H. The price is that the horizon is seen later: first at 0.18 M_H rather than 0.065. For flattened peaks the chosen scheme is to start in the soft gauge and switch to the standard one after the horizon appears.
This was the mandatory test for the 3D code: the round clump in 3D reproduced the spherical map, and the test is passed for horizon growth from 0.09 to 0.6 M_H. The moment of horizon birth itself cannot be compared yet. In 1D the first horizon is a core of 0.0025 M_H, while the smallest horizon the 3D grid can resolve is 0.06–0.18 M_H, depending on the gauge. The built-in horizon finder skipped over the small horizon. So the horizon is located with rays from the centre and with spheres along 96 directions: the horizon area agrees with the ray estimate to 2·10⁻⁵.
Warm particles in 3D. The same clump with a velocity dispersion of σ = 0.027 at horizon entry (this is q ≈ 90). In the horizon growth above 0.18 M_H the dispersion is barely visible: the mass is 0.4–1% lower, the time 0.1% later, and noise flattens the horizon by 0.2–0.5%. The 1D calculation says the same: such a dispersion changes only the very first, tiny horizon. That one cannot be made out in 3D, and here 1D remains the main tool.
Stage 1b: starting from horizon entry done
All the runs above started from the ready-made dust solution at turnaround, 0.5 t_C(0). But a flattened clump cannot be set up that way: there is no exact solution for it. So a start from the very beginning is needed, from the initial data of Yoo and colleagues, when the clump is still larger than the horizon (t_i, about 1800 times earlier than t_H).
This start turned out to be very sensitive. While the clump is larger than the horizon, the binding energy of each shell is a small difference between two large numbers. A relative error η in the mass that a shell “sees” changes its binding energy 100–1000 times more strongly. The usual deposition of particles onto the grid was off by 8·10⁻⁴, and the core came out 20% less bound. We had to:
Tune the particle masses
The masses are adjusted iteratively so that the density deposited on the grid matches the one required by the constraint equations. The error drops from 8·10⁻⁴ to 10⁻⁵.
Run the early phase on the simplest clocks
The shift gauge drags the particle lattice across the grid, which is unacceptable in the super-horizon phase. Up to 0.5 t_C(0) the calculation runs on geodesic clocks, and then the shells reproduce LTB to 10⁻⁴. The working gauge is switched on afterwards.
Fix two one-step lags
A test on a homogeneous Universe showed that the sources and the cosmological curvature lagged by one step. With a step of 3% of the age of the Universe, over hundreds of steps this made it “closed”. After the fix the homogeneous Universe holds to 10⁻⁶.
Result: a full run of the round clump from the t_i start (7 levels, 20 hours on a CCX33) reproduced the horizon map. It agreed with the late start to 0.007 t_C(0) and with exact dust to 0.010, and the horizon grew to 0.70 M_H. A side observation: Yoo and colleagues use one particle per cell and simple deposition on an 80³ grid, and such a scheme cannot hold the binding energy to the required accuracy of 10⁻⁴. This may be part of the reason their thresholds are higher than ours. Not yet quantified.
The first flattened clump in progress
Since 7 October we have been running e = 0.2 at μ = 0.3, with the same pipeline that passed all the checks on the round clump. It is at this ellipticity that Yoo and colleagues saw a horizon from μ = 0.05 upwards. We measure the time of the first horizon, its mass from its area and the growth of M_AH against proper time, all compared with the round clump.
What comes next
e = 0.2 at μ = 0.3
How much ellipticity delays and shrinks the horizon in full GR. In progress.
Flattened peaks at small μ
e = 0.13 and 0.19 at μ = 0.04–0.06, to narrow the threshold bracket 0.039…0.08.
A finer grid at the centre
The horizon is born across only 2–5 of the finest cells. To see it earlier and smaller, 2–3 more grid levels are needed, at roughly twice the cost.
What we plan to get
The first full-GR threshold for a typical flattened peak. Is it closer to the spherical 0.039 or to the Newtonian 0.06–0.08?
Whether the filament is compact enough for a horizon to form, or the final lump is needed. A Newtonian calculation cannot settle this.
Achieved: the round clump in 3D reproduced the spherical map and exact dust to about 1% in the proper time of matter, from 0.09 to 0.7 M_H, in two gauges and with two start options. Along the way, three bugs were found and fixed.
There is already a partial answer: at q = 45 the field bounces where particles produce a black hole. Wave pressure matters as long as the Kaup mass exceeds the mass of the first horizon.
In parallel, two requests remain open: to the Yoo group (V0), which has its own 3D code with particles, and for the GRAMSES code.
At the moment all 3D work runs on a rented CCX33. For the series of flattened peaks the plan recommends a more powerful AX102-1: 16 cores, 128 GB, €257.30 per month. The current two-core server is only good for 1D calculations and verification tests.
Checking the tools
Before trusting the new code with horizons, we separately verified the module that finds them (pilot PBH-VERIFY-01). The module computes the mass and the light-ray expansions and classifies spheres. One agent derived the physics, another wrote the code, and a third checked it. The module passed 34 of 35 independent tests and caught all 13 deliberately planted bugs. The most useful finding concerns work already done. The 1D criterion “2m/R ≥ 1 at one node or more” fires on a core one cell wide, where the direction of light rays is not yet resolved. The first reliable horizon appears about 1% later in time: at μ = 0.10 it is 28.43 t_H instead of 28.14. This affects the threshold maps at the percent level (not yet recomputed). The rule for the 3D stage that follows: a horizon counts only if the grid resolves it.
Why all this matters
How many black holes are born depends exponentially on the threshold. A factor-of-two shift in μ_th is the difference between “a negligible number of black holes” and “they make up all of the dark matter”. We are answering a question that nobody has yet solved in full general relativity: how much time does the early Universe need to turn barely noticeable density ripples into black holes?
The 3D stage is under way: for a round clump, the black hole from a particle swarm in full GR matched the spherical map to about 1%, and the first flattened clump is already running (section “The 3D stage”). Next, this calculation will decide where inside the 0.039…0.08 bracket the threshold for a typical flattened peak lies, and whether the threshold has a floor with real small-scale ripples.