On the validity of perturbative studies of the electroweak phase transition in the Two Higgs Doublet model

Abstract
 
 Making use of a dimensionally-reduced effective theory at high temperature, we perform a nonperturbative study of the electroweak phase transition in the Two Higgs Doublet model. We focus on two phenomenologically allowed points in the parameter space, carrying out dynamical lattice simulations to determine the equilibrium properties of the transition. We discuss the shortcomings of conventional perturbative approaches based on the resummed effective potential — regarding the insufficient handling of infrared resummation but also the need to account for corrections beyond 1-loop order in the presence of large scalar couplings — and demonstrate that greater accuracy can be achieved with perturbative methods within the effective theory. We find that in the presence of very large scalar couplings, strong phase transitions cannot be reliably studied with any of the methods.

Abstract: Making use of a dimensionally-reduced effective theory at high temperature, we perform a nonperturbative study of the electroweak phase transition in the Two Higgs Doublet model. We focus on two phenomenologically allowed points in the parameter space, carrying out dynamical lattice simulations to determine the equilibrium properties of the transition. We discuss the shortcomings of conventional perturbative approaches based on the resummed effective potential -regarding the insufficient handling of infrared resummation but also the need to account for corrections beyond 1-loop order in the presence of large scalar couplings -and demonstrate that greater accuracy can be achieved with perturbative methods within the effective theory. We find that in the presence of very large scalar couplings, strong phase transitions cannot be reliably studied with any of the methods.

Introduction
The great triumph of the Standard Model (SM) of particle physics was achieved in 2012 with the discovery of its last missing piece, the Higgs boson, by the ATLAS and CMS experiments at the Large Hadron Collider (LHC) [1,2]. Although the properties of the observed Higgs boson so far agree with the SM predictions, it may be just one member of an extended Higgs sector. Frameworks with non-minimal scalar sectors are amongst the best motivated beyond SM (BSM) scenarios, as they may provide solutions to many of the SM shortcomings, such as the origin of the observed baryon excess in the universe.
One of the most promising frameworks for producing this asymmetry is electroweak baryogenesis (EWBG), which produces the baryon excess during the electroweak phase transition (EWPT) at a temperature T ∼ 100 GeV. Although the SM contains all the required ingredients for EWBG [3][4][5], it is unable to explain the observed baryon excess due to its insufficient amount of CP violation [6][7][8][9][10] and the lack of a first-order EWPT. Nonperturbative lattice studies in the SM have revealed that the Higgs boson is too heavy to lead to a large potential barrier between the symmetric and broken phases, and the electroweak symmetry breaking in the SM is thus a smooth crossover [11][12][13][14][15].
The study of EWBG is well-motivated in BSM models with extended Higgs sectors, which allow for new sources of CP violation and could provide a strongly first order phase transition. Among the simplest non-minimal Higgs frameworks are the Two Higgs Doublet Models (2HDMs), where the scalar sector is extended with one scalar doublet that has the same quantum numbers as the SM Higgs doublet. Both CP conserving and CP violating 2HDM frameworks have been studied in detail in the literature [16][17][18][19][20][21][22][23][24][25][26][27].
A common feature of BSM models with strongly first-order EWPT is that the relevant new fields can be light and hence dynamically active during the phase transition. This setup potentially leads to a multi-step transition with a tree-level potential barrier between the intermediate minimum and the final Higgs phase [28][29][30][31][32][33]. Alternatively, radiative corrections from new fields strongly coupled to the Higgs boson can induce a large barrier between the origin and the Higgs phase and facilitate a strong single-step transition. In what follows, we will focus on the latter option and leave the discussion of multi-step phase transitions for future work. Strong phase transitions are interesting also because they can produce gravitational waves that may be observed in the near future [34,35]. For both baryogenesis and gravitational-wave predictions, precise knowledge of the equilibrium properties of the EWPT is crucial.
In the context of the EWPT, variations of 2HDMs have been considered where the phase transition is analyzed using the perturbative effective potential [36][37][38][39][40][41][42][43]. Generically, in these works, a strongly first order EWPT is achieved through scalar couplings of O (1) or larger, which raises concerns of the performance of perturbation theory already at zero temperature. Additionally, at finite temperatures perturbative expansions suffer from severe infrared (IR) divergences in the presence of massless bosons [44]. In particular, the symmetric Higgs phase is inherently nonperturbative and cannot be rigorously described by perturbative weak coupling methods, including the ring-improved perturbation expansions [45][46][47][48].
The IR problem can be overcome with lattice Monte Carlo simulations. However, it is not known how to implement lattice fermions with non-Abelian chiral gauge couplings, rendering simulations of the full electroweak sector of the theory impossible [49]. Hence, the predominant approach for EWPT simulations is to make use of dimensionally-reduced effective theories (EFTs) [50][51][52][53] (see however [14]). In short, the EFT is obtained by integrating out non-zero Matsubara modes including all fermions, which have effective masses of order 2πT in the heat bath and decouple from the long-distance physics governing the phase transition. The EFT is then effectively three dimensional (hereafter 3d EFT), simplifying both perturbative and nonperturbative computations 1 .
In the paper at hand, we present a state-of-art study of the equilibrium dynamics of the EWPT, focusing on the CP-conserving 2HDM for simplicity. We carry out simulations with two dynamical doublets, overcoming the limitations of the previous analysis of Ref. [59] where only a limited region of the parameter space could be studied nonperturbatively. Since nonperturbative methods are very time consuming, we limit our analysis to two phenomenologically-motivated benchmark (BM) points where we carry out nonperturbative simulations, and perform a thorough comparison with conventional perturbative approaches. We find that for the moderately strong transitions studied here, the nonperturbative effects from IR physics are small in comparison to inaccuracies arising from bad convergence of perturbation theory due to the large scalar couplings, even at 2-loop level. As very strong transitions are generally associated with even larger couplings, our results suggests that in these cases perturbation theory fails to qualitatively describe the EWPT, and purely nonperturbative methods are then required.
The remainder of the paper is organized as follows. In section 2 we review the scalar potential of the 2HDM and identify experimental and theoretical constraints applicable to our analysis. In section 3 we introduce the 3d EFT and discuss the basic ideas on which dimensional reduction is based. Our benchmark points for the lattice analysis are presented in section 4. Section 5 is devoted to describing the lattice simulations, while in section 6 we compare the perturbative and nonperturbative treatments of the model and justify the validity of our effective theory. Finally, in section 7 we draw our conclusions.

The Two Higgs Doublet model
We start by introducing the 2HDM. It consists of gauge and fermion sectors as in the SM, a scalar potential V (φ 1 , φ 2 ) for the two SU(2) L scalar doublets with hypercharges Y = 1 and their kinetic terms, as well as Yukawa interactions. We describe below the scalar potential and the structure of the Yukawa sector.

The scalar potential and the Yukawa sector
In general, models with multiple Higgs doublets which can couple to fermions are at the risk of introducing Flavour Changing Neutral Currents (FCNCs) at tree level, which are tightly constrained by experiment. In the case of the 2HDM, imposing a Z 2 symmetrywhich can be softly-broken -on the scalar potential and extending it to the fermion sector forbids these FCNCs. We therefore focus on a scalar potential of the form where µ 2 12 causes a soft violation of the Z 2 symmetry, and explicitly Z 2 -breaking terms of have been discarded. In general, µ 2 12 and λ 5 can be complex. We write the field composition of the two scalar doublets as where the vacuum-expectation values (VEVs) v 1 and v 2 could in principle be complex. The complex phases of λ 5 , µ 2 12 and the VEVs are connected to CP-violation in the scalar sector, relevant for baryon number violation during the EWPT. In the case a softly broken Z 2 symmetry, spontaneous CP violation can occur when Im(λ * 5 [µ 2 12 ] 2 ) = 0 [18,62] and there exist no basis in which λ 5 , µ 2 12 and the VEVs are real. An exact Z 2 symmetry forbids the soft-breaking term µ 2 12 , which in turn leads to a real λ 5 2 . However, CP-violating phases also contribute to the electric dipole moment (EDM) and are hence heavily constrained by the strong bounds, in particular, on electron EDM from the ACME collaboration [64]. As a result, baryogenesis does not seem feasible in a simple 2HDM framework. However, our goal is not to solve the full EWBG problem, but to study how accurately one can determine the main properties of the phase transition. Hence, in the paper at hand, we fix µ 2 12 , λ 5 as well as the VEVs to be real, and do not discuss CP violation further. In the absence of CP violation, the mass eigenstates {h, H, A, H ± } can be written in terms of the fields in Eq. (2.2) and mixing angles α, β as where s α = sin α, c α = cos α, and the angle β is related to the doublet VEVs by tan β = v 2 /v 1 . The Z 2 charge assignment to fermions classifies four independent types of Yukawa interactions which are known as Type-I, Type-II, Type-X and Type-Y 3 [65][66][67]. We shall study a Type-I 2HDM, where fermions are coupled only to the φ 2 doublet. Consequently, constraints from flavor physics are less stringent than in Type-II, where down-type quarks are coupled to φ 1 instead 4 .
Finally, if the vacuum respects the Z 2 symmetry the φ 1 doublet does not acquire a VEV, and the model is reduced to the Inert Doublet Model (IDM) [69][70][71], for which a detailed EWPT analysis using a full 2-loop resummed effective potential has been performed in [43].

Input parameters for the numerical analysis
For our analysis, we renormalize the theory in the MS scheme and treat tan β and α as input parameters, together with the softly Z 2 -breaking parameter µ 2 ≡ −µ 2 12 and the physical pole masses {M h , M H , M A , M H ± }. These are related to the Lagrangian parameters via a 1-loop renormalization procedure described in detail in Appendix A. To summarize, the renormalized parameters are solved by requiring that the poles of loop-corrected propagators match the pole masses. Input parameters for the scalar sector are thus where the electromagnetic fine-structure constant α EM essentially fixes the combination v 2 ≡ v 2 1 + v 2 2 at 1-loop level, and the gauge couplings g, g are fixed by loop corrections to the masses of W ± , Z bosons. The scheme-dependent parameters tan β, α and µ 2 are input at a fixed MS scale. However, we will also discuss the EWPT in the presence of an inert φ 1 (µ 2 = v 1 = 0). In this case, we input the MS couplings λ 1 and λ 345 /2 ≡ (λ 3 + λ 4 + λ 5 )/2, corresponding to self-interaction and portal coupling of the dark matter candidate H [69][70][71], instead of the mixing angles. In what follows, we identify h as the observed scalar with M h = 125.09 GeV.
The quantity cos(β − α) is important for 2HDM phenomenology as it controls the interaction strengths of the CP-even scalars h, H to electroweak gauge bosons. The case cos(β − α) = 0 corresponds to the alignment limit where h couples to SM particles exactly as the physical Higgs in the SM. In practice, 2HDMs are driven to the alignment limit by constraints from collider experiments [72][73][74]. Additional particles may introduce important radiative corrections to gauge boson propagators. Furthermore, electroweak precision measurements of the oblique parameters [75][76][77][78] are satisfied when the charged scalar H ± is close in mass with either H or A [78,79]. In our analysis, we consider the M H ± = M A case.
In the phase transition analysis, we take into account the top Yukawa coupling, y t , to φ 2 and neglect the Yukawa couplings of other fermions due to their subdominant contribution. As a result, our EWPT analysis is valid for 2HDMs of Type-I 5 . It is worth mentioning that the light Yukawa couplings are included in our renormalization procedure, where we assume Type-I Yukawa couplings, but have verified that their numerical effect on the self energies is negligible.

High-temperature dimensional reduction
The concept of dimensional reduction in thermal field theory is based on the observation that a quantum system in a heat bath possesses a natural scale hierarchy. In the Matsubara formalism, this hierarchy can be made explicit by Fourier expanding the fields with respect to the imaginary time. The theory can then be described in terms of 3d Fourier modes, where the modes with a non-zero Fourier frequency obtain a mass correction of order 2πT in the propagators. These non-zero Matsubara modes can be integrated out perturbatively in an IR safe manner, and the resulting 3d EFT for the IR-sensitive Matsubara zero modes can then be studied nonperturbatively on the lattice. In general, the EFT is purely bosonic due to the lack of fermionic zero-modes, and numerical simulations are straightforward.

Three-dimensional effective theory for the 2HDM
The dimensionally-reduced EFT for the CP-conserving 2HDM has been derived in Ref. [80], extending the previous derivations of Refs. [55,81], and has the form where the (three-dimensional) field-strength tensor F rs and the covariant derivative D r are understood to contain both the SU(2) and U(1) gauge fields with the gauge couplings denoted byḡ andḡ . We have suppressed the gluon and ghost terms which are irrelevant for the EWPT, as they do not directly couple to the scalars 6 . Having effectively integrated out the temporal direction, this theory is defined in a three-dimensional Euclidean space and contains only the zero Matsubara frequency modes of the original fields. We denote the fields with the subscript 3d to emphasize this fact. Furthermore, we absorb the factor 1/T multiplying the action into the field definitions, so that the fields have the dimension GeV 1/2 and all couplings have a positive mass dimension. Higher-order operators, such as (φ † 1 φ 1 ) 3 3d , have been dropped from Eq. (3.1); we will discuss these operators in section 6.4. This EFT is constructed perturbatively by matching the Green's functions in both theories at O(λ 2 i ) accuracy in the scalar couplings and O(g 4 ), O(y 4 t ) in the gauge and top Yukawa couplings. This corresponds to a 1-loop matching of four-point functions and 2-loop matching of the scalar two-point functions. We use high-T expansion in computation of the sum-integrals, which leads to additional contributions of order O(µ 2 i λ j ), O(µ 2 i g 2 ), O(µ 2 i y 2 t ) that are also contained in the matching relations. Let us note that the construction of the theory in Eq. (3.1) also involves integrating out the temporal components of the gauge fields, which generate effective masses of order gT due to Debye screening. This results in a small correction to the EFT parameters. A detailed derivation has been presented in Ref. [80]; in particular, see section 3.3 there for explicit matching relations for the couplings of the effective theory.
We emphasize that the 3d EFT approach is useful not only for lattice simulations, but also as a way of organizing perturbation theory. The reason is that thermal resummations beyond 1-loop order are automatically implemented in the parameter matching. Indeed, the renormalized massesμ 2 11 ,μ 2 22 are just the screened masses evaluated at 2-loop level [82]. Loop calculations are also simplified, as the 3d EFT is purely spatial and contains one mass scale less than the full theory (the temperature), and is furthermore super-renormalizable.

Lattice formulation of the effective theory
Our discretized action in three dimensions, corresponding to the 3d EFT in Eq. (3.1), reads where a is the lattice spacing and the dimensionless constant β G is given by In Eq. (3.2), U i (x) are the SU(2) gauge links and P ij is the standard Wilson plaquette.
Following previous lattice studies of the EWPT [11,52,[83][84][85], we have dropped the U(1) gauge field from the lattice action as its effect on the dynamics of the transition is small [13]. The masses m 2 11 , m 2 22 , m 2 12 and fields Φ 1 , Φ 2 in the lattice action are related to the corresponding continuum quantitiesμ 2 11 ,μ 2 22 ,μ 2 12 and (φ 1 ) 3d , (φ 2 ) 3d in the MS scheme by relations that can be found in Ref. [85] (see Appendix B therein). We emphasize that due to the super-renormalizability of the effective 3d theory, Eq. (3.1), these lattice-continuum relations are exact and not susceptible to perturbative errors [86,87]. Couplings in the 3d theory do not run, so their lattice-continuum relations are trivial. For the actual simulations, we find it convenient to make the fields and couplings dimensionless by scaling them with appropriate factors of T as in Ref. [85].

Choosing the benchmark points
Because of the computational effort required for lattice simulations, it is not possible to perform nonperturbative scans over the whole parameter space allowed by theoretical and experimental constraints. We therefore need to focus our analysis on a couple of points from which one can hope to draw more general conclusions about the performance of perturbation theory. Let us reiterate that in order to generate a potential barrier for a strong single-step EWPT, some of the Higgs sector couplings will necessarily have to be large. Unfortunately, in many strong-EWPT scenarios present in the literature, some of these couplings are so large that the convergence of perturbation theory is at best marginal (cf. section 6.4), and care needs to be taken when constructing the 3d EFT perturbatively.
In order to guarantee the accuracy of our 3d EFT, the couplings should be kept small enough so that loop corrections from the heavy Matsubara modes remain under control. Furthermore, the thermal scale hierarchy should be respected, so that all scalar degrees of freedom have to be lighter than 2πT in both phases near the critical temperature, where φ c ∼ T c for strong transitions.

2HDM scenarios to be studied on the lattice
With the above considerations in mind, we have chosen two phenomenologically-viable BM scenarios -described in Table 1 -where we expect the 3d EFT to accurately describe the EWPT. In addition to verifying boundedness from below of the scalar potential at 1-loop level (see section 6.2), we have checked that tree-level unitarity constraints [88,89] are satisfied and that the largest coupling λ 3 stays below 2π at scales relevant for the EWPT 7 (c.f. section 5.1). Although some of the vacuum masses in both benchmarks are heavy compared to the electroweak scale, the high-T expansion used in dimensional reduction can still be expected to converge well as the fields are lighter near the phase transition [43,93].  Table 1: Input parameters for our benchmark points. In BM1, the combination λ 345 ≡ λ 3 + λ 4 + λ 5 corresponds to a dark matter portal coupling in the IDM [39], and µ = −µ 2 12 represents the soft Z 2 -breaking parameter. Masses are assumed to be the physical pole masses, while the remaining parameters are input directly in the MS scheme at the initial renormalization scale Λ 0 .
BM1 is specific to the IDM and has been studied perturbatively in Refs. [39,43]. Our main motivation for studying this particular point on the lattice is to produce a quantitative comparison with the resummed 2-loop result of Ref. [43], where the 2-loop corrections to the effective potential were found to make the transition considerably weaker relative to a 1-loop calculation. To make this comparison, we have modified our renormalization procedure to match that of Ref. [43]; In particular, we neglect the U(1) sector by setting g = 0 in BM1, but have numerically verified that its inclusion in the 1-loop calculation has a negligible 7 The perturbativity bound λi < 2π is motivated by Refs. [90][91][92], where the breakdown of perturbation theory is demonstrated for couplings much smaller than the naive upper bound of 4π. In practice, the magnitude of our largest coupling λ3 in BM2 is roughly λ3 ∼ 4 at the scale of EWPT dynamics, with the other couplings being substantially smaller. effect on the renormalized parameters listed in Table 2. On the phenomenological side, BM1 provides a dark matter candidate H, which can constitute a fraction of the observed dark matter relic density [39].
BM2, in the softly Z 2 -breaking 2HDM, lies in the mass-hierarchy region where earlier studies based on the 1-loop effective potential report strong first-order phase transitions [37,41,94]. Our BM2 approaches the parameter-space points in which Refs. [34,40] predict gravitational-wave signatures in the sensitivity range of LISA. However, BM2 represents a more conservative EWPT scenario than what is shown in e.g. [40] with a large λ 3 , but where perturbation theory still converges reasonably well (modulo the usual IR problems), while also providing a moderately strong transition.
Phenomenologically, BM2 is motivated by possible collider signatures in the following processes, away from the alignment limit and in the small tan β region: the ratio of decay rates of h (the SM-like Higgs boson) to those of the h SM (the Higgs boson in the SM) in bb, τ + τ − , gg channels, as well as the signal strength of the SM-like Higgs boson in the τ + τ − decay mode 8 . Complementary to those, the hW + and HW + or AW + decays of the heavy charged Higgs produce interesting experimental signatures for this BM point [20]. Finally, the A → Zh channel acts as an extra probe of this BM point [38,95].
Experimental constraints on 2HDMs come from gauge bosons width [72], direct searches for charged scalars and lifetime of charged scalars [96], Higgs total decay width [73], Higgs invisible branching ratios and Higgs to γγ signal strength [97,98], A → Zh searches [99], direct searches for extra Higgs bosons at the LHC [100] and flavour constraints [101]. We have verified that our benchmarks satisfy all current experimental bounds arising from these sources.
In Table 2, we list the renormalized parameters, obtained from the input parameters, as described in Appendix A. We have chosen different input MS scales for BM1 and BM2 points. In BM1, the parameters in Table 2 are solved from the loop-corrected pole-mass conditions at scale M Z ≈ 91 GeV, in accordance with Ref. [43]. In BM2, however, we choose the initial scale to be the average of the pole masses, Λ 0 = (M h + M H + M A + 2M H ± )/5, in order to reduce the size of logarithmic corrections in the self energies. This choice for the input scale is justified by the numerical analysis of Ref. [102], where the scale dependence of different renormalization prescriptions is discussed. Table 2: Renormalized parameters corresponding to the input parameters in Table 1. The recipe for obtaining these is described in the main text and in Appendix A. In BM1, we have set the U(1) coupling to zero for the sake of comparison with the results of Ref. [43]. The SU(3) coupling g s is fixed at tree level.

Lattice simulations
The discretized theory, Eq. (3.2), is studied nonperturbatively by evaluating expectation values of quantum operators using Monte Carlo integration. Our simulations are performed with the same Monte Carlo code that was used in Refs. [11,84,85], and the practical procedure is similar to that of Ref. [85] where the Minimally Supersymmetric Standard Model (MSSM) was studied on the lattice. We emphasize that since gauge fixing is not needed on the lattice, the results obtained here are manifestly gauge invariant.

Obtaining lattice parameters
Starting from a set of given input parameters, our analysis proceeds as follows. In order to account for logarithmic corrections of thermal origin, we first run the parameters in Table 2 to the scale with γ being the Euler-Mascheroni constant, using 1-loop β functions found e.g. in Refs. [19,80]. Eq. (5.1) has the physical interpretation of corresponding to the average momentum of integration over the non-zero Matsubara modes [103]. We then apply the matching relations of Ref. [80] to obtain the parameters of the effective theory, Eq. (3.1), as functions of the temperature, which are converted to lattice parameters using the relations provided in Ref. [85]. When the RG scale is chosen as in Eq. (5.1), all logarithmic corrections in the matching relations vanish, apart from those corresponding to the (exact) RG running of the masses within the 3d theory [80,103]. RG scale in the effective theory is separate from that of the full theory and we fix it as Λ 3d = T , although the lattice results are insensitive to this choice as RG running is exactly contained in the lattice-continuum relations.

Finding the transition point on the lattice
From the simulations, we obtain the gauge-invariant expectation values of the operators in Eq. (3.2) for fixed volume V and β G which sets the lattice spacing. At the critical point, the metastability of the phases is so strong that normal simulation methods do not efficiently tunnel between the phases. Hence, we apply multicanonical simulations [104] to overcome the potential barrier suppressing mixed-phase configurations. In both of our BM points, the Φ 2 field (in lattice discretization) dominates the phase transition dynamics, while the other doublet Φ 1 is so heavy that its condensate, Φ † 1 Φ 1 , changes only slightly at the transition point. We can therefore treat the expectation value Φ † 2 Φ 2 as an effective order parameter that determines the transition point.
The composite operators are divergent in the ultraviolet (UV) and can hence obtain negative values when renormalization is applied. The behavior of the condensates is plotted in Fig. 1, where we have converted the lattice fields Φ † i Φ i to the respective quantities in the MS scheme, (φ † i φ i ) 3d , by substracting the lattice divergence [87]. We shall drop the field subscripts in the following discussion. In BM1, the change in the φ † 1 φ 1 condensate is a result of field fluctuations becoming more constrained due to the φ 2 field obtaining a VEV. In BM2, on the other hand, φ 1 will also develop a VEV due to the requirement that tan β is non-zero in the T = 0 vacuum. This change is compensated by decreased fluctuations around the new minimum, and as a result, the condensate φ † 1 φ 1 remains in practice constant in the transition.   : Gauge-invariant condensates of the two doublets as measured on the lattice with fixed volume and β G , converted to MS quantities using the relations in Ref. [87]. In both cases, the doublet φ 1 is heavy and almost inert at the phase transition. Instead, its fluctuations become more constrained as φ † 2 φ 2 changes due to the λ 3 term in the action. In BM2, φ † 1 φ 1 increases smoothly towards its T = 0 value as governed by our choice of tan β.
In a first-order phase transition, the probability distribution of φ † 2 φ 2 has a twopeak structure such as the one shown in Fig. 2, with the peaks corresponding to the symmetric and broken phases. The probability of field configurations between the peaks is strongly suppressed, and the separation of the peaks corresponds to a potential barrier that the system has to overcome in order for the phase transition to occur. At the critical temperature T c , the probability of finding the system in either phase is equal. Hence, our criterion for finding T c is that the integrated probability under the peaks in the histogram is the same. In practice, the simulations are carried out at a temperature close to T c , and the precise critical temperature is then conveniently found by reweighting [105] the multicanonical distribution to a temperature that minimizes the integrated probability difference of the two phases. Figure 2: Unnormalized probability distributions of the expectation value φ † 2 φ 2 3d in BM1, measured on the lattice with varying β G and converted to the corresponding continuum quantity. The histograms have been obtained by reweighting in T and minimizing the integrated probability difference between the peaks. As β G increases, the phases become more separated and thus the transition grows stronger. Configurations between the two phases are suppressed exponentially by the free energy carried by the phase interface, with the suppression increasing with the area of the interface.

Determining physical observables from the simulations
In both BM points, properties of the transition depend strongly on the lattice spacing, or β G , as well as on the volume to a lesser extent. Continuum results for the equilibrium quantities characterizing the phase transition are obtained by first extrapolating to an infinite volume, and taking the lattice spacing to zero afterwards (corresponding to 1/β G → 0). The simulated values of V and β G are listed in Table 3. Each simulation has been weighted by an appropriate multicanonical weight function and consists of 1.5 × 10 6 − 2.5 × 10 6 measurements, depending on the volume. Cylindrical lattices have been used for a precise measurement of the interface tension (see below).  Our results are collected in Table 3 along with statistical errors, obtained with jackknife sampling, related to the continuum extrapolations. The behavior of the extrapolations is qualitatively similar to the MSSM case of Ref. [85], however in our simulations, the latent heat and order parameter discontinuity contain substantial dependence on the lattice volume.
Critical temperature: For individual simulations with fixed volume and β G , the critical temperature is determined from equal-weight histograms as described above. T c in both BM points is insensitive to the volume, and the dependence on 1/β G ∝ a is linear over the entire range of the lattice spacings used. Extrapolations to continuum are shown in Fig. 4. Statistical errors from the Monte Carlo simulations -as well as those from reweightingare small for the critical temperature. Discontinuity in the order parameter: From probability distributions at the critical temperature, such as those shown in Fig. 2, it is straightforward to measure the discontinuity in doublet condensates. In both of our BM points, however, only the Φ 2 condensate can be used as an order parameter, related to the condensate in MS renormalization as described in Ref. [87]. Volume dependence of the dimensionless quantity 2∆ φ † 2 φ 2 3d /T is shown in Fig. 5. Unlike in the MSSM study of Ref. [85], there is a significant dependence on the lattice volume, and we find the dependence to be approximately linear in the dimensionless combination 1/(V T 3 ).
β G dependence of the infinite-volume results is plotted in Fig. 6. In BM1, a leastsquares quadratic extrapolation fits nicely to the data points, but we have also included a linear fit. The difference between the two extrapolations is roughly 3% in the continuum limit, with the quadratic fit resulting in a slightly larger error. We shall use the quadratic extrapolation from now on.
In BM2, however, extrapolation is more difficult due to the smaller range of β G used in the simulations. As a result, a quadratic polynomial does not fit the points well. We have therefore used a linear fit in BM2, but emphasize that missing higher-order terms may have a numerical impact. Reliably probing the effect of 1/β 2 G and higher terms would require the use of very small lattice spacings, which is computationally expensive due to the large lattices required. Given that the perturbative uncertainty from dimensional reduction and zero-temperature renormalization is possibly considerably large, of the order 20% as estimated in section 6.4, we shall be content with a linear extrapolation here.
Often, the quantity used for determining the strength of the EWPT is not the order parameter discontinuity ∆ φ † φ , but the discontinuity in the (gauge-dependent) doublet VEV    Figure 6: Continuum extrapolations of the order parameter discontinuity in the infinitevolume limit. Both quadratic (green) and linear (blue) fits are shown in BM1. A reliable higher-order polynomial fit in BM2 would require the use of very large β G -and consequently extremely large lattices -due to the very massive scalar degrees of freedom.
-or, in the case of multiple doublets, some combination of their VEVs -which can conveniently be measured from the effective potential. We define a gauge-invariant counterpart to the conventional φ c /T c as which is dimensionless due to the different normalization of fields in the effective theory.
Latent heat: By definition, the latent heat L is the discontinuity in the energy density of the system, and is a concrete physical quantity characterizing the strength of the phase transition. In terms of the partition function, which can be determined from changes in the expectation values of composite field operators. Following Ref. [85], we measure the above quantity directly on the lattice and extrapolate the results to the continuum. Unlike in the MSSM case of Ref. [85] where all couplings were very small, we find the temperature dependence of quartic couplings to be significant and hence include also the interaction-term expectation values in the evaluation. Schematically, where the ellipsis represent the remaining interaction terms. The above quantity is easily obtained by reweighting.
Extrapolations of the latent heat to infinite volume and zero lattice spacing are shown in Figs. 7 and 8. We find that the numerical behavior of the latent heat is very similar to that of the order-parameter discontinuity and have hence used the same fitting ansatzes as for ∆v/T . Again in BM1, we will use the quadratic fit as the final result. As pointed out already in Ref. [85], statistical errors are somewhat larger for the latent heat. This is a natural result as there are more condensates involved in the determination of L/T 4 .

Surface tension:
The final equilibrium quantity we measure is the tension σ of the phase boundary separating the symmetric and broken phases. The interface tension reduces the likelihood of mixed-phase configurations, visible in the probability distributions of Fig. 2 as a suppressed "valley" between the two phases. The suppression is proportional to exp(−σA/T ), with A being the area of the phase boundary, and can be measured from the probability distributions using the histogram method [106]. Specifically, the quantity where P max and P min denote the maximum and minimum probability distribution between the peaks, respectively, will tend to σ/T in the infinite-volume limit. In cylindrical lattices with L z L x = L y , the phase interface generally will form perpendicular to the z direction, as this configuration is energetically favored over other possibilities. As in Refs. [11,84,85], we apply a finite-volume scaling ansatz in order to reduce large-volume effects related to lattice geometry. For L x = L y , an appropriate ansatz is [106]    where G = 0 for cylindrical lattices. A periodic lattice will contain two interfaces, and it is assumed that their mutual interactions can be neglected in Eq. (5.7). In practice, this condition is fulfilled for long lattices where the interfaces form far enough from each other. Our lattices generally have L z ≈ 4L x and, in this regard, are more ideal than the lattice shapes previously used in EWPT simulations.   Figure 9: Extrapolation of the surface tension to an infinite surface area. Statistical errors are substantial in the β G = 32 case and could be improved with more simulations.
The extrapolations 1/(AT 2 ) → 0 for the dimensionless combination σ/T 3 are shown in Fig. 9. The area dependence in both BM points is linear for all of our lattice spacings, but the extrapolations for β G = 32 come with large statistical uncertainty. Improving the fits would require the use of very large lattices, being computationally expensive, and given the limited use of the surface tension in practical applications we have chosen not to pursue better accuracy here.   10 shows a continuum estimate of the surface tension. Due to the large uncertainty in the β G = 32 infinite-volume extrapolation, we again have chosen a linear fit. Somewhat surprisingly, we find the dependence on β G to be roughly as important as the dependence on the interface area. This is due to the heavy degrees of freedom which become dynamical only at small lattice spacings.

Comparison with perturbation theory
Having obtained the equilibrium characteristics of the EWPT nonperturbatively on the lattice, we now wish to evaluate the same quantities in perturbation theory and compare the results. We shall use the effective potential V eff (ϕ 1 , ϕ 2 ) for classical background fields ϕ i and assume that the background is only in the neutral component, i.e, the doublets φ i are shifted as The background fields modify mass eigenvalues and couplings of the theory. The fielddependent masses of the scalars are given by the eigenvalues of the mass matrix where V is the tree-level potential (Eq. (2.1)) and the indices i, j refer to the components of the doublets (see Eq. (2.2)). Correspondingly, gauge boson masses are obtained by diagonalizing Finally, the field-dependent mass of the top quark in a Type I 2HDM is and we neglect other fermions from the phase-transition analysis due to their small couplings to the scalars. At the critical temperature, the loop-corrected effective potential will have a symmetrybreaking minimum that is degenerate with the symmetric minimum at the origin, and the strength of the phase transition can be determined from the potential barrier separating the two minima. However, it is well known that this simple approach will result in explicit gauge dependence of the location of the minimum [107,108]. As such, the critical temperature and consequently all other thermodynamical quantities calculated from the potential are not gauge invariant. A solution to this gauge problem is to compute V eff for the composite operators φ † i φ i rather than for ϕ i , and the thermodynamical quantities can then be obtained from the potential in a gauge-invariant manner [109].
In Ref. [43] -where the Feynman-t'Hooft gauge was used -the authors argued that ambiguities related to gauge dependence are overshadowed by the uncertainty related to higher-order corrections from the large scalar couplings. We shall therefore be content with the practical approach described at the beginning of this section, as is frequently done in the literature, and use Landau gauge for the perturbative calculations.

Perturbative calculation in the effective theory
We start with the effective potential constructed within the 3d EFT, Eq. (3.1). This calculation is simpler than the conventional V eff in the full theory, both conceptionally and computationally, as thermal corrections have already been accounted for in the dimensional reduction procedure. As mentioned briefly in section 3.1, this also includes thermal resummations beyond 1-loop order.
The 3d V eff allows for a direct comparison with the results obtained from lattice simulations that is not affected by possible uncertainties related to dimensional reduction. In particular, the magnitude of nonperturbative effects related to the "ultrasoft" scale g 2 T , for which no resummation is possible, can be estimated by comparing the 3d V eff to the lattice results. As RG running in 3d starts only at 2-loop level, it is desirable to calculate the 3d V eff to two loops. We carry out the calculation in d = 3−2 spatial dimensions using the MS scheme. Since the U(1) subgroup has been left out from the lattice simulations, we choose to drop its contributions to the 3d V eff as well. However, we have also performed the analysis with the full U(1) contributions included at 2-loop level and verified that their effect on the phase transition is small in comparison to systematic uncertainties in the calculation.
For a 3d EFT containing only one Higgs doublet, the 3d V eff and a list of the relevant integrals have been presented in Ref. [103]. We extend this calculation to our EFT with two doublets. Having integrated out fermions already in the dimensional reduction, the 1-loop correction to the 3d V eff is given by the bosonic zero modes as where the sum runs over all scalar eigenstates, and the UV-finite integral J 3 is given by We emphasize that all parameters entering the 3d V eff are those of the 3d EFT, Eq. (3.1), which themselves are functions of the renormalized parameters of the full theory and the temperature, as dictated by dimensional reduction. Consequently, the masses in Eq. (6.5) are the thermally-screened masses. Note that when the U(1) sector is neglected (ḡ = 0), we have m Z = m W . In the Landau gauge, the SU(2) ghosts are massless and do not enter the 1-loop corrections. The 2-loop correction is obtained from vacuum diagrams shown in Fig. 11. Although the calculation is standard, and the integrals can be found in Ref. [103], the number of diagrams is large due to the many scalar fields present in our theory. For this reason, we choose not to give the 2-loop contribution in an explicit form, and shall simply present the numerical results obtained with the full 2-loop potential in section 6.3.
UV divergent contributions arise at 2-loop level, which can be cancelled by introducing mass counterterms in the tree-level potential. The divergence matches the counterterms previously obtained in an independent calculation in Ref. [80], which serves as a cross check of our 3d V eff . As a result of super-renormalizibility, the RG running of the mass parameters with respect to the MS scale Λ 3 in 3d is defined exactly by the 2-loop counterterms, while the couplings remain RG invariant.

1-loop effective potential in the full theory
Next, we discuss the predominant approach for studying the EWPT perturbatively, namely the 1-loop resummed effective potential calculated in the full theory. It is given by where V CW is the T = 0 1-loop Colemann-Weinberg correction to the tree-level potential V , and V T is the 1-loop finite temperature part. The Colemann-Weinberg part is [110] where the sum runs over the particle content of the model, N j is the internal number of degrees of freedom which is positive for bosons and negative for fermions, m j is the field dependent mass, Λ is the renormalization scale, and C j = 5/6 for gauge bosons and C j = 3/2 for fermions and scalars. The unimproved thermal part is given by the integral This expression can be improved by accounting for the thermal dispersion relations of the (quasi)particles. Originally Parwani [46] suggested to do this by replacing the field dependent masses of bosons by their thermal masses, m j → m j (T ) in the whole 1-loop part of the effective potential. This approach is not self-consistent however. It leads to T -dependent divergences at higher orders and the choice of the thermal mass is ambiguous [43]. In a more consistent approach by Arnold and Espinosa [47], one introduces thermal masses only in the cubic terms. This corresponds to screening only the IR-sensitive zero modes and results in the following ring-improved potential On the other hand, ring improvement relies on the high-T expansion, which fails when the bosonic zero modes are also heavy. For this reason, the Parwani resummation allows a smoother continuation to the nonrelativistic limit in theories which contain heavy degrees of freedom [36,43]. We shall perform our numerical analysis using both of these resummation methods.
Fermions and transverse gauge boson modes do not receive thermal corrections. For the longitudinal gauge bosons, the thermal masses are obtained by diagonalizing M 2 g + δM 2 g , where 11) and the thermal scalar boson masses are obtained from M 2 + δM 2 , where with c k = 1 for bosons and c k = −1/2 for fermions.

Numerical results in perturbation theory
The condition that the symmetry-breaking minimum becomes degenerate with the minimum at the origin determines the critical temperature T c and the critical field value φ c ≡ ϕ 2 1,c + ϕ 2 2,c . Discontinuity in the quantity φ c /T c then corresponds roughly to the order parameter discontinuity obtained from lattice simulations (Eq. (5.3)), aside from the ambiguities related to gauge fixing. In the 3d EFT where the fields are scaled to mass dimension GeV 1/2 , the corresponding quantity is φ 3d c / √ T c . In the thermodynamic limit, the value of V eff in its minimum coincides with the grand canonical free energy density. Hence, the latent heat can be obtained from the effective potential as In the 3d analysis, the factor 1/T multiplying the 3d action has been absorbed into the definition of V 3d eff , and so the latent heat is obtained from the 3d effective potential as (6.14) For simplicity, we shall not compute the surface tension perturbatively.
The effective potential bears and explicit dependence on the RG scale Λ. While the full effective action is Λ-independent, in perturbative expansions there is always an uncertainty related to the scale variation, higher in order than the one under consideration. Due to the large scalar couplings present in our analysis, this ambiguity in the choice of Λ is a significant source of uncertainty. We demonstrate this by varying the RG scale along a range of mass scales with dominant contributions to the effective potential. Clearly the most dangerous logarithms are those proportional to the scalar couplings. In the full theory, contributions of the form λ 2 i ln m 2 /Λ 2 , where m denotes a scalar mass, arise from the T = 0 loop corrections, Eq. (6.8). However, thermal fluctuations generate additional logarithms of the type ln(Λ/(πT )), making it difficult to ensure that all logarithmic corrections stay small simultaneously. We choose to vary Λ between 0.5πT and 1.5πT .
The situation is different in the 3d EFT, where the 3d V eff contains only logarithms of the type ln m 2 3 /Λ 2 3 , which furthermore arise only at 2-loop level at the earliest. Having integrated out the mass scale πT , we vary the RG scale of the 3d EFT, Λ 3 , between 0.5T and 2T . However, scale variations in the full theory affect the parameters of the EFT via the matching relations obtained from dimensional reduction. Although this uncertainty is less severe than the corresponding ambiguity in the full V eff as only thermal logarithms appear in dimensional reduction. For the sake of having a one-to-one comparison of perturbative and nonperturbative analyses in the 3d EFT, we fix the RG scale of dimensional reduction as in Eq. (5.1), but discuss variations of this scale further in section 6.4.
Turning to numerical analysis, we present our findings in Fig. 12 and Table 4. In Fig. 12, we have plotted the global minimum of the potential as a function of the temperature with different resummation implementations, as well as in the 3d EFT. For each scenario, the RG scale has been varied as described above and the results for the smallest and largest scale are plotted in Fig. 12. The colored bands depict the uncertainty related to the scale sensitivity. In all plots, the parameters given in Table 2 have been run to the final scale using 1-loop β functions in the full theory. In the 3d analysis, the corresponding parameters in the 3d EFT are first obtained from dimensional reduction, and then run to the final 3d scale using exact RG evolution.
For both Parwani and Arnold-Espinosa resummations in the full theory, there is significant uncertainty in determination of critical temperature, with Parwani resummation giving a considerably smaller T c . On the contrary, the magnitude of the jump φ c /T c is not very sensitive to scale variations. In BM2, scale variations lead to much larger uncertainty than in BM1 as the couplings and masses are larger. In 3d perturbation theory, scale uncertainties are significantly less alarming, which is not surprising as the 3d EFT is super-renormalizable and the 3d V eff is evaluated to two loops.
Thermodynamic quantities obtained from the perturbative calculations are collected in Table 4, along with the nonperturbative lattice results from section 5.3. Comparing first the perturbative and nonperturbative results within the 3d EFT, we see that in both BM points, the 3d V eff describes the EWPT quite well. The strength of the transition is slightly underestimated by the perturbative analysis in BM2 and there is a ∼ 5% discrepancy in the critical temperature. Qualitatively the behavior is similar to that of the MSSM case [85,111].   Table 4: Comparison of thermodynamic quantities as obtained either with the EFT approach, or with the resummed 1-loop potential in the full theory. Error bars indicate sensitivity to the RG-scale as described in the text. For the lattice results we show the statistical uncertainty. Numbers in the last column correspond to central values without scale uncertainties.
We emphasize that apart from perturbative corrections beyond two loops, the main difference between the 3d perturbative and nonperturbative approaches is the handling of the IR sensitive fields in the symmetric phase. As such, we conclude that these nonperturbative IR effects are, in fact, already suppressed for the marginally strong phase transitions considered here. This is reassuring, and suggests that the EWPT can be studied reliably with the relatively simple 2-loop 3d V eff , at least in the cases considered here.
Moving on, we find that the situation is not as good for the 1-loop potential in the full theory. In particular, this approach overestimates T c by a large margin compared to the analysis in the 3d EFT. This is because the critical temperature is particularly sensitive to thermal mass corrections and resummations, which are incorporated at 2-loop level in the 3d EFT. With additionally the scale uncertainty in T c being about 10 to 25 percent, we conclude that at least in our studied BM points, the full 1-loop V eff with thermal resummations is not a reliable tool for determining the temperature scale of the EWPT. On the other hand, L and φ c /T qualitatively match the lattice results, but the dimensionless combination L/T 4 c , important for gravitational-wave predictions, is underestimated due to the overly large T c . BM1 has previously been studied at the full 2-loop level without dimensional reduction in Ref. [43], where a Parwani-type resummation was applied. According to their Table 4, the 2-loop corrections to the potential make the transition slightly weaker and therefore shift the result away from what we find in our lattice analysis. As the authors point out, their resummation scheme fails to take into account O(λ 2 3 µ 2 i ) and O(λ 2 3 T 2 ) corrections to the thermal masses, which in turn are included in our 3d EFT and could explain the discrepancy. This ambiguity in resummation at 2-loop level is another reason why the 3d EFT approach is preferable over calculations in the full theory.

Validity of the effective theory
Let us reiterate that while our nonperturbative analysis in the 3d EFT is exact within statistical errors of Monte Carlo methods, the derivation of the EFT by dimensional reduction is still performed perturbatively. Therefore, estimating the perturbative accuracy of the dimensional reduction is a crucial part of our comparison with fully perturbative methods. In this work, we have applied the dimensional reduction of Ref. [80], which is performed at 1-loop level for the couplings, and at 2-loop level for the masses. In the SM, where all couplings are small compared to the scalar couplings in our BM points, this next-to-leading order (NLO) calculation (cf. dimensional reduction at leading order meaning tree-level matching for couplings and one loop for masses) is highly accurate and reproduces the equilibrium thermodynamics of the EWPT with errors of only 1% or less [53].
At NLO, the couplings do not obtain significant loop corrections, while the mass parameters in the EFT -the Debye screened masses -have the schematic form where µ 2 0 is the corresponding mass parameter in the full theory and the loop corrections Π 1-loop and Π 2-loop are of the order O(λT 2 ) and O(λ 2 T 2 ), respectively, with additional mass-dependent corrections coming from the high-T expansion. When running of the parameters is taken into account, dependence on the RG scale Λ can be shown to cancel exactly up to corrections formally of 3-loop order [53,80]. The phase transition occurs when at least one of the doublets in the EFT becomes very light, for which a cancellation between the tree-level mass µ 2 0 and the thermal loop corrections is necessary. A qualitative description may already be obtained at leading order with just the 1-loop thermal mass, but for quantitative results the 2-loop part in Eq. (6.15) can be significant, especially when some couplings are large. In BM1, the relative importance of the 2-loop corrections for mass parametersμ 2 11 ,μ 2 22 andμ 2 12 are 13, 17 and 2 percent, respectively, when the temperature is fixed to the critical value obtained from the simulations. In BM2, the respective numbers are 11, 19 and 2 percent.
Although the above numbers do not immediately signal bad convergence, it is possible for corrections at higher loop orders to be substantial. We can estimate the importance of higher-order effects of O(λ 3 ) and higher by varying the RG scale Λ in Eq. (6.15) and the other matching relations (see Ref. [80]). This leads to increased uncertainty in our results within the 3d EFT, and while it is difficult to quantitatively estimate this effect on the lattice due to the computational effort required, we may use the 3d V eff to address the issue. Hence, we have repeated the perturbative analysis of the previous section for the 3d V eff with different scales used in the dimensional reduction. We chose the same range for Λ which was used for the V eff in the full theory, 0.5πT to 1.5πT , and the results are shown in Table 5.  Table 5: Results in 3d perturbation theory with the renormalization scale of dimensional reduction varied from 0.5πT to 1.5πT , and the 3d scale varied from 0.5T to 2T as in Table 4. For corresponding scale uncertainties, we show the most pessimistic values in the error bars.
In BM1, the variation of Λ amounts to a slight increase of the uncertainty compared to the earlier case in Table 4, where Λ was fixed as in Eq. (5.1). In BM2, there is a clear shift in central values and the uncertainties are substantially larger than in the case of fixed Λ. We therefore estimate that higher order effects in the dimensional reduction procedure could lead to inaccuracies of a few percents in the thermodynamic quantities for the BM1 case, while for BM2, the uncertainty may be as large as few tens of percents, which can compromise quantitative predictions.
From the simple consideration above, it would clearly be preferable to include nextto-next-to leading order (NNLO) contributions to the dimensional reduction. For the parameters of the effective theory in Eq. (3.1), this means 2-loop contributions in the couplings and 3-loop contributions in the masses, and also higher-order terms originating from high-T expansions. While this calculation is highly non-trivial, important simplifications can be made by focusing only on the most dominant contributions from the large scalar couplings at order O(λ 3 ). Recent developments in calculating higher-order corrections to the dimensionally-reduced EFT of hot QCD (see [112,113] and the references therein) motivate applying similar techniques in the context of the EWPT. However, we will not consider such improvements in this study.
Finally, at NNLO it is no longer justified to neglect higher-dimensional operators from the 3d EFT. In our case, we expect the most dominant operators to be scalar operators of the type (φ † φ) 3 3d , which in principle are straightforward to include in the EFT by matching scalar 6-point functions. Unfortunately, these higher-dimension operators ruin the super-renormalizibility of the 3d EFT, turning it into "only" a renormalizable theory. Consequently, lattice analyses are complicated as the relations to continuum parameters are no longer exact. We may nevertheless estimate the effect of these operators by including them in the perturbative V eff . If the thermodynamic quantities obtained from this potential differ considerably from those in Table 4, we can then conclude that the higher-order operators' contributions are too significant to be ignored in the lattice simulations.
A systematic determination of the aforementioned NNLO contributions, as well as a numerical study of their effects in BM1 and BM2, will be published in a separate work. According to preliminary results of this work, one can estimate that the inclusion of higherorder operators weakens the transition by a few percentages in BM1 and about 20 percents in BM2, while the critical temperature is not affected significantly. This estimate has been obtained by performing a 1-loop dimensional reduction for all scalar operators of dimension six (ignoring operators containing derivatives) and including their effect in the 3d effective potential at a full 2-loop level.
Combined with the scale variation estimate above, we approximate that the accuracy of the dimensional reduction for thermodynamic quantities is within a few percentages in BM1. In BM2 with somewhat larger couplings, the accuracy is significantly worse, of the order 20%, and suggests that 3-loop corrections should be included if one seeks quantitative results. For even larger couplings -such as those used in Refs. [34,40] to produce very strong transitions in the heavy M A = M H ± regime -we believe that perturbation theory may fail to give even a qualitative picture of the EWPT. New techniques relying neither on the perturbative effective potential nor dimensional reduction are then required for reliable results.
We emphasize that the shortcomings related to NNLO corrections in the dimensional reduction are also present in perturbative analyses in the full theory, such as the 1-loop V eff , Eq. (6.7), in the form of missing higher-order corrections. However, due to the 2-loop mass corrections from dimensional reduction and the efficient handling of IR resummations, our EFT approach is still superior to the frequently-used 1-loop V eff approach.

Conclusions
In this paper, we perform a state-of-art study of the equilibrium properties of the electroweak phase transition in the Two Higgs Doublet model. The main analysis is based on a dimensionally-reduced effective theory [53,80], obtained by integrating out the heavy thermal field modes while simultaneously incorporating thermal resummations beyond the leading order. Using the effective theory, we perform nonperturbative lattice simulations in two benchmark points motivated by model phenomenology. As a comparison, we also apply conventional perturbative methods to compute the critical temperature and the strength of the phase transition at 1-loop level.
In simple beyond the Standard Model settings where only one phase transition occurs near the electroweak scale, it is necessary to have sufficiently strong interactions in the scalar sector in order to produce a strong transition. We demonstrate that this requirement of large couplings, substantially reduces the predictive power of traditional perturbative approaches, due to the significant corrections at two loops and beyond. A clear sign of such behaviour is the large uncertainty arising from residual dependence on the renormalization scale. Furthermore, in both of our benchmark points the ring-improved 1-loop effective potential fails to reproduce the nonperturbative critical temperature T c obtained from lattice simulations. We argue that this is due to higher-order resummations missing from the 1-loop potential.
On the contrary, resummations beyond the leading order are conveniently incorporated by performing dimensional reduction on the high-temperature theory. In section 6 we show that the effective potential evaluated within the high-T effective theory -for which a 2loop calculation is straightforward -displays better convergence than its counterpart in the full theory and agrees with the lattice results well within accuracy. This result suggests that the aforementioned shortcomings of the resummed perturbation theory are not due to nonperturbative effects related to the ultrasoft thermal modes, but rather a consequence of bad convergence caused by large scalar couplings.
It should be pointed out that for cosmological applications, the thermodynamic quantities in Table 4 should be obtained at the nucleation temperature T n instead of the critical temperature T c . The computation of T n nevertheless requires precise knowledge of the critical temperature, which, according to our results, is beyond the reach of 1-loop resummed perturbation theory, and the 2-loop correction calculated in Ref. [43] does not significantly improve the result for T c . For this reason, we advocate the use of the effective-theory approach even if the study is performed purely perturbatively, without numerical simulations. This strategy has already been recommended earlier in Refs. [43,111].
For both baryogenesis and gravitational-wave production, the phase transitions considered here are possibly too weak. A tempting resolution could be to move in the parameter space towards even larger couplings and consequently larger potential barriers -a lapse that would come with a high price, as perturbative convergence will then only worsen. In fact, the accuracy of many results in recent literature which invoke large scalar couplings could be in jeopardy due to the 1-loop potential simply being too inaccurate in describing the phase transition. We believe that more care needs to be taken in such scenarios, if one is willing to push perturbation theory to its limits. Unfortunately, for very large couplings even the effective-theory approach is likely to fail, and purely nonperturbative methods would then be necessary. Naturally, one would expect this conclusion to also hold for other models which rely on large couplings for a strong phase transition.
Finally, we highlight an alternative setup for strong phase transitions where more elaborate dynamics of multiple light fields are responsible for strengthening the transition. While such multi-step transitions typically call for a degree of fine tuning, the requirement of large couplings could be avoided in theses scenarios.  Table  4 of Ref. [34] A Renormalization in the MS scheme

Acknowledgments
As mentioned in section 2.2, the input parameters for our analysis are the pole masses of the scalars {M h , M H , M A , M H ± }, the scheme-dependent parameters {tan β, cos(β − α), µ 2 } and the gauge boson and fermion pole masses [72]. In order to relate these to parameters appearing in the Lagrangian, we apply a standard pole-mass renormalization at 1-loop level in the Minkowski space vacuum. As some of the couplings and masses in our theory are fairly large, loop effects can modify the relations between the two substantially. We emphasize that although this calculation is important for making accurate physical predictions, it is not directly related to the high-T behavior of the theory, which is the main focus of our paper. Starting with the CP-conserving 2HDM, we first require that the VEVs, v i (c.f. Eq. (2.2)), minimize the tree-level potential, and rotate the fields to a diagonal basis as in Eq. (2.3). The parameters {λ 1−5 , µ 2 11 , µ 2 22 } can be expressed in terms of the mass eigenvalues m i , the mixing angles and µ 2 ; the explicit relations can be found in e.g. Appendix B.2 of Ref. [80].
In the renormalized theory, dressed propagators are of the form where Π i (p 2 ) denotes a self energy evaluated at external momentum p and MS scale Λ. Because the minimization conditions (A.1) are imposed only at tree level, the VEVs generate one-particle-reducible tadpole contributions to the self energies that need to be accounted for, and are included in our study. The condition that the input masses M i correspond to poles of the loop-corrected propagators is then equal to having which is a system of four equations for the eigenvalues {m h , m H , m A , m H ± }. Similar polemass conditions are obtained for the gauge fields W ± , Z as well as for the top quark t (a detailed discussion can be found in Ref. [53]). Pole conditions for the light fermions are not needed in our case, as their Yukawa couplings are negligible in the phase-transition analysis.
Together with the direct input parameters tan β, cos(β − α) and µ 2 , these equations are sufficient to fix all but one parameter in the electroweak sector at some input scale Λ 0 . The final relation is obtained from charged particle scattering in the Thomson limit. This fixes the electromagnetic fine-structure constant α EM = 1 4π in the MS scheme viâ where α EM = 1/137.036 [72], and the photon self energy Π γγ and the (traverse) Zγ correlation function are computed at one loop. In the 2HDM, δα EM /α EM obtains a subdominant correction from H ± loops. Using the above relations and omitting terms beyond leading order in the correlation functions, one obtains loop-corrected expressions for the renormalized parameters: µ 2 11 (Λ) = µ 2 t β + will not be listed here. Explicit formulas can however be found in Ref. [114], and we have verified that our results match the expressions in their Appendix C. In practice, we have used the Feynman-t'Hooft gauge to simplify the calculation of the self energies. To 1-loop order, the UV counterterms required for renormalization are the same as those used in dimensional reduction. Apart from gauge-dependent terms, the counterterms can be found in Ref. [80]. Eqs. (A.6)-(A.15) form a non-linear system of equations for the renormalized parameters, as the correlations functions and δα EM /α EM themselves are functions of the MS couplings. The situation can, however, be simplified by substituting the renormalized parameters with the physical, scheme-independent parameters inside the loop corrections, i.e, making the replacements m i (Λ) → M i andα EM (Λ) → α EM , and the difference between the two prescriptions is formally of higher order in perturbation theory.
In the IDM limit (µ 2 = 0 and v 1 = 0), the calculation proceeds analogously, but instead of tan β and cos(β − α) we input the self-coupling λ 1 (Λ) and the combination λ 3 (Λ) + λ 4 (Λ) + λ 5 (Λ) /2. Detailed formulas in the IDM case (neglecting g contributions to the self energies) can be found in Ref. [43]. In Ref. [43], it is also discussed how higherorder corrections to Eqs. (A.6)-(A.15) can be partially resummed by solving the equations "self-consistently", without performing the linearization. In BM1, we adopt their approach and solve the parameters iteratively, dropping g terms inside the loop corrections. This ensures that our study of thermal effects in BM1 is directly comparable to the 2-loop results of Ref. [43].
In BM2 however, the iterative approach does not converge very well, whereby we take the simplified approach. We find corrections of order 10% to the couplings relative to their tree-level values. For the smallest coupling λ 2 , however, the correction is ∼ 25%, while the mass parameter µ 2 22 is modified by ∼ 50%. Such large corrections, and the bad convergence of the iterative approach, again indicate that our BM2 is already close to the border of applicability of perturbation theory.