Verification report: arXiv:2609.14218v1 (Yoo, Escrivà, Harada, Kohri)

Source used: the arXiv abstract page and the full PDF (https://arxiv.org/pdf/2609.14218v1, 29 pp., 13 figs, submitted 13 Sep 2026). The PDF was saved locally by WebFetch and I extracted its complete text with pdftotext and rendered the figure pages to images, so every quote below is from the paper itself, not from memory. Where I read a number off a plot or derived it from the paper's equations, I say so explicitly.

Ledger #

1. Bibliographic data and abstract (verbatim) #

Title: Simulation of PBH formation in a matter-dominated universe Preprint numbers: NU-QG-28, RUP-26-21. arXiv:2609.14218v1 [gr-qc] 13 Sep 2026. Comments: 29 pages, 13 figures. Categories: gr-qc, astro-ph.CO, hep-ph, hep-th.

Authors and affiliations: Chul-Moon Yoo¹·², Albert Escriv๷³·⁴, Tomohiro Harada⁵, Kazunori Kohri⁶·⁷·⁸·⁹ 1 Graduate School of Science, Nagoya University, Nagoya 464-8602, Japan 2 Kobayashi-Maskawa Institute for the Origin of Particles and the Universe (KMI), Nagoya 464-8602, Japan 3 Asia Pacific Center for Theoretical Physics, Pohang 37673, Republic of Korea 4 Department of Physics, Pohang University of Science and Technology, Pohang 37673, Republic of Korea 5 Department of Physics, Rikkyo University, Toshima, Tokyo 171-8501, Japan 6 Division of Science, NAOJ, and SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 7 Department of Astronomy, The University of Tokyo, Bunkyo-ku, Hongo, Tokyo 113-0033, Japan 8 Theory Center, IPNS, KEK, 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan 9 Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

Abstract:

"We investigate primordial black hole (PBH) formation during an early matter-dominated era using fully nonlinear numerical relativity. The initial condition is set by a functional form of the curvature perturbation including ellipticity, which makes the configuration triaxial. Two kinds of matter descriptions are considered: dust fluid and collisionless particles. In the dust fluid description, numerical computation crashes associated with the appearance of a singularity at which the fluid density diverges, unless the singularity is hidden well inside the apparent horizon. We found that, for the dust fluid description, to observe the horizon formation before calculations crash, the initial amplitude must be larger than the previous analytic estimation by a factor of 2. On the other hand, with the particle system description, calculations do not crash, and we may observe black hole formation after subsequent evolution of the system. Then the threshold of black hole formation is significantly smaller than the previous analytic estimation by an order of magnitude."

Section structure: I. Introduction; II. Basic setups in the form of 3+1 decomposition (A. BSSN, B. Gauge conditions, C. Scale-up non-Cartesian coordinates); III. Initial condition; IV. Spherically symmetric cases (A. LTB spacetime, B. Correspondence with long-wavelength solutions, C. Nakedness of the singularity and PBH formation, D. Numerical simulation of spherical cases with dust fluid); V. Non-spherical Collapse (A. Criterion from the Zel'dovich approximation and the hoop conjecture, B. Comparison with numerical simulation); VI. A collisionless particle system (A. …stress-energy tensor, B. Particle settings, C. Results); VII. Summary; App. A Linear perturbation equations; App. B Consistency between the BSSN simulation of collisionless particles and the framework of the exact LTB solution.

2. Curvature profile, metric ansatz, parameters #

Metric decomposition (Sec. II-A):

Eq. (2.1): ds² = −α²dt² + γ_ij(dx^i + β^i dt)(dx^j + β^j dt) Eq. (2.2): "γ_ij = a²ψ⁴γ̃_ij with det γ̃ = det f, where det f is the determinant of the 3-dim flat metric f_ij." Eq. (2.5): K_ij = a²ψ⁴Ã_ij + (1/3)Kγ_ij

Definition of Ψ (Sec. III):

"Specifically, using the function Ψ(x) := exp(ζ/2), we express the initial conditions for the geometrical variables as follows:" [Eqs. 3.4–3.9]

Leading-order metric (Sec. IV-B), Eq. (4.7):

"ds² = −dt² + a²Ψ⁴(dr² + r²dΩ²), where Ψ is given by a function of the radial coordinate r with spherical symmetry. This leading-order form of the metric commonly applies to the time-slicing conditions of the constant-mean-curvature slice, the comoving slice, and the uniform density slice for adiabatic systems. In this form, the radial coordinate is often called the isotropic coordinate, in which the spatial metric is described in a conformally flat form."

So γ_ij = a²Ψ⁴f_ij = a²e^{2ζ}f_ij at leading order; the brief's two forms are equivalent. Since ln Ψ(0) = μ/2, the peak of ζ = 2 ln Ψ is μ (derived from the stated definition, not written in the paper).

Spherical profile, Eq. (4.13): ln Ψ = (μ/2) exp(−k²r²/6).

Non-spherical profile (Sec. V), Eq. (5.1):

ln Ψ = (μ/2) exp[−(1/6)k²r²] [1 + (k²/6)(p(2X² − Y² − Z²) + 3e(Y² − Z²))], "where r² = X²+Y²+Z² and the ellipticity e and prolateness p characterize the non-spherical symmetry. This profile is realized as the typical profile [24] for the power spectrum: P(k̃) ∝ (3√6/√π)(k̃³/k³) exp(−3k̃²/(2k²)), (5.2) and reduces to Eq. (4.13) for e = p = 0."

No other definition of e or p is given. Their link to the deformation tensor is Eq. (5.13): (α̃, β̃, γ̃) = (2μ/15)(k²/(a²H²))(1 + 3e + p, 1 − 2p, 1 − 3e + p), where "α̃, β̃ and γ̃ are the eigenvalues of −∇D|_{r=0} with α̃ ≥ β̃ ≥ γ̃."

3. Horizon entry and time units (exhaustive) #

Eqs. (4.14)–(4.16), Sec. IV-C, verbatim:

"The radius of the maximum compaction function r_m can be given by (∂_r + r∂_r²) ln Ψ |_{r=r_m} = 0 ⇔ k²r_m² = 6. (4.14) Therefore, a typical scale is given by ar = √6 a/k. Estimating the horizon entry by using this typical scale, we find k²/(a²H²) |_ent = 6. (4.15) Then, since a³H³ = 2H_i²/(3t) = 50k²/(3t), the horizon entry time t_H is given by t_H = 100√6/k = 10√6 L. (4.16)"

Symbol: t_H. Numerically t_H ≈ 244.9/k ≈ 24.49 L. (With t_i = 2/(3H_i) = 2/(15k), t_H/t_i = 750√6 ≈ 1837; a(t_H) = 150 a_i. These ratios are my arithmetic from the paper's H_i, a_i.)

Time axes in figures:

• Fig. 3 (spherical dust, 80³), Fig. 5 (e = 0.1 dust), Fig. 7 (particle α₀), Fig. 10 (particle Hamiltonian constraint): x-axis label "cosmological time [t_H]". Fig. 3 and 5 run to 16 and 20 t_H; Figs. 7 and 10 run to 250 t_H.
• Fig. 8, 9, 11 (particle snapshots): "t = 102t_H, 122t_H, 163t_H, 176t_H, 204t_H, 245t_H"; Fig. 9 caption "Particle distribution at t = 245t_H."
• Fig. 1 (LTB null geodesics): axes "t × 10⁻³ [1/k]" and "r [1/k]" (the singularity t_s(0) for μ ≈ 0.32–0.33 sits near 2.05×10³/k ≈ 8.4 t_H, my reading).
• Fig. 12–13 (App. B, μ = 1.3, 40 grids): "t = 10L, 50L, 100L, 150L, 200L" (200L ≈ 8.2 t_H).

So "t = 250" in Fig. 7/10 means 250 t_H = 25000√6/k ≈ 6.1×10⁴/k ≈ 6.1×10³ L ≈ 4.6×10⁵ t_i (derived). The paper does not state how "cosmological time" is reconstructed from the coordinate time under the modified slicing (2.20). The only statement is for the simpler gauge (2.18): "In this gauge condition, at the boundary, since the value of α_c is fixed, the time coordinate is equivalent to the cosmological time of the background spacetime." For (2.20) they say it "initially resembles the e-folding number, and approaches to the cosmological time somewhat before the horizon entry."

Caution — two different "horizon entry" notions appear in the paper. Eq. (2.19) (gauge section) defines a_k = (k/(a_iH_i))^{−2/(1+3w)} a_i as "the background scale factor at the horizon entry of the scale 1/k" (i.e. k = aH, giving a_k = 25 a_i here), whereas Eq. (4.15) defines entry at k/(aH) = √6 (a = 150 a_i). t_H and all figure time units use the (4.15)/(4.16) definition.

4. Initial slice (Sec. III, verbatim) #

"K = −3H, (3.4) ψ = Ψ[1 − (1/6) q (1/(aH))²], (3.5) γ̃_ij = f_ij − (4/5) p_ij (1/(aH))², (3.6) Ã_ij = (2/5) p_ij H (1/(aH))², (3.7) α_c = 1 − (1/3) q (1/(aH))², (3.8) β^i = 0, (3.9) where q(x) = −(4/3) △Ψ/Ψ⁵, (3.10) p_ij(x) = (1/Ψ⁴)[−(2/Ψ)(D_iD_jΨ − (1/3)f_ij△Ψ) + (6/Ψ²)(D_iΨD_jΨ − (1/3)f_ij D^kΨD_kΨ)] (3.11) with △ := f^{ij}D_iD_j. The initial hypersurface is a constant-mean-curvature slice, since the value of K is constant. Once the geometrical variables are given, through the constraint equations (2.13) and (2.14), we can calculate E and J^μ. Then, we obtain ρ and V^μ from Eqs. (3.2) and (3.3). In the numerical simulations presented in this paper, we set H_i = 5k = 50/L and a_i = 1."

Also (3.1)–(3.3): T^{μν} = ρu^μu^ν, E = Γ²ρ, J^μ = EV^μ. Long-wavelength solutions attributed to Refs. [66, 67] (Shibata–Sasaki 1999; Harada, Yoo, Nakama, Koga 2015).

Note: the lapse actually initialized in the runs is Eq. (2.21), α = H_k/H_i − (H_k/(3H_i)) q (1/(a_iH_i))², i.e. (3.8) rescaled by H_k/H_i.

5. LTB section (Sec. IV-A, IV-B) — exhaustive #

Eq. (4.1): "ds² = −dt² + (∂_rR)²/(1 − k(r)r²) dr² + R(t,r)²dΩ², where k̃ is an arbitrary function of the radial coordinate r and the area radius R is a function of t and r." [Note: (4.1) prints k(r) but the text and all later equations use k̃; this k̃ is unrelated to the wavenumber k and to the k̃ in (5.2).] "The energy-momentum tensor is given by the dust fluid form as Eq. (3.1) with u^μ = (∂_t)^μ." Eq. (4.2): (∂_tR)² = −k̃r² + m(r)r³/(3R). Eq. (4.3): ρ = ((∂_r m)r³ + 3mr²)/(24πR²∂_rR). Eq. (4.4): M = (1/6) m r³ (Misner–Sharp). Eq. (4.5): R = r m^{1/3}(t − t_B(r))^{2/3} S(x), "where t_B is an arbitrary function of r, x = k̃((t − t_B)/m)^{2/3}" and S(x) is given by (4.6) (cosh/cos parametric forms), "S(0) = (3/4)^{1/3}. The function S(x) is analytic for x < (π/3)^{2/3}."

Sec. IV-B, verbatim:

"There are three arbitrary functions in the expression of the LTB solution. In this paper, we focus on the case t_B(r) = 0 since inhomogeneity induced by the nontrivial functional form of t_B(r) corresponds to decaying modes. Then, there are two functional degrees of freedom k̃(r) and m(r). These degrees of freedom correspond to growing modes and the gauge degree of freedom associated with the choice of the radial coordinate."

"rΨ²(r) = lim_{t→0} R(t,r)/a(t) (4.8) Ψ⁴(r) = lim_{t→0} (1/(1 − k̃r²)) (∂_r(R(t,r)/a(t)))² (4.9) We choose the background FLRW scale factor a = a_i(t/t_i)^{2/3} with a_i = a(t_i), where the time t = t_i is identified with the initial time considered in Sec. III. Thus we set a_i = 1, consistent with the parameter choice stated at the end of Sec. III. Noting lim_{x→0} S(x) = S(0) = (3/4)^{1/3}, from the first equation (4.8), we find m(r) = 4Ψ⁶/(3t_i²). (4.10) Taking the derivative of the first equation (4.8) with respect to r, we obtain Ψ² + 2rΨ∂_rΨ = ∂_r lim_{t→0} R(t,r)/a(t). (4.11) Substituting this equation into (4.9) and swapping the order of the r-derivative and the t → 0 limit that is guaranteed for long-wavelength solutions, we obtain k̃(r) = (1/r²)[1 − (1 + 2r∂_r ln Ψ)²]. (4.12) Therefore, the functional forms of k̃(r) and m(r) are uniquely given in terms of Ψ(r) in the present setting."

(I checked both algebraically: (4.8) with (4.5) at t_B = 0 gives m = 4Ψ⁶/(3t_i²); (4.9) gives 1 − k̃r² = (1 + 2r∂_r ln Ψ)².) App. B confirms usage: "The functional forms of k(r) and m(r) are fixed by setting the functional form of Ψ(r) through Eqs. (4.10) and (4.12)."

6. LTB nakedness results (Sec. IV-C) — exhaustive #

Method, verbatim:

"Now we are interested in whether the spacetime is globally naked or locally naked [71]. The spacetime singularity is located at R = 0, indicating t = 0 in the past and t = t_s(r) := (π/3) m k̃^{−3/2} after the time evolution. Since it can be shown that any future-directed outgoing geodesics cannot be emanated from the singularity t = t_s with r ≠ 0 (see, e.g., Ref. [72]), the singularity can be naked only at (t, r) = (t_s(0), 0). In general, multiple future-directed radial null geodesics can be emanated from this point. Therefore, we need to search for the most past null geodesic emanated from the point r = 0 and t = t_s(0). In practice, we solve the past-directed inward null geodesic equation from the radius (t, r) = (t_t, 1/k) with t_t being the trial value, and find the critical time t_cr above which the null geodesic hits the singularity. Then we solve the future-directed outgoing null geodesic equation from (t, r) = (t_cr, 1/k). Since the current density profile is a compensated profile, namely, the central over-dense region is surrounded by an under-dense region, once the outgoing null geodesic reaches the under-dense region, we determine the spacetime is globally naked; otherwise, the null geodesic is terminated at the singularity, and it is locally naked. We may safely say PBH forms if the singularity is locally naked but not globally naked [35, 73]."

Result, verbatim:

"In Fig. 1, we show the trajectories of the singularity t = t_s(r) (dashed line) and the most past null geodesic emanating from (t, r) = (t_s(0), 0) (solid line) for μ = 0.32 (blue) and μ = 0.33 (red). Since the solid and dashed red lines have two intersection points, the most past null geodesic emanating from (t, r) = (t_s(0), 0) terminates at the singularity. Therefore, the spacetime structure is locally naked, and the singularity is hidden by the event horizon from the asymptotic observers for μ = 0.33. On the other hand, since there are null geodesics emanating from (t, r) = (t_s(0), 0) and can reach the asymptotic infinity for μ = 0.32, the spacetime structure is globally naked."

Null geodesic eqs: (4.17) −ṫ² + (∂_rR)²/(1 − k̃r²) ṙ² = 0; (4.18) ẗ + (∂_rR)(∂_t∂_rR)/(1 − k̃r²) ṙ² = 0.

Figures: Fig. 1 (t_s(r) and most-past null geodesic, r ∈ [0, 2.5]/k, t ∈ [2, 5]×10³/k; no apparent horizon drawn). Fig. 2 (schematic Penrose-like diagrams; "Blue curves describe the trajectories of apparent horizons"; left panel annotated "Globally Naked μ ≲ 0.32", right "Locally Naked μ ≳ 0.33"). Fig. 6 caption: "The red arrows indicate the threshold between the globally naked (downward) and locally naked (upward) cases of the LTB solution." Fig. 13 (App. B): compactness 2M/R vs proper length; "the middle [intersection with 2M/R = 1] is the marginally outer trapped surface, which indicates the formation of a black hole."

Shell crossing: the LTB analysis considers only the shell-focusing singularity R = 0 at t_s(r); no check for ∂_rR = 0 (shell crossing) before t_s(0) is mentioned anywhere in Sec. IV. NOT FOUND.

The criterion is about the event horizon, not the apparent horizon: "the singularity is hidden by the event horizon from the asymptotic observers for μ = 0.33." The spherical 3D dust runs (Sec. IV-D, μ = 0.30–0.50, 80³) all find an apparent horizon at ≈ 11–12 t_H (Fig. 3, my reading); for μ = 0.3 "an apparent horizon is found in relatively late times. At the formation time, the size is too small, and the constraint is significantly violated outside the apparent horizon. This behavior may reflect the global nakedness of the spacetime, although the difference from the cases which have slightly larger initial amplitudes is not distinct."

7. Analytic threshold (Sec. V-A) #

Eq. (5.12): h(α̃, β̃, γ̃) := (2/π) (α̃ − γ̃)/α̃² · E(√(1 − ((α̃ − β̃)/(α̃ − γ̃))²)) "In Ref. [29], from the hoop conjecture [74], the criterion for PBH formation is given by h(α̃, β̃, γ̃) < 1 at the horizon entry." Eq. (5.13): (α̃, β̃, γ̃) = (2μ/15)(k²/(a²H²))(1 + 3e + p, 1 − 2p, 1 − 3e + p) Eq. (5.14): μ > ĥ(e, p) × 6a²H²/k² Eq. (5.15): ĥ(e, p) := (15/π) · e/(1 + 3e + p)² · E(√(1 − (e+p)²/(4e²))) "The contour map of ĥ(e, p) is depicted in Fig. 4." Fig. 4 caption: "The shaded region is excluded by the condition α̃ ≥ β̃ ≥ γ̃ ≥ 0." Contour labels visible: 0.1, 0.2, 0.3, 0.4, 0.42, 0.44, 0.46, 0.47, 0.48.

Also (5.11): D ≃ (4/5)(1/(a²H²)) ∇ ln Ψ, and (5.10) Ψ_lin ≃ ξ ≃ Ψ − 1 ≃ (5/6) ln Ψ? — the extraction reads "Ψ_lin ≃ ξ ≃ Ψ − 1 ≃ (5/6)?" ambiguous; I did not rely on it.

At entry (aH = k/√6) the criterion is μ > ĥ(e,p); Fig. 6's blue line is "μ = ĥ(e, 0)". No tabulated values appear. Evaluating (5.15) (E with modulus √(1−1/4) for p = 0, i.e. parameter m = 0.75): ĥ(0.01,0) = 0.055, ĥ(0.05,0) = 0.219, ĥ(0.1,0) = 0.342, ĥ(0.2,0) = 0.452, ĥ(0.25,0) = 0.472, ĥ(0.3,0) = 0.481 — consistent with the blue line in Fig. 6 and with the brief's 0.34 / 0.45. The text remarks "the line of the analytic criterion approaches μ = 0 at e = 0", consistent with ĥ ∝ e at small e (slope (15/π)E ≈ 5.78, my derivation); no "μ_th ≈ e S(t)" expression exists in the paper.

8. 3D dust results (Sec. IV-D, V-B) #

"We use 80 grids in each direction. The formation of the apparent horizon is monitored during the simulation." (IV-D) / "Here, let us focus on the p = 0 cases. We use 80 grids for each direction, as with the spherically symmetric cases." (V-B) "The results can be classified into three cases: stable evolution with horizon formation ("HF w/o crash"), crash after horizon formation ("crash w/ HF"), and crash without horizon formation ("crash w/o HF"). The smallest value of μ in the case "HF w/o crash" may give a sufficient condition for PBH formation. The real threshold may exist below the value giving the sufficient condition." "We perform the numerical simulations to find the classification of the cases for e = 0, 0.01, 0.03, 0.05, 0.1, 0.15, 0.2, and 0.25." Fig. 5 (e = 0.1, p = 0) legend: μ = 0.750 HF w/o crash; 0.725 HF w/o crash; 0.700 crash w/ HF; 0.675 crash w/ HF; 0.650 crash w/o HF; 0.600 crash w/o HF. (Horizons appear at ≈ 12–13 t_H in Fig. 5, my reading.) Sec. VI-C: "As an example, let us focus on the case e = 0.2. According to Fig. 6, the calculation breaks down before horizon formation for μ ≤ 0.95." "Although we cannot give a lower bound of the threshold for μ, since the analytic estimation with aH = k/√6 gives obviously smaller values compared to the boundary between the green triangles ("crash w/ HF") and the red crosses ("crash w/o HF"), the real threshold might be larger than the values indicated by the blue solid line." "Our results suggest that the threshold value might be significantly larger, so PBH formation would be harder than the analytic estimation within the range in which the fluid approximation is valid."

Fig. 6 (my reading of the markers): e = 0.2: red crosses up to ≈ 0.95–1.0, green triangles ≈ 1.0–1.05, blue circles ≈ 1.1; e = 0.25: red ≈ 1.05, green 1.1–1.4, blue ≈ 1.5; e = 0: red arrows at ≈ 0.32–0.33 with green triangles from ≈ 0.1 up and blue circles at ≈ 0.45–0.5.

Shell-crossing attribution, verbatim: Intro: "for a dust fluid, we suffer from the crash of the numerical computation associated with the shell-crossing singularity, which generically happens due to the intersection of fluid elements. Therefore, we cannot probe the spacetime after the moment of the shell crossing with the dust fluid description." Summary: "The breakdown of the simulation is caused by shell focusing or shell-crossing singularities that appear in the dust fluid description." Footnote 1: "The constraint violation occurs around when and where the particles intersect and cross each other. This situation corresponds to the shell crossing or shell focusing singularity for the dust fluid case."

9. Particle (collisionless) results — exhaustive #

Method (Sec. VI-A): geodesic 3+1 equations (6.2)–(6.4) from Ref. [75]; stress tensor (6.5) T^{μν} = −Σ_p m_p δ³(x − x_p)/(u_p^λ n_λ √γ) u_p^μ u_p^ν; (6.6) E = Σ_p m_pΓ_p δ³(x−x_p)/√γ; (6.7) J^i = Σ_p m_pΓ_pV_p^i δ³/√γ; (6.8) S^{ij} = Σ_p m_pΓ_pV_p^iV_p^j δ³/√γ. Deposition kernel (6.9), verbatim:

"Since the delta functions cannot be treated numerically, for each particle, we instead assign a cubic domain based on the coordinate spacing and consider uniform distributions inside the cubic box. Namely, the delta function is replaced by the uniform support function as δ³(x − x_p) → Θ(Δx/2 − |x − x_p|) Θ(Δy/2 − |y − y_p|) Θ(Δz/2 − |z − z_p|) · 1/(ΔxΔyΔz), (6.9) and √γ is evaluated at x = x_p."

Particle settings (Sec. VI-B), verbatim in full:

"In order to convert the fluid distribution into a particle distribution, we consider a regularly aligned initial distribution of particles. In addition, we assume the numerical region is initially filled with no space and no overlap between particles. Since the particles are properly aligned initially, we cannot introduce inhomogeneity by an inhomogeneous particle distribution. Instead of an inhomogeneous distribution, we introduced different proper masses for individual particles so that the density distribution may have the desired inhomogeneity. That is, we set each particle mass m_p as m_p = (E(x_p)/Γ_p) √γ(x_p) ΔxΔyΔz (6.10) on the initial hypersurface."

Number of particles per cell and total N: not stated. Fig. 7 caption: "Time evolution of α₀ for each value of μ with e = 0.2 and N = 80." (N is introduced in VI-A as the particle count, but 80 is evidently the grid number per direction; the paper does not resolve this.) Fig. 9 caption: "The number of particles has been reduced to 1/80" (for plotting only). Initial velocities: no explicit statement in VI-B; App. B calls the spherical test "a spherically symmetric system of collisionless particles free from velocity dispersion".

Results (Sec. VI-C), verbatim:

"First, as is described in Appendix B, we performed the comparison between the numerical simulation of the spherical case and the corresponding LTB solution for μ = 1.3. We can find good agreement for this parameter setting until the time of horizon formation. In this section, we are mainly interested in the situation where the simulation breaks down for the dust fluid case. As an example, let us focus on the case e = 0.2. According to Fig. 6, the calculation breaks down before horizon formation for μ ≤ 0.95. However, for the particle simulation, the calculation does not crash even for much smaller values of μ. In Fig. 7, we show the value of α₀ as a function of the cosmological time for each value of μ. We also plot the case for e = 0 and μ = 0.3 as a reference. As is shown in Fig. 7, the value of α₀ starts to steeply decrease from a certain time and an apparent horizon is formed for μ ≥ 0.05. In contrast, for μ ≤ 0.045, the decrease of α₀ stops, and it keeps a middle value without horizon formation. Let us check the overall dynamics following snapshots of the density profile. It can be found that the density perturbation δ = 8πE/(3H²) − 1 first increases around the center and the tip of the over-dense distribution along the longest distribution axis, which is the y-direction in the present case. This behavior around the tip of the over-dense region is similar to the behavior reported in Refs. [77, 78]. In the case μ = 0.045, the collapse halts, and it may describe a halo formation below the threshold amplitude. This halo would be supported by the effective pressure generated by the velocity dispersion of the particles. To explicitly see that, we plot the particle distribution in Fig. 9 together with the arrows presenting the velocity of each particle. One can find that the velocities of the particles are not coherent in the central region, and velocity dispersion is associated with the random motion." "The calculation does not break down for μ ≤ 0.045 and will continue to run unless manually stopped. If one accepts this result, the threshold of PBH formation is at around μ ∼ 0.045, which is much smaller than the values indicated by Fig. 6 and even smaller than the analytic estimation (5.14) with aH = k/√6. However, there is a big caveat. Although the calculation does not break down, constraints are significantly violated, as shown in Fig. 10 and 11. The origin of the localized numerical constraint violation requires further investigation and limits the quantitative interpretation of the particle results.¹ Nevertheless, it would be suggestive enough showing an example of the numerical simulation in which the mean absolute value of the Hamiltonian constraint is small enough (see the lower panel of Fig. 10). That is, the violation of the constraint equations remains localized."

Footnote 1 (verbatim):

"The constraint violation may be regarded as a mismatch between the stress-energy tensor described by the particle system and the geometry. Then we might expect that this mismatch originates from the failure of the collisionless particle description and can be resolved by introducing more physically realistic matter fields. The constraint violation occurs around when and where the particles intersect and cross each other. This situation corresponds to the shell crossing or shell focusing singularity for the dust fluid case. Since the matter density and curvature are divergent at the shell crossing and shell focusing singularity, a more realistic description of the matter field than dust fluid and particle system would be needed to properly describe the small-scale dynamics around the high-density regions."

Runs in Fig. 7/10 (e = 0.2 unless noted): μ = 0.300 (e = 0.0, reference), 0.300, 0.150, 0.100, 0.075, 0.060, 0.050, 0.045. Fig. 7 plots the unnormalized α₀ (values ≈ 9–13.5 before collapse, in contrast to α₀/α_L in Figs. 3, 5). Horizon-formation times are not stated in the text. From the α₀ drops in Fig. 7 (my reading): μ = 0.30, e = 0: ≈ 12 t_H; μ = 0.30, e = 0.2: ≈ 22 t_H; 0.15: ≈ 47 t_H; 0.10: ≈ 73 t_H; 0.075: ≈ 108 t_H; 0.060: ≈ 148 t_H; 0.050: staircase drop from ≈ 180 to ≈ 215 t_H; 0.045: dips to ≈ 8 near 200 t_H then flattens at ≈ 9 through 250 t_H. Fig. 10 (upper): max-norm Hamiltonian violation rises to O(1) by ≈ 50–100 t_H for all runs and reaches ≈ 10 near the horizon-formation times; (lower) mean absolute value stays ≲ 2×10⁻³. Why the μ ≥ 0.05 runs terminate shortly after the α₀ drop is not stated (no excision statement for particle runs; excision is mentioned only for the spherical dust μ = 0.5 case).

Gauge (Sec. II-B), verbatim:

"(∂_t − β^i∂_i)α_c = −2α_c(K + 3H), (2.18)" — described as "the following modified version of the 1+log slicing condition is often used". "we adopt the following time slicing condition: (∂_t − β^i∂_i)α = −2α(K + 3H) + (3/2) αH_k exp(−C²α_b²/a_k²), (2.20) where α_b is the value of the lapse function at the corner of the numerical box x = y = z = L, and we fixed the constant C as 5 after some numerical trials. The initial profile of α is set as α = H_k/H_i − (H_k/(3H_i)) q (1/(a_iH_i))². (2.21)" Shift: "(∂_t − β^i∂_i)β^i = (3/4)B^i, (2.22); (∂_t − β^i∂_i)B^i = ∂_tΓ̃^i − 3HB^i, (2.23) where Γ̃^i = −D_jγ̃^{ij}." a_k = (k/(a_iH_i))^{−2/(1+3w)} a_i (2.19); "1/H_k … the timescale at the horizon entry … characterized by the scale of the inhomogeneity k and is time-independent." Ref. [65] (Ning, Cai, Wang, Yoo 2026) cited for "a similar approach".

Box and coordinates (Sec. II-C): "we consider the region 0 ≤ X^i ≤ L"; non-Cartesian coordinates X^i = x^i − (η/(1+η))(L/π) sin(πx^i/L) (2.24), "first implemented in Ref. [56]", "The value of η is set to 10." L = 10/k follows from "H_i = 5k = 50/L". Boundary conditions: only "This functional form is compatible with the boundary conditions adopted in this work and satisfies x^i = 0 at X^i = 0 and x^i = L at X^i = L." — never stated as periodic/reflective.

App. B (particle vs LTB test): μ = 1.3, "we introduce only 40 grid points for each direction in this demonstration. Even for this low resolution, we can find good agreements, which clearly demonstrate the robustness of the simulation as long as the fluid description remains valid. At the time t = 200L, we could not find an apparent horizon due to the low resolution." No other resolution study exists.

10. Matter-era duration / velocity dispersion #

Reheating appears only in the Intro: "Such early matter-dominated phases can arise from the oscillation of inflaton or other scalar fields after inflation before reheating completes." No discussion of a finite matter era, deadline, or dependence of the threshold on its duration. Velocity dispersion: Intro cites "[36] for effects of inhomogeneity and velocity dispersion"; the only other discussion is the qualitative halo remark quoted in §9. No parametric study.

11. Citations #

• Carr, Iovino, Perna, Vaskonen, Veermäe 2026, arXiv:2601.06024 = Ref. [12], "Riv. Nuovo Cim. 49, 225 (2026)". Cited once, in the Intro list "Useful reviews in different aspects and a summary of observational constraints can be found in, e.g., Refs. [4–12]." No asteroid-mass-window statement anywhere.
• Shapiro & Teukolsky 1991 (PRL 66, 994) = Ref. [77], with Yoo–Harada–Okawa 2017 [78]: "This behavior around the tip of the over-dense region is similar to the behavior reported in Refs. [77, 78]." (Sec. VI-C)
• Harada, Yoo, Kohri, Nakao, Jhingan 2016, arXiv:1609.01588 = Ref. [29]: Intro ("[29–32]"; "ellipticity [29]"); Sec. V-A the Zel'dovich/hoop criterion "following Ref. [29]", h defined "Following Ref. [29] (see also Ref. [31])", "In Ref. [29], from the hoop conjecture [74]…".
• Harada, Kohri, Sasaki, Terada, Yoo 2023, arXiv:2211.13950 = Ref. [36]: Intro "Ref. [36] for effects of … velocity dispersion"; also in "[29–36]" analytic works "with not necessarily valid perturbative approximations and hypotheses".
• Kokubu, Kyutoku, Kohri, Harada 2018, arXiv:1810.03490 = Ref. [35]: Intro "[33–35] … effects of inhomogeneity"; Sec. IV-C "We may safely say PBH forms if the singularity is locally naked but not globally naked [35, 73]."
• de Jong, Aurrekoetxea, Lim 2021 (arXiv:2109.04896) = [57]; de Jong, Aurrekoetxea, Lim, França 2023 = [58]: Intro "3+1 dimensional numerical investigation into black hole systems in an expanding background has been developing rapidly in recent times [25–28, 57–61]."
• Ebrahimian et al.; Milligan et al. 2025: NOT CITED (no match in the 78-item reference list).
• Other relevant: [30] Harada+2017 spins, [31] Saito+2024, [32] Ye+2025 (spin discussion in Summary); [62] Munoz & Bruni 2023 (dust in NR); [65] Ning, Cai, Wang, Yoo 2026 (slicing); [46] Escrivà+2026 N-body GW; [61] Escrivà 2026.

12. Code #

No code name appears (no "COSMOS" or any other name), no statement of public availability, no timing/cost figures. The only cost-related sentence: "An efficient slicing condition that can save computational time may be given by …" (Sec. II-B). Acknowledgements list only JSPS KAKENHI grants JP25K07281, JP24K07027, JP26K17141.

13. Conclusions and outlook (Sec. VII, verbatim) #

"We have investigated the formation of primordial black holes in an early matter-dominated universe using fully nonlinear 3+1 numerical relativity. We constructed initial data for the dust fluid using a long-wavelength curvature perturbation with ellipticity, allowing us to study both spherical and triaxial collapse configurations. After reviewing spherically symmetric dust collapse and the connection with long-wavelength solutions, we examined the limitations of the dust fluid description in a matter-dominated background. The dust fluid simulations show that sufficiently large amplitude perturbations can form apparent horizons, while smaller amplitudes fail to form a black hole before the simulation breaks down. The breakdown of the simulation is caused by shell focusing or shell-crossing singularities that appear in the dust fluid description. To address this breakdown of the fluid approximation in a concrete collisionless realization, we also implemented a collisionless particle system in full numerical relativity. The particle-based simulations avoid the immediate crash associated with fluid shell crossing, and they demonstrate horizon formation for amplitudes larger than a certain value. The resulting collapse threshold is significantly smaller than earlier analytic estimates based on simplified matter-dominated criteria, although the late-time evolution is accompanied by significant Hamiltonian constraint violations after the particle intersection happens. If one interprets the failure of the fluid approximation as indicating that singularity formation associated with shell focusing or shell crossing suppresses PBH formation, then the threshold may be significantly larger than analytic estimates, and PBH formation may be harder to realize. On the other hand, if one accepts the particle simulation results shown here, PBH formation becomes much easier than previously expected. Thus, our results forcefully suggest that there is a large theoretical uncertainty in the PBH formation criterion during a matter-dominated epoch. Here, we should remark on the impeding mechanism for PBH formation in a matter-dominated era associated with black hole spin. Throughout the analyses in Refs. [30–32], the conclusion has evolved, but the latest work suggests that spin effects are negligibly small compared with the non-spherical collapse suppression mechanisms studied here. However, these spin estimates were obtained using analytic models based on the Zel'dovich approximation and the hoop conjecture, so the change of the threshold value in the non-spherical collapse studied numerically in this paper may also play an important role in the spin effect. Then continued analyses of spin effects remain necessary. Future work should refine the treatment of the high-density regime, improve constraint preservation, and extend the analysis to quantify the PBH abundance and the possible gravitational wave signatures associated with matter-dominated collapse. These developments are necessary to make definitive predictions for PBH formation in realistic early-universe scenarios."

14. Points where the paper differs from (or is weaker than) the brief's reading #

1. Gauge: the brief's condition (∂_t − β^i∂_i)α = −2α(K+3H) is Eq. (2.18), which the paper calls the commonly used form; the runs use Eq. (2.20) with the additional (3/2)αH_k exp(−C²α_b²/a_k²) term (C = 5) and a rescaled initial lapse (2.21). Any reproduction that uses (2.18) is not the paper's gauge.
2. Resolution: consistent with "80³ single resolution" — the only other resolution anywhere is the 40³ spherical particle test (μ = 1.3) in App. B; no convergence study for any 3D e ≠ 0 dust or particle run. The Fig. 7 caption "N = 80" is ambiguous (N is defined as the particle number in VI-A).
3. Horizon-formation times for the particle runs are never given numerically; the only source is Fig. 7, where the μ = 0.05 drop occurs at ≈ 180–215 t_H and μ = 0.045 shows no drop up to 250 t_H. The brief's "100–250 t_H" is a plot-reading range, not a stated result.
4. LTB numbers (0.32 globally naked / 0.33 locally naked) are exactly as in the brief; but note the criterion is event-horizon based, and no shell-crossing check is reported.
5. Two "horizon entry" definitions coexist: (2.19) uses k = aH (a_k = 25 a_i) for the gauge; (4.15)/(4.16) use k/(aH) = √6 for t_H. All "[t_H]" axes refer to the latter.
6. Particle count, particles per cell, boundary conditions, code name, cost: none stated.
7. ĥ tabulated values and any small-e formula: not in the paper; the brief's 0.45 and 0.34 are reproduced by evaluating (5.15) and agree with the Fig. 6 blue line, but they should be cited as computed from (5.15), not as quoted.
8. Abstract's "factor of 2" (dust) and "order of magnitude" (particles) refer to comparison with the analytic ĥ(e,0) line at aH = k/√6 (≈ 0.45 at e = 0.2 vs. dust HF at μ ≳ 1.0 and particle threshold μ ≈ 0.045–0.05).

Local copies of the extracted text and figure-page images are in /tmp/claude-0/-root-PBH/a34e6be4-0ee0-4598-bfd0-88626255d779/scratchpad/ (yoo_raw.txt, yoo.txt, pg-10.png … pg-19.png); the PDF itself is at /root/.claude/projects/-root-PBH/a34e6be4-0ee0-4598-bfd0-88626255d779/tool-results/webfetch-1790190337945-w8z33i.pdf.