Primordial black holes as dark matter: Interferometric tests of phase transition origin

We show that primordial black holes - in the observationally allowed mass window with $f_{\rm pbh}=1$ - formed from late nucleating patches in a first order phase transition imply upcoming gravitational wave interferometers will see a large stochastic background arising from the bubble collisions. As an example, we use a classically scale invariant $B-L$ model, in which the right handed neutrinos explain the neutrino masses and leptogenesis, and the dark matter consists of primordial black holes. The conclusion regarding the gravitational waves is, however, expected to hold model independently for black holes coming from such late nucleating patches.


Introduction
Putting aside the possibility being observationally misled by a theory of modified Newtonian dynamics [1][2][3][4][5][6], the microscopic nature of the cosmic fluid labelled dark matter (DM) is an open question.One possibility is for all of the DM to be made up of primordial black holes (PBHs) [7][8][9].In this paper, we adopt the standard notation of denoting the PBH-to-DM density ratio, and assume throughout Observationally, this is allowed for approximately monochromatic PBH distributions with mass in the range [10], 10 −16 M ⊙ ≲ M PBH ≲ 10 −10 M ⊙ . (2) Here the lower limit comes from the would-be-flux of particles coming from Hawking evaporation -which would have altered the CMB [11,12], produced too much extragalactic background light [13][14][15], been detected by Voyager [16] or SPI/INTEGRAL [17][18][19] -and the upper limit comes from microlensing observations [20].
In order for PBHs to be produced in the early universe, physics beyond the standard model (SM) of particles and cosmology is required.One possibility is for a feature in the inflaton potential to give an enhancement in the amplitude of overdensities at small scales [21,22].When the overdensities re-enter the Hubble horizon, these will collapse into PBHs, provided the density contrast is sufficiently large for the given equation-of-state and shape of the curvature power spectrum [23][24][25][26][27][28][29][30][31][32][33].As the evidence for inflation is rather compelling -albeit circumstantialthis possibility has garnered a good deal of attention.In such scenarios, gravitational waves (GWs) are produced at second order in perturbation theory, when the enhanced curvature perturbations at small scales re-enter the Hubble horizon.Such GWs may be detected at future interferometers [34][35][36][37][38][39].
An alternative is for the overdensities which collapse into PBHs to be produced at a later stage in the evolution of the universe.This could occur if there was a sufficiently strong firstorder phase transition (PT) between the end of slow roll inflation and big bang nucleosynthesis.This possibility has also been the topic of a number of papers which have appeared in the last four decades, albeit more sporadically than those regarding the inflation scenario, and a number of mechanisms allowing for PBH formation during a PT have been identified [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57].
Among the PT mechanisms that trigger PBH formation, the formation of large false vacuum remnants from stochastically late-nucleated patches in supercooled phase transitions stands as an intriguing possibility [42,43,47,48,[52][53][54][55][56][57].This is an unavoidable process if sufficiently strong PTs take place.The remnants lead to Hubble sized overdensities being present in the radiation following post-PT reheating, most closely emulating the inflation mechanism.The concurrent production of smaller PBHs has also recently been considered in [55], using a Hoop conjecture criterion, but the spectrum remains relatively peaked. 1he reheating following bubble percolation in the PT leads to an initially inhomogenous state, not modelled in General Relativistic numerical simulations of PBH formation [24][25][26][27][28][29][30][31][32], but the large number of bubbles present and the fact that an O(1) overdensity is achieved, leads us to assume that collapse is rather plausible.Therefore, for the scope of this paper, we operate under the assumption that this mechanism functions as expected, along the lines of [47, 48, 52-54, 56, 57].Nevertheless, we warn the reader that the collapse criterion in this scenario, which will necessarily feature a certain level of anisotropy, turbulent inhomogeneity, non-sphericity, and initial velocity vectors for the radiation fluid differing from the inflationary overdensity scenario, is far from certain.
Strong first order PTs also produce a stochastic background of GWs [64,65].The production of the PBHs in the PT will lead to GWs from the bubble collisions, and these GWs are observable in upcoming interferometers, as we explain next using a simple argument.In the rest of the subsequent paper we provide a more detailed calculation.

A simple argument
In the delayed patch mechanism, the mass of the PBH is set by the energy contained inside a Hubble volume during the PT, where M Pl ∼ 10 19 GeV is the Planck mass, and T RH is the reheating temperature following the completion of the PT (we assume rapid reheating).The observationally allowed mass window for which f PBH = 1 is possible, Eq. ( 2), constrains the reheating temperature to lie between [56] 10 TeV ≲ T RH ≲ 10 4 TeV.
The abundance is controlled by the time the PT takes to complete compared to the Hubble rate.To achieve an abundance of f PBH = 1, we require the inverse timescale of the transition to be [56], where Γ bub is the bubble nucleation rate per unit volume and H is the Hubble rate during the PT.
The amplitude of GWs from bubble collisions during the PT scales as ∝ R 2 ∼ v 2 w /β 2 , where R is the bubble radius and v w is the wall velocity.For the PTs of interest here, it is a good approximation to set v w ≃ 1, and we do so throughout the rest of this paper.The precise amplitude and spectral shape of the GWs has been addressed in a number of works [66][67][68][69][70][71][72][73][74][75][76][77].In our calculations we shall use the recent determinations, relevant for our type of PTs, in which the peak amplitude redshifted to today is [73,75,77] Said GW estimates still feature uncertainties for the very strong PTs we will be interested in, due to the expansion of the universe during the PT itself, see [78,79].So in this regard our results should be considered preliminary, albeit promising, until the GW predictions can be further refined.The peak frequency of the GW signal is similarly set by the inverse bubble radius at the collision time, which is then redshifted to today, From the above, we see that in such a f PBH = 1 from a PT scenario, we will have a peak amplitude Ω GW ∼ 10 −8 , i.e. well above stochastic astrophysical foregrounds, somewhere in the frequency range 10 −3 Hz ≲ f peak ≲ 1 Hz.The GWs at lower peak frequencies are detectable by LISA [80].Those at higher peak frequencies are similarly detectable by the Einstein Telescope (ET) [81].A mid-frequency interferometer such as BDECIGO [82] would be sensitive to the entire allowed PBH mass range and could in any case help confirm and characterise the signal.Similar conclusions hold if we replace a detector with a variant of broadly similar equivalence, such as ET with Cosmic Explorer [83] or BDECIGO with AEDGE [84].Note we consider only the GWs from the bubble wall collisions, the calculation and addition of GWs sourced from the enhanced curvature perturbations at small scales in the context of the PT mechanism [85], is left for future work.
The rest of the paper is dedicated to calculating the PBH production and GW signal within a simple beyond the SM scenario, motivated by neutrino masses and leptogenesis, for which the PBHs will act as a DM candidate.We will scan over the parameters giving M PBH in the allowed window.The model we choose is a classically scale invariant, gauged B − L, extension of the SM.Such close-to-conformal models can provide the necessary supercooled PTs [86][87][88].Even though fine-tuning is required to get the cosmological constant close-to-zero today, this seems to be the case almost generically when dealing with the vacuum energy density in a cosmological context.We therefore put questions of fine-tuning aside, which may anyway similarly be present in other models triggering PBH formation, and simply explore the scenario of achieving the observed ρ DM .
We find it beneficial to conduct an explicit calculation within the aforementioned beyond the SM physics scenario.The calculation in such a model serves to illuminate several facets of the subject matter.For instance, many previous studies (e.g.[47,56]) assume a model independent approximation for the bubble nucleation rate where t n is the bubble nucleation time, by definition when Γ bub = H 4 , and β acts as the coefficient of the first order Taylor approximation of the full nucleation rate.One may wonder, however, whether the second and possibly higher derivatives will effectively play a role in setting f PBH , because the probability of a Hubble patch collapse is very sensitive to the precise behaviour of Γ bub .Nevertheless, our results below show the approximate form yields accurate predictions for the required β, at least in close-to-conformal potentials.Thus confirming the applicability of the GW signal based on the requirements on β for achieving f PBH = 1, as outlined in our discussion around Eqs. ( 3)- (7).
Let us now also mention some work which previously headed along this path, but with different overall scope, e.g.[47,48,53,89,90]. 2 Indeed, already in 1984, Witten was discussing the concomitant aspects of a strong QCD phase transition, GWs, production of PBHs, and future pulsar timing array measurements [64].Compared to the aforementioned papers, we provide a thorough scan over the expected PBH mass range coming from the PT, giving a more detailed picture of the expected GW signal.We also further develop the late nucleating patch formalism [47,56], showing the calculation can be recast from using time as a variable to a conveniently chosen temperature, here that of the false vacuum plasma.
Our focus here is on f PBH = 1 which results in a scale of new physics beyond collider reach.The interplay of micro-lensing hints for 0 < f PBH < 1 at larger PBH masses, GW signals from supercooled PTs, and collider physics has instead recently been studied in [91].
3 The Model

Particles and charges
As our example, we take a classically conformal B − L extension of the SM, in which we add three right handed neutrinos, N i with Q B−L = −1, to cancel off anomalies [86,92,93].To break the B − L symmetry, we add a complex scalar, ρ, with Q B−L = −2.Ultimately, the radiative breaking of the B − L symmetry induces a negative mass term for the electroweak Higgs field via a portal coupling λ ρh , thus eventually breaking electroweak symmetry.The kinetic term for the new scalar reads where D µ ≡ ∂ µ − 2ig B−L Z ′ µ is the covariant derivative with gauge coupling g B−L .This will generate a mass for the gauge boson once ρ gains a vacuum expectation value (vev) radiatively, ⟨ρ⟩ = v ρ / √ 2. At the classical level, the scalar potential reads where H is the SM Higgs doublet.In addition to the SM Yukawas we also have where l Li are the SM lepton doublets and H ≡ iσ 2 H * .The second term in Eq. ( 11) will give Majorana masses to the N i .These can then decay in a CP violating manner via the first Yukawa coupling realising leptogenesis, on which we will briefly comment later.The neutrino oscillation data [94][95][96][97][98][99][100][101] is explained through a type-I seesaw [102][103][104][105] (here at relatively low scales ≲ 10 8 GeV).In the regime we are interested in, there is a significant hierarchy between v ρ and the electroweak vev v ew ≃ 246 GeV, which implies a small cross quartic coupling, where m h ≃ 125 GeV is the SM-like Higgs mass.This hierarchy situates the scenario beyond the reach of current collider constraints.

Radiative symmetry breaking at finite temperature
Our treatment of the field theory in this section, leading to the approximation of the effective potential, V eff , follows [87].At one loop the beta function of the B − L Higgs quartic is where the sum over the y N i is understood (i = 1, 2, 3).Consistency requires β λρ > 0, so that λ ρ turns negative when running down from high scales, triggering radiative symmetry breaking.
Making the here justifiable approximation of the symmetry breaking as dominantly in a single field direction, following the approach of Gildener and Weinberg [106], we write the one loop zero temperature potential as where from now we denote v ρ as the vev at the zero-temperature minimum and ρ denotes the classical field value.At the B − L symmetry breaking scale, λ ρ ≃ 0, because the coupling is switching signs.We also require a small λ ρh , as mentioned above.Therefore we can approximate the beta function as The mass of the physical scalar after symmetry breaking is m ρ ≃ √ β λ v ρ .Its field dependent mass gives only a small contribution to the thermal corrections in the effective potential.The dominant thermal corrections instead come from the field dependent mass of the B − L gauge boson, and the heavy Majorana neutrinos, Their thermal contributions are where the thermal functions are for bosons and fermions respectively.To improve the perturbative analysis, we include the corrections due to the resummation of daisy diagrams [107], where the thermal mass of the longitudinal component of the B − L gauge boson is The full effective potential we consider is therefore after the aforementioned approximations have been made.To be consistent with observations, showing a small vacuum energy in the present day, a cosmological constant term, needs to be added to the potential, so that V 0 (v ρ ) = 0.The gravitationally important terms, Λ vac , together with the Planck mass, M Pl , must be considered external to our classically scale invariant sector.An example showing the temperature evolution of the potential is shown in Fig. 1.Nucleation is possible once the false and true minima become degenerate at the critical temperature (time), T c (t c ).

Bubble nucleation and percolation
The bubble nucleation rate per unit volume can be approximated as [108][109][110] where T is the temperature, A is a mass dimension one pre-factor we elucidate later, and S is the Euclidean bounce action.At zero temperature the configuration minimizing the action, S ≡ S 4 , is O(4) symmetric and In this work we can limit ourselves PT dominated by single field directions.The action S 4 corresponds to quantum tunneling through the potential barrier.The equation of motion is with boundary conditions dρ dr r=0 = 0, and lim At finite temperature the field instead becomes periodic in 1/T in the time coordinate.The configuration minimizing the action is O(3) symmetric.Furthermore, at sufficiently high temperatures, the minimum action configuration becomes constant in the time direction and This represents bubble formation through classical field excitation over the barrier.The corresponding equation of motion is then with the same boundary conditions as above.The solution with a non-trivial periodic bounce in the time coordinate, corresponding to quantum tunneling at finite temperature, is more difficult to evaluate in practical calculations.It is eventually, however, well approximated at low T by S 4 but calculated with V eff including any finite temperature corrections.For simplicity, we therefore take a rather standard estimate for our nucleation rate where R c ∼ 1/T is the bubble radius in the low T limit of radiative symmetry breaking, and we have used standard estimates of the pre-factor (for discussion regarding uncertainties see [111]).
For the PTs we study, we have verified that Γ bub is given by the S 3 action.The nucleation temperature in an average Hubble patch, T n , or the equivalent nucleation time, t n , is defined as when the nucleation rate equals the Hubble rate, In practice, we compute the nucleation temperature using the Hubble rate in the false vacuum, where T is the false vacuum temperature, and g * are the radiation degrees-of-freedom.In our scenario, g * ≃ 116 prior to the PT, and g * ≃ 113 just after the PT.Henceforth, any instances of the Hubble parameter and the scale factor presented without an index should be understood as these quantities evaluated in the false vacuum phase.We will be interested in strong supercooled phase transitions for which the universe becomes vacuum dominated.This occurs at a temperature For the case we will be interested in, namely efficient reheating and negligible changes in degreesof-freedom between the two phases, we have where T RH is the reheating temperature following the completion of the phase transition, and T p is the bubble percolation temperature, which we now review how to calculate.
In the vacuum dominated regime, the temperature (time) at which the bubbles percolate, T p (t p ), may be appreciably different from T n (t n ).It is useful to introduce the comoving radius of a bubble at false vacuum temperature T , nucleated at some higher temperature T ′ , where a is the scale factor, we have assumed v w ≃ 1 shortly after nucleation, and have used dT = −T Hdt.The physical radius is then where we have assumed a constant g * and so used a(T )T = a( T ) T .The expected volume of true-vacuum bubbles per comoving volume (double counting overlapping regions and including fictitious nucleations in true vacuum) is given by [112][113][114][115][116] To find the probability of a given point in the comoving volume to be in the false vacuum, we need to exclude fictitious nucleations and avoid double counting overlapped regions, and the appropriate expression is given by [112][113][114][115][116] The change in the physical false vacuum volume, V false = a 3 P (T ), normalized to the Hubble rate, is captured by the equation [116] 1 We define the percolation temperature T p as the highest temperature for which both I(T ) > 1 and d log V false /dt < −H hold (sometimes slightly less stringent conditions are employed, e.g.see discussion in [117], the precise choice does not affect our qualitative results below).Numerically, we find the first condition, I(T ) > 1, occurs later and hence determines the percolation condition in our parameter space.
To characterise the strength of the phase transition, we find the free energy difference between the true and false vacuua normalized to the radiation density, where T * is set to either the nucleation or percolation temperature. 3The inverse timescale of the phase transition normalized to the Hubble rate reads By numerically calculating the bubble action S for a sufficiently large number of points over a suitable temperature range, we obtain a good approximation of Γ bub (T ), which allows us to find all the above quantities.Phase transition parameters as a function of g B−L are shown in Fig. 2. Clearly in the limit of a small β λρ we have strong supercooling, with β H ≲ 8 and α ≫ 1, which allows for cosmologically significant PBH production, e.g. as found in [56].We now turn to calculating the PBH abundance.
4 Primordial black hole production

Background evolution
We first consider an average Hubble patch.The total energy density, ρ, consist of the vacuum energy density, ρ vac , and the radiation, ρ rad .The latter is made up of the plasma of relativistic particles together with the bubble walls moving at v w ≃ 1 [56].The Friedmann equations are given by where H bkg denotes the Hubble rate of the average background patch.We can relate ρ vac to the bubble nucleation rate through To solve the Friedmann equations, we find it useful to switch coordinates, from time to false vacuum temperature.We then have where T is the temperature in the false vacuum, and which we solve to find ρ rad (T ).Note ρ vac , through the quantity I which contains the bubble radius, effectively also contains the scale factor.Hence the equations should be solved self-consistently by evaluating I by taking into account the scale factor at each time step, which in turn can be found through H bkg .In order to solve the equations, however, a simplification is possible.First note that before percolation, H bkg is for practical purposes given by its corresponding false vacuum value.After percolation, the ρ vac source term for the radiation is, by definition, approximately negligible.We therefore first evaluate I numerically using the scale factors and Hubble rate in the false vacuum, H, and use this to solve the Friedmann equations including H bkg of the background patch in the first equation.We then check whether the approximation is a good one by using the resulting scale factor in an updated determination of I.The resulting shift in T p is tiny -we will quantify it below -and so we confirm the approximation is justified.
For our calculations, we will also need the scale factor as a function of the false vacuum temperature.Ignoring changes in degrees-of-freedom, for the scale factor in a false vacuum patch we have the standard relation, For the scale factor in the background patch that has nucleated, we instead solve The equation is solved starting from a sufficiently high temperature where H bkg = H and a bkg (T ) = a(T ).Note, for bookkeeping purposes, the temperatures which appear here should be understood as the temperature in the false vacuum patch, not the temperature in the nucleating background patch, which contains nucleated bubbles and partially reheated plasma.Using our solution, it is convenient to define which changes from unity to take into account deviations from the false vacuum relation, Eq. ( 44), for the background scale factor.In terms of the false vacuum temperatures and false vacuum H, the bubble radius in the background patch is then more precisely given by It can also be of interest to consider the bubble density as a function of radius [116,117], which rewritten as a function of the false vacuum temperature is given by where here T ′ ≡ T ′ (R) is understood as being the false vacuum temperature at which a bubble of physical size R bkg was nucleated, found by numerically inverting Eq. ( 47), from which one also finds the corresponding derivative.

Late patch evolution
The above describes the evolution of an average background nucleating patch.We are now interested in the evolution of a late nucleating patch.We assume no bubble nucleates in the patch until some temperature T i < T n .By modifying the high temperature terminal in Eq. ( 36), we find the vacuum energy in the late patch, Again we use the same approximation as before, in assuming H = H false , when calculating I late .The probability of finding a point in the false vacuum in the late patch is then Thus the vacuum energy density in the late patch is simply We then also solve the Friedmann equations for the late patch Eventually we will also need the scale factor in the late patch, found by solving Similarly to what we did for the background patch, it is convenient to define to parametrize deviations from the false vacuum scale factor/temperature relation.The bubble size in the late patch is then given by The bubble spectrum of the late patch is given by where the Θ function takes into account that no bubbles formed at T ′ > T i in the late patch.

Collapse and fractional abundance
For each choice of T i we can also calculate the contrast density in radiation This reaches a maximum shortly after late patch percolation, as the energy density in the background patch has began to become diluted a little earlier, while the late patch energy density is still constant due to the vacuum.After late patch percolation, δ decreases again as H late > H bkg leading to a faster redshifting.We define T δmax as the temperature at which δ is maximized, δ max ≡ δ(T δmax ).The smaller T i the larger δ max .We assume a patch collapses if δ max reaches a threshold value, δ c , for critical collapse.Calculations, some based on full general relativistic simulations, indicate 0.4 ≤ δ c ≤ 0.66 in the context of overdensities from inflation re-entering the Hubble horizon [23][24][25][26][27][28][29][30][31][32][33].The precise value of δ c depends on the shape of the overdensity, but its effects can be captured by considering the maximum of the so-called compaction function and its curvature, and thus be related to the power spectrum of curvature perturbations.Note the calculations have been performed assuming spherical symmetry and typically also isotropic pressure (although see [118]).We expect departures from these assumptions in the PT scenario, due to the nature of the Hubble patch just after bubble percolation.Note non-sphericity is expected to increase δ c [119][120][121][122][123]. The value of δ c will have an effect on which input parameters return f PBH = 1, but the strong GW signal will not be sensitive to our precise choice.In this work, we follow previous PT literature [56], and take δ c = 0.45 to aide comparison.The late patch nucleation temperature, T i , is thus from now on fixed by requiring the maximal density contrast, δ max , be equal to the critical collapse threshold, δ c .
We assume the collapse occurs at T δmax with corresponding time t δmax .The probability of the horizon size patch, which eventually collapses into a PBH, having no bubbles at T i is given by [47] P no bub = Exp − where the temperature and Hubble rate are those of the false vacuum, and the volume factor, represents the comoving volume of the Hubble sized patch at the start of the collapse.However, the calculation of δ for the critical collapse also assumes no larger bubbles -from background patches surrounding the late nucleating patch -enter into the collapsing volume before δ max is attained.Thus, it has been advocated [56], that the collapse probability is better estimated as where the volume factor is Here, in the evaluation of r(t δmax , t ′ ), Eq. ( 34), we use the background value for the scale factor.The eventual PBH mass is estimated as the energy inside the sound horizon of the collapsing patch at t δmax (e.g.see [124]) where c s = 1/ √ 3 is the sound speed (although included, ρ vaclate is negligible at t δmax ).The current study is limited to the above monochromatic estimate, we leave for future work the full derivation of a detailed mass spectrum taking into account the full critical collapse phenomenon [125][126][127][128][129]. Some preliminary estimates, showing the monochromatic estimate is a good one, are given in App. A.
We now wish to find, f PBH , the PBH-to-DM fraction at late times.To do this we note that PBH formation occurs a little after bubble percolation, around time t δmax when the false vacuum temperature is T δmax , and we denote the corresponding background temperature as T bkg form (note T bkg form ≈ T RH ≈ T infl with small differences found numerically and taken into account in our results).The ratio of densities at this time is where the ratio of Hubble rates is a correction which takes into account that M PBH is set by the late rather than background patch density.Later, at matter radiation equality, T eq ≃ 0.8 eV, we have ρ rad = ρ m = ρ DM (Ω DM +Ω B )/Ω DM ≃ 1.2ρ DM .Between PBH formation, at T ≃ T bkg form , and T eq , one has ∝ 1/a 3 dilution so From entropy conservation we have where g * s are the entropic degrees-of-freedom.The radiation instead redshifts as ∝ 1/a 4 , but it also receives reheating contributions, so Combing all the above we have where we have used g * (T bkg form ) = g * s (T bkg form ), and we remind the reader that g * s (T eq ) = 3.91 and g * (T eq ) = 3.36.

Application to example model
Having developed the above machinery, we now calculate the PBH fraction for our example B − L model, and search for parameter space in which f PBH = 1.In Fig. 3, we characterise PBH formation and the strongly first-order phase transition for a benchmark point for which v ρ = 10 3 TeV, y N i = 0.01, and g B−L = 0.2958.These values were judiciously chosen so as to return f PBH = 1.The PBH mass is M PBH ≃ 4.5 × 10 −12 M ⊙ .In the top left plot we show our full determination of Γ bub (T ) and compare it with the approximate form, Eq. ( 8), rewritten in terms of temperature Deviations between the full and approximate expression are visibly present at T = T i ≈ 54 GeV which returns the required contrast for PBH formation.We will eventually show the resulting modest differences in the β H required to achieve f PBH = 1 below.In Fig. 3 top right we first show the evolution of energy densities with the false vacuum temperature T .The dashed lines in the plot represent the vacuum, radiation, and total energy densities of an average Hubble patch, while the solid lines illustrate these same quantities for a late-nucleating patch.At sufficiently high false vacuum temperatures, the radiation energy densities for both types of Hubble patches align, and the same is observed with the vacuum energy densities.As nucleation starts within an average Hubble patch, its radiation energy density begins to increase relative to that of the late-nucleating patch.Once the percolation threshold of the true vacuum phase is achieved, there is a drastic drop in the vacuum energy density, leading the average background patch to become dominated by radiation.Concurrently, the energy density within this patch begins to dilute due to expansion, whilst the late-nucleating patch remains vacuum dominated.Once the percolation threshold is reached within the latenucleating patch, an excess of radiation energy density relative to the surrounding background Hubble patches follows, triggering the collapse into a PBH.
In the middle row of Fig. 3, we show the dependence of δ max on T i (right) and the dependence of the density contrast with the false vacuum temperature T (left).We observe how we start off from a homogeneous state at sufficiently high false vacuum temperatures.The density contrast drops to its minimum value approximately when bubbles start nucleating in the background Hubble patches and increases again when nucleation starts in the delayed patch.Ultimately, it reaches a maximum shortly after late patch percolation, as the energy density in the background patch has began to become diluted a little earlier, while the late patch energy density is still constant due to the vacuum.
Finally, in the bottom row we show the bubble distribution of the average background patch and the late patch with δ max = 0.45 at their respective percolation temperatures.The late patch features approximately twice the number of bubbles as the background patch.Due to the first bubbles being delayed, the distribution of bubbles in the late patch features a cut in the distribution, here at roughly half the Hubble length.In the right plot, we weight the distribution by R 3 to show the volume occupied by the bubble.Therein we see a significant volume is occupied by rare large bubbles, especially for the background patch, which explains why the late patch is filled instead by roughly twice as many small bubbles.Now we turn to Fig. 4. For a given choice of v ρ , we can scan over g B−L to find where f PBH = 1 (we always set y i = 0.01 as an example).The required g B−L to achieve f PBH = 1 is shown in Fig. 4 (top left) together with the corresponding value of the β λρ function (top right).
In the middle panel of Fig. 4 we show various temperatures characterizing the PT.The percolation temperature and also T δmax remains above QCD confinement, so we do not need to take into account QCD enhancement of Γ bub , present at lower temperatures [130][131][132][133].A smaller β λρ corresponds to a smaller vacuum energy difference between the two phases, which strengthens the phase transition [134].Furthermore, smaller values of g B−L imply a weaker running of couplings which, in turn, delays the transition time and makes the derivative of the bubble action S smaller [86].
As mentioned above, before solving the Friedmann equations, we first fix H = H false in our calculation of I(T ).To check this is a good approximation, we then iterate, using our H bkg determination to calculate a new I(T ), and use this to find an updated percolation temperature T p new .The ratio of the two is displayed in Fig. 4 middle right, showing only very small changes in T p , thus justifying the approximation.In Fig. 4 bottom we shown the resulting M PBH and various measures of the bubble radius at collision.
In the bottom row on the left, we show how the PBH mass decreases as we go to larger vevs, which corresponds to smaller physical Hubble volumes during the PT.In the right panel, we show the measures of the bubble radius at percolation, normalized to the Hubble length (mean, approximated mean, and radius corresponding to the mode of the bubble volume distribution).As the inverse timescale of the transition does not need to vary much from β H ≈ 8 in order to give f PBH = 1 over the entire allowed window, the typical bubble sizes when normalized to Hubble are also always similar.

Gravitational waves 5.1 The spectra
We next turn to the spectra we employ for the gravitational wave signal defined as where ρ GW is the energy density in GWs and ρ c is the critical density.
As spectra from dedicated studies are often given at production rather than today, we briefly review how to take into account the redshifting.We note that GWs redshift like radiation but The resulting PBH mass (left) and the bubble radii at collision (right), namely the mean radius R mean , the approximate radius, R approx = π 1/3 /β(T n ) [115], and the radius at which the bubble energy (and volume) distribution is maximized, R max [117].
are not reheated, and Ω rad (T RH ) ≡ ρ rad /ρ c ≃ 1 just after the PT.The above implies an amplitude today [66,73] where Ω GW * is the GW signal just after the PT, T 0 ≃ 0.23 meV is the CMB temperature today, and H 100 ≡ 100 km/s/Mpc.The frequency of the signal is redshifted as where f * is the frequency just after the PT.Concretely, for the reheating temperature, we take which follows directly from the first Friedmann equation.We are interested in supercooled PTs with highly relativistic bubble walls.We will use three determinations of the spectrum, showing our conclusions are hardly affected by the precise choice.
• The first estimate comes from (3 + 1) dimensional lattice simulations of thick wall bubble collisions by Cutting, Escartin, Hindmarsh, and Weir [75].The spectrum is found to be where the spectral shape is given by and the peak frequency is fthw = 32 mHz g * (T RH ) 100 In the above, we have used the central values of the fit provided in [75] for the thickest walls of their simulations, and taken the bubble diameter to be D bub H ≃ (8π) 1/3 β −1 H [115] in converting the spectrum in [75] to be in terms of β H . From Fig. 4 we see this diameter is in good agreement with the mean diameter extracted from the bubble distribution when using β H (T n ), which we subsequently use in calculating the GW spectrum.Note if we were to instead use 2R max , the diameter at which the bubble energy distribution is maximized, as advocated in [117], our GW signals would be stronger.
• The second estimate comes from the hybrid simulations by Vaskonen and Lewicki [77].
Therein, the anisotropic stress induced in a bubble collision is first determined in a (1+1) dimensional simulation.This is then used as a source at points at which walls collide in a (3 + 1) dimensional lattice simulation in the thin walled limit.The advantage here is that the lower dimensional simulation allows one to study the effect of non-trivial scalar gradients and associated gauge field production during the bubble collision.The simulations find differences in the spectra between the non-gauged [76] and gauged cases [77], of U(1) symmetry breaking.Of course we use the gauged example here, where the spectrum is with the shape, and peak frequency fhyb = 22 mHz g * 100 • Reassuringly also semi-analytic methods which avoid the need for running lattice simulations have been developed [72,73].We use the results of the bulk flow from Konstandin [73] with amplitude and spectral shape with peak frequency at fbulk = 21 mHz g * (T RH ) 100 Note all the above results are given in the limit α ≫ 1, v w ≃ 1, and the spectral shape functions have been normalized to return unity at their respective peak frequencies.As the finite cosmological horizon is not taken into account in the above simulations, we also impose the correct scaling Ω GW ∝ f 3 , for super-horizon modes at the time of the PT [135][136][137][138][139].These correspond to the frequenices at IR tail of the spectrum, below as measured today.In Fig. 5 we show the resulting GW spectra for two example parameter points with f PBH = 1.For practical purposes here, we see the GW spectra from the three estimates are all rather similar.They are all above foregrounds and detectable at some upcoming interferometer.
One word of caution is in order.The above estimates do not take into account the Hubble expansion during the transition itself.The issue has been studied in Ref. [78], which extended the analysis to an expanding background, but only for PTs in a radiation dominated background and using the older envelope approximation.The results indicated a suppression of an orderof-magnitude in the signal for PTs with β H ≈ 10.Note even with such a suppression, our GW signals would still be easily detectable, although this effect adds considerable theory uncertainty.To better gain a handle on the GW signal, the bulk flow analysis or the simulations will also have to be modified to take into account expansion, in particular for a vaccum dominated background.Our analysis further motivates such efforts.

Signal-to-noise ratio
We now wish to estimate the signal-to-noise ratio and also show foregrounds will not mask our signal.A similar method to what we describe in this section, up to some small updates, has also been used in [140].The signal-to-noise ratio is given by [141][142][143][144][145] where t obs is the observation time and Ω sens the sensitivity of the interferometer.For these quantities we take without deviation the choices made in [140] and use this to also calculate the power-law-integrated sensitivity curves shown in Fig. 5.Note the second and third terms in the denominator of the integral in Eq. ( 84) means the SNR is saturated for large GW signals.We do not want to count signals below astrophysical foregrounds as detectable, so as a simple and naive estimate, we define a foreground limited signal-to-noise ratio in which we only count GW signals above the astrophysical contamination.Namely, we compute, As estimates of the foregrounds we take: • The component of the GW galactic white dwarf binaries which are not subtractable after four years of LISA observations, with approximate form given in [146,147].The same quantity has also been estimated in [148], with a higher peak value, but at lower frequencies.For this latter estimate we use the analytic fit from [149].To be conservative we simply include both estimates.
• The upper value of the extra-galactic white-dwarf binary estimate from [150], a broadly similar estimate is derived in [151].
• The central value binary neutron star and binary black hole signal inferred from the observations in [152], extrapolated to lower frequencies assuming the f 2/3 scaling, as is appropriate for an stationary ensemble of binaries with circular orbits losing energy solely through GW emission [153].• The continuous super-massive binary black hole signal estimate of [151] will not mask our signal, but we display it in our plots.
To summarize, in Fig. 6 we show the latent heat α, β H , peak amplitude Ω GW (f peak ), the peak frequency f peak , SNR, and SNR FGL as a function of M PBH (the vev v B−L and coupling are fixed by M PBH and the requirement f PBH = 1).
In the top row of Fig. 6 we see, as expected, PBH formation occurs for α ≫ 1 and β H (T n ) ≲ 8, with slightly slower transitions required for more massive PBHs, in order to compensate the smaller enhancement in ρ PBH /ρ rad between PBH formation time and matter-radiation equality.We also compare to the required β H found for f PBH = 1 in Ref. [56], which used the approximate ansatz Eq. ( 8) for the nucleation rate.We obtain the required abundance with a slightly higher β H , because the bubble nucleation rate is somewhat suppressed in the full calculation compared to the approximate method.Note the collapse probability, P coll , is very sensitive to β H , thus a small shift in the latter can compensate for changes in Γ bub between the full and approximate expressions.
We thus confirm the approximate method is valid, at least in close-to-conformal potentials as studied here, which do not have a zero-temperature barrier.(If a zero temperature barrier exists, then there is the parameter dependent risk of having the field become permanently trapped in the false vacuum [154,155].Then PBH formation would not occur simply because the background patches themeselves do not percolate.Such models may require a more careful examination.Nevertheless, bubble sizes and hence GW signals are still necessarily large in successully percolating models which produce PBHs and feature more dramatic decreases in Γ bub below T n , see [55] for work along these lines.) In the middle row, we see the amplitude in Ω GW is large and almost constant over the range of M PBH , as β H is almost constant.The peak frequency covers three orders of magnitude, ranging over the sensitivies for LISA, BDECIGO, and ET, as anticipated through Eqs. ( 3) and (7).
In the bottom row of Fig. 6, we observe that the typical supercooled phase transitions in the classically conformal B − L model explaining all the dark matter in our universe in the form of PBHs will give rise to extremely strong GWs signals detectable by LISA, BDECIGO, and ET.These experiments will be able to probe the entire parameter space of the model producing f PBH = 1.Current limits do not constrain this PBH formation mechanism in the given range of M PBH [156,157].

Validity of assumptions 6.1 Rapid decay of the condensate excitations
In the above, we have assumed rapid reheating, so that T RH ≃ T infl .To check whether this is a valid assumption in our model, we can compare the decay rate of the scalar ρ with the Hubble rate at the end of the PT.The first partial width we consider is of ρ decaying to the EW Higgs doublet Here we take into account the thermal mass of the Higgs, which begins to become present as the bath is reheated, The above width is tiny, due to the suppression of λ 2 ρh , and possible kinematic factors as m H (T ) becomes large compared to m ρ .But we can also consider the partial width to the N i , As M Z ′ ≈ 10T RH , the thermal masses of the N i are negligible following the PT.For standard assumptions regarding the typical Yukawa coupling size -as inferred in from the type-I seesaw and the atmospheric mass squared difference of the active neutrinos -we can assume the N i decay rapidly compared to the Hubble rate, i.e. the strong washout regime [158].We can thus simply compare Eq. ( 87) to the Hubble rate to determine whether we can avoid a long lived ρ or other matter domination following the PT.In Fig. 7 left we show both the partial widths compared to the Hubble rate for our parameter choice, showing our assumptions of rapid decay, giving T RH ≃ T infl is justified.Note the Z ′ will also rapidly decay to SM fermions.

Pressure on the bubble wall
For supercooled transitions, such as those studied here, leading order pressure on the wall as particles gain their mass [159], is generically insufficient to stop the wall from accelerating.The presence of gauge bosons which change their mass across the wall, however, leads to a NLO retarding pressure which increases linearly with γ, the bubble wall Lorentz factor [160].We therefore check whether the vacuum energy density is converted into the wall energy, or whether sufficient soft-gauge bosons are produced at the wall, for the energy to instead be stored in the latter before wall collision.Such a situation has been modelled in [161].Whether this would affect the GW signal appreciably is a separate question, macroscopically both cases seem to be captured by the bulk flow model (e.g.see [161] and discussion in [162]).Soft gauge boson production leads to a retarding pressure on the wall P NLO ∝ γT 3 p M Z ′ [160,163].Here we use the leading-log determination for this term given in [164], where κ ≈ 4, (Q eff B−L ) 2 = 40 is an effective charge factor summing over the bath degrees-offreedom with 3/4 (1) weighting for fermions (scalars), α B−L ≡ g 2 B−L /4π, and we have used the thermal mass of the Z ′ as an IR cutoff in the logarithm.The maximum Lorentz factor is obtained at collision for an effectively run-away wall, i.e. approximately zero retarding pressure, and is given by [164] γ ≈ 1 3 where R nuc ≈ 10/T n is the bubble radius at nucleation and R coll ≈ π 1/3 /β is the bubble radius at collision (we have checked both analytic estimates through numerical determination).For consistency of the above γ estimate we require P NLO < Λ vac .We calculate these quantities and display the results in Fig. 7 right, showing the wall is in the effective run-away regime at collision.

Reliability of the effective potential
We now briefly discuss the reliability of the approximations we have made in deriving V eff .To answer this question requires the application of more refined techniques such as the renormalization group improved effective potential [111,165].Indeed, in the context of the particle-DM dilution mechanisms, there are known examples in which going to the RGE improved V eff changes the qualitative outcome of the mechanism itself [166].This is due to changes to the details of the nucleation, percolation, and reheating temperatures, although the PTs remain very strong.Such an investigation is left for future work.

Comments on leptogenesis
In Fig. 8 we show the dilution factor due to the entropy production following the PT.As the dilution factor is quite large, ∼ 10 7 − 10 9 , it seems unlikely that we could use a mechanism which relies on the dilute false vacuum plasma as a source of the baryon asymmetry.For example, the mass gain mechanism when the N i suddenly gain their mass during the PT and then decay to generate the BAU [167][168][169][170][171]. Similar conclusions most likely also hold for the related mechanism in [172].Electroweak style baryogenesis mechanisms suffer from the same problem in the current context, together with a completely negligible yield for v w ≃ 1, due to lack of diffusion back into the symmetric phase [173].
If instead, heavy particles are produced in the reheating process itself, following along the lines of [174,175], then these may have a larger abundance compared to the entropy density, and act as a source of the asymmetry.But note our bubble collisions are expected to be inelastic, which suppresses heavy particle production [176].
Anyway, for our choice of example Yukawa coupling y i = 0.01, the asymmetry will be generated by the usual thermal leptogenesis around T ≈ M N after the PT (for our parameter choices T RH ≃ 10M N i ).In either case, resonant leptogenesis is generically required [177][178][179][180], as the heavy neutrinos are much below the Davidson-Ibarra bound [181] (also see [182]).Note the heavy neutrinos are efficiently produced in the early universe due to their Yukawa coupling to ρ and their gauge interactions.

Conclusion
In this work we have studied the GW signal from PBH formation during a PT.We have assumed f PBH = 1, which limits the range for the peak frequency of the GWs, and calculated the required underlying parameters of the model to obtain the PBH abundance.We used the properties of the PT to find the expected GW signal from bubble collisions over the allowed parameter space, showing it can be detected by BDECIGO or a combination of LISA and ET, independently of the precise M PBH , over the entire range allowed for f PBH = 1.Substituting one of the experiments by a broadly equivalent one will, of course, not affect the conclusions.An important caveat is that the GW predictions for strong phase transitions may overestimate the signal because expansion during the transition is not taken into account.This motivates further refinement of such estimates, along the lines of [78], but extended to vacuum dominated transitions and beyond the envelope approximation.We note a suppression in Ω GW by an order of magnitude for β H ≈ 10, would still return SNRs of 10 2 and above over the f PBH = 1 parameter space, so prospects remain very promising.
The classically scale invariant B − L model we considered does not contain an automatically stable DM candidate, although given a small enough Yukawa coupling, N 1 could play the role of the possibility of a ∼ keV scale sterile neutrino (provided an appropriate mechanism was also in play in order to set its relic abundance correctly, as it carries gauge charge).Such sterile neutrino scenarios are by now by now heavily constrained [183].Additional field content can of course be included to act as DM [184].We have instead proposed using the PBH generation mechanism from the strong PT realised with the close-to-conformal potential to obtain the observed Ω DM .Leptogenesis in this model should occur sometime after the PT, through resonant enhancement of the CP violation, as in the usual low-scale thermal leptogenesis [177][178][179][180]. Furthermore, it would also be of interest to re-examine models of particle DM, in which the relic abundance relies on supercooled PTs [87,140,166,[185][186][187], to check compatibility with PBH production.
Our calculations in a specific model, taking into account the full temperature dependence of the bubble nucleation rate, have shown only modestly small changes to the inverse timescale of the transtion β H required to achieve f PBH = 1, compared to model independent calculation which used an approximate ansatz for the nucleation rate.As β H sets the typical bubble size at collision and therefore also controls the properties of the GW signal, our results regarding the sensitivities of future detectors to such a late patch mechanism are independent of the underlying particle physics model.
Models of slow roll inflation producing PBHs also result in detectable GWs, in the same frequency range as the PT, due to the enhanced power spectrum on small scales leading to significant anisotropic stress [34][35][36][37][38][39].The additional power also present on small scales from the PT, would also induce anisotropic stress and a further GW source.Accordingly, this paper justifies the further development of computational techniques in order to pin down the expected PBH spectrum in the late patch mechanism (beyond the monochromatic approximation primarily used here), the bubble collision signal, together with any additional GWs at second order in perturbation theory, in light of testing f PBH = 1 in upcoming interferometers.Along with the GW signal, the PBHs can be searched for using improved lensing studies, through future MeV telescopes [188,189], and the 21-cm absorption signal [190,191].

A Towards a non-monochromatic mass distribution
In the main text, we used a monochromatic approximation for the PBH mass, We also found where (Ω DM + Ω B )/Ω DM ≃ 1.2.We then used where M Hor = M 2 Pl /(2H bkg ) is the energy inside the Hubble horizon of a background patch.In the monochromatic approximation P coll (T i ), given in Eq. (60b), is evaluated at T i chosen so that δ max = δ c .
We now assume instead critical collapse, as is known from numerical codes relevant for the inflation scenario, in which for δ m > δ c the PBH mass is found to be [125,192] where 1 ≲ K ≲ 10 is an efficiency factor, and γ c ≃ 0.36 is the critical exponent assuming radiation domination (both found from simulations).In the inflationary scenario, δ m is the overdensity as measured at the maximum of the compactification function, at cosmological horizon crossing time, e.g.see [27].We now make assume something similar holds in the PT mechanism and replace δ m → δ max in the above equation.We warn we are making a big leap making such an assumption for critical collapse in the PT mechanism, and do so simply out of curiousity regarding the resulting PBH spectrum, mainly as an academic exercise.Note in using Eq. ( 94) we also do not include the small correction taking into account the different Horizon masses of the collapsing and background patches as we have done for our monochromatic approximation in Eq. ( 91).From Eq. ( 94) we see that patches with δ max = δ c fail to form a black hole, but instead a whole spectrum of PBH masses is generated from patches with larger δ max , coming from progressively smaller T i .The probability of keeping large regions bubble free also falls as T i is reduced.(Roughly speaking, P coll is the probability of keeping a sufficiently large region bubble free sufficiently long, it therefore includes the probability of staying bubble free even longer.In other words it acts as a cumulative probability function.)A plot of P coll (T i ) is shown in Fig. 9.
We now wish to calculate the resulting PBH mass spectra assuming Eq. ( 94) holds.The first spectrum typically considered in the literature is normalized to the PBH density and is conventionally defined as [192] Ψ(M PBH ) ≡ In our scenario we have where the normalization factor is given by Note dP coll /dM PBH < 0 given our definition of P coll in Eq. ( 60b), but the factor appears both in the numerator and the denominator, N , so Ψ is positive.Integrating Ψ(M PBH ) over the entire PBH mass range yields unity.In the above, we have related the spectrum to the collapse probability calculated in Eq. (60b), and we have used that M Hor and T bkg form are approximately independent of M PBH for PBHs originating from overdensities with a preferred length scale (in the current context of the late patch mechanism this is simply the Hubble length during the PT).This is further justified by the strongly falling P coll with decreasing T i , as shown in Fig. 9.
The second commonly used spectrum is instead normalized to the observed DM density and is defined as [192] so that In our scenario we find where we have again used that M Hor and T bkg form are approximately independent of M PBH for PBHs originating from overdensities with a preferred length scale.We now include an overall negative sign because dP coll /dM PBH < 0. For either of the spectra of interest, we can evaluate PBH respectively [192], coming from values of T i close to which δ max = 0.45.These IR scalings can be understood from Eqs. ( 96) and (100c).
the last two differential factors numerically as functions of M PBH , and thus find the desired quantity.
We now return to our example parameter point used in Fig. 3 and calculate the PBH mass spectra assuming the critical collapse phenomenon holds.The resulting PBH spectra are shown in Fig. 10.Clearly the spectra are strongly peaked at M PBH ≈ c 3 sound M Hor ≃ 0.2M Hor , the value assumed for the monochromatic approximation.Furthermore, the total integrated PBH density normalized to the observed DM density remains largely unchanged.Thus the monochromatic approximation is justified.We can also compare the resulting spectra to existing limits as done in Fig. 11, but only to gain a rough idea, as the limits are strictly only valid in the monochromatic approximation.Still it is reassuring that the spectra are strongly enough peaked that the IR tails do not intersect with the low mass constraints, as long as ⟨M PBH ⟩ ≳ O(10 −16 ) M ⊙ .The UV cut-off is much steeper and so even less of a concern.
Figure 11: The PBH spectra for two example parameter points, both with y N i = 0.01 and K = 1, compared to limits on monochromatic spectra from the CMB [11,12], extra galactic background light [13], galactic gamma ray background [15], Voyager e ± [16], and Subaru Hyper Suprime-Cam (HSC) lensing [20].The GWs from the same two parameter points are shown in Fig. 5.
Even assuming the critical collapse relation holds for the late patch mechanism, the above is still not a complete estimate of the full PBH spectrum.This is because: (i) we have ignored the possibility of PBHs being formed at slightly different times with different M Hor , and (ii) the collapsing volume is not limited to Eq. (61b), but may take a whole range of values (for us, horizon and superhorizon at the time of the PT, but one may also imagine that subhorizon patches with larger densities can also collapse, as in Ref. [55]).Given the theoretical uncertainties involved in applying the critical collapse to the PT scenario, however, together with the computational time required for the numerical evaluation, we have not yet evaluated this full expression.Nevertheless, the peaked nature of the PBH spectrum is expected to hold also in this more general case, because the probability of obtaining larger V coll will also fall rapidly.

Figure 1 :
Figure 1: The evolution with temperature of the effective potential for a benchmark point with g B−L = 1, y N i = 1 and v ρ = 100 TeV.The barrier separating the two minima means the evolution from the false to the true vacuum corresponds to a first order phase transition.The barrier is always present when the finite temperature corrections are added to the radiatively induced symmetry breaking potential.

Figure 2 :
Figure 2: Left: Temperatures characterising the phase transition as a function of the gauge coupling.We display the following temperatures: critical, T c , nucleation , T n , percolation, T p , and the temperature at which thermal inflation starts, T infl .Right: The dimensionless phase transition latent heat, α, and inverse timescale normalised to the Hubble rate, β H .

Figure 3 :
Figure 3: Characterisation of several important quantities for PBH formation and the phase transition dynamics for a bechmark scenario with v ρ = 10 3 TeV, y N i = 0.01, and g B−L = 0.2958 that predict f PBH = 1, with M PBH ≃ 4.5 × 10 −12 M ⊙ .We depict the behaviour of the nucleation rate comparing the full calculation of the nucleation rate with the approximation often used in the literature [47, 56] (top left), energy densities (top right) and density contrast (middle left, with T i ≃ 53.6 GeV) as a function of the false vacuum temperature T .We also show the density contrast as a function of T i (middle right).Finally we display the bubble distribution at percolation (bottom left), where RH is the bubble size as a fraction of the Hubble length.Weighting the distribution by R 3 we see the differences in the large bubbles between the background and late patches (bottom right).

Figure 4 :
Figure 4: Top: The required gauge coupling (left) and associated beta function (right) giving f PBH = 1 as a function of the B − L vev, v ρ .Shaded regions are excluded by PBH constraints.Middle: Various temperatures characterizing the PT as a function of v ρ .Note the new percolation temperature T pnew following our iterative check shows only a per mille correction with respect to T p (right).Bottom:The resulting PBH mass (left) and the bubble radii at collision (right), namely the mean radius R mean , the approximate radius, R approx = π 1/3 /β(T n )[115], and the radius at which the bubble energy (and volume) distribution is maximized, R max[117].

Figure 5 :
Figure 5: GW spectra for two PTs returning f PBH = 1.For each of the PTs, we show three estimates of the spectra -given the macroscopic PT parameters -available in the literature.Also shown are power-law-integrated sensitivity curves (SNR = 10) for current and future interferometers together with astrophysical foregrounds.Other key parameters associated with these example PTs can be read off Fig. 4.

Figure 6 :
Figure 6: Top: The latent heat α and inverse timescale of the transitions β/H returning f PBH = 1as a function of M PBH .Also shown is a comparison to the results of[56], which uses the approximate Γ bub ansatz, so some differences are expected.Middle: The peak amplitude and peak frequency of the three GW spectrum models as a function of M PBH .Bottom: the signal-to-noise ratio (SNR), and signal-to-noise ratio above foregrounds (SNR FGL ), as a function of M PBH .The signal can be detected by BDECIGO or a combination of LISA and ET over the entire allowed range of M PBH with f PBH = 1.The shaded band encompasses the three different determinations of the GW signal.Additional theory and experimental forecasting uncertainty means the true error band is larger.

Figure 7 :
Figure 7: Left: The decay rates compared to the Hubble rate.Right: The LO and the γ enhanced NLO retarding pressure on the wall compared to the driving pressure, i.e. the vacuum energy difference.

Figure 8 :
Figure 8: The entropy production (dilution) factor following the PT.

Figure 9 :
Figure 9: Left: The cumulative probability P (T i ) for the example parameter point, v ρ = 10 3 TeV, y N i = 0.01, and g B−L = 0.2958.Right: For convenience we again show δ max (T i ) for the same parameter point, where δ max = 0.45 for T i = 53.6GeV (for a larger range of T i see Fig. 3, middle row, right).

Figure 10 :
Figure 10: The PBH spectra for the example parameter point, v ρ = 10 3 TeV, y N i = 0.01, and g B−L = 0.2958, now taking into account that δ max is a function of T i , and assuming the critical collapse scaling with K = 1.Upon integration we find ⟨M PBH ⟩ ≃ 5.4 × 10 −12 M ⊙ and an integrated total f PBH ≃ 1.05, compared with M PBH = 4.5 × 10 −12 M ⊙ and f PBH ≃ 1.02 using the monochromatic approximation.The grey dashed lines show the expected IR scalings, Ψ(M PBH ) ∝ M 1/γc PBH and f PBH (M PBH ) ∝ M 1+1/γc s (T eq ) g * (T eq ) * s (T eq ) g * (T eq )