Real effective potentials for phase transitions in models with extended scalar sectors

The effective potential obtained by loop expansion is usually not real in the range of field values explored by its minima during a phase transition. We apply the optimized perturbation theory in a fixed gauge to singlet scalar extensions of the Standard Model in order to calculate a one-loop effective potential that is real by construction. We test this computational scheme by comparing such a potential obtained in Landau gauge to that derived based on the Higgs pole mass. We carry out the latter construction by imposing physical renormalization conditions, which yields a potential without residual regularization scale dependence. We use our effective potential to study the parameter dependence of the critical temperatures in a two-step phase transition of the form $(0,0)\to (0,w')\to (v,w)$ that occurs for decreasing temperature in scalar extensions of the SM with two vacuum expectation values $v$ and $w$.


Introduction
We know from measurements of anisotropies in the cosmic microwave background [1,2], Big Bang nucleosynthesis [3], and direct searches [4] that the Universe contains baryons and trace amounts of anti-baryons as ordinary matter. From the value of the baryonto-photon ratio η 6 · 10 −10 [5] we know that in the early Universe there should have been approximately one extra baryon for a billion baryon and anti-baryon pairs. These pairs could annihilate leaving the extra baryons to survive and form the Universe we know today. This observation cannot be explained by the Standard Model of particle physics (SM), as it does not predict sufficiently large violation of the CP symmetry [6,7] required by baryogenesis [8]. Nevertheless, it is an intriguing question whether the known baryon asymmetry can be explained by a suitable extension of the SM.
There are two main directions to explain generation of baryon asymetry in particle physics. One is the mechanism of baryogenesis [9,10] in which sphalerons provide the necessary baryon number violating processes. Baryogenesis requires a first order phase transition as part of Sakharov's conditions [8], which is strong enough to prevent the washout of the created baryon number. Another one is the mechanism of leptogenesis [11,12], which seeks to create a lepton number asymmetry by involving heavy neutral leptons (singlet states under the SM gauge group), then turns this lepton asymmetry into a baryon asymmetry by standard processes [13]. For leptogenesis one needs to know the temperature range where the relevant phase transition or transitions happened and follow the thermal evolution of the masses within this range.
Beyond the SM theories may provide a solution to the baryon asymmetry as they generally incorporate new degrees of freedom. In particular, models with extended scalar sectors have interesting phenomenology regarding phase transitions in the early Universe [14][15][16][17]. A first order phase transition requires a potential barrier between the symmetric and the symmetry breaking ground states that can either be generated already at tree level by a suitable Lagrangian [14], or at loop level by the thermal contributions of bosonic degrees of freedom [18]. It is known that in the SM, the thermal effects do not produce a sufficiently strong first order electroweak phase transition (EWPT) as the physical Higgs boson is too heavy [19,20]. In extensions of the SM, additional bosons coupled to the Brout-Englert-Higgs (BEH) field can strengthen this transition [21].
The simplest scalar extension of the SM adds one gauge singlet scalar. For a Z 2 symmetric scalar a tree-level potential barrier is excluded between the electroweak minimum and an additional minimum with a non-vanishing singlet vacuum expectation value (VEV) [14], while thermal effects are generally too weak to strongly affect the thermally generated symmetry breaking. In these extensions a strongly first order phase transition is only possible if the singlet VEV vanishes in the T = 0 broken phase (see e.g. Ref. [15,22] for a real singlet and Ref. [16] for a complex one). Hence, baryogenesis is ruled out if the second VEV is also nonzero. Nevertheless, leptogenesis is still possible in such models if heavy neutral leptons (HNLs) exist. For example, in the superweak extension of the SM (SWSM) [23] the HNLs are the right-handed neutrinos, which receive their masses from the finite VEV of a singlet scalar through the BEH mechanism. To study such scenarios of leptogenesis, we have to find the critical temperatures of the phase transition of the model.
In quantum field theory the effective potential is the relevant tool for thermodynamic calculations and the study of phase transitions. Evaluated at homogeneous field configuration from the Legendre transform of the generating functional, it is convex, infrared (IR) finite, and real 1 . Some functional approaches used to compute the effective potential indeed satisfy these expectations [25]. However, for models with non-convex (double-well) classical potentials the effective potential obtained by loop expansion at a fixed order is not convex in general [26]. Additionally, beyond leading order the effective potential is complex for small field values, and even worse, in the presence of massless modes it be-comes infrared divergent starting at three-loops (and its derivatives already at lower loops [27,28]). In this case the elimination of IR divergences from the effective potential can be achieved by resummation of the Goldstone self-energy [28,29]. In the SM the calculation of the effective potential was carried out recently at two-loop [30] and even at three-loop order [30][31][32][33]. Furthermore, the infrared instabilities can also be treated with dimensional reduction, where the heavy modes in the theory are integrated out and the problematic low energy theory is studied nonperturbatively on the lattice [34,35]. For recent applications of the method in scalar extensions of the SM see Refs. [36,37].
The origin of the imaginary part of the effective potential arising in between its minima within the conventional loop expansion is a consequence of the non-convexity of the classical potential [38]. In Ref. [39] the authors interpreted the imaginary part of the one-loop result as the decay rate of a specific Gaussian state centered about field values where the classical potential was concave. Indeed, there are no stable homogeneous equilibrium states in this region [40]. However, the loop expansion of the effective potential is reliable only when the defining path integral is dominated by a single saddle point, which is the case for convex classical potentials. For models with a non-convex classical potential the usual expansion fails to correctly approximate the exact effective potential as competing saddle points are usually not taken into account [41,42]. Mathematically the issue can also be formulated as the non-interchangeability of the loop expansion and the analytic continuation of a parameter that controls the convexity, as illustrated in Sec. 13.5 of Ref. [42]. A modified loop expansion, which takes into account two saddle points, was considered in Ref. [38]. This procedure leads to an effective potential that is convex and real. Also, a generalized effective potential that is real and that can be computed perturbatively was constructed in Ref. [43] by coupling the external current to a polynomial of the field. Motivated by the above considerations, our aim in this paper is to study the temperature-driven phase transition of scalar extended models with Z 2 symmetry as a first step in a phenomenological investigation of leptogenesis. In particular, we plan to find the ranges of the critical temperatures in a two-step phase transition where the singlet acquires a non-zero VEV at high temperature, and remains finite even at temperatures below the EWPT (for vanishing singlet VEV, see e.g. Ref. [44]).
In order to calculate the one-loop effective potential, which is enough to estimate the range of the critical temperature available in a given model, we employ the optimized perturbation theory (OPT) of Ref. [45]. This computational scheme was applied mainly in scalar models where it was shown to preserve Ward identities (Goldstone theorem) and renormalizability even at higher orders in the loop expansion [46] as it preserves the symmetries of the Lagrangian. For an application of the method to describe the BEC-BCS crossover in a non-renormalizable fermionic model, see Ref. [47]. The OPT scheme yields an infrared finite and real effective potential in a sufficiently wide, relevant region of the scalar fields involved. The reality of the OPT potential is ensured by the convexity of the classical potential, which in turn ensures that the conventional loop expansion is reliable. Thus the application of the OPT solves the problem of the conventional loop expansion which does not give a real effective potential at vanishing field value, even at high temperature when it becomes the stable ground state of the system.
The paper is organized as follows. In Sec. 2 we deal with the SM one-loop effective potential at zero temperature. Using Landau gauge, we construct an IR finite effective potential which is used as a benchmark for the effective potential calculated with the OPT method. We also discuss the physical parametrization of these potentials. In Sec. 3 we discuss the tree-level potential and construct at zero temperature the one-loop OPT effective potential in a complex scalar extension of the SM. The parametrizations of the one-loop OPT potential of this and the SWSM model are discussed in Sec. 4. In Sec. 5 we provide the one-loop thermal corrections to the OPT potential and investigate the phase transition of the SWSM model. We present several auxiliary computations in the Appendices: the Higgs self-energy needed for the benchmark potential in App. A and steps of the construction of the benchmark potential in App. B, while the derivation of the thermal masses is discussed at length in App. C.

Standard model effective potential at one loop
The SM scalar sector includes an SU(2) L doublet, the BEH field, that in the spontaneously broken phase is parametrized in terms of real scalar fields as: .
Here φ i (x) and h(x) are the Goldstone and Higgs boson fields and v is the constant background field about which the Higgs field oscillates. The associated scalar potential depends on the mass parameter µ 2 and self-coupling λ. The classical potential is obtained upon neglecting all oscillating fields in the BEH doublet, and taking the potential to be a function of the homogeneous background field v: Here we explicitly indicated the dependence on the mass squared parameter for later convenience. In order for this potential to have a non-trivial minimum, we require that µ 2 < 0 and λ > 0. These tree-level parameters are determined by the VEV v 0 and the mass M h of the Higgs boson (c.f. Table 1): . The positivity of λ also implies that the resulting classical potential is bounded from below.
The full effective potential is a sum of the classical potential and loop corrections, which we write formally as In particular, the first order correction is given by the sum of various one-loop contributions where each individual mode i with multiplicity n i has the Coleman-Weinberg (CW) form Here m 2 i (v; µ 2 ) is the tree-level background-dependent mass squared, and the overall sign is s i = −1 for fermions and s i = 1 for bosons. The CW form (2.6) is obtained using dimensional regularization (DR) and modified minimal subtraction scheme (MS), with Q as the regularization scale. The UV divergence proportional to D MS = −ε −1 + γ E − ln(4π) is canceled by appropriate counterterms. The value of the constant c i appearing in (2.6) is 3/2 for fermions and scalars and 5/6 for massive gauge bosons 2 .
In the SM all massive particles obtain their masses through the BEH mechanism. In principle every particle with a mass contributes to the effective potential, but we can safely neglect all but the heaviest ones. Working in Landau gauge (ξ = 0), in which the ghosts are massless, the background-dependent tree-level masses of the relevant particles are where g L and g Y are the SU(2) L and U(1) Y couplings and y t is the Yukawa coupling of the top quark. The multiplicities of these particles are n G = 3n h = 3, n W = 2n Z = 6 and n t = 12. There are two issues regarding the effective potential that need to be addressed. First, the second derivative of the one-loop effective potential, which is the curvature mass of the Higgs boson, is IR divergent in the minimum v 0 due to the massless Goldstone modes. Secondly, as µ 2 < 0, the Goldstone mass squared becomes negative for v < −µ 2 /λ ≡ v 0 (c.f. Eq. (2.7c)). This makes the effective potential complex for field values v < v 0 . A complex potential would indicate instability, which we know is not the case, thus this feature is understood to be an unphysical artifact of the loop expansion of the effective potential. We discuss these issues in detail in the following two sections.

Effective potential with parametrization based on the Higgs pole mass
If one tries to construct the one-loop effective potential using renormalization conditions that preserve the tree-level values of the minimum and the curvature mass of the Higgs boson [21], one encounters a logarithmic IR divergence caused by the massless Goldstone modes, collectively denoted by G. This is because, although the one-loop potential (2.6) is finite for a vanishing mass, its second derivative is not: the bubble integral at vanishing external momentum B(p 2 = 0; m G ) is IR divergent for m G (v 0 ; µ 2 ) = 0. The standard way to avoid the appearance of this logarithmic divergence is to impose a renormalization condition on the pole mass M p ≡ M h (see e.g. [48][49][50]) instead of the curvature mass of the Higgs boson. At the minimum v 0 of the effective potential the pole mass is given in terms of the self-energy Π of the Higgs boson as [51,52] With this procedure, the bubble integral appearing in the curvature mass is replaced by B(p 2 = M 2 p ; m G ), which is finite for a vanishing mass because p 2 = M 2 p acts as an IR regulator.
In order to formulate the pole mass based parametrization of the one-loop effective potential, we add the following term to the tree-level Lagrangian of the Higgs field which contains through the decomposition δZ = δ Z + ∆Z, δµ 2 = δ µ 2 + ∆µ 2 , δλ = δ λ + ∆λ, (2.9b) infinite counterterms and finite corrections to the tree-level couplings, both nonvanishing and scheme-dependent beyond tree level. Then, the usual relations between the bare couplings λ 0 , µ 2 0 and the renormalized ones are µ 2 0 = µ 2 R + δ µ 2 , λ 0 = λ R + δ λ , with µ 2 R = µ 2 + ∆µ 2 and λ R = λ + ∆λ. Here µ 2 and λ are finite couplings that satisfy the tree-level relations given below (2.3) by construction and enter into the expression of the tree-level masses, while δ Z and ∆Z represent the infinite and finite part of the wave function renormalization factor.
Working in a strict perturbative expansion, we choose the finite corrections to the treelevel parameters, ∆Z, ∆µ 2 , and ∆λ, such that they remove all finite one-loop quantum corrections from the Higgs one-and two-point functions. This treatment ensures that the pole mass, the residue at the pole, and the minimum of the potential do not change compared to their respective tree-level values. The counterterms δ Z , δ µ 2 and δ λ , which remove UV divergences in the MS scheme, are determined based on the Higgs self-energy given in App. A. They cancel the divergences in the potential introduced in (2.5).
The renormalized one-loop effective potential obtained with the above procedure, pre- where Ω = − i =G n i m 4 i (v 0 ) 128π 2 is an uninteresting finite constant, while the functions l i (r) are with F (r) = −1 + 4/r and A(r) = arctan 1/F (r) . These functions encode what remained from ∆µ 2 and ∆λ after using them to bring the contributions of the modes into the form seen in the last two terms of (2.10). Without finite wave function renormalization (∆Z = 0), the functions were computed in Ref. [48] for regularization scale Q = M t . We provide a comparison of our formulae to those in Ref. [48] in App. B. The requirement of the unit residue of the Higgs propagator leads to an improved oneloop effective potential in Eq. (2.10) as compared to the one presented in Ref. [48], because the function V [1] (v) is explicitly independent of the regularization scale. We also mention that we obtained V [1] (v) through IR finite counterterms, as opposed to those appearing in Eq. (91) of [53] where the problem posed by the massless Goldstone mode was not addressed.
The Higgs curvature mass can be determined from the potential given in (2.10). It turns out that the values given in Table 1 result in a Higgs curvature mass that is about 5 % smaller than the pole mass. This gives an idea on the size of the possible deviation if one parametrizes the effective potential based on the curvature mass instead of the pole mass as we do in the OPT approach in the next subsection.

Optimized perturbation theory approach
In this subsection we apply the OPT procedure [45] to the SM scalar potential in order to prevent the one-loop effective potential from becoming complex at small values of the background field, which occurs due to the presence of µ 2 < 0 in the tree-level mass squared formulae of the scalars. We do not discuss the renormalization of the OPT scheme, which can be done as in Ref. [45]. Here, we simply drop the divergent piece D MS in the expression of the CW potential given in Eq. (2.6).
Following the approach outlined in Ref. [45], we change the classical potential by adding and subtracting a mass term in the Lagrangian where m 2 > 0. The second term is treated as an interaction term (or a finite part of the counterterm) [54] that contributes first at one-loop order. The new scalar potential takes the same expression as given in (2.2) with µ 2 changed to m 2 . Therefore, within the OPT scheme we employ the replacement µ 2 → m 2 in the tree-level mass squared formulae for the scalars in Eq. (2.7), making the Goldstone masses positive for any value of the background field v. When writing the effective potential at one-loop order in the OPT scheme, the interaction term (µ 2 − m 2 )|φ| 2 is added to the classical potential that contains now m 2 instead of µ 2 , resulting in V Note that the original µ 2 mass parameter is restored in the classical part, but in the genuinely one-loop contributions only the new m 2 mass parameter appears. Thus, if we find a physical parametrization with m 2 > 0, then the potential will be real for all v.
The one-loop effective potential has three unknown parameters: µ 2 , m 2 , and λ. In order to determine their values we need to impose three conditions on the form of the effective potential. In the SM, the two unknown parameters µ 2 and λ are determined by fixing the value of the VEV (minimum of the potential) and the value of the Higgs mass: Condition 2: Here we used the approximation that the Higgs pole and curvature masses are equal. In Sec. 2.1 we saw that the latter is smaller by about 5 %. For the VEV and Higgs boson mass we take the values given in Table 1. To fix the third parameter, we employ the principle of minimum sensitivity (PMS), namely that the potential at the minimum (physical point) should not depend on the value of m 2 : Choosing Q = M t for the regularization scale, with the value of M t given in Table 1, the solution of the three conditions is Since m 2 > 0, it is guaranteed that the OPT one-loop effective potential is real for all values of v, as shown in Fig. 1. As a comparison we have also plotted the classical scalar potential and the real part of the potential in Eq. (2.10) obtained in the previous subsection.

Effective potential in the singlet extension
Scalar extensions of the SM have been popular in the past decades [55][56][57][58] as they provide an interesting phenomenology from a particle physics as well as from a cosmological point Effective potential [GeV 4 Re[V [1] (v)] Figure 1. Comparison of SM effective potentials. The OPT potential (black, solid) is compared to the real part of the potential obtained from parametrization via the pole mass (red, dotted), and the classical potential (blue, dashed). The OPT potential matches the one obtained from the parametrization via the pole mass to good accuracy: less than 1 % relative difference below the minimum.
of view. The differences between particular models lie in the number of new scalars introduced, as well as in their transformation properties. Two very well studied models are the two Higgs doublet model (2HDM) [59] and the scalar singlet extension of the SM [17,57].
In this paper we are going to investigate the one-loop effective potential of the latter.

Classical potential
We consider the model where we add a single complex scalar field χ to the SM: Similarly to the BEH doublet, we allow the singlet scalar to have a non-vanishing vacuum expectation value w. The fluctuations around the constant background are described by the scalar fields s(x) and χ 2 (x) (the new Goldstone boson).
To allow for non-trivial phenomenology (i.e. having a not completely decoupled scalar), we allow mixing between the BEH doublet φ and the singlet χ, and we consider the potential from which the classical scalar potential is obtained via averaging over all fluctuating fields: As the classical potential is now a function of two background fields, the conditions for having a potential that is bounded from below, and has a non-trivial minimum are more involved. For λ φ,χ > 0 the requirement is [60] λ The boundedness condition in itself does not constrain positive values of λ . The background-dependent Goldstone masses are obtained as Due to the mixing of scalar fields, the states h(x) and s(x) are not mass-eigenstates. The curvature masses for the Higgs boson and the singlet scalar are obtained by diagonalizing the Hessian of the classical potential: The general formulae for the eigenvalues m 2 cl,± (v, w; µ 2 φ , µ 2 χ ) are complicated, but they can be simplified at the minimum (v 0 , w 0 ) of the potential. In this case the Goldstone masses in Eq. (3.5) vanish and we can eliminate the µ 2 φ,χ parameters in favour of the vacuum expectation values: Exploiting the relations in Eq. (3.7), the classical mass squared eigenvalues of Eq. (3.6) in the minimum of the potential are Requiring that both classical scalar masses are real leads to a constraint on the mixing λ : which extends the already existing boundedness condition of Eq. (3.4) to positive values of λ as well.
In our investigation of the singlet extension of the SM we set the smaller eigenvalue to be equal to the observed Higgs mass squared, and the larger one to be the yet unknown mass of the singlet: Experimental constraints exclude the possiblity of the opposite assignment [61,62]. The expressions for the masses in Eqs. (3.8) and (3.10) are satisfied by two separate parameter points {λ φ , λ χ } for any physical λ . This double covering of the same physical observables is removed by requiring that the Higgs mass is dominantly due to the vacuum expectation value of the φ field. This requirement is satisfied when Indeed, at λ = 0 the lighter mass is then given by m 2 cl,− = 2λ φ v 2 0 as in the SM. The limiting value λ max is determined for a given set {M s , w 0 } by the condition λ φ v 2 0 = λ χ w 2 0 , which gives There are no solutions to the parametrization if |λ | > λ max . This constraint on the mixing λ is stricter than that of Eq. (3.9), thus any parametrization in terms of non-negative squared scalar masses that results in real parameter values will yield a real potential bounded from below.

One-loop corrections
The tree-level treatment described in the previous subsection is general, that is, it holds for any theory in which the extension includes a singlet scalar. However, the corrections to it incorporate loops involving all fields in the model. We describe the general approach to the one-loop parametrization in this subsection, and present the formulae obtained in the minimal singlet extension and in a more complete model in the following section.
The singlet extended scalar potential in Eq. (3.2) has 5 unknown parameters: µ 2 φ , µ 2 χ , λ φ , λ χ , and λ . In order to fix them, we extend the procedure applied to the SM in Sec. 2.2. In addition to the known mass M h of the Higgs boson and the vacuum expectation value v 0 of the φ field, we pick a fixed value for the mass M s of the s scalar and another one for the vacuum expectation value w 0 of the χ field. These give 4 conditions for the 5 parameters of the potential. We investigate the parametrization by treating λ as a free parameter.
In general the one-loop effective potential is where the sum is over all fields present in the given model. The background-dependent masses in the one-loop correction depend in general on both vacuum expectation values. The singlet extended potential suffers from the same problem as that in the SM, it becomes complex in a certain region of {v, w}. To deal with this problem, we proceed as in Sec. 2.2 and introduce new mass parameters for both scalars (c.f. Eq. (2.12)) denoted with m 2 φ and m 2 χ : 14) The terms in square brackets are treated as interactions, first appearing in a one-loop level calculation of the effective potential. These terms restore the original mass parameters µ 2 φ/χ in the classical part of the one-loop level OPT effective potential, as in Eq. (2.13), while the genuinely one-loop part added to this contains the new mass parameters m 2 φ/χ : The conditions to fix the parameters in the effective potential Eq. (3.15) are similar to those used in the SM. First of all, v 0 and w 0 are imposed as minima of the bivariate potential: Secondly, we require that the curvature masses at the minimum of the potential are equal to the Higgs boson mass M h and the new scalar mass M s . Due to mixing between the physical states, one has to look at the eigenvalues of the Hessian 18) and relate those to the masses. Let these eigenvalues be given by m 2 We assumed that the Higgs boson is the lighter scalar, thus the conditions for the masses are: Thirdly, we have two stationarity conditions based on the PMS, Condition 6: ∂V We use these six conditions to fix the values of m 2 φ/χ , µ 2 φ/χ and λ φ/χ , while λ is kept as a free parameter. For a given model, the existence of a real solution to the parametrization conditions is not guaranteed due to the presence of logarithmic contributions involved [63]. While an estimate for the region of validity can be given in any model based on tree-level considerations, in general the loop contributions may be significant enough that the parameter regions where a solution exists are only available through trial and error.

Parametrization of SM extensions with a singlet scalar
In this section we are going to showcase the parametrization outlined in the previous section in two models. In both cases the potential V(φ, χ) in the Lagrangian is the same, given in Eq. (3.2). First, we are going to present the parametrization of a model that is the SM plus a complex singlet scalar only. Second, we take a look at the so-called super-weak extension of the standard model that includes a complex singlet scalar among other non-scalar new degrees of freedom.

SM with a singlet scalar
This model involves no new fields apart from χ(x), parametrized as in Eq. (3.1), thus the loop corrections to the effective potential will be the same as those in the SM, plus the contribution of the scalar s(x) and the Goldstone χ 2 (x). In this model the backgrounddependent tree-level masses of the gauge bosons and the top quark remains those given in Eqs. (2.7a)-(2.7b). All Goldstone and scalar masses are modified due to the mixing. The Goldstone masses were given in Eqs. (3.5a)-(3.5b), while the scalar masses are understood as the eigenvalues of the Hessian in Eq. (3.6).
We show an example parametrization in Fig. 2 with a scalar mass M s = 260 GeV, χ(x) vacuum expectation value w 0 = 10v 0 , and regularization scale Q = M t . The figure shows all 6 fitted parameters as functions of the quartic mixing λ . Solutions exist for a narrow range of λ , symmetric around λ = 0 in rough agreement with the conclusion we drew in Eq. (3.12) from a tree-level analysis. The characteristic feature of the figure is the elliptic shape of all parameter curves: there are two solutions at any value of λ that can be related to each other by a simple reflection. As discussed around Eq. (3.11), this doubling is expected and we can select the physical solution as that with the lighter scalar mass (i.e., the Higgs boson) mostly given by the vacuum expectation value of the φ field. As the scalar sectors completely decouple at λ = 0, the left-hand side plots in Fig. 2 should reproduce the SM result in Eq. (2.17) (shown with red stars). The physical selection in Eq. (3.11) is equivalent to choosing the half of the elliptic curve on which the SM point lies.
The reason for parametrizing the model using the OPT approach was to obtain an effective potential that is real for any value of v, w. In the SM this was shown to be possible in Sec. 2.2, however the SM is a theory with enough observables to uniquely fix all of its parameters. Any extension to the SM will have new degrees of freedom, thus yet unobserved values for certain physical quantities. However, there is no guarantee in the parametrization conditions that the potential will remain fully real for any solution. In Fig. 2 the parameter value m 2 φ becomes negative above 3 a certain value of λ indicated by the dotted blue line. At these parameter values, the potential is complex in a given region centered at the origin of {v, w}, however it is real at the minimum.
The parametrization depends on the values of M s and w 0 . For a particular choice, a generic solution is shown in Fig. 2. In fact, a real parametrization does not always exist. The characteristic features of the parametrization can be understood mostly using treelevel relations detailed in Sec. 3.1. In what follows we explore the effects of changing the (mostly) arbitrary input values for M s and w 0 . From Eq. (3.12) we determine that the non-trivial sector for these parameters is when M s and M h do not coincide, and w 0 is not much larger than the scales set by the masses.
Changing the vacuum expectation value w 0 while keeping the mass M s fixed will result in parameter ranges shown in Fig. 3. In all plots we see the changing of the physical λ domain, as described by Eq. (3.12). For µ 2 φ and m 2 φ (and less obviously for µ 2 χ and m 2 χ ) we see a change in the tilt of the elliptical curves with changing VEVs. From tree-level relations one can find the difference between the values of µ 2 φ/χ for the two extreme values for λ = ±λ max , As a base assumption we have M s > M h , thus for an increasing w 0 we expect the solution curves for µ 2 φ and µ 2 χ to rotate clockwise, which is confirmed by Fig. 3. The solutions for λ χ in Fig. 3 shift upwards with decreasing vacuum expectation values w 0 . The border between the physical and unphysical solutions at tree-level satisfies the simple relation λ φ v 2 0 = λ χ w 2 0 (see Eq. (3.11)). In particular, the minimum value of the physical λ χ is expressed as which describes the effect we see in the last plot of Fig. 3. Next, we look at the parametrization of the singlet scalar extended SM for a fixed VEV w 0 with varying the M s scalar mass. The resulting parameters are shown in Fig. 4 for w 0 = 10v 0 and singlet scalar masses in the domain M s ∈ [200, 300] GeV. We see an interesting new feature when M s is varied: m 2 χ can become negative for a symmetric region around λ = 0, as indicated by the dotted lines in Fig. 4. This is not necessarily an issue, as the minimum point of the potential is still real (i.e., the parametrization exists). The largest singlet scalar mass for which the potential remains real at λ = 0 (i.e., at the completely decoupled scalar sector limit) is approximated by setting the values of µ 2 χ and λ χ to their tree-level expressions and solving the PMS condition, to obtain m 2 χ . The upper bound for having a real solution is independent of the vacuum expectation value w 0 , and we find the condition for having a real potential at λ = 0 to be  Larger M s values are still compatible with the parametrization, but the resulting potentials will become complex for an increasing range in λ . The parametrization breaks down (i.e., we would require complex parameters to satisfy the conditions at the minimum of the potential) for M s /Q 1.74.
The choice of the regularization scale Q influences the existence of real solutions to the parametrization, but their actual value is unphysical. From Eq. (4.5) we see that larger values for the unphysical scale Q can accommodate heavier singlet scalars while the potential remains real. In Fig. 5 we compare two parametrizations done with different regularization scales, Q = M t and Q = √ 2M t . The value of the quartic coupling depends logarithmically on the regularization scale, which is a standard result, while the mass squared parameter µ 2 φ/χ is independent of it. The shifted mass parameters m 2 φ/χ introduced in the OPT scheme only appear in the one-loop part of the effective potential, thus their values depend heavily on the regularization scale. It is seen in Fig. 5 that a parametrization done at a larger Q results in larger positive mass squared parameters m 2 φ/χ , thus a real potential even for heavier scalars, in accordance with Eq. (4.5).

Super-weak model
In this section we apply the OPT formalism to a model which includes extra massive fermionic degrees of freedom. The model chosen here is the SWSM [23]. This model is an economic U(1) z extension of the SM focusing on experimentally reachable light new physics, and probing the parameter space for various BSM scenarios. The spectrum of the model extends that of the SM by a complex singlet scalar field χ, 3 generations of righthanded neutrinos N i (RHN), and a massive gauge boson of the new U(1) z symmetry Z . The U(1) z symmetry is spontaneously broken by the non-zero VEV of χ, generating mass for Z and the RHNs through the Higgs mechanism. In Refs. [64,65] it was shown that the model is capable of explaining MeV scale sterile neutrino dark matter if the coupling g z ∼ O(10 −5 ) is tiny, and the VEV of χ is comparable to the SM VEV, 2v 0 < w 0 < 100v 0 . We use this region of the parameter space in the following.
The effective potential depends mostly on the heavy particles (m heavy M W ) in the spectrum. For the parameter space outlined above Z and the DM candidate RHN N 1 are light, and their small contribution to the effective potential is neglected. The masses of the heavier neutrinos only depend on the VEV of the χ field (SM contribution is neglected due to the see-saw mechanism): Each RHN carries n N j = 2 degree of freedom (particle and antiparticle), and the Yukawa couplings y N j are free parameters of the model. The masses of the remaining SM degrees of freedom remain approximately unchanged: the gauge boson masses were given in Eq. (2.7a) (an O g 2 z /g 2 Z 0 correction to the Z boson mass is neglected), and the top quark mass is (2.7b). The scalar sector is the same as that detailed in Sec. 3.1, in particular see Eq. (3.5) and Eq. (3.6) for the Goldstone and scalar masses. The addition of RHNs to the spectrum only affects the singlet direction in the effective potential, as their masses only depend on w. Additional fermions decrease the one-loop value of the mass squared parameter, but increase that of the quartic coupling. Assuming that in the SWSM the two heavy RHNs have approximately the same mass, the one-loop corrections to the parameters at λ = 0 are approximately (neglecting all contributions apart from the RHNs) (4.7) As w 0 is constrained to be larger than the BEH VEV, the contribution of RHNs is largely negligible unless their masses lie well above the electroweak scale. An example parametrization is shown in Fig. 6 for various RHN masses for a given value of the scalar mass M s and VEV w 0 . The figure shows that the SM sector parameters (µ 2 φ , m 2 φ , and λ φ ) are indeed independent of the choice for the RHN masses. The RHNs shift the singlet scalar sector parameters as described by Eq. (4.7). The SWSM does not change the parametrization significantly as compared to the singlet scalar extension, and we refer to the previous subsection, Sec. 4.1, for the qualitative picture of the effects of varying the scalar mass M s or the VEV w 0 .
In this section, we have shown the OPT parametrization of the SWSM effective potential. The tiny gauge coupling renders the description shown here essentially equivalent to any extension which includes extra fermions and a complex singlet scalar.

Effective potential at one loop and finite temperature
The finite temperature corrections to the effective potential are given by [66] where m 2 i is the background dependent tree-level mass. The sum involves in principle all particles in the theory. Similarly to the zero temperature case, we can safely neglect the contribution of the light modes, and only consider the heavy degrees of freedom.
To cure the IR divergency produced by the static (n = 0) Matsubara mode, it is costumary to perform a daisy resummation on the bosonic modes, which then adds a thermal Debye mass to their zero temperature one [18]. Note that for the gauge bosons only the longitudinal masses are modified, since at high temperature the leading order of the transverse part of the self-energy is ∝ mT whereas for the longitudinal part it is ∝ T 2 [67]. Following standard notation, we denote the mass eigenvalues including the thermal contributions (i.e. the thermal masses) as m i . We give the thermal masses in Appendix C. In the Arnold-Espinosa approach [68] the thermal functions J ± (m 2 , T ) are given by (5.2) For example, in the SM the longitudinal modes would correspond to W L , Z L , and γ L , the scalars are the Higgs and the Goldstone bosons, while W T , Z T are the transverse modes. Even though the photon is massless at zero temperature, its longitudinal component (which does not even exist in the vacuum) gets a non-zero contribution at finite temperature. We introduced I ± as the fermionic (+) or bosonic (−) integral The numerical evaluation of this integral for any T and background field value is computationally expensive, so we use the semi-analytic method described in Sec. 2.2 of Ref. [69] to approximate its result. The idea of this approximation is to expand the integral for small and large positive values of a 2 [21], then connect the two series continuously at some particular intermediate point a 0 where the derivatives match. Using sufficiently many terms in both expansions, the analytic result matches that of the numerical one to great accuracy. With increasing temperature the minimum position of the potential changes, and at sufficiently large temperatures the system eventually settles in its symmetric phase, i.e. all VEVs become zero. The order of the phase transition from the broken phase to the symmetric one and the transition temperature T c can be estimated by analyzing the potential at finite temperature. As described in the introduction, a strongly first order phase transition from tree level effects is excluded in a singlet scalar extension of the SM where the broken phase has a non-zero singlet VEV. Then only thermal effects create a barrier between the degenerate symmetric-and broken-phase minima, and two competing effects have to be considered for the EWPT. On one hand, requiring large couplings of the singlet to the Higgs boson results in heavy (M s T c ) singlets that do not have any significant effect at EWPT due to decoupling. On the other hand, a lighter scalar (M s = O(T c )) is not decoupled, but consequently it has a small coupling to the Higgs and its contributions are negligible. A strongly first order EWPT is thus not expected in our investigations of the singlet extension, and we will not consider it. We emphasize however, that the OPT approach can be employed straightforwardly in models displaying a strong first order EWPT, such as singlet extensions of the SM with vanishing singlet VEV.
In order to be able to track the evolution of the temperature dependent minimum of the effective potential, we need the potential to be real at any value of the background fields. In the symmetric phase (at high temperature) the vacuum expectation values of the scalar fields vanish and the minimum of the potential is at the origin. However, in conventional calculations the mass squared parameter is negative and the potential is complex below the zero temperature minimum (see Eq. (2.10)), and minimization of the potential is not defined here unless the real part is taken. Without such a procedure the symmetric phase cannot be recovered.
In the OPT scheme the masses of the scalars and Goldstone bosons could be made real for all values of the background fields at T = 0. Consequently the effective potential is generally real and its minimzation is well defined. In the leading order of the high temperature expansion (HTE) the thermal corrections to particle masses are given in Appendix C. For fields where the mass eigenstates differ from those originally in the Lagrangian (gauge eigenstates), thermal masses need to be calculated by diagonalizing the full temperaturedependent mass matrix. However, the thermal contributions to the scalar mass squared may be negative if λ < 0, leading to an imaginary mass (and effective potential) at sufficiently high temperatures. Imaginary thermal masses are unphysical and the requirements for their reality can be viewed as constraints on the parameter space. These constraints are tenuous as there is nothing inherently wrong with having negative contributions in the thermal mass [70]. The appearance of the imaginary masses can be interpreted in two ways: either the given parameters are unphysical and should be excluded, or the model in question is not sensible at these temperatures (which sets the energy scale).
We present the thermal evolution of the minimum of the potential in the SWSM in Fig. 7. A benchmark point was chosen with M s = 200 GeV, M N = 150 GeV, and w 0 = 5v 0 where the absolute value of the mixing was fixed to approximately its maximal value, |λ | = 0.0394. Mixing indicates that the scalar fields that appear in the Lagrangian are not mass eigenstates, and we have to diagonalize the thermal mass matrix at every temperature to obtain the real masses as eigenvalues. The sign of the mixing also matters, as the evolution at late times (at low temperatures) is qualitatively different: positive mixing (left panel) implies that w(T ) approaches its T = 0 value from above, while for negative mixing (right panel) w 0 is approached from below as T → 0 [17]. With these parameters and within the OPT scheme the effective potential is real. The two panels differ in the sign of the mixing: in the left panel λ > 0 and w 0 is approached from above, while in the right panel λ < 0 and w 0 is approached from below as T → 0. Note that apart from the qualitative difference at low temperatures, there is also a quantitative difference between the critical temperatures.
While the overall shape of the thermal evolution of the minimum of the effective potential remains the same as that shown in Fig. 7, the value of the critical temperature for each phase transition depends on the model parameters. The critical temperatures in the SWSM depend on four parameters: the masses of the singlet scalar and the RHN M s and M N , the singlet scalar VEV w 0 , and the scalar mixing λ . The λ parameter takes values in a finite interval given by Eq. (3.12). Thus by fixing the masses we can study the critical temperatures as functions of only w 0 while varying λ ∈ [0, λ max ] (we focus on positive λ for simplicity to avoid possible imaginary thermal contributions). We show one such dependence in Fig. 8 for M s = Q = 300 GeV and M N = 150 GeV. The phase transitions happen within a temperature band depending on the value of λ . The EWPT critical temperature T EW c depends weakly on w 0 even for larger values of the scalar sector mixing λ . On the contrary, the SWSM phase transition T SWSM c changes roughly linearly with increasing w 0 . The observation of phenomenological relevance depicted in the right panel of Fig. 8 is that the SWSM phase transition had to happen at or above temperatures of O(1) TeV even for a light (but heavier than M h ) new scalar.

Conclusions
In this paper we applied the OPT scheme for the construction of a one-loop effective potential that is real for all values of the background field or fields. We first showed the viability of this approach on the SM effective potential by comparing it to a benchmark potential, then explored BSM scenarios. In particular, we studied in detail the parametrization of the SM with an additional singlet scalar, and looked at the finite temperature behavior of the effective potential in the SWSM. It is straightforward to extend the OPT method for other models with multiple scalar fields.
We constructed the benchmark SM potential in Landau gauge by revisiting the method presented in Ref. [48] to tackle the infrared divergence in the second derivative of the SM one-loop effective potential at zero temperature. In our approach we not only imposed the value of the pole mass of the one-loop Higgs propagator at the minimum of the potential, but also required that the pole should have a unit residue. As a result of the finite wave function renormalization, this procedure provides an IR finite one-loop effective potential without residual dependence on the regularization scale of the dimensional regularization scheme used, as opposed to the case when the unity of the residue is not imposed. The OPT potential reproduces this benchmark potential to high relative accuracy (within 1 %) even for temperatures where the EW phase transition takes place. We expect that after fixing the gauge, it is safe to use the OPT approach for the calculation of the effective potential in order to estimate the critical temperatures in phase transitions occurring also in scalar extensions of the SM.
We have explored phase transitions with the finite temperature effective potential of the SWSM via the OPT approach. By restricting our search to two-step phase transitions (0, 0) → (0, w ) → (v, w) with w = 0, we found that for a light (but heavier than M h ) new singlet scalar the EWPT depends only slightly on the BSM scalar sector parameters, whereas the singlet phase transition happens at temperatures at or above T = O(1) TeV. Moreover, the EWPT is first order (as in the SM) but too weak to provide baryogenesis in any region of the parameter space of the SWSM, as expected from general considerations of singlet extensions. However, the relatively high scale of the SWSM phase transition may provide ample time for leptogenesis. This will be investigated in a forthcoming work.
The p-independent divergence of Π bub is proportional to v 2 , hence it can be removed by 2v 2 δ λ (see Eq. (B.2)), with the coupling counterterm One can check using Eqs. (A.3) and (2.7c) that δ λ consistently cancels also the v-dependent divergence of Π tad . Its remaining divergence is removed by the mass counterterm Using the correspondence D MS → 1 + ln(M 2 /Λ 2 ) (see e.g. [71]) where M and Λ are the renormalization and regularization (here cutoff) scales, and Eqs. (12.50) and (12.53) of Ref. [72], namely one can verify that the above counterterms give the correct expression of the one-loop γ and β functions, appearing for instance in Eqs. (A34), (A35), and (A36) of Ref. [28].

B Construction of the SM potential based on the Higgs pole mass
In this Appendix we work out the procedure outlined in Sec. 2.1, giving the ingredients needed to reproduce the effective potential (2.10).
With the notation introduced in Eq. (2.9), the one-and two-point functions of the Higgs field have the following form at one-loop level where we used that the one-loop self-energy is a sum of tadpole and bubble diagrams. We also used that the contribution of the 1PI tadpole diagrams appearing in the twopoint function multiplied by v is the corresponding tadpole contribution to the one-point function. The infinite counterterms are used such that each individual square bracket is finite, which can be explicitly checked using the formulae of Appendix A. The finite part of any contribution to the self-energy is denoted in what follows by the superscript F, e.g., Such a decomposition is scheme dependent and we use the MS subtraction scheme, as in Appendix A.
The requirement that the minimum of the one-loop potential remains where it was in the tree-level one is equivalent to the requirement that v 0 , satisfying µ 2 + λv 2 0 = 0, is the nontrivial solution of (B.1), which implies that Two more relations between the finite parameters can be obtained from the Higgs propagator evaluated at the minimum v 0 of the potential. We use Taylor expansion withΠ F bub defined by the equality, and exploit the tree-level relation M 2 h = 2λv 2 0 to obtain Then, the requirement that the expressions in the two parentheses vanish, that is guarantees that p 2 = M 2 h is a pole of the Higgs propagator with unit residue due to the limit lim 3) and used the definition of δ Z in Eq. (A.5). Conditions similar to those in (B.7) were used in Ref. [73] for the renormalization of the Abelian-Higgs model. We note that while δ Z receives contribution only from the top quark and gauge bosons, all modes contribute to the finite wave function renormalization factor ∆Z.
The renormalized one-loop potential can be constructed using ∆µ 2 and ∆λ, given in Eqs. (B.4) and (B.7), and the infinite counterterms δ µ 2 and δ λ , but bringing the result in the form presented in Eqs. (2.10) proves difficult. A better strategy is to work with bare couplings. The classical potential written in terms of bare couplings reads where we separated a (bare) constant Ω 0 for convenience. We fix Ω 0 from the condition V [1] (v 0 ) = 0, which results in Our goal is to express the bare couplings µ 2 0 and λ 0 and bring the sum of (B.9) and (2.5) into the form (2.10). The equations determining µ 2 0 and λ 0 are The first equation is (B.1), written in terms of bare couplings. We separated in it the contribution of the Goldstone bosons. The contribution of all the other modes is written using the comment below (B.2) and is denoted with a bar on the quantity they contribute to, e.g.
CW (v, m 2 i (v, µ 2 )). The second equation is obtained from iG −1 h (M 2 h ; v 0 ) = 0, using the first equation of (B.7) and the relations in (and below) Eq. (2.9b) on the right hand side of Eq. (B.2).
In order to relate (B.10b) to the effective potential without introducing IR divergences, one should separate the contribution of the Goldstone modes. Using that the one-loop contribution of each non-Goldstone mode to the Higgs curvature mass squared is equal to its contribution to the one-loop self-energy evaluated at vanishing external momentum, one has the identity in which all terms are IR finite. Note that the contribution of the tadpoles cancels in the square bracket: . We denote this finite difference by ∆Π bub (M 2 h ; v 0 ). Using the relation M 2 h = 2λv 2 0 and Eq. (B.11), one can solve the equations in (B.10) for the bare couplings and find Having determined the bare couplings, one can use the explicit expression of the Higgs self-energy given in Appendix A to perform the final step of the calculation. Adding (B.9) to (2.5) and using the contribution of each mode to µ 2 0 and λ 0 given in (B.12), a straightforward calculation leads to (2.10), showing that all UV divergences and dependencies on the regularization scale cancel. In the absence of a finite wave-function renormalization a dependence on the regularization scale Q would remain from the contributions of the top quark and gauge bosons to ∆Π bub (M 2 h ; v 0 ). The derivativesV (1) andV (1) appearing in the expression of µ 2 0 and λ 0 are to be rewritten as mass derivatives for each mode involved. They combine withV (1) , giving the expression appearing in the last term of (2.10). This finite contribution of a non-Goldstone mode corresponds to the subtraction of the first three terms of the Taylor expansion of ln(p 2 + m 2 i (v)) about ln(p 2 + m 2 i (v 0 )) in the defining integral of V CW (m 2 i (v)). The functions listed in (2.11) are obtained from the term proportional to C in (B.12). Considering only the divergent part of the self-energy when taking the derivative in (B.12c) corresponds to ∆Z = 0, which is the case considered in Ref. [48]. A detailed comparison shows that for the top quark and the gauge bosons the result given there in (A1) is not in line with the Higgs self-energy computed e.g. in Ref. [51], with which we agree. Namely, the last term in the contribution of the top quark appearing there should be 4r, while concerning the real part of the gauge boson's contribution only the first term seems correct. 4 A final remark concerns the imaginary part of the self-energy that appears when contribution of individual modes is considered. As shown in [51], the imaginary parts produced by bubble integrals with vanishing mass add up to zero in the self-energy evaluated for v 0 and p 2 = M 2 h . In other words, the coefficient of B(M h ; 0) coming both from the Goldstone boson contribution and the last term in the contribution of gauge bosons in (A.2c) vanishes: Other bubble integrals do not produce an imaginary part when evaluated with physical masses at v 0 and p 2 = M 2 h . Therefore, only the real part of the bubble integral has to be considered when working with the contribution of individual modes to the self-energy.

C Thermal masses
In this appendix we give the thermal masses of various types of fields, calculated at leading order in HTE. First, general formulae are given for the different contributions, then these formulae are applied to the SM fields. Similar discussions can be found in Refs. [16,74]. For relations regarding the SU(N ) algebra we take the notation and formulae presented in Ref. [75].
We take the following general Lagrangian for a fermion ψ that transforms under some local gauge group: where T a are the generators of the group in the fundamental representation and g is the corresponding generic coupling between the fermions and the gauge fields. The SM is a gauge theory with SU(3) c ⊗SU(2) L ⊗U(1) Y local symmetries where the gauge covariant derivative of a field f (x) is where σ denotes the vector of Pauli matrices, λ is the vector of the 8 Gell-Mann matrix, i, j = 1, 2, 3, and a, b = 1, 2, . . . , 8. The heaviest particle in the SM is the top quark, for which the Yukawa term is In the terms involving the fields φ 3,4 the bottom quark field b appears (we neglect mixing present in the d-type quarks as we focus only on the top-bottom quark doublet). The mass 4 We obtain 3 4 RefW (r) = (12 − 4r + r 2 )F (r) arctan (1/F (r)) − 12 + 5r − 3r ln of the top quark m t = y t v/ √ 2 is generated by the Higgs mechanism when φ 1 acquires a non-zero vacuum expectation value v (c.f. first term in the Yukawa Lagrangian).
The scalar potential of the BEH field φ (parametrized as in Eq. (2.1)) in the SM involves a quartic term in the scalar fields

C.1 Scalar thermal mass
The thermal mass of scalar fields is obtained from the self-energy diagrams with fermion, scalar, and gauge boson loops [18]. For a real scalar field with a Yukawa coupling to a fermion such as that in Eq. (C.1), the thermal contribution is In the SM, the Yukawa Lagrangian for the top quark is given in Eq. (C.3). It is seen, that φ 1,2 and φ 3,4 have different vertices, nevertheless their thermal contribution is the same Given a scalar potential for a real field φ of the form V (φ) = λφ 4 /4, the two point function gets contributions from scalar self-interactions and yields the thermal contribution to the mass: In the SM the potential is given in Eq. (C.4). The self-energy for any φ i (i = 1, 2, 3, 4) is a sum of a tadpole involving the same field φ i and the 3 tadpoles involving φ i =j in the loop. The thermal contribution is The gauge boson contributions to the self-energy are due to the kinetic term of the bosons (c.f. Eq. (C.2)). The couplings between the scalar field and the gauge boson is proportional to the gauge boson mass M . Let M 2 (φ) = g 2 φ 2 /4, then In the SM, only the W ± and Z 0 gauge bosons are massive, thus In total, the SM Higgs gets the following thermal contribution to its mass: In the SWSM there exists a mixing term between the scalar fields L mix = λ |φ| 2 |χ| 2 which induces new diagrams to the self-energy. These contributions are calculated the same way as the ones in the SM where we put a φ i =j into the tadpole of the one-loop φ i propagator. The BEH field is a doublet, it involves 4 real scalar fields, whereas the complex singlet has 2 real scalars, thus the thermal contributions to the self-energy are Π (mix) φ = 2 · λ T 2 24 = λ T 2 12 and Π (mix) χ = 4 · λ T 2 24 = λ T 2 6 .
(C. 12) In order to calculate the thermal masses of the physical scalar states, we need to add the one-loop self-energy to the mass matrix and find the eigenvalues. We mention that the thermal contribution of the Z gauge boson is negligible due to M Z M Z , or put differently due to the feeble coupling O(g 2 z /g 2 Z 0 ) 1.

C.2 Gauge boson thermal mass
The gauge boson thermal masses are calculated in e.g. Refs. [18,67,76]. The gauge boson self-energy gets contributions from self-interactions (non-Abelian gauge fields only), fermion loops, and scalar loops. In the following we present the Debye masses (corresponding to taking the limit |q| → 0 at q 0 = 0) of the longitudinal modes of gauge fields [77], obtained at leading order in the HTE. The contribution due to self-interaction (derivative cubic and quartic vertices) and ghosts is Π (gb.) 00 (g) = where C A is the quadratic Casimir of the given SU(N ) group for the adjoint representation, and g represents either g L or g s for the SU(2) L and SU(3) c groups respectively. For a group with N ≥ 2 the Casimir is C A = 2T F N , with T F being the normalization of the generators in the fundamental representation, for which we use the conventional value T F = 1/2 for any N . For the fermionic contribution we have to sum over all fermion states that can enter the loop. We will relate the self-energy contributions of fermion fields to the well-known QED result Π QED 00 = e 2 T 2 3 . (C.14) Since QED is a non-chiral theory, the above self-energy is twice the corresponding contribution of an explicit chirality field. The interactions between gauge fields and fermions can be read from the covariant derivative in Eq. (C.2). The only fermions that transform non-trivially under the SU(3) c group are the quarks. Given N G = 3 fermion generations (quarks and leptons) there are N f = 2N G = 6 quark and lepton flavours. The quark contribution to the gluon self-energy is Π SU(3) 00,ab = N f tr(t a t b )Π QED 00 (e 2 → g 2 s ) = δ ab g 2 s T 2 , (C. 15) where t a = λ a /2 (a = 1, 2, ..., 8) are the half Gell-Mann matrices, and tr(t a t b ) = 1 2 δ ab .
To the U(1) Y gauge boson self-energy all fermions with non-zero hypercharge contribute. We denote the hypercharges of any field i by Y i . Left-and right-handed fields differ under the U(1) Y group, thus we sum over the specific modes (each separate chiral mode will have a contribution half that of Eq. (C.14)) of all fermions and find Π U(1) 00 Here Y Q (Y L ) is the hypercharge of the left-handed quark (lepton) doublets, while Y e , Y u , Y d denote the hypercharges of the respective right-handed fields.
The scalar contributions only appear for the SU(2) L and U(1) Y gauge fields, as the BEH field is a color singlet. The self-energy due to scalars is Π (φ) 00 (g) = In order to obtain the full thermal Debye masses for the longitudinal mode of the gauge bosons, we have to diagonalize their mass matrix with the temperature dependent selfenergies included: In the SWSM, the U(1) z gauge boson B weakly mixes with W 3 and B. This mixing modifies the SM thermal masses for Z L and γ L only at the level of O(g 2 z /g 2 Z 0 ) 1. Furthermore, the new mass eigenstate Z receives a thermal mass correction of order g 2 z , and thus its contribution to the effective potential is negligible.

C.3 Fermion thermal mass
For completeness we also present the fermion thermal masses. For chirally invariant gauge theories they were calculated at high temperature in Ref. [78].
For a general SU(N ) gauge group with a coupling g to the fermions (c.f. Eq. (C.1)), the contribution of gauge fields to the thermal mass is where C R is the eigenvalue of the quadratic Casimir operator in the representation R of the fermion field. In the SM, all fermion fields transform in the fundamental (F) representation. For a U(1) group C F = 1, and for SU(N ) it is given by C F = (N 2 − 1)/2N . Note that leftand right-handed fields couple differently to gauge fields, hence their respective thermal masses will differ. The scalar contribution from a single real scalar field φ with a Yukawa interaction given in Eq. (C.1) is The left-(L) and right-handed (R) top quark thermal masses in the SM are