It is important to stress the nonperturbative nature of the non-Gaußianity in this work. In cosmology, the commonly used parameters like $f_{\mathrm{NL}}$ or $g_{\mathrm{NL}}$ describe an expansion around Gaußianity. While it might well be that for a given rarity, such as $6\sigma$, the mentioned expansion might actually be sufficient for the derivation of the studied effects, the exact exponential tail, however, cannot be obtained from any finite expansion based on Gaußianity.

We appreciate that most of the literature discusses continuous situations in which a peak is an extended feature (unless a delta-peak is assumed). On the contrary, in our discrete study we refer to a peak as an excess in a single pixel (representing an initial Hubble patch at the time of horizon crossing of an overdensity).

More advanced estimates utilize the so-called compaction function (see, e.g., Ref. [43])—an approach which unavoidably induces some degree of smoothing, hence—given the discrete (pixelated) approach used in this work—it is a relevant question to ask whether and how the use of the smoothing due to the use of compaction functions alters the results. This could certainly have large implications on subhorizon scales, but in our work, each pixel is thought to represent an entire initial Hubble patch, so we do not expect a strong effect.

Observe that the broadening of the PBH mass function does not only extend toward higher masses, due to the absorption of clustered peaks, but also toward lower masses, since neighboring (partly overlapping) peaks generically lower the PBH formation threshold (see Ref. [40]).

Note that this effect cannot directly compare to the Gaußian random field case, since Gaußian random fields do not possess a heavy tail. Specifically, if the threshold for PBH formation is set to 20 standard deviation ($6\sigma$ rarity) of the exponential random field, then a Gaußian random field of the same variance would yield essentially no threshold exceedance.