Primordial Black Holes from Conformal Higgs

Scale-invariant extensions of the electroweak theory are not only attractive because they can dynamically generate the weak scale, but also due to their role in facilitating supercooled first-order phase transitions. We study the minimal scale-invariant $U(1)_{\rm D}$ extension of the standard model and show that Primordial Black Holes (PBHs) can be abundantly produced. The mass of these PBHs is bounded from above by that of the moon due to QCD catalysis limiting the amount of supercooling. Lunar-mass PBHs, which are produced for dark Higgs vev $v_\phi\simeq 20~\rm TeV$, correspond to the best likelihood to explain the HSC lensing anomaly. For $v_\phi\gtrsim 400~\rm TeV$, the model can explain hundred per cent of dark matter. At even larger hierarchy of scales, it can contribute to the $511~\rm keV$ line. While the gravitational wave (GW) signal produced by the HSC anomaly interpretation is large and detectable by LISA above astrophysical foreground, the dark matter interpretation in terms of PBHs can not be entirely probed by future GW detection. This is due to the dilution of the signal by the entropy injected during the decay of the long-lived $U(1)_{\rm D}$ scalar. This extended lifetime is a natural consequence of the large hierarchy of scales.

In this letter, we investigate whether the electroweak phase transition (EWPT), minimally extended by new physics, can generate PBHs.Apart from its second derivative around ⟨H⟩ ≃ 246 GeV measured in 2011 [32,33], the shape of the Higgs potential remains unconstrained [34].In 1976, Gildener and S. Weinberg [35] suggested that the Higgs mass could vanish at tree level and be dynamically generated when quantum loop corrections, introduced by Coleman and E. Weinberg (CW) in 1973 [36], drive the Higgs quartic to negative values in the infrared.The original model has been ruled out by the large top mass [37], necessitating the introduction of additional bosonic fields into the Standard Model (SM) to make it viable.One of the simplest scale-invariant extensions of electroweak theory involves a complex scalar field ϕ, neutral under the SM but charged under an additional U (1) D gauge group [37].This model has attracted considerable attention due to its small number of new parameters and potential implications for dark matter, leptogenesis and neutrino oscillations [38][39][40][41][42][43][44][45][46][47].Phase transitions driven by a Coleman-Weinberg potential plus thermal corrections: are known to be strongly first order (1stOPT) and supercooled [48,49].We have Taylor-expanded thermal corrections around the symmetric minima ϕ = 0. We have introduced the beta function β λ of the quartic coupling λϕ 4  and the coupling g with particles of the plasma.For the U (1) D model with gauge coupling g D studied in this letter, one has β λ ≃ 3g 4 D /8π 2 and g 2 ≃ g 2 D /12.Spontaneous breaking ⟨ϕ⟩ ̸ = 0 is communicated to the Higgs field H through the mixing λ|H| 2 ϕ 2 .The scalar ϕ can be interpreted as the "Higgs of the Higgs" [50].Thanks to the logarithmic potential, the temperature T n at which bubbles nucleate and percolate is exponentially suppressed for β λ ≪ 1: with respect to the critical temperature T c when the 1stOPT becomes energetically allowed.S c is the bounce action at percolation.Coleman-Weinberg dominance in the potential in Eq. ( 1) additionally implies that the bounce action cost S 3 /T -which controls the bubble nucleation rate per unit of volume Γ V = T 4 exp(S 3 /T ) -decreases only logarithmically with the temperature.This implies that the completion rate of the PT: can be relatively small for β λ ≪ 1.
We investigate the minimal U (1) D scale-invariant extension of the electroweak theory.We study the tunneling dynamics both numerically [51] and analytically taking into account QCD catalysis and collider constraints.We find that the two necessary conditions for the formation of PBHs in observable quantities derived in the companion paper [21] -significant supercooling (T n ≪ T c ) and a slow rate of completion (β/H ≲ 7) -can be satisfied as soon as the Dark Higgs Vacuum Expectation Value (VEV) is larger than v ϕ ≳ 10 TeV.The gravitational wave detection prospects are calculated in all the parameter space, accounting for astrophysical foregrounds.We account for the matter era following the phase transition resulting from the dark Higgs becoming long-lived at large VEV.

II. SCALE INVARIANT EXTENSION OF THE ELECTROWEAK THEORY
The scalar potential -We extend the SM with a complex scalar field Φ charged under a hidden U (1) D dark photon V µ with coupling constant g D [37]: with D µ = ∂ µ + ig D t a V a µ .Supposing Bardeen's conjecture that scale-invariance is a fundamental symmetry of Nature at tree level [52], the only allowed operators in the scalar potential are: where λ hϕ is the mixing coupling with the SM Higgs H.In the unitarity gauge, we can parameterize the scalar fields as: Scale-invariance is violated by quantum interactions encoded in Renormalized Group Equations (RGE), see e.g.[43].In App.B, we show that upon assuming the ordering g 2 D ≫ λ ϕ , λ hϕ , the 1-loop corrected potential is well approximated by its leading-log term: We use the full RGE improved coupling λ ϕ (ϕ) in our numerical computations.The Coleman-Weinberg potential in Eq. ( 7) generates a VEV for the U (1) D field ⟨ϕ⟩ = v ϕ , which through the mixing λ hϕ h 2 ϕ 2 generates a VEV for the Higgs ⟨h⟩ = v EW .Fixing the Higgs VEV to v EW = 246 GeV and its mass to m h = 125 GeV implies: (8) At finite temperature the potential for ϕ receives thermal corrections which under the ordering g 2 D ≫ λ ϕ , λ hϕ are dominated by the dark photon V µ .The 1-loop and Daisy contributions read [53]: with m 2 V (ϕ) = g 2 D ϕ 2 and Π V = g 2 D T 2 /6 [54].QCD catalysis -In App.B, we show that the motion of the Higgs field can be safely neglected during the tunneling if v ϕ ≳ TeV.An exception arises if the universe super-cools below the temperature of QCD confinement T n ≲ Λ QCD .Then the top chiral condensate induces a negative linear slope in the Higgs direction − yt √ 2 ⟨t L t R ⟩ h [49], see also [55][56][57][58][59][60][61].The Higgs acquires a VEV: which through the mixing λ hϕ induces a negative mass for the hidden scalar ϕ: Since both the Higgs quartic λ h and the top Yukawa y t diverges below Λ QCD , it is not possible to predict the precise value of ⟨h⟩ QCD [57,60].Waiting for further studies, we approximate the temperature dependence of QCD effects with a Θ function and fix ⟨h⟩ QCD = 100 MeV.
Bubble nucleation -The thermal corrections in Eq. ( 9) generates a barrier between the minima in ϕ = 0 and ϕ = v ϕ .Thermal fluctuation can drive the field over the barrier with a rate per unit of volume [62]: where S 3 is the O 3 -symmetric bounce action.We wrote a C++ undershooting-overshooting algorithm CCBounce [51] to find the tunneling trajectory and compute the bounce action S 3 .We found it to be well approximated by the thick-wall formula: T QCD is the temperature when the QCD contribution in Eq. (11) becomes larger than the thermal barrier of ϕ and below which the EW-PT must necessarily complete.We refer to App.B for more details on the bounce action calculation and for a comparison of the thick-wall formula with the proper numerical treatment.
Nucleation happens when the decay rate in Eq. (B27) multiplied by the causal volume becomes comparable to the Hubble expansion rate: Γ V (T n ) ≃ H 4 (T n ).Plugging Eq. ( 13), we deduce the nucleation temperature: where we neglected the QCD contribution for clarity.We use a proper root solving algorithm to find T n in our plots.In the absence of QCD, effects, the bounce action S 3 /T flattens at low temperature, causing the absence of solution for Γ V (T n ) ≃ H 4 (T n ) below the temperature [63][64][65]: where H n ≡ H(T n ) and M given by Eq. ( 13).Instead, QCD confinement generates a sharp decrease of S 3 /T down to zero below the temperature T QCD .
We can now evaluate the rate of PT completion β ≡ d log Γ V /dt at nucleation.Above QCD confinement T n ≫ T QCD , it decreases logarithmically with the temperature β/H n ≃ −4 + 127 Instead, around  The PT completion rate β/H becomes lower as the amount of supercooling increases, until QCD effects kick in around Tn ≃ TQCD in Eq. (13).The lines are obtained from varying gD.The parameter α along the horizontal axis is the latent heat fraction.
The latent heat ∆V ≃ β λ v 4 ϕ /16 dominates the energy density of the universe, generating an inflationary period, below the temperature T eq : Reheating after matter domination -Supercooled PTs proceed with ultra-relativistic bubble wall Lorentz factor γ ≫ 1 [66][67][68][69].Bubble growth convert the latent heat into scalar field gradient stored in the walls and plasma excitation resulting from the friction operating on the walls.The fraction stored in the former at percolation is given by: where γ coll is the wall Lorentz factor at collision accounting for friction [68], while γ max is the same quantity neglecting friction.We refer to App.C for the details.The main dark Higgs decay channel is into two Higgs ϕ → hh.However the decay rate is suppressed by the small mixing λ hϕ : where m ϕ = β λ v ϕ is the scalar mass and is the Hubble rate at percolation.The universe after percolation contains a matter component occupying a fraction κ coll in Eq. ( 18) of the total energy density.According to Eq. ( 19), this matter component is long-lived for v ϕ ≳ 10 3 TeV.
It dominates the universe below the temperature T dom = κ coll T eq , hence generating an early matter-domination (MD) era, which ends with the decay of ϕ at the temperature T dec = 1.2 M pl Γ ϕ /g 1/4 * .The decay of ϕ reheates the SM and dilutes the abundance of any relics decoupled from the SM with a factor: As discussed in the next section, this will have a minor impact on the PBH detectability but a major impact on the GW one.In an adiabatic universe, supercooled phase transition produce gravitational waves within the frequency window of future interferometers LVK run O3 [70], run O5 [71], ET/CE [72][73][74] and LISA [75][76][77][78] with an amplitude larger than astrophysical foregrounds from white dwarfs binaries and compact binaries [70,76,[79][80][81].Instead, in the minimal scale-invariant extension of the SM, the long lifetime of the scalar field dilute all the GW signals at large dark Higgs VEV v ϕ ≳ 10 4 TeV.Instead the PBH abundance is only slighly impacted by the matter era.See App.C for more details on the experiments reach and expected astrophysical foregrounds.The purple region shows the microlensing constraints (ML) from HSC/subaru telecope [82,83].The blue U-shaped region can explain the HSC event [84][85][86].The boundary of the brown region can explain dark matter.
Non-observation of energy injection from Hawking radiation in CMB [96][97][98] and BBN [90,99,100] exclude the yellow and green regions.See [8] for a review of PBH constraints.The GW constraints shown in light blue disappear at large v ϕ ≳ 10 4 TeV due to entropy injection following the decay of the long-lived dark Higgs.The associated dilution factor D is indicated with blue lines labelled "D".Part of the PBH DM region can not be probed with GWs.On the left of the dashed black line, the PT is catalysed by QCD confinement.In the blue region in the top-left corner, the universe can not tunnel without "QCD effects", due to ΓV = H 4 having no solution Tn.However, the QCD catalysis temperature given in Eq. ( 13) is smaller than the Hubble scale TQCD ≲ Hn, implying that the tunneling dynamics receives large correction from "gravity effects" [101,102], which we leave for future research.The gray region displays collider constraints on additional Higgs boson [103][104][105][106], see Fig. 10 in App.B for more details.On the right of the very thin blue line, bubble wall friction convert most of the the latent heat into the plasma, which reduces the energy fraction stored in the scalar field and therefore reduces the dilution factor D.

III. PHENOMENOLOGY
Primordial black holes -Due to the stochastic nature of thermal tunneling, distinct causal patches percolate at slightly different times.Patches which percolate later hold the latent heat longer in the form of vacuum energy while their neighbors already released it into radiation-like energy density.When the latest patches finally percolate, they are overdense which, beyond a given energy density threshold δ c , lead to their collapse into PBHs.We set the over-density threshold to δ c = 0.45, which assumes a Mexican hat density profile during radiation domination [107][108][109][110]. Different profile shapes would result in different values within the range δ c ∈ [0.40, 0.67] [111,112].The determination of the density profile shape as well as the precise equation of state right after percolation are left for future studies.In particular, in the runaway bubble wall regime with κ coll = 1, PBH formation could occur during a matter era with a significantly lower critical threshold [113], leading to an enhanced PBH abundance.The PBH mass is given by the mass within the sound horizon at the time of percolation: where M d ≃ 3.7 × 10 −8 M ⊙ is the moon mass.The calculation of the PBH abundance has been completed in the companion paper [21] and reported in App. A. Normalised to the dark matter relic density, it gives: where P coll is the fraction of causal patches collapsing into PBHs, see App. A. We added the dilution factor D in Eq. ( 20) to account for the entropy injection resulting from the matter era after the PT.In Fig. 2, we show that supercooled PTs with β/H ≲ 7 would produce PBHs detectable by on-going experiments, and that the effect of the dilution factor D -opaque vs transparent -is very limited.PBHs formation in the scaleinvariant U (1) D extension of the SM is shown in Fig. 3 together with astrophysical and cosmological constraints [8].
For comparison, in Fig. 9 we present the same figure with QCD catalysis neglected.This is equivalent to setting ⟨h⟩ QCD in Eq. ( 10) to 0 instead of 100 MeV.
Gravitational waves -As the supercooled PT completes, the large latent heat ∆V is converted into scalar field gradients with an energy fraction κ coll .This drives bubble walls to ultrarelativistic Lorentz factors [66][67][68][69].A fraction 1 − κ coll is converted into extremely thin shells of ultra-relativistic particles [67][68][69][114][115][116] or ultra-relativistic shock-waves [117][118][119][120], which trail bubble walls [67][68][69][114][115][116].The GW spectrum from thin shells of stress-energy momentum tensor can be described by the "bulk flow" model [121][122][123][124][125]. We refer the reader to App.C for the details.For large VEV v ϕ ≳ 10 3 TeV, corresponding to small Higgs mixing λ hϕ ≲ 10 −5 , the PT is followed by a matter domination era, further redshifting the GW signal according to: where D is given by Eq. ( 20) and with additional spectral distortion in the causality tail [126][127][128][129][130]. We refer to Fig. 11 in App.C for more details and for a visualisation of the GW spectra.The impact of the dilution factor on the GW detectability is visible in Figs. 2 and 3. We can see that the existence of the matter era -which is a consequence of the smallness of the Higgs mixing λ hϕ ≪ 1 in Eq. ( 8) -limits the GW detectability to v ϕ ≲ 10 4 TeV, corresponding to PBH masses M PBH ≳ 10 −14 M ⊙ .Part of the region explaining DM with PBHs might never be probed with GWs, in contrast to the naive expectation [131].As a comparison, the same figure with matter era effects neglected is provided in Fig. 8 in App.B.

IV. CONCLUSION
In this letter, we have demonstrated that particle physics models featuring scalar potentials dominated by logarithmic terms, similar to the Coleman-Weinberg potential in Eq. ( 1), are promising theoretical frameworks for producing PBHs from 1stOPTs.By focusing on the U (1) D scale-invariant extension of electroweak theory, we have found that PBHs with lunar masses can be produced for dark Higgs vev v ϕ ∼ 20 TeV.This provides a plausible interpretation for the candidate event from the HSC microlensing search [85,86,132].At larger vev v ϕ ∈ 4 × [10 2 , 10 5 ] TeV, PBHs can explain 100% of the observed Dark Matter (DM) relic den-sity [6][7][8][9].Another consequence of the larger hierarchy of scales is the suppression of the dark Higgs decay width as Γ ϕ ∝ (v EW /v ϕ ) 4 .The decay after a matter era can generate a tremendous amount of entropy dump which dilutes the gravitational wave (GW) signal.GW constraints disappear for v ϕ ≳ 10 4 TeV, leaving PBHs as the sole detectable evidence for a scale-invariant extension of the standard model across a wide range of energy scales.This includes areas that could account for 100% of DM in terms of PBHs.As pointed in [133][134][135][136][137], part of the PBH DM region will be explored with future telescopes [138][139][140].The detection of PBH dark matter without GW counterpart could be a smoking-gun for a conformal Higgs scenario.However, it would be degenerate with standard PBH formation scenarios if they are followed by entropy injection.This degeneracy could be resolved by measuring the mass distribution of PBHs [22,27] or by using cosmic string archaeology to decipher the equation of state, duration and epoch of the entropy injection phase [141][142][143][144][145][146][147].Two remarks are important to highlight before closing.Firstly, in contrast to logarithmic potentials, traditional polynomial potentials typically do not satisfy the two necessary conditions for PBH formation which are latent heat domination (T n ≪ T c ) and slow rate of completion (β/H ≲ 7).This is because the bounce action decreases too rapidly with temperature, following a power-law [148] instead of the logarithmic decrease in Eq. ( 13).Secondly, it is important to emphasize that logarithmic potentials are not exclusively predicted by weakly-coupled Coleman-Weinberg models.They also arise in nearly-conformal strongly-coupled models featuring a light dilaton [58,63,64,114,[149][150][151][152][153][154] and their weakly-coupled 5D holographic dual [56,59,[155][156][157][158][159][160][161][162][163][164][165][166][167][168].This letter opens new avenues for investigating the existence of PBHs in our universe using nearly scale-invariant particles physics models.We report here the derivation of the PBH abundance during supercooled PTs which has been proposed in the companion paper [21].It accounts for fluctuation in the nucleation time of the first bubble in the full past lightcone of the collapsing patch and consider the collapse to occur after percolation, in a radiation-dominated universe.In contrast, other approaches proposed in Refs.[14,19,20,26] are restricting the collapsing patch to remain 100% vacuum dominated until collapse, leading to a lower PBH abundance, while in Ref. [16] nucleation is not accounted in the entire past light-cone of the collapsing patch, leading to a larger PBH abundance.The derivation in the companion paper [21] has been reproduced in details in Ref. [22].While in the companion paper [21] a monochromatic mass distribution is considered, Ref. [22] extends the analysis by calculating the PBH mass distribution.Ref. [27], which appeared after the first version of this work, improves the treatment of [21,22] by accounting for fluctuations in the nucleation time of the first j c bubbles in a given causal patch, with j c being arbitrary.Refs.[21,22] can be viewed as the case j c = 1.We leave the quantitative analysis of the importance of those corrections to the PBH abundance for future works.
During a supercooled PT, the latent heat ∆V ≃ β λ v 4 ϕ /16 can dominate the energy density of the universe, leading to an inflationary stage below the temperature T eq : which also corresponds to the maximal reheating temperature after the PT, up to ratio of g * .Inflation ends when bubbles percolate around a temperature T n given by: It is convenient to approximate the tunneling rate with the bounce action Taylor-expanded around the nucleation  time t n at first order: where the expression for Γ 0 follows from Eq. (A2).The validity of the β-parametrization in Eq. (A3) is discussed below Eq. (B31).
During the nucleation phase, the vacuum energy ∆V is converted into kinetic energy of bubble walls which redshifts akin radiation.The time and position at which a bubble nucleates are stochastic variables.Hence the percolation time of each causal patch fluctuates with respect to the background.While the average universe has already percolated and is cooling like a radiation-dominated universe, late-bloomers -patches which start nucleating later -maintain an almost constant energy density until they finally percolate and convert their vacuum energy into radiation.Hence, causal patches which nucleate the latest inflate the longest and become the densest.Late-blooming region ("late") collapse into PBHs if the contrast radiation density with respect to the background ("bkg") is larger than: where t ni is the time at which a late-bloomer nucleates its first bubble.It is also the time until which it remains 100% vacuum-dominated.The overdensity generated by a supercooled 1stOPT is expected to source an adiabatic curvature perturbation ζ to the FLRW universe: where we have assumed a spherically-symmetry perturbation.Indeed, a connection between energy density contrast and curvature perturbation arises from the linearized Einstein equations [110,169,170]: The density contrast in Eq. (A4) is the density contrast in Eq. (A6) averaged over a Hubble sphere.This motivates to apply the results on PBH formation from curvature perturbations generated by inflation to the curvature generated during a 1stOPT. 1 Large curvature peaks collapse into PBHs if the averaged associated density contrast is larger than a given threshold value δ c [171].For collapse occurring during radiation domination era, the threshold has been found to take values within the range δ c ∈ [0.40, 0.67] depending on the profile shape [111,112].Specifically, δ c = 0.40 corresponds to the broadest profiles, while δ c = 0.67 corresponds to the sharpest profiles.Awaiting future studies to determine the shape of the density profile generated by 1stOPT, we assume δ c = 0.45, corresponding to a Mexican hat profile shape in radiation-domination.Such choice has been widely used in the literature [107][108][109][110].
In presence of bubble growth, the vacuum energy density decreases as: where F (t; t ni ) is the volume fraction of remaining false vacuum at time t [172]: where a(t) is the scale factor of the universe evolving with energy density ρ tot = ρ V + ρ R .Energy-conservation in an expanding universe implies: The critical time of nucleation postponing t PBH ni beyond which a causal patch reaches the threshold in Eq. (A4) can be found from numerically solving for ρ R (t; t ni ) in Eq. (A9) as a function of t ni for late-bloomers and at t ni → 0 for the background.The fraction P coll of causal patches which collapse into PBHs is given by the probability that causal patches remain vacuum-dominated until t PBH ni .Denoting by t max the time when the density contrast reaches its maximal value, set equal to the threshold in Eq. (A4), we have: where V (t ′ ; t max ) is the comoving volume at time t ′ which is in causal contact with the collapsing patch at t max : with r H ≡ (aH) −1 the comoving horizon and r(t; t ′ ) = t t ′ d t/a( t) the light travel distance between t ′ and t.In Ref. [21], a ready-to-use formula fitted on the numerical calculation is proposed, valid in the supercooled limit T n ≪ T eq : with a ≃ 0.5646, b ≃ 1.266, c ≃ 0.6639 and δ c ≃ 0.45.Since the collapse occurs during radiation, the PBH mass is given by the mass inside the sound horizon c s H −1 at the time of the collapse, e.g.[173]: where c s = 1/ √ 3 was used and where M d ≃ 3.7 × 10 −8 M ⊙ is the moon mass.We refer to [27].Counting the number of collapsing patches in our past light-cone and weighting them by M PBH , we obtain the PBH abundance in unit of the DM relic density: T eq 1 TeV . (A14) As we show in Fig. 4, values of β/H ≲ 7 are needed to produce PBHs in observable abundance f PBH ≳ 10 −5 .

B. Scale-invariant SM extension 1. SM only
In the SM, the electroweak scale can not be explained by radiative correction alone and a mass scale must be added at tree level.This is because the one-loop contribution to the Higgs potential reads [36]: Due to the large top mass, the constant β λh < 0 has the wrong sign and h = 0 is the only minima of the potential.

Effective potential
To generate the weak scale with a scale-invariant set-up, we complete the SM with a complex scalar Φ = ϕe φ/ √ 2 charged under U (1) D .The tree level potential in the unitarity gauge reads: The Renormalized Group Equations (RGE) of the theory are given by, e.g.[43]: Integrating the first equation, we obtain: The field excursion ϕ * during tunneling is of the order of the nucleation temperature ϕ * ∼ T n , see Eq. (B26) below.
The interval 1 ≲ T n /v ϕ ≲ 10 −10 characteristic of supercooled phase transition implies −23 ≲ t ≲ 0. Hence, we can safely neglect the running of g D : Due to collider bounds, the Higgs-mixing is small λ hϕ ≲ 10 −2 , so we can safely neglect its RGE running in Eq. (B4).The integration of Eq. (B5) gives: where λ 0 ϕ ≡ λ ϕ (v ϕ ) is a boundary value.The potential at 1-loop with counterterm δλϕ 4 /4 is: We impose the renormalization condition that ϕ = v ϕ is the minimum of the potential: which leads to: Plugging Eqs.(B8) and (B11) into Eq.(B9) and Taylor-expanding as a function of the logarithm t, we obtain: with: where the two right-hand sides assume g 2 D ≳ λ 0 ϕ , λ hϕ .We conclude that we can safely truncate the quantum loop correction to first-order in t (leading-log approximation) as long as: In this study, to minimize the number of parameters, we assume λ 0 ϕ = λ ϕ (v ϕ ) ≪ √ g D , which in practice leads us to set λ 0 ϕ = 0. Instead, despite its small value given in the main text, we keep λ hϕ finite in our equations.Computing the second derivative of the potential around its minimum, we get the scalar mass which is loopsuppressed with respect to the vector boson mass: The total potential is given by: The second piece contains the finite temperature correction at one loop and the Daisy contributions which accounts for the resummation of the Matsubara zero modes [53]: with vacuum and thermal masses m 2 V (ϕ) = g 2 D ϕ 2 and Π V = g 2 D T 2 /6 respectively [54].The third piece contains the Higgs contribution below the temperature of QCD confinement, see main text

Bounce action
The cost for nucleating a bubble through thermal tunneling is given by the action [62] where the tunneling trajectory reads We calculate the bounce action S 3 in Eq. (B20) both numerically and analytically.Numerical method -We wrote an over-shooting/under-shooting C++ algorithm CCBounce [51] using the full potential in Eq. (B17) and the quartic coupling being given by Eq. (B8).We checked that the error is below the percent level as explained in [65, chap. 6].Outputs of CCBounce are displayed in Fig. 5.The bounce action after slight smoothening is shown with thick lines in Fig. 6.
Thick-wall method -We numerically find that the field value ϕ * at the center of a nucleated bubble is of the order of the temperature ϕ * ∼ T n ≪ v ϕ .This allows two approximation schemes.First of all, we can Taylor-expand the potential V (ϕ) = V T =0 + V T + V QCD at first order in ϕ/T ≪ 1:   Second of all, we can use the thick-wall formula [64], see for instance the textbook [65, chap. 6]: In terms of the physical parameters, S 3 becomes: and the field value at the tunneling exit reads: Comparison between O 3 and O 4 -symmetric bounce trajectories -The decay rate of the false vacuum due to O(3)and O(4)-symmetric tunneling trajectories is given by [62,[174][175][176]: where A 3 = T 4 (S 3 /T /2π) 3/2 and A 4 = R −4 n (S 4 /2π) 2 are prefactors and R n is the bubble radius at nucleation.As shown in bottom-left panel of Fig. 5, the bubble radius at nucleation is about R n ≃ 10/T n , hence log where S crit is the critical bounce action (either S 3 /T or S 4 ) when nucleation becomes efficient.The O 4 -symmetric bounce for quantum tunneling reads in the thick-wall limit [64,65]: Comparing with Eq. (B24), we conclude that the O 4 tunneling rate is always smaller than the O 3 -symmetric tunneling rate when: Therefore, in the parameter space of interest tunneling is always dominated by the O 3 -symmetric contribution.Fig. 6 shows that the thick-wall formula compares reasonably well with the numerical formula.We can safely use it in the rest of this work.The nucleation temperature T n and phase transition rate β/H are plotted in Fig. 8, and also Fig. 1.
Validity of the β-parametrization of the tunneling rate -We now discuss the validity of truncating the Taylorexpansion of the tunneling rate logarithm at 1st order in t − t n in Eq. (A3).Assuming Γ V (t) ≃ T 4 e −S3/T , the Taylor-expansion at 2nd order reads This error is evaluated either at the nucleation delay time t PBH ni or at the time t max when the density contrast δ reaches its maximum value set to δ c ≃ 0.45, both of which are determined numerically.The black line represents the values of the PT completion rate β/H, while the brown region indicates the range that produces observable PBHs.
In the large supercooling limit we can set Ḣ = 0, and the ratio of the 2nd order term over the 1st order term becomes In Fig. 7, we show with blue lines the relative error ϵ β in Eq. (B33) evaluated at two characteristic times.First, we consider the time t PBH ni , which corresponds to the time until which nucleation is delayed within the past light cone of the PBH-forming region.Second, we consider the time t max , at which the density contrast reaches its maximum value, set equal to the threshold δ c ≃ 0.45.Both times are determined numerically following App. A. We conclude that neglecting the second-order correction in the Taylor expansion of the tunneling rate logarithm is a reasonable approximation when nucleation occurs through the Coleman-Weinberg potential alone, without QCD catalysis.Specifically, for values β/H ∼ 6 associated with efficient PBH production, the error ϵ β is only 5% at t PBH ni and 15% at t max .However, the error ϵ β can reach 50% at t PBH ni and 30% at t max when nucleation is catalyzed by QCD confinement.The relative error ϵ β is expected to increase at smaller values of v ϕ , where the effects of QCD are more pronounced.Given the uncertainties in modelling QCD catalysis, we consider this level of precision sufficient for the purposes of this work and leave a more detailed analysis for future studies.
For future reference, we note that the time t PBH ni can been estimated analytically from Ref. [21, App.C2] as where P coll is the collapse probability in Eq. (A12).We find that Eq. (B34) overestimate its numerically calculated value by less than 50%.
Validity of the single-field approximation -Apart from the QCD contribution, in the above discussion we have neglected the motion of the Higgs field h during tunneling of the U (1) D scalar ϕ.We now check that this is a valid approximation.In the supercooling limit T n ≪ T eq , the exit point of tunneling ϕ * in Eq. (B26) is hierarchically small: The motion of ϕ during tunneling does induce a motion of h if the negative induced mass − 1 2 λ hϕ ϕ 2 * is larger than   3 in the main text under the assumption that the QCD catalysis does not arise and that the universe follows an adiabatic evolution after the PT.The later point would arise if the dark scalar decays right after percolation, see e.g.[22].In contrast to the minimal case presented in the main text, all the regions producing observable PBHs are testable with future GW interometers, above astrophysical foregrounds, including the region explaining DM as made of PBHs.The thin orange region in the bottom-left corner can explain the six ultra-short lensing events observed with OGLE telescopes after five years of observation of the Galactic bulge and which could be attributed to PBHs with masses 10 −6 M ⊙ > M PBH > 10 −3 M ⊙ [84,177].This region is absent in Fig. 3, which implies that effects from QCD catalysis prevent the scale-invariant EWPT to explain OGLE events.In principle, OGLE events could be explained if QCD corrections to the EW phase transition are weakened for some reasons.
the positive thermal mass m 2 h (T ) ≃ ( 3g where we used λ hϕ = 2λ h (v EW /v ϕ ) 2 and λ h ≃ 0.13.We conclude that we can safely neglect the motion of the Higgs during tunneling (except for QCD catalysis discussed in the main text).

Collider constraints
Mixing angle -Within a view to applying collider bounds to the SM extension we consider, we review the derivation of the mixing angle θ hϕ .The Lagrangian after spontaneous symmetry breaking (with quantum corrections included) reads: with and Φ = vϕ+ φ √ 2 the gauge eigenstates.The minimum of the potential at (v ϕ , v EW ) implies: The mixing angle θ hϕ is defined as the rotation angle between the gauge eigenstates and the mass eigenstates: with associated mass eigenvalues: We find: where where m ϕ = β λ v ϕ is the scalar mass and H 2 n = ∆V /3M 2 pl is the Hubble rate at percolation.The other decay channel, through Higgs mixing, is further suppressed by (v EW /v ϕ ) 2 : where Γ h ≃ 4.1 MeV is Higgs decay width [178].
Higgs searches -We use Eq.(B43) to recast existing collider constraints on additional Higgs boson on the parameter space of the minimal scale-invariant U (1) D extension of the SM, see Fig. 10.For g D ∼ 0.5, we find that the tightest constraints are given by LEP searches v ϕ ≳ 1.5 TeV, the mass of the dark Higgs being m ϕ ≳ 80 GeV., resonant production at LHC [104] (green), electroweak fit of the W mass [104] (red), Higgs signal shape [105] (purple), Higgs decay to new particle [106] (gray).In the brown region, the scalar mass m ϕ is induced by the Higgs VEV.

C. Bubble dynamics 1. Energy budget
The motion of bubble walls is subject to the driving vacuum pressure ∆V and to the retarding friction pressure induced by the bremsstrahlung emission of vector bosons V by scalar particles ϕ entering bubbles.After resumming leading-logs quantum corrections (LL), the friction pressure reads [66][67][68][69]: where γ is the bubble wall Lorentz factor.We suppose the leading-order pressure P LO to be sub-dominant since it does not grow with γ [179].In the absence of friction P LL ≪ ∆V , all the vacuum energy E vac = 4πR 3 ∆V /3 of a bubble of radius R is converted into kinetic energy of the wall E wall = 4πR 2 σγ, where σ is the wall surface tension.
In that case, the wall Lorentz factor at collision reaches the maximal value allowed by energy conservation: The mean bubble radius at collision is given by R coll ≃ π 1/3 /β [180].We introduced R crit = 2σ/∆V which is the critical radius for nucleating a bubble, found as the saddle point of the free energy F (R) = 4πσR 2 − 4πR 3 ∆V /3 in the thin-wall approximation.As bubble walls accelerate, the retarding pressure P LL in Eq. (C1) grows linearly with γ.The walls stops accelerating as soon as P LL ≃ ∆V , with the associated Lorentz factor: Whether the latent heat of the PT is dominantly converted into the plasma or into the kinetic energy of the wall depends on whether Eq. (C2) or (C3) dominate.The Lorentz factor at bubble collision is given by the formula: The fraction of the latent heat converted into wall kinetic energy reads (see also [181]):

Reheating after a matter era
At the time of bubble percolation, the latent heat of the universe has been converted into two fluids, the radiationlike energy density ρ shock of the ultra-relativistic shock waves generated by the work of the friction pressure, and the radiation-like energy density of the scalar field gradient ρ ϕ .Both are highly peaked distribution of energy-momentum tensor.After that walls pass through each other, the highly peaked energy stored in the wall E wall ≃ γ coll σR 2  coll is converted into a broadly-distributed oscillating scalar field condensate with energy E osc ≃ ∆V R 2 coll ∆R where ∆R is the distance between the two peaks in energy-momentum tensor after their collision.The radiation-like energy stored in the collided walls is converted into the matter-like oscillating scalar condensate E wall = E osc after the time: If the scalar field is long-lived Γ ϕ ≪ H, then the scalar field energy density starts redshifting like matter with an energy density after percolation given by: The radiation energy density reads: The scalar field decays after the time approximately given by its inverse decay width: where we also gave the associated temperature.If long-lived enough, this would lead to a matter-domination era starting at the time and temperature: obtained from equating Eqs.(C7) and (C8) with a ∝ t 1/2 .We conclude that if t dom < t dec then the PT is followed by a matter era starting at t dom and ending at t dec .A transition from a matter to radiation inject an amount of entropy in the universe, e.g.[182][183][184], whose magnitude is given by the dilution factor: where S i and S f is the total entropy in the universe just before and just after the decay of the scalar field, assuming an instantaneous decay.The second equality comes from evolving matter and radiation from T dom to T dec .The entropy injection dilutes the abundance of any relic present in the universe before t dom .This is the case of PBHs and GWs.For PBHs, the abundance is simply rescaled by: For GWs, the spectrum is rescaled by: where the D factors trivially correct for the redshift of the peak amplitude and frequency [130].Instead, the effects behind S M (f ) are more subtle and consist to change the spectral slope to ∝ f 1 [126][127][128][129] instead of f 3 [185][186][187] for modes with wavenumber k = 2πf which are super-horizon k < H n /a at percolation and enter the causal horizon k = H/a during the matter era.In other words, S M (f ) is designed to change the spectral slope to: assuming continuity with the standard spectrum at f dom .We have introduced the frequencies corresponding to the causal horizon H/2π at the beginning and end of the matter era red-shifted up to today: with H ≃ 0.3g 1/2 * T 2 /M pl given by Friedman's equation and: T dom 100 GeV g * (T dom ) 100 where the dilution factor D is only applied to f dom .The GW spectrum Ω GW (f ) from bubble dynamics is presented in the next section.

GW signal
Scalar field gradient -A fraction κ coll given by Eq. (C5) of the latent heat is stored in the scalar field gradient localised in bubble walls.Initially, the GW spectrum has been calculated in the "envelop" approximation where walls are infinitely thin and collided parts are removed from the calculation [192][193][194][195][196]. The collided part were studied later both analytically [121] and numerically [122][123][124][125] in what is known as the "bulk flow" model.It is found that the collided part continue to source GW after collision, generating the slow decreasing IR slope Ω GW ∝ f 1 .The GW spectrum in the bulk flow model gives (v w = 1) [122]: Figure 11 : Gravitational waves spectrum produced by the supercooled first-order phase transition predicted in the minimal U (1) D extension of the SM.Solid lines shows the GW spectra accounting for the matter-dominated period after the PT, which suppresses the signal at large dark Higgs VEV v ϕ ≳ 10 4 TeV.In comparison, dotted lines show the same GW spectra without the matter era.The colorful regions show the existing constraints from Earth-based laser interometers LIGO-Virgo-Kagra (LVK) run O3 [70], assuming Signal-to-Noise Ratio (SNR) detection thresholds of 2, and prospective constraints from LIGO run O 5 [71], their follow-ups Einstein Telescope (ET) [72,73] and Cosmic Explorer (CE) [74], space-based laser interferometer LISA [75][76][77], Earth-based atom interferometer AION [188,189] in its space version AEDGE [190].For LISA, ET and CE, we optimistically assume a low SNR threshold of 1 (equivalent to a 1σ deviation from noise) after 20 years of observation [191].The gray regions show the expected astrophysical foreground from galactic white dwarf binaries (WD B.) according to model 1 [79,80] and 2 [76], and extragalactic compact (neutron stars and black holes) binaries (Compact B.) fitted on LIGO O3 data [70].Instead, the contribution from extragalactic supermassive black holes binaries lies at lower frequencies [81].All the GW constraints in Fig. 9 and Figs of the main text are derived assuming that the GW spectra is above the astrophysical foregrounds inside the frequency window of experiments.
with the redshift factor between percolation "n" and today "0": T eq 100 GeV g * (T eq ) 100 and the spectral shape S(f ) peaking on f p : S(f ) = 3(f /f p ) 0.9 2.1 + 0.9(f /f p ) 3 , f p = a n a 0 0.8 β 2π . (C21) The dilution factor D ≥ 1 given by Eq. (C11) accounts for the additional redshift due to entropy injection in presence of an eventual early matter era if the scalar ϕ is long-lived.We added the correction factor: with c * = O(1) to impose the scaling Ω GW ∝ f 3 for emitted frequencies smaller than the Hubble factor H * /(2π) and entering during radiation [129,[185][186][187].We set c * = 1 and defer the determination of c * for future studies, see e.g.[197,198].The factor S M (f ) corrects the causality tail f 3 → f 1 for modes entering during the matter era instead of the radiation era, as explained in Eq. (C14).Ultrarelativistic shells -We know discuss GW generation from the fraction 1 − κ coll of latent heat stored in terms of fluid motions.In the ultra-relativistic limit, one expects bubble walls to be followed by ultra relativistic shells of particles which are generated by plasma/wall interactions [67-69, 114-116, 118] which are extremely thin [114] and whose impact on the GW spectrum is not yet understood [119,120] at the time of writing.Waiting for future studies to investigate the precise GW spectrum from ultrarelativistic shells, we posit that from the point of view of gravity, the GW production from extremely thin shells of scalar field gradient should be indistinguishable from extremely thin shells of ultra-relativistic particles. 2A difference could arise however due to ultra-relativistic shocks being potentially long-lived [118] while scalar field gradient having a reduced lifetime given by Eq. (C6).One expect such difference to enhance the GW spectrum by a factor β/H. Since this is a small factor in the regime of interest for PBH production, we do not include it in our plot.However, we checked that the conclusion of this paper do not depend on this detail.Additional enhancement due to turbulence [199][200][201][202][203] and second order GWs [204,205] are foreseeable and left for future works.
Hence, we suppose that the bulk flow model holds in all the parameter space.However, the friction still has an important impact on the GW spectrum through the dependence of the dilution factor on the latent heat fraction into scalar field gradient, see Eqs. (C11), (C10) and (C5).We show a family of GW spectrum for different values of PT scale v ϕ in Fig. 11 together with LISA reach and astrophysical foregrounds.The constraints on the scale invariant U (1) D extension of the SM are shown in Figs. 3 and 9.For large dark Higgs VEV v ϕ ≳ 10 4 TeV, the extended dark Higgs lifetime generates entropy injection and suppresses all GW signals, leaving PBHs as the only possible signatures.

Figure 1 .
Figure 1.The PT completion rate β/H becomes lower as the amount of supercooling increases, until QCD effects kick in around Tn ≃ TQCD in Eq. (13).The lines are obtained from varying gD.The parameter α along the horizontal axis is the latent heat fraction.

Figure 2 .
Figure 2.  In an adiabatic universe, supercooled phase transition produce gravitational waves within the frequency window of future interferometers LVK run O3[70], run O5[71], ET/CE[72- 74]  and LISA[75][76][77][78] with an amplitude larger than astrophysical foregrounds from white dwarfs binaries and compact binaries[70,76,[79][80][81].Instead, in the minimal scale-invariant extension of the SM, the long lifetime of the scalar field dilute all the GW signals at large dark Higgs VEV v ϕ ≳ 10 4 TeV.Instead the PBH abundance is only slighly impacted by the matter era.See App.C for more details on the experiments reach and expected astrophysical foregrounds.

5 Figure 3 .
Figure 3. PBH and GW predictions in the minimal scale-invariant U (1)D extension of the SM.The model has only two parameters: the hidden U (1)D breaking scale v ϕ and its gauge coupling gD.The purple region shows the microlensing constraints (ML) from HSC/subaru telecope[82,83].The blue U-shaped region can explain the HSC event[84][85][86].The boundary of the brown region can explain dark matter.Hawking evaporation overproduce cosmic rays in the red region[87][88][89][90][91][92][93][94][95], but could contribute to the 511 keV excess in the green spot[87][88][89].Non-observation of energy injection from Hawking radiation in CMB[96][97][98] and BBN[90,99,100] exclude the yellow and green regions.See[8] for a review of PBH constraints.The GW constraints shown in light blue disappear at large v ϕ ≳ 10 4 TeV due to entropy injection following the decay of the long-lived dark Higgs.The associated dilution factor D is indicated with blue lines labelled "D".Part of the PBH DM region can not be probed with GWs.On the left of the dashed black line, the PT is catalysed by QCD confinement.In the blue region in the top-left corner, the universe can not tunnel without "QCD effects", due to ΓV = H 4 having no solution Tn.However, the QCD catalysis temperature given in Eq. (13) is smaller than the Hubble scale TQCD ≲ Hn, implying that the tunneling dynamics receives large correction from "gravity effects"[101,102], which we leave for future research.The gray region displays collider constraints on additional Higgs boson[103][104][105][106], see Fig.10in App.B for more details.On the right of the very thin blue line, bubble wall friction convert most of the the latent heat into the plasma, which reduces the energy fraction stored in the scalar field and therefore reduces the dilution factor D.

Figure 4 :
Figure 4 : PBH abundance as a function of the amount of supercooling T n /T eq and rate of completion β/H.

Figure 5 :
Figure 5 : Outputs of CCBounce [51]: the scalar potential (top), the bounce profile at nucleation (bottom left) and the bounce action for different values of the temperature (bottom right).We have fixed the U (1) D gauge coupling g D = 0.45 and the U (1) D scalar vev v ϕ = 17 TeV.

Figure 6 :
Figure 6 : Bounce action S 3 /T (left) and rate β/H of the phase transition (right).The thickwall analytical formula is a good approximation of the numerical results.

Figure 7 :
Figure 7 : The blue lines depict the relative error ϵ β = β 2 (t − t n )/2β in Eq. (B33) resulting from neglecting the second-order correction in the Taylor expansion of the tunneling rate logarithm ln ΓV (t) ≃ β(t − t n ) + β 2 (t − t n ) 2 /2.This error is evaluated either at the nucleation delay time t PBH

Figure 8 :
Figure 8 : Parameter space of the minimal U (1) D scale-invariant model with large amount of supercooling T n /T eq < 1 (purple lines) and slow completion rate encoded by β/H (orange lines).See Fig. 3 for the rest of the legend.

Figure 9 :
Figure9: Same as Fig.3in the main text under the assumption that the QCD catalysis does not arise and that the universe follows an adiabatic evolution after the PT.The later point would arise if the dark scalar decays right after percolation, see e.g.[22].In contrast to the minimal case presented in the main text, all the regions producing observable PBHs are testable with future GW interometers, above astrophysical foregrounds, including the region explaining DM as made of PBHs.The thin orange region in the bottom-left corner can explain the six ultra-short lensing events observed with OGLE telescopes after five years of observation of the Galactic bulge and which could be attributed to PBHs with masses 10 −6 M ⊙ > M PBH > 10 −3 M ⊙[84,177].This region is absent in Fig.3, which implies that effects from QCD catalysis prevent the scale-invariant EWPT to explain OGLE events.In principle, OGLE events could be explained if QCD corrections to the EW phase transition are weakened for some reasons.

Figure 10 :
Figure10: Collider constraints on mixing angle θ hϕ in scale-invariant U (1) D extension of the SM (black line) as a function of the extra Higgs mass m ϕ .We show LEP searches for 4-fermion final states from Higgs strahlung e + e − → Zϕ[103] (blue), resonant production at LHC[104] (green), electroweak fit of the W mass[104] (red), Higgs signal shape[105] (purple), Higgs decay to new particle[106] (gray).In the brown region, the scalar mass m ϕ is induced by the Higgs VEV.
We conclude that the PT rate β/H is minimized for coupling g D associated with the longest supercooling stage while not being catalyzed by QCD Λ QCD ≲ T n ≪ T eq .The nucleation temperature T n and the PT rate β/H n in scale-invariant SM+U (1) D are shown in Fig. 1.

Table of contents
Whenever Eq. (B36) is verified, during tunneling the Higgs field acquires the value h exit ≃ ϕ * λ hϕ /λ h .The motion of h backreacts on the tunneling if the negative induced mass − 1 2 λ hϕ h 2 exit is comparable to the thermal mass 2EW and λ ϕ = 11β λ /6.The dominant dark Higgs decay channel is into two Higgs ϕ → hh: