We leave a proper de-convolution from the top-hat window function $W(r,R)=\mathrm{\Theta }(R-r)/V$ followed by inverse-Fourier transformation for future work.

In each patch, we randomly sample the nucleation times $t_j$ and positions $d_j$ of the first $j_c=50$ bubbles and assume a deterministic evolution for subsequent nucleation dynamics [20]. The probability distributions for $t_j$ and $d_j$ within a patch of radius $R$ are derived analytically in Appendix A of Ref. [20] and are summarized in the SuM [52].

As known from textbooks, e.g., Ref. [65], and rederived in SM [52], the approximation ${\mathrm{\Phi }}^{\prime }\simeq 0$ holds for superhorizon modes only in the absence of nonadiabatic pressure perturbations, ${\delta }_{p,\mathrm{nad}}=0$. In our scenario, bubble dynamics generate ${\delta }_{p,\mathrm{nad}}\ne 0$, leading to departures from ${\mathrm{\Phi }}^{\prime }\simeq 0$ on superhorizon scales. Nevertheless, ${\mathrm{\Phi }}^{\prime }\simeq 0$ remains a good approximation near horizon crossing, $k\simeq \mathcal{H}$, justifying its use in evaluating the PBH collapse threshold. The time evolution of $\mathrm{\Phi }(t)$ is shown in SM [52].

The relation $\mathcal{R}\simeq -9{\delta }^{(C)}/4$ follows from Eqs. (16) and (17) for $k=\mathcal{H}$, $\omega =1/3$, and ${\mathrm{\Phi }}^{\prime }\simeq 0$.