Strong first-order phase transitions in the NMSSM — a comprehensive survey

Motivated by the fact that the Next-to-Minimal Supersymmetric Standard Model is one of the most plausible models that can accommodate electroweak baryogenesis, we analyze its phase structure by tracing the temperature dependence of the minima of the effective potential. Our results reveal rich patterns of phase structure that end in the observed electroweak symmetry breaking vacuum. We classify these patterns according to the first transition in their history and show the strong first-order phase transitions that may be possible in each type of pattern. These could allow for the generation of the matter-antimatter asymmetry or potentially observable gravitational waves. For a selection of benchmark points, we checked that the phase transitions completed and calculated the nucleation temperatures. We furthermore present samples that feature strong first-order phase transitions from an extensive scan of the whole parameter space. We highlight common features of our samples, including the fact that the Standard Model like Higgs is often not the lightest Higgs in the model.


Introduction
One of the enduring problems in modern physics is the origin of the baryon asymmetry of the Universe (BAU) [1][2][3]. This asymmetry cannot be an initial condition in any cosmology that includes inflation, as that would wash out any initial asymmetry. 1 Therefore baryon asymmetry must be produced; however, as yet there is no experimental confirmation of any production mechanism. Any mechanism that produces the BAU must satisfy three criteria [5]: 1. charge (C) and charge-parity (CP) violation, 2. baryon number (B) violation, and 3. departure from equilibrium.
The Standard Model (SM) has the ingredients to satisfy all three criteria: there is a CP violating phase in the CKM matrix, B is violated through sphalerons which are unsuppressed at high temperature and there could be departures from equilibrium following two phase transitions (PTs) that occur in the SM vacuum as it cools -the electroweak (EW) and the QCD transition. Quantitatively, however, the CP violating phase in the CKM matrix is far too feeble to produce enough baryon asymmetry. Furthermore the two transitions that occur in the SM at high temperature are both crossover transitions rather than first-order phase transitions (FOPTs) and therefore do not provide a large enough departure from equilibrium (see e.g., Ref. [6]). As such one has to look beyond the SM for explanations.
While the origin of the baryon asymmetry is a mystery, its measurement is on a firm foundation. During big bang nucleosynthesis, the baryon asymmetry is an input to the set of Boltzmann equations that govern the production of primordial light elements. Since we can measure some of these primordial abundances (deuterium in particular) with high accuracy, this constrains the baryon asymmetry 2 to be [7] Y B ≡ n B s = 8.2 -9.4 × 10 −11 (95% CL).
Furthermore the baryon asymmetry produces acoustic oscillations in the power spectrum of the cosmic microwave background (CMB) [8]. Observing these oscillations gives an even tighter bound on the BAU, The fact that there is a concordance between these two unrelated measurement approaches is a triumph of cosmology. Along with dark matter and inflation, the origin of the BAU is a powerful cosmological argument for physics beyond the SM. 1 For an exception to this rule of thumb see Ref. [4]. 2 We convert measurements of the photon-baryon ratio to Y B by Ref. [2] Y B ≡ n B s ≈ 1 7.04 Electroweak baryogenesis is a minimal and natural explanation for the origin of the baryon asymmetry in the Universe [3,. It utilizes the electroweak phase transition (EWPT) which is known to have occurred in our cosmic history providing the reheating temperature was not unnaturally small. Although this transition is a crossover in the SM, its character may be modified by the introduction of new weak scale bosons such that the transition becomes a strongly FOPT (SFOPT) and proceeds by bubble nucleation. Such a phenomenon is all the more interesting because it might directly be probed by future gravitational wave detectors [44][45][46][47][48][49].
This mechanism can be in principle realized within supersymmetry. In the Minimal Supersymmetric Standard Model (MSSM) a barrier between the EW symmetric and broken vacuums arises from thermal corrections from stops; however, one requires light stops to catalyze the PT such that it is sufficiently strongly first order [50,51]. This is all but ruled out by LHC constraints on stop masses [52]. Much more attractive is the possibility of the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [53,54] where a light singlet scalar can catalyze a strongly first order EWPT [16,25,55]. Unlike the stop which catalyzes the PT through thermal effects, the singlet can change the potential such that there is a barrier even at zero temperature.
Electroweak baryogenesis was recently considered within the NMSSM [15,[56][57][58][59][60] and it was found that the baryon asymmetry can vary by an order of magnitude depending on whether the singlet acquires a vacuum expectation value (VEV) before or during the EWPT (with a simultaneous transition providing more efficient baryon production) [58]. Furthermore the baryon yield is proportional to the maximal variation of the ratio of the two Higgs VEVs, ∆β, and it was shown in Ref. [15,30,[61][62][63][64] that ∆β can be an order of magnitude larger in the NMSSM compared to the MSSM.
In this work we explore the plausibility of EW baryogenesis within the NMSSM, focusing on the PT and leaving CP violation to future work (see [56,[65][66][67][68][69][70][71][72][73][74][75][76][77] for various approaches to generating CP violation). We consider the case where the superpartners are all heavy enough to have their thermal contributions Boltzmann suppressed during the transition. Thus we can match our model to a two Higgs doublet model plus a singlet (THDMS). We sample the parameter space to find points with an EW SFOPT. For such points, we investigate the phase structure, that is the evolution of the minima of the effective potential as the Universe cools. This investigation includes determining whether the singlet acquires a VEV during or before the EWPT and it also involves calculating the strength of the PT.
As we focus on the third Sakharov condition (a departure from thermal equilibrium), we do not consider explicit or spontaneous CP violation in the Higgs sector. We instead assume that CP violation enters the Higgs sector radiatively, though remain agnostic about the exact source of CP violation and do not examine constraints on complex phases (such as electric dipole moments). This simplification allows us to focus only on PTs between the ground states of CP-even fields, easing the numerical problem of finding vacua of a multifield scalar potential.
The structure of this paper is as follows. In Sec. 2 we introduce the NMSSM and the THDMS, fixing the notation we will use in the paper. Following this, in Sec. 3 we describe the radiative and finite temperature corrections that we include in our analysis. Then in Sec. 4 we outline the procedure we use to determine if a point in the parameter space has a FOPT or not, and if so calculate the critical temperature and transition strength. The results of our scan are presented in Sec. 5 and finally our conclusions are given in Sec. 6.

NMSSM
The NMSSM extends the MSSM particle content by adding one singlet superfield,Ŝ. Here we work in the Z 3 symmetric NMSSM where the µ-term of the MSSM is forbidden and instead an effective µ-term, µ eff = λ S , is generated when the singlet develops a VEV, thus solving the µ-problem of the MSSM. The superpotential is given by where a hat is used for superfields, i, j ∈ {1, 2, 3} are family indices, and we have introduced the SU (2) L dot product, A·B = A 1 B 2 −A 2 B 1 . The discrete Z 3 symmetry is spontaneously broken when the Higgs fields or singlet obtain a VEV. We assume that following the strategies of Ref. [78][79][80] domain wall problems can in principle be avoided without impacting any phenomenology.
Under the SM gauge group where the first two entries inside the parentheses give the representation under SU (3) C and SU (2) L , respectively, while the third entry gives the U (1) Y hypercharges without GUT normalization. There are three contributions to the tree-level Higgs potential of the NMSSM: Here the F -and D-term contributions are where g and g are respectively the SU (2) L and U (1) Y gauge couplings without GUT normalization. Finally, the soft-breaking terms are The couplings λ and κ and the corresponding trilinears, A λ and A κ , are in general complex. Three of the four phases, however, may be removed through field redefinitions of H u , H d and S. Since current LHC limits and the 125 GeV Higgs mass measurements require squarks and gluinos to be TeV-scale, the mass spectrum of the NMSSM must contain a large hierarchy between the SM particles and colored sparticles. Furthermore the states with the largest couplings include both heavy sparticles and light SM particles, i.e., stops and the top quark. Therefore higher-order corrections will always include large logarithms since one cannot choose the renormalization scale Q to simultaneously minimize ln m t /Q and ln M SUSY /Q. This makes it challenging to perform precise calculations when working in the full theory. To improve the precision of our calculations we will integrate out the heavy superpartners and use an effective field theory (EFT) which contains only the light states. This makes it possible run to Q = m t and perform calculations in the EFT which are free from large logarithms.

Matching to the THDMS
Since we want to consider scenarios in which all superpartners are too heavy to impact the PT, we match the NMSSM to a two Higgs doublet model plus a singlet (THDMS), which in this context is an effective field theory of the full NMSSM theory valid below M SUSY . 3 The tree-level potential of a Z 3 symmetric THDMS model is where the couplings λ 7 , m 4 and m 5 may be complex. Two of the three phases, however, may be removed by redefinitions of H u , H d and S, leaving a single complex phase, as in the NMSSM. In (10) we follow the conventions in Ref. [15,[81][82][83]; in particular the We match the NMSSM to the THDMS at the scale M SUSY by identifying the tree-level conditions We furthermore included a dominant one-loop threshold correction to the matching for λ 2 , Although we stated the potential and matching conditions for λ 7 , m 4 and m 5 without loss of generality, we later consider only real, CP conserving parameters. As discussed in Sec. 1 we assume that the CP violation demanded by Sakharov's first condition originates in a different sector of the NMSSM, e.g., the squark sector. Although CP violation must enter the Higgs sector through loops, since we only consider the dominant one-loop corrections in the matching, CP violating phases that may appear outside of the Higgs sector do not enter our calculation. At higher orders, however, we would be forced to consider complex parameters and consequently (as later discussed) PTs involving CP-odd fields. An examination of the potential impact this could have is left for future study. Since we match the NMSSM to a THDMS, our results are also applicable to a subspace of the THDMS, which is well-motivated even in the absence of supersymmetry.

Effective potential at zero temperature
In the R ξ gauge the one-loop corrections to the potential, ∆V , are given by [84] where Q is the renormalization scale, m i are field dependent MS masses and the n i are the numbers of degrees of freedom for field i. The first term sums fluctuations of scalar fields, which at the EW breaking minimum can be separated into physical Higgs bosons and Goldstone bosons, the second term sums transverse and longitudinal massive gauge bosons, the third one scalar gauge boson fluctuations, and the final one fermions. We neglect contributions to the vacuum energy. The numbers of degree of freedom for the particles that we include are n t = n b = 12, n τ = 4 (16) for the real scalar, vector and Dirac fermion fields in our model, where A 0 i , H + i and H − i include the physical Higgs states and the Goldstone bosons. At zero temperature, the minimum of the one-loop potential lies at non-zero values for the Higgs fields, which we refer to as VEVs, and assume may always be written as where v u , v d and v S are real, i.e., we do not consider charge or CP breaking VEVs. 4 As we assume that the VEVs are CP conserving, a tadpole condition forces CP violating phases in the potential to vanish. To construct the field dependent masses appearing in (13), we consider the potential as a function of the fields corresponding to the VEVs, i.e., we consider the h u , h d and s components of the fields, where h u , h d and s are real. The field dependent masses are functions of h u , h d and s. In principle, we could consider variation of the charged and CP-odd fields which cannot all be eliminated by gauge fixing. However, because we consider PTs only between charge and CP conserving vacua, we set charged and CP-odd Higgs fields to zero in the field dependent masses. The expressions for the field dependent masses are given in App. A. The effective potential also contains explicit dependence on the gauge parameter ξ. The physical, gauge-independent content of the effective potential may be found through Nielsen identities [87], which express the fact that at extrema, h, the gauge dependence of the effective potential vanishes, since and thus The location of the extrema, however, are gauge dependent, i.e., ∂h/∂ξ = 0. See e.g., Ref. [84,88] for further discussion of this issue. We work in the ξ = 1 (Feynman) gauge.
The effective potential furthermore depends on a choice of renormalization scale, which could in fact have greater impact than gauge ambiguities [89].

Effective potential at finite temperature
To describe the conditions of the early Universe we need to take into account temperature corrections. We calculate one-loop finite temperature corrections including daisy terms using the Arnold-Espinosa method [90] in the ξ = 1 (Feynman) gauge. The effective potential can be written as a sum of zero temperature and finite temperature pieces The one-loop thermal corrections in the R ξ gauge are [84] where the field dependent masses are the same as those appearing in (13) in the previous section, and the expressions for them are given in App. A. The degrees of freedom, n, are as in (14); we again neglect contributions to the vacuum energy and the thermal functions are Here the upper/lower signs are for bosons/fermions. For m 2 T 2 the thermal functions are exponentially suppressed by a Boltzmann factor. This ensures that the massive supersymmetric particles that we integrated out do not impact the finite temperature corrections.
The daisy terms are where we sum over the Higgs fields (including Goldstone bosons) and massive gauge bosons, andm 2 are field dependent mass eigenvalues that include Debye corrections to the tree-level masses in the mass matrices. The Debye corrections add additional T dependent terms of the form c Φ T 2 |Φ| 2 for all complex scalar gauge eigenstates and c A T 2 A µ A µ for all gauge bosons associated with the original gauge symmetries before EWSB. For the THDMS we find, where the couplings g , g, y t , y b and y τ are as in (41). The corrections for the gauge bosons are in the gauge basis before symmetry breaking and every component of a gauge representation receives the same Debye correction. The coefficients are gauge independent, as they originate from a high-temperature expansion of (22) in which the dependence on ξ cancels, We cross-checked our results in (25) against general expressions in Ref. [91] and results for the coefficients in the THDM in Ref. [92]. Thus we have described the full finite temperature potential, which is a function of the fields h u , h d and s and the temperature, T .

First-order phase transitions
Having constructed the finite temperature effective potential, we investigated whether there was a FOPT in which the vacuum of the potential changed abruptly as the Universe cooled. For such a transition to occur, the potential must exhibit two minima separated by a barrier. The temperature at which the two minima are exactly degenerate is known as the critical temperature. That is, at the critical temperature, where caligraphic fonts, h u etc, indicate a minimum of the scalar potential, i.e., Below the critical temperature, the potential develops a minimum that is deeper than the other minima. The system may tunnel through the barrier to the new vacuum state with the lower minimum [93][94][95]. As discussed below, however, the transition might not complete.
We developed a C++ program, PhaseTracer, to map the temperature dependence of the minima of the effective potential and to find potential PTs between them. It enhances the algorithm that was developed in CosmoTransitions [96] to map out the phase structure, and to find out possible PTs between every phase. The numerical method coded in PhaseTracer is briefly described in App. B. This method is different from the one applied in the code BSMPT [91] and previous works on SFOPTs in the NMSSM [60], which may only find a single PT between the EW symmetric vacuum and the observed EWSB vacuum. Our method is able to map out a more complicated phase structure and find multiple PTs in it. Of equal importance, by analyzing the phase structure obtained by PhaseTracer, we confirmed that not all potential tunnelings actually take place in the early Universe. This may happen because the tunneling rate is too slow or because the PT is located on a branch of the phase structure that the system never utilizes because it evolved in a different direction.
To exhibit spontaneous EWSB as the Universe cooled, the vacuum of the finite temperature effective potential (21) should respect EW symmetry at high temperature, which is 1 TeV in this work, and should violate it at zero temperature. Thus at high temperature the global minimum should lie at the origin, h u = 0 and h d = 0, and at zero temperature the deepest minimum should lie at the observed EWSB VEV. We can use this information to fix the boundaries of the phase structure by finding all minima of the potential at T = 0 and T = 1 TeV and checking that spontaneous symmetry breaking occurs. Starting from T = 0 then we can use PhaseTracer to find all possible PTs.
The strength of such a transition is described by an order parameter. For baryogenesis, we consider the order parameter The singlet VEV is not included here because it does not affect EW sphalerons. Order parameters of about γ EW 1 are considered strong and could catalyze baryogenesis. The Nielsen identities in (20) imply that the critical temperature is gauge independent, since the effective potential is gauge independent at extrema. Our one-loop truncation of the effective potential, however, means that it is gauge independent only at the tree-level extrema. Thus the critical temperature, which we find from the effective potential at the one-loop minima, is gauge dependent. See Ref. [84] for further discussion and a procedure that may enforce gauge independence. The location of the minima, furthermore, and thus the order parameter, always depend on the gauge parameter ξ.
A first order transition occurs through bubble nucleation and there is a finite probability per unit time and volume for tunneling to a new phase. The new phase dominates once the following condition is satisfied [97,98], where S E stands for the Euclidean bubble action, and T N is the so-called nucleation temperature. If there is no solution, we conclude that the transition cannot complete. During the scan, we identify all possible PTs without checking whether they successfully nucleate. After classifying phase structures, we check nucleation temperatures for a subset of our samples using CosmoTransitions [96].

Parameter space, constraints and sampling strategy
To explore all possible PTs in the NMSSM, including strong EWPTs, we sampled the parameter space of the model within the ranges shown in Tab. 1. The first four parameters, λ, κ, A λ and A κ are from the tree-level NMSSM potential and enter the matching conditions at tree-level (11), while the fifth parameter, the stop trilinear A t enters at the one-loop level (12). These parameters are all defined at the matching scale m SUSY which we also take as an input and represents the geometric mean of the left and right soft SUSY breaking masses of the stops, which have been integrated out, i.e.
The final two parameters are the ratio of the Higgs VEVs tan β ≡ v u /v d and the singlet VEV, v S , which are defined at the top quark mass, m t = 173.1 GeV. Therefore our model has eight free parameters. From these inputs the parameters of the THDMS at m t are obtained using Flexible-SUSY-2.1.0 [99,100], coupled with 5 SARAH-4.12.3 [103][104][105][106], which implements the matching and running procedure described in Sec. 2.1, with (11) specified as a boundary condition in the FlexibleSUSY model file. 6 Since all running and effective potential calculations are performed in the THDMS it is not necessary to specify any further soft-breaking masses in the NMSSM. Because the quartic coupling λ can always be made positive through field redefinitions, we do not consider negative values for it, but we do consider both negative and positive values for the soft trilinears, κ and v S . Lastly, as discussed earlier, for self-consistency we only consider real parameters.
The field dependent masses which enter the one-loop corrections to potential are calculated with FlexibleSUSY, and the thermal functions are evaluated using the implementation described in Ref. [107]. We use PhaseTracer to find the phases and critical temperatures by exploring the finite temperature potential between T = 0 and T = 1 TeV, as described in Sec. 4. Since this involves varying the field values that enter the fielddependent masses, in principle it is possible that this could re-introduce large logarithms and lead to perturbativity problems, therefore we do not consider VEVs greater than 1.6 TeV. In practice in all our results the VEVs are significantly smaller than this, and are less than 300 GeV in all but one very special category of points, therefore this restriction does not have an impact on our results. 7 The main experimental constraints on the parameter region of interest come from LEP chargino searches and the observed Higgs properties. The Higgs sector of our model must be compatible with observations of an SM-like Higgs boson with a mass close to 125 GeV. The observed Higgs, however, could correspond to any one of the three neutral Higgs Parameter Range Metric λ 0, π/2 flat |κ| 0, π/2 flat |A λ | 0, 10 TeV hybrid |A κ | 0, 10 TeV hybrid |A t | 0, 10 TeV hybrid m SUSY 1, 10 TeV log |v S | 0, 10 TeV hybrid tan β 1, 60 log Table 1: Ranges and metric of parameters that we scanned in the NMSSM at the SUSY scale. We considered positive and negative κ, v S and trilinear couplings. The "hybrid" metric is flat below 10 GeV, and logarithmic elsewhere. The top mass was fixed to its measured value 173.1 GeV [7]. bosons in our model. We calculated tree-level reduced couplings between the neutral Higgs bosons and SM fermions by taking into account mixing between the neutral Higgs bosons. We furthermore calculated one-loop reduced couplings between the Higgs bosons and photons and gluons using FlexibleSUSY routines developed in Ref. [108]. By passing this information and the Higgs masses to Lilith-1.1.4_DB-17.05 [109], we find a chi-squared, χ 2 Higgs , for our Higgs sector from Run I and II measurements of the Higgs boson at the LHC.
We penalized points in tension with LEP bounds on charginos [7] by introducing a chi-squared for the effective µ-parameter We constructed this function to guide our sampling algorithm towards acceptable solutions with mχ± 1 100 GeV, rather than precisely reflect experimental constraints from LEP. We furthermore penalized points without an SFOPT by the chi-squared The role of this term is to focus our sampling algorithm on SFOPT with γ EW 1; it is in fact equivalent to a Gaussian penalty log 10 γ EW = 0 ± 0.2.
Since the parameter space shown in Tab. 1 is eight-dimensional we sampled points from our model using MultiNest-3.10 [110][111][112] with a chi-squared We saved and considered all points evaluated by MultiNest, i.e., we disabled the cuts ordinarily placed on saved points by the MultiNest algorithm. To be consistent with the LHC Higgs measurements and LEP bounds on charginos [7], and to satisfy our SFOPT requirement, we select points with χ 2 Higgs − min χ 2 Higgs ≤ 6.18, µ eff ≥ 100 GeV and γ EW ≥ 1, where min χ 2 Higgs = 22.3 was the minimum χ 2 Higgs found in our scan. After that, we further required that remaining points satisfied LHC and LEP bounds on BSM Higgs bosons using HiggsBounds-5.3.2beta [113][114][115][116][117], which we interfaced via NMSSMCALC [118].

Classification of phase transitions
After collecting more than three million valid points, we found that the possible phase structures in the NMSSM harbored rich and novel phenomenology. To reflect this phenomenology, we classify these points into three categories that differ by the nature of the first possible PT in the cosmological history: 1. Type-H-and-S: EW symmetry is spontaneously broken such that at least one Higgs field and the singlet field obtain non-vanishing VEVs simultaneously.
2. Type-Only-H: EW symmetry is spontaneously broken by one or both Higgs fields obtaining VEVs, but the singlet VEV remains zero.
3. Type-Only-S: EW symmetry remains unbroken, but the singlet field obtains a VEV.
The Higgs obtain non-vanishing VEVs in a SFOPT afterwards, during which the sign of singlet VEV may be maintained or flipped. Thus we further classify this type into two subcategories: • Type-Only-S(maintain): the strongest PT maintains the sign of singlet VEV.
• Type-Only-S(flip): the strongest PT flips the sign of singlet VEV.
It is important to understand that at this stage we do not have the means to ensure that a PT is definitely part of the cosmological history. More precisely, for such an extensive sample of parameter points, we are not in the position to calculate nucleation temperatures, actions, decay rates, etc. for each potential transition in the phase structure. For this reason, unless specified otherwise when we say 'PT' we typically mean 'possible PT'. To simplify our discussion of this non-trivial structure, we introduce the following shorthand notation: • We denote the minimum value of the potential in a given direction with a calligraphic font. For example, s is a value of singlet field s at a minimum of the scalar potential.
• At a critical temperature, two vacua are degenerate. However, we always define the true vacuum to be the deepest of these vacua just below the critical temperature, and the other one is the false vacuum in our notation.
• In case of multiple SFOPTs we refer to the SFOPT with the greatest γ EW as the strongest one.
as the signed geometric mean of the Higgs fields.

Benchmark points
In Fig. 1, we present a phase history for a typical point in each category. For these benchmark points, we checked our results with CosmoTransitions and calculated the nucleation temperature for every possible FOPT. The corresponding input parameters, Higgs properties and transition information are shown in Tab. 2. On each panel, the lines show the signed geometric mean of the Higgs fields (left) or the singlet field (right) at a minimum of the potential as a function of temperature. 8 Two phases linked by an arrow at a given temperature are degenerate and thus a FOPT could occur in the direction indicated by the arrow (i.e., below the critical temperature, the phase at the end of the arrow contains a deeper minimum). When there is more than one possible sequence of FOPTs that leads from the origin at T = 1 TeV to the observed vacuum at T = 0, we show the FOPTs that belong to the sequence that includes the strongest FOPT by black arrows, and PTs that are not part of that history by gray arrows. Note though that other possible FOPTs between phases that are never degenerate are not marked. For example, in the upper left panel, the minima in phase 2, which appears at about T = 88 GeV, always lies shallower than that in phase 3. A FOPT between them is possible, although there is no critical temperature. From Fig. 1 we can see that in all categories at high temperature, T > 400 GeV, the true vacuum is always at the origin (as described in Sec. 4 this is a requirement in our scan). In the upper left panel, the first (and only) PT occurs at T 145 GeV between (0, 0, 0) and (106,117,276) with γ EW = 1.09 and nucleation temperature T N = 116 GeV. Thus it is classified as Type-H-and-S.
In the upper right panel, only one of the Higgs fields, h u , develops a VEV during the first crossover transition at T = 155 GeV. The first transition in the cosmological history was never first order in our Type-Only-H samples. As the Universe cools, however, a deeper minimum exists between T = 151 GeV and T = 124 GeV at about (0, 0, 450), which belongs to phase 2. The FOPT to this deeper minimum would (temporarily) restore EW symmetry; however, we find that it cannot complete as (34) cannot be satisfied. If it completed, EW symmetry would subsequently be permanently broken by another SFOPT at T 123.6 GeV which would complete, from (0, 0, 463) to (91,162,274) with γ EW = 1.5 and T N = 119 GeV. Indeed, in all the Type-Only-H samples that we found, EW symmetry was broken, possibly restored and finally broken again, and the final FOPT would be the strongest, just as in this example. However, these sequences of transitions are impossible, as the actions for the transitions that restore EW symmetry are always so large that bubbles cannot nucleate properly. Thus although there appear to be SFOPTs with γ EW > 1 and nucleation temperatures in the Type-Only-H samples, they cannot explain the observed baryon asymmetry of the Universe, as a previous transition in the cosmological history would not complete.

Figure 1: Phase structures for typical points in the categories Type-H-and-S (upper left), Type-Only-H (upper right), Type-Only-S(maintain) (lower left) and Type-Only-S (flip) (lower right). The lines show the field values at a particular minimum as a function of temperature. The arrows indicate that at that temperature the two phases linked by the arrows are degenerate and thus that a FOPT could occur in the direction of the arrow. The dots in the lower panels represent transitions that do not change the corresponding field values. The black arrows and dots show a path that includes the strongest EW FOPT, while the gray ones are not in that path.
interesting gravitational wave signatures.
Finally, we consider a Type-Only-S(flip) point (lower right panel). The singlet field develops a negative value during the first transition at T = 368 GeV, which is first-order and completes at T N = 143 GeV. At T 368 GeV, just below the critical temperature of the first transition, a phase with positive s develops, which is approximately symmetric with respect to the phase with negative s. Eventually, a second first-order transition at T = 104 GeV breaks EW symmetry and flips the sign of the singlet by transitioning to this approximately symmetric phase. Although this is the strongest PT, it cannot complete,  as the barrier between the phases means that the tunneling action is too large for (34) to be satisfied. This phenomenon appears in a large fraction of our Type-Only-S(flip) samples. Phase histories of types Type-H-and-S and Type-Only-S(maintain) were previously investigated in Ref. [15,[58][59][60]; however, the richer phase histories in Type-Only-H and Type-Only-S(flip) have not been discussed in the literature as far as we are aware. Note that the barrier between the minima in the Type-Only-H and Type-Only-S(flip) are usually so high that the tunneling may not happen. This shows the importance of studying phase structure as well as calculating the transition strength.

Type-H-and-S
We also checked the robustness of our results against the change of the renormalization scale. For the Type-H-and-S benchmark point in Tab. 2, we found a mild (1% -2%) variation of the critical temperature and the transition strength as the renormalization scale changes in the (m t /2, 2m t ) range. We furthermore checked gauge dependence by repeating our calculations for our benchmark points in the ξ = 0 (Landau) gauge. We found, as anticipated, that gauge dependence was present but typically mild, especially for the critical temperatures. The gauge dependence could, nevertheless, motivate the application of gauge independent techniques in future works.

Reaching the observed SM vacuum
During the scan we required that the deepest minimum at zero temperature agreed with the observed VEV, h = 246 GeV. We call the phase associated with the observed VEV the SM vacuum. We split our samples by two ways of reaching the SM vacuum. First, in Sec. 5.4.1 we consider samples for which the strongest SFOPT ends in the SM vacuum, which changes smoothly to h = 246 GeV at T = 0. Second, in Sec. 5.4.2 we consider samples for which the strongest FOPT does not end in the SM vacuum. As we discuss, such samples must feature at least one further FOPT that ultimately ends in the observed vacuum at T = 0. In both cases, the Type-Only-H scenario was by far the rarest, with noticeably few samples shown in the following scatter plots.

The strongest FOPT ends in the SM vacuum
We selected samples in which the strongest FOPT ended in the SM vacuum. For our samples, it was sufficient to check that h u > 0 GeV and h d > 0 GeV for the true vacuum of strongest FOPT. All of our benchmark points in Tab. 2 are in this category. In Fig. 2, we present the true and false minima of the strongest FOPT at the critical temperature. It demonstrates some features of each of the types of point that we described above. For Type-H-and-S, the first transition, in which the Higgs and singlet fields acquire VEVs, is usually also the strongest FOPT. There are however three exceptional points where the singlet field values at the false minimum are non-zero. They have similar phase structures to the upper right panel of Fig. 1 except that the minima of phase 1 is always deeper than the minima of phase 2 in all three cases. Thus there is no critical temperature between these phases, and so the strongest FOPT for these three points is not in the cosmological history.
According to the definition of Type-Only-H, only the Higgs fields obtain VEVs in the first transition in the history, while the upper right panel of Fig. 2 shows that the false vacuums of strongest FOPT have zero Higgs VEVs, h u = h d = 0, but a non-vanishing singlet VEV, s = 0. This means that there must be an intermediate transition that restores EW symmetry and generates a singlet VEV. Since the number of Type-Only-H scenarios that we found are quite small, we checked each one in detail. We found that this intermediate transition exists for all Type-Only-H samples, but the corresponding tunneling probabilities are too small. Nonetheless it is possible that there exist scenarios of this type where the transition does complete. The lower panels of Fig. 2 display samples of Type-Only-S where the strongest FOPT maintains (left) or flips (right) the sign of the singlet VEV. We see that the singlet VEV can evolve to up to 1.6 TeV after the first transition, and then shifts to about 150 GeV to 650 GeV during the strongest FOPT. The singlet VEV s of the true vacuum can be both positive or negative, because the input v S includes both signs. We have checked that the singlet VEV at the true vacuum has the same sign as the input v S .
In all scenarios, the spread in the possible true vacuum for the Higgs fields at the critical temperature is small, and typically it matches and rarely exceeds the input EWSB vacuum, i.e., h 246 GeV. This can be further seen in Fig. 3 We now delineate the regions of the NMSSM parameter space in which our four scenarios occur. We checked that in all scenarios the stops were truly decoupled by checking stop mixing, X t = A t − µ eff cot β, which could potentially split the stop mass eigenvalues making one of them light. We found that most samples were actually concentrated within the range −m SUSY ≤ X t ≤ m SUSY and no particular value of m SUSY was preferred by our samples.
In Fig. 4 we show that the Higgs sector parameters (µ eff , tan β) are severely constrained. Indeed, the Type-H-and-S and Type-Only-H scenarios require tan β 3, whereas the Type-Only-S(maintain) and Type-Only-S(flip) scenarios permit tan β 7 and tan β 17, respectively. For all types, the upper limit of tan β decreases with µ eff increasing. The effective µ-parameter, and thus the higgsinos, are always light, |µ eff | 300 GeV. Thus we find further motivation for scenarios with small µ eff 1 TeV, which are also motivated by naturalness, and we anticipate that the searches for higgsinos at the LHC could be sensitive to our models. Samples with µ eff < 0 were extremely rare in the Type-H-and-S and Type-Only-S(maintain) scenarios, and not present in the Type-Only-H samples.
We see, furthermore, in Fig. 5, that quartic couplings of around λ ≈ 0.6 and κ ≈ 0.2 could result in an SFOPT in all our scenarios, though a broad range of couplings result in SFOPTs in Type-Only-S(flip) scenario, including couplings with values far above the limits that would be set if we required perturbativity up to the GUT scale. The constraints strongly prefer that λκ > 0, a combination that is invariant under the field redefinition S → −S. Since we worked in a λ > 0 convention, the inequality λκ > 0 is equivalent to κ > 0. In the Type-Only-S(maintain) scenarios, however, we find a few solutions for which κ < 0. Fig. 6 shows the trilinear couplings (A λ , A κ ) with the quartic coupling κ shown by the color bar. The trilinears play an important role. As different types of sample require different sign of singlet VEV at low temperatures, the parameter space of each type shows distinguishable tendency. The samples in Type-H-and-S, Type-Only-H and Type-Only-S(maintain) scenarios are concentrated at negative A κ with positive κ or positive A κ with negative κ, as well as a horizontal slice of points at A κ ≈ 10 GeV for Type-H-and-S and Type-Only-S(maintain). On the other hand, A λ is typically positive but 500 GeV. The one point with negative A λ in Type-H-and-S and the two points with negative A λ in Type-Only-S(maintain) correspond the point of negative µ eff in Fig. 4. The distinction between Type-H-and-S and Type-Only-S(maintain) is that Type-Only-S(maintain) favors smaller A κ and A λ . Finally, Type-Only-S(flip) shows two approximately symmetric regions that were previously identified in Fig. 4 by the sign of µ eff . The region of positive (negative) A λ and A κ corresponds to positive (negative) µ eff .
We emphasize again that the parameter spaces shown in Fig. 4, Fig. 5 and Fig. 6 can only ensure the existence of a SFOPT with γ EW 1. Establishing whether this SFOPT is definitely part of the cosmological history requires further investigation, which we only present for our benchmark points.

The strongest FOPT does not end in the SM vacuum
Other than the scenario discussed above, we have plenty of samples in which the strongest FOPT does not end in the SM vacuum, as shown in Fig. 7. In these samples, in the true vacuum for the strongest FOPT, h is always negative and s is either zero or has a different sign to µ eff , so this almost certainly does not belong to the SM vacuum in which h = 246 GeV. The spread in the possible true vacuum for the Higgs fields at the critical temperature is substantial, and could differ considerably from the observed EW vacuum. Because of this, we no longer find that h ≈ 246 GeV, allowing enhancement or suppression of the strength of the PT in Fig. 8, which differs markedly from Fig. 3. Indeed, in the Type-Only-S(maintain) scenario, SFOPTs are possible for substantial critical temperatures of up to T C 500 GeV.
At first glance, these points might seem uninteresting, as they do not end in the correct zero-temperature vacuum. They may be especially interesting, however, as this means that in order for such samples to achieve the correct zero-temperature vacuum, there must be another EW FOPT transition or sequence of transitions that complete and end in the correct vacuum. Thus in Fig. 9 we histogram the number of possible FOPTs with γ EW 1 for each sample. Let us stress that strictly speaking, we count the number of temperatures at which two vacua are degenerate. This differs from the number of FOPTs that can take place in one cosmological history, since only particular routes through the phases are possible. For example the upper right panel of Fig. 1 exhibits one or two FOPTs from phase 1 to phase 3, but we would count this though as three. Furthermore FOPTs may also occur between phases that were never degenerate, but such possibilities are not included in our count.
For the samples that end in the SM vacuum (left panel), there is usually a single FOPT with γ EW > 1, except in the Type-Only-H scenario, in which there are often two FOPTs with γ EW > 1. For the samples that do not end in the SM vacuum (right panel), almost all of Type-Only-H and Type-Only-S(maintain) samples and about half the Type-H-and-S samples have more than one EW SFOPTs. We also checked that for most of them the second strongest FOPT does end in the SM vacuum.
Thus, without further calculations, the samples for which the strongest FOPT does not end in the SM vacuum could still potentially explain the observed baryon asymmetry. We display the parameter spaces in Fig. 10, Fig. 11 and Fig. 12. Compared to the scenario in which the strongest FOPT ends in the SM vacuum, the parameter spaces of Type-H-and-S and Type-Only-H are roughly unchanged, while Type-Only-S(maintain) and Type-Only-S(flip) exchange parameter spaces with each other. This is because here the Type-Only-S(maintain) (Type-Only-S(flip)) requires a minimum on the singlet axis with s < 0 (s > 0), opposite to the Type-Only-S(maintain) (Type-Only-S(flip)) scenarios in which the strongest FOPT ends in the SM vacuum.
From Fig. 10 we see that the constraints on the effective µ-parameter are stricter than they are in the scenario in which the strongest FOPT ends in the SM vacuum, especially for small tan β. The Type-H-and-S, Type-Only-H and Type-Only-S(flip) scenarios require an effective µ-parameter smaller than about 200 GeV, whereas the Type-Only-S (maintain) permits µ eff 400 GeV. The slender bar in the Type-Only-S(maintain) scenario at tan β 1 and µ eff ∈ [200, 400] GeV corresponds to samples with T C 200 GeV for the strongest FOPT, displayed in the lower left panel of Fig. 8.
In Fig. 11, a visible difference appears in Type-Only-S(maintain) compared to Fig. 5. The parameter space of λ and κ splits into two separate regions, and relatively large λ 0.5 is favored. For instance, when κ 1.4, here λ is always larger than 1, whereas in lower left panel of Fig. 5 λ can be as low as 0.5.
On the trilinear couplings (A λ , A κ ) plane, there are two additional regions in the Type-Only-S(maintain) scenario (lower left, Fig. 12) compared to the Type-Only-S(flip) samples for which the strongest FOPT end in the SM vacuum (lower right, Fig. 6). First, there is an additional region at A λ 0 GeV. This region corresponds to the previously mentioned region at tan β 1 and µ eff ∈ [200, 400] GeV, with T C 200 GeV for the strongest FOPT. Second, there is an additional region at A κ −50 GeV and A λ > 0. This region is similar to one in the Type-Only-S(flip) scenario (lower right, Fig. 12). Indeed, for this region, as well as the strongest FOPT that maintains the sign of singlet, there is another weaker FOPT that flips the sign of singlet.
In summary, the scenario in which the strongest FOPT does not end in the SM vacuum introduces new interesting regions of parameter space that were not covered by the scenarios in which the strongest FOPT ends in the SM vacuum. These scenarios may be especially interesting because they could be followed by additional FOPTs. However, at the same time there is an additional requirements to ensure that the subsequent transitions actually lead to the EW breaking phase we observe today, which we have not checked.

Properties of the Higgs bosons
As shown by our benchmark points, although our points pass experimental constraints from LEP and the LHC, our scenarios are not in a decoupling regime in which Higgses other than the 125 GeV one are heavy. This is not surprising since it is well known that in the NMSSM a light singlet Higgs state plays an important role in generating a FOPT that breaks EW symmetry [16,25,55], without the need for light stops which are heavily constrained by LHC searches [119,120]. In fact, in all our benchmarks, all Higgs bosons are lighter than about 400 GeV, while there are always at least two CP even Higgs states with masses below 600 GeV in the samples from our scan, with the SM-like Higgs being either h 1 or h 2 .
In Fig. 13 we show the masses of the non-SM-like CP even neutral Higgs bosons in our four scenarios by plotting the mass of h 3 , which is never SM-like, against the mass of the Higgs (either h 1 or h 2 ) that did not play the role of the SM-like Higgs. Samples that are allowed by experimental constraints are shown by green points. We also show excluded samples to aid explanations (gray and blue points).
For the samples where the strongest FOPT ends in the SM vacuum we see that the SM-like Higgs is actually the next to lightest CP even Higgs for almost all allowed samples (green points) in Type-H-and-S, Type-Only-H and Type-Only-S(maintain), with just three exceptions that all appear in the Type-H-and-S samples. In contrast, in the Type-Only-S(flip) scenarios, the SM-like Higgs can be either the lightest Higgs or the next to lightest Higgs. The samples where the strongest FOPT does not end in the SM vacuum show very similar results, but as usual the patterns of the Type-Only-S(maintain) and Type-Only-S(flip) scenarios are exchanged. The reason we see so few samples where the SM-like Higgs is the lightest state for the categories mentioned above, seems to the constraints from the observed Higgs signals, which we have implemented with Lilith. We note that, although it is not clear in the plot, for these types of scenarios there are already a significantly larger number of gray points where the SM-like Higgs is the second lightest CP even Higgs boson than there are for the case where it is the lightest. Applying the constraints for the SM-like Higgs from Lilith-1.1.4_DB-17.05 then reduces the samples where it is the lightest to almost zero. The samples excluded by HiggsBounds-5.3.2beta in Type-Only-S(flip) scenario for which the strongest FOPT ends in the SM vacuum and Type-Only-S(maintain) scenario for which the strongest FOPT does not end in the SM vacuum shown by blue points, are mainly restricted by the constraints from an LHC search for a scalar resonance decaying to a pair of Z bosons [121]. It is also worth noting that even without the requirement of an SFOPT, similar observations have been made in the NMSSM previously. A preference for the SM-like Higgs being the next to lightest one was also found in a global analysis of the NMSSM [122] that did not consider PTs.

Conclusions
Motivated by EW baryogenesis and gravitational wave experiments, in this article we investigated the properties of PTs in the NMSSM. We employed an effective field theory approach to calculate the finite temperature effective potential by matching the NMSSM to the THDMS. By tracing the change in the minima of the effective potential with temperature, we mapped out the phase structure and computed the strengths of any EWPTs, γ EW . By scanning the parameter space of the NMSSM, we obtained millions of samples that featured an SFOPT with γ EW > 1 and satisfied the constraints from LHC Higgs measurements and LEP bounds on charginos.
We classified them into three categories, Type-H-and-S, Type-Only-H and Type-Only-S, based on the nature of the first PT in their cosmological histories. The Type-Only-S samples were further divided into Type-Only-S(maintain) and Type-Only-S (flip) according to whether the singlet VEV changed sign during the strongest EWPT. In the Type-H-and-S samples, the first PT in the cosmological history breaks EW symmetry and gives the singlet a VEV at the same time. This transition is usually the strongest one. The strongest FOPT does not end in the SM vacuum The Type-Only-H samples, on the other hand, go through a series of PTs that break, restore and break again EW symmetry. The first one is a crossover transition during which only the h d field obtains a non-vanishing VEV, and the last one is the strongest EW FOPT. This scenario was by far the rarest in our scan. For the Type-Only-S(maintain) samples, during the first transition EW symmetry remains unbroken, but the singlet field obtains a non-vanishing VEV. Then EW symmetry breaks through a subsequent FOPT. Both of the transitions can be SFOPTs, which could give interesting gravitational wave signatures [123] as well as triggering an EW baryogenesis mechanism. The first PT of the Type-Only-S(flip) samples is usually a FOPT with very small γ, and the following strongest FOPT flips the sign of the singlet VEV. We found, however, that the tunneling rates in Type-Only-H and Type-Only-S(flip) scenarios could be problematic. For our benchmarks, the SFOPT in the Type-Only-H scenario did not complete, and in the Type-Only-S(flip) scenario, a preceding PT required to reach the SFOPT did not complete. Thus, unfortunately, these scenarios might not help EW baryogenesis.
The regions of NMSSM parameter space in which the four scenarios occur show different features. In the samples for which the strongest FOPT ends in the observed zero temperature phase: • The observed 125 GeV Higgs is often the second lightest Higgs in the model, not the lightest one.
• All of the input parameters are severely constrained, except the SUSY scale m SUSY and stop trilinear A t .
• Quartic couplings of around λ 0.6 and κ 0.2 could result in an SFOPT in all our scenarios, though a broad range of couplings result in SFOPTs in the Type-Only-S (flip) scenario, including couplings far away from limits on perturbativity.
• The scenarios predict different trilinear couplings, i.e., they are distinguishable on the (A λ , A κ ) plane. The A λ and A κ of the Type-Only-S(flip) samples always have the same sign, while in the other scenarios the samples are concentrated in the quadrant of negative A κ and positive A λ . Compared to the Type-H-and-S scenario, the Type-Only-S(maintain) scenario favors smaller |A κ | and A λ .
In addition we found substantial samples for which the strongest FOPT does not end in the SM vacuum. The regions of parameter space are similar to the samples for which the strongest FOPT ends in the SM vacuum, except that Type-Only-S(maintain) and Type-Only-S(flip) exchange parameter spaces with each other. There are, furthermore, two additional regions that appear in the Type-Only-S(maintain) scenario, and one of them results in critical temperatures higher than 200 GeV.
In summary, we mapped out and classified intricate patterns of symmetry breaking that are possible in the NMSSM, and checked which scenarios could in principle help provide a viable theory of EW baryogenesis or potentially lead to a gravitational wave signal. We found viable scenarios in which the Higgs fields and singlet or only singlet first acquired VEVs. We checked that the sequences of required PTs actually nucleated, contained a SFOPT, and that the model satisfied constraints from LEP and the LHC. The combination of constraints lead to the predictions that λ 0.6, κ 0.2 and that the observed 125 GeV Higgs tends to be the second lightest Higgs in the model.

A Field Dependent Masses
When exploring the potential away from minima we need to account for the Higgs field dependence of the MS mass eigenstates in (13). Therefore here we present the so-called field dependent masses of the THDMS.
The field dependent masses of the gauge bosons and top, bottom and tau fermions are given by the simple tree-level expressions where the gauge couplings are without GUT normalization, and the y t , y b and y τ are the bosons from physical Higgs bosons. In the ξ = 1 gauge, however, they are treated on an equal footing and we do not need to identify Goldstones.
Finally, the charged Higgs mass matrix is We apply an algorithm developed in CosmoTransitions [96] to trace phases in steps no greater than ∆T : 0. We select a minima m ≡ (h u , h d , s) at temperature T to trace.
1. We use a local minimum finding algorithm, such as Nelder-Mead [124], to find the minimum m at T = T + ∆T .
2. We check that the new minimum m lies close to that expected from a shift caused by thermal corrections.
We calculate the difference 3. If R ≤ max R, where max R governs the maximum acceptable changes in the field, we accept that the minima m at temperature T belongs to the same phase as the minima m at temperature T . We continue to trace the phase by returning to step 1 with m → m and T → T , and we reset any changes to ∆T .
If R > max R, we assume that the change in temperature dramatically changed the potential. We reduce the change in temperature by a factor of two, ∆T → ∆T /2, and return to 1.
If, however, R > max R and |∆T | < min ∆T , where min ∆T governs the smallest permissible step in temperature, we conclude that the phase must have ended, as the minima appears to change abruptly with a small change in temperature.
We save the sequence of minima and temperature found through this process -this is a phase. We find all the phases by tracing all T = 0 minima up to at most 1 TeV (the phase may end earlier) and all T = 1 TeV minima down to T = 0 (in which case ∆T < 0). After removing degenerate phases, we denote the i-th phase by m i (T ).
If any two of the phases, e.g., the i-th and j-th phase, coexist between temperatures T 1 and T 2 , and if there must exit a critical temperature, T C , between temperatures T 1 and T 2 at which they are degenerate, We calculate the critical temperature using bisection, and find properties of the transition, e.g., the strength of transition from (33).