Nucleation is More than Critical -- A Case Study of the Electroweak Phase Transition in the NMSSM

Electroweak baryogenesis is an attractive mechanism to generate the baryon asymmetry of the Universe via a strong first order electroweak phase transition. We compare the phase transition patterns suggested by the vacuum structure at the critical temperatures, at which local minima are degenerate, with those obtained from computing the probability for nucleation via tunneling through the barrier separating local minima. Heuristically, nucleation becomes difficult if the barrier between the local minima is too high, or if the distance (in field space) between the minima is too large. As an example of a model exhibiting such behavior, we study the Next-to-Minimal Supersymmetric Standard Model, whose scalar sector contains two SU(2) doublets and one gauge singlet. We find that the calculation of the nucleation probabilities prefers different regions of parameter space for a strong first order electroweak phase transition than the calculation based solely on the critical temperatures. Our results demonstrate that analyzing only the vacuum structure via the critical temperatures can provide a misleading picture of the phase transition patterns, and, in turn, of the parameter space suitable for electroweak baryogenesis.


Introduction
Cosmological observations suggest that our Universe has a large dark energy component, and that its matter component is dominated by an unknown form of dark matter [1]. Only 5 % of the energy budget of the Universe consist of the particles of the Standard Model (SM), mainly its baryons. Extensive tests at particle accelerators and other laboratory experiments have found no (unambiguous) deviations from the SM predictions. However, whereas the SM accurately describes the behavior of the particles making up the ordinary matter, it fails to give an explanation of how they came to be.
Under the assumption that particles and anti-particles are produced in equal number in the early Universe, the SM predicts that they would have long annihilated each other without leaving any remnant matter today. As first enunciated by Sakharov [2], producing a baryon asymmetry, i.e. more matter than anti-matter, requires baryon number violation, C and CP violation, and out-of-equilibrium processes to all occur at the same time. The SM does provide for C, CP, and baryon number violation through the electroweak interactions and sphalerons, respectively. The Electroweak Phase Transition (EWPT), however, is a smooth crossover in the SM and, thus, is not giving rise to sufficient deviations from thermal equilibrium [3]. In addition, the amount of C and CP violation in the SM is too small to generate the observed baryon asymmetry [4] even if the EWPT were to provide out of equilibrium conditions.
In order to generate the observed baryon asymmetry, sources of CP violation and out-of-equilibrium processes beyond those found in the SM must be realized in nature. One interesting possibility to achieve the latter is via a Strong First Order Electroweak Phase Transition (SFOEWPT), yielding promising conditions for electroweak baryogenesis. Accommodating a SFOEWPT demands modifications of the Higgs potential. Such modifications may be induced predominantly by thermal effects, as it happens e.g. in the Minimal Supersymmetric extension of the Standard Model (MSSM) [5][6][7][8][9][10][11], or by zero-temperature effects that have a lasting consequence after thermal effects are taken into account. The latter situation naturally occurs in models of new physics containing additional light scalar particles with sizable couplings to the Higgs.
To study the phase transition patterns of models with extended Higgs sectors, most previous works solely rely on analyses of the temperature-dependent vacuum structure via the computation of the critical temperature, T c , at which two (distinct) local minima of the effective potential become degenerate. While the critical temperature is indicative of the thermal history since it is the temperature at which the role of the global minimum passes from one vacuum phase to another, this calculation does not account for the probability of the associated phase transition actually taking place. First order phase transitions proceed via bubble nucleation, and the probability of the system transitioning from the false vacuum to the (new) true vacuum is computed via the bounce action, the Euclidean space-time integral over the effective Lagrangian, see, e.g., ref. [12] for a review.
Heuristically, bubble nucleation becomes difficult if the barrier separating two local minima becomes too high, or if the distance (in field space) separating the minima is too large. These conditions occur most readily if multiple scalar fields participate in the phase transition. For the EWPT, the possibility of a SM gauge singlet field participating in the phase transition is particularly interesting. While electroweak precision data tightly constrains the couplings and vacuum expectation values (vevs) of any fields charged under the electroweak symmetry, such constraints do not apply to gauge singlets. Since its couplings are free parameters, a gauge singlet field can radically alter the shape of the effective potential, enabling a SFOEWPT. On the other hand, a gauge singlet may induce large barriers separating local minima and acquire a large vev during the EWPT, increasing the distance between the local minima and reducing the nucleation probability. Therefore, a careful analysis of these effects is necessary in order to determine the region of parameter space leading to a successful SFOEWPT.
The Next-to-Minimal Supersymmetric extension of the Standard Model (NMSSM) [13,14] is a well-motivated example of physics beyond the SM that may solve the hierarchy problem of the electroweak scale [15][16][17][18][19] and provide a dark matter candidate [20][21][22][23][24][25][26][27][28][29][30]. Its scalar sector contains a (complex) gauge singlet and two SU (2) doublets, thus, it is wellsuited for a case study of the comparison of the phase transition patterns suggested by the critical temperature calculation and those obtained from calculating the nucleation probabilities. Moreover, the NMSSM provides a range of possibilities for C and CP violation beyond what is found in the SM. For example, CP violation can occur in the Higgs sector, or between the superpartners of the SM particles. Assuming the latter, CP violation in the Higgs sector is induced only via (small) quantum corrections, and one can study the EWPT in the CP-conserving limit of the scalar potential.
The EWPT in the NMSSM has been studied previously in the literature. To the best of our knowledge, Pietroni [31] was the first to consider electroweak baryogenesis in the NMSSM, noting that the dimensionful coupling of the singlet to the Higgs doublets, A λ , allows for shapes of the scalar potential suitable for a SFOEWPT at tree level. This is to be contrasted with the situation in the MSSM, where a barrier between the trivial and the physical minimum necessary for a SFOEWPT arises only from thermal effects. Subsequent work on the EWPT in the NMSSM includes refs. [32][33][34][35][36][37][38][39][40][41][42][43][44][45], and work on closely related models can be found in refs. [46,47]. Many of these papers focused on numerical scans of the NMSSM parameter space, aiming at identifying regions of parameter space suitable for realizing a SFOEWPT. Analytic studies have been carried out in refs. [31, 34-36, 38, 39]. A common idea in these works was to use parameters shaping the potential in the singlet-only direction to characterize the EWPT.
In the NMSSM, in general, there are ten degrees of freedom in the Higgs sector. In practice it suffices to consider the three-dimensional subspace spanned by the CP-even neutral scalar degrees of freedom. Nevertheless, computing the bounce action in this threedimensional field space is still numerically expensive, and, until now, results for the phase transition based on the nucleation calculation have only been presented for a few benchmark points in parameter space, see, e.g., refs. [33,36,40,43,45]. These studies mainly reported small-to-moderate supercooling, i.e. nucleation temperatures not much smaller than the corresponding critical temperatures for their benchmark points. More importantly, the thermal histories indicated by the critical temperatures agree with the ones obtained by the nucleation calculation. The notable exception is the recent work of Athron et al. [45], where results for the nucleation temperatures of four benchmark points were presented: For two of those four points, the authors reported small-to-moderate supercooling, while for the two remaining points the authors found that the nucleation condition could not be satisfied and, hence, the transition pattern indicated by the calculation of the critical temperatures was not a good indicator of the thermal history.
In this work, we present results for the EWPT in the NMSSM based on the nucleation calculation for a broad scan of the parameter space. We use CosmoTransitions [48] for the calculation of the bounce action, and support our results with analytic studies. We focus on the region of parameter space where alignment-without-decoupling is realized in the Higgs sector, and on small-to-moderate values of tan β, the ratio of the vevs of the scalar SU (2) doublets. This is motivated by the phenomenology of the 125 GeV Higgs boson observed at the Large Hadron Collider (LHC). In the NMSSM, a mass of 125 GeV of the SM-like Higgs boson can be achieved in the low-to-moderate tan β 5 regime without the need for large radiative corrections. The couplings of this state to SM particles are SM-like if it is (approximately) aligned with the interaction eigenstate that couples like the SM Higgs boson to other SM particles. In the NMSSM, there are two ways to achieve such alignment: i) the decoupling limit, that requires the non-SM-like interaction eigenstates to have masses much larger than the SM-like interaction state, and ii) the alignment-without-decoupling limit, where the parameters of the theory conspire to suppress the mixing of the SM-like interaction state with the non-SM-like interaction states [49]. The latter is of particular interest for realizing a SFOEWPT in the NMSSM: in the alignment-without-decoupling limit the non-SM-like states can have masses comparable to that of the SM-like Higgs boson, and hence, they can easily alter the shape of the scalar potential in ways relevant for the EWPT.
The null-results from searches for superpartners at the LHC suggests that the squarks and gluinos are heavy and decoupled from the EWPT. We use an effective field theory approach, integrating out all superpartners except for the neutralinos and charginos. This leaves the full SM particle content, an augmented scalar sector consisting of two SU (2) doublets and a complex singlet, and the electroweakinos (composed of the superpartners of the photon, the Z-and W -bosons, the two Higgs doublets, and the scalar singlet) as dynamical degrees of freedom; similar approaches have been taken in refs. [40,45,47]. In order to maintain the location of the physical minimum in field space, the mass of the SMlike Higgs boson, and the alignment of the singlet-like and SM-like interaction eigenstates after including the radiative corrections to the effective potential from these remaining dynamical degrees of freedom, we add a set of (finite) counterterms, see refs. [41][42][43] for similar schemes.
The outline of our work is as follows: We begin by discussing the scalar sector of the NMSSM in section 2. In section 2.1 we discuss the radiative corrections to the scalar sector of the NMSSM, and in section 2.2, the thermal corrections. After analyzing the zerotemperature vacuum structure of the NMSSM in section 2.3, we discuss the phase transition behavior of the NMSSM in section 2.4, in particular, we identify the relevant characteristics of the transition patterns for a SFOEWPT, and develop some analytical intuition for the regions of parameter space where phase transitions can successfully nucleate. In section 3, we present our numerical results. In section 3.1 we study the region of parameters in which the proper physical minimum is obtained. We compare the results for the phase transitions obtained from the nucleation calculation with the transition patterns suggested by the critical temperature analysis in section 3.2. In section 3.3 we comment on the collider and dark matter phenomenology in the region of parameter space where we find SFOEWPTs. We summarize and present our conclusions in section 4. Explicit formulae for the field-dependent masses, the finite temperature corrections to the masses, and the equations we use to fix the counterterms are listed in appendices A, B, and C, respectively.
Let us here already highlight our main result: We find that the phase transition patterns of given parameter points vary substantially between the critical temperature analysis and the nucleation calculation. Thus, calculating only critical temperatures is not enough to identify the regions of parameter space favorable for electroweak baryogenesis.

The Next-to-Minimal Supersymmetric Standard Model
The Next-to-Minimal Supersymmetric Standard Model augments the particle content of the MSSM by a SM gauge-singlet chiral superfield S, see refs. [13,14] for reviews. The best-studied version of the NMSSM is the Z 3 -NMSSM. In this model, an additional discrete symmetry is imposed, under which all left-handed chiral superfields transform as Φ → e 2πi/3 Φ and all gauge superfields transform trivially. An interesting consequence of the Z 3 symmetry is that it renders the superpotential of the NMSSM scale invariant; in particular the Higgsino mass parameter µ arises from the vacuum expectation value (vev) of the scalar component of the singlet superfield, S. Thus, the NMSSM alleviates the MSSM's µ-problem.
Of greater phenomenological interest is that the NMSSM can accommodate a 125 GeV SM-like Higgs boson without the need for large radiative corrections to its mass. Furthermore, the presence of the scalar gauge singlet makes a SFOEWPT easily achievable in the NMSSM [31-45, 47, 50]. This should be contrasted with the situation in the MSSM, where, in the presence of a 125 GeV SM-like Higgs, the scalar potential is constrained such that a SFOEWPT is only possible if the stops are very light [5][6][7][8][9][10][11]. Such stops have been virtually ruled out by the LHC, not only via direct searches but also by the fact that such light stops would lead to a variation of the Higgs production cross section and decay branching ratios that are in conflict with current Higgs precision measurement data [51][52][53][54][55][56][57]. This places severe pressure on the possibility of electroweak baryogenesis in the MSSM. In the NMSSM, the presence of the singlet S, the bosonic component of S, allows for radically different shapes of the scalar potential, which make a SFOEWPT possible in the NMSSM without the need for light stops.
The superpotential of the Z 3 -NMSSM is given by are the usual SU (2)-doublet Higgs superfields, we use a dot-notation for SU (2) products and W Yuk indicates the Yukawa terms which are identical to those in the MSSM [58]. Including F -, D-and soft SUSY-breaking terms, the scalar potential reads where m 2 i and A i are soft SUSY-breaking parameters of dimension mass-squared and mass, respectively, and g 1 and g 2 are the U (1) Y and SU (2) L gauge couplings.
The Higgs fields have large couplings amongst themselves, to the electroweak gauge bosons, and to third generation (s)fermions. These couplings lead to sizable radiative corrections to V 0 , to which we return in section 2.1. However, many of the properties of the scalar potential can already be seen from the tree level potential, eq. (2.3).
In order to be compatible with phenomenology, the NMSSM must preserve charge. While in the MSSM the scalar potential is sufficiently constrained to make charge-breaking minima very rare (see, e.g., ref. [59]), the additional freedom of the NMSSM's scalar potential makes such minima a much larger problem. However, ref. [60] demonstrated numerically that, while charge-breaking minima may be present in the NMSSM, they are virtually always accompanied by additional charge-conserving minima, and the tunneling rate from the metastable physical minimum to these charge-conserving minima is larger than to the charge-breaking minima. Hence, we can neglect such charge-breaking minima; in the following we will assume that for all phenomenologically relevant vacua the vevs can be rotated to have the form Without loss of generality, one can furthermore take all vevs to be real-valued: While the Z 3 -NMSSM does allow for stationary points in the scalar potential which spontaneously break CP, at tree level such points are either saddle points or local maxima [61]. In summary, it suffices to allow the neutral real components of H d , H u , and S to take non-trivial vevs 1 when studying the vacuum structure of the NMSSM. This reduction from a ten-dimensional to a three-dimensional field space makes the task considerably more tractable.
In order to ensure that the scalar potential has a stationary point at the physical minimum, we use the minimization conditions replacing the squared mass parameters m 2 H d , m 2 Hu , and m 2 S with the vevs v d , v u , and v S in eq. (2.3). In practice, it is convenient to re-parameterize the vevs, The observed mass of the electroweak gauge bosons is reproduced by fixing v = 174 GeV, removing one of the NMSSM's free parameters. In order to account for the constraints on the NMSSM imposed by the SM-like couplings of the observed 125 GeV Higgs boson, it is useful to write the Higgs fields in the extended Higgs basis [49,[62][63][64][65][66][67][68] H SM , H NSM , and H S are the three neutral CP-even interaction states of the Higgs basis, A NSM and A S are the CP-odd states, and H ± is the charged Higgs. The neutral and charged Goldstone modes are denoted by G 0 and G ± , respectively, and we used a shorthand notation In this basis, the couplings to pairs of SM particles take a particularly simple form. Focusing on the CP-even states, the couplings to pairs of down-type and up-type fermions and pairs of vector bosons (VV) are (2.11) where m 2 Z = v 2 g 2 1 + g 2 2 /2 is the Z-boson mass. While m h 125 receives sizable radiative corrections via the stops, see section 2.1, it is interesting to note that the term proportional to λ 2 v 2 allows one to obtain m h 125 = 125 GeV already at tree level for small values of tan β 3 if λ takes values 0.7 λ 1. Thus, there is no need for large radiative corrections to the mass of the SM-like Higgs, i.e. no need for heavy stops, in the NMSSM. Including moderate corrections from the stops, the required value for the mass of the SM-like Higgs boson is obtained for 0.6 λ 0.8 in the small-to-moderate tan β 5 regime.
In Here, the M 2 S,ij again are the entries of the squared mass matrix for the CP-even states in the basis H SM , H NSM , H S . There are two possibilities to achieve such (approximate) alignment of h 125 with H SM : either, the entries of the squared mass matrix corresponding to such mixing are small, or, the right hand sides of eqs. (2.15) and (2.16) become large. The latter option is the so-called decoupling limit. Realizing alignment in this way implies {m H , m h S } m h 125 . As we will see below, a relatively light singlet-like state gives the scalar potential a favorable shape for SFOEWPT. Thus, the former option, the so-called alignment without decoupling limit, is more interesting for electroweak baryogenesis.
At tree-level, alignment between the two states originating from the Higgs doublets, eq. (2.15), is achieved for It is convenient to instead rewrite this condition as because this form is robust against radiative corrections [49]. Identifying M 2 S,11 = m 2 h 125 , one obtains the alignment condition For small to moderate values of tan β, this condition yields 0.6 λ 0.7. It is interesting to note that, for moderate values of tan β 5, this range of λ coincides with the range for which one obtains m h 125 = 125 GeV without the need for large radiative corrections. Suppressing the mixing of H SM with H S , eq. (2.16), yields a second alignment condition from demanding M 2 S,13 → 0, namely where we introduced the parameter A is the (squared) mass parameter of A NSM and controls the mass scale of the mostly doublet-like CP-even and CP-odd mass eigenstates as well as the mass scale of the charged Higgs boson. The alignment condition eq. (2.20) gives rise to a mass spectrum where, provided κ < λ, the doublet-like mass eigenstates have approximate masses m H , m A , m H ± ∼ 2µ/ sin 2β [69,70].
In the remainder of this work, we will consider the NMSSM in the alignment limit, choosing parameters to satisfy eqs. (2.19) and (2.20). While current data [71,72] allow for some deviation from perfect alignment, the phenomenological impact of such departures on the EWPT in the NMSSM is small. Note also that in refs. [69,70] it was demonstrated that, in random parameter scans where the alignment conditions are not a priori enforced, requiring compatibility with the phenomenology of the observed 125 GeV Higgs boson selects the region of parameter space where eqs. (2.19) and (2.20) are (approximately) satisfied.
The NMSSM parameter space is constrained by a number of additional arguments. Let us briefly discuss two of them here, while we derive constraints arising from the stability of the electroweak vacuum in section 2.3. It is well known, that large values of the dimensionless parameters λ and κ lead to Landau poles. Avoiding the appearance of Landau poles below the GUT scale [Q GUT ∼ O(10 16 ) GeV] entails constraining the values of the NMSSM's couplings, at the electroweak scale, to [14] λ 2 + κ 2 0.7 . (2.22) As discussed above, both the SM-like nature of the observed Higgs boson and its mass value lead to a preference of sizable values of 0.6 λ 0.7 in the NMSSM. Hence, avoiding Landau poles below Q GUT limits the value of |κ| 0.3 in the alignment limit. Note that the NMSSM with larger couplings (and Landau poles between the TeV and the GUT scale) is known as λ-SUSY, see, for example, refs. [73][74][75].
The parameter space is also constrained by avoiding tachyonic masses. The most relevant constraint arises from the singlet-like CP-odd mass eigenstate a S . Taking into account first-order mixing effects, its mass is approximately [49] (2.23) Recalling that alignment requires M 2 A 4µ 2 / sin 2 (2β), we can deduce the condition the NMSSM parameters must satisfy to keep a S from becoming tachyonic: For small-to-moderate values of tan β and in the alignment limit, where 0.6 λ 0.7, the right-hand side of eq. (2.24) is approximately κ × O(1). Hence, equation (2.24) implies µA κ v 2 for κ > 0, while for κ < 0 the condition becomes µA κ v 2 ; in particular, disfavoring sgn(µA κ ) = −1 for κ < 0.

Radiative Corrections
The scalar potential receives sizable radiative corrections from the large couplings between the Higgs bosons themselves as well as from their large couplings to the electroweak gauge bosons and the (s)fermions, in particular the (s)tops, see, for example, refs. [14,[76][77][78]. Since the precise interplay between the higher-order corrections to the Higgs mass and the mass values of the SM particles and their superpartners does not play a relevant role in our study of the EWPT, we shall take only the dominant one loop corrections into account in this work. The null-results from SUSY searches at the LHC suggest that all squarks as well as the gluinos have masses 1 TeV. LHC constraints on new states neutral under QCD are less stringent. Furthermore, to yield a scalar potential sufficiently different from that of the SM to accommodate a SFOEWPT, the Higgs bosons' masses should not be much larger than the electroweak scale. These considerations motivate studying a scenario in which all sfermions 3 and the gluinos are heavy and can be integrated out, yielding an effective theory where the remaining dynamical degrees of freedom are the SM particles, the new Higgs bosons {H, h S , A, a S , H ± }, the five neutralinos χ 0 i , and the two charginos χ ± i ; see refs. [40,45,47] for similar approaches. The parameters of this effective theory are obtained by matching onto the full theory (containing all the NMSSM's degrees of freedom) at an intermediate scale. The leading operator one obtains from this procedure is where h t is the top Yukawa coupling determined from the (running) top quark mass m t = h t v sin β, M S is the geometric mean of the stop masses, and A t is the soft trilinear stop-Higgs coupling. We note that for small to moderate values of tan β, the top quark superfield has a sizable coupling only to H u in the superpotential. After the singlet acquires a nontrivial vev, an effective µ-term is generated and additional effective quartic couplings, which involve not only H u but also H d , arise via stop loops. However, these contributions are suppressed by powers of µ/M S . We shall work in a region of parameter space where |µ| M S and, hence, the dominant contribution induced by integrating out the stop sector is given by eq. (2.25). At higher loop orders, the exact relation between ∆λ 2 and the parameters in the stop sector is modified, but, for small values of |µ|, the stop radiative corrections can still be effectively parametrized by ∆λ 2 (see, for example, refs. [80,82,83]).
The scalar potential of this effective theory is then given by This new contribution gives sizable corrections to the Higgs mass matrix. In particular, the mass of the SM-like Higgs state is given by where F i = 0 for bosons and F i = 1 for fermions. The constant C i takes values C i = 3/2 for scalars, longitudinally polarized vector bosons, and fermions, while for transversal vector bosons C i = 1/2. We denote the field-dependent masses computed from V eff 0 by and work in the Landau gauge; explicit expressions for the m 2 i are collected in appendix A. The bosonic fields entering eq. (2.29) are B = h i , a i , H ± , G 0 , G ± , Z, W ± with n B = {1, 1, 2, 1, 2, 3, 6} degrees of freedom, respectively. Here, h i and a i denote the three neutral CP-even and two CP-odd Higgs bosons, H ± the charged Higgs, G 0 and G ± the neutral and charged Goldstone modes, and Z and W ± the electroweak gauge bosons. The fermionic fields entering the Coleman-Weinberg potential are 4 i and χ ± i denote the five neutralinos and two charginos, respectively, and t is the top quark. We have chosen m t as the renormalization scale, implying that the parameters are defined at such scale. In order to guarantee the one-loop renormalization scale independence and preserve the supersymmetric relations, the parameters at the scale m t must be related with those at higher energies, up to the supersymmetry breaking scale, by including all particles in the effective theory in the running to higher energies.
Note that since the Goldstone modes' masses vanish at the physical minimum, their contributions to the Coleman-Weinberg potential lead to divergent contributions to physical masses and coupling coefficients computed from derivatives of the loop-corrected effective potential. This divergence is an artefact of the perturbative calculation [85,86] and can be dealt with by shifting the masses of the Goldstone modes by an infrared regulator, m 2 G → m 2 G + µ 2 IR . In our numerical calculations, we use a value of µ 2 IR = 1 GeV 2 ; note, however, that in numerical calculations numerical errors on m 2 G typically suffice to "regulate" the logarithmically divergent contribution from m 2 G → 0, even before including an explicit infrared regulator.
Including the Coleman-Weinberg contributions, the (effective) scalar potential at zero temperature is given by The Coleman-Weinberg corrections alter the location of the minima as well as the physical masses. We include a set of counterterms to keep the location of the physical minimum at H SM , H NSM , H S = √ 2 {v, 0, µ/λ}, ensure M 2 S,13 → 0, preserving alignment, and maintain m h 125 = 125 GeV. Note that these counterterms correspond to a redefinition of the soft SUSY-breaking terms 5 , see refs. [41][42][43] for similar approaches. We list equations for the fixing of the counterterms in appendix B.
The input parameters for our model are thus All other parameters are fixed by the various conditions we impose on the model, namely, λ and M 2 A are determined by alignment, ∆λ 2 by setting m h 125 = 125 GeV, and the counterterms are fixed by the conditions discussed in the previous paragraph.

Thermal Corrections
So far, we have discussed the scalar potential at zero temperature. At finite temperatures, thermal corrections to the potential have to be taken into account. The one-loop finite temperature potential is given by where analogously to our definition of the Coleman-Weinberg potential, eq. (2.29), the sum runs over bosonic and fermionic degrees of freedom, n i counts the degrees of freedom of species i, and F i = 0 (F i = 1) for bosons (fermions). We denote thermal (field-dependent) masses with a tilde, m 2 i . Compared to the field-dependent masses, which we denote with a hat, m 2 i , the thermal masses include the so-called Daisy corrections re-summing hard thermal loops, The Daisy coefficients c i are only non-zero for bosonic fields. Furthermore, only the longitudinal polarization states of vector bosons receive non-zero Daisy corrections, gauge symmetry protects the transversal degrees of freedom. We list the Daisy coefficients for the relevant fields in appendix C. The thermal functions are defined as Following ref. [87], we improve the calculation of the thermal corrections by replacing the field-dependent masses with the thermal masses in the Coleman-Weinberg potential, Including the Coleman-Weinberg and the thermal corrections, the temperature-dependent effective potential at one-loop order is given by

Zero-Temperature Vacuum Structure
While the NMSSM's scalar potential is subject to radiative as well as thermal corrections as discussed in sections 2.1 and 2.2, one can already learn much about the possibility of a SFOEWPT from considering the effective potential, V eff 0 , obtained after integrating out all sfermions and the gluinos and prior to including the Coleman-Weinberg and thermal corrections. In this section, we derive the most interesting regions of NMSSM parameter space for realizing a SFOEWPT from V eff 0 . As we shall show later on, these regions of parameter space are only mildly affected by radiative corrections. Recall that in order to study the vacuum structure of the NMSSM, it suffices to consider the three-dimensional field space spanned by the neutral CP-even fields H SM , H NSM , H S , where V eff 0 is the potential given in eq. (2.27). As discussed above, the singlet plays a special role for realizing a SFOEWPT. Its coupling to the Higgs doublets, λ, and its self-coupling, κ, are free parameters, while the quartic couplings between the Higgs doublets are governed by the gauge couplings (and ∆λ 2 ). Furthermore, as a consequence of This residual Z 2 symmetry ensures that any extrema in the singlet-only direction, i.e. where H SM = H NSM = 0, are also extrema (or saddle points) of V eff,3 0 . In the alignment limit (or, more specifically, as long as the second alignment condition, eq. (2.16), is satisfied) the scalar potential in the singlet-only direction is given by This potential has extrema at The first of these field values corresponds to the trivial minimum of the scalar potential H SM = H NSM = H S = 0, and the second value coincides with the vev of H S at the physical minimum v S = µ/λ. The third field value marks a new special location in H S space, which, in the following, we refer to as Recall that since we used the minimization conditions, eq. (2.5), to replace the m 2 i parameters in the scalar potential with v, tan β, and µ, the physical minimum H SM , H NSM , H S = √ 2 {v, 0, µ/λ} is also guaranteed to be a stationary point of the scalar potential. Hence, in the alignment limit, all first-order derivatives of V eff,3 0 vanish at The potential V eff,3 0 may have additional stationary points; we will return to the possibility of such minima below.
In order to constrain the allowed parameter space, we consider the value of the potential at the field values given in eq. (2.43) and demand the physical minimum to be the global minimum. As we will see, the |µ| vs. v S /v S plane is a useful projection of the parameter space. In the alignment limit, the potential at the physical minimum takes the value where we used eqs. (2.28) and (2.42) for the second equality. We can derive a first constraint on the parameter space by demanding the physical minimum to be deeper than the trivial minimum. The scalar potential vanishes at the trivial minimum, V eff,3 0 (0, 0, 0) = 0. Thus, in the alignment limit, demanding Comparing with eq. (2.45), we see that this stationary point of the potential is never deeper than the physical minimum; H SM , H NSM , H S = 0, 0, √ 2µ/λ is a saddle point of the scalar potential in the alignment limit.
On the other hand, at H SM , H NSM , H S = 0, 0, √ 2v S , the scalar potential (in the alignment limit) takes the value Demanding this minimum to be shallower than the physical minimum, V eff,3 defining a range of v S /v S for which the physical minimum is deeper than the minimum at H SM , H NSM , H S = 0, 0, √ 2v S .
As we noted above, the potential may feature additional stationary points beyond those listed in eq. (2.43). In particular, minima deeper than the physical minimum can easily appear in the NMSSM for field configurations where H NSM and H SM take non-zero vevs. Such minima break the electroweak symmetry, and, unless H NSM = 0 and H SM = √ 2v, do not lead to electroweak physics compatible with observations. In general, V eff, Such field configurations are very challenging to identify analytically, thus, we resort to numerical techniques to infer the constraints on the NMSSM parameter space arising from demanding the physical minimum to be deeper than any minima where H NSM = 0. 6 Finally, the parameter space of the NMSSM is also constrained by avoiding tachyonic 6 We use the package HOM4PS2 [88]   600 masses. As discussed in section 2, the most relevant constraint arises from avoiding the singlet-like neutral CP-odd state, a S , becoming tachyonic. In terms of v S /v S , the constraint arising from eq. (2.24) can be rewritten as . Correspondingly, these constraints come from avoiding the trivial minimum, the minimum at H SM , H NSM , H S = 0, 0, √ 2v S , or minima with H NSM = 0, becoming deeper than the physical minimum. We also show the region where the singlet-like CP-odd mass eigenstate a S becomes tachyonic, eq. (2.50), with the gray shade. Note that overlapping regions are marked by dashed lines of the corresponding colors. In all figures, we truncate the x-axis at |µ| = 100 GeV; smaller values of |µ| are disfavored by null results of chargino searches at LEP. Since we imposed alignment (without decoupling), the scalar potential is uniquely specified by v S /v S (see eq. (2.42)), µ, tan β, and κ/λ, and the potential is insensitive to the sign of µ. As we can see from eqs. (2.46) and (2.49) (the green and blue shaded regions, respectively), the conditions stemming from the trivial minimum and the minimum at H SM , H NSM , H S = 0, 0, √ 2v S becoming deeper than the physical minimum do not depend on the sign of κ and are relatively insensitive to the value of |κ|.
For large values of κµ 2 , eq. (2.46) implies that the physical minimum is deeper than the trivial minimum for v S /v S < 1/2. Equation (2.49) on the other hand implies, for large κµ 2 , that v S /v S > −1 to avoid the minimum at H SM , H NSM , H S = 0, 0, √ 2v S becoming deeper than the physical minimum. These constraints relax for smaller values of 600 With λ ∼ 0.65 in the alignment limit, we find λm h 125 v ∼ 120 GeV. In figures 1-3 we can observe the corresponding change in the blue and green shaded bounds for |µ| 120 GeV/ |κ/λ|.
Finally, as discussed above, the region of parameter space where minima with H NSM = 0 are deeper than the physical minimum can only be inferred by numerically investigating the vacuum structure. From figures 1-3 we see that such constraints become more stringent with larger |κ| and depend on the value of tan β. Furthermore, the constraints arising from avoiding such minima are sensitive to the sign of κ; for κ < 0, avoiding minima with H NSM = 0 effectively sets a lower limit on the value of v S /v S , while for κ > 0, avoiding these minima sets an upper bound on the value of v S /v S .

Thermal History: Analytical Understanding
In this section, we explore the possible phase transition patterns in the NMSSM. We first discuss the effective potential at very high temperatures, which gives guidance on the starting point of the thermal evolution. Then, we discuss the requirements a phase transition must satisfy to provide favorable conditions for electroweak baryogenesis via a SFOEWPT. We continue by discussing specific phase transition patterns which appear in the NMSSM, and fix a shorthand notation we will use to identify them. We close this section by discussing the regions in parameter space where we expect to observe different transition patterns, in particular, the regions in which we expect the nucleation probabilities of first order phase transitions to be sufficiently large for such transitions to complete.
Let us start with the vacuum structure at very high temperatures. In the limit T 2 m 2 i , and neglecting the Daisy coefficients, the finite temperature potential, eq. (2.34), can be written as The ellipsis [. . .] in eq. (2.51) indicates terms which are independent of the field values. It is straightforward to see that in this limit, the field-dependent terms of the thermal potential are parameterized by the Daisy coefficients, (see eq. (C.2)) where the ellipsis now includes both the field-independent and higher-order terms. Explicit expressions for the c ij can be found in appendix C. Note that the symmetries of the NMSSM enforce this particular form of the high-temperature potential. In particular, the A phase transition must satisfy certain requirements in order to give rise to favorable conditions for electroweak baryogenesis: In order for a baryon asymmetry to be produced in the transition, and such asymmetry not to be subsequently washed out in the low temperature phase, electroweak sphalerons must be active in the high-temperature phase and suppressed in the low temperature phase. Estimating the rate of the sphaleron suppression is a notorious problem in the perturbative approach to the phase transition calculation, see, for example, refs. [3,89,90], and even more so if the electroweak symmetry is broken in multiple steps, see, for example, ref. [91].
We shall demand as conditions for a SFOEWPT. Here, Φ hT ( Φ lT ) is the value of Φ in the high (low) temperature phase at the temperature T where the phase transition occurs. The first condition ensures that electroweak sphalerons are inactive in the low-temperature phase, while the second condition requires the sphalerons to not be unduly suppressed in the high temperature phase. We stress that while the numerical thresholds for the order parameters we chose in eq. (2.53) are indicative for the possibility of generating the baryon asymmetry through a SFOEWPT [89], obtaining the exact conditions would require a gauge-invariant evaluation of the sphaleron profile through the bubble wall which is beyond the scope of this work.
In the remainder of this paper, we use a shorthand notation to classify the phase transition patterns we observe in the NMSSM: • We use an integer (1, 2, . . .) to denote the number of steps in the transition patterns.
• For 2-step transitions (we don't observe transition patterns with more than 2 steps in our data) we use a roman number to classify the intermediate phase: • We use a lower case letter to denote the strength of any transitions in which electroweak symmetry is broken in the low-temperature phase, -"a" denotes a SFOEWPT, -"b" denotes a first order phase transition that is not a SFOEWPT, i.e does not satisfy one (or both) of the conditions in eq. (2.53), -"c" denotes a second order phase transition.
Thus, for example, "1-a" denotes a direct one-step SFOEWPT from the trivial phase to the electroweak phase. "2(I)-b" denotes a two-step transition pattern, where the first step is from the trivial phase to a singlet-only phase (since electroweak symmetry is not broken in this intermediate phase, we do not differentiate the pattern with respect to the strength of this first transition), and the second step is a first order (but not SFOEWPT) transition from the singlet-only to the electroweak phase. "2(II)-ca" on the other hand denotes a twostep phase transition pattern, where the first transition is a second order phase transition into a phase in which electroweak symmetry is broken (but which is distinct from the electroweak phase), and the second transition is a SFOEWPT from this intermediate phase to the electroweak phase. We can get some intuition about the different regions of parameter space suitable for the respective phase transition patterns from the shape of the effective potential. While thermal effects alter the shape of the potential at finite temperatures, the zero-temperature vacuum structure still indicates the relative importance of the different possible local minima for the thermal history. Thus, we expect the results from section 2.3 to be indicative for the transition patterns suggested by the critical temperature calculation. For example, we can expect direct one-step transition patterns to most prominently be realized in the parameter region close to where the trivial minimum becomes the global minimum at zero temperature (green shade in figures 1-3). Similarly, we can expect "2(I)" transition patterns to appear in the parameter regions adjacent to where H SM , H NSM , H S = 0, 0, √ 2v S becomes the global minimum at zero temperature (blue shade), and "2(II)" transitions are expected to appear in regions close to those where the global minimum has non-trivial vev of H NSM = 0 (orange shade).
The vacuum structure gives however little information about the tunneling probability from one local minimum to another, i.e. if a first order phase transition suggested by the critical temperature calculation can actually nucleate. The tunneling rate is controlled by the height of the barrier and the distance (in field space) between the respective local minima. The higher the barrier, and the larger the distance between the minima, the lower the nucleation probability. Although the shape of the potential is modified by thermal effects, we can learn some lessons from the zero-temperature potential. As discussed above, the trivial minimum is the global minimum of the effective potential at very high temperatures. Thus, any phase transition pattern starts at H SM = H NSM = H S = 0. The distance between the trivial and the physical minimum (at zero temperature) is given by 2v 2 + 2µ 2 /λ 2 . Since the values of v = 174 GeV and λ ∼ 0.65 are fixed by electroweak precision data and the alignment conditions, respectively, the distance between the trivial and the physical minimum is controlled by |µ|. The distance increases with the value of |µ|, hence, nucleation proceeds more easily for small |µ|.
The height of the barrier around the trivial minimum can be inferred from the squared mass parameters of the fields H SM , H NSM , and H S around the trivial point, i.e. the fielddependent masses given in appendix A at H SM = H NSM = H S = 0. In order for a phase transition to occur, the smallest of the eigenvalues of the squared mass matrix should be approximately zero, implying a flat direction around the trivial point at zero temperature. If the smallest eigenvalue is too large, the barrier around the trivial minimum is large, and hence the tunneling rate will be too small to allow for successful nucleation. If the smallest squared mass eigenvalue is negative, the trivial minimum is a saddle point of the potential (at zero temperature). Finite temperature effects can still give rise to a barrier between the trivial and the physical minimum required for a SFOEWPT in this situation, but only if the absolute value of the smallest squared mass parameter is not too large, such that thermal effects can overcome the zero-temperature shape of the potential. Hu is the smaller of the doublet-like eigenvalues. In the alignment limit, The alignment conditions enforce sizable values of λ ∼ 0.65, thus, the value of m 2 S is controlled by κ 2 µ 2 (v S /v S ). Since the temperature corrections to m 2 S , eq. (C.7), are of order 0.2 T 2 , one would expect that at the characteristic temperature of the EWPT of order 100 GeV, the tunneling rate could only be large enough for successful nucleation if the squared mass parameter controlling the barrier m 2 S (100 GeV) 2 . This condition can be achieved in two ways: either, |v S /v S | 1, or |κµ| 100 GeV. Note that the conditions m 2 Hu ≈ 0 or m 2 S ≈ 0 are indicative for the possibility of a first order phase transition to successfully nucleate at finite temperature since they imply the presence of an approximately flat direction around the trivial minimum at zero temperature. However, this analysis does not predict the transition pattern, which is determined by the shape of the potential away from the trivial minimum (at the transition temperature). The bounce solution of the fields (the trajectory in field space connecting the local minima) is, in general, not a straight line in field space; in particular, m 2 S ≈ 0 does not necessarily lead to "2(I)" transition patters, and m 2 Hu ≈ 0 does not directly imply "2(II)" patterns.

Numerical Results
In order to explore the EWPT in the NMSSM, and, in particular, find which regions of parameter space give rise to phase transition patterns suitable for electroweak baryogenesis, we perform an extensive numerical study using CosmoTransitions v2.0.5 [48]. In this section, we first describe our implementation of the NMSSM in CosmoTransitions and sketch the steps of the calculations CosmoTransitions performs. As discussed in section 2, in the alignment limit, the Higgs sector of the NMSSM can be described by the four parameters {tan β, κ/λ, µ, v S /v S }, and we perform random scans in this parameter space. We show the results of our numerical scans in figures 4-10. In section 3.1 we discuss the regions of the parameter space where points satisfy the boundary conditions we implement in our CosmoTransitions calculation. In section 3.2 we discuss the phase transition patterns suggested by the critical temperature calculation and we compare these results with the thermal histories obtained by calculating the nucleation rate. As we shall see, the phase transition patterns obtained from the nucleation calculation differ substantially from those indicated by the critical temperature calculation, and thus, computing only the critical temperatures provides a misleading picture of the regions of parameter space favorable for electroweak baryogenesis. In section 3.3 we comment on the collider and dark matter phenomenology in the region of parameter space promising for baryogenesis via a SFOEWPT.
The CosmoTransitions package provides a framework for calculating phase transitions in single-and multi-field models (in the perturbative approach). The implementation of a model into CosmoTransitions proceeds via the specification of the effective potential. We have described the effective (temperature-dependent) potential of the NMSSM in section 2, in particular, it consists of the terms is the tree-level potential of the effective theory obtained after integrating out the sfermions and gluinos, V CW 1−loop is the Coleman-Weinberg potential (including counterterms as shown in eq. (2.32)), and V T =0 1−loop contains the thermal corrections to one-loop order, see section 2.2. Explicit formulae for the field-dependent masses, the counterterm coefficients, and the Daisy coefficients are collected in appendices A, B, and C, respectively.
The calculation of the phase transition pattern with CosmoTransitions proceeds in multiple steps: • First, we compute the locations of the local minima at zero temperature 7 .
• Second, the phases, i.e. the temperature-dependent locations in field space and values of the effective potential at the local minima, are computed from the list of zerotemperature minima. Note that if a phase ends at some temperature, i.e. ceases to be a local minimum, CosmoTransitions tries to find other local minima nearby in field space and then traces the corresponding phases as well. Thus, CosmoTransitions attempts to include phases which cannot be obtained from the list of zero-temperature minima because they exist only at finite temperatures.
• Third, using the phases as input, CosmoTransitions analyzes the temperaturedependent vacuum structure of the potential. The most relevant output from this step is a list of critical temperatures, the temperatures at which two distinct local minima of the potential have the same potential value. At the critical temperatures, the role of the global minimum of the effective potential passes from one phase to another, suggesting the phase transition pattern.
• Finally, for possible first order phase transitions indicated by the analysis of the vacuum structure, CosmoTransitions allows to compute the probability of the transition taking place. First order phase transitions proceed via bubble nucleation, and the nucleation rate is commonly parameterized via the bounce action S E , the Euclidean space-time integral over the (effective) Lagrangian density. In practice, it typically suffices to compute the three-dimensional effective Euclidean action, S E S 3 /T . The technically most challenging part of this computation is finding the bounce solution for the scalar fields, i.e. the trajectory in field space connecting the two local minima which minimizes the Euclidean effective action. 8 The bubble nucleation rate per unit volume at finite temperature T is given by Γ/V ∝ T 4 e −S 3 /T ; requiring the nucleation probability (for the EWPT) to be approximately one per Hubble volume and Hubble time leads to the nucleation condition [95] (see, e.g., ref. [12] for a review), The nucleation temperature T n is the (highest) temperature for which S 3 /T 140. If S 3 /T > 140 for all T > 0, the corresponding transition does not occur because the tunneling probability through the barrier separating the respective local minima is too small. Typically, this is caused by a too high barrier and/or a too large distance (in field space) between the local minima.
Since the calculation of the nucleation temperature (involving the computation of the bounce action) is numerically expensive, to date such calculations have only been presented for a few benchmark points in the NMSSM, see refs. [33,36,40,43,45]. Here, we present results based on the full nucleation calculation for a broad scan of the parameter space.
We focus our study on the region of parameter space where alignment without decoupling is realized, i.e. the region of parameter space for which the NMSSM features a Higgs mass eigenstate which (at tree-level) couples to SM particles like the SM Higgs boson. As discussed in section 2, the alignment conditions fix the values of λ and M 2 A (or, equivalently, A λ ), leaving {tan β, µ, κ, A κ } as the four free parameters which control the effective potential. We fix the mass and mixing parameters of the stop sector (parameterized by the threshold correction ∆λ 2 in V eff 0 , see section 2.1) to obtain m h 125 125 GeV for the mass of the SM-like Higgs boson. As discussed in section 2.3, we use v S /v S to re-parameterize A κ . Here, v S = µ/λ is the vev of the CP-even singlet interaction state at the physical minimum, H S = √ 2µ/λ, and v S = − (µ/λ + A κ /2κ) is the location of an extremum of V eff 0 in the singlet-only direction, H SM , H NSM , H S = 0, 0, √ 2v S . In summary, we use as input parameters for our numerical evaluation. Note that throughout our calculations, we fix the bino and wino mass parameters, which enter the radiative corrections from the charginos and neutralinos (see eqs.  We will return to this issue in section 3.3. Note that searches for neutralinos and charginos at the LHC do not constrain the parameter space for |µ| 100 GeV in a relevant way, see, for example, ref. [98]. 9 Admixtures of H NSM and H S of this size modify the production cross sections and branching ratios of h125 by 10 % compared to the SM prediction. The currently best-measured production cross section of the observed Higgs boson is via the gluon-fusion mode with a 1σ uncertainty of ∼ 15 % [71,72]. Similarly, the largest branching ratios of the observed Higgs bosons are measured with ∼ 15 % uncertainty [71,72].
• We check that, at zero temperature, the physical minimum is the global minimum of the effective potential. 10 For each point satisfying all boundary conditions, we compute the phase transition pattern with CosmoTransitions as discussed above. Figures 4-9, to be discussed in detail in sections 3.1 and 3.2, show the results from our parameter scans for tan β = {1.5, 3} and κ/λ = {−0.1, 0.1, 0.3} in the |µ| vs. v S,CW /v S plane; these are the same slices of parameter space for which we have shown constraints from the zero-temperature vacuum structure of the effective tree-level potential, V eff 0 , in figures 1-3. In order to compare the results of the respective calculations, we color-code the points according to the transition patterns indicated by the critical temperature calculations in the left panels of figures 4-9, while in the right panels, points are color-coded according to the thermal history obtained from the full nucleation calculation; see section 2.4 for our shorthand notation of the phase transition patterns. Points violating the boundary conditions described above are labeled "failed BC" in figures 4-9. Points which satisfy all boundary conditions, but for which CosmoTransitions fails to return a phase transition pattern starting from the trivial minima at high temperature and ending in the physical minimum at zero temperature are labeled "no transitions". Note that the left and right panels show the same set of points in parameter space, the only difference is the color-coding of the points.

Boundary Conditions
Let us begin the discussion of the results of our parameter scans with the regions of parameter space where points fail to satisfy the boundary conditions. The boundary conditions are independent of the thermal calculation, hence, the same points are labeled "failed BC" in the left and right panels of figures 4-9.
We This behavior can largely be understood from the discussion of the zero-temperature vacuum structure in section 2.3, see also figures 1-3. The analysis of the vacuum structure in section 2.3 was based on V eff 0 , the potential of our effective model after integrating out all sfermions and the gluinos, but prior to including the Coleman-Weinberg corrections. We indicate the region of parameter space for which, per the analysis in section 2.3, the physical minimum is the global minimum of V eff 0 at zero temperature with the thin black contours in figures 4-9. Since these contours are derived from V eff 0 , the y-axis for these contours is v S /v S , where v S = − (µ/λ + A κ /2κ) is the tree-level value. We see that, although these contours are derived from V eff 0 , they describe well many of the features of the boundary conditions seen in our parameter scan, which incorporates radiative corrections. The largest deviations appear for κ/λ = −0.1, see figures 4 and 5. While the contours here allow only a narrow range of v S /v S values, we see that the points from our parameter scan satisfy the boundary conditions for a much wider range of values of v S,CW /v S than what the contours suggest. Comparing with figure 1, we see that this discrepancy occurs in regions of parameter space where the analysis of V eff 0 suggested that a minimum with H NSM = 0 was the global minimum of the potential (indicated by the orange shade in figure 1). This constraint was derived numerically in section 2.3, and hence is challenging to understand quantitatively. However, it is not surprising that the region of parameter space disfavored by vacua with H NSM = 0 becoming the global minimum of the potential changes considerably after including the Coleman-Weinberg corrections: the potential is subject to larger radiative corrections in the doublet-like directions of the effective potential than in the singlet-like direction, and furthermore, the H NSM direction is affected by the counterterms we have included to maintain the location of the physical minimum, H SM , H NSM , H S = √ 2 {v, 0, µ/λ}. Before moving to the discussion of the phase transition patterns we observe for points satisfying the boundary conditions in section 3.2, let us briefly mention a few features visible in figures 4-9. First, we can see a gap in the points around v S,CW /v S ≈ 1, which widens for small values of |µ|. This gap is due to numerical difficulties in our algorithm to find v S,CW if v S,CW ≈ v S . Identifying the value of v S,CW is particularly challenging for small |µ|, because |µ| controls the size of v S = µ/λ.
Second, an arc of points failing the boundary conditions crosses the region of parameter space consistent with the physical vacuum being the global minimum at zero temperature, starting at small values of |µ| and negative v S,CW /v S and ending at larger values of |µ| and positive v S,CW /v S . This feature is particularly pronounced for tan β = 1.5, and is due to the mass parameters of the interaction states H SM and H S becoming approximately degenerate for those points. As discussed below eq. (3.5), in this situation, even small deviations from the alignment conditions lead to a sizable H S component of h 125 , and thus, these points are forbidden by our requirement C S h 125 < 0.1. Neither of these issues is related to the thermal history of a given parameter point, and these issues do not occur in regions of parameter space which are of special interest for the phase transition calculation. Hence, we ignore them in the following.
We also note that in the left panels of figures 4-9, where we show the results of the critical temperature calculation, points labeled "no transition" appear. As discussed in section 2.4, the trivial minimum is guaranteed to be the global minimum of the potential at high temperatures, and for any point passing the boundary conditions, the physical minimum is the global minimum at zero temperature. For points labeled "no transition", CosmoTransitions failed to return a transition pattern starting in the trivial minimum at high temperatures and ending in the physical minimum at zero temperatures. This is due to numerical errors arising in the second step of the numerical calculation described above, i.e. the step in which CosmoTransitions attempts to trace the local minima of the effective potential with changing temperatures. We have investigated these numerical issues, and have not found any indication that they bias our results towards particular regions of parameter space. Thus, we expect that our scanning over a large number of points throughout the parameter space gives an accurate picture of the regions of parameter space suitable for electroweak baryogenesis.

Comparison of Critical Temperature and Nucleation Results
In this section, we compare the phase transition patterns obtained from the nucleation calculation with the ones suggested by the analysis of the temperature-dependent vacuum structure at the critical temperatures. In figures 4-9, the color-coding of the points in the left panels shows the phase transition patterns suggested by the critical temperature calculation. In the right panels of figures 4-9, we color-code the points according to the thermal histories obtained from the nucleation calculation. Comparing the left and right panels, we see that the thermal histories obtained from the nucleation calculation differ significantly from those the critical temperature analysis suggests, leading to a marked shift in the regions of parameter space which allows for a SFOEWPT.
Let us begin by discussing the results for tan β = 1.5 and κ/λ = −0.1, shown in figure 4. For the critical temperature results, shown in the left panel, we observe that onestep SFOEWPT patterns ("1-a", dark green points) occur at the upper range of the values of v S,CW /v S allowed by the boundary conditions, and that the range of v S,CW /v S for which we find such "1-a" transition patterns becomes wider for smaller values of |µ|. For smaller values of v S,CW /v S and larger values of |µ|, we find two-step transition patterns where the intermediate phase is in the singlet-only direction ("2(I)", blue points). However, except for a few "2(I)-a" points at values of µ 250-300 GeV and small values of v S,CW /v S , the EWPT for these points is weakly first order ("2(I)-b") or a second order transition ("2(I)-c") as indicated by the lighter blue shades of the points.
Qualitatively, the patterns suggested by the critical temperature calculation can mostly be understood from the discussion of the zero-temperature vacuum structure in section 2.3. The left panel of figure 1 shows the different constraints on the zero-temperature vacuum structure (at tree level) for the same slice of parameter space as figure 4. At large values of v S /v S , the trivial minimum is deeper than the physical minimum, indicated by the green shade in figure 1. Thus, towards large v S,CW /v S , we expect the trivial minimum to play a large role in the thermal history, and accordingly, we find one-step transitions from the trivial to the physical minimum in this region of parameter space in the left panel of figure 4. Similarly, for small values of v S /v S , the minimum in the singlet-only direction is deeper than the physical minimum (blue shaded region in figure 1), hence, the singlet-only phase plays a larger role in the thermal history, explaining the appearance of "2(I)" transition Focusing now on the results of the nucleation calculation, we should recall that electroweak baryogenesis requires a SFOEWPT, i.e. one of the phase transition patterns labeled with an "a" in our shorthand notation. The only such patterns we observe for tan β = 1.5 and κ/λ = −0.1 in the right panel of figure 4 are direct one-step transitions ("1-a", dark green points), that occur for a narrow range of values v S,CW /v S ∼ 0. At small values of |µ|, the range of values of v S,CW /v S for which we find SFOEWPTs widens slightly, before being truncated by the boundary conditions. For values of v S,CW /v S just below the "1-a" patterns, we find one-step transitions from the trivial to the physical minimum which are not strong first order ("1-b" and "1-c", lighter green colors). For even smaller values of v S,CW /v S , we find two-step transitions where the intermediate phase is in the singlet-only direction and where the second transition step, in which electroweak symmetry is broken, is weakly first order or second order ("2(I)-b" or "2(I)-c", light blue points). Note that outside of these bands in v S,CW /v S , we do not find points for which the nucleation calculation indicates thermal histories ending in the physical minimum. This should be contrasted with the phase transition patterns suggested by the critical temperature calculation, where we observe "1-a" patterns at much larger values of v S,CW /v S . The The reason for the mismatch between the critical temperature and nucleation results was discussed in section 2.4: While the behavior of the critical temperatures can be understood from the zero-temperature vacuum structure, the nucleation probability is controlled by the height of the barrier separating the local minima, and the distance in field space between the local minima. For all parameter points, the thermal evolution starts in the trivial minimum at high temperatures. For large values of v S /v S , the barriers around the trivial minimum are large, making the tunneling probability prohibitively small. Hence, for larger values of v S /v S , the fields are "stuck" at H SM = H NSM = H S = 0, even if at zero temperature the trivial minimum is no longer the global minimum of the potential as required by the boundary conditions. For v S /v S → 0, the zero-temperature effective potential becomes flat in the singlet direction around the trivial point, and for v S /v S < 0 the trivial point turns into a saddle point of the potential, see eq. (2.55). For small values of |v S /v S |, thermal effects can still give rise to a barrier around the trivial minimum at finite temperatures, while for large negative values of v S /v S , thermal effects can no longer overcome the zero-temperature shape of the potential to give rise to the barrier required for a SFOEWPT. This behavior of the barrier explains why the nucleation calculation singles out the region around v S,CW /v S = 0 for a SFOEWPT in the right panel of figure 4.
For tan β = 3, shown in figure 5, we find similar results as for tan β = 1.5. Beginning with the critical temperature results (left panel), the main difference is that for the larger values of tan β, we observe that two-step transition patterns ("2(II)", orange and magenta points) appear at small values of |µ|. This is somewhat difficult to understand from the analysis in section 2.3. The constraints coming from local minima in the doublet-like directions (orange shade in figure 1) are the only vacuum structure constraints depending on the value of tan β. However, as mentioned in section 3.1, the doublet-like directions are subject to large radiative corrections, explaining the mismatch between the region where "2(II)" patterns appear in our numerical results and the orange shaded region of the treelevel vacuum structure analysis in figure 1. The appearance of the "2(II)" patterns can however be understood from the region of parameter space for which m 2 Hu < 0, eq. (2.54). In section 2.4, this condition was discussed in the context of the zero-temperature barrier in the H u -direction disappearing for m 2 Hu 0, leading to large tunneling rates at finite temperature. To understand the critical temperature results, it is more relevant to note that for m Comparing the nucleation calculation results for κ/λ = −0.1 and tan β = 1.5 with those for tan β = 3, shown in the right panel of figures 4 and 5, respectively, we see that the preferred region of parameter space for a SFOEWPT is almost independent of the value of tan β. The main difference is that for tan β = 3, points with smaller values of |µ| satisfy the boundary conditions, leading to the band of points around v S,CW /v S = 0 for which we find SFOEWPTs ("1-a", dark green points) extending to lower values of |µ| than for tan β = 1.5. For tan β = 3, we also see the emergence of two-step transition patterns, where electroweak symmetry is broken in the intermediate phase, ("2(II)", orange and magenta points) for positive values of v S,CW /v S and small values of |µ|. As discussed around eq. (2.54), for small values of |µ|, the barrier around the trivial point in the H u direction disappears. Note however that these points (except for one parameter point at v S,CW /v S ∼ 0) do not feature a SFOEWPT step, but both steps are weakly first order or second order.
Let us now discuss the results for κ/λ = 0.1, shown in figures 6 and 7 for tan β = 1.5 and tan β = 3, respectively. Comparing the κ/λ = −0.1 critical temperature results (left panels) with those for κ/λ = 0.1, we find that many of the features remain the same. The two main differences are that the boundary conditions relax for small values of |µ|, allowing a larger range of values for v S,CW /v S , and that for tan β = 3, "2(II)" patterns appear even more prominently in the low |µ| region. The behavior of the boundary conditions is discussed in section 3.1, hence, we focus on the latter difference here. As for the κ/λ = −0.1 case, the appearance of "2(II)" patterns can be understood from the region of parameter space where m 2 Hu < 0. From eq. (2.54), we find that, for tan β = 3 and κ/λ = −0.1, the mass parameter for H u becomes tachyonic for |µ| 230 GeV, while for κ/λ = 0.1, this critical value increases to |µ| 320 GeV. Accordingly, we see that "2(II)" patterns appear for larger values of |µ| for tan β = 3 and κ/λ = 0.1 (left panel of figure 7) than for κ/λ = −0.1 (left panel of figure 5).
Let us now concentrate on the nucleation results for κ/λ = 0.1. For tan β = 1.5, see  For tan β = 3, we likewise find similar behavior for κ/λ = 0.1 (right panel of figure 7) and for κ/λ = −0.1 (right panel of figure 5). Here, the main difference is that for κ/λ = 0.1, two-step transition patterns where electroweak symmetry is broken in the intermediate phase ("2-II") play a larger role than for κ/λ = −0.1, restricting the values for which we find SFOEWPTs to a narrower band of values of v S,CW /v S . This can again be understood from the range of values for which H u becomes tachyonic around the trivial point. Note that the presence of this tachyonic direction in the effective potential (at zero temperature) makes it more difficult to achieve transition patterns favorable for baryogenesis, which we see reflected in the absence of "1-a" transition patterns for |µ| 200 GeV in the right panel of figure 7.
Considering finally the critical temperature results for κ/λ = 0.3 (left panels of figures 8 and 9), we find that compared to the results for smaller values of κ/λ, two-step transition patterns play a much larger role. Comparing eq. (2.45) with eq. (2.48), we see that the depth of the singlet-like minimum is much more sensitive to the value of κ/λ than the depth of the physical minimum, and thus, the minimum in the singlet-only direction plays a larger role in the thermal history for larger values of κ/λ, leading to "2(I)" patterns (blue points) appearing more prominently for κ/λ = 0.3 than for κ/λ = −0.1 and 0.1. Likewise, we see "2(II)" patterns (orange and magenta points) appearing more prominently in the region of parameter space not ruled out by the boundary conditions. For tan β = 1.5 and κ/λ = 0.3, we find from eq. (2.54) that m 2 Hu < 0 (at zero temperature) for |µ| 180 GeV, while for κ/λ = 0.3, the critical value is |µ| 840 GeV.
Regarding the nucleation results, for tan β = 1.5 and κ/λ = 0.3, shown in the right panel of figure 8, we find SFOEWPTs for small values of |µ| and v S,CW /v S . The scaling of the depths of the respective local minima with v S /v S becomes faster the larger the value of |κ/λ|, making the change in phase transition behavior with the value of v S,CW /v S more rapid for this larger value of κ/λ than what we have observed for lower values of κ/λ. Thus, the range of v S,CW /v S leading to (one-step) SFOEWPTs is smaller for all values of |µ| than what we found for κ/λ = ±0.1. Furthermore, we observe that "2(II)" transition patterns appear for small values of |µ| due to the disappearance of the barrier in the H u direction. This behavior is even more pronounced for tan β = 3 and κ/λ = 0.3, see the right panel of figure 9. In this slice of parameter space, m 2 Hu < 0 (at zero temperature) for |µ| 840 GeV, and we do not find any parameter points with a SFOEWPT.
We stress that for all slices of parameter space shown in figures 4-9, the region providing favorable conditions for electroweak baryogenesis via a SFOEWPT differs markedly when the thermal history is inferred from the nucleation calculation instead of the simpler calculation of studying only the vacuum structure at the critical temperatures. While the critical temperature results can be explained from the zero-temperature vacuum structure, the regions of parameter space where SFOEWPTs actually nucleate can only be understood when considering the barriers of the effective potential. We find that SFOEWPTs can only nucleate if v S,CW /v S 1 and |κ/λ| is not too large, leading to a small barrier in the singlet direction. If |κµ| is significantly smaller than the weak scale, larger val-  ues of v S,CW /v S can still lead to a small barrier in the singlet direction and a successful SFOEWPT. For larger values of κ/λ and tan β, the barrier in the H u direction disappears in the small |µ| region, leading to multi-step phase transition patterns where electroweak symmetry is broken in the intermediate phase, and typically, no SFOEWPT is realized.
In figure 10, we collect the results of our scans over the different slices of parameter space shown separately in figures 4-9. As before, we classify points based on the thermal histories suggested by the critical temperature calculation in the left panels, while in the right panels, parameter points are color-coded according to the results of the nucleation calculation. In order to highlight the region of parameter space for which the respective calculations indicate a SFOEWPT, we show only the points falling in one of the transition patterns "1-a", "2(I)-a", "2(II)-aa", "2(II)-ab", "2(II)-ac", "2(II)-ba", or "2(II)-ca" in figure 10. In the upper panels, we show results in the |µ| vs. v S,CW /v S plane. Comparing the left and the right panels, it is evident that the critical temperature calculation gives a misleading picture of the parameter space favorable for electroweak baryogenesis. We note also that a one-step SFOEWPT ("1-a", green points) is by far the most generic possibility to realize a SFOEWPT in the NMSSM. While multi-step transitions including a SFOEWPT step can occur in the NMSSM, our results suggest that such transition patterns require very particular combinations of parameters, making them rare in a (random) parameter scan.

Collider and Dark Matter Phenomenology
In this section we discuss the prospects for collider searches to cover the region of parameter space where we find SFOEWPTs and comment on the possibility of realizing a dark matter candidate in this parameter space.
In the lower panels of figure 10, we show the points from our parameter scans for which we find a SFOEWPT in the plane of the masses of the two non-SM-like neutral CP-even Higgs bosons. Recall that we denote the state with the largest H S component by h S , and the state with the largest H NSM component by H. Comparing the left and the right panels, we observe that, similar to what we saw in the |µ| vs. v S,CW /v S plane, the results based on the full nucleation calculation lead to a considerably tighter relation between m H and m h S for points with SFOEWPTs than the results of the critical temperature calculation, as well as a significant shift of the preferred region of parameter space. As we have seen above, SFOEWPTs occur in the region of parameter space where |v S /v S | 1, or |κµ| 100 GeV. In this limit, the mass of the singlet-like mass eigenstate (at tree-level and in the alignment limit), is approximately given by while the mass of the doublet-like mass eigenstate is approximately  While we leave a study of the collider phenomenology of the region of parameter space where we find a SFOEWPT for future work, we can make some broad statements. As we have seen in section 3.2, see also figures 4-10, SFOEWPTs can be realized in the NMSSM for small values of |κ/λ| and tan β, and not too large values of |µ|, leading to relatively light non-SM-like Higgs bosons. From eq. (3.7) we find that the doublet-like state can be as light as m H ∼ 200 GeV for tan β = 1.5 and |µ| ∼ 100 GeV, as shown in the lower right panel of figure 10. Similarly, the singlet-like state can be as light as m h S ∼ 70 GeV for tan β = 3, κ/λ = 0.1, and |µ| ∼ 100 GeV. Despite the relatively small masses, this region of parameter space is challenging to probe at colliders: The direct production cross section of the singlet-like state is suppressed by its small doublet component, we find C NSM h S 10 % for the points featuring a SFOEWPT. The doublet-like state H has sizable production cross sections. However, its decay patterns make it challenging to probe for the small values of tan β preferred by a SFOEWPT. Considering the decays into pairs of SM fermions, due to the small value of tan β, the decay mode into top-quark pairs will dominate if kinematically accessible. Thus, for m H 350 GeV, the branching ratio into pairs of top quarks will be large and this final state is very challenging to probe at hadron colliders such as the LHC [99][100][101][102][103][104][105][106]. For m H 350 GeV on the other hand, the H → τ + τ − channel could provide some sensitivity. However, due to the preference for small values of |µ| and |κ/λ|, the Higgsinos and singlinos will be relatively light; their mass parameters are µ and 2κµ/λ, respectively. Thus, decays of H into pairs of neutralinos will be kinematically allowed in the parameter region preferred by a SFOEWPT, and the associated branching ratios will be sizable, suppressing H → τ + τ − decays. The final states arising from decays of H into neutralinos are challenging to probe at the LHC, see, for example, refs. [98,104,[107][108][109][110]. Out of the di-boson final states, decays of H into two SM(-like) states, e.g. h 125 h 125 , ZZ, and W + W − will be strongly suppressed due to alignment [49,111]. However, the branching ratios into final states containing at least one singlet-like boson, such as h 125 h S or a S Z, will be sizable if kinematically allowed [49,69,70,[111][112][113][114][115][116][117], making these channels a promising means to explore the region of parameter space preferred for a SFOEWPT.
Considering the neutralino sector, we find that the region of parameter space where a SFOEWPT is realized features light singlinos. However, a singlino-like neutralino is only a good dark matter candidate if its spin-independent cross section is suppressed by the so-called blind-spot cancellations, see, for example, refs. [20,22,28]. For a singlinolike dark matter candidate, the blind-spot condition in the NMSSM is 2κ/λ ≈ sin 2β, requiring larger values of κ/λ or tan β than those for which we find SFOEWPTs. On the other hand, the value of the bino mass parameter M 1 has practically no influence on the SFOEWPT 11 . Thus, the most promising dark matter scenario in the region of parameter space where we find SFOEWPTs is a bino-like lightest neutralino. The interaction cross sections of such a bino-like neutralino can be sufficiently small to be compatible with the null results from direct detection type experiments without requiring additional (blindspot) cancellations [28,118]. However, its couplings are too small to provide the correct dark matter relic density via standard thermal production. For |M 1 | m t , the correct relic density for a bino-like lightest neutralino can be achieved via co-annihilation with the singlino-like neutralino in the so-called new well-tempered scenario, where |M 1 | ∼ |2κµ/λ| [28]. The bulk of the region of parameter space where we find SFOEWPTs features smaller values of |µ|. There, the correct relic density for a bino-like lightest neutralino could be achieved via resonant annihilation through the singlet-like CP-even or CP-odd states, h S or a S , requiring the mass of the lightest neutralino χ 1 to satisfy m χ 1 m h S /2 or m χ 1 m a S /2, respectively. Alternatively, the NMSSM neutralinos may be unstable (on cosmological scales) and the dark matter may be provided by particles not included in the NMSSM, like axions and axinos (see, for example, ref. [119]).

Summary and Conclusions
Electroweak baryogenesis is a compelling scenario for the generation of the baryon asymmetry of the Universe. It relies on the presence of a Strong First Order Electroweak Phase Transition (SFOEWPT). The Standard Model (SM) of particle physics does not provide appropriate conditions for electroweak baryogenesis, hence, new physics is needed for realizing a SFOEWPT. Calculating the phase transitions in models of new physics is numerically expensive, and hence, most studies in the literature content themselves with studying the vacuum structure at the critical temperatures. At the critical temperature, the role of the global minimum of the potential passes from one local minimum to another, hence, this calculation ensures that a necessary condition for a first order phase transition is met. However, the critical temperature calculation does not ensure that the (quantum-mechanical) tunneling rate through the barrier separating the false from the true vacuum is large enough for such a first order phase transition to occur. In this work, we have investigated if a more complete calculation including the computation of the nucleation probability is necessary to understand the phase transition patterns in models of new physics. As an example model, we chose the Next-to-Minimal Supersymmetric extension of the Standard Model (NMSSM). We focused our case study of the NMSSM on the region of parameter space where alignment-without-decoupling is realized. For the purposes of the phase transition, the remaining four-dimensional parameter space is well described by the set of parameters κ/λ, tan β, |µ|, and v S /v S , where v S is the vev of the singlet H S at an extremum of the effective potential in the singlet-only direction, and v S is the vev of H S at the physical minimum. 12 Using extensive parameter scans, we have demonstrated that successful nucleation of a SFOEWPT occurs mostly in a narrow region of parameter space where |v S /v S | takes small values, and that the range of v S /v S leading to a SFOEWPT becomes increasingly narrow for larger values of κ/λ, tan β, and |µ|. This region of parameter space differs markedly from what one would have inferred from the critical temperature calculation alone, that, in general, suggests a SFOEWPT for much larger values of v S /v S . The difference between the two results can be understood from the shape of the effective potential. In the region of the parameter space suggested by the critical temperature calculation, the barriers around the trivial minimum, where the thermal evolution of the model begins at very high temperatures, are large, leading to prohibitively small tunneling rates. However, the barrier in the singlet direction diminishes for small values of |v S /v S |, enabling tunneling from the trivial minimum. As we have shown, the requirement on the values for v S /v S loosens for values of |κµ| far below the weak scale. The dependence of the parameter region where we find a SFOEWPT on the value of tan β arises mainly from the disappearance of the barrier in the H u -direction, triggering a phase transition which tends to lead to thermal histories incompatible with electroweak baryogenesis. The barrier in the H u -direction disappears for small values of |µ|, and the range of values of |µ| for which this occurs is broader for larger values of tan β and κ/λ.
Note that our findings are obtained in a perturbative expansion of the effective potential (to one loop, improved by relevant resummations), and, hence, may be affected by the well-known shortcomings associated with this expansions [7,[120][121][122][123][124][125][126][127][128]. Nonetheless, our results stress the relevance of computing the nucleation probability to obtain the regions of parameter space promising for electroweak baryogenesis via a SFOEWPT. Our computations strongly rely on the accuracy of CosmoTransitions, thus, they would profit from corroboration with an independent calculation of the tunneling rate.
While we have focused on the phase transitions, the region of parameter space where a SFOEWPT occurs also leads to interesting collider and dark matter phenomenology. We find masses of the singlet-like state 70 GeV m h S 200 GeV. The mass of the new doubletlike Higgs H, on the other hand, depends more strongly on tan β. At moderate values of tan β, we find m H 350 GeV, and hence, H decays prominently into pairs of top quarks.
For smaller values of tan β ∼ 1.5, H can be lighter than the top pair production threshold. Although in principle this enhances the branching ratio into tau leptons, collider searches in conventional SM decay modes of these non-SM-like Higgs bosons are quite challenging due to the presence of decays into light non-standard Higgs, neutralino, and chargino states. The most promising search channels seem to arise via the so-called Higgs cascade decays, e.g., H → h 125 + h S . We reserve a more detailed study of the collider phenomenology for future investigation.
The preference for small values of κ/λ for a SFOEWPT implies the presence of a light singlino in the spectrum. While the spin-independent cross section of such a singlino is too large to be compatible with the null results from direct detection experiments in the region of parameter space where we find a SFOEWPT, a viable dark matter candidate could be realized via a bino-like lightest neutralino, if its annihilation cross section is enhanced through co-annihilation or resonant annihilation.
In closing, we would like to stress that arguably the most important result of this work is that the nucleation calculation yields qualitatively different results for the phase transition patterns in the NMSSM than what the simpler analysis based only on the vacuum structure at the critical temperatures suggests. While our numerical results are obtained in the NMSSM, we expect similar behavior to appear in other models where multiple scalar fields participate in the EWPT. Our results emphasize that, in order to infer the regions of parameter space where electroweak baryogenesis can be realized, it is critical to compute the thermal histories based on the nucleation probabilities.

A Field-dependent masses
In this appendix, we present explicit expressions for the field-dependent masses after inclusion of the leading stop corrections, but without corrections from the Coleman-Weinberg potential. As argued in section 2, it suffices to study the potential as a function of the three neutral CP-even degrees of freedom H SM , H NSM , H S . Let us begin by presenting the expression for the field-dependent (squared) masses in the scalar sector. These can be directly obtained from the scalar potential, (A.1) The entries involving the CP-even interaction states are or in words, the charged Goldstone mode G ± is massless and decouples from the charged Higgs H ± . The remaining entries of the (symmetric) (10 × 10) matrix of the m Φ i ,Φ j not listed above vanish due to CP-and charge conservation.
The field-dependent masses for the electroweak gauge bosons are given by with the weak mixing angle cos θ W = g 2 / g 2 1 + g 2 2 = m W /m Z . The masses of the vector bosons at the physical minimum are related to the gauge couplings as For the 5 neutralinos, the (symmetric) matrix of field-dependent masses in the basis B, W 3 , H 0 d , H 0 u , S can be written as Finally, the field-dependent mass of the top quark is given by where the Yukawa coupling h t is related to the (running) top quark mass m t via h t = m t /s β v.
We compute the contributions to the Coleman-Weinberg potential as well as to the thermal potential in the Landau gauge. This is useful since in the Landau gauge the ghosts decouple and we do not have to include them in our calculations. The quantities entering the Coleman-Weinberg and the thermal potential are the eigenvalues of the respective mass matrices. Recall that the number of degrees of freedom are n i = 1 for the three neutral CP-even and three neutral CP-odd states, n i = 2 for the two charged Higgs states, n i = 6 for the W ± bosons, and n i = 3 for the Z-boson. Out of the fermions, the top quark has n i = 12 and the five neutralinos have n i = 2 each. Since we wrote the chargino mass matrix, eq. (A.23), in terms of four Majorana states (which combine to two physical Dirac fermions), the four eigenvalues of eq. (A.23) enter with n i = 2 each.

B Counterterm coefficients
In order to maintain the location of the physical minimum at H SM , H NSM , H S = √ 2 {v, 0, µ/λ}, preserve m h 125 = 125 GeV, and M 2 S,13 → 0 (i.e. alignment of H S and H S ) after including the Coleman-Weinberg corrections, we include the counterterms given in eq. (2.32). The solutions for the counterterms to satisfy these conditions are is the effective potential including the Coleman-Weinberg corrections V CW 1−loop at zero temperature, all derivatives are evaluated at the physical minimum, H SM , H NSM , H S = √ 2 {v, 0, µ/λ}, and m h 125 is an input parameter which sets the mass of the H SM interaction eigenstate of the Higgs basis.

C Daisy coefficients
The Daisy coefficients c i for the thermal masses can be obtained from the high-temperature limit of the thermal corrections to the effective potential, Note that for the derivation of the Daisy coefficient, V T =0 1−loop = V T =0 1−loop ( m 2 i ) is computed with the temperature independent field-dependent masses m 2 i , while when computing the temperature-dependent effective potential, the Daisy-resummation improved thermal masses m 2 i are inserted in V T =0 1−loop as well as in the Coleman-Weinberg potential. Note also that while we gave explicit expressions for the m 2 i as a function of the three neutral CP-even Higgs boson interaction states, H SM , H NSM , and H S , in appendix A, when computing the Daisy coefficients via eq. (C.2), the field-dependent masses must be inserted as a function of all bosonic fields, i.e. The non-vanishing coefficients involving the neutral Higgs bosons are and those involving the charged Higgs bosons are Note that the photon gets a temperature-dependent mass. In order to properly account for this appearance of the longitudinal degree of freedom of the photon, the Daisy resummation improved thermal masses of the neutral electroweak gauge bosons must thus be included as the eigenvalues of mass matrix, −g 1 g 2 −g 1 g 2 g 2 1 +T 2 5g 2 2 /2 0 0 13g 2 1 /6 . (C.12) After removing the contribution from the neutralinos and charginos to V T =0 1−loop ( m 2 ), these results agree with the results in ref. [45] (where the neutralino and chargino contribution were neglected).