Detectable Gravitational Waves from Very Strong Phase Transitions in the General NMSSM

We study the general NMSSM with an emphasis on the parameter regions with a very strong first-order electroweak phase transition (EWPT). In the presence of heavy fields coupled to the Higgs sector, the analysis can be problematic due to the existence of sizable radiative corrections. In this paper we propose a subtraction scheme that helps to circumvent this problem. For simplicity we focus on a parameter region that is by construction hidden from the current collider searches. The analysis proves that (at least) in the identified parameter region the EWPT can be very strong and striking gravitational wave signals can be produced. The corresponding gravitational stochastic background can potentially be detected at the planned space-based gravitational wave observatory eLISA, depending on the specific experiment design that will be approved.


Introduction
First-order phase transitions can establish testable links between cosmology and particle physics. This is particularly interesting for the electroweak phase transition (EWPT) that relates to the properties of the Higgs sector which is currently tested at the LHC. Well known examples for this link are electroweak baryogenesis [1] and gravitational wave (GW) production [2][3][4][5].
In the Standard Model the EWPT is not of first order [6][7][8][9] but this behaviour can change in extensions of the standard model (SM). For instance, a strong EWPT is possible in ultraviolet (UV) embeddings that delay the SM-like electroweak symmetry breaking to temperatures below O(10 GeV) [10][11][12]. In most of the cases, however, the EWPT turns to be strong because of new electroweak-scale fields that extend the Higgs sector or induce radiative modifications to it.
Among the plausible UV completions of the SM, supersymmetric theories -which have both an extended Higgs sector and new fields coupled to it -are candidates, that, in principle, can naturally lead to a strong EWPT. In the Minimal Supersymmetric SM (MSSM), detailed studies of the EWPT have been carried out [13][14][15][16][17][18][19][20][21][22][23][24]. It turns out that the measurement of the Higgs mass, on top of LHC constraints on new physics, is cornering the scenario into a parameter region that is in tension with naturalness and collider data [25][26][27][28]. It is hence worth studying the electroweak symmetry breaking dynamics in supersymmetric extensions where the 125-GeV Higgs mass can be accommodated more naturally.
When looking for supersymmetric theories (and their parameter regions) that possibly have a strong EWPT, it is a good guiding principle to start from non-supersymmetric extensions of the SM with a strong EWPT. These non-supersymmetric theories may be viewed -1 -JCAP03(2016)036 as low energy approximations to the supersymmetric ones. The scenario we have in mind is a supersymmetric extension of the SM with most superpartners in the TeV region. The difficulty lies in rephrasing the successful low-energy (electroweak scale) parameter regions in terms of the high-energy (TeV-scale) ones. In practice, one has to ensure that the successful low-energy parameter configurations satisfy the supersymmetric constraints at the high scale. These relations can be strongly modified by large radiative corrections induced by heavy superpartners.
A strategy to deal with these two issues is to analyze the EWPT after employing a "match and run" procedure (cf. e.g. [24,29]). The disadvantage of this technique is that it is precise only when the heavy fields are much above the TeV scale. Since present LHC bounds on new physics do not yet require such large masses [30][31][32][33], in this paper we propose an approach useful to determine the EWPT in scenarios where the radiative corrections of the heavy fields are sizable but not so large as to require a renormalization group (RG) resummation. This method is based on a renormalization prescription that shares some similarities with the on-shell scheme. In particular, we apply this method to determine the EWPT and the corresponding GW production in an illustrative supersymmetric theory with heavy particles. The choice of this theory is motivated by the simple non-supersymmetric SM extensions exhibiting a very strong EWPT.
In many of the SM extensions with a strong EWPT, the barrier that is required between the electroweak symmetric and broken minima of the Higgs potential, is produced by thermal effects or quantum corrections. In these cases, the EWPT typically exhibits small supercooling and latent heat, i.e. is not very strong. The contrary tends to occur in models where the electroweak breaking minimum is close to metastability with the unbroken phase, and the height and width of the barrier between the electroweak symmetric and broken phases is generated by tree-level terms. These terms in fact can easily induce a large jump of the order parameter during the transition and lead to a regime of large supercooling. The most simplistic model along these lines is the SM extended by a real scalar singlet. In its minimal form with a Z 2 symmetry, the tree-level scalar potential of the Higgs and singlet fields, h and s, reads For a suitable range of parameters, this potential can be written as Forλ m andμ 2 s small, the potential has an electroweak broken minimum at {h, s} $ $ EW = {v h , 0} as well as electroweak symmetric minima at {h, s} EW = {0, ±v s } withv s v h /α. The parameterμ 2 s controls the potential difference between the broken and symmetric phase while the parameterλ m controls the height of the potential barrier between the two phases. For a suitable set of parameters, the model exhibits a very strong phase transition while particle phenomenology is in accord with all collider constraints [34][35][36][37][38][39][40][41][42][43][44].
In fact, for some parameter values the potential (1.2) predicts a two-stage phase transition: at very large temperatures the ground state of the system breaks neither the electroweak symmetry nor the Z 2 symmetry. At temperatures around the electroweak scale, the Z 2 symmetry is first broken by a vacuum expectation value of the singlet field ( h EW = 0, s EW = -2 -JCAP03(2016)036 0). At even lower temperatures the electroweak symmetry is broken whereas the Z 2 symmetry is restored ( h $ $ EW = 0, s $ $ EW = 0). It is this second phase transition that is very strong and can lead to cosmological implications.
Motivated by this feature of the singlet extension of the SM, in the present work we study the EWPT in one of its possible UV embeddings: the general next-to-minimal supersymmetric extension of the SM (NMSSM) [45,46]. The aim is to identify the parameter region with a very strong EWPT along the lines of the potential (1.2). 1 That this is possible is not guaranteed, since, as stated above, supersymmetry implies constraints between different couplings of the model, and predicts new particles whose phenomenology may be in conflict with present experimental limits. To accommodate these limits, some fields need to be heavy and their radiative corrections must be kept under control to avoid the destabilization of the tree-level results. We circumvent this issue by employing the subtracting-scheme approach mentioned above.
The paper is organized as follows. In section 2 we introduce the general NMSSM and we identify a parameter region that is promising for a two-step EWPT. Since this region involves some heavy particles (around the TeV scale), in section 3 we review how to deal with the sizable radiative corrections induced by these heavy fields. In section 4 we describe an approach useful to study the EWPT in the presence of heavy fields. Section 5 and section 6 discuss the particle phenomenology of the parameter region identified in section 2, and provide some parameter configurations that are safe from plausible forthcoming LHC constraints. Section 7 contains the analysis of the EWPT performed for some benchmark points, whereas section 8 explains the corresponding gravitational backgrounds and their detection perspectives at the forthcoming space-based gravitational wave observatory eLISA. Section 9 is devoted to some final remarks and conclusions.
2 Tree-level analysis of the phase transition A supersymmetric theory that potentially reproduces the singlet extension of the SM at low energy, is the general NMSSM [45,46]. In this supersymmetric theory all renormalizable interactions respecting gauge symmetry and R parity are allowed. The superpotential involving the superfields of the singlet and the Higgs doublets,Ŝ,Ĥ u , andĤ d , is 1 By a very strong phase transition we mean a phase transition that requires a large supercooling and hence generates GWs with amplitudes detectable at experiments such as eLISA [47]. In sensible supersymmetric models a very strong EWPT was found only before the LHC bounds [48,49]. Further findings on strong (but not very strong) EWPTs in singlet extensions of the MSSM have been reported in refs. [50][51][52][53][54].
We assume all parameters to be real as we are not interested in effects of CP violation. Our notation follows the definitions employed in the public code SARAH [55] which we use to check some of our results. To reproduce at low energy the singlet Z 2 -symmetric extension of the SM, we push the parameters towards the regime where the potential (2.2) resembles (1.2) once the heavy fields have been decoupled. We hence impose the mass of the heaviest CP-even eigenstate to be above the electroweak scale, and the singlet VEV in the electroweak broken phase to be vanishing, s $ $ EW = 0. The latter is generically achieved by a shift-redefinition of S, which helps to identify the relevant terms breaking the Z 2 symmetry.
At this stage, it is useful to fix some of the parameters of the potential by specifying the electroweak breaking minimum. We fix the quantities by rephrasing the electroweak broken phase as with v h = 246 GeV. This implies whereB µ = B µ + λL 1 .
With this set of parameters it is more transparent how to reproduce the regime where the heaviest CP -even and CP -odd Higgses decouple. This is achieved by sendingB µ → ∞ while keeping tan β constant. In the original parameters, for fixed v h this limit implies that also m 2 H d , m 2 Hu → ∞. In general, one cannot impose a Z 2 symmetry on the potential (2.2) with the electroweak minimum (2.4). This would imply µ = 0 which is in conflict with realistic chargino masses (cf. section 6). 2 However, if the Z 2 symmetry is only imposed in the electroweak symmetric phase (i.e. h u EW = h d EW = 0), one obtains the constraints Interestingly, away from the symmetric minimum, the resulting scalar potential is Z 2 symmetric up to a term of the form

JCAP03(2016)036
In the limitB µ → ∞, the light linear combination of h u and h d corresponds to the excitation along (h u , h d ) ∝ (sin β, cos β), while the orthogonal combination is heavy and its VEV goes to zero. Correspondingly, at low energies the term in (2.9) vanishes and the effective scalar potential for the light degrees of freedom displays an (approximate) Z 2 symmetry. The potential in the symmetric phase is then characterized by the coefficients of the quartic and quadratic terms only. In order to keep as many physical parameters as possible fixed, it is also advisable to set the singlet VEV in the electroweak symmetric phase. We adjust m S by parametrizing the electroweak symmetric minima as Using this parametrization, the potential difference between the broken and symmetric phases can be written as Hencev 2 s κ 2 should not be too large in order to ensure that the electroweak broken minimum is the global minimum of the potential. This will be less constraining when considering the one-loop potential, since the loop corrections seem to lower the broken phase (and thereby lift the Higgs mass). At the same time, adjusting the parameterv 2 s κ 2 allows to tune the model to be close to metastability.
Of course, above Z 2 symmetry is not exact. In particular, the heavy scalar and pseudoscalar Higgses still break this symmetry explicitly. This means that quantum corrections will not respect the Z 2 symmetry and neither thermal corrections will (even though the thermal corrections of the heavy particles are in general rather small). To preserve the Z 2 symmetry order by order in perturbation theory will hence require a certain tuning in the parameters.
In any case, since we are interested in the model in the limit of heavy scalar Higgs, very different scales are involved in the calculations. One way to deal with this problem is to use dimensional reduction (DR) regularization, to integrate out the heavy degrees of freedom and then run the parameters down to the electroweak scale using the RG equations of the effective theory. The disadvantage of this method is that it is precise only when the heavy fields are well above the TeV scale (unless also higher-dimensional operators are taken into account, which makes the analysis very cumbersome). Since in our case we do not deal with overly large logarithms, we here propose a different approach: we adopt a regularization scheme where the EWPT quantities are rather insensitive to radiative corrections. Our specific choice will be discussed after briefly reviewing the perturbative problem caused by heavy fields in minimally subtracted renormalization schemes.

Heavy fields in a toy model
As previously explained, it is important to have a reliable strategy to study the EWPT in scenarios with heavy fields. These fields can in fact induce large radiative corrections that  completely spoil the tree-level results, or even generate perturbative problems. 3 The main problem compared to the usual analysis used in collider phenomenology is that we require a scheme in which the perturbative expansion is under control not only close to the broken phase but also in the symmetric phase. Both regions are equally important in the phase transition analysis. In this section we discuss some ideas on how to deal with this issue.
For simplicity we focus on a toy model with one heavy and one light real scalar field. The Lagrangian of the model is where the parameters are renormalized quantities and counterterms are omitted. The field φ is light, whereas the field Φ is heavy, m 2 M 2 . Both m 2 and M 2 are assumed positive such that the fields do not develop any VEV.
It is educative to analyze how the renormalized mass parameter m 2 is related to the mass observable m 2 obs in different subtraction schemes. We consider the one-loop correction to the propagator of the light field φ. Since the one-loop energy in the current toy model is momentum independent, this only amounts to studying the mass shift from the diagrams depicted in figure 1.
Using DR, the naive one-loop relation between the renormalized and the observed mass parameter reads Notably, this relation is problematic when M 2 is large. In particular, for very large values of M 2 m 2 obs ∼ Q 2 the equation has no solution for m 2 at all. This problem can be avoided when Q is chosen such that the loop corrections are small. However, such a strong dependence on the choice of Q is undesirable and finding such a peculiar Q is a problematic task if there are several light and/or heavy fields in the theory.

JCAP03(2016)036
Notice also that this problem is not strictly related to large logarithms. The problem already arises when the loop corrections are larger than the observed mass, which does not necessarily lead to large logarithms in the sense of Only if this relation holds, the resummation using the RG evolution is relevant. Finally, notice that using the naive relation (3.2) in the "matching and run" approach leads to the same problems. Indeed the matching condition between above toy model and the low-energy effective theory (where only the light-field interactions are present) reads where the renormalization scale Q is of order M and the parameterm is the light mass parameter in the effective theory. Since the left-hand side of this matching condition has the size of m obs , the problems above discussed arise also in this case. One way to remove these issues is to reorganize perturbation theory by writing the normalized mass parameter as The sum of the two terms in the brackets is interpreted as the tree-level mass that is used in the propagator of φ. The last term is interpreted as a one-loop counterterm that is grouped with the true one-loop quantum corrections. The resulting relation for the mass turns out to be Hence for the choice one obtains and therefore m 2 obs = m 2 + ∆m 2 Φ + small one-loop correction .  perturbation theory using counterterms that the renormalization conditions are still fulfilled. The difference to the reorganization of perturbation theory using (3.5) is that the tree-level mass parameter is still m 2 rather than m 2 + ∆m 2 . Even though in the current example both procedures are equivalent, the reorganization of perturbation theory in (3.5) has advantages. For example, resummation at finite temperature has to be done using this reorganization in order to avoid temperature dependent bare mass parameters.
That reorganizing perturbation theory resolves the issues indicates that the problem actually stems from the breakdown of perturbation theory. The class of diagrams that have to be resummed is depicted in figure 2. The corresponding contribution is proportional to where we used again the definition ∆m 2 Φ for the loop contribution of the heavy field as in (3.7). Accordingly, the relation (3.2) turns again into Similarly, using the resummed one-loop expression for the light mass, the matching condition (3.4) in the effective theory becomes almost trivial 4 Notice that reorganizing perturbation theory is equivalent to this resummation. In fact, using counterterms does not only correspond to the Daisy resummation on the on-loop level. The counterterm in the reorganized theory will remove all the Daisy diagrams we have resummed in the last paragraph. Resumming relevant (hand-picked) contributions is only one option. After all, it is not guaranteed that the resummed diagrams really contain all contributions that are relevant to improve the convergence of the perturbative series. A more sophisticated approach would be to use the two-particle-irreducible effective action [57]. The corresponding diagrams are depicted in figure 3. This time the propagators are understood to be the full propagators. In the current toy model, the only relevant diagram for resummation is the mixed light/heavy diagram while the other two diagrams only contribute small shifts in the masses. Only considering the mixed diagram then leads to the same result as the explicit resummation above.

JCAP03(2016)036
In summary, only removing divergences and adjusting the mass parameters will not lead to reliable results once the loop contributions to the self-energies approach the order of the physical masses. In this regime, the convergence of the perturbative series is much better in renormalization schemes that use counterterms and are closer to on-shell renormalization. Alternatively, the Daisy diagrams of perturbation theory can be resummed. Since it is easier to implement, we will use the former approach in the following.
A similar discussion applies also to the cases where the light scalar particle obtains a VEV. In these scenarios, the VEV is often used as renormalized (input) parameter instead of the squared mass parameter. The VEV, which is obtained by solving the tadpole equation, is then kept fixed order by order in the perturbative expansion. The loop corrections of this equation contain the mass of the light scalar in terms of its mass parameter and its VEV. Also the heavy-field masses expressed in terms of the VEV, enter in the loops and they induce large contributions to the tadpole equation. To compensate these contributions, the mass parameter has to be drastically modified, with the consequent dangerous effect in the loop involving the light field. On the contrary, if the VEV is kept fixed by introducing a counterterm that cancels the large radiative contributions to the tadpole equation, the issue is mitigated. In particular, such a counterterm amounts to the contribution that keeps the mass parameter fixed (cf. eq. (3.7)) up to tiny radiative corrections as e.g. the one to λ φ (which is small when eq. (3.3) holds). In this regime, therefore, applying counterterms to fix the VEV in the tadpole equations is perturbatively as good as the resummation above. This is the scheme usually employed in NMSSM phenomenology [58][59][60][61][62].
As mentioned at the beginning of the section, when studying the phase transition, it is not enough to ensure the convergence of the perturbation theory in the broken phase. The symmetric phase is equally important and one needs a subtraction scheme in which observables (in particular the VEVs and masses) are stable against quantum corrections also in the symmetric phase. Using counterterms for the three tadpole equations only in the broken phase will not achieve this and lead to sizable corrections to the effective potential in -9 -

JCAP03(2016)036
the symmetric phase at one-loop level. Our strategy will be to introduce counterterms for all soft supersymmetry breaking parameters. These will be fixed by conditions on the effective potential that also involves the symmetric phase. In particular, we would like to preserve (at least approximately) the Z 2 symmetry on loop level. The precise renormalization conditions are discussed in the next section.

One-loop construction
In the one-loop analysis we use the scalar potential (2.2), supplemented by the one-loop effective potential obtained using DR scheme, where n i is the number of degrees of freedom of each species (with a negative sign for fermions) whose squared mass m 2 i is expressed as a function of the h u , h d , s backgrounds fields. Besides, we need to add the effect of the different counterterms These counterterms are considered to be of one-loop order and hence do not change the masses used in the one-loop contributions (4.1). Moreover, they are in accordance with softly broken supersymmetry.
Since we want to study the phase transition of the system, we choose renormalization conditions that do not only hold the VEVs in the broken phase fixed but also the VEVs in the symmetric phase. The counterterms are determined imposing on the full one-loop potential the following tadpole conditions in the electroweak broken phase (2.4) 4) and the following tadpole condition in the two electroweak symmetric phases (2.10) Furthermore, we impose in the symmetric phase This leaves us with one free counterterm that we could in principle use to minimize the mixing between the light Higgs particle and the scalar singlet in the broken phase. In practice this mixing is small due to the Z 2 symmetry on tree level and we use ∆B µ = 0. At first sight, the above choice for the renormalization conditions seems arbitrary and unphysical. However, it turns out that these choices keep almost all masses relatively stable -10 -

JCAP03(2016)036
against loop corrections (cf. section 3). This of course helps in identifying the parameter region where the low-energy theory resembles the singlet extension of the SM with a very strong EWPT. Besides, also the properties of the phase transition are observable, and it is desirable that they are not heavily influenced by loop corrections. Moreover, the above construction reduces the dependence on the renormalization scale Q tremendously compared to the plain (i.e. not resummed) DR scheme.

Tree-level spectrum
The main features of the spectrum can be captured at tree level after imposing the minimization conditions explained in the previous section.
In the basis {h 1 , h 2 , h 3 } with h 1 = s, h 2 = h d cos β +h u sin β and h 3 = h d sin β −h u cos β, the symmetric CP-even squared mass matrix is given by with m 2 A 2 =B µ /(sin β cos β) and m 2 Z = (g 2 1 + g 2 2 )v 2 h /4. A similar procedure can be applied also to the pseudoscalar and charged Higgs squared masses. After a rotation leading to A 2 = h uI cos β + h dI sin β and G 0 = −h uI cos β + h dI sin β (with h uI = Im 2), and omitting the entries corresponding to the Goldstone boson G 0 , the CP-odd squared mass matrix in the basis {A 1 = Im[S]/ √ 2, A 2 } turns out to be Performing the same rotation on the squared mass matrix of the charged Higgses, one obtains a massless linear combination, which corresponds to the charged Goldstones, and the physical charged Higgs, H ± , with mass Notice that in the decoupling limit m A 2 → ∞, the off-diagonal entries in M 2 S and M 2 P are negligible and correspondingly the fields h 1 , h 2 , h 3 , A 1 and A 2 are a good approximation to the (tree-level) eigenstates fields whose masses are the diagonal elements of the matrices. In this case a judicious choice of some parameters, e.g. B S , establishes the (tree-level) mass hierarchy between the singlet-like CP -even and CP -odd Higgses. Moreover, due to the negligible mixing between the singlet-like and SM-like Higgses, radiative corrections coming from sfermions are not expected to change such a hierarchy. These corrections are instead -11 -

JCAP03(2016)036
important for the SM-like Higgs, h 2 . In fact they can efficiently lift m 2 h 2 to the experimental value ∼ 125 GeV. We use this property to set the squark and slepton masses, assuming vanishing trilinear terms and masses all degenerate (although in the regime tan β 10 we will consider, only stops are relevant for the Higgs boson mass and the EWPT). In practice this implies all sfermions are at the TeV scale.
In concrete cases, m A 2 does not necessarily need to be much above the electroweak scale to lead to the decoupling scenario characterized by a low-energy Z 2 symmetry. A reasonable hierarchy of m 2 A 2 versus m 2 Z and λ 2 v 2 h is sufficient to guarantee a SM-like Higgs h 2 without peculiar cancellations in M 2 S,23 . Also the states h 1 and h 3 have a tiny mixing if µv h is small enough compared to m 2 A 2 . Finally, if A 1 is light, its mixing with A 2 is negligible for M S sufficiently small (although A 1 is not required to be small to the aims of the EWPT). In conclusion, A 2 (and A 1 ) can be just above the electroweak scale without peculiar cancellations in the entries M 2 S,23 and M 2 S,13 (and M 2 P,12 ) if the Higgsino mass term µ is (and the singlino mass term M S are) small. This choice would have implications for the neutralino and chargino mass spectra, which can be deduced from the mass matrices the usual NMSSM neutralino and chargino bases are employed and the electroweak breaking minimum of eq. (2.4) is assumed. Since gauginos do not play any relevant role in the EWPT, for simplicity we assume all of them to be degenerate at 1 TeV.

Phenomenological constraints
In the previous sections we have identified a region where the tree-level scalar Higgs sector of the general NMSSM has an approximate Z 2 symmetry at low energy. In this section we analyze the main experimental constraints on this parameter region. This will guide us to select some benchmark points to prove that in some simple supersymmetric models, the present and foreseeable future experimental bounds do not forbid a very strong EWPT with a sizable stochastic GW background. As previously mentioned, for concreteness we fix the stops and gauginos at around 1 TeV. This is in agreement with stop-gluino searches [30,31] and implies enough radiative corrections [63][64][65] to make m h 2 compatible with the SM-like Higgs mass measurements [66].
Further constraints on h 2 come from the measurements on the 125-GeV Higgs signal strengths [67,68]. In general the field h 2 may depart from the SM-like behavior because of two reasons: existence of new decay channels and presence of sizable mixing with h 1 and h 3 . Both issues can be avoided in the identified parameter region. Indeed no dilution of the h 2 branching ratios occurs if m h 1 , m A 1 , m χ 0 1 > m h 2 /2. Moreover, as discussed in section 5, h 2 almost does not mix with h 3 since m A 2 is set to be in the decoupling limit.

JCAP03(2016)036
On the other hand, at tree level h 2 has no mixing with h 1 since M 2 S,12 is vanishing by construction. Radiative corrections generate such a mixing 6 but typically well within the 95% C.L. experimental limit, which requires sin 2 γ < 0.23 [73,74], where γ is defined as . (6.1) As the numerical results of the next section show, in the identified parameter region, tan 2γ turns out to be very small also at one-loop. The tiny mixing of h 1 and A 1 with the Higgs doublets (as well as the fact that the singlet does not acquire a VEV) is essential to overcome also the experimental bounds on h 1 and A 1 . The literature has broadly investigated direct and indirect searches sensitive to non-standard CP-odd and CP-even Higgs bosons [41][42][43][44][75][76][77]. Besides the decays h 2 → A 1 A 1 , h 1 h 1 which are kinematically closed in our case (note that h 3 is heavy and thus rarely produced), the h 1 and A 1 decays into SM particles are tightly constrained by the BaBar, LEP and LHC analyses in a wide mass region of m h 1 and m A 1 . Nevertheless, none of these bounds seem to rule out the identified parameter region with m h 1 , m A 1 80 GeV when M 2 S,13 is kept reasonably below its critical value corresponding to sin 2 γ 0.23 [77]. Finally, unless m A 2 is extremely heavy, the requirement of not too large Higgs-singlet mixing without tuning in the parameters, and the necessary boost of the tree-level mass m h 2 , implies a limit on µ and tends to favor light Higgsinos close to the electroweak scale (cf. eqs. (5.1) and (5.2)). For the same argument, also the singlino mass term M S should not be larger than m A 2 . This however allows for a hierarchy between the singlino and neutralino mass. In particular, if the singlino is sufficiently heavy, as well as the gauginos, the Higgsino states are almost degenerate in mass. Consequently the lighter chargino and the two lightest neutralinos often exhibit multibody decays into pions, which are challenging at the LHC but might be detected at the ILC [78]. In this sense the benchmark points with the Higgsino as lightest supersymmetric particle (LSP) are not problematic with respect to current bounds.
However, from a cosmological point view, the Higgsino as LSP is not fully satisfactory since it cannot account for the observed DM relic density. Among several further DM possibilities (e.g. axion, axino or gravitino), assuming a light singlino can solve the problem. For instance for M S ∼ µ the LSP can yield the correct DM relic density while being safe from direct and indirect detection constraints [79]. The collider constraints on this chargino-neutralino configuration are weak: due to the relatively small mass splitting between χ ± 2 ,χ ± 1 ,χ 0 2 and χ 0 1 , charginos and neutralinos decay via off-shell gauge/Higgs bosons whose products are soft and hence difficult to trigger [80]. Notably, although we will not study numerically the case M S ∼ µ, a priori it seems compatible with the requirement of a very strong EWPT.

Very strong EWPT
In this section we present some numerical results on the strength of the EWPT in the general NMSSM. We are not interested in providing a comprehensive analysis of the ample parameter space but we content ourselves with several benchmark points with a first-order EWPT. In order to avoid the phenomenological bounds discussed in section 6 we impose most of the supersymmetry breaking parameters that are not relevant for the EWPT to be around the TeV scale (in fact exactly 1 TeV if not noted otherwise). This includes gaugino masses, squark and slepton masses. We also choose B µ large to decouple some of the scalar degrees of freedom. For concreteness the parameter M S is assumed sizable, which lifts the singlino mass well above the Higgsino mass with µ = 300 GeV (cf. eq. (5.4)), whereas the squarks trilinear parameters and L 1 are set to zero for simplicity. The slepton and squark masses are assumed degenerate and adjusted to reproduce m h 2 ≈ 125 GeV at one loop. The parameters m S and B S establish the hierarchy between h 1 and A 1 . Since the EWPT is not very sensitive to B S one can fix it to avoid any tension between m A 1 and collider data. Moreover,

JCAP03(2016)036
Hu , m 2 S , T λ , T κ and ξ 1 are fixed to enforce the approximate tree-level Z 2 symmetry as described in section 2.
Concerning the dimensionless parameters, we choose a combination of λ and tan β that provides a relevant boost to the Higgs mass without spoiling the perturbativity of the theory below the grand unified scale [81]. Moreover, κ is set at a small value in order to have the electroweak phase as the global minimum of the potential even in the presence of largev s (cf. eq. (2.12)).
With these choices, the only light particles in the scalar sector (meaning masses much smaller than 1 TeV) are one neutral Higgs, the singlet, one pseudoscalar, two neutralinos and one chargino. The remaining particles, which are heavy, can lead to large one-loop contributions to the potential of the light degrees of freedom. We keep these radiative corrections under control by means of the counterterms discussed in sections 3 and 4.
In the following we analyze four benchmark points in parameter space with phase transitions ranging from weakly first order to very strongly first order and hence close to metastability. Table 1 displays the parameters we use in each scenario, and table 2 (left) shows the corresponding tree-level masses of the light degrees of freedom. Out of these, the bosonic fields are potentially strongly affected by the radiative corrections. However, the Higgs and singlet masses depend on the parameter combination (B S + m 2 S ) while the pseudoscalar mass depends on (B S − m 2 S ). Accordingly, we use this freedom to adjust the pseudoscalar mass freely and its precise one-loop mass is therefore not very relevant for the EWPT. The scalar one-loop masses in turn are not free to chose and are given in table 2 (right). These are obtained (not on-shell but at zero external momentum) from the one-loop effective potential. It turns out that the 125-GeV Higgs has a tiny singlet component (sin 2 γ = 10 −3 ), much JCAP03(2016)036 125.6 sin 2 γ 10 −3 Table 2. Tree-level spectrum (left) and one-loop spectrum and mixing (right) of the four benchmark scenarios quoted in table 1. In the four scenarios the light masses and the mixing angle γ are the same within our calculation uncertainties. below its experimental upper bound, sin 2 γ 0.23 [77]. Also the spectrum is in accord with the collider constraints discussed in the last section and is rather insensitive to changes in v s . This is due to the fact that the dependence of the scalar masses onv s in (5.1) and (5.2) is suppressed by factors of κ, which is small as explained above.
The calculation of the tunneling action S 3 /T is in the current model a multi-field problem. In order to simplify the analysis, we assume a smooth path in field space that passes through the two minima and the saddle point of the scalar potential. The nucleation temperature T n , the ratio of the vacuum energy to the thermal energy of the plasma, α, the quantity β/H = ∂(S 3 /T )/∂ log T , and the strength of the EWPT v h (T n )/T n (where v h (T n ) is the VEV of the Higgs at the temperature T n ) are then calculated by the usual one-dimensional over-/under-shooting procedure [82][83][84]. The resulting tunneling action can be quite large due to the fact that the scalar fields have to cover a relatively long distance in field space. Therefore the scalar potential only requires a rather small potential barrier at the nucleation temperature (see figure 4).
As anticipated, the model displays a very strong first-order phase transition. By increasing the singlet VEV in the symmetric phase,v s , we dial the strength of the phase transition. According to (2.12), changingv s has a direct impact on the potential difference between the symmetric and broken phases. At a certain valuev s 325 GeV the system will not be able to overcome the potential barrier and reside in the symmetric phase. Close to metastability (scenario D), we observe large Higgs bubbles (β/H 10) and sizable latent heat (α 0.1).
In the next section, we discuss the expected GW signal of the benchmark scenarios.

Gravitational waves
In this section we describe the calculation of the GW production in the EWPT of the benchmark points of table 1. For this we have to take into account that the phase transition does not occur in an empty universe, but rather in the presence of a hot plasma of relativistic particles ('fluid'). The calculation will be performed in two ways. Firstly, we will treat the expanding bubbles and the fluid they drag with them as a system of vacuum bubbles only [85]. We model the phase transition as collision of these 'vacuum bubbles', and compute the resulting gravitational wave signal. In a second stage, we model the fluid in a more detailed manner. The phase transition then leads to the creation of sound waves, which in turn produce gravitational waves [86,87]. Finally the resulting gravitational waves signals will be compared to the sensitivities of the experiment that will be possibly launched in the forthcoming ESA Gravitational Wave Mission, (preliminarily) called eLISA.   Table 3. Characteristics of the EWPT in the benchmark points of table 1.

Gravitational wave signal from bubble collisions
The calculation of the production of GWs from bubble collisions will follow [47]. For this we only need three characteristics of the phase transitions in each benchmark point: α, β/H and v w . The former two quantities are calculated in the previous section and quoted in table 3. The last parameter is the expansion velocity of the bubbles. The dynamics of the bubble walls are quite complicated and have been studied in [88][89][90][91][92][93][94][95]. These calculations have been also employed to determine v w in some NMSSM scenarios where the EWPT is not very strong [54]. Fortunately there is a simple criterion to see if the bubble wall approaches the speed of light: it occurs when the effective potential in the mean-field approximation at T = T n is larger in the unbroken phase than in the broken phase [96]. Since this condition is fulfilled for our benchmark points, we can consider v w 1.
If the bubble walls can run away and reach very large gamma factors in the case of very strong phase transitions is still under debate. In a conservative estimate, we assume that bubble acceleration stops long before the bubbles collide. So we consider the case of fast detonations. Let us stress that in this case almost the entire released latent heat will go into the fluid, either into heating the plasma or setting it into motion. The scalar field itself is then completely irrelevant when it comes to the production of gravitational waves, but the envelope approximation might still apply.
Assuming that the dynamics of the fluid can be modelled as that of a scalar field in the envelope approximation leads to a gravitational wave signal with maximum amplitude and peak frequency as [47,85] h 2Ω env = 1.67 · 10 −5∆ κ 2 H β 2 α α + 1 where the parameters∆, κ and f env /β are approximated as Consequently, the spectrum of the GW stochastic background has the shape given by [85] Ω env (f ) =Ω env 3.8(f /f env ) 2.8 2.8 + (f /f env ) 3.8 . Beyond the peak frequency the signal falls as f −1 . The GW spectra corresponding to the four benchmark scenarios are reported in figure 5.

Gravitational wave signal from sound waves
Numerical simulations [86,87] show that the description of the fluid as a scalar field is too simplistic. 7 While colliding scalar field bubbles behave very non-linearly, the fluid dynamics turns out to be very linear and can be described as an ensemble of sound waves. These waves are mostly generated when the bubbles collide, but they do not disappear when the phase transition is completed. They are rather expected to be damped mostly by the Hubble expansion. As a result, the gravitational wave amplitude goes as (H/β) rather than (H/β) 2 . This effect should not be viewed as an enhancement compared to another source. It is just a more accurate description of one and the same source: the fluid.
Notice that the contribution from sound waves has been only simulated up to values α 0.1. Besides, it is not certain how long the sound waves persist in the plasma due to weak shock formation. We consider the envelope approximation and the sound wave approximation as conservative and somewhat optimistic cases of GW production.
The peak amplitude of GW radiation from sound waves is given by [47,86,87] which is larger than the result from the envelope approximation by a factor β/H. The peak frequency isf The scalar field-like description would only apply to the case of strict runaway until bubble collision. and a reasonable fit to the numerical spectrum is given by Beyond the peak frequency the signal falls off as f −4 . The GW spectra corresponding to the four benchmark scenarios are reported in figure 6.

Probing the signals at eLISA
The first GW experiment that might be able to probe the stochastic background in the mHz range is eLISA. This laser interferometry is under discussion in the context of the ESA-L3 Gravitational Wave Mission. The precise experiment architecture is still under debate and several designs, with their corresponding detection performances, are under investigation. In the present analysis we consider three possible experimental scenarios, whose sensitive curves correspond to: two arms of 2-Gm length and five years of data taking (Design 1); three arms of 1-Gm length and five years of data taking (Design 2); three arms of 5-Gm length and five years of data taking (Design 3). It is important to stress that the sensitivity curves for a stochastic background are different from the ones calculated for isolated sources. For the GW signals generated by a first-order phase transition it is more appropriate to consider sensitivity lines based on "power-law integrated curves" [97]. The power-law sensitivity curves corresponding to the three designs above have been determined within the eLISA working groups [98] and applied to study the eLISA capabilities to probe the phase transitions [47]. They are reprinted in figures 5 and 6 (blue lines). 8

JCAP03(2016)036
As figure 5 shows, in the considered benchmark scenarios the bubble collisions during the EWPT produce sizable GW stochastic backgrounds. However the detection of the signal depends on the specific case. In Scenario D, which is very close to metastability, the signal is detectable if eLISA is approved with Design 3. Scenario C is border line with the sensitivity curve of Design 3. Establishing its detectability by adopting Design 3 is delicate as the conclusion depends on the uncertainties/assumptions the sensitivity curves are based on. On the other hand, it seems plausible that improvements in the data analysis can make Scenario C detectable. Scenario A and B are instead too weak for eLISA, independently of the experiment design.
Similarly, figure 6 shows that also the sound waves during the EWPT produce sizable GW stochastic backgrounds. In this case, even Designs 1 and 2 are able to detect the projected signals in some of the scenarios. The scenarios A, B and C are then within the sensitivity of Design 1, while Scenario D is borderline.
In summary, the GW energy released by the EWPT in our benchmark scenarios is large. Still in general it seems unlikely that eLISA can probe the EWPT of the general NMSSM if the (more economic) option Design 1 is adopted.
As final remark, let us emphasize the differences between the present analysis and the ones in [48,49] that also study GW signals in the singlet extension of the MSSM (namely the nMSSM). First of all, in [49] the parameter scenarios with a strong EWPT have been identified via a vast scan over the parameter space (without quadratic and cubic terms in the superpotential) while here we construct them using an approximate Z 2 symmetry as guideline. This approximated Z 2 symmetry, with the singlet having negligible VEV in the electroweak broken phase, avoids problems with the LHC Higgs measurements and constitutes a crucial difference from refs. [48,49] for what concerns collider constraints. Moreover the scenarios with the strongest EWPT found in [49] have somewhat larger α and also larger β than those we produce in our benchmark points. We stress that in our setup we could reach even bigger values of α at the cost of some tuning but we refrain from doing so. A further difference is that ref. [49] is based on two estimates for gravitational waves caused by turbulence (namely [101] and [102]). However, recent numerical analyses [86,87] show that sound waves dominate the dynamics of the plasma and turbulence is probably negligible, at least for frequencies around the peak of the GW spectrum. Last, we use improved sensitivity curves, as published in [47,98], that take into account the stochastic long-lasting nature of the GW signal and the capabilities of the possible eLISA configurations. This sensitivity improvement is the main reason why our conclusion is more optimistic than in ref. [49] for what concerns the observation of GWs in the singlet extension of the MSSM.

Conclusion
We studied the electroweak phase transition in the general NMSSM. We found new regions in parameter space with very strong first-order electroweak phase transitions. In these regions the electroweak minimum is close to metastability of the symmetric phase. The strong phase transition occurs in a parameter domain where one of the Higgs doublets is rather heavy and the scalar potential of the light degrees of freedom displays an approximate (and accidental) Z 2 symmetry. The squarks and sleptons were assumed heavy to easily fulfill all experimental constraints.
The presence of heavy particles leads to large one-loop corrections that have to be renormalized. In the present analysis, we did not only impose renormalization conditions -20 -in the broken phase of the effective potential but also in the symmetric one. In this way we ensured that the approximate Z 2 symmetry is still present at the one-loop level. In general, it guarantees that the properties of the phase transition are not too strongly affected by the loop contributions. Although we applied this renormalization method to a specific supersymmetric model, its application looks promising also for more generic theories where the EWPT depends radiatively on heavy fields.
We also checked the collider phenomenology of the parameter space under consideration. With our choice of parameters, the singlet is the unique new light degree of freedom with a mass around 100 GeV. One pseudoscalar and some of the (Higgsino-like) neutralinos/charginos have masses of about 300 GeV. All the other degrees of freedom beyond the Standard Model have masses of TeV range. After all, the arising spectrum is in accord with present collider constraints.
As anticipated, the model displays a very strong first-order electroweak phase transition in the identified parameter space. In part of this parameter region, the phase transition is so strong that the electroweak broken phase becomes metastable and the system does not tunnel during the course of the universe. Close to this situation, we found very strong phase transitions with large Higgs bubbles and sizable latent heat.
Moreover we analyzed the expected gravitational wave signal. We found that close to the metastable parameter regime (where the latent heat is of order of the radiation energy) the gravitational wave signal is strong. Under some circumstances, the corresponding gravitational wave stochastic background falls into the sensitivity ballpark of the eLISA interferometer. However, the detection seems feasible only if the three-arms eLISA architecture is finally approved.
We stress that the analysis did not aim to be a comprehensive study of the full available parameter space. The aim was to provide a proof of principle showing that in supersymmetric models, striking gravitational wave signals coming from the electroweak phase transition are possible. For this reason a more complete study of the parameter space would be worthy. For instance, we imposed an approximate Z 2 but a priori small departures from this configuration could still lead to sizable gravitational waves. This would also open the possibility of having the heaviest Higgs below the TeV scale, provided one of the MSSM-like Higgses is aligned to the SM one. Moreover, for simplicity we assumed the singlino to be somewhat heavier than the Higgsinos. By modifying this assumption and allowing a singlino close in mass to the Higgsinos (which in principle does not seem in tension with the strong EWPT requirement) would be interesting in order to tackle the dark matter puzzle. Last but not the least, the discovered parameter region should have appealing implications for electroweak baryogenesis. These interesting lines of research are left for future studies.