Cosmological Phase Transitions and their Properties in the NMSSM

We study cosmological phase transitions in the Next-to-Minimal Supersymmetric Standard Model (NMSSM) in light of the Higgs discovery. We use an effective field theory approach to calculate the finite temperature effective potential, focusing on regions with significant tree-level contributions to the Higgs mass, a viable neutralino dark matter candidate, 1-2 TeV stops, and with the remaining particle spectrum compatible with current LHC searches and results. The phase transition structure in viable regions of parameter space exhibits a rich phenomenology, potentially giving rise to one- or two-step first-order phase transitions in the singlet and/or $SU(2)$ directions. We compute several parameters pertaining to the bubble wall profile, including the bubble wall width and $\Delta\beta$ (the variation of the ratio in Higgs vacuum expectation values across the wall). These quantities can vary significantly across small regions of parameter space and can be promising for successful electroweak baryogenesis. We estimate the wall velocity microphysically, taking into account the various sources of friction acting on the expanding bubble wall. Ultra-relativistic solutions to the bubble wall equations of motion typically exist when the electroweak phase transition features substantial supercooling. For somewhat weaker transitions, the bubble wall instead tends to be sub-luminal and, in fact, likely sub-sonic, suggesting that successful electroweak baryogenesis may indeed occur in regions of the NMSSM compatible with the Higgs discovery.


I. INTRODUCTION
The origin of the matter-antimatter asymmetry in the Universe remains one of the key open problems at the interface of cosmology and particle physics. While a primordial asymmetry is not observationally ruled out, the relatively small size of the observed asymmetry -roughly one part in 10 billion -combined with the success of the inflationary paradigm in diluting away primordial relics, suggest a dynamic origin of the asymmetry. Many such baryogenesis mechanisms have been proposed, occurring at widely differing points in the history of the Universe, and involving correspondingly different physical scales [1].
Among baryogenesis scenarios, those invoking physics at the electroweak phase transition (electroweak baryogenesis, or EWB) have attracted much interest. This is because such models generically require new physics at scales that are within the reaches of current-generation colliders, intensity-frontier precision measurements, and, potentially, planned gravitational wave detectors at the cosmic frontier (for a recent review of electroweak baryogenesis see e.g. [2]). One framework which, in principal, contains all of the ingredients necessary for successful electroweak baryogenesis is the minimal supersymmetric extension of the Standard Model, or MSSM [3].
The discovery of a Higgs particle at the Large Hadron Collider (LHC), and the lack of any signals for supersymmetric particles, have started to put significant pressure on the corner of MSSM parameter space compatible with EWB [4]. In particular, it was shown in a recent analysis [5] that requiring both a strongly first-order phase transition and a Higgs mass around 125 GeV forces one of the two stop squark masses to be in the 10 6 TeV range, with the other one around 120 GeV. This in turn requires a light neutralino, with a mass lower than 60 GeV, contributing to the invisible decay width of the Higgs and thereby reducing the enhancement in gluon-gluon fusion Higgs production from the light stop [6].
Additionally, searches for stop squarks at the LHC have all but ruled out the light stop scenario in the MSSM directly [7,8], albeit with some possible loopholes. Apart from the requirement of a strongly first order phase transition, the CP-violating sources in the MSSM are also highly constrained by current limits from electric dipole moment experiments [9][10][11].
In light of these recent developments, MSSM EWB appears quite unnatural at best, if not plainly contrived.
An alternative supersymmetric framework for electroweak baryogenesis is one in which the MSSM Higgs sector (consisting of two Higgs doublets) is enlarged to include an additional singlet superfield -a setup known as the next-to-minimal supersymmetric extension to the Standard Model, or NMSSM. It has long been appreciated [12,13] that tree-level cubic terms in the Higgs potential might provide a strongly first-order electroweak phase transition, relaxing the model-building challenges one faces in the MSSM. Nevertheless, relatively little is known about the details of how the phase transition proceeds in this case. Many previous studies of the NMSSM and related non-minimal SUSY models have exclusively considered one-step electroweak phase transitions [12][13][14][15][16][17][18][19][20][21][22], while others have considered the possibility of a more complicated pattern of symmetry breaking [23], albeit prior to the Higgs discovery and with several assumptions about the spectrum that are no longer justified. Also lacking are detailed predictions for the quantities pertaining to the bubble wall profile and which enter the calculation of the baryon asymmetry in the context of EWB. To the best of our knowledge, amongst the previously mentioned studies only Ref. [14] attempts this task.
However, they considered a slightly different scenario, and, more critically, did not have knowledge of the Standard Model Higgs discovery. A major aim of our study is to bridge this gap.
The baryon asymmetry produced during electroweak baryogenesis depends crucially on the SU (2) and singlet wall thickness between the broken and unbroken electroweak phase (L w and L s ), the relative variation of the neutral SU (2) Higgs fields across the wall (a quantity dubbed ∆β), and the velocity of the wall relative to the plasma (v w ). In many cases the resulting baryon asymmetry is directly proportional to ∆β and inversely proportional to L w [24,25]. The dependence on the wall velocity is more complicated [24,25]: if v w is too small, the baryon asymmetry will be diluted by interactions with the plasma in front of the bubble wall, while for quickly moving walls, the SU (2) sphalerons do not have sufficient time to act on the diffusing chiral current, leading again to a suppression of the asymmetry.
In the literature, typical ranges for L w and ∆β have been given for the case of the MSSM [26] and the general (non-scale-invariant) NMSSM [14], but no detailed studies exist taking The present study also addresses a few technical aspects of the calculation of the electroweak phase transition that had previously been glossed over in the NMSSM. Most importantly, to deal with potentially large logarithms arising in the one-loop effective potential, we integrate out both stop squarks (which we assume are heavy) and work in an effective field theory with two Higgs doublets, a singlet, neutralinos/charginos, and the remaining SM spectrum. We have implemented this setup in the widely-used CosmoTransitions software package [27] and detail our set-up for use in future studies.
Our approach will be to focus on a well-defined corner of the NMSSM parameter space which is both phenomenologically viable and demonstrative of a broad variety of outcomes for the patterns of electroweak symmetry breaking. By "viable" we mean that every one of the points we consider features a Higgs sector entirely compatible with results from the LHC, a sparticle spectrum compatible with LHC searches, a lightest neutralino with a thermal relic abundance matching the observed dark matter density, and with direct and indirect detection rates in accordance with current limits.
Focusing on a specific region of the NMSSM parameter space, we produced a set of points exhibiting several different patterns of electroweak symmetry breaking. For each point, we calculated the various phase transition properties, including the relevant friction coefficients needed for a detailed estimate of the bubble wall velocity, as well as the predicted spectrum of gravitational waves from bubble collisions in the early universe. One of our key findings is that even this small region of parameter space can exhibit a broad variety of outcomes for the quantities relevant for EWB. Parameters pertaining to the bubble wall profile and expansion, which are critical inputs for computing the baryon asymmetry in EWB, can vary by up to an order of magnitude, even across the small slice of parameter space near our benchmarks, and take on favorable values for successful EWB. For example, the SU (2) wall width varies in the broad range 5 L w T 20; the parameter ∆β also spans a large range of values, 0.01 ∆β 0.5; the bubble wall velocity is found to be roughly ∼ O(0.01 − 0.1) for sub-sonic bubbles. We do not find any detectable level of gravitational wave emission for any of the points considered.
The rest of this manuscript is organized as follows. Sec. II introduces the model under consideration (the scale-invariant NMSSM), the associated effective potential, and our computational strategy. In particular, we explain how we deal with large logarithms, and detail how we implement the effective potential in CosmoTransitions. Sec. III presents the parameter space we consider, focusing on the phenomenology of the relevant particle sectors (neutralinos/charginos, sfermions, Higgses). Subsequently, Sec. IV describes how we study the electroweak phase transition and its associated properties, and presents our key findings. The computation of the wall velocity is then given in Sec. V

II. THE SCALE-INVARIANT NMSSM AND EFFECTIVE FIELD THEORY
A. Overview of the Strategy Before delving into the details of our analysis, we outline the key steps comprising our strategy for studying phase transitions in the NMSSM: 1. Begin by choosing NMSSM parameters at M t , the scale set by the stop squark masses.
These can be fed into a spectrum calculator (we use NMSSMTools [28], as discussed further below) for a precise determination of the spectrum to check against phenomenological constraints.
2. Integrate out the stop squarks from the spectrum, since they are assumed to be heavy. 4. Solve the renormalization group equations (RGEs) for the various quartic and dimensionful parameters in the potential to evolve them from M t to m t . We will do so approximately using a fixed order approximation, but this can also be done exactly.
5. Include the 1-loop zero temperature contributions to the potential. Impose counterterms to minimize the dependence of the effective potential on the renormalization scale Λ.
6. Add the finite temperature contribution to the potential. The full 1-loop effective potential can then be fed into CosmoTransitions [27] to compute the phase transition properties, as discussed in more detail in Sec. IV A.
Readers not concerned with the details of the above procedure should proceed directly to our discussion of the particle spectra in Sec. III. The remainder of this Section is devoted to elucidating the aforementioned steps of our analysis.

B. The Model
Our starting point is the so-called scale-invariant NMSSM. In this incarnation, a Z 3 symmetry forbids dimensionful parameters in the superpotential W , which is here given by The hatted quantities above represent the chiral SU (2) superfields S is a gauge singlet chiral superfield. The dot represents the usual antisymmetric SU (2) product. Supersymmetry is softly broken, with The tree-level potential is then given by where g 2 ≡ (g 2 1 + g 2 2 )/2, g 1 and g 2 denote the U (1) and SU (2) gauge couplings, respectively, and the sum over the F -terms is over H 0 u,d , S, with F ≡ ∂W/∂φ i . For details regarding the rest of the spectrum, we refer the Reader to Refs. [21,28], whose notation and conventions we follow here unless otherwise stated (see also Ref. [29] for a more detailed review).

C. NMSSM Effective Potential
Anticipating electroweak symmetry breaking, we write the Higgs and singlet fields as In terms of these fields, the portion of the tree-level potential relevant for the electroweak phase transition becomes The vacuum expectation values (VEVs) of the scalar fields are assumed to be real at all temperatures, and we do not consider charged vacua (although we do ensure that the potential is stable in the charged and imaginary directions). After electroweak symmetry breaking, and at zero temperature in the physical vacuum, we define h u ≡ v u , h d ≡ v d , and s ≡ v s as usual. Note that the non-zero vacuum expectation value of s generates an effective µ term in the superpotential Eq. 1, µ ≡ λv s / √ 2.
The zero temperature potential in Eq. 5 receives quantum corrections from all fields which couple to h u,d and s. These corrections are given by the well-known Coleman-Weinberg expression Here m 2 i are the field-dependent mass-squared values, n i are their associated number of degrees of freedom, and Λ is the renormalization scale. The constants c i depend on the renormalization scheme. We choose to work in the M S prescription, whereby c = 1 2 for the transverse polarizations of gauge bosons, while c = 3 2 for their longitudinal polarizations and for all other particles. The plus and minus signs in Eq. 6 are for bosons and fermions, respectively. The sum over the relevant particles i includes all Standard Model particles (we ignore fermions lighter than the bottom quark since their Yukawa couplings are small), the physical Higgs and other scalar particles, their associated Goldstone bosons, the neutralinos and the charginos. We work in Landau gauge so that the ghost bosons decouple and need not be included in the spectrum.
Note that the one-loop potential contains explicit gauge-dependence [30][31][32][33], which cancels against the implicit gauge-dependence of the VEVs at every order in . As is common practice, we do not consider the effects of the implicit gauge-dependence, and so our results will contain gauge artifacts. However, our primary purpose in examining the effective potential is to estimate whether or not a first-order phase transition is possible and to infer its general properties in comparison with e.g. the MSSM, for which purpose a calculation with gauge-dependence is acceptable. Additionally, since the transitions are driven primarily by the singlet contributions to the potential, Landau gauge would appear to be a reasonable choice, since the gauge-dependent contribution to the finite temperature effective potential in this gauge is small. Defining a gauge-independent version of the one-loop effective potential including the full bosonic thermal contribution is an open problem and is beyond the scope of the present study.
We will be primarily interested in regions of the NMSSM with moderately heavy (∼ 1 − 2 TeV) stops. The stops enter the 1-loop zero-temperature effective potential through Eq. 6.
The top Yukawa coupling controls both the coupling of the top quarks and stop squarks to the Higgs fields, and since there is a large hierarchy in scale between the tops and stops, large logarithms will inevitably arise from a naive application of Eq. 6. This situation is quite unacceptable for calculating the phase transition properties using the one-loop potential, since the large differences between tree-level and one-loop expressions signify a substantial dependence on the renormalization scale Λ and the need to go to higher orders in the loop expansion. To avoid this issue, we employ an effective field theory approach, which we describe below. It involves integrating out the stops and working in the resulting model with two Higgs doublets, a singlet, neutralinos/charginos, and the rest of the SM spectrum with a renormalization scale near the electroweak scale.

D. Dealing with Large Logarithms
In keeping with the usual convention, we will use as inputs the NMSSM parameters defined at the stop mass scale, M t . To avoid the large logarithms arising from the stops in Eq. 6, we will integrate the stops out at M t (see e.g. Refs. [34][35][36][37][38]). Below this scale, the tree-level potential can be described by the following general two Higgs doublet + singlet (2HD+S) potential [39][40][41] Comparing the above potential with Eq. 5, at the scale Λ = M t one can define Similarly, for the mass terms, we have Integrating out the stops results in threshold corrections to the parameters in Eq. 8.
Keeping only renormalizable terms in the potential, the only relevant threshold correction at M t is to the h u quartic, Then, at the scale M t , the parameters in Eq. 7 are given by and where ∆λ i = 0 for i = 2.
Ultimately, in analyzing the phase transition structure of the theory using the effective potential, we would like to work at a renormalization scale Λ ∼ m t to reduce the logarithmic contribution from the top quark in the physical minimum. To do so requires the various parameters at the scale Λ, which can be obtained by solving the relevant renormalization group equations. The most important contributions to the RGEs are those from the top quarks, gauge bosons, Higgs and singlet bosons, Higgsinos, and singlinos. We do not include the gaugino contributions, since in the benchmark points we study the wino is always rather heavy and the bino does not significantly affect the running. We list the relevant RGEs, along with more details about which contributions we consider, in Appendix A.
A complete resummation of the large logarithms requires solving the full set RGEs.
However, for our purposes it is sufficient to consider the lowest order solutions, given by The above approximation corresponds to keeping only the first term in the loop expansion.
With the parameters defined at m t , one can take derivatives of the effective potential of Eq. 7 and obtain the mass matrices for the various Higgs bosons. We will denote the lightest CP-even (odd) Higgs mass eigenstate as h s (a s ), since it is very singlet-like for all of our benchmarks. The second-lightest CP-even Higgs mass eigenstate is denoted by h and will be Standard Model-like for all points considered. We use v = 246 GeV, as well as the values of tan β and v s (obtained from µ and λ) as input parameters at Λ = m t instead of m 1,2,3 , solving for the latter by minimizing the tree-level potential. We find that the above procedure yields good agreement with more complete calculations, the masses and mixing matrix entries for the scalars and pseudoscalars falling within a few percent of those calculated using NMSSMTools [28]. For our purposes, this is sufficient and it encapsulates the sizable corrections from the stop sector to the SM-like Higgs mass, as well as the other dominant 1-loop contributions.
The one-loop effective potential given by Eq. 6 has non-zero derivatives at the tree-level minimum, so the one-loop and tree-level minima do not coincide. We follow the strategy of Ref. [42] and employ counterterms to cancel the one-loop effects, ensuring that the electroweak minimum of the one-loop plus counterterm-corrected effective potential is the same as that given by the parameters β, µ, and λ. The counterterms can be written as In principle, one could also impose counterterms to ensure that the masses do not change between the tree-level and one-loop potentials. However, we find that the masses at oneloop do not differ drastically from their tree-level counterparts, and since it is the tree-level masses that enter the one-loop finite temperature contribution to the effective potential, we do not include these terms.
Finally, using the improved effective potential, we can compute the various masses for the particles and use them to calculate the 1-loop zero-temperature effective potential via Eq. 6. The full one-loop contribution at finite-temperature is given by where and again the upper (lower) signs correspond to bosons (fermions). At high temperature, the validity of the perturbative expansion of the effective potential breaks down. Quadratically divergent contributions from non-zero Matsubara modes must be re-summed through inclusion of thermal masses in the one-loop propagators [43,44]. This amounts to adding thermal masses to the longitudinal gauge boson degrees of freedom and to all of the scalars.
That is, the bosonic mass matrices M 2 ij , from which the individual eigenvalues m 2 i are calculated, receive extra thermal mass contributions M 2 ij → M 2 ij +Π ij T 2 . The thermally corrected eigenvalues are then re-input into Eq. 14, yielding the re-summed finite-temperature effective potential.
The full one-loop effective potential at finite temperature is then given by where the masses m 2 i are field-dependent (calculated from the 2HD+S potential with counterterms) and include thermal mass corrections. This potential can then be used to determine the phase structure of phenomenologically viable NMSSM parameter space points.
We emphasize that the application of this approach to the NMSSM is a novel feature of our analysis and allows for a correct treatment of the stop sector in studying phase transitions close to the electroweak scale. For a similar strategy in the MSSM at two loops, see Ref. [37].
With our strategy laid out, we now turn to the parameter space we wish to explore.

III. THE PARAMETER SPACE
Even in its Z 3 symmetric incarnation, the NMSSM parameter space is very large. Instead of performing numerical scans over all of the parameter space, our approach is to study the phase transition properties in a particular region of the NMSSM motivated by Higgs physics, SUSY searches, dark matter constraints, and naturalness arguments. As we will see, even in the small parametric window we consider, the phase structure exhibits a rich phenomenology.
Most recent studies [20,22] have focused on regions with very light Higgs and neutralino states. In light of the Higgs discovery and the non-observation of non-SM particles at the LHC, we will instead consider parameter space which accommodates: • A significant tree-level contribution to the SM-like Higgs mass. This corresponds to regions with sizable λ and low tan β.
• Very little mixing between the singlet-like and SM-like Higgses, the latter in good agreement with LHC observations.
• Moderately heavy stops and other sfermions, in the 1-2 TeV range. This is motivated by naturalness arguments [45][46][47] and the non-observation of superpartners at the LHC.
• A viable neutralino dark matter candidate which saturates the observed relic abundance and is compatible with direct-and indirect-detection experiments.
• Higgs and chargino/neutralino spectra compatible with current LHC limits.
Additionally, we will focus on regions with relatively small κ ∈ [0.1, 0.15]. This is appealing from the standpoint of electroweak baryogenesis, since κ governs the quartic couplings involving the singlet, and a smaller quartic coupling tends to strengthen the phase transition along the corresponding field direction (this behavior is familiar from the SM case [48]).
Taken together with our other parametric choices, this results in the singlet-like CP-even Higgs state being lighter than the SM-like Higgs for the points we consider. We will comment further on this feature below.
Our study will be centered around several representative benchmark points in line with the considerations outlined above. The details of the benchmarks and their associated phenomenology are given in Table I below. These points are in good agreement with all relevant experimental observations, involve relatively low fine-tuning, and will be shown to exhibit dramatically different symmetry breaking patterns in the early universe. Note that the spectra in Table I are quite insensitive to changes in A κ . To a good approximation, varying A κ amounts simply to varying the singlet-like Higgs masses. This can be seen by comparing the spectra of BM 1 and 2, or BM 3 and 4. Both pairs feature an identical choice for the remaining parameters, with the individual points differing only in their values for A κ .
To illustrate the variation of the phase transition properties around the benchmark points, we will vary A κ (and hence m hs ) around the values listed in Table I. For this purpose we define three sets of points, Sets I, II, III, corresponding to the parameter values (except A κ ) for BM 1/2, BM 3/4, and BM 5. We will scan over A κ for these three sets in Sec. IV.
Before moving on to our analysis of the phase transition properties, in the remaining portion of this section we discuss the various phenomenological features of these points, explaining why they are good candidates for beyond-the-Standard Model physics and promising for electroweak baryogenesis. The spectra are computed using NMSSMTools 4.2.1 [28], with the DM properties calculated by MicrOmegas 3.3 [49]. We additionally check the Higgs sector using the HiggsBounds 4.1 package [50]. We require all of our points to pass all relevant experimental constraints implemented in NMSSMTools, MicrOmegas, and HiggsBounds, except those arising from the muon g − 2 and from requiring perturbativity up to the GUT scale. The former can be ameliorated by e.g. including lighter sleptons and the latter by demanding that some new physics enter below the GUT scale [51,52]. Neither will affect any of our analysis of the phase transition or the particle phenomenology. BM 1 and 2 have a Landau pole below M GU T as a result of the large values of λ and small tan β considered in these cases. We discuss the other relevant constraints in more detail below.

A. Neutralino/Chargino Sector
One of supersymmetry's virtues is that it can provide a viable dark matter candidate in the lightest supersymmetric particle (provided R-parity conservation), which in many cases is the lightest neutralino. Consequently, we require all of our benchmarks to have a lightest neutralino LSP in agreement with both direct and indirect detection experiments and consistent with the observed relic density of dark matter [53,54] 0.091 ≤ Ωh 2 ≤ 0.138, (17) where h here is the local Hubble expansion parameter in units of 100 km/s/Mpc. The interval quoted above corresponds to the 2σ limits from the WMAP 9-year data including 10% theoretical uncertainty, and it also encompasses the range suggested by PLANCK data.  A viable neutralino DM candidate is not a requisite feature of SUSY; it is possible to instead have e.g. an axion [55] or gravitino [56] dark matter particle. Requiring the LSP to make up 100% of the observed dark matter density is not essential to our study and can be dropped in favor of another dark matter mechanism 1 . However, one feature that we do find to be important for EWB is relatively low values of µ 300 GeV [21] 2 . This is because µ = λv s / √ 2 and v s must generally be near the electroweak scale for the singlet to significantly impact the EWPT. For typical strongly first-order transitions involving the singlet (in the absence of large supercooling), T n ∼ v s (T n ), and so if v s v, one will typically find v(T n )/T n 1 (for a one-step phase transition). Low values of µ imply the existence of at least two relatively light neutralinos and one light chargino pair. Higgsino-like neutralinos are not viable dark matter candidates, since they are always under-abundant. Thus, requiring a neutralino LSP DM particle necessitates introducing either light bino or singlino states as the LSP.
For the points we consider, the lightest neutralino is bino-like and close in mass to the MeV. This allows for efficient co-annihilation to suppress the thermal relic density, since pure bino dark matter is systematically over-abundant. Consequently, our benchmarks all feature moderately light neutralinos and charginos, with the exception of the wino-like particles, which we take to have masses of 600 GeV so that their effects are largely decoupled from the standpoint of current experimental searches. Since the LSP is bino-like, the cross-sections for spin-independent and spin-dependent scattering of χ 0 1 with nucleons, denoted σ SI , σ SD , respectively, are in agreement with LUX [57], XENON100 [58], and with other direct detection experiments for all the benchmark points we consider. The zero-temperature annihilation rates for the various benchmark points are also compatible with limits from the Fermi large area telescope [59].
Of course, dropping the requirement of a viable neutralino dark matter candidate, the signals predicted for direct-and indirect-detection experiments will be significantly weakened, since the LSP will only make up a fraction of the dark matter abundance.
This set-up is also comfortably compatible with current LHC limits on electroweak-ino production. While the winos are heavy, the Higgsino-like charginos and neutralinos in all of our benchmarks have masses between 200-300 GeV, and decay with almost 100% branching ratio to final states involving χ 0 2 . Since the mass splitting between χ 0 1 and χ 0 2 is small, we have typically BR(χ 0 2 → χ 0 1 γ) ≈ 100%. The resulting photon is very soft and although the decay can result in a displaced vertex, the decay products will simply be counted as missing energy, since such soft photons fall well below the trigger thresholds in the ATLAS and CMS electromagnetic calorimeters. Thus, from the standpoint of LHC searches, the parameter space we explore can be thought of as an effective light singlino-Higgsino scenario, as considered in e.g. Refs. [60,61]. The most relevant LHC constraints are those arising from ATLAS [62] and CMS [63] searches for trileptons with missing transverse energy. A simple application of the limits from Ref. [60] shows that Benchmarks 1-5 lie well within the allowed region of parameter space, although they may be probed at the 14 TeV LHC [61]. have not yet been validated for use with CheckMATE). While the collider phenomenology of the regions we consider here deserves a more detailed study, it is beyond the scope of this work. Here, we simply emphasize that the benchmark points listed in Table I all satisfy current LHC constraints on chargino and neutralino searches.

B. Sfermion Sector
Heavy superpartners are not fundamentally required from the standpoint of electroweak baryogenesis, but rather by their non-observation so far at the LHC. For m χ 0 1 ∼ 100 GeV, as in our benchmarks, ATLAS and CMS currently exclude stops with masses less than about 700 GeV [71-73], albeit with some caveats and potential loopholes. Nevertheless, we take the stop masses to be at or above 1 TeV. All other sfermions are also assumed to have larger masses, set to 1.5 TeV for all benchmarks. Although stop masses above a TeV are already in some tension with naturalness [45][46][47], all of our benchmarks fall within the 10 − 30% fine-tuning range as computed in NMSSMTools 4.2.1 [28,74] and are arguably quite natural in this sense. To quantify the amount of tuning for each of our benchmarks, in Table I we show the value of ∆ max for each point, defined by [74] where p GUT i are the input parameters of the theory at the GUT scale. The amount of finetuning is given approximately by 1/∆ max . For all benchmarks, p GUT max = m 2 hu , except for BM 2, which features p GUT max = A λ . For the sfermion mixing parameters, we choose all tri-scalar couplings to be less than 2 TeV. This, again, takes its cue from naturalness, but also avoids vacuum stability issues that can arise for large A t [75,76]. From the standpoint of the electroweak phase transition in this case, the precise values of the sfermion masses and mixing parameters (other than those for the stops) are immaterial: their effects are decoupled since they are far above the EW scale and because of their small couplings to the Higgs sector. We choose specific values only to concretely demonstrate the phenomenology of the various benchmark points.
The gluino mass is set to 1.5 TeV for all benchmarks. This choice is again motivated by the lack of evidence for a light gluino at the LHC [77, 78] and by naturalness considerations.
The precise value of the gluino mass has little impact on the phenomenology we consider here.

C. Higgs Sector
Of course all of our benchmarks should also be consistent with the Higgs discovery, i.e. feature a Higgs boson in the range 124 GeV m h,SM 127 GeV [79,80] with couplings very similar to those predicted by the Standard Model [81, 82]. More precisely, we require that all of our benchmarks predict signal strengths that match those observed by the ATLAS and CMS collaborations in the various production and decay channels. To this end, one usually defines the signal strength parameter µ(X, Y ) as where σ h,X is the total production cross-section of h via the process X, where X = ggF (gluon-gluon fusion), ttH (associated production with a top quark pair), VBF (vector boson fusion), or VH (associated production with a gauge boson). Here h is taken to be the SM-like Higgs in the NMSSM (the second-lightest CP-even Higgs for our benchmarks), and h SM is where Y = γγ, V V , bb, or τ τ , and similarly for h SM . Typically, the production modes are grouped and analyzed together as X 1 ≡ ggF+ttH and X 2 ≡ VBF+VH (see e.g. Ref. [83] for further discussion). Measurements from ATLAS and CMS, together with those from the Tevatron, can then be used to derive constant likelihood contours in the µ(X 1 , Y )−µ(X 2 , Y ) planes by means of a global fit. This has been done recently in Refs. [83,84].
To check whether or not our points agree with experimental data on the Higgs signal rates, we compute µ(X, Y ) for all relevant channels using NMSSMTools 4.2.1 [28] and verify that they lie within the 95% C.L. regions derived in Ref. [83]. The results for the various benchmarks are shown in Fig. 1, superposed on the likelihood contours obtained from the χ 2 in Ref. [83]. The figure shows good agreement between our points and the results from ATLAS, CMS, and the Tevatron. All benchmarks lie within the 68% C.L. regions and are very Standard Model-like; if the enhanced h → γγ rate continues to decrease in statistical significance, all points will move into even better agreement with the data. We have also cross-checked all of our benchmarks with HiggsSignals 1.1 [85], which takes into account the various correlations and systematics that enter the fit. This can be important in cases when there are multiple Higgs bosons near 125 GeV, as is the case for several of our benchmarks. We again find good agreement with current observation for all points considered.
We must also ensure that the rest of the Higgs sector does not violate the current constraints from LEP, the Tevatron, or the LHC. This is precisely what is checked by HiggsBounds. In all of our benchmarks, the scalar closest in mass to the SM-like Higgs is singlet-like, with couplings of order 10% or less of those for a SM Higgs boson with the same mass. These suppressed couplings make it difficult to detect these states and allow them to be even lighter than the SM-like Higgs, as will be the case for all the points we consider. In fact, in the NMSSM, scalars and pseudoscalars can be extremely light and still compatible with current collider and meson decay limits, provided the mixing is small [86] (see e.g. Ref. [87] for a more detailed discussion of possible strategies to search for these additional states at the 14 TeV LHC). Additionally, all of our benchmarks feature low tan β, making searches for charged Higgs states difficult as well. For all the cases we consider, the strongest bounds on the Higgs sector come from the couplings of the SM-like Higgs and all other Higgs constraints are satisfied by a comfortable margin.

IV. THE PHASE TRANSITIONS AND THEIR PROPERTIES
With the relevant particle phenomenology in hand, we now turn to analyzing the phase transition properties for the various benchmarks described above. First, we define and discuss the parameters of interest from the standpoint of electroweak baryogenesis and cosmology, and detail how we compute them. We then move on to present our results for the various points considered in Sec. IV C.

A. Studying Phase Transitions in the NMSSM
The cosmological phase transitions predicted by the NMSSM are of interest primarily because they may be able to support successful electroweak baryogenesis. If a first-order phase transition is strong enough, it can provide the out-of-equilibrium dynamics necessary to quench the processes which wash out the baryon asymmetry after it is produced, namely the SU (2) sphalerons. The requirement that the SU (2) sphaleron rate be sufficiently suppressed inside the bubble is usually phrased in terms of the order parameter v(T n )/T n , where T n is the nucleation temperature of the bubble (defined in more detail below). Aside from the gauge-dependence inherent in this quantity (discussed in Sec. II C), the correct baryon number preservation condition depends on how large the baryon washout can be. This in turn depends on the strength of the CP-violating sources generating the chiral current, as well as the details of the diffusion of the various charge densities in front of the bubble wall [30]. For example, as a result of these uncertainties, the correct baryon number preservation condition in the Standard Model can range from v(T n )/T n 0.4-1.4, depending on the details of the CP-violation and transport [30].
Following the usual convention, we will define a "strongly first-order phase transition" as a first-order transition such that ∆φ Here, we have defined the slightly more general quantity ∆φ = Still, such a situation may be of interest from the standpoint of cosmological signatures such as gravitational radiation, or extended electroweak baryogenesis scenarios. For the same reason, we will also consider transitions with ∆φ SU (2) = 0 but ∆φ s ≥ 1.
For each of our benchmark points, we therefore calculate the high-and low-temperature VEVs and the difference between these VEVs in both the SU (2) direction (∆φ SU (2) ) and the singlet direction (∆φ s ). In addition to the phase transition order parameter, we will also compute the change in energy density (∆ρ) at the transition, and the change in pressure (∆p). Since the finite-temperature effective potential V (φ, T ) is equal to the free energy density, the change in total energy density is given by where φ 0 and φ n are the (three-dimensional) field values in the false and true vacua, respectively. Note that if there is no supercooling, ∆ρ is the same as the transition's latent heat.
The quantity ∆ρ provides another measure of the strength of the transition, with larger ∆ρ corresponding to more strongly first-order transitions. The pressure difference between the phases is simply the change in the finite-temperature effective potential from the high-to low-temperature VEV: Also of crucial importance from the standpoint of electroweak baryogenesis (and any other microphysical calculation involving the bubble wall) are the details of the bubble wall profile, φ(z) (here z is a spatial coordinate in the frame comoving with the bubble wall and φ is a vector in field space). Several of the CP -violating sources which enter the microphysical calculations of the baryon asymmetry are proportional to dβ/dt [24], and the variation of the Higgs VEVs in the bubble wall is crucial for obtaining non-vanishing CP-violating sources for the chiral current. In the literature, dβ/dt is typically estimated as [24] dβ/dt ∆βv w /L w .
Here ∆β is defined as the change in the angle β = arctan(h u /h d ) from the high-temperature VEV to the low-temperature VEV, L w is the instanton bubble wall widths in the SU (2) direction, and v w is the velocity of the bubble wall in the frame of the plasma far away from the wall outside the bubble. We will thus calculate ∆β, L w , and v w for our various benchmark points, as well as L s , the wall width in the singlet direction, as this quantity will enter the CP-violating sources involving the singlet VEVs.
In addition to the bubble wall properties, we also compute the spectrum of gravitational radiation produced by the strongest first-order transitions. For this we require the relative change in energy density over the transition, α, and the inverse of the duration of the transition 3 , ζ. From Refs. [31,88], we have where ρ rad = g * π 2 30 T 4 n and g * is the number of relativistic degrees of freedom, taken here to be 100. Meanwhile, where S 3 is the Euclidean action, and H is the Hubble expansion rate during the transition.
These quantities allow us to calculate the gravity wave overall amplitude and peak frequencỹ and K is an efficiency factor, given 4 as a function of α in e.g. Ref. [88].
We employ the CosmoTransitions package [27] to study the phase structure and thermal tunneling properties of each of the benchmark points. The full one-loop finite-temperature effective potential (Eq. 7) is directly input into the program along with the zero-temperature electroweak VEV. The minimum at the VEV is traced upwards in temperature (the VEVs are generally temperature dependent) until it either disappears or merges into a distinct hightemperature phase via a second-order transition. Either way, a new minimum is found at the temperature at which the first phase ends, and the structure of the new phase can likewise be found. In this way, the CosmoTransitions package determines the theory's complete temperature-dependent phase structure. In all of our benchmarks, the low-temperature phase is just the electroweak phase determined by the electroweak VEV, and the hightemperature phase resides at the origin of field space. However, different benchmarks exhibit qualitatively different intermediate phases. They can generally lie along either the singlet or h u directions, a mixture of several directions, or be missing entirely (that is, the system can transition directly from the high-T to electroweak phase). We detail the features of individual benchmarks in Sec. IV C below.
Once the phase structure is found, it is straightforward to calculate the critical temperatures T c at which any two temperatures have degenerate minima (corresponding to equal free energy and equal pressure). The CosmoTransitions package then finds the nucleation temperatures T n below T c such that there is unit probability to nucleate one bubble of lowerenergy phase within a higher-energy background phase per horizon volume per Hubble time.
Numerically, to determine the nucleation temperature we use the rough criteria [48] that S 3 (T n )/T n ∼ 140. The bubble action S 3 is found via the CosmoTransitions path deformation routine, and represents the energy of a bubble whose surface tension exactly balances the pressure gradient across its bubble wall.
If the critical bubble is very thick walled, then its center may be significantly displaced from the low-temperature VEV. Once it starts growing, its center will quickly roll down to the bottom of the minimum. Therefore, the wall profile of an expanding bubble may not match that of the incipient critical bubble. To account for this possible discrepancy, we calculate the quantities L w , L s , v w and ∆β for bubble walls moving at constant velocity with constant friction; we dub these "late time" bubble profile parameters for this reason.
Details of this calculation are given in Sec. V, below. An important first test for wall runaway has been suggested in Ref. [89]. Physically, there are effectively two competing forces acting on the bubble wall: the vacuum free-energy density difference between the phases, V (φ 0 , T = 0) − V (φ n , T = 0) ≡ ∆V (T = 0), acting as a driving force, and the pressure exerted by the plasma particles on the scalar field background, F p /A (A here is the area of the wall, neglecting curvature). A necessary, but not sufficient [90,91], condition for a runaway wall is that F p /A < F vac /A = ∆V (T = 0).
The pressure exerted by the plasma in the relativistic limit is given by [89,92] where n i is the number of degrees of freedom for species i, m i (φ) are the field-dependent masses (including the thermal masses for the bosons), E i,p (φ) = p 2 + m 2 i (φ), and where f eq i,p are the equilibrium distribution functions. Since the wall's motion is assumed to be ultrarelativistic, the passage of the wall changes the masses sharply but leaves the distribution functions as they were in the symmetric phase to leading order in 1/γ As pointed out in Ref. [89], the above expression is equivalent to the free-energy density difference between the minima in the so-called mean field T = 0 thermal effective potential, V T [89,92]. V T is simply given by a Taylor expansion of V 1,T around the symmetric minimum in field space, truncated at quadratic order in the field-dependent masses, i.e. [89,92] V T (φ, With this definition, the condition that must be satisfied for the wall to run away can be and defining the full 1-loop mean field effective potential This suggests the following criterion for determining whether the wall is safe from runaway [89]: if tunneling to the broken minimum φ n is not energetically favored in the mean field potential, the wall cannot run away. Thus, for all of our benchmarks, it suffices to compute both the full effective potential and the mean field potential along the tunneling direction; if the symmetry-breaking minimum disappears or is raised above the symmetric minimum in the mean field limit, the wall will remain sub-luminal. We will check against this criterion for all of the points we consider.

C. Results
We can now move on to our main results, summarized in Table II Clearly, the broken minimum is raised above the origin in the mean field limit, and so, by the results of Sec. IV B, the bubble wall must approach a stationary state.
We to the critical bubble profile in Ref. [14]. We also differ significantly from Ref. [14] in our treatment of the stops 5 .
To illustrate the possible range of the one-step electroweak phase transition strength, as well as the parameters of ∆β and L w , we perform a scan over values of A κ and for the three sets of points, Set I, II, III, corresponding to BM 1/2, 3/4, 5, respectively, keeping the rest of the parameter values fixed 6 . As discussed in Sec. III, the rest of the spectrum does not depend sensitively on A κ and so this amounts to varying the singlet-like CP-even/CP- 5 It is also worth noting that Ref. [14] considers a different region of parameter space than we do 6 Although the only benchmarks we consider with a strongly first order one-step EWPT are BM 1 and 2, all three sets of points we considered exhibit one-step electroweak phase transitions for some range of A κ   Fig. 4. The circle-, diamond-, and square-shaped points correspond to Sets I, II, III.
From the standpoint of a one-step electroweak phase transition, the effects of varying A κ are clear: larger |A κ | results in a larger tree-level contribution to the barrier by the κA κ s 3 term in the effective potential. Meanwhile, larger |A κ | results in a smaller m hs for the rest of the parameters fixed. This also explains the differences between the three curves: for a given singlet mass, Set I has the largest |A κ |, while Set II has the smallest. The effect of increasing |A κ | on L w is also rather straightforward to understand: a larger barrier results in a thinner wall (parametrically, L w ∝ ∆φ/ √ ∆V , where ∆V represents an overall rescaling of both the barrier height and the difference in pressure between the two VEVs). From Fig. 4, ∆β is larger for lighter singlet masses/stronger phase transitions. This effect is less simple to understand, as it ultimately results from the complicated interplay of the various parameters in the potential. However, we find a clear correlation between ∆β and the strength of the phase transition. The effect on the wall velocity will be discussed in the next section.
The black points in Fig. 4 do not have a runaway solution for the bubble wall equations of motion, and so may be viable for electroweak baryogenesis, as is the case for BM 1.
However, as advertised in Sec. IV B, for very strong first order transitions, the bubble wall might not be slowed sufficiently by the plasma and can accelerate without bound. Such is the case for the cyan-colored points in Fig. 4. As an illustrative example of this runway case, we can consider BM 2, which corresponds to A κ = −129 GeV in Set I. The transition again proceeds in one step to the broken phase, as with BM 1. The late-time bubble wall parameters are not well-defined in this case, since the wall may never enter a regime with constant velocity and friction.
The large value of |A κ | leads to a very strongly first-order transition, and inspection of the mean field and full effective potentials for this point in Fig. 5 clearly shows the existence of a runaway solution. This suggests that BM 2 likely cannot result in successful electroweak baryogenesis 7 Nevertheless, fast moving bubble walls can be interesting from the standpoint of gravitational wave production. Table III lists the values of α and ζ/H, and the resultant amplitude h 2 Ω and peak frequencyf of gravitational waves produced by this strongly first-order transition, assuming that the bubble does in fact run away (see Sec. IV A for explanation of these quantities). Unfortunately, the relatively low temperature and small α result in a spectrum much too faint to be observed by Big Bang Observatory 7 The existence of a runaway solution does not necessarily imply that the wall will in fact run away. For example, Ref. [90] showed that there can be sources of hydrodynamic obstruction which prevent the bubble from ever reaching the runaway regime. We have checked against the criteria outlined in Ref. [90], and find that the hydrodynamic obstruction is negligible for BM 2. This is because of the relatively large amount of supercooling for this point (T c = 142.0 GeV). Thus, we expect that the bubble wall will indeed run away (although there may be other possible exceptions; see Ref. [91]). (BBO) [93] or eLISA [94]. However,the spectrum of gravity waves we consider is only that coming from collisions of the bubbles and neglects other possibly important contributions from turbulence and other hydrodynamic effects [95][96][97][98]. Future work is required to assess whether our conclusions hold once these additional contributions are properly accounted for.
In addition to the one-step EWPT case, there are several other patterns of symmetry breaking possible in the NMSSM. Although not all are promising from the standpoint of electroweak baryogenesis, they may still yield interesting cosmological signatures. These possibilities are illustrated by Benchmarks 3-5.
Benchmark 3 features a strongly first-order transition in only the singlet direction. This occurs at a nucleation temperature T n = 186.5 GeV, while the critical temperature is T c = 188.7 GeV. The amount of supercooling is significantly less than in the transitions involving the SU (2) field directions. The SU (2) transition, which occurs at T n = 169.6 GeV, is second-order and so is not shown in Table II. This is benchmark has the smallest value of |A κ |, and correspondingly a rather thick wall, with L s = 49.3/T n . As can be seen from the third panel of Figure 5, no runaway solutions exist. We find this to be a generic feature of singlet-only transitions in the parameter space we have investigated: singlet-only transitions tend to have smaller total pressure difference than transitions in all three directions, since only one field is changing its value.
Another novel possibility in the NMSSM is that the phase transition can proceed in two steps: a first transition results in a non-zero singlet VEV, while a second transition breaks electroweak symmetry. This is exemplified by Benchmark 4. Here, the singlet transition occurs at T n = 165.7 GeV, while the SU (2) transition is at a slightly lower temperature, T n = 165.4 GeV. For this particular benchmark, only the singlet phase transition is strongly first-order by our definition (∆φ/T n > 1), while the transition in the SU (2) directions is weakly first-order (has a value of ∆φ/T n < 1). From the discussion in Sec. IV A concerning the baryon number preservation condition, future work is required to determine whether or not such a weak transition can lead to successful electroweak baryogenesis. Regardless, this point serves as a proof of principle that the phase transition in the NMSSM can proceed in two steps. Regarding the bubble profile parameters, the wall width for the first transition is slightly thinner than that of BM 3, corresponding to the larger value of |A κ |. The second transition yields the thickest wall, simply because the phase transition is weak and the barrier between the two minima relatively low. Note that ∆β is significantly smaller for this transition, with a value in the range typical for the MSSM [26]. This is expected, since the singlet is not participating in the transition. Note that because of the relatively low pressure differences for each transition, neither permits runaway solutions, as shown in Fig. 5.
Finally, Benchmark 5 produces a one-step strongly first-order transition at T n = 103 GeV, and is unique among the benchmarks we considered in that h u has a nonzero hightemperature VEV. The symmetry is first broken by a second-order transition in the h u  direction, and is then further broken in the remaining two directions when the first-order phase transition occurs at T n = 103 GeV. Since electroweak symmetry is already partially broken in the phase outside the bubble, sphalerons will already be suppressed in this region, and so it is unlikely that successful electroweak baryogenesis will occur. Also, from the sixth panel of Figure 5, we see that a runaway solution exists, making this transition, like that of BM 2, even less attractive for electroweak baryogenesis, but possibly interesting from the standpoint of detectable gravitational radiation (once again the hydrodynamic obstruction [90] is negligible due to the large supercooling, T c = 141.2 GeV). Table III shows that, like BM 2, BM 5 produces gravity waves with a peak frequency on the order of several mHz and an amplitude on the order of 10 −17 , which unfortunately lies below the range observable by BBO and eLISA [88,93,94].
To summarize the results of this section, we have found that the very narrow region of parameter space investigated features a rich phase transition phenomenology. Transitions can proceed in either one-or two-steps and can occur from either the origin in field space or from a phase with a non-zero VEV prior to tunneling. Some transitions, namely those with a large amount of supercooling, produce bubble walls that can accelerate without bound; the gravitational radiation produced by these bubbles is however too faint to be observed by BBO or eLISA. Most other transitions do not produce runaway bubble walls, and thus can give rise to successful electroweak baryogenesis. The parameters governing the wall profile, such as ∆β and L w , typically take on values more promising for electroweak baryogenesis than in the MSSM, making the NMSSM an even more attractive framework for simultaneously explaining the Higgs mass, dark matter, and baryogenesis.
The only other important parameter not yet computed is the bubble wall velocity. This is the task we turn to next.

V. ESTIMATING THE WALL VELOCITY
We have seen in the previous section that many cases with a strongly first order electroweak phase transition predict bubble walls that approach a finite steady-state velocity.
However, the transport processes important for electroweak baryogenesis typically depend quite sensitively on the precise value of v w . We would thus like to go beyond the analysis above and obtain a quantitative estimate of the v w for the scans shown in Fig. 4.
Determining the wall velocity requires computing the drag force on the bubble wall, which is in general a difficult problem. However, the situation is simplified in two limiting cases: the ultra-and non-relativistic (or "slow-wall") limits. The former is the simplest, since the drag does not depend on the wall velocity or on the deviations from equilibrium of the various species in the plasma at lowest order in 1/γ. The friction saturates and if the driving force is greater than this value, the bubble wall runs away. We considered this limit in our analysis of runaway solutions in Sec. IV B. In this Section we are interested in the opposite case: we will assume that the wall is propagating with velocity such that γ ≈ 1. In this regime we can estimate the friction force microphysically, as has been done for the SM [99] and MSSM [100] some time ago. This will allow us to determine the value of a phenomenological friction parameter Γ that reproduces the approximate wall velocity, provided v w is not too large.
We emphasize that we start with the assumption of a non-relativistic wall in our calculation. Specifically, we assume a simple form for the friction coefficient in the bubble wall equations of motion that does not match on to the solution in the relativistic regime. Despite the shortcomings of this parametrization, it is nevertheless expected to provide a good estimate for the case of non-relativistic bubble walls [91]. Consequently, we only compute the wall velocity for points that do not have runaway bubble wall solutions. Furthermore, our results in this limit may be useful for more in-depth future studies of the wall velocity in the NMSSM, providing preliminary reference values to match on to the techniques of e.g.
Refs. [91,101,102]. In particular, our computation of the various interaction rates of the particles in the plasma is quite general and can be used in future studies independent of the simplifying assumptions outlined above. The reader should interpret our results with the above provisos in mind.
In the sub-sections below we discuss our computation of the stationary-state solution to the bubble wall equations of motion (Sec. V A), as well as the microphysical friction coefficients (Secs. V B, V C; see also App. B). Readers primarily interested in the results may proceed directly to Sec. V D.

A. Simplified equations of Motion
To find the bubble wall velocity in full generality, one must couple the field equations of motion to the radiation, and then also couple the radiation fluid equations to the field. The radiation will tend to slow down the bubble wall's expansion, whereas the forward motion of the wall will tend to heat the fluid and change its velocity. The effective potential is a function of both the field value and the temperature, so heating the fluid will change the free-energy (or pressure) difference between the two phases, thereby changing the bubble wall's acceleration. Additionally, the fluid velocity feeds back upon the wall by changing the magnitude of the drag. The interplay between these effects leads to a rich phenomenology: bubble walls can expand in steady detonations or deflagrations [103], as runaway events [89], or, if reheating is large enough, the two phases can reach pressure equilibrium and there can be an adiabatically varying period of phase coexistence [104].
Solving the full set of relativistic hydrodynamic equations is a difficult problem. However, we are primarily interested in rough estimates of the wall velocity, and are particularly interested in how the velocity compares to that in the MSSM light stop scenario. For this purpose, we ignore fluid velocity and temperature differentials across the bubble wall and their associated effects upon the bubble wall, and we treat the drag upon the wall as fieldindependent. The field equations of motion are then where u µ is the fluid four-velocity, and Γ is the drag coefficient. The bubble will become thin-walled as it grows so we can approximate its motion as planar, and we assume that it reaches a steady velocity and a steady profile. The wall is static in its rest frame, and the equations of motion become The velocity v w is the wall's velocity relative to the fluid. This is the quantity we wish to estimate. As usual, γ = 1/ 1 − v 2 w . If we center the bubble wall at x = 0 in its rest frame, then the boundary conditions are φ(x = −∞) = φ T (the field should be in the true vacuum deep inside of the bubble) and φ(x = ∞) = φ F (the field should be in the false vacuum far outside of the bubble). Also, d is very similar to the equation governing the shape of the initial critical bubble, the only difference being that here the friction term is constant instead of inversely proportional to the bubble's radius. In one field dimension, we can solve for the boundary conditions using an overshoot/undershoot method similar to the one used for the critical bubble. When solving for the critical bubble, the friction term is fixed and the initial position φ(r = 0) is varied to satisfy the boundary conditions (φ(r = ∞) = φ F ). When finding the velocity v w , the initial position is fixed and the overall coefficient Γγv w can be varied to satisfy the boundary condition. If the coefficient is too low, the field will move too quickly (|dφ/dx| will be too large) and it will overshoot the false vacuum. Conversely, if the coefficient is too high, the field will not reach the false vacuum before x = ∞. Instead, dφ/dx will go to zero at finite x, and the field will eventually oscillate about the potential barrier maximum separating the two phases: an undershoot. By iterating between overshoots and undershoots one can converge upon the correct coefficient. We perform this iteration using CosmoTransitions, but modified to vary the drag coefficient instead of φ(r = 0). When the potential has multiple fields, the problem changes in exactly the same way for the critical bubble as for the drag coefficient calculation. We use the pathDeformation module in CosmoTransitions to find the correct path through field space.
With the above approach we can calculate the overall coefficient Γγv w for the various points of interest. These values are shown for the benchmarks without runaway solutions in Table IV. However, we would like to go beyond determining Γγv w and find the wall velocity itself. To do this, we need to calculate Γ. We do so by matching onto a (simplified) microphysical calculation of the friction force in the non-relativistic limit, as described below.

B. Sources of Friction
Friction on the bubble wall arises from interactions of the wall with the various species in the plasma which dissipate the wall's energy. Consequently, to determine Γ in the non- and singlet directions, calculated from the late-time bubble wall profile, for the benchmarks without a runaway bubble wall solution. These quantities enter the calculation of the wall velocities.
relativistic regime, one must know not only the bubble properties, but also the out-ofequilibrium distributions of the various particles in the plasma that interact with the wall. This is generally a very difficult problem, however it has been solved in several different approximation schemes for the SM and MSSM cases in the past (see e.g. Refs. [99,100]).
In this Section, we will apply the techniques laid out in Refs. [99,100] to calculate Γ for our various benchmarks starting from the various interaction rates in the plasma with several simplifying assumptions. Our goal is not a precise calculation of Γ, but rather a quantitative estimate allowing us to compare the wall velocities across our benchmarks and to determine whether the bubble is likely to be sub-sonic or not.
We begin with the equations of motion for a set of infrared (IR) scalar field condensates in a plasma. The scalar field Lagrangian and the conservation of energy result in a set of Klein-Gordon equations with damping terms for the fields φ i . In the bubble wall frame, ignoring the curvature of the wall, we have [99] − where the index j runs over all particles that couple to the scalar field φ i , and δf j connotes the deviation from the thermal equilibrium distribution of the particle species j in the plasma, with in the fluid frame. In the above expression, the δ j are generally spacetime-dependent perturbations from equilibrium generated by the interactions of the particle species j with the bubble wall. Provided that the species j satisfies the WKB condition p j L −1 w (we will deal with soft excitations below), the distributions satisfy Boltzmann equations (semi-classical versions of the Louville equations 8 ), given in the fluid frame by where C[f ] is a collision term which depends upon the various interaction rates of j in the plasma: In the fluid approximation, the quantum mechanical field perturbations are assumed to take the form of a perfect fluid, where the bg subscript corresponds to the perturbations of the background fluid, assumed to comprise the degrees of freedom with small couplings to the Higgs fields. Inserting this form of the perturbations back into the Boltzmann equation yields These equations can be recast in the form [99,100] A where A is a matrix with entries ∼ c 2,3,4 v w , Γ is a matrix involving the various interaction rates, F is the source term, and δ is a vector comprising the various perturbations δ i .
After solving the above set of Boltzmann equations for the perturbations, the (space-timedependent) solutions for the perturbations in the bubble wall can then be plugged into the Higgs equations of motion. Using δf j f 0 δ j , the equations of motion for the background Higgs fields become In general, one needs to solve the coupled set of Boltzmann equations represented by Eq. 41 to determine the perturbations δ i . However, the situation is simplified by noting that, when the wall velocity is small and the wall is not too thin, the terms involving δ i ∼ δ i /L w are multiplied by c j v w and can thus be significantly smaller than the terms involving the δ i (the latter are multiplied by the various rates which are typically of O(10 −2 T ) or larger).
Thus, in this regime, we can approximate the perturbations as roughly constant in the wall, δ = 0. Of course this approximation breaks down for faster moving, thinner walls, but we use it here to obtain a rough estimate to compare between our benchmarks and the MSSM. Ref. [100] compared the friction coefficients found using this approximation with the full numerical solution for the light stop MSSM scenario and found discrepancies up to a factor of 3 for the case of light stops. When discussing the wall velocities for each of our benchmarks we address the implications of these possible differences and find that our overall conclusions remain unchanged.
Under this simplifying assumption, the Boltzmann equations can be easily inverted for the perturbations δ j , resulting in δ = Γ −1 F. Plugging in these solutions, the equations of motion become where the φ i arises from the ∂ t m 2 (φ) term on the RHS of Eq. 40. The η i are (constant) viscosity coefficients which characterize the friction from the plasma on the field direction i of the expanding bubble wall. These coefficients depend on the interactions of all species present in the plasma with one another and with the bubble wall and their calculation is quite difficult. We dedicate Sec. V C below to their computation.
The form of Eq. 44 does not yet match that of the fluid equation (34). In particular, the damping term in Eq. 44 carries additional space-time dependence by virtue of the φ 2 term multiplying the derivative. However, Γ and η can be easily related as follows [106]: let us consider only one field direction for simplicity, with viscosity coefficient η in Eq. 44 and Γ in Eq. 34 (the generalization to multiple field directions is given below). Multiplying Eq. 34 by φ and integrating over x results in Γv w σ = ∆V (45) where ∆V ≡ V (φ 0 , T n ) − V (φ n , T n ) is the pressure difference between the phases and is the surface tension of the bubble wall. Performing the same integration on Eq. 44 after multiplying by φ yields where φ n is the field value in the broken minimum at the nucleation temperature. The last equality follows from assuming assuming a simple shape for the bubble wall profile Finally, combining Eqs. 45 and 47 yields the desired relation between Γ and η for one fielddimension: Thus, given the values of η determined microphysically from the theory, along with the phase transition order parameter, one can determine the values of Γ that enter into Eq. 34.
In our case, we require the generalization of Eq. 49 to multiple fields. In solving the simplified fluid equations analogous to Eq. 34 we assume a common Γ for all field directions.
Of course the friction coefficients η u,d,s are different, since each field couples to different degrees of freedom. However, we can determine an approximate value of Γ that should produce the same wall velocity as that found by solving Eqs. 44. This is done by carrying out the same procedure as for the single-field case, multiplying each equation by φ i , integrating over x, and adding the three equations together, noting that Here ∆V is the pressure which must be balanced out by the friction force for a steady state bubble as before. Setting the wall velocities from both resulting equations equal to one another, we find If one assumes a simple hyperbolic tangent profile for the fields, Eq. 51 can be approximated by Here φ i,0 , φ i,n are the field values before and after tunneling, respectively, L i are the corresponding wall widths, and ∆φ i ≡ φ i,n − φ i,0 . For the simplified cases of a singlet-or SU (2)-only transition from the origin, Eq. 53 simplifies further to where β here is understood to be defined at the T = T n minimum, and we have assumed L u = L d for the second case. To obtain the above result we have assumed that tan β is constant along the tunneling path (i.e. ignored ∆β), while in general it is space-time-dependent.
For our calculations of the wall velocity, we use the full late-time wall profile (via Eq. 51) as computed by CosmoTransitions. Comparing these results with the simple hyperbolic tangent approximation via Eq. 53, we find that the hyperbolic tangent approximation tends to underestimate the full result by a few percent but otherwise is in rather good agreement with the results obtained by using the full profile. This approximation can thus be useful in future studies should the full profile not be computed.
In addition to the damping from particles with p L −1 w described above, infrared bosons will also contribute to the friction [107]. At low momenta, the gauge fields can be treated by classical field theory and be shown to undergo over-damped evolution [108] (scalar fields are not over-damped and so their contribution is numerically much smaller; we ignore these contributions). This gives rise to a viscosity coefficient that can be comparable to those obtained from the WKB contribution above 9 . The relevant Γ terms have simple analytic expressions for the case of one field dimension [108], where m D is the Debye mass of the SU (2) gauge bosons. Note that this expression is valid at leading log order, and hence there is an undetermined O(1) constant that will be added to the log (all such non-logarithmic terms are dropped in this approximation). In keeping with previous work, we neglect this constant term. Eq. 55 generalizes easily to the case of three field directions as before: where x solves L w m W (x ) = 1 (this cuts off the logarithmic divergence in the numerator of Eq. 56 and corresponds to the breakdown of the kinetic theory [108]). Assuming a simplified hyperbolic tangent profile and neglecting ∆β, Eq. 56 simplifies to Here again tan β is the ratio of SU (2) Higgs VEVs at T n . Note that we did not have to solve any Boltzmann equations to determine Γ IR and so it is free of the uncertainty associated with our δ = 0 approximation. In what follows, we take the total friction coefficients Γ to be sums of the Γ W KB and Γ IR , which should be correct up to O(α w ) [103,107].
Armed with expressions (49), (57), and (58), it remains to determine the friction coefficients η u , η d , η s for the fields h u , h d , and s, respectively. This is the task we turn to next.

C. The Friction Coefficients
We would like to evaluate the rates for various particle interactions in the plasma. In particular, the fields most relevant for the drag on the bubble wall are those with large couplings to the scalar fields. In the Standard Model, the fields with the largest couplings to the Higgs are the top quarks, Higgs and SU (2) gauge bosons. In the NMSSM, the couplings of the neutralinos and charginos to the Higgses are also significant. All fields with sizable couplings to h u,d,s can contribute substantially to the friction on the bubble wall.
To compute the friction precisely involves evaluating a large number of interaction rates which are, in general, space-time dependent, due to the changing VEVs. This results in a complicated network of coupled Boltzmann equations represented schematically by Eq. 41.
A full analysis of the drag force on the wall is beyond the scope of this paper; recall that we are already making several simplifications in the treatment outlined above. Thus, in what follows we will make several additional simplifying assumptions and approximations that will allow us to estimate the friction force.
For the SU (2) Higgs fields h u,d we will retain the simplification employed in previous calculations [99,100] and neglect the friction exerted by Higgs bosons. This was justified for the SM case because there there is only one Higgs degree of freedom (not counting the Goldstones) and for the EWPT to be strong required a Higgs lighter than the W mass.
While neither is true in the NMSSM, including the Higgs contribution will only increase the friction. Since we are mostly concerned with the bubble wall moving too quickly for efficient electroweak baryogenesis, ignoring the Higgs bosons in the friction for the h u,d fields will yield a conservative estimate of the drag, which is sufficient for our purposes. Similarly, we will also ignore the neutralino and chargino contributions to the friction in the SU (2) directions.
These approximations are not necessary but simplify our calculations significantly.
We will take the friction on h u,d as arising from the top quarks and SU (2) gauge bosons, neglecting the U (1) contribution. All other species are treated as a common background (see Ref. [99] for details). The SU (2) gauge bosons are treated as one species W with common chemical potential, as in Ref. [99], with the masses, couplings, and rates averaged over the three physical fields. Similarly, the left-and right-handed components of the tops, as well as anti-tops of both helicities, are considered as one species. This corresponds simply to adding the relevant Boltzmann equations together, averaging, and multiplying by the number of degrees of freedom. Thus, the rates simply add to one another. Our set-up is precisely that of Refs. [99,100] to which we refer the Reader for further clarification and discussion.
A novel feature of phase transitions in the NMSSM is that the singlet field s is typically involved. Thus, we need to compute the drag on this field direction as well. Since s is a singlet under all SM gauge groups, its only interactions are with the Higgs bosons, Higgsinos, and singlinos. Once again, a proper calculation would involve many interactions with complicated matrix elements. However, we can once again neglect the Higgs contribution and analyze only those coming from the Higgsino/singlino sector. This will result in a conservative estimate of the friction (within our approximation scheme), since the drag provided by the Higgs bosons will add to that from the Higgsinos and singlinos. The field-dependent masses for the Higgsinos/singlinos (in the gauge eigenstate basis) are given by The contribution of these fields to the friction on s is proportional to dm 2 i (s)/ds and so the contribution of the singlet and singlino to the friction will be suppressed by κ 2 /λ 2 , which is small for our benchmarks. Therefore we can drop these contributions and consider only friction arising from the Higgsinos.
We need to calculate the various interaction rates for the tops, gauge bosons, and Higgsinos which enter the matrix Γ discussed in the previous subsection. The precise definition of Γ in terms of the various interaction rates is given in Ref. [100]. The quantities we need for each species i are given by integrals of the collision term in Eq. 38: Thus, we require the various matrix elements for the processes contributing to C[f ].
The relevant processes include contributions from diagrams with the t-channel exchange of a potentially soft gauge boson. At finite temperature, the contributions from soft bosonic degrees of freedom must be treated with some care, since, naively, they result in infrared divergences in the massless limit. In reality, the divergence is cut off by the thermal self-energy in the propagator, which results instead in a finite, large logarithm ∼ log 1/α. Consequently, the first analysis [107] of the Standard Model case utilized a "leading logarithmic expansion" of the interaction rates, in which only t-channel processes are kept and the external particles were treated as massless. Additionally, Ref. [107] replaced the full thermal self-energies of all massless exchanged particles with the corresponding thermal/Debye masses. The phase space integrals of Eq. 38 were then performed analytically. In Ref. [100], a similar approximation was made, with the additional improvement that the integrals were performed numerically. Subsequent analyses in different contexts [109][110][111] have improved upon these methods, and we draw on these analyses in our treatment of the rates.
Specifically, we consider all leading log order contributions to the various interaction rates, systematically dropping terms of order m/T , as in Refs. [100,107]. However, in the case of the Higgsinos, there are also important s-channel contributions which are e.g. enhanced by N f , the number of fermions or by the top Yukawa coupling and N c . These contributions, although formally higher in the logarithmic expansion, will contribute comparably to the tchannel pieces (see Appendix B). We therefore include these contributions as well. Following the results of Refs. [109][110][111], the matrix elements entering the collision integral are computed in terms of Mandelstam variables at zero temperature, and subsequently 'translated' into the appropriate form for the finite temperature results, i.e. including the thermal self energies in the propagators, using the dictionary provided in Refs. [109][110][111]. The relevant vacuum matrix elements are given in Appendix B. For the evaluation of the phase space integrals in Eq. 38, we use the hard thermal loop approximation for the full thermal selfenergies for all exchanged particles. These are again detailed e.g. in Ref. [110]. We evaluate the phase space integrals numerically using the Cuba Monte Carlo integration package [112].
The evaluation of these integrals is straightforward and is detailed in Refs. [109][110][111], along with the relevant expressions for the thermal self-energies. Further details can be found in Appendix B.

D. Approximate results for v w
With the interaction rates evaluated, we can now compute the wall velocities for all points in the scans shown in Fig. 4, given our assumptions. We only consider points without runaway solutions. The results are shown in Fig. 6. The curves of Fig. 6 demonstrate the expected parametric behavior: stronger transitions occur for smaller values of m hs , which therefore tend to have larger pressure and hence larger wall velocities. The difference between the three sets of points is due to variations in the quantity Γγv w , which is related to the pressure difference driving the wall expansion (see Eq. 45). For a given m hs , Set III has the smallest Γγv w and hence the smallest wall velocity (the drag coefficients η u,d,s are similar for the three sets considered).
We emphasize that the velocities we have computed are estimates that likely underestimate the full result, given by including the δ term in Eq. 41. This was pointed out in Ref. [100] for the case of the MSSM, in which case the δ = 0 approximation underestimated 98  the wall velocity by a factor of 3 for the case of light stops. A full treatment of the Boltzmann and wall equations is necessary to accurately determine v w , and we hope to address this in future work. Additionally, the results shown in Fig. 6 do not match on to the velocities in the relativistic regime. It is thus quite possible that some of the points with smaller m hs are significantly higher, since they should match on to the points with runaway bubble walls.
However, comparing our results to the δ = 0 case in Ref. [100], preliminarily suggests that the effect of friction in the NMSSM is comparable to that in the light stop scenario of the MSSM, and hence that the bubble walls can move sub-sonically (in the non-runaway case).
This in turn indicates that successful transport-driven electroweak baryogenesis can occur in the NMSSM parameter space considered. Even an order of magnitude increase in the wall velocities shown in Fig. 6 would not change this conclusion. On the other hand, neglecting the friction from the Higgsinos typically increases the wall velocities by more than an order of magnitude, highlighting the crucial role played by the Higgsinos in slowing down the bubble wall.

VI. SUMMARY AND CONCLUSIONS
In this paper we analyzed the nature and properties of the electroweak phase transition in the next-to-minimal supersymmetric Standard Model in light of the Higgs discovery at the LHC. We honed in on a region of parameter space featuring significant tree-level contributions to the Higgs mass, a viable dark matter candidate in the lightest supersymmetric particle (the lightest neutralino), stops in the TeV range, and the rest of the particle spectrum compatible with LHC searches. We employed an effective field theory approach to carefully compute the finite temperature effective potential, which was then fed to the CosmoTransitions software package to study the details of symmetry breaking and compute the values of several parameters crucial for the calculation of the baryon asymmetry of the universe.
We showed that the phase transition structure, for phenomenologically viable parameter space points, is very rich. There can be one-step electroweak phase transitions from the origin in field space or from some non-zero VEV prior to tunneling, one-step singlet-only transitions (which might give rise to observable gravitational radiation in some other regions of parameter space), as well as two-step phase transitions. Although some of these possibilities seem not to lead to successful electroweak baryogenesis, such a major event in the history of the universe is certainly interesting cosmologically. We believe it would be worthwhile to study the potential observational probes of such scenarios in the hopes of disentangling the NMSSM (or other theories with additional gauge singlets) from more minimal models.
In addition to computing the patterns of symmetry breaking, a major aim of our study was investigating the microphysical parameters that enter into any realistic calculation of the baryon asymmetry via electroweak baryogenesis. To this end, we studied the bubble wall profile, and in particular the quantities L w , L s , and ∆β. These parameters can vary to up to an order of magnitude even across the narrow slices of parameter space around our benchmarks. Notably, we found that these parameters tend to take on values more promising for electroweak baryogenesis than in the MSSM, further suggesting the viability of NMSSM electroweak baryogenesis.
A crucial part of our study comprised our analysis of the bubble wall velocity for realistic parameter space points. We found that ultra-relativistic solutions to the bubble wall equations of motion exist typically when the phase transition is very strong (when there is substantial supercooling). For weaker transitions the bubble wall velocity tends to be sub-luminal, and thus potentially compatible with successful electroweak baryogenesis. We provided an approximate estimate of the bubble wall velocity in the non-runaway case, hinging upon a microscopic treatment of the various sources of friction acting on the expanding bubble wall. Although this estimate is rather rough in many respects, our results suggest typical values for the wall velocity are in the O(0.01 − 0.1) range, comparable to that of the MSSM light stop scenario. Future work is required to improve this estimate. However, we stress that our computation of the various interaction rates of the particles in the plasma was quite thorough, and can be used in future studies beyond the simple framework we have employed here.
We believe this study can serve as a jumping off point for more detailed investigations of NMSSM electroweak baryogenesis, especially as the LHC, dark matter searches, and intensity frontier experiments continue to clarify what physics might exist beyond the Standard Model. Overall, our results suggest that the NMSSM might not only explain the observed Higgs mass, the nature of dark matter, and alleviate the hierarchy problem, but also explain the origin of the baryon asymmetry of the universe via electroweak baryogenesis.

NOTE ADDED
While the final stages of this work were being completed, Ref. [113] appeared, which investigates phase transitions in a similar region of the NMSSM. They reach conclusions similar to ours in terms of the possible patterns of symmetry breaking, although they do not discuss two-step scenarios with both a first-order singlet and SU (2) transition. Their study differs substantially from our investigation in several respects, especially in their treatment of the stop contributions to the effective potential. Also, they perform large numerical scans of the parameter space, while our study is much more focused on the precise details of the phase transition and bubble wall profile, which they do not attempt to address. In this sense, their study is largely complementary to ours, but quite disjoint.

(A3)
In the above equations, λ and κ are the parameters appearing in the neutrlalino/chargino mass matrices. We do not include their running, since they enter only in the one-loop contributions to the effective potential, and thus their running is formally a higher-order effect. The same reasoning holds for the Yukawa and gauge couplings. It should also be noted that the NMSSM parameters we match onto at the stop scale are technically DR running parameters (appropriate for a supersymmetric theory), while the parameters we use in the effective theory are defined in the MS scheme. Converting between the two schemes results in a small threshold correction to several of the quartic couplings. However, this only affects the quartics at the O(1%) level, and so we neglect these corrections in our calculations.
In practice, to match onto the spectrum calculated by NMSSMTools, when computing the corrections to the 2HD+S potential parameters we consider only the top, gauge boson, Higgsino, and singlino contributions, since NMSSMTools does not include Higgs boson corrections to the lightest CP-even state [28]. We have modified NMSSMTools accordingly so that our spectra exhibit good agreement.
W f → W f : It should be noted that some of the amplitudes listed above do not simplify to those found in Ref. [99]. This is due to some errors in the treatment of Ref. [99] as was subsequently pointed out in Ref. [109]. We have taken these into account and checked that our results for the squared matrix elements do match up to those listed in Refs. [110,114], where applicable.