Thermodynamics of a two-step electroweak phase transition

New field content beyond that of the Standard Model of particle physics can alter the thermal history of electroweak symmetry breaking in the early universe. In particular, the symmetry breaking may have occurred through a sequence of successive phase transitions. We study the thermodynamics of such scenario in a real triplet extension of the Standard Model, using nonperturbative lattice simulations. Two-step electroweak phase transition is found to occur in a narrow region of allowed parameter space with the second transition always being first order. The first transition into the phase of non-vanishing triplet vacuum expectation value is first order in a non-negligible portion of the two-step parameter space. A comparison with 2-loop perturbative calculation is provided and significant discrepancies with the nonperturbative results are identified.


I. INTRODUCTION
In the Standard Model (SM) of particle physics, electroweak (EW) gauge symmetry is spontaneously broken by the vacuum-expectation value (VEV) of the Higgs field. Thermal corrections to the Higgs potential restore this symmetry in the early universe. For the physical Higgs mass this transition is a smooth crossover rather than a true phase transition [1][2][3], i.e., there is no distinction between the symmetric and broken "phases". In many beyond the Standard Model (BSM) scenarios, the introduction of additional scalar fields can result in a scalar potential having vastly different thermal behavior from that of the SM. In particular, these extensions may yield a bona fide electroweak phase transition (EWPT) that is first order, with cosmological consequences that include conditions needed to generate the cosmic matterantimatter asymmetry through electroweak baryogenesis (EWBG) [4][5][6] and production of gravitational waves (GW). A conclusive test of this possibility could result from present and future high energy collider experiments [7] and GW probes [8][9][10].
An extended scalar potential may admit a richer thermal history than in the SM. The new fields may have phase transitions of their own, and the universe may undergo several symmetry-breaking transitions before settling down to the present EW vacuum. While such a thermal history would be interesting in itself, multi-step EW symmetry breaking could have important implications for cosmology. Specifically, EWBG could be realized in a sequence of symmetry-breaking transitions around the EW scale [11][12][13][14][15]. This setup also leads naturally to a strong first-order transition into the final EW phase through a tree-level potential barrier. Furthermore, a non-minimal pattern of EW symmetry breaking can produce topological solitons, such as monopoles and domain walls, with potentially interesting properties. Such defects are absent in the SM, but are generic in grand unified theories [16]; many analogues also exist in condensed matter systems [17].
The simplest extension of the SM scalar sector admitting distinct phases of broken EW symmetry in the early universe is the real triplet model with three BSM degrees of freedom, collectively denoted by Σ. In the resulting "ΣSM", EW symmetry breaking may occur directly in a single step from the unbroken phase O to the Higgs phase φ, or in two steps, O → Σ → φ, where EW symmetry is broken in both the Σ and φ phases. A delineation of the model parameters leading to either possibility is given in the perturbative analysis in [12]. Analogous studies in other models containing new scalars either charged or neutral under the SM gauge symmetries indicate that multi-step transitions may arise generically [11,14,[18][19][20][21][22][23][24][25]. Thus, a more thorough investigation of the thermal history and phase diagram of the ΣSM is well-motivated.
A robust determination of the phase diagram is a nontrivial task even for theories that are weakly coupled at zero temperature. The EWPT is driven by infrared (IR) bosonic fields, the Matsubara zero modes, whose mutual interactions are boosted by Bose enhancement. This results in a poor convergence of perturbation theory and ultimately renders the momentum scale ∼ g 2 T nonperturbative, g being a gauge coupling [26]. This problem affects gauge bosons in the symmetric high-temperature phase and scalar fields near a phase transition where their correlation lengths can grow large. Indeed, perturbation theory incorrectly predicts a first-order EWPT in the minimal SM. There is no a priori reason to trust the perturbative description in BSM settings either, unless one is interested solely in properties of the Higgs phase, where the VEV provides a perturbative mass for most excitations. Large couplings in the scalar sector may further aggravate the IR problem [27,28].
For the EW theory, a solution to the IR problem is known: the thermodynamics are well described by a three-dimensional (3d) effective field theory (EFT) for which nonperturbative lattice simulations can be carried out [29][30][31]. This "dimensional reduction" amounts to perturbatively integrating out nonzero Matsubara modes, and the resulting theory describes thermal fluctuations of the bosonic zero modes.
Here, we report on a nonperturbative study of the ΣSM using the 3d EFT. The results are used to obtain a realistic picture of the two-step EWPT scenario. We also assess the performance of the perturbative treatment in light of our nonperturbative results.

II. MODEL
The color neutral scalar field Σ carries no hypercharge, transforms under the adjoint representation of SU(2) L , and does not couple to SM fermions. For simplicity, we further require invariance under the Z 2 transformation Σ → −Σ, which allows for the VEV v Σ to vanish at T = 0. Doing so ensures consistency with bounds on the EW ρ parameter while enabling the neutral field Σ 0 to contribute to the dark matter relic density [32,33]. Recent studies of the corresponding collider and dark matter phenomenology appear in [34,35]. The most general, renormalizable scalar potential then reads where a = 1, 2, 3 is the adjoint index, with √ 2Σ ± = Σ 1 ∓ iΣ 2 and Σ 0 = Σ 3 .
For µ 2 φ > 0, the potential has a symmetry-breaking minimum in the Higgs direction, φ † φ = 1 2 v 2 with v Σ = 0. This corresponds to the standard EW minimum with three BSM excitations from the Σ field, whose masses are degenerate at tree level [32]. Following [36], we relate the Lagrangian parameters to EW observables through polemass renormalization at one-loop level, taking the mass M Σ of Σ 0 as an input parameter. We treat the couplings a 2 and b 4 as input parameters directly at the MS scale M Z . The one-loop correction is necessary to match the accuracy of our EFT construction below.
If µ 2 Σ > 0, a second minimum of V (φ, Σ) appears in the Σ direction, with v = 0. Physics in this Σ vacuum resembles that of the broken phase of the SU(2) Georgi-Glashow model [37]: SU(2) breaks to a U(1) gauge group distinct from that of the usual electromagnetic interaction.
In the Σ vacuum phase, the system admits 't Hooft-Polyakov monopole excitations [38][39][40]. These are topological soliton solutions of the field equations, carrying a magnetic charge under the residual U(1) gauge group. When the system crosses to the Σ vacuum, these monopoles can freeze-out as a result of existing longwavelength thermal fluctuations [41]. They may also play a role in the dynamics of the finite-temperature phase transition itself [42].
Thermal corrections can modify the vacuum structure. In the high-T limit, the leading effect is a Tdependent reduction of the squared mass parameters: Here Π φ,Σ are O(g 2 ) constants, where for notational convenience g 2 denotes a general quartic coupling. The thermal correction turns µ 2 φ negative at T φ ∼ 100 GeV, relaxing the Higgs VEV to zero. Two-step EWSB occurs if the thermal corrections drive µ 2 Σ negative at a higher temperature T Σ > T φ . The universe then resides in the symmetric phase O at high temperatures before transitioning into the Σ phase (O → Σ) at T Σ , followed by another phase transition (Σ → φ) into the final Higgs phase at T φ . The presence of a tree-level saddle point separating the φ and Σ minima suggests a first-order transition in the second stage.

III. EFFECTIVE THEORY AT HIGH TEMPERATURE
We derive the 3d EFT in the imaginary time formalism by integrating out modes with a nonzero Matsubara frequency, including all fermions, leading to the Euclidean space Lagrangian: Here F a ij is the SU(2) L field strength tensor. Thermal corrections from the hard scale πT are included in the barred parameters, whose matching was worked out to O(g 4 ) accuracy in [36] and includes corrections from temporal components of the gauge fields that generate a Debye screening mass and can be integrated out [31,36]. These couplings in (2) are dimensionful, and the fields are scaled by T −1/2 . We have neglected the U(1) Y gauge field and the SU(3) C sector as they have only a small effect on the EWPT [43] and do not couple to Σ.
The EFT is formally valid in the high-T limit m πT . By construction, its region of validity overlaps with that of the consistent daisy resummation of [44] as required to correctly describe physics at the "soft" scale gT . The EFT systematically includes these corrections.
To probe the parameter space for a two-step EWPT, we have scanned the parameters using the effective potential V eff calculated to two-loop order in the EFT. Evolution of the different minima is tracked using the gaugeinvariant approach described in [45,46]; details of the  Colored regions correspond to different types of EW symmetry-breaking transitions. The allowed parameter space is dominated by direct transitions into the Higgs phase (regions IV and V). Regions II and III lead to a two-step symmetry-breaking history distinguished by whether the O → Σ transition is a crossover; a first order EWPT; or a second-order EWPT (a line somewhere in the grey region). Both O → φ and O → Σ transitions grow stronger as the quartic portal coupling a 2 increases. In region I the EW minimum is not the global minimum at zero temperature, according to the one-loop, T = 0 effective potential. Our lattice benchmarks are marked with a cross; points BM1 and BM2 are discussed in detail below. calculation can be found in the Appendix. Two-step transitions occur in a narrow band separating the parameter space of one-step EWPTs (O → φ) from that where the EW minimum is metastable at T = 0. For b 4 = 0.25, this is illustrated in Fig. 1.
In a vast region of parameter space, the EWPT is driven solely by the Higgs doublet, which becomes parametrically light near the critical temperature,μ 2 φ ∼ (g 2 T ) 2 , due to a cancellation between vacuum and thermal masses. This allows us to integrate out Σ as an UV mode near the critical temperature T c , resulting in a simpler EFT for which the nonperturbative phase diagram is known [1,47]. This approach was taken in [36] to identify where the O → φ transition is a crossover: region V of Fig. 1. Deep in region IV, integrating out Σ is no longer justified, but our simulations verify that the transition remains first order here. The line separating regions IV and V corresponds to second order transitions. Its location is only accurate within ∼ 10% due to neglect of higher-dimensional operators [36].
It is interesting to ask where in the two-step EWPT region the O → Σ transition is first order. As in the SM case, perturbation theory does not provide reliable guidance; genuine nonperturbative input is needed. Qualitatively, at temperatures near a O → Σ transition we expect the IR physics to resemble that of a Georgi-Glashow type theory containing just gauge fields and Σ. The corresponding phase transition terminates at a finite value of the scalar self coupling [42]. Our simulations confirm this expectation: the O → Σ transition is crossover in region II; first order in region III; and terminates somewhere in the grey region between. We have not attempted a more precise determination of the end line. The IR behavior also suggests that the O → Σ transition grows stronger at small b 4 , which we have verified with simulations using b 4 = 0.15, 0.20. However, the two-step region itself becomes narrower due to a decrease in the Σ minimum vacuum energy, which goes as at tree level. There is no two-step EWPT if the T = 0 potential is deeper in the Σ direction than in the Higgs minimum (region I).

IV. SIMULATIONS
Simulations in the full ΣSM are not practical due to the chirally coupled fermions. A systematic method for implementing the fermionic corrections (which are significant) is provided by the dimensionally-reduced EFT (2). To discretize it, we employ the (unimproved) Wilson action for the gauge links and couple these to the scalars through gauge-invariant hopping terms. Parameters in the lattice action are related to the continuum parameters in Eq. (2) by expressions given in [48]. These relations become exact in the continuum limit as a consequence of super-renormalizability of the 3d EFT.
In lattice simulations, we determine probability distributions of gauge-invariant operators by generating field configurations in the canonical ensemble. For the EWPT, the observables of interest are scalar condensates, particularly φ † φ and Σ a Σ a , whose probability distributions in a first order transition develop a two-peak structure. The peaks correspond to the bulk phases and have equal integrated probabilities at T c [47].
In the region separating the bulk phases, the ensemble is dominated by mixed-phase configurations where the two phases exist simultaneously on the lattice [49]. The phase interface carries free energy proportional to its surface area, and the probability of tunneling between phases is thus exponentially suppressed. On large lattices, this makes it difficult to obtain the probability distributions using conventional update algorithms for the canonical ensemble. Instead, we apply multicanonical simulations to boost the probabilities of the mixed configurations relative to the bulk phases [50]. The ensemble is modified by a suitable weight function W as where Φ multi is typically an order parameter-like quantity that distinguishes the phases. The canonical distributions are then ob- tained by reweighting the measurements [51]. While W itself can be calculated recursively [52,53], the efficiency of multicanonical simulations depends on the choice of Φ multi . For O → Σ transitions, we choose the volume average of Σ a Σ a , and the simulations proceed analogously to those of Refs. [42,47]. Consistent with perturbation theory we found that, for all the cases we studied (crosses in Fig. 1), the Σ → φ stage is a first-order transition with strong suppression of the mixed configurations. We have not found a simple choice of Φ multi that would efficiently take the system both ways between the two broken phases. In either phase, one of the scalar condensates develops large bulk fluctuations, and an even larger fluctuation is required to start the tunneling process. Instead, we determine T c by restricting the simulation to sample the mixed-phase configurations only. At T c , neither phase is preferred over the other, and probability distributions of order parameters in the allowed range become approximately flat.
The simulation results carry mild dependence on lattice volume and spacing a, but extrapolations V → ∞ and a → 0 can be taken in a controlled fashion [47,54,55]. For the first-order transitions studied here, O(a) errors appear negligible for 4/(aḡ 2 ) ≥ 20, with both T c and condensate values changing by less than 1% if a is decreased. Volume dependence appears to be even smaller, suggesting that our finite-size effects are well under control. Below we quote results from only the largest lattices.
Our code for simulating the SU(2) theory with fundamental and adjoint Higgses has been cross-checked by reproducing histograms in Refs. [42,47]. We employ con-ventional heatbath updates for the gauge links [56] and a mixture of Metropolis and over-relaxation updates [47] for the scalars.

V. RESULTS & DISCUSSION
The condensates require additive renormalization, but their discontinuities across a phase transition are renormalization group invariant and directly related to the latent heat L [48,57], a physical quantity characterizing transition strength. Fig. 2 shows the condensate evolution for two benchmark (BM) points giving a two-step EWPT, together with perturbative estimates. In the high-T phase the condensates stay close to zero, while at low temperatures φ † φ obtains a large value. The existence of an intermediate Σ phase is clearly visible. The condensates can be negative because of the additive renormalization.
In BM1, the O → Σ stage is a crossover: we find no evidence of phase coexistence, ruling out a first-order transition. To investigate the possibility of a second order transition we studied finite-size scaling of the dimensionless Σ 2 susceptibility, where the subscript denotes volume averaging. As shown in Fig. 3, χ(Σ 2 ) peaks at T ≈ 142 GeV but converges to a finite value as V → ∞, consistent with crossover behavior. By contrast, for a second-order transition the susceptibility diverges with a critical exponent. The first tran- sition in BM2 is first order and, depending on the criterion for baryon number preservation within the Σ phase, could be strong enough to support two-step EWBG [12].
To assess the reliability of perturbation theory, we compare the nonperturbative results to those obtained from the two-loop V eff (solid lines in Fig. 2). In the gaugeinvariant treatment used here, the potential is minimized by expanding the VEVs around their tree-level values (see the Appendix). Near the O → Σ transition, this approach breaks down due to the absence of a small expansion parameter, and the potential encounters an IR divergence [45,57]. This is the reason for the spiking of Σ a Σ a /T in fig. 2. Consequently, the crossover in BM1 is not visible perturbatively.
Outside the temperature range of O → Σ transitions, perturbation theory provides some rough qualitative guidance but performs poorly in quantitatively describing both T c and the "strength" (condensate discontinuities). For the Σ → φ transition this finding is, perhaps, surprising, as the two minima are present already in the tree-level potential. Nevertheless, scalar loops can significantly alter the transition dynamics because of the large a 2 coupling necessary for a two-step EWPT. On dimensional grounds, the high-T expansion parameter for φ − Σ interactions is of the form a 2 T × (scalar mass) −1 . At the second transition, the masses are bounded from below by the non-zero VEVs but are still numerically small compared to a 2 T on both sides of the transition. Hence, a low-order perturbative description is not necessarily reliable. There may also be additional nonperturbative effects from magnetic monopoles that exist in the Σ phase.
After extrapolating V → ∞ and a → 0, the latent heat is L/T 4 c = 0.4109(2) for the Σ → φ transition in BM1; in BM2 the value for the first (second) transition is 0.151(2) (0.5895 (9)). Errors are the stastistical uncertainties. The perturbative values, where applicable, are smaller by 30% in BM1 and larger by 40% in BM2. The discrepancy is dominated by the error in T c . The two-loop potential is crucial for even a qualitative agreement with the nonperturbative results: at one loop, the jumps in condensates are more than 50% smaller than at two loops, while the temperatures differ only by a few percent.
Applicability of these results to the full 4d ΣSM depends on the overall accuracy of our 3d EFT. Dimensional reduction produces operators of dimension six (in 4d units) that we have neglected here. We anticipate that the operators c φ (φ † φ) 3 /T 2 and c Σ (Σ a Σ a ) 3 /T 2 yield the largest contribution, with a potentially significant effect in the presence of a non-vanishing condensate. Following [31], we estimate their effects on scalar VEVs at tree level. For T > 50 GeV, the operators cause relative shifts of less than 1% in the VEVs in both BM1 and BM2, suggesting that the performance of our dimensional reduction is comparable to the SM case, despite the relatively heavy scalar excitations in the Higgs phase; indeed, top quark contributions still dominate c φ .
Overall, our results for the ΣSM phase diagram (Fig. 1) validate the expectations from purely perturbative studies that the early universe could have undergone successive EWSB transitions. To our knowledge, this work provides the first non-perturbative demonstration of this possiblility. At the same time, a robust determination of the character of these transitions and a quantitative determination of their properties (T C , latent heat, and model parameter-dependence) requires a nonperturbative treatment. Looking ahead, we anticipate that future nonperturbative studies will be essential for obtaining dynamical properties (e.g., rates for nucleation [49], sphaleron transitions [59], and monopole-catalyzed processes [60]) necessary for a complete picture of the associated thermal history in the ΣSM and other extended scalar sector scenarios. In this context, we consider the present study as the first step in a exciting program aimed at building a rigorous understanding of nonminimal electroweak symmetry breaking.

Appendix
Appendix A: Gauge-invariant effective potential to two loops Here we collect details of the perturbative calculation that was used for comparison with the nonperturbative results in the main text. The goal is to compute thermal corrections to the effective potential V eff and extract from it values for T c , latent heat and the condensates, and we do this at two-loop level. The perturbative expansion of V eff in terms of quartic couplings has a peculiar structure at finite temperature, with fractional powers such as λ 3/2 appearing as a consequence of Debye screening. A consistent inclusion of these effects requires daisy resummation in the high-T approximation [44], but as discussed in the main text, it is easier to work directly in the 3d EFT given in Eq. (2), where these resummations are incorporated automatically. We take this approach, generalizing the calculation of [61] to include a background field for the triplet.
Parameters of the EFT are related to those in the full theory by matching relations presented in [36]. Here we simplify the notation by dropping the overline from the EFT parameters. In what follows, all parameters are assumed to be those of the 3d theory (2) and therefore temperature dependent. For completeness we also include the U (1) Y hypercharge field B i , so the covariant derivatives read For comparison with the nonperturbative results we have set g = 0, as the U (1) Y field is not included in our lattice simulations. We parametrize the scalars as where v and x are real background fields. The Euclidean Lagrangian (2) becomes Here L (2) 3d and L (I) 3d contain quadratic and interaction terms respectively. Terms linear in φ i or Σ i do not contribute to the effective potential, which is defined (at a finite volume V) through .
The symbolic measure Dφ denotes functional integration over all dynamical fields, and the expectation value is to be calculated perturbatively. As discussed in [46], the value of V eff in its minimum is guaranteed, by Nielsen identities, to be gauge invariant order-by-order in the loop-counting parameter . Expanding the potential and its minima as and generalizing the analysis of [45,46] to the case of two background fields gives the " expansion" .
All derivatives are to be evaluated at the tree-level minimum (v 0 , x 0 ). Note that corrections to the VEVs contribute only at O( 2 ). This form of V eff (v min , x min ) is gauge invariant, and we shall calculate it in Landau gauge ξ = 0. With this choice, ghost fields remain massless after symmetry breaking and decouple from Goldstone modes.
, while the O( ) part can be calculated by diagonalizing the quadratic Lagrangian in momentum space. In the V → ∞ limit, the result is the familiar Coleman-Weinberg correction in d = 3 − 2 Euclidean dimensions: Here the integral is finite in 3d, and the field-dependent masses read The O( 2 ) correction consists of the 2-loop potential evaluated at a tree-level minimum, as well as 1-loop corrections to locations of the minima. The latter is obtained from Eqs. (A8)-(A10), while the former requires computation of one-particle-irreducible vacuum diagrams depicted in Fig. 4. Including the minus sign from exp − 1 d 3 x L (I) 3d in the diagrammatic vertex rules, V 2 is given by minus the sum of diagrams in Fig. 4.
The calculation of V 2 (v, x) at general field values is somewhat complicated as one needs to introduce a fielddependent mixing angle for the neutral scalars. A simpler approach is to perform the computation directly at a tree-level minimum (v 0 , x 0 ), which is all that is needed for the O( 2 ) correction. In the case of ΣSM, there is then no mixing between the mass eigenstates of φ and Σ as guaranteed by the Z 2 symmetry. Consequently, the masses m 2 ± in Eq. (A12) reduce to m 2 3 and m 2 4 . Below we present results for the different diagram topologies at two loops, expressed in terms of master integrals, but emphasize that these results are not applicable if v 0 and x 0 are simultaneously non-vanishing. Contributions from the Σ field are collected in curly brackets.
(VVV) = 1 2 The loop integrals are defined, in dimensional regularization with MS scale Λ, as Some special cases of the vector "sunset" integrals have been calculated previously in [61]. In the presence of the hypercharge gauge field, the following generalizations are needed.
Many of the expressions above utilize integration-by-part techniques developed in [62,63] (for thermal sum-integrals, see [64]). We are grateful to Philipp Schicho for providing particularly simple expressions for the special cases of D V SS , D V V S and D V V V . The two-loop diagrams are UV divergent, but are regulated (apart from the vacuum divergence) by mass counterterms in the tree-level part (A4). These are given by g 4 − 5 16 g 4 − 9 8 g 2 g 2 + 3λ(3g 2 + g 2 ) − 12λ 2 − 3 2 a 2 2 + 6a 2 g 2 (A38) which were also obtained independently in Ref. [36]. Due to super-renormalizability, there are no further corrections to the counterterms at higher loop orders. Apart from contributions from the triplet and the hypercharge field, the two-loop expressions in Eqs. (A14)-(A21) agree with those given in [61] for an SU(2)-Higgs theory.
To study the phase structure, we evaluate V eff (v 0 , x 0 ) separately in the three phases, and varying the temperature (which is now encapsuled in the 3d parameters). Not all of the above minima exist simultaneously at a given temperature; this needs to be checked separately. The condition for T c is that the value of V eff in any two minima is degenerate, e.g. V Σ eff (T c ) = V φ eff (T c ) for Σ → φ transitions. Latent heat is calculated from the 3d potential (which has units GeV 3 ) as Finally, the scalar condensates are given by [57] As discussed in the main text and in Refs. [45,57], the two-loop potential constructed here is not useful for studying thermodynamic properties near the critical temperature for transitions out of the symmetric phase (v 0 , x 0 ) = (0, 0). The issue lies in Eq. (A6), which assumes that the true minimum is related to the tree-level one through small perturbations. This assumption breaks down at temperatures close to O → φ or O → Σ transitions, for which the tree-level condition for T c is that the thermally-corrected mass parameter vanishes, µ 2 φ (T c ) = 0 or µ 2 Σ (T c ) = 0. Given that the high-T expansion parameter is ∼ g 2 T × (mass scale) −1 , perturbation theory is unreliable near T c . In particular, there is an explicit divergence at two-loop order due to the vanishing scalar mass [45]. This problem does not arise for transitions between two broken phases (Σ → φ) as the tree-level masses need not vanish for such transitions.
One may hope to regulate the problem by giving up on the expansion of (v min , x min ) altogether and solve for the minimum of V eff "exactly", as is frequently done in the literature. This automatically incorporates higher-order corrections of the VEVs into the potential. The downside is that these corrections also include uncancelled gauge dependence, and estimating the effects of this residual gauge dependence on final results is not a well-defined endeavor. Even if the resulting potential is free of spurious IR divergences, there is still no guarantee that the perturbative description near T c is reliable.