Gravitational waves from supercooled phase transitions: dimensional transmutation meets dimensional reduction

Models with radiative symmetry breaking typically feature strongly supercooled first-order phase transitions, which result in an observable stochastic gravitational wave background. In this work, we analyse the role of higher order thermal corrections for these transitions, applying high-temperature dimensional reduction to a theory with dimensional transmutation. In particular, we study to what extent high-temperature effective field theories (3D EFT) can be used. We find that despite significant supercooling down from the critical temperature, the high-temperature expansion for the bubble nucleation rate can be applied using the 3D EFT framework, and we point out challenges in the EFT description. We compare our findings to previous studies, and find that the next-to-leading order corrections obtained in this work have a significant effect on the predictions for GW observables, motivating a further exploration of higher order thermal effects.


Introduction
The observation of a stochastic gravitational wave (GW) background from a primordial first-order phase transition would unravel information about underlying particle physics beyond that of the Standard Model (SM).A very interesting beyond-the Standard-Model (BSM) scenario is the case of a supercooled first-order phase transition, which typically arises in models with classical scale invariance (or nearly conformal dynamics) [1][2][3][4][5].In such a case, the phase transition completes at a temperature much below the critical temperature.As a result, the potential energy difference between the high-temperature and low-temperature phases becomes very large, and the amount of energy released -relative to the radiation energy density -is orders of magnitude larger than in scenarios without significant supercooling.Large energy release results in a strong GW signal sourced by the sound waves in the plasma or the collisions of the bubble walls [6][7][8][9][10][11]. Predictions of the GW spectrum for models with classical scale invariance [4,5,[11][12][13][14][15][16][17][18][19][20][21][22][23] indicate that the signal could be readily observed by the Laser Interferometer Space Antenna (LISA) [24] and other next-generation GW detectors [25][26][27].This makes models with classical scale invariance and strong supercooling an interesting theoretical playground, and accurate predictions of the GW spectrum in terms of the free parameters of such models are essential to determine if a potentially observed GW signal was caused by a phase transition in such a model.
Predicting the GW signal requires a determination of thermal parameters describing the phase transition, such as the percolation temperature T p , the strength α, the (inverse) time or length scale of the transition, β or R * and the wall velocity v w .In many studies, the phase transition parameters are obtained from the one-loop effective potential at finite temperature, with so-called daisy resummation accounting for a resummation of a class of diagrams enhanced in the infrared (IR) due to thermal screening.In recent years, it has become clear in the context of Higgs portal models [28] and the Standard Model Effective Field Theory [29] that this approach might not predict the thermal parameters with sufficient precision, and the corresponding uncertainty of the GW signal can be several orders of magnitude.In the work at hand, we apply similar higher-order thermal corrections to models with classical scale invariance.
The reason for the poor convergence of the computation at finite temperature, is that bosonic low-energy modes become highly occupied in a thermal plasma.This results in a breakdown of the usual loop expansion [30][31][32].Indeed, the standard one-loop procedure suffers from an incomplete treatment of the perturbative expansion, which reveals itself as an uncancelled dependence on the renormalisation scale [28,32].The escape out of this distress is the use of a dimensionally reduced effective field theory (EFT) [33][34][35][36], that is constructed to account for thermal scale hierarchies and consistently incorporates the required thermal resummations, which significantly reduces the uncertainty of the GW signal predictions [28,29] (c.f. also ref. [37]).This method allows one to construct an EFT for only the degrees of freedom that are driving the phase transition at IR length scales.The heavy ultraviolet (UV) modes are integrated out, and their effect is captured in the parameters of the EFT via matching.Constructing the EFT can be a technically challenging endeavour compared to the use of mere one-loop thermal functions with minimal daisy resummation that encode the leading behaviour of the effective potential, but this obstacle has been largely removed by DRalgo [38], which has automated the matching procedure and the computation of the effective potential in the EFT for generic models.Furthermore, the formulation in terms of an effective field theory combined with strict perturbative expansions has been shown to provide a theoretically sound setup for computations, that is free of residual gauge dependence, imaginary parts, spurious IR-divergences or double counting contributions [39][40][41][42][43][44][45].Indeed, in the terminology used in ref. [45] we implement the mixed method in the computation of the bubble nucleation rate, which is based on the strict expansion for the action around the leading order bounce solution.
So far the dimensionally reduced EFT approach has not been applied to models with classical scale invariance.At first glance, the approach might not even seem suitable for the study of supercooled phase transitions, as the construction of the dimensionally reduced EFT relies on scale hierarchies in a high-temperature (HT) expansion, assuming that the field-dependent masses are small compared to the temperature.This assumption seems not at all appropriate for a phase transition in a scale-invariant model: the position of the minimum of the potential of the transitioning field exceeds the temperature by multiple orders of magnitude.This suggests that applying the dimensionally reduced EFT to models with classical scale invariance might do more harm than good: does the inclusion of higherorder corrections in the effective potential come at the cost of applying the HT expansion in a regime where it is not at all valid?
In this work, we will argue that the EFT relying on the HT expansion can be used for parts of the computation.The crux is that the transitioning field does not transition directly to the minimum of the potential, but remains in the regime of validity of the HT approximation.Therefore, along the path formulated in refs.[39,[43][44][45][46], we compute the thermal contributions coming from the so-called hard scale (c.f.section 4) and construct an EFT for the bubble nucleation at the soft scale.This EFT can be used for the determination of the nucleation and percolation temperature, and the typical length scale of the transition.Other parameters, such as the phase transition strength, do depend on the value of the potential at its minimum.These quantities have to be determined without the hightemperature expansion, but follow from the zero-or low-temperature potential (see e.g.ref. [11] for the details on how to compute the reheating temperature, or the potential energy difference ∆V ).For concreteness, we will demonstrate the approach explicitly in the SU(2)cSM model [4,47], a conformal extension of the SM.
We find that the next-to-leading (NLO) corrections included in the EFT modify the predictions for the properties of the phase transition significantly, as compared to earlier results based on daisy resummation [11].For example, the percolation temperature can change by 100%, whereas the changes in the length scale, given by the normalised bubble radius R * H * , reach 50%.Since the signal is expected to be well visible with LISA, it would be possible to reconstruct the values of R * H * and the reheating temperature with good accuracy [48].This clearly shows the importance of providing the most precise theoretical predictions possible.
Interestingly, the modification of the potential at NLO accounts only for a part of the large correction described above.A correction of the kinetic term in the action, only appearing at NLO, is responsible for a significant shift in the results.This kind of correction is not straightforward to include within the conventional daisy-resummed approach, which shows the importance of using the EFT framework.On the other hand, this correction is also a main source of uncertainty in our computation, and it could indicate the breakdown of the mass hierarchies at the root of the applied EFT.These non-trivial issues in the construction of the EFT for classically conformal theories set the stage for further studies.
This article is organised as follows.In section 2 we review the previous knowledge on phase transitions in models with classical conformal symmetry and discuss the applicability of high-and low-temperature approximations.In section 3 we introduce our concrete BSM model, the SU(2)cSM, for which in section 4 we construct an effective description for bubble nucleation at high temperature, using an EFT at the soft scale.In section 5 we present our numerical results and we summarise our findings in section 6.For the convenience of the reader, we provide the expressions for the running couplings in appendix A and our implementation of dimensional reduction using DRalgo in appendix B.

Supercooling at high temperature
The equilibrium properties of a high-temperature plasma can be described in Matsubara's imaginary time formalism [49].In this formalism, fluctuations of several mass scales arise: modes with non-zero Matsubara frequency have masses of the order πT or higher at temperature T , while zero modes can have masses which are parametrically smaller m ∼ gT [50,51].Here, the dimensionless coupling g 2 (4π) 2 ≪ 1 parametrises the hierarchy.Such a hierarchy allows for a HT expansion with respect to m T ∼ g, and suggests an EFT picture, where an EFT for long-distance IR physics for phase transitions is constructed, by integrating out short-distance non-zero modes in the UV.These UV modes screen the modes in the IR and generate thermal mass corrections [52].Capturing such effects requires resummation of perturbative expansions.The scalar field zero modes undergo the phase transition, and since they are static and live in three spatial dimensions, this procedure is called high-temperature dimensional reduction.We describe this in more detail in section 4.
In models with classical scale invariance, all fields are massless at classical level, and massive modes are generated radiatively at loop level by quantum corrections.This is called dimensional transmutation [53].Physical masses depend on the vacuum expectation value of the scalar field which is typically much larger than the nucleation temperature, thus m T ≫ π, seemingly invalidating the use of HT expansion.Hence, at first sight, supercooling seems incompatible with the formalism of dimensional reduction as the latter relies on the HT expansion.In this section, we delve into this seeming contradiction to formulate a consistent prescription for treating supercooled phase transitions with due accuracy.For the current purpose, we will phrase our discussion in terms of the one-loop effective potential with Arnold-Espinosa -or daisyresummation [32]: a framework which is familiar to most readers and corresponds exactly to dimensional reduction at leading order (see section 4.2), where thermal corrections to the masses are computed at one-loop, and resummed to all orders.We start with a brief review of the temperature-dependent effective potential and the parameters characterising phase transitions.

Perturbative description of a phase transition
In perturbation theory, the effective potential is the central object for computing the properties of the phase transition.It is given by a loop expansion as where V (0) corresponds to the tree-level potential, V (1) is the one-loop Coleman-Weinberg correction, including counterterms to remove divergences, and V T contains thermal corrections, at one loop.Herein, two-loop corrections are not considered.In this work, we will use Landau gauge and the MS renormalisation scheme.
The one-loop zero-temperature correction is given by the well-known formula [53] being a sum of contributions from different fields with n a counting the number of degrees of freedom as where s a denotes the spin of a given particle, Q a = 1(2) for neutral (charged) particles, and N a = 1(3) for uncoloured (coloured) particles.C a = 3/2 for scalars and fermions, and C s = 5/6 for vector bosons.Here, for simplicity, we assume the effective potential to be a function of a single field but it can be straightforwardly generalised to the multi-field case.
The one-loop thermal correction is given by where the thermal function is defined as with the "+" sign for fermions and "−" sign for bosons [51].When the temperature is high with respect to the mass scale, M a /T ≪ 1, the thermal functions can be expanded as and in the opposite regime M/T ≫ 1, the expansion is which is the same for bosons and fermions.Above, γ E denotes the Euler-Mascheroni constant, ζ is the Riemann zeta function and Γ the gamma function.
At high temperature, diagrams that involve IR-sensitive zero modes get enhanced due to thermal screening, and these diagrams need to be resummed.For now, we will employ the Arnold-Espinosa daisy resummation method [32], and further sections will be devoted to a detailed discussion of applying the dimensional reduction scheme.The one-loop potential with daisy resummation reads where and M 2 i,th (φ, T ) denotes the thermally corrected mass squared, which is the sum of the squares of the zero-temperature and the thermal mass.Diagrammatically, the first term here is a sum of an infinite number of diagrams of a one-loop zero-mode diagram with nonzero-mode one-loop diagrams attached around it ✿, hence the name daisy.The physics behind this construction is clear: the non-zero UV modes screen the zero mode living at the IR scale, and the first term in eq.(2.9) is nothing but a result of a one-loop diagram of the zero-mode, with resummed mass.The second term merely removes the double counting, since the zero mode contribution with unresummed mass is already included in the cubic term of eq.(2.6).We emphasise that eq.(2.9) should only be added to the effective potential whenever the HT expansion is valid, as it is an essential assumption in its derivation.Using the 3D EFT approach, it is straightforward to derive eq.(2.9), and we compute it explicitly in section 4.
The temperature evolution of the effective potential determines the details of the phase transition.A supercooled phase transition typically proceeds as follows.At high temperature, the scalar field fluctuates around the symmetric minimum, and as the temperature decreases, another minimum is formed.At the critical temperature T c the two minima become degenerate and at lower temperatures the symmetry-breaking minimum becomes energetically favourable.It is characteristic of supercooling that the transition does not proceed right after it has become energetically favourable.First, at temperature T V , the Universe enters a stage of thermal inflation induced by the large amount of energy stored in the false vacuum.Then, at some point, the field transitions to the true ground state, due to getting kicked by thermal fluctuations (typically quantum tunnelling is much less probable [11,16]).The nucleation temperature, T n , at which at least one bubble of the true vacuum is nucleated per Hubble volume is considered as the onset of the transition.Later the bubbles percolate at the percolation temperature, T p .To consider the phase transition complete, not only the fraction of the volume turned into the true vacuum has to be big enough, but also the volume of the false vacuum should be shrinking at T p [54], see also ref. [55].This condition is not trivially satisfied for transitions taking place during a phase of thermal inflation and thus it constrains the available parameter space.The size of the bubbles at the moment of collision, R * , can be used to estimate the length scale of the transition.It can be used interchangeably with the (inverse) time scale of the transition, β * , given by the derivative of the decay rate of the false vacuum.During the phase transition latent heat is released, which in the case of supercooling is tightly related to the strength of the transition α given by where ρ rad is the energy density stored in radiation and ∆V is the potential energy difference between the false and true minima.The released energy is partially converted to gravitational waves.In the case of strong supercooling, two production mechanisms -via bubble collisions and sound waves in the plasma -can be effective (see e.g.refs.[7,11]).For the calculation of the terminal Lorentz factor of the bubble wall one needs to consider the pressure difference across the wall, and for the NLO pressure contribution we use γ-scaling [56][57][58] (for further details, see [11]).We will use the efficiency factor for production via sound waves, κ sw to determine the dominant source, as the efficiency for production via bubble collisions is given by κ col = 1 − κ sw .However, most of the energy goes back to the plasma, reheating it back to T V . 1 The reheating temperature and the length/time scale of the transition (evaluated at T p ) are the parameters most relevant for the determination of the resulting GW spectra.Since the GW spectrum depends on α/(α + 1) and α ≫ 1 for strongly supercooled phase transitions, the exact value of α becomes irrelevant.More detailed definitions of the relevant parameters listed above can be found in ref. [11], the approach of which we follow here.
The decay of the false vacuum is controlled by the rate given by2 where A is T -dependent pre-factor.The three-dimensional Euclidean action is evaluated at the so-called bounce configuration which corresponds to the solution of the bounce equation with boundary conditions dφ b dr = 0 for r = 0 and φ b → 0 for r → ∞.

High-and low-temperature regimes
The ratios of the field-dependent masses to the temperature determine whether the HT or LT limit should be considered.Large field values correspond to large masses and thus LT, while small field values correspond to the HT limit.In models with classical scale invariance, which feature supercooled phase transitions, the scales associated with the global minimum of the potential (the location of which determines the strength of the transition) and with the location of the barrier (where the tunnelling takes place) are widely spread.Therefore, we cannot use just one of the limits, either LT or HT, to have the full picture of the transition.
In classically scale-invariant models, the potential around the global minimum, for temperatures below the critical temperature, is in the low-temperature regime.This means that we can use the one-loop thermally corrected potential, without daisy resummation to compute the temperature at which thermal inflation starts, T V , and to compute the vacuum energy close to the nucleation temperature we can even neglect the thermal corrections.
In the presence of various energy scales, we should use the renormalisation group (RG) improved effective potential to resum the field-dependent logarithmic terms and make the potential perturbative over a wide range of field values.For theories in which the one-loop corrections to the potential are dominated by a single mass scale M (φ) this is straightforward to attain by the field-dependent choice of RG scale, µ = M (φ).Accordingly, all couplings are run to this scale.However, going from the high-field regime to lower field values, which are relevant for the tunnelling, at some point the ratio M (φ)/T becomes small, which signals the onset of the high-temperature regime.
In the high-temperature regime we can use the expansion of eq.(2.6) in the potential but also resummations of higher-order terms are obligatory.The one-loop potential in the high-temperature limit reads (for simplicity we consider here models with bosons only) Note that the dependence on the logarithm of the field-dependent mass cancelled out. 3ow the only logarithm present is of the ratio of the temperature and the renormalisation scale.It is thus clear, that in order to preserve perturbativity of the computations one should fix the renormalisation scale to be proportional to the temperature, with some O(1) proportionality factor.The most natural choice is µ = 4πe −γ E T , which cancels the logarithmic term entirely.Nonetheless, any choice of µ ∼ κT , where roughly κ ∈ (1, 2π) is acceptable.
This choice builds a bridge between the HT and LT regimes.In the HT regime we have µ = κT , whereas in the LT µ = M (φ).In the intermediate regime, we thus have µ ≈ κT ≈ M (φ): deep in the HT regime we should have M (φ)/T ≪ 1 and we expect that the breakdown of the applicability of the HT approximation occurs for M (φ)/T ∼ O(1).These observations teach us how to treat the RG-improved potential for the sake of phase transition-related computations: at large field-values the scale should follow the field, whereas at low field-values the scale should be set by the temperature as4 µ = max(M (φ), κT ). (2.15) The thermal cutoff on the running of the couplings prevents them from reaching the Landau poles of e.g. the top Yukawa coupling at small field values so it regulates the behaviour of the potential around φ = 0.One should note, however, that the HT effective potential of eq.(2.14) is not renormalisationscale independent, in contrast to the zero-temperature effective potential [28,32].While the running of the parameters in V (0) cancels the explicit logarithms appearing in V (1) , there is still RG-scale dependence leftover.The implicit running of the term that corresponds to the thermal correction to the mass at one-loop, is of the same order as the running of V (0) , yet there are no T 2 -dependent logarithms of the form of in the potential for compensation.Indeed, at HT such terms only appear at two-loop order, due to the relatively slower convergence of perturbation theory, induced by enhancement due to thermal screening.The running of the other thermal contributions, in particular from the daisy resummation, is of higher order.In perturbation theory, renormalisationscale dependence can be used to estimate the size of the missing corrections, 5 and omission of such corrections is the source of one of the largest uncertainties in predictions of GW signals originating from phase transitions [29].Dimensional reduction will be the technique advocated in this work to be used in the HT regime to include these missing large corrections.In section 4 we discuss how these missing corrections are included and further resummations are performed with care.In summary, due to different properties of the theory in the UV (large field values, low temperature) and in the IR (small field values, high temperature), the field space naturally divides into two parts: 1.The low-temperature (LT) regime, where no resummations in the thermal part of the potential are needed.The RG-improved potential with running couplings and fields should be used and the renormalisation scale should follow the value of the field.
2. The high-temperature (HT) regime, where the scale at which computations are performed is set by the temperature.This is the region where thermal resummations are inevitable and the dimensionally reduced theory can be used to cancel renormalisation-scale dependence and systematically include higher-order corrections.
In the case of scale-invariant potentials, which we treat as models for supercooled transitions, there is also a natural division of parameters which need to be computed in different field-or scale-regimes: 1. LT regime: the location of the symmetry-breaking minimum, the energy of the true vacuum (needed to determine the phase transition strength α), the temperature at which thermal inflation starts T V , the reheating temperature, T r .
2. HT regime: the percolation and nucleation temperatures, the size of bubbles at collision R * H * , the inverse time scale of the transition β/H * (normalised to the Hubble length/time).
Therefore, one can compute the quantities in the first category using a finite-temperature RG-improved potential without resummations, while for the quantities in the second category, which are related to solving the bounce equation, the use of the dimensionally reduced theory is required.Below, when we consider a concrete model, we will demonstrate that indeed the escape point for the bounce trajectory lies within the HT region.In the literature the HT approximation has been used for computations related to the phase transition in supercooled transitions, see e.g.refs.[60][61][62].However, using dimensional reduction in the HT regime and relating it to the low-temperature limit through the RG running and using both of them is a novel approach which we present in this work.
We will demonstrate how this construction works by applying it to a model with classical scale invariance and an extra SU(2) X gauge group in the proceeding sections.

Introducing the model
In this section, we will apply the discussion of the previous section to a concrete BSM model, the so-called SU(2)cSM model [4,47].It is an extension of the conformal version of the SM (without the explicit mass term for the Higgs field) with a new, "dark" SU(2) X gauge group and a scalar that is a doublet under this new symmetry, while transforming as a singlet under the SM gauge group.This model has been studied extensively in the literature [4,11,16,17,47,[63][64][65][66][67][68][69], in particular in ref. [11] a thorough analysis of the thermal history of the Universe within this model has been performed using RG improvement and daisy resummations.In this work, our aim is to improve these results by including higher-order thermal corrections obtained with the dimensionally reduced theory.

Model at zero temperature
The model contains two complex scalar doublet fields.We exploit the symmetries to rotate the fields such that the vacuum expectation values are only non-zero in one direction of the field space for each doublet.Then we can write the tree-level potential for the scalar background fields as In principle, there are two independent field directions, however, as was discussed in the literature [11,16,17], the tunnelling proceeds along the direction of the new scalar field φ and subsequently the Higgs field h rolls to the true vacuum.Therefore, in our analysis we will focus solely on the φ-direction.
In the one-loop correction to the effective potential, the dominant contribution comes from the dark gauge bosons, X µ .Therefore, the one-loop potential along the φ direction at zero temperature reads with M X (φ) = 1 2 g X φ.The scalar one-loop contributions can be neglected since they are subdominant [11].

Thermal effective potential at leading order
The thermal one-loop correction to the effective potential along the φ direction reads and the daisy correction in the Arnold-Espinosa scheme is given by [32] V daisy (φ, where [17] Since scalar loops are negligible, the only contribution affected by the daisy resummation is that of the zero Matsubara mode of the dark gauge field temporal component, which acquires a Debye mass.The full one-loop potential with daisy corrections is given by the sum of the contributions listed above The label "LO" stands for leading order, as we will see that this potential coincides with the leading order effective potential obtained in section 4.1, as long as the matching is performed at leading order only (see section 4.2).
As explained in the previous section (see also ref. [11]) in models with classical scale invariance we should use the RG-improved effective potential, since vastly different scales are present in the model.We improve the potential by replacing the field and couplings with their running versions as φ → Z(t)φ (see the discussion in sec.3.3), λ φ → λ φ (t), g X → g X (t), where t = log µ µ 0 and µ 0 corresponds to the Z boson mass.The β functions and the anomalous dimension for the φ field are listed in appendix A.
Next we choose the scale µ as stated in eq.(2.15), with no running included in the field dependent mass and g X defined at the scale This choice ensures that at large field values the field-dependent logarithmic term is (almost) cancelled by the field-dependent renormalisation scale, 6 while at lower field-values, for φ < 2κT /g X , the field-dependent logarithms cancel between the zero-temperature and finite-temperature contribution, and the scale is fixed to κT , cancelling a remaining Tdependent logarithm.

Tunnelling and normalisation of the field
When defining the potentials we start from setting the values of the mass of the X boson and its coupling g X at µ = M X .Following the procedure described in ref. [11] we recover the values of all the couplings at the electroweak scale set by M Z and we define the theory at that scale.This is reflected in the choice of the reference scale µ 0 , see the discussion above eq.(3.7).This means that the field φ is defined at µ = M Z , i.e. at this scale it is canonically normalised (for a comprehensive discussion of the scale-dependence of scalar fields, see ref. [70]).As we RG-improve the effective potential, we evolve the couplings and fields along their RG flows.This means that at other scales the field is not canonically normalised and the field is rescaled by the field renormalisation constant, Z(t).In the HT regime the running is frozen at µ = κT and the normalisation of the field is given by Z(log κT M Z ).At the same time, the usual bounce equation is derived from an action containing a canonically normalised scalar field, see eq. (2.13).This means that we cannot simply use the bounce equation with the RG-improved potential evolved down to the thermal scale.We could rederive the bounce equation in terms of the field defined at µ = M Z but we choose to rather redefine the field for the purpose of solving the bounce equation and computing the action -we introduce a new field that is defined at the scale µ = κT and is thus canonically normalised.It is related to the old field via the rescaling Remembering about this subtlety, we will not introduce extra subscripts on the field symbol φ for simplicity of notation. 7he factor of √ Z may seem unimportant, as Z stays close to 1, as long as we are within the perturbative regime of the theory, however, it turns out that it affects the results visibly and we therefore take this Z into account in the LO computations of section 5.This issue has not been appreciated in the literature, see for example ref. [11].

High-temperature effective theory
As mentioned in the introduction, thermal field theory suffers from poor convergence of the perturbative expansion, which can hamper the precision with which the properties of the phase transition are determined.Even though daisy resummation (see eq. (2.9)) resums a leading set of IR-sensitive diagrams, and is hence correct at O(g 3 X ), a problem persists: parametrically large O(g 4 X ) contributions are still missing.In particular, there is an uncancelled RG-scale dependence due to the omission of two-loop thermal masses, and furthermore additional resummations are required at the same order. 8The root of both problems lies in the Bose-enhancement of the low-energy modes, resulting in an enhancement of the effective parameters of these modes.

Effective field theories
A systematic way to deal with the above-mentioned problems, is to construct a series of effective field theories describing the thermodynamics at the different relevant energy scales.See ref. [45] for a recent discussion of the possibly relevant scales.For our purposes, it suffices to distinguish the following two energy scales (see table 1):

Name Energy scale Scaling of expansion parameter
Hard πT Soft gT g π Table 1: Relevant energy scales for the SU(2)cSM model (see section 3) at finite temperature.Conventionally, g denotes the largest relevant coupling in the theory, in our case g X .
The hard scale For the construction of the EFTs, we make use of the partition function given by Here, τ is the imaginary time coordinate and its periodicity is set by the reciprocal of the temperature.The functional integration DΦ is performed over all fields.L E denotes the Euclidean space Lagrangian density.The fields can be written as a sum over momentum modes, with momenta P = (ω n , p), with the Matsubara frequency ω n = 2πnT for bosons and ω n = (2n+1)πT for fermions.The theory described by eq.(4.1) contains all momentum modes, but modes with masses larger than πT get Boltzmann-suppressed, so we see that the largest relevant energy scale in the HT regime 9 in the problem is the so-called hard scale of O(πT ).
The soft scale We can obtain the effective theory at the soft scale by formally integrating out all n ̸ = 0 momentum modes with ω n ⩾ πT .The resulting partition function for the theory containing only the scalar fields and gauge bosons is given formally by Here the path integral is over the zero modes only.The 1/T factor coming from the integral over τ is absorbed by the fields in the 3D Lagrangian such that the exponent is dimensionless, see the discussion below eq.( 4.4) and e.g.ref. [51].The fields in the effective Lagrangian L soft 3 are static and three-dimensional; they carry no momentum in the imaginary time direction.f 0 is the coefficient of the unit operator, related to the pressure in the symmetric phase [36].
As we will see explicitly below, the zero modes of the temporal components of the gauge fields obtain a Debye mass from the screening by the hard modes.This mass is of the order m D ∼ gT , with g the relevant gauge coupling.This mass scale defines the so-called soft scale, which is the largest energy scale of the EFT constructed by integrating out the hard modes.The spatial components of the gauge fields do not get screened, and are thus lighter than the soft modes. 10The screening of the hard scale also generates a mass ∝ gT for the scalar field.
Let us now get more explicit.The action of the 3D EFT -the action obtained after integrating out the hard modes -separates into S soft 3 = S soft, dark 3 + S soft, SM 3 + . .., where the ellipsis denotes the contribution from a portal that couples the two sectors together.From now on, for simplicity, we focus solely on the dark sector part within the EFT as that is what we need for the computation of the bubble nucleation rate.Yet, it is good to keep in mind that we still capture the SM contributions coming from the hard scale in the matching relations.For the dark SU(2) sector, the action reads ϕ now denotes the scalar field and F a ij is the gauge field strength tensor of the spatial gauge field X a i with spatial Lorentz indices i = 1, 2, 3 and SU(2) isospin index a = 1, 2, 3.The gauge coupling of the EFT is denoted by g X, 3 .The temporal components of the gauge field, X a 0 , are Lorentz scalars in the EFT and transform in the adjoint representation of SU (2).The covariant derivatives for scalar doublet and triplet are , respectively.The scalar potential reads Note that within the EFT, fields have dimension T 1 2 , but we do not give them an explicit label "3", and the mass terms m 3 , m D,X and the couplings λ 3 , κ 3 , h 3 have dimension T .We include marginal operators for the doublet ϕ with coupling constants c n , which correspond to the terms containing the ζ-terms in eq.(2.6).We use these marginal operators only as an indicator of a breakdown of the HT expansion: when their effect becomes non-negligible at low temperature, or more importantly, at large field values, the HT expansion starts to break down.
The parameters of the 3D theory are obtained by a matching procedure, as is common in the construction of any EFT.We use DRalgo [38] for the determination of the parameters of the soft scale EFT, and list the result in appendix B.2.We highlight that the momentum-dependent field normalisation contributions of the hard modes are absorbed into the parameters of the EFT for all fields, rather than including Z-factors in the kinetic terms of the soft scale EFT action [35].For illustration of the matching procedure, see e.g.refs.[35,71], or appendix B.1 in ref. [29].We emphasise that the construction of the 3D EFT is performed in the symmetric phase, and relies on the high-temperature expansion for the matching, which assumes that m/T ≪ 1 for all the fields.In conformal models all fields are massless in the symmetric phase, yet deep in the broken phase the assumption about the HT expansion is no longer valid, as discussed in section 2.2.
We now assign a background field to the scalar field resulting in masses for the spatial and temporal gauge bosons respectively The tree-level potential for the field v 3 at the soft scale is then given by Note that as explained above, the marginal operators with coefficients c 2n are suppressed in the HT expansion and hence of higher order, and included only for inspecting departure from the HT regime at large field values.In order to obtain a cubic term, we have to integrate out the gauge field modes, resulting in a new EFT expansion for the effective potential and the effective action.This is a generic feature of models where the tree-level potential does not include a barrier11 required for first-order phase transition [39,[43][44][45]72].Since all the masses in this theory are formally soft, it is not immediately obvious that we can integrate out the gauge modes.We return to the issue of mass hierarchies in section 4.6.
Integrating out the gauge field modes at one-loop order results, together with the tree-level term of the soft scale EFT, in the leading order contribution of the final EFT expansion12 6(m 2 X,3 ) Here the last, field-independent, term has been added to normalise the potential to zero for a vanishing field.The validity of this EFT expansion will be further discussed in sections 4.5 and 4.6, and higher-order corrections are determined in section 4.3.

Relation to one-loop thermal effective potential
One can confirm straightforwardly that eq.(3.6) in the HT expansion is reproduced exactly by eq.(4.8) when the matching relations (c.f.appendix B.2) are truncated as follows ) T 2 (see eq. (B.10)).We note that while we simply took eq.(2.9) for leading daisy resummation from the literature, in eq.(4.8) we actually derived it in the 3D EFT approach and we can clearly see the physics behind it: this term originates from the fact that hard modes screen the soft zero mode, and once this hard scale screening is accounted for by the soft mode mass at one-loop order, the EFT automatically resums this contribution to all orders.In the EFT, a one-loop computation with the resummed propagator is easy and, furthermore, two-loop diagrams are also straightforward, and we exploit this fully in the next section.

Higher-order corrections within the EFT
Now that we have illustrated what is behind the EFT approach, let us highlight what happens when we do not truncate the matching relations: in this case, we account for thermal corrections of the hard scale screening in not just the mass, but also in all couplings, at one-loop order.This includes momentum-dependent contributions, i.e. we account for field renormalisation factors, see ref. [35].Even more importantly, we include two-loop corrections in the thermal mass, and together with field renormalisation factors this guarantees that 3D EFT parameters are renormalisation scale invariant to the order we compute.
The construction of the NLO contribution to the potential follows the prescription of refs.[39,41,46].We include two-loop contributions of the soft modes, which yield the NLO correction to the effective potential [38,45] The last line is independent of the field value, and ensures that the potential is zero at the origin.We will denote the sum of the LO and NLO potential as

.15)
In addition to the corrections to the effective potential, integrating out the gauge modes also results in a field normalisation term for the effective action (note that we normalise the potential such that it is zero at the origin, which also implies discarding f 0 in eq.(4.3)) with Z 3 (v 3 ) given by (see also refs.[73][74][75]) Within the EFT, the term with the Z 3 -factor is an effective derivative operator generated by the gauge modes.The role of this Z 3 -factor within the EFT is different from field renormalisation in the parent theory, as it is not related to UV running.
Let us pay close attention to the soft 3D EFT RG-scale, µ 3 .We can split the action into two parts S EFT

3
, where where NLO is suppressed compared to LO by the soft expansion parameter, formally O( g X π ).The LO action depends explicitly on the 3D EFT mass parameter, m 3 , which runs (see eq. (B.6)). 13The NLO part of the action is independent of the scalar mass, yet contains explicit logarithms of the RG-scale.In ref. [41] it has been shown that the action is RG invariant at the order considered (O(g 4 X )).For completeness, we reproduce this argument here.We simply have as the running of the scalar mass is governed by the beta function which exactly cancels the µ 3 -dependent terms in V EFT,NLO

3
. In our computation, we use the freedom to choose the RG scale and use µ 3 = κ RG m X,3 , where the coefficient κ RG is of the order unity.Should one aim to optimise the choice of RG-scale, one could follow e.g.refs.[59,76].

Thermal parameters at NLO
To compute the exponential part of the nucleation rate in the EFT, we follow refs.[39,41].We use a method based on the strict expansion of the action in order to obtain results that are gauge invariant.
In the strict expansion, the critical bubble configuration v 3,B is formally expanded as 2 ) and the LO bounce v 3,b is found only using the LO effective potential: with boundary conditions The NLO action is then simply evaluated at the leading order bounce, and higher-order corrections to the bounce result in contributions that are formally beyond our accuracy goal [41,74].It has been shown in refs.[40,41] that despite the singularity in the Z-factor at zero field value, its contribution to the action is finite.
The bounce solution depends on the scalar mass and hence inherits RG-scale dependence through it.Note, however, that the implicit running of the bounce does not contribute at the order considered, because the LO action is extremised by the bounce, i.e. when applying the chain rule in eq.(4.20), we see that the term µ 3 For illustration, let us use a simplified expression for the nucleation condition, i.e. that the nucleation rate equals the Hubble parameter, H A(T n )e −S(Tn) Here we will use i.e. we estimate the prefactor simply using dimensional analysis since we do not compute this contribution properly, but only assume it is suppressed compared to the exponential part.M pl is the reduced Planck mass, and ξ g = 30/(g * π 2 ) with g * the number of degrees of freedom in the plasma.For convenience, we define the shorthand notation We can determine the nucleation temperature by the mixed method14 (c.f.ref. [45]) where we compute the action in the strict expansion, yet we directly solve for the nucleation temperature from the condition This method is not only gauge invariant, but also invariant with respect to the soft 3D EFT renormalisation scale, since the sum S 0 + S 2 is invariant.In this sense, only using the sum S 0 +S 2 describes the full NLO soft scale corrections.In our numerical analysis we apply this method to compute the nucleation and percolation temperatures, i.e. we evaluate the NLO action at the LO bounce solution and use the resulting S 0 (T ) + S 2 (T ) with the standard formulas for nucleation and percolation temperatures, following the approach described in ref. [11].
In analogy to ref. [45], we could use the strict expansion for the nucleation temperature as well, by formally expanding in δ (which is set to 1 at the end of the computation): S = S 0 + δ 2 S 2 + O(δ 3 ) and A = δ 3 A 3 + O(δ 4 T 4 ) [39], and furthermore expand T n = T n,0 + δT n,1 + δ 2 T n,2 + O(δ 3 ).Here δ denotes suppression with respect to S 0 .We can then expand the simplified nucleation condition, and by equating the first two non-vanishing orders we find that T n,0 is given by and the first non-zero correction to it reads where all quantities are evaluated at T n,0 and in the second expression we assumed again A 3 = T 4 .However, we do not find this expansion useful as the sum T n,0 + T n,2 appears not to be scale invariant, and neither are the individual terms and we therefore choose to use the mixed method. 15.5 Power counting and validity of the EFT Now that we have given our expressions for the LO and NLO effective potential, as well as the approach to determine the nucleation rate and temperature, let us understand better in which sense eq.(4.8) describes the LO behaviour, and what kind of power counting this implies for perturbation theory.
For starters, let us consider the scaling of different contributions to the potential in the presence of a radiatively generated barrier.The existence of the barrier requires [32] i.e. all terms in the potential are of the same order.As discussed in detail in ref. [45], one possible realisation of this occurs in the temperature regime where the scalar mass parameter is parametrically lighter than soft, in particular m 2 3 ∼ g In such case, the gauge field modes are soft m X,3 , m X 0 ,3 ∼ gT , and eq.(4.8) is interpreted as the construction of the supersoft EFT below the soft scale.In ref. [45] such effective description exists, since the mass parameter has the schematic form m 2 3 = m 2 0 + Cg 2 T 2 , where C is a positive constant and m 2 0 is a negative zero temperature mass parameter.The partial cancellation of these two terms, which are individually of order (gT ) 2 , makes it possible that the effective scalar mass is parametrically lighter than soft in some temperature range.The effective description for the phase transition is constructed in such a temperature window, and furthermore the background field is assumed to scale as v 3 ∼ √ T in the broken phase, such that the gauge modes are soft.The potential itself follows the scaling where corrections are due to marginal operators with c 2n ∼ g 2n π 2(n−1) [29].At parametrically smaller field values v 3 ∼ √ gT the effective description is compromised.Indeed, in the terminology used in ref. [39] gauge field modes are scale-shifters: at field values close to the symmetric phase, in the bubble tail, the effective description based on a derivative expansion of the effective action alone fails, yet can be used at the bulk of the bubble to compute the nucleation rate [39,41,74,77,78].We will return to this issue below.
In our case of a dimensionally transmutated theory, the previous discussion becomes more subtle and needs to be modified, since the scalar mass parameter cannot be lighter than soft m 2 3 ∼ (gT ) 2 , since there is no zero temperature mass m 2 0 .In this case, eq.(4.28) implies λ 3 ∼ T . 16Such a huge value of the background field pushes the gauge field modes formally to the hard scale m X,3 , m X 0 ,3 ∼ πT and furthermore signals a possible breakdown of the HT expansion, as the scaling of the potential becomes and marginal operators are no longer strongly suppressed, but could contribute at leading order.Does this mean that the effective description based on the effective potential that we derived cannot be used after all?We argue that this is not the case, for the following reasons: Even though the aforementioned formal power countings for masses and the background field can bring clarity about how to organise the perturbative expansion, it is not clear how strictly they should be followed.As scale-shifters, the gauge field modes vary over multiple mass scales, from the hard to the soft to eventually the non-perturbative ultrasoft scale ( g 2 π 2 T ).As the bounce solution interpolates between the two phases, the bubble nucleation rate obtains contributions from different scales.Intuitively, the EFT description could capture most of the effects reliably, provided that 1) Non-zero Matsubara modes are much heavier than zero modes.
2) Gauge field zero modes are much heavier than the nucleating scalar field zero mode.
In the following subsection we study the mass hierarchies in our model in more detail and explain how we treat the scale-shifters.

Mass hierarchies and scale-shifters
Let us study the typical mass hierarchies of the problem, which we encounter when using the potential of eq.(4.8) to find the bounce solution using the schematic figure 1.In this figure, the masses of the gauge field modes m X,3 , m X 0 ,3 are given as a function of the critical bubble radial coordinate r, together with soft mass parameters m D and m 3 .17Indeed, both m D and m 3 are parametrically of the order gT , i.e. soft, yet note that due to group theory factors in the LO dimensional reduction matching relations they differ approximately by a factor 2, c.f. eqs.(4.12), (4.13).Close to the center of the bubble at small r, m X,3 , m X 0 ,3 indeed become very large, yet they are still below the lightest bosonic non-zero Matsubara mode with mass 2πT at the escape point.In this case, one needs to be very cautious with the HT expansion.However, we demonstrate in sec.5.1 that indeed up to the escape point, the HT expansion converges well and as a consequence also marginal operators are suppressed, providing support that the EFT picture is reliable in the small r regime.
On the other hand, we know that the EFT picture fails at the bubble tail r > r t , where r t is defined at the radial distance where the spatial gauge mode mass becomes comparable with the nucleating scalar mass, suggesting that it is not possible to integrate out the gauge modes.Therefore, we can trust the EFT that we constructed for field values above Finally, we emphasise the following: as long as points 1) and 2) stated at the end of section 4.5 are valid, the higher-order corrections to eq. (4.8) given in eq.(4.14) are of the same form regardless of the assumed formal power counting for the gauge field modes and the background field.Indeed, the EFT expansion has the same functional form as long as 2πT ≫ m X,3 , m X 0 ,3 ≫ m 3 ≫ g 2 π 2 T , yet we do not need to fix the formal power counting for these in-between scales between the hard and ultrasoft scales, and indeed we cannot, since the gauge field modes are scale-shifters [39].
In ref. [39] it is explained in detail how to treat the nucleation EFT construction with a scale shifter.In essence, the one-loop contribution from the gauge fields is still resummed to the LO effective potential to provide a barrier, and this affects the LO bounce solution.For the action however, their contribution is computed without derivative expansion, i.e. they contribute in the prefactor in analogy to the soft scalar modes.Then, one needs to subtract this gauge field contribution from the exponential part of the rate, to avoid double counting.In our analysis, we do not compute the prefactor, and hence we stick to the procedure described earlier in section 4.4.We check the accuracy of this approach by estimating the contribution of the tail of the bounce, where the assumed mass hierarchies are violated, for details see section 5.4.2.

Numerical results
In this section we present the results of a scan of the entire allowed parameter space using the mixed method described in section 4.4 for computing the action with the NLO effective potential of eq. ( 4.15) and the Z 3 of eq.(4.17).We use this action to compute the nucleation and percolation temperatures, as well as the normalised radius of bubbles at the moment of collision (R * H * ) and the efficiency factor for bubble collisions following the procedures described in ref. [11].We also evaluate the expected observability of the signals in terms of the signal-to-noise (SNR) ratio, using the spectra for supercooled phase transitions from ref. [10].Moreover, we compare the NLO predictions to the LO ones, obtained from eq. (3.6), updating them with respect to ref. [11] by including the thermal cutoff on the running also in the zero-T part of the effective potential (see section 2.2 for details) and by redefining the field to be canonically normalised at the thermal scale (see section 3.3).

Effective potential at LO and NLO
Let us start with comparing potentials computed with different approximations, all evaluated at T n computed from the NLO potential.In the left panel of figure 2 we focus on the low-field value or high-temperature regime, the range for the plot is chosen such that the barrier is well visible.The blue solid line represents the full one-loop potential of eq.(2.8) (see also eq.(3.6)), with the daisy term included.It agrees very well with the high-temperature approximation of eq.(2.14) (long-dashed light green).The NLO potential computed within the EFT of eq.(4.15) (dotted red line) differs from the LO result mildly, while exclusion of the daisy term (dashed green line) modifies the result significantly, which indicates that the daisy diagrams are indeed very relevant for the shape of the potential. 18The right panel of figure 2 shows the large-field or low-temperature behaviour of the effective potential.The full one-loop LO potential of eq.(2.1) (solid blue line) is closely approximated by the low-temperature approximation with only the first term in the sum in eq.(2.7) (long-dashed light green line).
From figure 2 we learn that the HT expansion works perfectly around the barrier but we should ask whether it is valid all the way up to the escape point for the bounce trajectory and beyond.This is shown in figure 3. The left panel presents the potential evaluated using different approximations around the escape point, which is marked by the vertical grey line (obtained with the NLO action).The solid blue line shows the full one-loop LO effective potential, while the dashed lines indicate the usage of HT expansion of the thermal function, eq.(2.6) (see also eq.(4.8)): long-dashed, green is the first approximation without the sum containing the ζ-terms in the second line, short-dashed light green includes the first term in the sum, while the dotted red line includes the first three terms from the sum.It is clear that all the approximations agree very well in the vicinity of the barrier and beyond, only on the verge of the displayed region small differences between the curves can be noticed, as the first approximation to the potential deviates slightly from the full solution.For larger field values, shown in the right panel of figure 3, the HT expansion is quickly invalidated which is clear from the fact that the approximations with more terms from the HT expansion included behave worse than the ones with fewer terms.At the same time, the LT approximation (dot-dased orange curve) works like a charm (it overlaps with the solid blue line representing the full potential).This is an explicit confirmation of our earlier claims that the high-temperature expansion can be used for the field values relevant to the tunnelling, while the low-or zerotemperature potential describes accurately physics associated with the minimum of the potential.

Phase transition and GW signal
For strongly supercooled phase transitions the only parameters that are relevant for the determination of the gravitational wave spectrum are the length or time scale of the transition evaluated at the percolation temperature and the reheating temperature.
The process of reheating is controlled by the decay rate of the scalar field, which measures its ability to transfer the energy to the SM plasma.When it is larger than the Hubble parameter, reheating can be considered instantaneous.Then the Universe reheats to the temperature at which thermal inflation started, which is controlled by the potential in the low temperature/large field regime.In ref. [11] it was shown that this is the case for most of the parameter space of the SU(2)cSM model.Only for low g X and high M X , the decay rate of the scalar field becomes smaller than the Hubble parameter.In this case the Universe cannot reheat immediately after the phase transition.This results in a period of matter domination when the scalar field oscillates around the minimum until the decay rate becomes large enough to transfer the energy to the SM sector [7,79,80].In this scenario of inefficient reheating, the final temperature is lower and the GW spectrum is modified by the modified expansion history.In ref. [11] the region of inefficient reheating was excluded by the percolation criterion.As we will see, this region opens up by including the NLO effects, however, it is still very small and we will not analyse it in detail since it is beyond the main focus of the present work.Thus, we will assume that the reheating temperature is given by T V and will not be changed compared to the analysis of ref. [11].Therefore, we will not show the results for T r here.
Below we present the results of the scan of the parameter space obtained using the theory given by eqs.(4.8), (4.14), (4.17).In figure 4 (upper panel) we present the values of the percolation temperature and the average radius of a bubble at the moment of collision (length scale of the transition) obtained using the EFT NLO potential with µ 4 = πT .We exploit the invariance of the potential with respect to the 3D RG scale and choose a fielddependent value for it as µ 3 = m 2 X,3 (v 3 ) (see eq. (4.6)).The excluded regions are marked in shades of grey.The leftmost part, corresponding to small M X , does not reproduce the electroweak vacuum correctly.The upper right part of the parameter space is excluded as there the g X coupling becomes nonperturbative at some scale between M Z and the QCD scale µ ≈ 0.1 GeV.In the lower part, the phase transition is triggered by the QCD condensate [7,14,16,68,[81][82][83][84].We assume this happens for T p < 0.1 GeV and this region is beyond the scope of the present work.Finally, in the light-grey triangle-shaped region we cannot assure the completion of the phase transition as the percolation criterion is not fulfilled.It was shown that supercooled phase transitions can produce primordial black holes (PBHs) [85,86], and for sufficiently slow transitions i.e. β/H * ≲ 6 − 8 it would cause PBHs overabundance.We indeed find such values, but only in the right corner of the non-percolation region, for large masses M X .The former two constraints do not depend on the temperature, therefore they will be identical in all the plots.On the other hand, the latter two (regions where QCD effects become relevant and the percolation criterion is not satisfied) depend on our predictions for the phase transition and can change depending on the approach.In the lower panel of figure 4 the same quantities computed from the LO potential of eq.(4.8) (with matching conditions truncated as in eqs.(4.9)-(4.13) to match the usual daisy-resummed approach) potential are presented.
The values of the percolation temperature obtained at NLO are between 0.1 GeV and about 380 GeV.As a general trend, it can be noticed that the percolation temperature goes up, as compared to the LO prediction.This extends the available parameter space to lower values of the gauge coupling g X .Moreover, the region of non-percolation is pushed to higher values of the X mass.This opens up a small region of the parameter space where percolation is possible, but reheating is inefficient.We do not study this effect in detail as it is beyond the scope of the present paper.
To better evaluate the differences between the LO and NLO predictions we present the relative differences between them (normalised to the NLO ones) in figure 5. We can see that the change in the percolation temperature between the two approaches is significant, ranging from O(50%) in the low-mass, large-coupling region, up to O(100%) in the smallcouplings, large-mass corner.This seems somewhat counter-intuitive, as one expects the largest corrections between the two methods for large couplings, but note that this concerns the coupling at the thermal scale.It should be stressed here that the coupling and mass displayed in figure 5 are defined at the scale µ = M X .They need to be RG-evolved to the thermal scale.In the large-mass corner, the coupling becomes significantly larger at the thermal scale, which explains why the difference between the LO and NLO approaches is largest in this part of parameter space.Let us point out again that the value of T p does not directly affect the GW spectrum, so the large corrections we find in the NLO description are not reflected in a strong modification of the GW signal.However, it signals that the differences between the descriptions at different orders in perturbation theory are non-negligible.
The right panel of figure 5 displays the difference in the typical length scale of the transition, R * H * .Here we see that the differences are smaller than for T p , but again the largest differences are observed in the large-mass corner, reaching O(55%).Even though the relative difference between LO and NLO for R * H * is smaller than for T p , these differences do modify the GW prediction.It can be seen in refs.[10,11] that both the GW amplitude and the peak frequency depend on R * H * and we, therefore, expect the predicted spectra to be shifted compared to the earlier results.Given the values of R * H * we expect the GW signals to still be well visible in LISA.This will be verified by computing the SNR.
One can also note the change in the overall allowed region indicated in figure 5 by the white dashed and dotted lines.The region of non-percolation is significantly shifted and also the region where the phase transition is expected to be sourced by QCD effects is pushed to lower values of g X at NLO.Therefore, the predictions for the GW signal in this region could be significantly altered by using the LO or NLO approach.
To formulate the predictions for the GW spectra, we should ask the question about the source of the gravitational waves -for which region of the parameter space are they produced via collisions of bubbles and where by sound waves? 19Our predictions for the efficiency factor for producing GWs via sound waves, κ sw based on the LO and NLO potential are shown in figure 6.It is clear that, depending on the region of the parameter space, GWs can be sourced by either of the mechanisms, and the difference in the predicted source between the two methods is limited to a rather narrow range of the parameter space.This is interesting, since in ref. [11] for scans performed at a fixed renormalisation scale no region was found where the dominant source would correspond to bubble collisions (as opposed to predictions based on the RG-improved potential).In our current approach, the scale is fixed in the HT regime, where the tunnelling takes place.However, it is fixed to κT so effectively it is different for every point in the parameter space and is proportional to the percolation temperature.This suggests, that allowing the scale to change is crucial for seeing bubble collisions. 20aving checked the predictions for the parameters of the transition, let us check what are the implications for the signal-to-noise (SNR) ratio predicted for LISA.The results of the LO and NLO approach are presented in figure 7.As expected, the values of SNR are very high throughout the entire allowed parameter space, implying that a first-order PT in the SU(2)cSM model should be well visible at LISA. 21At NLO we observe a slightly lower SNR, around 10, at the edge of the parameter space corresponding to large M X .The reason is that the peak frequency of the spectrum is higher for higher reheating temperatures.The latter grows with M X so for the largest values of M X allowed at NLO (excluded at LO by the percolation criterion), the signal moves out of the sensitivity range of LISA.The very high values of SNR imply that, if a signal from a phase transition in SU(2)cSM is observed, we will be able to reconstruct the values of T r and R * H * with very good precision, for very strong signals even better than 1% [48].Because of this prospect, the differences in predictions between the LO and NLO approaches, see figure 5, will be much greater than the experimental uncertainties and this highlights the importance of efforts to include higher-order thermal corrections.

Implications for dark matter abundance
The SU(2) X model contains DM candidates -the new gauge bosons are stable due to a residual SO(3) symmetry and can in principle attain the correct relic density in different parts of the parameter space via the freeze-out mechanism [11, 16, 23, 63-66, 68, 69] or the super-cool DM mechanism [11,16,68,69].The DM phenomenology is not the focus of the present paper but we will comment on how the NLO modifications of the allowed parameter space affect the existing results.In the lower-mass regime, the DM relic abundance is produced via the thermal freeze-out mechanism and a rather large gauge coupling is required.Reference [11] showed that the correct relic abundance, in agreement with the current direct-detection limits, can be attained for 1.2 TeV ≲ M X ≲ 1.8 TeV and 0.82 ≲ g X ≲ 0.96 (similar results were found in ref. [23]).The NLO results extend the allowed parameter space, where a first-order phase transition happens independently of QCD effects, to lower values of g X so we do not expect additional constraints on the region with the correct thermal DM abundance presented in refs.[11, 16, 23, 63-66, 68, 69].
On the other hand, the mechanism of super-cool DM requires inefficient reheating, such that the temperature of the Universe after the phase transition is below the decoupling temperature of the X particles.This is realised for larger M X , approximately above 3000 TeV, and lower g X around 0.7 (see e.g.figure 8 of ref. [11]).This region was excluded by the LO analysis of ref. [11], however, the NLO results re-open a small part of this regime, as discussed in the previous subsection.Therefore, a small region with the correct relic density obtained via the super-cool DM mechanism (supplemented by a subthermal population produced via scattering) may be possible.However, to give a definite answer requires a dedicated computation of the reheating temperature and the resulting DM abundance.

Evaluation of the uncertainties
In section 4.5 and 4.6 we have discussed the challenge of constructing an accurate EFT in a theory with scale-shifting fields, and the inaccuracy associated with the contribution from the bubble tail.Moreover, there is an uncertainty associated with the omission of higher-order corrections, which we can study by varying the RG scale.In this section, we quantify the uncertainties associated with our computation of the thermal parameters.

Dependence on the renormalisation scale
As was explained earlier, the common approach of using the one-loop effective potential with daisy resummation suffers from an uncalled RG-scale dependence that can be cured by the inclusion of certain two-loop level diagrams.This is achieved in our NLO effective potential (with the matching also performed at two-loop level).There, the RG-scale dependence cancels up to terms of order higher than the order to which we compute.Moreover, the potential (and the full action) are independent of the 3D scale, up to higher-order corrections.As was shown in the literature [28,29], omission of significant perturbative corrections (revealed by the scale dependence) is the main source of uncertainty in predicting the GW signals.We are now in a position to check the RG-scale sensitivity of the NLO predictions and to contrast it with the LO result.Note, that with the RG-improvement procedure implemented in this work, when we say that we change the 4d scale, in fact we mean that we change the thermal cutoff in the running in eq.(3.7).So changing the scale from πT to 2πT means changing κ in eq.(3.7) from κ = π to κ = 2π.
Figure 8 (left panel) presents the relative difference in T p obtained from the NLO action at two different scales, µ 4 = πT and µ 4 = 2πT .We observe a mild dependence on the 4D scale, the result for T p changing, between the two RG-scales, by at most 10%.The changes in R * H * are much smaller, and they never exceed 2%, therefore we do not show the plot illustrating this difference.We have seen that the predicted SNR for LISA for this model is large, implying that thermal parameters can be reconstructed with very good precision.It therefore needs to be determined if the 2% uncertainty in R * H * leads to an observable difference.
For comparison, figure 8 (right panel) presents the dependence on the scale of the LO results.We can see that the change in results for T p is much larger -the relative difference is approximately between 15% and 30%.This confirms our earlier claims that the inclusion of the NLO corrections cancels the residual scale dependence present at LO.The RG-scale dependence of the bubble radius is again milder and is of the order of 5% in the whole parameter space, which is again larger than the uncertainty in the NLO result.
In both approaches, LO and NLO, the overall allowed region is only slightly modified by changing the 4D scale as indicated in figure 8 by the white dashed and dotted lines.

Importance of the Z 3 -factor and the bubble tail
Looking at the rather good agreement between the different potentials in figure 2, the sizeable differences in T p and R * H * observed in section 5.2 might come as a surprise.It turns out, that the largest cause of the difference between the LO and NLO descriptions is the Z 3 -factor multiplying the kinetic term.As observed from its explicit form in eq.(4.17), Z 3 diverges as v 3 → 0. This corresponds to the regime where the derivative expansion of the action is no longer valid, and which our description hence cannot capture, see ref. [77].Since along the bounce trajectory, as v 3 (r) → 0 also ∂ i v 3 (r) → 0, and S 2 thus remains finite.However, by comparing different contributions to S 2 shown in figure 9, we can see that the contribution of Z 3 is still dominant in S 2 .To quantify the effect of Z 3 we have computed the thermal parameters in the mixed method with the NLO contribution to the potential, but without the Z 3 for a representative set of parameter choices.In all cases, we observe that the NLO approach gives a correction to T p only of the order of ∼ 25% compared to the LO result, which is much smaller than the corrections observed in figure 5 and this confirms explicitly that the main correction comes from Z 3 .
A related question, which was already discussed in section 4.1 and 4.4, is the fraction of the action coming from the region with r > r t .A large contribution from this region signals several problems.First, the solution gets a large contribution from the Z 3 -factor in a region where its expression is not valid.Second, the EFT breaks down, due to the scale shifting nature of the fields.Last, large contributions coming from the kinetic term suggest that the derivative expansion, which allowed us to compute the bounce in a momentum-independent potential background, breaks down.In section 4.4 an approach was suggested to estimate the contribution from the region r > r t .For the 3D RG-scale given by µ 3 = m X,3 , we find for a representative set of benchmark points that the fraction of the action given by r > r t is of the order 30 − 40%.For the benchmark point considered in figure 9, r t = 0.29 GeV −1 and the tail's contribution to the action is 31%.This is a significant fraction, and it should motivate us to investigate the validity of our expansion further.However, this estimate should not be taken too literally, as it depends on the choice of µ 3 .The fraction can be made smaller, at the cost of increasing S 2 with respect to S 0 .
The way forward is to include the contributions of the gauge and scalar modes in the functional determinant of the nucleation rate, following the approach of refs.[77,88].This will allow us to assess better the validity of the derivative expansion, and we leave it for future work.

Summary and outlook
This work is devoted to the accurate theoretical description of supercooled phase transitions, with the SU(2)cSM model [4,47] as a concrete example.The main motivation for this study is the great prospect for observability of a GW signal from supercooled firstorder phase transitions [4,5,[12][13][14][15][16][17][18][19][20][21][22], which would be so strong that spectral parameters could be reconstructed with a very good accuracy [48].This exciting possibility calls for an increased accuracy in the theoretical description, which is not accessible with the popular machinery of the daisy-resummed one-loop effective potential.A perfect tool for increasing the accuracy of phase-transition-related predictions is dimensionally reduced effective field theory [33][34][35][36] which allows us to perform resummations systematically.Therefore, we pursued the task of reconciling DR, which is based on a high-temperature expansion, with the description of supercooled phase transitions.This work resulted in several new findings which are summarised below.
In section 2.2 we demonstrated that the relevant quantities, describing the phase transition and setting the GW signal, can be divided into two groups: large-field-related (∆V , T r ) and small-field-related (T n , T p , R * H * ).The former ones correspond to the low-temperature limit of the effective potential and no resummations are needed to compute them and we can follow the approach described e.g. in ref. [11].The latter ones are related to the hightemperature regime.To compute them accurately we need a high-temperature effective field theory which takes into account the hierarchies between different energy scales in the presence of high temperatures.
As we studied the relation between the HT and LT regimes, we also elucidated certain points in the computation of the thermal parameters of the phase transition in the 4D approach.First, by studying the interplay of the RG-improvement of the potential, and RG-scale cancellations between the HT limit of the thermal contribution to the effective potential and the zero-temperature part, we came to the conclusion that there is a preferred scale for the phase transition computations in the 4D theory which is the thermal scale, µ ∼ T , see section 2.2.This may sound like a trivial observation, but it is not commonly employed in the computations found in the literature.Moreover, we explained the role of the normalisation of the field in the phase-transition-related computations in section 3.3.
Furthermore, we have checked that the bounce solution, corresponding to the tunnelling trajectory of the field, is always within the HT regime.However, as we move along the bounce, the considered masses change substantially, which causes difficulties in treating mass hierarchies.The constructed EFT is expected to break down in the tail region, which is the region where the gauge field mass becomes smaller than the scalar field mass.In this region, contributions of the gauge field modes to the action should be accounted for without derivative expansion.These issues were discussed in sections 4.5, 4.6 and 5.4.2.
In sections 4.1 and 4.3 we have constructed an EFT in the HT regime, with matching at two-loop level and going to the NLO level in the couplings.This is the first time that these higher-order corrections are taken into account for a classically conformal model.It should be emphasised that at NLO, besides new contributions to the effective potential, there are also new effective operators that contribute to the kinetic part of the action, described by Z 3 .This contribution can straightforwardly be included in the EFT framework, but is absent in the typical daisy-resummed approach.This is a serious shortcoming of the standard approach, as we find that the effect of Z 3 is significant.The significance of the Z 3 -factor modifying the kinetic term as well the presence of the scale-shifting [39] fields suggests that the derivative expansion of the effective action might not be fully reliable.This should be studied by computing the functional determinant in the action prefactor, using the methodology developed in refs.[74,77,78,88].
We formulated predictions based on the theory at NLO, by implementing a gaugeinvariant and 3D RG-scale-invariant approach in section 4.4, based on [39].In such a setting, we performed a scan of the full parameter space of a BSM model, which has not yet been done in the context of GW production. 22We show in section 5.2 that the differences in the percolation temperature T p between the LO and NLO approach become as large as O(1), with the largest corrections occurring in the large-mass, small-coupling corner.The differences in the predicted length scale of the transition R * H * are more moderate, but also more relevant for the GW prediction, which depends strongly on R * H * .
We thoroughly study the scale dependence of the NLO predictions and compare them to the LO results in section 5.4.1.We find that the dependence on the RG-scale of the 4D theory, which is a measure for inaccuracy associated to missing higher-order corrections (see e.g.[28,29]), becomes reduced in the NLO prediction compared to the LO prediction, indicating that higher-order corrections indeed are required to reduce this source of uncertainty.
To sum up, we have demonstrated that higher-order corrections in the computation of thermal parameters in theories with supercooling can and should be included.We have found that the higher-order corrections have a significant effect on the GW signal, and that further studies into the contribution of the scale-shifting nature of the gauge fields is required.and the masses by  with γ E the Euler-Mascheroni's constant and A the Glaisher's constant.The effective masses depend on the effective couplings between the scalar and the several temporal gauge modes, given by λ V L5 = g 2 1 T λ hφ 8π 2 , (B.11) where λ V L5 , λ V L6 denote couplings between φ and the SM gauge fields, and λ V L8 the coupling between φ and the dark gauge field.

3 Figure 1 :
Figure 1: Schematic mass hierarchies we encounter along the bounce solution as a function of the radial coordinate.

Figure 2 :
Figure 2: Comparison of different approximations to the effective potential at low field values, around the barrier (left) and at large field values, around the minimum (right) at T n = 14.59GeV, for g X = 0.8, M X = 10 TeV.

5 Figure 3 :
Figure 3: Illustration of the validity of the high-temperature expansion at T n = 14.59GeV for g X = 0.8, M X = 10 TeV, for field-values in the vicinity of the escape point (vertical grey line) and beyond (left panel), and for larger field values (right panel).

Figure 4 :
Figure 4: Upper panel: Percolation temperature (left) and the transition length scale (right) obtained from the NLO potential with µ 4 = πT and µ 3 = m 2 X,3 (v 3 ).Lower panel: Percolation temperature (left) and the transition length scale (right) obtained from the LO potential with κ in eq.(3.7) set to κ = π to match the choice of µ 4 for the NLO potential.The grey regions are excluded for reasons explained in the main text.

Figure 5 :
Figure 5: Absolute value of the differences in the predictions for T p (left) and R * H * (right) between the NLO and LO potentials, normalised by the NLO quantities.The excluded regions are the same as in figure 4 (lower panel).The white dashed and dotted lines indicate the excluded regions obtained at NLO (as presented in the upper panel of figure 4).

Figure 6 :
Figure 6: The efficiency factor for sound waves κ sw computed from the LO potential.White lines correspond to excluded regions of parameter space obtained in the NLO setup.The red lines correspond to contours of κ sw = 0.99 (upper) and κ sw = 0.01 (lower), obtained with the NLO potential.

Figure 7 :
Figure 7: The values of SNR predicted at NLO (left panel), and at LO (right panel).

Figure 8 :
Figure 8: Left panel: Absolute value of the differences in T p obtained from the NLO action with µ 4 = πT and µ 4 = 2πT , normalised by the result with µ 4 = πT .Right panel: Absolute value of the differences in T p obtained from the LO action with κ = π and κ = 2π, normalised by the result with κ = π.

Figure 9 :
Figure 9: Different contributions to the integrand of the NLO action S 2 of eq.(4.19) for g X = 0.8, M X = 10 TeV, and T n = 14.18GeV.The blue solid line represents the full NLO contribution, the red dot-dashed line the contribution from the NLO contribution to the potential and the dashed green line the contribution from the kinetic term with Z 3 .