Verification report: PBH formation from an oscillating scalar field (checked 2026-09-23)

Sources were downloaded and grepped locally (arXiv HTML with LaTeX alt-text, plus PDF→text) under /tmp/claude-0/-root-PBH/a34e6be4-0ee0-4598-bfd0-88626255d779/scratchpad/papers/ (2504.02600v1/v2, 2109.04896v2, 2306.11810v2, 2403.02878v1, 2501.13046v1/v2). Quotes below are verbatim from those files (math is the source's LaTeX alt-text). IOPscience is behind a bot-wall, so journal refs come from the arXiv pages / GRTresna reference list.

Ledger #

2025–2026 follow-ups (same groups) #

No paper by either group on PBH formation from a directly perturbed oscillating scalar in 3D was found. Closest items:

Verbatim quotes #

Paper 1 — arXiv:2504.02600 (v1 3 Apr 2025; v2 26 Aug 2025, "21 pages, 9 figures. Accepted for publication in JCAP") #

(a) Title/authors/abstract (v2). Title: "Primordial Black Hole Formation in a Scalar Field Dominated Universe". Authors: Ethan Milligan, Luis E. Padilla, David J. Mulryne, Juan Carlos Hidalgo. Abstract:

"We present a numerical code that solves the Misner-Sharp system for a spherically symmetric cosmological model containing both a scalar field and a perfect fluid. While the code is capable of exploring general scenarios involving a minimally coupled scalar field and perfect fluid, we focus on the regime where the scalar field dominates the dynamics, particularly in the post-inflationary scalar field-dominated scenario, where the universe is governed by a rapidly oscillating scalar field for a period lasting a few e-folds. We analyse the threshold for PBH formation under quadratic and quartic potentials, evolving configurations from superhorizon scales. Our results confirm that a quartic potential behavior is similar to the radiation-dominated universe, resulting in a PBH formation threshold close to the well-established value in radiation backgrounds. Conversely, in the quadratic case, we observe a significant deviation from the expected dust-like behaviour, due to wave-like effects opposing the gravitational collapse. While numerical limitations prevent us from evolving a wide range of initial conditions to determine a precise threshold for PBH formation, our findings suggest that PBH formation may be suppressed with respect to the pure dust scenario, allowing the formation of stable solitonic structures instead. This study highlights the importance of properly accounting for wave dynamics in oscillating scalar fields when characterising PBH formation."

(v1 abstract identical except "an minimally coupled".)

(b) Setup.

• Sec I, Eq. (1): "V(φ)=λ_n/(2n) φ^{2n} ... where λ_n is a free parameter of the model under consideration. In the case n=1, λ_1≡μ² can be interpreted as an effective mass parameter, whereas for the case n=2, the parameter λ_2≡λ is interpreted as a self-interaction term."
• Sec I: "we hereby introduce a code with the capability of exploring general cosmological scenarios involving the combination of a scalar field minimally coupled with a perfect fluid. Specifically, we extend the numerical formalism presented in [76] for the perfect fluid case and produce a numerical code that solves the Misner-Sharp system for a universe containing both a scalar field and a perfect fluid."
• Sec II, Eq. (3): "In this paper, we use the Misner-Sharp formalism [74], which uses the following line element of the spacetime ds²=−e^{2φ}dt²+e^{λ}dA²+R²dΩ². The variables φ, λ and R are functions of A and t only."
• Sec III.3.2: "Given that our code is set in a coordinate system that moves with the perfect fluid, we cannot, in general, simply set the fluid contributions to zero. Instead, we consider an extremely diluted, homogeneous, and isotropic fluid".
• Sec IV: "When simulating the scalar field collapse, we adopt a radiation-like reference fluid with an equation of state parameter w_pf=1/3 ... we consider a subdominant radiation component, maintaining an energy density ratio of \tilde{\rho}{\varphi}/\tilde{\rho}{pf}\simeq\mathcal{O}(10^{8})."
• Sec IV: "We numerically integrate the evolution equations along the logarithmic time coordinate ξ using 4^{th} order accurate Runge-Kutta method. The spatial derivatives are performed using \mathcal{O}(\Delta\bar{A}^{4}) finite difference methods ... our code contains regridding procedures to adaptively increase the number of grid points by interpolating data ... We stabilise the method by implementing a N=3 Kreiss-Oliger [80] dissipation".
• Code name / public availability: NOT FOUND (no name, no URL, nothing in Acknowledgments).

(c) Mass vs Hubble.

• Sec V.1: "if the oscillation of φ is sufficiently undamped, which is the case if 2\mu\gg 3H, then the friction term in the Klein-Gordon equation can be neglected. Therefore, the background dynamics of the scalar field evolves according to (see appendix A): \tilde{\varphi}{\rm b}=\tilde{\varphi}{\rm 0,b}\cos\left(\frac{2}{3}\tilde{\mu}e^{\xi}+\delta_{0}\right). (56) We see that φ oscillates with angular frequency that grows exponentially in ξ."
• Appendix A, after Eq. (A7): "where \tilde{\mu}=R_{H}\mu."
• Sec VI.1: "In all our simulations, we set the self-interaction parameter to \tilde{\lambda}=10."
• Sec VI.2: "as the initial amplitude is reduced, the scalar field undergoes an exponentially increasing number of oscillations before collapse."
• Sec I: "the challenges associated with tracking scalar field oscillations over a number of Hubble times."
• Fig. 1 caption: "Numerical and analytical dimensionless frequency of the background universe in a quadratic dominated scenario."
• Numerical value of μ̃, m/H, H_i, or oscillations per Hubble time: NOT FOUND.

(d) κ range, horizons, bounds, reason.

• Sec VI.2 (v2): "we further analyze the evolution of different perturbations, considering various initial conditions for κ within the range \kappa\in{0.06,0.2}. As shown in Fig. 7, we observe that smaller initial amplitudes collapse at later times, with the time of collapse defined as the moment when the apparent horizon forms."
• Sec VI.2 (v2): "For initial amplitudes below \kappa=0.06, we did not obtain simulations sufficiently reliable, according to our consistency tests, and within the limitations of the computational power available. One of the key challenges in this regime is that, as the initial amplitude is reduced, the scalar field undergoes an exponentially increasing number of oscillations before collapse. This significantly increases the runtime and imposes stringent demands on both spatial resolution and numerical stability. In particular, tracking the highly oscillatory evolution of the scalar field over extended timescales complicates maintaining control over the Hamiltonian constraint. Despite applying regridding techniques and increasing the resolution (from N=800 up to N=10,000), the magnitude of our normalized \mathcal{H} eventually exceeds \sim 10^{-2} for very small amplitudes. This saturates the upper limit of acceptable numerical error. Therefore, in this work we restrict our analysis to the robust simulation regime \kappa\in[0.06,0.2], where we can ensure the numerical accuracy and physical reliability of our results."
• v1 wording of the same: "For initial amplitudes below 0.06, we have not obtained simulations with sufficient confidence in their physical validity. As the initial amplitude decreases, the runtime increases exponentially due to the correspondingly exponentially increasing number of scalar field oscillations. Our current study focuses on the robust simulation regime between initial amplitudes of \kappa=0.06 and \kappa=0.2, where we have high confidence in the results."
• v1: "In all cases we evolved, we observe similar behavior: the perturbation initially collapses, and at a certain point, the central density begins to flatten before an apparent horizon forms. This recurring feature indicates that, within the framework of a collapsing quadratic scalar field, we can establish only an upper bound on the threshold necessary for PBH formation, which is associated with the minimum κ value we evolved numerically. Specifically, our results indicate: \mathcal{C}{\rm th}\leq 0.2,\quad\delta{\rm th}\leq 0.077. (61)" — v2 identical except "Across this range, we observe similar behavior" and Eq. (62).
• Sec VI.2 (v2), after Eq. (62): "although our simulations appear to predict the formation of a central soliton + NFW envelope structure, we find that these configurations still collapse into PBHs even for perturbation amplitudes smaller than those predicted by Newtonian estimates. This implies that the second reported threshold value in Eq. (61) is not accurate enough to reliably estimate PBH formation abundances in these scenarios, allowing PBHs to form at lower amplitudes than previously expected." (Eq. 61: δ^{sol}{th}=0.019, δ^{sol+NFW}{th}=0.238, from Newtonian ref [68].)
• Fig. 9 caption: "Top-panel: Evolution of the maximum of the compaction function for the particular case of \kappa=0.06. Bottom panel: Final time slicing, with a vertical line marking the location where the outgoing radial null geodesic vanishes".

(e) Soliton / NFW / flattening.

• Sec VI.2: "we evolve the same type of initial conditions (with \kappa=0.2) for both dust-like and scalar field perturbations. In Figure 5 ... although the scalar field density initially appears to evolve similarly to dust perturbations (see the top panel of the figure), the bottom panel reveals that, once the quantum pressure kicks in, the central region of the density in the quadratic case becomes flattened."
• Sec VI.2: "the central region of the density profile can be well approximated by the solution of a central soliton (depicted by the red dotted-dashed line), with deviations not exceeding 2% when compared to the analytical soliton profile. For this soliton, we use the profile reported in [87], given by: \rho(R)=\frac{\rho_{c}}{(1+0.091(R/R_{c})^{2})^{8}}, (60) ... Moreover, we observe that the decline of the central profile does not follow that of the soliton for large values of R. In fact, the profile decreases more gradually, transitioning from a falloff of \propto R^{-3} (which closely resembles the NFW decay at large radii) [fn 3: The transition from the soliton core to the polynomial envelope was set at R\simeq R_{c}, in accordance with the behavior observed in Newtonian simulations.] to a subsequent power-law decline \propto R^{-1.78} with oscillations, before asymptotically settling into the cosmological background. Let us emphasize that, to our knowledge, although the formation of solitons with NFW-like envelopes has been previously reported as an outcome of cosmological scalar field collapse (both, in spherical and 3D symmetries), this is the first time such a structure has been reported within the full general relativistic description."
• Ref [87]: "H.-Y. Schive, T. Chiueh, and T. Broadhurst, 'Cosmic structure as the quantum interference of a coherent dark wave,' Nature Physics, vol. 10, pp. 496–499, July 2014."

(f) Initial data, ε, duration.

• Sec III.3: "We require the initial conditions to be cosmologically consistent with being generated from inflationary perturbations. In particular, the initial perturbations should consist solely of the growing component on scales larger than the cosmological horizon." ... "The method for constructing initial conditions that selects the growing component in the linear regime was first developed in [48]. This formalism was clarified and further developed for the perfect fluid dominated case by Bloomfield et al [76]."
• Sec III.3.1: "for the scalar field-dominated case, we truncate our initial conditions to first-order as the complexity of the resulting system of equations makes higher-order corrections impractical".
• Sec III.3.2: "we start by adopting the approach of generating appropriate cosmological initial conditions based on an initial curvature profile [51, 42, 43], which was claimed to pick correctly the growing mode." ... "To construct the initial data, we assume that at superhorizon scales, the curvature profile is a time-independent quantity, K(A,t)=K(A). Under this assumption, we expect that both \delta_{U} and \delta_{m_{\varphi}} have a time-dependence of \propto\text{exp}[2(2n-1)\xi/3n]." ... "We simplify our description by assuming that at \xi=0 (the initial time in our numerical simulations), we have \tilde{\varphi}{\rm b0}=0 and then \tilde{\Pi}{\rm b0}^{2}/2=1. Our major simplification will be to assume that at time \xi=0, {\tilde{\varphi}{0}}=\tilde{\chi}{0}=0, so all the information of the initial profile that we will use is contained in the kinetic term of the scalar field." (Eqs. 48–51 give δ_{Π0}, δ_{U0}, δ_{R0} in terms of δ_{m0φ}.)
• Sec IV, Eq. (54): "\delta_{m_{T}}=\kappa e^{-A^{2}/2\sigma^{2}}, where, for simplicity, we used in all our simulations the value \sigma=2R_{H}, [fn 2: Note that this value is just slightly larger than the cosmological horizon. We can make this simplification due to the fact that in the previous section we showed that our initial conditions correctly select only the growing mode of the solutions.] leaving κ as the only free parameter in our simulations." Also: "In all the simulations, we fix the value a(t_{0})=1, \xi_{0}=0 defined by t=t_{0}e^{\xi} ... Automatically, we obtain R_{H}=t_{0}/\alpha and H_{0}=\alpha/t_{0}."
• ε = k/aH: NOT FOUND (in Sec III.3.1 "ε is an order counting parameter for the derivative expansion different from ϵ denoting perturbative ordering").
• Duration: Fig. 7 caption only: "Time interval from horizon entry until apparent horizon formation, as a function of different initial perturbation amplitudes κ. The left vertical axis displays logarithmic time \xi_{AP} and the right shows the number of e-folds N_{AP}." Radiation test (Sec V.3): "evaluated at the time of horizon crossing, \xi_{\rm HC}=1.25".

(g) Convergence / constraints.

• Sec V.1: "we compute the relative errors of different variables with respect to the background analytical FRW solutions ... we compute the L^{2} norm at each time step" ... "The background dynamics of the scalar field are successfully replicated numerically by our code as presented in Fig 1, where the relative error between the numerical and analytical solutions remains below 10^{-3}."
• Sec V.2: "A key diagnostic for the reliability of our simulations is the consistency with the Hamiltonian constraint, Eq.(28j). Since we do not evolve \bar{\Gamma}^{2} using an independent evolution equation, but rather compute it iteratively at each time step from the constraint itself and then use it within the evolution system, we cannot directly compare it to a dynamically evolved quantity. Instead, to assess whether the constraint is being consistently satisfied, we differentiate Eq.(28j) with respect to the spatial coordinate \bar{A} and monitor whether the resulting relation holds throughout the simulation." ... "we define a normalized version of \mathcal{H} and consider the simulation to be numerically consistent if the magnitude of the Hamiltonian constraint remains below 10^{-2} throughout the evolution."
• Sec V.3: "the threshold value for PBH formation in our simulations is \kappa=0.1735, corresponding to a maximum compaction function value of \mathcal{C}_{\rm th}=0.50 ... Our results closely match Figure 12 in Ref. [76]."
• Resolution-convergence study (Richardson etc.): NOT FOUND.

(h) Conclusions and future work (Sec VII, v2, verbatim).

"In this paper, we have developed and presented a fully relativistic numerical code based on the Misner-Sharp formalism for spherically symmetric cosmological scenarios containing both a scalar field and a perfect fluid. Our focus has been on cases where the scalar field dominates the dynamics, allowing us to explore the evolution of perturbations in scenarios such as slow-reheating after inflation.

We validated our numerical implementation by reproducing known results in the (perfect fluid) radiation-dominated scenario. Our results closely match the threshold for PBH formation in this regime, with a critical compaction function of approximately \mathcal{C}_{\rm th}=0.5, consistent with previous works. The ability of our code to reproduce these results confirms the reliability of our approach.

We explored the case of a quartic potential for the scalar field, which is known to mimic the behaviour of a radiation-dominated universe in perturbative scenarios. As expected, we found that the PBH formation threshold is, up to the current resolution of our code, equal to that of the radiation-dominated case. This confirms the consistency of the analogy between a quartic scalar field potential and a radiation-dominated universe, even in a non-linear, strong gravity regime.

Interestingly, we find that the quadratic potential case deviates significantly from dust-like behaviour in the non-linear regime. Our simulations indicate that overdensities experience an effective quantum pressure that counteracts gravitational collapse. However, due to limitations in computational power, we were unable to evolve fluctuations with sufficiently low amplitudes, which in turn prevented us from determining a precise threshold for PBH formation in this case, leaving us only with an upper bound. Nevertheless, the behaviour of the perturbations we evolved suggests that PBH formation may differ significantly from the pure dust case for which we have not managed to find the formation of an apparent horizon, with a reasonable possibility that stable solitonic/virialised structures could emerge as the final outcome of gravitational evolution for some initial conditions.

A key limitation preventing us from exploring smaller-amplitude perturbations is the difficulty in maintaining numerical control over the highly oscillatory dynamics of the scalar field for small inhomogeneities in this case. The exponential increase in the number of field oscillations, combined with the steepening of metric gradients, requires high spatial resolution and computational resources. Even after applying regridding techniques and using high-resolution grids, we find that the Hamiltonian constraint control expression eventually exceeds acceptable bounds (on the order of 10^{-2}), marking the limit of numerical reliability. This constraint violation reflects not a breakdown of the physical dynamics, but rather the onset of numerical instability. As a consequence, we restrict our analysis to a regime where the simulations remain stable and accurate, enabling a robust assessment of PBH formation.

In summary, our study provides new insights into PBH formation in scalar field-dominated scenarios using a fully relativistic numerical approach. While PBHs readily form in a quartic potential scenario, the quadratic case suggests the possible emergence of solitonic structures instead of PBHs. These findings highlight the importance of considering wave-like effects and quantum pressure in scalar field cosmologies.

In previous works, the results of analytic approximations to the scalar field collapse, in terms of associated threshold amplitudes for PBH formation, have been employed in constraining inflationary models with two characteristics: extended reheating phases, and features in the potential towards the end of inflation [69, 70, 72]. The implications for probes of the extended reheating scenario could be recast in light of the results of the present paper. On one hand, we are now confident that the constraints of the standard instantaneous reheating are readily applicable to the \phi^{4}, self-interacting case. On the other hand, the \phi^{2} case is still not conclusive and deriving constraints may require information from the number of e-folds and energy scale of the reheating period. Future work should explore these effects in greater detail, looking into more general scalar field models, as well as considering simulations with more than one matter component, phase transitions, and isocurvature modes sourcing the inhomogeneities. All these issues shall be explored elsewhere.

Furthermore, our results may have implications beyond the study of PBHs. Scalar fields are frequently considered as viable dark matter candidates, particularly in the context of ultra-light or self-interacting scalar field models. Although resolving the small-amplitude perturbations relevant for structure formation in these scenarios remains a major numerical challenge, our fully relativistic simulations offer valuable insights into the nonlinear dynamics of scalar fields under gravitational collapse. In particular, they may shed light on the possible formation of massive or even supermassive black holes in scalar field dark matter models, where solitonic cores and wave-like effects play a crucial role in the formation of these objects (see for example [85])."

v1's fourth paragraph instead reads: "...due to numerical limitations, we were unable to evolve fluctuations with an initially low amplitude ... suggests that PBH formation may differ significantly from the pure dust case, with a reasonable possibility that stable solitonic/virialised structures could emerge as the final outcome of gravitational collapse for some initial conditions." Also Sec VI.2 (v2) ends: "Future work is needed to determine with more accuracy the threshold value for the collapse of inflaton stars."

Paper 2 — de Jong, Aurrekoetxea, Lim, arXiv:2109.04896, JCAP03(2022)029 #

(a) Footnote 1 (Sec I, attached to the sentence on ξ diluting):

"In principle, we could use a single massive scalar \phi. 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."

(b) Setup.

• Sec I: "Since the field \xi has no potential, it will only influence dynamics via its gradients. Furthermore, it will dilute much more rapidly than dark matter, and hence not affecting the long term dynamics of the system once its initial job of sourcing a perturbation is done" ... "Meanwhile, the massless scalar field \xi provides the energy density perturbation that will trigger BH formation. In this paper, we exclusively consider initially static spherically symmetric perturbations ... \xi(t=0,r)=\Delta\xi~\tanh{\Big[\frac{r-R_{0}}{\sigma}\Big]}, (9) where \Delta\xi, \sigma and R_{0} are the amplitude, width and the initial size of the perturbation respectively. The mass of the initial perturbation scales roughly as R_{0}^{2}. We emphasise that this perturbation is non-linear, despite its moniker." ... "The background scalar field \phi starts from rest, so that \dot{\phi}=0".
• Appendix A.2, Eqs. (25)–(26): "V_{\phi}(\phi)=\frac{1}{2}m^{2}\phi^{2},\qquad V_{\xi}(\xi)=0 ... \phi(t=0,x^{i})=\phi_{0}, \xi(t=0,x^{i})=\Delta\xi\tanh[(r-R_{0})/\sigma_{0}], \partial_t\phi=\partial_t\xi=0."

(c) Accretion when hoop fails.

• Abstract: "In particular, for the latter case, the initial perturbation does not have to satisfy the hoop conjecture for a black hole to form."
• Sec II.1, Eq. (18): "applying the hoop conjecture [fn 6] suggests that if the condition \sigma(a_{*})<2GM_{\mathrm{infall}} is satisfied ... then a black hole will form."
• Sec II.2: "On the other hand, if \sigma(a_{*})>2GM_{\mathrm{infall}}, a black hole does not form directly. In this case, the energy density of the perturbation \rho_{\xi} disperses after reaching the centre and becomes locally sub-dominant to the background energy density \rho_{\mathrm{DM}}. Nevertheless, the presence of \xi generates a gravitational potential well in the centre, which seeds accretion of the background DM and eventually causes a collapse into a black hole."

(d) Parameters, code, solver, cost.

• Sec II: "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 \phi to be \phi_{0}=7.8\times 10^{-3}M_{\mathrm{Pl}}, with m\approx 10^{2}H_{0}. ... R_{0}H_{0}\in[0.575,~1.6] ... \Delta\xi M_{\mathrm{Pl}}^{-1}\in[0.075,~0.12], whilst keeping the initial width fixed to \sigma_{0}=0.15H_{0}^{-1}, such that the ratio between the maximum gradient energy density to dark matter energy density is \rho_{\xi}/\rho_{\mathrm{DM}}\sim 1."
• Appendix A.2: "In our simulations, mH_{0}^{-1}=62.6." (Note the main-text "≈10²H_0" vs appendix 62.6 discrepancy.)
• Appendix A.1: "This work was written based on simulations run using \mathtt{GRChombo} [133], with the CCZ4 formulation of the Einstein equations [134]."
• Sec I: "We choose a conformally flat ansatz for the 3-metric ... We choose an initially expanding spacetime with K=-3H_{0}, so that the periodic integrability condition is satisfied [fn 3: ... K^{2}=24\pi V(\phi_{0})] ... 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 \psi numerically [Eq. (14), radial ODE]". Appendix A.2: "Momentum constraints are trivially satisfied, and we solve the Hamiltonian constraint (14) numerically to find the conformal factor \psi."
• Appendix A.2: "we reduce the computational cost of evolution by simulating one eighth of the system using symmetric boundary conditions." Appendix B: "three different base grid resolutions, namely N_{\mathrm{LR}}=80, N_{\mathrm{MR}}=96, N_{\mathrm{HR}}=128 ... showing convergence to 1%."
• Sec III: "we were unable to track the growth of PBH beyond a few factors of their initial mass, as the numerical cost becomes prohibitive." Sec III: "M_{\mathrm{BH}}H_{0}\sim 10^{-2}M_{\mathrm{Pl}}^{2} – see Fig. 6."
• Box size L: NOT FOUND in this paper. Core-hours: NOT FOUND; Acknowledgements list only facilities (SuperMUC-NG PRACE 2018194669, JUWELS PRACE 2020225359, COSMA7, DiAL DiRAC ACTP238, CSD3).

Paper 3 — de Jong, Aurrekoetxea, Lim, França, arXiv:2306.11810 (v2 7 Jul 2023; JCAP 10 (2023) 067 per GRTresna ref [49]) #

(a) Sec II.1: "For all simulations we choose \phi_{0}=7.8\times 10^{-3}m_{\mathrm{Pl}} and the mass m\approx 62H_{0}. We find that these values accurately model a matter-dominated universe expansion on average." Footnote 1: "Planck units \hbar=c=1, such that G=m_{\mathrm{Pl}}^{-2}, where m_{\mathrm{Pl}} is the non-reduced Planck mass." Eq. (5): "H_{0}^{2}=\frac{4\pi m^{2}}{3m_{\mathrm{Pl}}^{2}}\phi_{0}^{2}." Sec I: "If \phi oscillates coherently with a period considerably smaller than a Hubble time, 2\pi/m\ll 1/H, its pressure averages to zero".

(b) Sec I.1, Eqs. (6)–(7): "\xi(t,r,\theta,\alpha)=R(r)\left[1+\sin{(\theta)}\Phi(t,\varphi)\right], (6a) \Pi_{\xi}=R(r)\sin{(\theta)}\partial_{t}\Phi(t,\varphi), (6b) where R(r)=A\exp{\left[-\frac{(r-R_{0})^{2}}{\lambda^{2}}\right]}, (7a) \Phi(t,\varphi)=B\cos{(k\varphi-\omega t)}. (7b) The constants A, R_{0}, \lambda, B, k, \omega represent the perturbation shell's initial radial amplitude, radius, width, spin amplitude, spin wavenumber and spin angular velocity, respectively". Sec II.1: "For the initial size of the superhorizon perturbation we use a range of R_{0}\in[1.2,1.5]H_{0}^{-1} with initial width \lambda=0.15H_{0}^{-1}. The spin amplitude and wavenumber are fixed to B=0.5 and k=2, respectively." ... "Consequently, for the perturbation's radial amplitude we use a range of A\in[0.0825,0.09]m_{\mathrm{Pl}}. We vary \omega\in[0,12]H_{0}^{-1} to parameterize the amount of angular momentum in the system."

(c) Sec I.1: "As is well known, all initial data in general relativity must obey a coupled system of Hamiltonian and momentum constraint equations, which we solve using the CTTK method [129]. We choose an initial spatially flat metric \gamma_{ij}=\delta_{ij} and solve for the trace K and traceless parts A_{ij} of the extrinsic curvature tensor. This choice is equivalent to choosing an initially homogeneous cosmological scale factor, where the rate of local expansion is determined by the matter distribution."

(d) NOT FOUND. All six footnotes checked (units; φ vs ϕ notation; x/y angular momentum negligible; "The angular momentum transfer between the fields \xi and \phi is minimal if they only couple via gravity."; accretion-rate caveat; AH area). The only characterisation is Sec I: "We simulate the collapse of superhorizon non-linear perturbations sourced by a massless scalar field, on a matter-dominated expanding background driven by an oscillating massive scalar field."

(e) Appendix C: "three different base grid resolutions, N_{1}=80, N_{2}=96 and N_{3}=128, which correspond to base grid spacings h of h_{1}=6.25\times 10^{-2}H_{0}^{-1}, h_{2}=5.21\times 10^{-2}H_{0}^{-1} and h_{3}=3.90\times 10^{-2}H_{0}^{-1}. We track Hamiltonian and momentum constraint violation, as well." Sec III: "Computational cost limits us from tracking the PBH evolution beyond a few e-folds after formation". Core-hours, AMR levels, box size: NOT FOUND in this paper (facilities: DiRAC@Durham ACTP238/ACTP316, DiRAC Leicester).

Paper 4 — de Jong PhD thesis, arXiv:2403.02878 (166 pp., supervisor E. A. Lim) #

(a) Abstract: "Independent of the formation mechanism, the PBH forms within an efold after collapse is initiated and with a small initial mass compared to the Hubble horizon, M_{\mathrm{BH}}H_{0}\sim 10^{-2}m_{\mathrm{Pl}}^{2}. Finally, we find that PBH formation is followed by extremely rapid growth M_{\textrm{BH}}\propto H^{-\beta} with \beta\gg 1, during which the PBH acquires most of its mass." Sec 3.5: "In the cases of both direct and accretion collapse, the initial mass of the PBH formed is small compared to the Hubble horizon, M_{\mathrm{BH}}H_{0}\sim 10^{-2}m_{\mathrm{Pl}}^{2} – see Fig. 3.9(a). Once the initial PBH has formed, the PBH accretes DM from its surroundings in the growth phase ... 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." Sec 3.6: "In both the direct collapse and accretion collapse formation cases, the initial mass of the PBH is roughly M_{\mathrm{BH}}\sim 10^{-2}H^{-1}m_{\mathrm{Pl}}^{2}, but formation is followed by an extremely rapid growth M\propto H^{-\beta} where \beta\gg 1."

(b) Ch. 5 intro: "In this chapter, we describe an extension of the modified cartoon method from Cook:2016soy, which involves the addition of matter fields. We discuss the dimensional reduction of the initial condition solver used in chapter 4 and the evolution equations used in chapters 3 and 4. In this chapter we limit ourselves to the case of a real scalar field". Sections: "5.2 Cartoon expressions", "5.3 Initial conditions", "5.4 Cartoon expressions for a scalar field in twist-free axisymmetry", "5.5 Vacuum bubble collision with dimensional reduction". Sec 6.1.2: "We extend the cartoon formalism presented in Cook:2016soy by adding matter, specifically a real scalar field, and we give cartoon expressions for dimensional reduction of the CCTK method for solving the Einstein constraints. ... It should be straightforward to import these expressions into working codes ... Furthermore, our own cartoon codes for both the initial constraints and the evolution will become public in the near future". Sec 1.1: "we run our simulations using the code \mathtt{GRChombo} ... which is open-source and publicly available". Sec 3.3.2: "the spherically symmetric scenarios described in this chapter can be studied in a computationally more efficient manner using a dimensionally reduced 1+1D code. However, we choose a 3+1 setup to make subsequent generalisation to scenarios with f[ewer symmetries]...".

Directly relevant to the brief (Sec 6.2 Future work): "Moreover, the dynamics of the collapse may change if the only field present is the background field, which is perturbed itself to provide the overdensities collapsing to PBHs. We find that this can have large effects on the local expansion rate of the universe and how this influences the collapse and PBH formation requires further study."

(c) Core-hours/node/wall-clock: NOT FOUND. Sec 3.3.3: "The length of the simulation box is L_{\textrm{box}}=2/H_{0} (L_{\textrm{box}}=2.5/H_{0}) for initially subhorizon (superhorizon) perturbations. ... The grid size of the coarsest level ... is 48^{3} and the refinement factor ... is 2. The Courant factor \Delta t/\Delta x=0.0384 ... The maximum number of levels achieved for these simulations is 4, including the base level." Sec 3.4: "with m=62.6H_{0}. For larger values of \phi_{0}, the \phi-oscillation amplitude decays rapidly during the first Hubble time ... value for m makes sure that \phi performs around 10 oscillations during the first Hubble time, which we find is sufficiently rapid to cause a matter-like universe expansion. At the same time, it is not so fast that we need to reduce the size of the Courant factor to a value so small that it would make the simulation unnecessarily computationally heavy." Sec 4.3.2: "We use periodic boundary conditions in the x- and y-directions and reflective boundary conditions in the z-directions ... L_{\textrm{box}}=4/H_{0} (L_{\textrm{box}}=5/H_{0}) ... coarsest level ... is 96 in the x- and y directions, and half in the z-direction ... maximum number of levels ... is 4". Sec 5.5.2: "Had the grid been three-dimensional, the computational cost would increase by roughly a factor 250 for the lightest simulation". Appendix: "We have not run new convergence tests with the corrected method, because these runs are computationally expensive."

Paper 5 — GRTresna, arXiv:2501.13046 (v2 26 May 2026; JOSS, DOI 10.21105/joss.09057) #

(a) Header: "GRTresna is a multigrid solver designed to solve the constraint equations for the initial data required in numerical relativity simulations. In particular, it is focussed on scenarios with fundamental fields around black holes and inhomogeneous cosmological spacetimes. The code is based on the formalism in Aurrekoetxea, Clough & Lim [1] and can be found at https://github.com/GRTLCollaboration/GRTresna". Sec 1: "GRTresna implements two variations of the CTT method recently introduced in Aurrekoetxea, Clough & Lim [1]: the CTTK and CTTK-Hybrid methods, which are particularly well-suited to cases with fundamental fields in the matter content." Sec 2: "Flexibility: ... It currently supports cosmological-type periodic spacetimes and a superposition of two boosted and/or spinning black holes (Bowen-York initial data), with fully general scalar field matter source configurations and the flexibility to adapt to other setups. While scalar fields are the only matter sources included in the current version of the code, the templated methods allow users to easily replace them with other matter types". "Methods: GRTresna incorporates the CTTK and CTTK-Hybrid methods to solve the Hamiltonian and momentum constraints." "Boundary conditions: The code implements extrapolating, reflective, and periodic boundary conditions, compatible with those in the NR evolution code GRChombo". Repo README: "GRTresna is licensed under the BSD 3-Clause License." (v1, Jan 2025, has the same sentences.)

(b) How K is solved: NOT described in the GRTresna paper. CTTK paper (arXiv:2207.03125, CQG 40 (2023) 075003) abstract: "instead of solving the Hamiltonian constraint as a 2nd order elliptic equation for a choice of mean curvature K, we solve an algebraic equation for K for a choice of conformal factor ... we show that this method provides rapid convergent solutions for several initial conditions ... namely (i) periodic inhomogeneous spacetimes with large random Gaussian scalar field perturbations and (ii) asymptotically flat black hole spacetimes with rotating scalar clouds." Examples in the GRTresna paper (Sec 4): "The robustness of inflation to inhomogeneities in the scalar field [47, 48]. The formation of oscillons during inflationary preheating [12]. Formation of spinning primordial black holes [49]. The effect of scalar dark matter environments around binary black holes [9, 10, 11]. The general relativistic evolution of polarized Proca stars [50]. Solving the initial conditions problem for modified gravity theories [13]." Forthcoming: "non-conformally flat metric data, new matter types including vector fields ... and dimensional reduction to 2D using the cartoon formalism". A growing-mode cosmological perturbation of the dominant massive scalar: NOT FOUND in the paper. Repo Examples/ScalarFieldCosmo/params.txt (not the paper) contains: phi_0 = 1e-1, dphi = 1e-1, pi_0 = 1e-1, dpi = 1e-1, scalar_mass = 1e-1, is_periodic = 1 1 1, sign_of_K = -1, deactivate_zero_mode = 1 ("the Garfinkle trick"); no growing-mode construction is described there.

Key takeaways for the brief #

1. Paper 1 is not a 3+1 EKG/BSSN code; it is a Misner–Sharp scalar+fluid code, unnamed and (as far as the text says) not public, and it never states the value of m/H — only that oscillations are "rapid" (2μ ≫ 3H) and that their number grows exponentially as κ decreases.
2. The claim "every amplitude in κ∈[0.06,0.2] formed a horizon" and the bounds C_th ≤ 0.2, δ_th ≤ 0.077 are supported; the limiting factor is the Hamiltonian-consistency test exceeding 10⁻² despite N up to 10,000.
3. The only explicit statements about the difficulty of perturbing the massive field directly are Paper 2 footnote 1 and thesis Sec 6.2 (large effect on the local expansion rate); Paper 3 has none.
4. Paper 2 states m/H_0 two ways (≈10² in text, 62.6 in appendix); the thesis says 62.6 and "around 10 oscillations during the first Hubble time".
5. No core-hour figures exist in any of the GRChombo papers or the thesis; grid/AMR parameters are only in the thesis.
6. No 2025–26 paper from either group on a directly perturbed oscillating scalar in 3D was found; nearest are 2509.10431 (spherical, quartic critical collapse) and 2509.26470 (existence of constraint solutions for cosmological overdensities).