Electroweak Phase Transition in the $Z_3$-invariant NMSSM: Implications of LHC and Dark matter searches and prospects of detecting the gravitational waves

We study in detail the viability and the patterns of a strong first-order electroweak phase transition as a prerequisite to electroweak baryogenesis in the framework of $Z_3$-invariant Next-to-Minimal Supersymmetric Standard Model (NMSSM), in the light of recent experimental results from the Higgs sector, dark matter (DM) searches and those from the searches of the lighter chargino and neutralinos at the Large Hadron Collider (LHC). For the latter, we undertake thorough recasts of the relevant, recent LHC analyses. With the help of a few benchmark scenarios, we demonstrate that while the LHC has started to eliminate regions of the parameter space with relatively small $\mu_\mathrm{eff}$, that favors the coveted strong first-order phase transition, rather steadily, there remains phenomenologically much involved and compatible regions of the same which are yet not sensitive to the current LHC analyses. It is further noted that such a region could also be compatible with all pertinent theoretical and experimental constraints. We then proceed to analyze the prospects of detecting the stochastic gravitational waves, which are expected to arise from such a phase transition, at various future/proposed experiments, within the mentioned theoretical framework and find them to be somewhat ambitious under the currently projected sensitivities of those experiments.

The customary measure of BAU, Y B , is the ratio of the difference between baryon and antibaryon densities (n B and nB, respectively) and the entropy density (s), i.e., Y B = (n B − nB)/s. Its most precise value to date (Y B = 8.65±0.09×10 −11 ) comes from the measurement at the Planck experiment [4] of the baryon acoustic oscillations that it gives rise to in the power spectrum of the cosmic microwave background (CMB).
For baryogenesis to take place, the much celebrated set of following three Sakharov criteria [5] are to be necessarily met: (i) baryon number non-conservation ( / B), (ii) C and CP violations ( / C,¨C P ) and (iii) departure from thermal equilibrium. EWBG is no exception. However, reference [5] was found to be rather prescient about scenarios based on Grand Unified Theories (GUTs) [6][7][8].
BAU can also be realized in some other motivated extensions of the Standard Model (SM) of particle physics where the same arises from a (s)lepton asymmetry (i.e., via leptogenesis) [9][10][11] in a supersymmetric (SUSY) framework or via the Affleck-Dine mechanism [12,13] or even with the help of gravitational effects [14]. Among all these, EWBG has attracted special attention as it necessarily invokes physics beyond the Standard Model (BSM) down at around the electroweak (EW) scale which is being (and will be) intensely probed at various experiments including at the colliders. Naturally, EWPT, as an essential trigger for EWBG, has continued to be an area of active research [8,[15][16][17][18][19][20].
However, at temperatures as low as around the weak scale, while / B (as anomaly effects [21], via finite-temperature 'sphaleron' transitions) and C and¨C P (induced by the CKM phase) could be present in a scenario with electroweak interactions like the SM, it is difficult to find a departure from thermal equilibrium [18]. EWPT can salvage the situation if it is a first-order phase transition (FOPT) and that also of a 'strong' nature. Such a transition proceeds in steps starting with the nucleation of bubbles of the broken phase in the cosmological plasma of the symmetric phase, followed by their expansions and eventual collisions and mergers until the whole space is engulfed by the broken phase. The process is violent enough to trigger local departures from thermal equilibrium in the vicinity of the walls of the rapidly expanding bubbles in the plasma.
Unfortunately, however, an FOPT (and hence EWBG) cannot be realized in the SM given that the observed SM-like Higgs boson is too heavy (m h SM ≈ 125 GeV [22,23]) for the purpose [24,25]. This is since such a value of m h SM signifies a large enough Higgs quartic coupling (λ H ) which virtually suppresses the term cubic in the Higgs fields in the 'effective' (higher-order) Higgs (scalar) potential. This deprives the potential of a crucial bump (as it varies with the field(s)) which is essential for an FOPT. Also,¨C P from the CKM phase in the SM is proven inadequate for generating enough chiral asymmetries [26][27][28][29] for / B to occur. Hence SM, as such, cannot lead to EWBG.
It is well known that popular SUSY extensions of the SM, viz., the Minimal SUSY SM (MSSM) and its next-to-minimal incarnation (NMSSM), a priori, provide the right setup [30,31] for EWBG. In the presence of an extended Higgs sector in these scenarios and other scalar degrees of freedom (in particular, the top squarks) in their spectra, an effective Higgs potential of the right kind for an FOPT to occur can be found. Furthermore, some new Lagrangian parameters could now be the sources of additional¨C P that triggers / B thus facilitating the generation of BAU. However, since the MSSM parameter space favoring a strong first-order electroweak phase transition (SFOEWPT) has now got highly disfavored (as it requires rather light top squarks [32] which are constrained by the LHC searches), the NMSSM (and its variants) has stolen the limelight. EWPT in the NMSSM is rather appealing because of the presence of a gauge singlet scalar field which helps generate a barrier between the symmetric and the broken electroweak phases of the Higgs potential that is required for an FOPT. Unlike in the MSSM, such a barrier may now arise even at the tree level and at zero-temperature thanks to the presence of cubic terms in the Higgs potential. Thus, the dimensionful couplings in these cubic terms in the soft Lagrangian involving the singlet scalar field and the doublet Higgs fields could play important roles in altering the barrier in favor of an SFOEWPT [31]. In addition, thermal effects including the so-called daisy contributions can also play an important role in the process [33].
Naturally, there has been a continued activity over the past decades exploring myriad aspects and possibilities of EWPT in the NMSSM. In particular, some of these shed light on how SFOEWPT, the experimental constraints on the dark matter (DM) observables and the spectrum of the singlet-and/or doublet-like scalars are connected across the NMSSM parameter space [34][35][36], its region over which simultaneous compatibility of SFOEWPT and the Galactic Centre Excess (GCE) can be found [37] while some others present analyses of EWBG in the presence of SFOEWPT [36,[38][39][40][41][42].
Further, a recent study [43] has undertaken a detailed probe into the patterns of phase transitions, based on calculations of critical temperature, that are possible over an extended region of the Z 3 -NMSSM parameter space using the package PhaseTracer [44] which is designed specifically for the purpose. Subsequently, in the context of such a scenario (in the socalled 'alignment without decoupling limit' in the Higgs sector), it has been demonstrated [45] with the help of the package CosmoTransitions [46] that ensuring a successful nucleation of a bubble of the broken electroweak phase is more crucial than just confirming the presence of a critical temperature for an FOPT. In both the studies mentioned above [43,45], only the Higgs-related constraints from the LHC and the bounds on the chargino-neutralino (electroweakinos) sectors from the LEP experiments are considered. Stringent bounds on the latter sector from the recent LHC studies and those on the DM observables, viz., the DM relic abundance and the DM direct detection (DMDD) rates for both the spin-independent (SI) and the spin-dependent (SD) cases, are, however, not imposed in either of these works. As mentioned there, such considerations are perfectly justified in dedicated studies of EWPT in which neither the properties of these electroweakinos in general, nor those of the DM are of much practical concern.
Going beyond, our goal in this work is to examine the prospects of SFOEWPT in the Z 3 -NMSSM once the latest constraints from the LHC and the DM sector are included in the analysis and their implications thereof. Given that these constraints are already known to be intricately connected over the Z 3 -NMSSM parameter space, such an exercise takes off by throwing the physics of EWPT into the mix. Together, these are likely to shed more light on the viability of EWBG in such a framework. We would, however (as is customary in such studies), remain agnostic about the extra sources of / C,¨C P or how / B is achieved with the understanding that those could always be arranged optimally.
It is perhaps straightforward to imagine [38] that the issues in the DM, the LHC and the EWPT sectors are all connected via the higgsino 'portal'. For, the effective higgsino mass parameter (µ eff ) of the scenario could affect all those sectors significantly, especially, intricately as there are a few other model parameters that appear both in the electroweakino and the Higgs sectors. As we will see, such a connection gives rise to a tantalizing possibility that relatively light higgsinos with masses under a few hundred GeV, in the presence of an even lighter singlino and/or a bino, with or without accompanying singlet-like scalar(s), might have been, somewhat comfortably, escaping their searches at the LHC even at this matured stage of the experiment. In particular, we seek to explore how small a µ eff could still be viable in view of the current experimental constraints given that it is somewhat motivated by 'naturalness' and, at the same time, is preferred by SFOEWPT and hence by EWBG.
On the other side of the proceedings, the dynamics of the nucleated bubbles could generate gravitational waves (GW) [47][48][49][50][51][52][53][54]. These would be stochastic in nature and could be detected by dedicated ground-based and space-borne experiments. Note that in the SM, EWPT is of a cross-over type. Hence it does not produce any GW in the early Universe. That is why the detection of such a stochastic background would likely to hint physics beyond the SM (BSM). In the context of the NMSSM, the production of such GW has recently been studied [38,39,55,56]. In this work, we also present, for a chosen set of scenarios, the prospects of detection of such a GW background in future experiments.
The present work is organized as follows. In section 2 we briefly discuss the Z 3 -NMSSM scenario with a focus on its scalar and electroweakino sectors which the present study is particularly sensitive to. Section 3 summarizes the generalities of EWBG by stressing SFOEWPT as its prerequisite. A schematic details of EWPT is then presented in the Z 3 -NMSSM scenario, matched to the THDSM, in terms of the finite-temperature effective (scalar) potential and the target region of the parameter space for our present study is outlined. A brief discussion on the mechanism of GW production in FOPT follows where we collect its basic theoretical ingredients. In section 4 we present our results where we delineate the relevant region of the parameter space, choose a few benchmark scenarios for our purpose that meet all primary constraints, show that some of these do survive explicit recasts of some recent, relevant LHC analyses (first of its kind, in the current context) and demonstrate in some detail how SFOEWPT is realized in each such case which together underscores an overall preference for a relatively small µ eff . Prospects of detecting the GW at future experiments in these viable scenarios are then presented. In section 5 we conclude with an outlook for the future. A three-part appendix outlines the key details of the implementation of our scenario in CosmoTransitions.
2 The theoretical framework: the Z 3 -NMSSM In this section we discuss the theoretical framework, i.e., the Z 3 -NMSSM with conserved R-parity, by outlining the superpotential, the soft SUSY breaking Lagrangian of the scenario followed by a brief description of its Higgs (scalar) and electroweakino sectors, that are relevant for the present work, at the tree-level.
The superpotential is given by [57] where W MSSM | µ=0 is the MSSM superpotential with its higgsino mass term (the µ-term) dropped, H u , H d and S are the SU (2) Higgs doublet superfields and the gauge singlet superfield, respectively, and 'λ' and 'κ' are dimensionless parameters. The (real) scalar component of the singlet superfield S assumes a non-zero vacuum expectation value (vev) v s during EWPT thus generating an effective µ-term as µ eff = λv s / √ 2. Correspondingly, the soft SUSY-breaking Lagrangian is given by where m S is the soft SUSY-breaking mass of the singlet scalar field, 'S', H u and H d are the doublet Higgs fields and A λ and A κ are the NMSSM-specific trilinear soft couplings with mass dimension one.

The Higgs sector
The tree-level Higgs (scalar) potential of the Z 3 -NMSSM takes the following form: where V F , V D and V soft represent contributions from the F -and the D-terms and the soft SUSY-breaking terms, respectively and are given by where g 2 = (g 2 1 + g 2 2 )/2 and g 1 and g 2 are, respectively, the U (1) and the SU (2) gauge couplings. In our present study, we consider the Lagrangian parameters λ, κ, A λ and A κ to be real. The complex scalar fields can be expressed as where h u = v u , h d = v d and s = v s are the vevs of the real components (CP -even) of the neutral scalar fields that refer to the tree-level scalar potential at zero temperature. Note On electroweak symmetry breaking (EWSB), the doublet and the singlet scalars could mix and the physical Higgs states arise. The tree-level mass-squared matrices for the CPeven, the CP -odd and the charged scalars in the bases {h d , h u , s}, {a d , a u , σ} and respectively, are obtained by expanding the scalar potential of equation 2.3 around v d , v u and v S (see equation A.1) and taking its double derivatives with respect to the scalar fields of the involved types. Diagonalizations of these mass-squared matrices lead to three CP -even, two CP -odd and two charged physical Higgs states. One of the lighter CP -even states has to be the observed SM-like Higgs boson, h SM . Thus, there is the interesting phenomenological possibility that one scalar state from each of the CP -even and the CP -odd sectors is light and is singlet-like (h S and a S ) and hence might have managed to escape detection at various collider experiments. Their heavier counterparts (H and A) would then be similar to those found in the MSSM. Note that these Higgs masses depend on the cubic couplings in which the parameters A λ and A κ appear and these are found to play important roles in achieving FOEWPT. Also, for an FOEWPT (and hence for an EWBG), of particular interest is the effective scalar potential at finite-temperature. We discuss its salient aspects in section 3.

The electroweakino sector
The neutralino sector of the Z 3 -NMSSM consists of five neutralinos which are mixtures of bino ( B), wino ( W 0 3 ), two higgsinos ( H 0 d , H 0 u ) and a singlino ( S) which is the fermionic component of the singlet superfield S appearing in the superpotential of equation 2.1. The symmetric, real 5×5 neutralino mass-matrix, M 0 , in the gauge (weak) basis where M 1 (M 2 ) is the soft SUSY-breaking mass for the bino (wino). The [5,5] element of M 0 is the singlino mass term, m S = √ 2κv S . M 0 can be diagonalized by an orthogonal 5 × 5 matrix 'N ', i.e., when the neutralino mass eigenstates, χ 0 i , are given in terms of the weak eigenstates, ψ 0 j , by and χ 0 i 's are ordered in increasing mass as 'i' increases. In this study, we set M 2 large. Thus, the heaviest neutralino (χ 0 5 ) is almost a pure wino and is indeed heavy with a mass m χ 0 5 ≈ M 2 . Hence the wino practically decoupled when M 0 effectively reduces to a (4 × 4) matrix. The scenario conserves R-parity which is odd for the SUSY excitations. Thus, the lightest SUSY particle (LSP) which is taken to be the lightest neutralino (χ 0 1 ) in this work turns stable and can be a good DM candidate.
The chargino sector of the Z 3 -NMSSM is exactly the same as in the MSSM but for µ → µ eff . The 2 × 2 chargino mass-matrix, M C , in the gauge bases ψ As in the MSSM, M C can be diagonalized by two 2 × 2 unitary matrices 'U ' and 'V ', i.e., where, in the present work, χ ± 1 (χ ± 2 ) is higgsino-like (wino-like) given that we set M 2 large. As we will find, a relatively light singlino-dominated neutralino, which, at times, can be the LSP, has a special context in this work [58]. The latter requires | √ 2κv S | < |µ eff |, |M 1 |. Given v S = √ 2µ eff /λ, this then requires 'κ' to be on the smaller side (with |κ| < 2λ) which is just what an SFOEWPT prefers. Furthermore, the mutual hierarchy among these electroweakinos would have important implications for their phenomenologies at the LHC.

EWPT in the NMSSM: a prerequisite to EWBG and its implications
In this section we take a quick tour into the generalities of EWBG and its association with (FO)EWPT followed by a brief discussion of the latter in the Z 3 -NMSSM. Some relevant analytical details which have gone into our implementations of the scenario in CosmoTransitions are deferred to the appendices.

Generalities of EWBG
Like any successful model of baryogenesis, EWBG also requires the three Sakharov criteria, as mentioned in the Introduction, are to be fulfilled. As noted there, EWBG exploits FOEWPT which triggers electroweak symmetry breaking (EWSB) at the characteristic energy-scale (∼ 100 GeV). In the process, important roles are played by the radiative [59] and finitetemperature [29] corrections to the Higgs potential. For the FOEWPT, the latter ensures an optimal evolution of the potential as the Universe expands and cools down from an early, hot (radiation-dominated) epoch where the electroweak symmetry was still intact [60,61].
A possible FOEWPT is envisaged when there appear (at least) two distinct local minima of the finite-temperature effective Higgs potential, separated already by a barrier when the temperature (T ) of the Universe is such that T > T c , where T c is the so-called 'critical temperature', i.e., the temperature at which the two minima become degenerate, still separated by a barrier. One such minimum is a trivial one with a vanishing potential for null values of the participating scalar fields where the electroweak symmetry is (trivially) preserved. For T < T c , the true (global) minimum emerges with a smaller value of the potential for finite field-values in the broken phase and the field(s) at the trivial (local) minimum (the false vacuum) naturally tries to tunnel to the true one [62][63][64].
The tunneling process is efficiently modeled in terms of a bubble of the broken electroweak phase nucleated locally in the cosmological plasma (in which the electroweak symmetry is intact) that starts growing as the rate of nucleation (Γ B , per unit volume) exceeds the same for the Hubble expansion. A bubble, once formed, continues to expand, collide and coalesce with other bubbles growing in the plasma until a giant one, formed this way, engulfs the whole space thus making EWSB permeate all over. At finite-temperatures (T ), in the semi-classical approximation, Γ B ∝ T 4 exp(−S 3 (T )/T ) [65][66][67], where S 3 (T ) is the effective three-dimensional Euclidean action evaluated at the ("bounce") solution of the classical field equation. The minimal requirement for a successful completion of an EWPT requires the bubble nucleation rate to be one per Hubble volume per Hubble time. This is met when 140 [64,68,69]. The corresponding nucleation temperature T n ( T c ) is the highest temperature for which S 3 (T ) T 140 is satisfied, as the Universe cools down. We use CosmoTransitions [46] to calculate this bounce solution by employing path deformation method.
Along the way, for EWBG to take place, the three Sakharov conditions are met and play their roles [1,8] in the following manner.
• The SM can give rise to / B [21] thanks to the triangle anomaly [70,71]. This is described in terms of the vacuum configurations of the static gauge field of the unbroken SU (2) L gauge theory in which alternating degenerate vacua with integer-value Chern-Simons numbers carry different baryon numbers and are separated by potential barriers whose constant height (E sph ) is given by static solutions that are known as "sphalerons" [72][73][74]. At finite-temperatures / B occurs via sphaleron transitions (hopping of the barriers) [1] from one vacuum to another. The transition rate (per unit volume per unit time) in the symmetric phase scales as T 4 , while in the broken phase, the same is suppressed exponentially as exp( ) [75][76][77]. The same mechanism works in the SUSY extensions of the SM, including the NMSSM [78,79].
• Complementary sphaleron-induced processes would generate similar excesses in baryons and antibaryons thus leading to a null baryon asymmetry. When the underlying theory possesses¨C P , preferential scattering of fermions with a specific chirality in the symmetric phase with the expanding bubble wall could generate both CP and C asymmetries in the particle number densities in that phase thus biasing the sphalerons there to generate more baryons than antibaryons [80,81]. As noted in the Introduction, while SM does not have a strong enough source of¨C P , SUSY extensions like the NMSSM have new sources of¨C P in the form of phases in the extended Higgs sector and/or in the gaugino masses etc., which make up for the deficit.
• Even in the presence of / B, / C and¨C P , the equilibrium average of net baryon number vanishes as a consequence of CP T -invariance [16,17]. Thus, to create a maintainable baryon asymmetry in the front of the bubble wall, the cosmological plasma in its vicinity should depart from thermal equilibrium. Such a departure is generically realized under FOEWPT when the nucleated bubble rapidly expands through the plasma.
• Some fraction of this baryon asymmetry thus generated in the symmetric phase subsequently diffuses into the broken phase [8,19,82] thanks to the motion of the bubble wall. For T < T c (more precisely, for T < T n ), is large in the broken phase and the exponential suppression in the rate of sphaleron transitions, as mentioned under the first item above, kicks in. Quantitatively, for φn Tn ≡ γ EW 1, 1 i.e., for a "strong" FOEWPT, where φ n = φ Tn in the broken phase [68,83,84], this rate per unit volume falls out of equilibrium thus rendering the rate of / B too slow to wash out the baryon-asymmetry that has sneaked into the broken phase. This completes the process of successful baryogenesis.
Given that a FOEWPT is central to the process of EWBG, we briefly review the same in the next subsection in the context of Z 3 -NMSSM.

Study of EWPT in the Z 3 -NMSSM
In this section we outline the formulation of EWPT in the Z 3 -NMSSM and discuss the viability of SFOEWPT that facilitates EWBG over the model parameter space. We assume that there is no spontaneous or explicit¨C P in the Higgs sector.

Effective Higgs potential at finite-temperature
To study the viability of SFOEWPT in the Z 3 -NMSSM, we start with the description of the effective potential for the Higgs sector. The zero-temperature radiatively corrected (at one-loop) effective potential for the (CP -even) Higgs sector is given (in the MS scheme and in the Feynman gauge) by [85] where V tree is the tree-level potential for the CP -even Higgs fields h u , h d and s (see appendix A) and V CW is the well-known Coleman-Weinberg [59] one-loop correction to V tree whose form is shown in the second and the third lines of the equation. There, m j and n j are the field-dependent (MS) masses (see Appendix C) and the degrees of freedom for the species 'j', respectively, and the n j 's are found to be as follows: Note that the scalar states, A 0 i and H ± i , include the Goldstone bosons and that the wino-like states are taken to be decoupled (as pointed out in section 2.2). At finite-temperatures, the 1 In the context of the NMSSM, γ EW = φn Tn = ∆SU (2) (CP -even) Higgs-sector potential receives additional contributions that are (in the Feynman gauge) given by [60,61,86] where the thermal function J B (J F ) captures the relevant thermal contribution from the bosons (fermions), and is given by with the upper (lower) signs appearing for bosons (fermions). This reveals that for m 2 i >> T 2 , i.e., for large |y 2 |, these thermal functions are exponentially (Boltzmann-) suppressed. Hence any massive new physics excitation that has been integrated out from the theory could never have a finite-temperature implication. In the reverse limit, i.e., at high temperatures with |y 2 | << 1, J B/F can be approximated as , γ E being the Euler-Mascheroni constant (≈ 0.577). The term − π 6 y 3 appearing in the high-temperature expansion of J B in equation 3.5a gives rise to a negative contribution cubic in the bosonic field in the finitetemperature effective potential V T . As pointed out earlier, the presence of this term can generate an energy barrier between two degenerate vacua, thus facilitating an SFOPT. Note that such a cubic term appears only for bosonic degrees of freedoms as it comes from the (Matsubara) zero mode propagator which exists only for them. 2 At high temperatures, the perturbative approximations at one-loop suffer from large temperature-dependent contributions from additional higher-order processes given by the so-called "daisy" (or "ring") diagrams. Their dominant contributions to the scalar masses obtained from the resummation of these diagrams are captured in the daisy potential given by [43,[88][89][90] where M 2 h and M 2 V are the eigenvalues of the thermally improved (i.e., Debye-corrected) masssquared matrices of the Higgs and the gauge bosons which are presented in Appendix C [88]. Note that only the longitudinal mode of each of the gauge bosons contributes to the daisy potential of equation 3.6. The one-loop finite-temperature effective potential thus becomes which is then used in the study of EWPT where one tracks its minima as a function of temperature. Its profile for T T c is important for the purpose. However, the locations of the extrema of V T , as well as the ratio φ c (T c )/T c , are both gauge-dependent [60,[91][92][93][94][95] 3 see, for example, references [85,[96][97][98][99]. We have checked the minimization of V T using both Landau and Feynman gauges and have found that the gauge-dependencies of both φ c (T c ) and T c are not significant for the benchmark scenarios we present.
It should be noted here that even for moderately heavy top squarks, which couple intensely to the doublet Higgs fields with the top quark Yukawa coupling, y t , their presence would give rise to large logarithms in V CW (in equation 3.1) because of a large enough hierarchy between m t and mt 1,2 . Such large corrections to tree-level potential point to significant dependence of the results on the renormalization scale 'Q' and a reliable study of phase transition would thus call for treating the potential at higher orders.
To circumvent the problem, one could adopt the effective field theory (EFT) approach in which the top squarks are integrated out from the theory thus resulting in a scenario with two Higgs doublets, a singlet scalar, the electroweakinos and the entire SM spectrum. Thus, the scalar sector of this scenario matches with the one known in the literature [100][101][102] as the (Z 3 -symmetric) Two Higgs Doublet Model with a Singlet scalar (THDMS) extension of the SM. Hence we adopt the tree-level scalar potential of the THDMS in the present work to study EWPT and make use of the relevant results obtained in references [38,43,[100][101][102] where a similar consideration is made. The model parameters of the tree-level scalar potential of the Z 3 -symmetric THDMS are derived in terms of those appearing in the corresponding potential in the Z 3 -NMSSM at the scale M SUSY where the latter is matched onto the former. This correspondence is discussed in Appendix A. We, thus, adopt the following steps [38] to compute the effective potential appropriate for our present study.
• The NMSSM model parameters are taken to be the DR ones at the scale M SUSY (= m Q 3 m U 3 ) following the convention of the spectrum generator NMSSMTools [57] which we use for generating the particle spectrum. • We assume that except for the additional Higgs bosons and the higgsino-, the singlinoand the bino-like electroweakinos, all new physics excitations are rather heavy and hence decoupled. Thus, we use the appropriate set of renormalization group equations (RGEs) which now include contributions from all the states in the THDMS scenario, along with those from these lighter electroweakinos, to obtain the respective THDMS parameters at a reference renormalization scale m t at which the logarithmic contribution from the top quark to the physical minimum of the potential is minimized, and which, also closely resembles the energy scale for the EWSB. V tree (at zero temperature) is then expressed in terms of these parameters of the THDSM at the scale m t . To make the present work self-contained, we present the set of relevant RGEs [38] in Appendix B.
• We then evaluate the zero-temperature one-loop contribution V CW of equation 3.1. Further, we make use of a specific remormalization condition to ensure the dependence of V CW on the renormalization scale (Q) is minimized [103] (see appendix C). Finally, the finite-temperature effective potential, V T , for the CP -even scalar fields, are obtained as described earlier.
We have used the package CosmoTransitions [46] to track the evolution of the finitetemperature effective potential V T and to find T c . Further, the evolution of the potential for T T c has also been studied in order to determine if successful bubble nucleation could occur for our benchmark scenarios. As has been done in some recent studies [43,45], we also study in detail the patterns of phase transitions for some of these scenarios. These pertain to issues like the number of steps taken for the transition to complete, whether it is of a first or a second-order type and the field directions along which a multi-step transition occurs.

Target region of the NMSSM parameter space
In this section we take a brief overall look into what the possibility of an efficient SFOEWPT would imply for the Z 3 -NMSSM parameter space when experimental constraints, in particular, from the observed Higgs sector and from the DM-sector, are also factored in. This leads to our target region of the parameter space from which we choose a few benchmark scenarios to examine their viability against recently reported LHC results on searches of electroweakinos.
For the purpose, it would be instructive to take a quick look into the tree-level NMSSM potential, V NMSSM tree , of equation A.1. Considering only the singlet field, a suitable barrier in the potential profile that makes an SFOPT possible develops when the relative contribution from the trilinear term ∼ κA κ s 3 increases in comparison to the quartic term ∼ κ 2 s 4 . 4 The strength of the transition (parametrized by the ratio of the cubic to the quartic term) increases for a reduced 'κ' since the latter term diminishes faster. Generically, for a given FOPT, increasing A κ (i.e., enhancing the cubic term above) strengthens the same. Furthermore, the term trilinear in the singlet and the doublet scalar fields (∼ λA λ h d h u s) in V NMSSM tree could further reinforce the SFOPT (which can now take place in all field directions) for a suitable A λ . Futhermore, it has been noted [34,36] that an SFOEWPT prefers relatively light singletand doublet-like scalars as these enhance higher order effects in the effective potential. It is also found that thermal effects, including the daisy contributions, could turn crucial in giving rise to coveted barriers between the involved minima.
The upshot is the following. A smaller 'κ' that an SFOPT already prefers leads to a lighter h S . At the same time, this causes a singlino-like state to turn lighter which has crucial implications for the DM and the LHC phenomenologies. On the other hand, to find a S on the lighter side, A κ needs to be so optimally small that it does not make its contribution to the cubic soft term, ∼ κA κ s 3 , insignificant. A near-parallel argument holds for the requirement on the size of A λ which controls the masses of the doublet-like Higgs states and has a somewhat similar role to play for the potential profile via the trilinear term ∼ λA λ h d h u s as does A κ via the terms cubic in 's', as mentioned above.
It may, however, be noted that since h S could mix with h SM on EWSB, a light h S quickly attracts stringent bounds from the experimental studies of h SM . Furthermore, a light h S is also somewhat disfavored by the observed upper limits on the DMDD-SI rates from various DM experiments unless in the presence of a so-called blind spot [104][105][106][107][108] occurring due to a destructive interference among the diagrams with CP -even Higgs states appearing in their propagators. Hence settling for a lone, light a S with a sizable coupling with h SM is a safer option when looking for an SFOEWPT in the Z 3 -NMSSM. Note, however, that regions of parameter space over which h SM could have on-shell decays to h S and/or a S (i.e., when m h S ,a S < m h SM /2) would be highly constrained by the latest LHC data on h SM [109][110][111][112].
As we have discussed earlier, opting for smaller values of 'κ' would, in turn, results in a light singlino-like LSP (m S = √ 2κv S ) which could be the viable DM candidate of the scenario. Apropos of this, as mentioned earlier, a light higgsino-triplet (comprised of a pair of neutralinos and a chargino) resulting from a smaller µ eff = λv S / √ 2 is very much in the context of the present work which, under circumstances, could as well provide the LSP. A larger value of v S can re-introduce the problems with large logarithms from one-loop corrections since the field-dependent masses depend on v S . Keeping this in mind, we consider v S ≤ 2 TeV. Thus, relatively small values of µ eff (O(100) GeV) is achievable for reasonably large values of 'λ'. This, in conjunction with relatively small values of tan β (< 10), helps find m h SM in the right ballpark, even for not-so-heavy top squarks thus letting m h SM appear somewhat 'natural' [113][114][115][116]. We, however, have not restricted ourselves very strictly to this regime and allowed for somewhat larger values of soft masses (m Q 3 and m U 3 ) for the squarks and trilinear coupling (A t ) from the third generation.
For smaller values of tan β, on the other hand, some extra regions of the parameter space could now find compliance with the DMDD-SI constraints by exploiting the so-called 'coupling blind spot' condition g h SM χ 0 5 when |µ eff | tends to approach the LSP mass. This allows us to study a rather nontrivial setup within the NMSSM with a large possible mixing of the higgsinos with the singlino or with the bino. Note that M 1 is not expected to influence the physics of the phase transitions in any drastic way since it enters the calculation of the finite-temperature effective potential via radiative corrections. Hence we have chosen its values (around the electroweak scale) to suit our purpose on DM and collider physics grounds. Thus, an involved situation might arise when all of 'κ', µ eff and M 1 are on the smaller side such that any of the lighter electroweakinos can be dominantly of a particular type or even mixed states. As we will soon find, its implications for the phenomenology of the electroweakinos at the LHC are rather subtle in connection to the physics of both DM and EWPT. It must, however, be noted that since the DMDD-SD rate has the dependence σ SD ∝ 1/µ 4 eff , lowering µ eff beyond a point would quickly attract stringent bounds from the relevant DMDD experiments.
Furthermore, given that we, by now, find that the optimal setup prefers smaller values of tan β, we could afford to consider doublet-like heavy Higgs bosons ('H' and 'A') of the scenario to be on the lighter side and still passing the latest relevant constraints on them from the LHC experiments in the form of bounds on the m H ± -tan β [118] and m A -tan β [119] planes. This is of some importance since relatively light doublet-like Higgs bosons could potentially render the FOEWPT stronger, provided such a light 'H' survives the DMDD-SI constraints. The stage is now set for a brief but important discussion on the phenomenology of such (relatively) light electroweakinos. Dedicated LHC searches for these states over the past years have put stringent lower bounds on their masses and those are becoming even stronger with time. However, these analyses are generally restricted to simplified MSSM scenarios in terms of the spectrum/hierarchy of these states and their consequent patterns of cascades leading to the final states of interest. In general, a scenario like Z 3 -NMSSM could easily invalidate such assumptions in the presence of possible new, light states (for example, the light singlet-like scalars and the singlino). These could then diminish the sensitivities of various target final states and/or tailored signal regions to the experimental analyses thus weakening the lower bounds on the masses of such electroweakinos. In fact, there are myriad such possibilities in our current Z 3 -NMSSM setup which could lead to such a situation [117,120]. On top of that, the LHC experiments mostly assume (at least, the analyses that are are relevant for the present work) the electroweakinos produced in the hard scattering are of wino type for which the relevant cross sections are the largest.
Given the wino decouples from our analysis, for any specific mass the next largest cross section is for the higgsino-pairs which is already about half of that for a corresponding winolike pair. This further reduces the sensitivity of various final states to the experimental analyses. As has been already pointed out, given the central role that µ eff plays in the DM and EWPT sectors, the search for light higgsino-like states has now become of special significance. The bottom line is that the published lower bound on the electroweakino masses are bound to get more relaxed for these higgsino-like states under a situation different from what the experiments assumed for their analyses. However, it is not a straightforward exercise to come up with the relaxed bounds for a given new situation and any such attempt requires thorough recasts of the existing analyses which we will attempt in this work.
On a conservative note, we do not consider possible situations which could result in weakened bounds on the masses of the electroweakinos when these have a compressed spectrum. In our case, for a light higgsino-triplet with a higgsino-like LSP, the lower bound on µ eff could go down to a value as small as ∼ 220 GeV, even with 139 fb −1 of data [121,122].
Guided by the above understanding, we lay down our strategy in section 4 for numerical exploration of the scenario before presenting there our results.

Production of GW from first-order phase transition
Given that the NMSSM could provide us with an ideal setup for an FOPT that might have taken place in the early Universe, a study of GW originating from such an FOPT, in the context of our present work, is in order. As noted in the Introduction, GW from an FOPT would exist in the form of a stochastic background and has been proposed to be searched for using the so-called "cross-correlation" method [123][124][125][126][127]. The salient mechanisms via which GW could arise from an FOPT and their corresponding contributions to the GW energy density (scaled by the critical density ρ c for the Universe with a vanishing cosmological constant Λ) are as follows.
• Collisions of the expanding bubble walls release stress energy located at their walls, as well as lead to possible subsequent shocks, in the intervening plasma made up of relativistic particles [128][129][130][131][132][133]. However, for a phase transition occurring in a thermal plasma, their contributions to GW energy density are believed to be negligible [134] and hence can be ignored.
• Bulk motion (velocity perturbations) of the plasma generates sound (acoustic) waves (longitudinal modes) that propagate in the same during the time interval between collisions of bubbles and the expanding new phases dissipating their kinetic energy in the plasma [135][136][137][138]. These sound waves contribute to the GW energy density as [139], with H 0 standing for the present-day (red-shift z = 0) value of the Hubble parameter, also known as the Hubble constant. Such acoustic contributions, when accumulated over the said duration, are expected to dominate.
The overall GW energy density can be approximated as a linear combination of the latter two contributions, i.e., A few key FOPT parameters, in addition to the bubble nucleation temperature, T n , that can be obtained from the particle physics models and which control these two contributions can be categorized as follows.
• The parameter 'α', which relates to the energy budget of the FOPT, is given by [146] where T * = T | t * with t * being the instant of time when the FOPT completes. In the absence of significant effects from reheating, is the difference between the potential energies at the false and the true minima and ρ * rad = g * π 2 T 4 /30 where g * is the number of the relativistic degrees of freedom at T = T * , taken here to be ∼ 100, following recent literature.
• The parameter 'β', which gives the inverse time-duration of the FOPT, can be derived in terms of the effective 3-dimensional Euclidean action (S 3 (T )/T ) as [131] where H * = H| T * . For a stronger GW signal, the EWPT should occur over a larger duration of time, i.e., it should be a slow process and hence the ratio β/H * needs to be on the smaller side.
• The parameter v w , which pertains to the bubble dynamics, i.e., the wall-velocity of the expanding bubble, needs to be larger for a more intense GW emission, although an optimally strong EWBG is known to be favored only for a tiny, subsonic v w instead.
The sound wave contribution to the GW energy density, Ω sw h 2 , as a function of the above FOPT parameters and the frequency 'f ' of the GW, is then given by [135,147,148] (3.11) where κ v is the fraction of the energy from the phase transition that gets converted into the bulk motion of the plasma which leads to GW and is of the form [54,149] κ f sw is the present day peak frequency for the sound wave contribution to GW energy density given by (with the approximation T ≈ T n ) [133] f sw = 1.9 × 10 −5 Hz Υ(τ sw ) is the parameter that brings in the effect of a finite lifetime of the sound waves which suppresses their contributions to the GW energy density and is given by [150,151] where the lifetime τ sw is considered as the time scale when the turbulence develops and is given by κvα 1+α is the root-mean-squared (RMS) fluid velocity obtained from a hydrodynamic analysis [51,153,154]. Note that as τ sw → ∞, Υ → 1, asymptotically. On the other hand, for all our benchmark scenarios presented in section 4.3.2, τ sw H * < 0.1 when Υ → τ sw H * . Furthermore, there is a growing realization [155] that v w might not enter the calculation of the EWBG. Then, to maximize the strength of the GW, it is assumed that the expanding bubbles attain a relativistic terminal velocity in the plasma, i.e., we consider v w 1.
Note that in the above calculation, the estimation of the portion of energy transferred to the fluid motion is based on the so-called bag model [146]. A recent work [156] proposes a model-independent approach (by going beyond the bag model) to obtain this quantity. In that work the parameter αθ quantifying the strength of the phase transition is given by where the subscript s (b) corresponds to the symmetric (broken) phase,θ is the difference in energy (e) and pressure (p) in a given phase and known as the pseudo-trace, Dθ ≡θ s (T s ) − θ b (T s ) is the difference in its value in the symmetric and the broken phases, c s b being the speed of sound in the broken phase which is defined as while ω s = (e + p) s is the enthalpy density in the symmetric phase.
The GW power spectrum due to sound wave from beyond the bag model can then be obtained from equation (3.11) by just carrying out the following replacement: Subsequently, the GW spectra, within and beyond the bag model, are compared in section 4.3.3. Furthermore, the MHD turbulence contribution to the GW energy density is given by [123] , (3.18) where k turb is not precisely known but is expected to be in the range of 5%-10% of k v [138]. We set k turb = 0.1k v in our calculation. The present-day peak frequency f turb of the GW spectrum from the turbulence contribution is given by with h * = 16.5 × 10 −6 Hz T n 100 GeV g * 100 1 6 .

Results
In this section we start by presenting the ranges of various input parameters of the scenario that we adopt to carry out a scan over the theory space. This is followed by a brief discussion on the pertinent constraints coming from relevant DM and collider experiments including the

Constraints from various sectors
In this work, we take into account constraints from various sectors, both theoretical and experimental. The theoretical ones include ensuring the spectra to be free from tachyonic states, the scalar potential not develop an unphysical global minimum and the evolutions of various pertinent couplings of the theory with energy not encounter Landau poles, etc. The experimental constraints include those coming from the Higgs, the DM and the flavor sectors and from various searches for new physics at the colliders. We further ensure the occurrence of SFOEWPT that facilitates EWBG and that such a transition does end up in the physical vacuum. To impose these constraints and for our general numerical analysis, we employ publicly available packages like NMSSMTools (v5.5.3) [157,158], HiggsBounds (v5.8.0) [159], HiggsSignals (v2.5.0) [160], CheckMATE (v2.0.34) [161], SModelS (v2.1.1) [162] and CosmoTransitions(v2.0.6) [46]. Below we briefly point out some of the important constraints that are obtained from these packages.
• NMSSMTools is used to compute and constrain various relevant observables from the Higgs, the DM, the flavor and the collider sectors. We impose the 2σ upper limit on the DM relic abundance, i.e., Ωh 2 ≤ 0.131 as reported by the Planck experiment [4,163]. The most recent (and improved) upper bounds on the DMDD-SI [164] and -SD [165,166] rates are taken into account after a commensurate downward scaling of these crosssections (as the relic abundance drops below the Planck-allowed band) is done. This helps the computed DMDD-SI  also takes into account, albeit simplistically, the constraints from the CMS analysis on the electroweakino searches in the 3 + / E T final state with 35.9 fb −1 worth of data [168].
• Using HiggsBounds and HiggsSignals, we retain only those parameter points which pass the thorough checks of the Higgs sector. With the help of the latter package, we allow for Higgs signal-strengths which are consistent with the experimental findings at a 2σ level. To take into account the theoretical uncertainties in the computation of m h SM , we consider m h SM over the range 122 GeV < m h SM < 128 GeV.
• A few representative (benchmark) scenarios out of the resulting set are then subjected to thorough recasts, via the packages CheckMATE and SModelS, of a multitude of relevant LHC analyses that include several recent ones with 139 fb −1 of data. These analyses and their availabilities in these two packages are indicated in table 2. Together, these are expected to provide us with the most stringent lower bounds on the masses of the electroweakinos under diverse circumstances which are pointed out while discussing those. In addition, there are a few more rather recent LHC analyses [121,184,185] which are expected to be sensitive to the scenarios we study but yet not available in the public versions of either of these two packages. We will get back to these in section 4.3.2.
• Parameter points that pass the previous set of constraints are subjected to analyses via CosmoTransitions to check for SFOEWPT that results in the physical EW vacuum.
We, however, do not consider the recent experimental finding on muon (g − 2) [186,187] since the dust is yet to settle over its BSM implications. We, thus, have set the masses of the smuons, along with all the sfermions, at a multi-TeV range. When using CheckMATE, we have generated, for each such analysis, Monte Carlo events for the leading order productions of all pertinent pair and associated productions of various electroweakinos at the 13 TeV LHC, i.e., for pp → χ j χ k , (χ j,k ∈ {χ 0 i , χ ± 1 }, with i ∈ {1 − 4}), with up to two additional partons, using MadGraph5 [188]. These events are then passed through PYTHIA8 [189] for generating parton showers, hadronization and decays of the unstable particles. The additional partonic jets from the matrix elements are then matched to those from the parton showers (the so-called ME-PS matching) using the MLM prescription [190] built-in in MadGraph5. The resulting events are passed through DELPHES [191] to include the detector effects. For an analysis using SModelS, we just provide the package with the SLHA file along with the MadGraph5-generated cross sections of various production processes as mentioned earlier. To account for the significant NLO+NLL contributions, all production cross sections have been multiplied by a flat k-factor of 1.25 [192]. Both the recast packages calculate a r-value for a given theory point, where r = (S − 1.64∆S)/S95, with 'S', ∆S and S95 signifying the predicted number of signal events, the associated Monte Carlo error and the experimental limit on 'S' at 95% confidence level, respectively. Nominally, r < (>)1 indicates the scenario to be allowed (disallowed).

Choice and study of benchmark scenarios
As pointed out earlier, we now look for a few benchmark scenarios from those that pass the selections of NMSSMTools, HiggsBounds and HiggsSignals. In figure 1 we present scatter plots of parameter points that pass those selections in the plane of |µ eff | − m χ 0 1 . The choice of the said plane is motivated by the physics of the relatively light electroweakinos that are in the context given the recent LHC searches and from the viewpoint of SFOEWPT. Presenting the bino (left plot) and the singlino (right plot) contents of the LSP (in the palettes) further clarifies the situations from the involved angles.
In both plots, scenarios having a higgsino-dominated LSP arise, by construct, along the diagonals (m χ 0 1 ≈ µ eff ). Points along the two horizontal streaks appearing at low m χ 0 1 correspond to a bino-or a singlino-dominated LSP DM that find h SM and Z-boson as funnels in their mutual annihilation. The sparse occurrence of a singlino-dominated LSP over these streaks points to some amount of tuning that is needed among the NMSSM parameters to comply simultaneously with the constraints from the Higgs and the DM sectors, an issue which is not of much concern for a bino-dominated LSP since M 1 could be altered practically freely without affecting the Higgs sector. m χ 0 1 30 GeV is disfavored since as a DM candidate χ 0 1 would require a relatively light Higgs boson (a S or h S ) below ∼ 60 GeV for an efficient (funnel) annihilation which, in turn, attracts severe constraints from the studies on h SM decays. Furthermore, the DMDD constraints are rather severe for such m χ 0 1 .
In each of these plots, another densely populated region appears along the edge of the diagonal where efficient coannihilations of the DM with closely lying electroweakinos, backed by favorable mixing among these states, are possible. In the rest of the (less populated) regions, compliance with the upper bound on the DM relic abundance is facilitated mainly by various Higgs boson funnels. Derth of points over the region bounded roughly by 100 GeV < m χ 0 1 < 200 GeV and 100 GeV < |µ eff | < 250 GeV is due to the constraints derived from the CMS search for electrweakinos in the final state 3 + / E T with 35.9 fb −1 of data [168]. A similar observation was made in reference [193] which finds further support in subsequent studies [120,194,195]. A low population of points at higher |µ eff | and for intermediate values of m χ 0 1 is mostly since the DM tends to be over-abundant due to its sub-optimal conannihilation rate and/or for a lack of suitable annihilation funnels.
In the subsequent subsections we settle for a few benchmark scenarios out of these allowed set which are representative of various situations of interest. We study their properties related to phase transitions at finite temperatures to ensure that an SFOEWPT (i.e., γ EW 1) could occur. It is important to note that while relevant LHC analyses with ∼ 36 fb −1 of data would continue to constrain our scenarios, the entire region of the NMSSM parameter space indicated in figure 1 can now be sensitive to some of the recent LHC searches for the electroweakinos with 139 fb −1 of data which all are listed in table 2. Hence we subject the benchmark scenarios to these analyses via their recasts using CheckMATE and SModelS. Furthermore, we check the future experimental sensitivity of the GW produced during the time of phase transition for a few of such allowed scenarios.

Studying the benchmark scenarios
In this subsection, we discuss how searches for the lighter electroweakinos at the LHC could restrict the region of parameter space which otherwise favors SFOEWPT and satisfy all other experimental bounds. The value of µ eff is in direct reference since a small value of the same, while favors SFOEWPT, draws substantial constraint from the above-mentioned searches. Our goal is to first identify the lowest values of µ eff (or, for that matter, the smallest values of the higgsino-like electroweakino masses, except when these form a triplet which contains the lightest of all the electroweakinos) that would be allowed under different circumstances, followed by a discussion of allowed scenarios with optimally light higgsinos. In the process, we highlight the role of different signal regions that play crucial roles.
For a scenario with a decoupled wino, lower bounds on the masses of the lighter electroweakinos to be derived from the LHC experiments would crucially depend on the values of two quantities, µ eff and κ/λ. This is all the more so for a smaller κ (which favors SFOEWPT) that renders the singlino lighter. With µ eff not so large, electroweakinos, with their dominant contents, could exhibit altered hierarchies in their masses which result in contrasting patterns in their cascades. These, depending on their mutual mass-splits, result in their altered sensitivities to different final states and/or signal regions at the LHC experiments.
A smaller M 1 could add further intricacies to the collider phenomenology [117,120,194,195] by placing the bino-like neutralino in the vicinity of the light singlino and the higgsinos while aiding compliance with various constraints from the DM sector. We also stick to small values of tan β ( 5) which favors SFOEWPT. As noted earlier, we further consider relatively large values of 'λ' ( 0.5) which, in conjunction with small tan β values, aid compliance with the observed value of m h SM in a more 'natural' way. Note that we seek to allow for relatively small values of m H as well since those are what is preferred by SFOEWPT. For small values of tan β (1 tan β 5) that we would like to restrict ourselves to, stringent lower bounds on m H ± ,H,A come dominantly from their searches in the tb (for H ± ) and τ τ (for H, A) final states [118,119]. While in the general scan of the parameter space, these constraints have eliminated some scenarios, the benchmark scenarios that we work with happen to lie outside the constrained regions. However, for the latter, as and when these become sensitive to similar future analyses, the presence of light electroweakinos in the spectrum could help evade those if the Higgs states could also decay to these electroweakinos. In this regard, searches for H ± is expected to be of immediate relevance and hence we mention its altered branching fraction to tb final state for our benchmark scenarios. Also given that the patterns of vacuum transitions get to be rather involved, we adopt the following convention to describe those in the upcoming discussions. The total number of steps involved in a given phase transition process is denoted by the roman numerals (i.e., I, II, etc.) while the type of the phase transition, i.e., whether it is of a first or a second-order kind, is denoted by the arabic numerals (i.e., 1, 2, etc.). On the other hand, for multi-step phase transitions, the various field directions along which the phase transitions occur are indicated by 'S', for the singlet direction and 'D', for the SU (2) field directions. For example, a direct, i.e., a one-step, FOPT along all three directions is denoted by 'I(1)', whereas a two-step FOPTs in which the first transition takes place along the singlet direction and the subsequent one along the SU (2) field directions is labeled as 'II-S(1)-D(1)'.

Disallowed scenarios with low µ eff
Benchmark points presented in table 3 are chosen with the following considerations. We seek to get an idea of how large a value of µ eff which is still consistent with SFOEWPT but is expected to be ruled out by the electroweakino searches at the LHC. We employ CheckMATE for the purpose by putting all its currently implemented set of LHC analyses in action. In case we find some such benchmark scenarios to be barely allowed, we subject the same further to SModelS (see table 2) in which a host of relevant LHC analyses (with 139 fb −1 of data) are incorporated to check if that is still the case.
In BP-D1, the LSP is singlino-dominated with m χ 0 1 ∼ 60 GeV. The higgsino-like electroweakinos are the immediately heavier states with their masses governed by |µ eff | (∼ 275 GeV) and range over 280 GeV -310 GeV. We set M 1 ∼ 480 GeV such that the binodominated neutralino effectively decouples on both collider and cosmology considerations. The DM relic is under-abundant thanks to the presence of the h SM -funnel which, 'λ' being large (=0.68), is rather efficient. This effectively scales down the reported upper limits of the DMDD-SI and -SD rates (as a function of m χ 0 1 ) thus aiding compliance of the scenario with these constraints.
The calculation for T c in BP-D1 suggests the possibility of a direct (type-I(1)) SFOEWPT (∆ SU (2) /T C = 1.14) from the trivial false minimum at {h d , h u , h s } ≡ {0, 0, 0} to the broken, true (global) minimum at {h d , h u , h s } ≡ {25.5, 145.6, −474.4} GeV at T c = 129.5 GeV. Note, however, that this is the only benchmark point that we present for which the phase with the true minimum does not nucleate successfully and hence the system would remain trapped at the metastable false minimum ({0, 0, 0}). This could render much of the parameter space (that otherwise favors SFOEWPT) cosmologically nonviable [45]. Nonetheless, we retain this point as a benchmark to demonstrate some characteristic collider-aspects, as discussed below, which could be equally instrumental in a scenario that does not have this shortcoming.
A CheckMATE analysis rules out BP-D1 (with r=1.12) via a CMS analysis [170] of 35.6 fb −1 worth data in the 3 + / E T final state where an opposite-sign, same-flavor (OS-SF) lepton (e or µ)-pair originates in the decay of an on-shell Z-boson coming from the decay of a heavier neutralino. Such a scenario, with µ eff as small as 275 GeV, would anyway be excluded more convincingly (i.e., with a larger r-value) by the recent LHC analyses in references [121,184] for the same final state which exploit 139 fb −1 of data.
The benchmark point BP-D2 is somewhat similar to that in BP-D1 in terms of the phenomenological features that are relevant for our discussion, i.e, the LSP is still singlinodominated with a very similar mass as before (≈ 60 GeV), the higgsino-like states are again the next heavier excitations with masses not very different (though on a little higher side) from those in BP-D1. Like BP-D1, BP-D2 also possesses relatively light singlet-like scalars. The DM phenomenologies, in both qualitative and quantitative terms, are rather similar in these two cases.
However, in contrast to that in BP-D1, in BP-D2, it is a two-step phase transition (type-II-S(1)-D(2)) as is suggested by the calculations of T c . First, a broken phase ({0, 0, 539.9} GeV) appears only along the singlet direction at T c = 151.5 GeV with a possibility of a first-order phase transition. This is followed by the appearance of another configuration at T c = 112.7 GeV for which SU (2) is now broken in the true minimum. This triggers the possibility of a second-order phase transition in which the scalar field could move from the evolved false minimum to the said true minimum. On the other hand, the calculation for T n now suggests that the tunneling process corresponding to T c = 151.5 GeV is so slow that what takes place instead is a strong (∆ SU (2) /T n = 2.2), one-step first-order transition along all three directions simultaneously (type-I(1)) from the trivial to the physical phase ({67.0, 197.8, 774.8} GeV) at T n = 96.2 GeV.
On the collider front, BP-D2 yields somewhat smaller production cross sections for the higgsino-like states than what BP-D1 gives because these states are a little heavier in BP-D1. A CheckMATE analysis indicates that in the CMS analysis in reference [170] that uses 35.9 fb −1 of data, the same final state (3 + / E T ) with an identical signal region as in BP-D1 becomes the most sensitive of the searches while, this time, barely disallowing (r = 1.01) the parameter point. A subsequent SModelS study indicates that the analyses in references [179][180][181], all involving 139 fb −1 of data, are even less sensitive. As for BP-D1, BP-D2 is also likely to be ruled out convincingly by the analyses of 3 + / E T final state with 139 fb −1 of data presented in references [121,184]. The resulting r-values would hint at how big a µ eff could thus be excluded in such a setup.
The point BP-D3 contains a somewhat heavier (∼ 91 GeV) singlino-like neutralino state where SFOEWPT (∆ SU (2) /T n = 3.8) remains viable, 'κ' being still small (∼ 0.071). The relic for such a singlino as a DM candidate is bound to be over-abundant in the absence of a suitable funnel, as is the case with BP-D3. The possibility of an efficient coannihilation, say with a bino-like state, requires a small mass-split between them which then tends to make the DD rates way too large to be acceptable. Instead, a bino-like LSP having a smaller mass and    possessing a suitable annihilation funnel via Z-or h SM could qualify as a DM. This is what we find in BP-D3 (with an h SM funnel, with |m χ 0 1 | ≈ 60 GeV). Spectrum-wise, this constitutes its basic difference from BP-D2. The pattern of phase-transition in BP-D3, as obtained from the critical temperature calculation, is also of a two-step kind (type-II-S(1)-D(1)) but is slightly different from what occurs in BP-D2, as can be seen in table 3. Although the bubble nucleation calculation indicates that both benchmark points have one-step SFOPTs in all three directions (type-I(1)).
However, the hierarchy among the lighter neutralinos now triggers important effects with major implications for the searches of the electroweakinos at the LHC. The higgsino-like states predominantly decay to the singlino-like NLSP and the Z-boson and/or h SM thanks to an enhanced higgsino-singlino mixing for a value of 'λ' which is on the larger side (∼ 0.57) [117]. Subsequently, the NLSP neutralino would undergo dominant decays to off-shell Zboson or h SM . Such cascades result in strengthened multi-lepton (more than three leptons) final states which now become far more sensitive to the recent LHC analyses when compared to the trilepton final states. The reason behind this is a much suppressed SM background for the former [170]. Indeed, a dedicated CheckMATE analysis confirms this effect and rules out BP-D3 rather emphatically (r = 2.13) by getting sensitive to the right (dedicated for finals states with more than 3 leptons) signal region ("G05") of the CMS analysis in reference [170] which considers data worth 35.9 fb −1 only. This needs to be contrasted with the verdict on higgsinos of very similar masses in BP-D2 in which those masses appear to be barely disallowed (r = 1.01) with the same set of data. Further, given the heightened sensitivity of the analysis to the multi-lepton finals states, it could eventually rule out even heavier higgsino-like states in such a setup.
The exercise undertaken in this section indicates how different types of spectrum for the light higgsino-like electroweakinos (i.e., smaller µ eff ), which otherwise comply with all relevant bounds including those from the DM sector and which allow for SFOEWPT, get ruled out by the LHC analyses with ∼ 36 fb −1 of data even when the latter's sensitivities to the targeted final states deteriorate significantly. The benchmark scenarios are so chosen that we end up with r 1. Such a value nominally reflects how light such electroweakinos could get before they start attracting bounds from the LHC analyses. Of course, more recent LHC analyses [121,184] with 139 fb −1 of data are expected to push these mass-bounds (and hence µ eff ) upwards but these are yet to be implemented in a recast package.

Allowed benchmark scenarios with successful nucleation
In this section we present a few benchmark scenarios that have all the good qualities of those listed in table 3 but now also pass the lower bounds on the electroweakino masses coming from some of the recent LHC analyses. Naively, this pushes µ eff up which impedes an efficient SFOEWPT with successful nucleation. The SFOEWPT now tends to proceed in two steps the details of which are presented in table 5 for our benchmark points. This is somewhat typical when the trivial and the global minima have a large separation between them in the field space [45]. This is since a larger µ eff corresponds to a larger v S at zero temperature for a given λ, a feature that governs the field-separation at T c .
The benchmark points in table 4 are picked up keeping in mind the following issues. While our goal is to find compatible points with smaller values of µ eff , we like to see the resulting scenarios have LSPs with different dominant admixtures. Allowing for this has a considerable bearing on both the DM phenomenology and searches for the electroweakinos at the LHC. Furthermore, these benchmarks possess a light singlet-like scalar below 100 GeV. This is since we set both 'κ' and A κ small which is preferred by SFOEWPT. Note that such a light singlet state inevitably affects both DM and collider phenomenologies, more so since the nature of the lighter electroweakinos, including the LSP, could get altered, simultaneously. As has been noted in section 3.2.2, to facilitate SFOEWPT we look for relatively light doublet-like Higgs bosons (by choosing A λ suitably) which are still allowed by the LHC Higgs searches.
In BP-A1 we have a higgsino-like lightest triplet with masses in the range ∼ 400 − 430 GeV with µ eff ∼ 422 GeV. Thus, the LSP and the NLSP are both higgsino-like (with their higgsino contents at 70% and 98%, respectively) while the lighter chargino is a nearly pure higgsino. As far as the DM sector is concerned, such a higgsino-like LSP DM is naturally under-abundant (Ωh 2 = 3.78 × 10 −4 ). This, in turn, generically helps satisfy the DMDD-SI and -SD constraints via downward scaling of the respective cross-sections.
As for the pattern of EWPT in BP-A1, the calculation for T c suggests that this is a two-step process of the type II-S(1)-D(1) as indicated in table 5 where first, at T c = 946.7 GeV, a broken phase appears in the singlet-direction followed by another in the SU (2) field directions at T c = 91.1 GeV. Subsequently, successful nucleations take place closely below the respective T c 's, down at T n = 946.6 GeV and 90.2 GeV. Note that in this particular case, calculations for both T c and T n suggest that the first phase transition (in the singlet-only direction) is just of a first-order kind while the second one, in the all-important SU (2) field directions that breaks the electroweak symmtery, is of a 'strong' first-order type (γ EW = 1.1) which is a crucial requirement for EWBG.
Searches for the lighter electroweakinos in BP-A1 effectively amounts to those for the higgsinos only where these states appear as the lightest triplet of electroweakinos which includes the LSP. This is since the heavier neutralinos, χ 0 3 and χ 0 4 , are singlino-and bino-like, respectively, whose productions are coupling-suppressed. In contrast to scenarios in which the higgsinos do not form the lightest triplet, here one loses out on the cascade of one of the neutralinos (which is the LSP in the present case). This restricts their abilities to contribute to diverse final states. On top of that, χ ± 1 and χ 0 2 , once produced in such a scenario, decays to the LSP, which is not far away in mass, via off-shell gauge and Higgs bosons (as indicated in table 4) thus resulting in associated leptons/jets to be generically soft. Both these issues have negative impacts on the experimental sensitivities of such a scenario. This is clearly reflected in the LHC analyses of such scenarios [121,122] which report much relaxed lower bounds (down to ∼ 220 GeV, conservatively) on the masses of such higgsino-like electroweakinos as a function of their mass-split with the LSP.
A CheckMATE analysis that includes all readily available analyses in its repository results in a 'r' value far below 1 for the point BP-A1 thus marking its total insensitivity to the LHC searches and hence allowed by the same. The relevant analysis and the most significant signal region therein are also indicated. Note that the higgsino masses for this benchmark point are     BM No.  Table 5. Phase transition characteristics of the benchmark points presented in table 4. For each benchmark point, presented are the T c 's and T n 's, the corresponding field values, the transition types ('FO' for first-order and 'SO' for second-order) and the strengths of the phase transition along the SU (2)-direction (γ EW ). See text for details. way above their current lower bounds for such a scenario as mentioned above. In passing, we note that in the future runs of the LHC such a scenario would likely attract bounds from the searches of the doublet-like heavy Higgs bosons sooner than from the direct searches for such electroweakinos.
Benchmark point BP-A2 is almost the same as BP-A1 except for M 1 now being brought down below µ eff . Thus, the LSP DM is now highly bino-dominated and its relic abundance (Ωh 2 = 0.107) now falls within the Planck-observed band. Towards this, the required depletion in the relic is again facilitated by the coannihilation of the bino-like DM with the higgsino-like chargino and neutralinos. Note that the sign on M 1 (with respect to that of µ eff ) ensures compliance with the experimentally observed latest upper bound on the DMDD-SI cross-section by setting up a so-called 'coupling blind spot' as discussed in section 3.2.2.
As in BP-A1, EWPT in BP-A2 is also of the type II-S(1)-D(1). The only notable difference that is found with respect to BP-A1 is in the delayed nucleation for the crucial phase transition in the SU (2) field directions (T n = 86.2 GeV, as opposed to 90.2 GeV in BP-A1) as shown in table 5. This results in a stronger FOEWPT (γ EW = 1.5, compared to γ EW = 1.1 in BP-A1). Its implications for the GW physics will be discussed in section 4.3.3. The delayed nucleation can be explained by the altered M 1 which modifies the thermal correction to the effective potential via terms that are only quadratic and quartic in m(φ)/T given that the bino is a fermion (see equation 3.5b). For our present benchmark scenario, successful nucleation would then require some appropriate modification in the term cubic in m(φ)/T which we achieve by a minor tweaking of A κ . In the process, the potential barrier gets modified in a way that leads to delayed nucleation compared to BP-A1.
On the LHC front, unlike in BP-A1, in BP-A2 cascades of both higgsino-like neutralinos (χ 0 2,3 ) will be important for the relevant final states. Although, just as in BP-A1, χ 0 2,3 and the higgsino-like χ ± 1 would undergo off-shell decays to Z, h SM and W ± , the corresponding branching fractions for χ 0 2,3 get suppressed in the presence of their significant on-shell branchings to a photon (for χ 0 2 ) and to a light a S (for χ 0 3 ). The relevant lower bound from the LHC on the masses of lighter electroweakinos with such mass-splits is presented in reference [121] for a wino(NLSP)-bino(LSP) system which can be conservatively taken as ∼ 300 GeV. In a scenario like BP-A2, such a bound would get weakened not only because of the suppressed off-shell branching fractions of the neutralinos as mentioned above but also, as described in section 3.2.2, since the collective production cross-sections for the higgsino-like electroweakinos are known to be smaller than if they were wino-like, for any given mass. This is corroborated by our CheckMATE analysis which indeed allows BP-A2. Compressed scenarios like BP-A1 and BP-A2 would, however, be sensitive to the HL-LHC. Also, as for BP-A1, BP-A2 is likely to be probed first in the searches for doublet-like heavy Higgs bosons at future LHC runs.
The benchmark point BP-A3, to start with, differs from BP-A2 in having a light singlinolike (NLSP, χ 0 2 ) state in-between the bino-like LSP (χ 0 1 ) and the higgsino-like chargino (χ ± 1 ) and neutralinos χ 0 3,4 . This is achieved by lowering the ratio κ/λ. The split between χ 0 2 and χ 0 1 is tailored to be rather small (∼ 5 GeV). Expectedly, the abundance of the highly binodominated LSP DM depletes via its coannihilation with the singlino-dominated NLSP. The DM relic abundance is found to lie within the Planck-observed band. Note that the proximity in the masses of these two states could, apriori, infuse a significant singlino component within the LSP thus pushing up the DMDD-SI cross sections dangerously. For the current benchmark scenario, such contamination has been tamed by requiring a relative sign between M 1 and m S (= 2κµ eff /λ) [117]. Achieving the coveted relative sign between these two quantities through a relative sign between M 1 and µ eff has an additional advantage since, as in BP-A2, this further helps restrict the DMDD-SI cross-section below its experimentally observed upper limit. The pattern of EWPT in BP-A3 is pretty similar to those in BP-A1 and BP-A2, i.e., this is a two-step process of type II-S(1)-D(1). However, the first transition along the singlet direction occurs somewhat later in time at around 644 GeV (in place of 945 GeV, as in BP-A2).
On the collider front, the higgsino-like χ 0 3,4 (χ ± 1 ) preferentially decay to singlino-like NLSP, χ 0 2 (thanks to their enhanced coupling given 'λ' is reasonably large [117]) along with an on-shell Z-boson and Higgs bosons (W ± boson). In turn, it is found that χ 0 2 dominantly decays to χ 0 1 γ (∼ 98%) as its decays to off-shell Z-and Higgs bosons are much suppressed due to a small mass-split between χ 0 2 and the LSP. Conservatively, when such photons go undetected due to their softness, cascades of the higgsino-like states via χ 0 2 would be effectively equivalent to their direct decays to LSP thus resulting in canonically sensitive final states like 3 + / E T and 1 + 2b-jets + / E T . Hence the reported bounds on the masses of the winolike electroweakinos from such final states, after correcting (relaxing) for the higgsino-like ones, would hold straightaway. Our CheckMATE analysis shows that BP-A3 survives this bound and is expected to be probed at the HL-LHC via the above standard searches for the electroweakinos as well as in the hunt for doublet-like heavier Higgs bosons.
The benchmark point BP-A4 differs from BP-A3 in the flipping of the nature of the LSP and NLSP, i.e., the LSP (NLSP) becomes singlino-dominated (bino-dominated). Furthermore, this is the only benchmark point where we find the CP -even singlet Higgs boson to be the lightest of the scalars (h S ∼ 74 GeV). Also, BP-A4 contains the smallest |µ eff | (∼ 335 GeV) among all four benchmark points presented in this table. The DM is found to be underabundant in the presence of multiple funnels (a S and h SM ). Hence the DMDD bounds are again satisfied thanks to the downward scaling of the DD cross sections.
The EWPT still takes place in two steps but is of type II-S(1)-D(2) as is suggested by the calculations of T c . The first of these is of the strong first-order type occurring along the singlet-direction at T c = 165 GeV. The subsequent transition occurs along the SU (2) direction at T c = 136.5 GeV and second-order in nature. However, the nucleation calculation indicates that the tunneling rate corresponding to the first transition is too small. Consequently, the actual nucleation from the trivial phase to the physical phase takes place directly (type I-(1)) at a later time at T n = 116.9 GeV. The possibility of such kind of a phase transition has already been pointed out in reference [45].
As in BP-A3, heavier higgsinos, χ 0 3,4 , decay to a bino-dominated NLSP (χ 0 2 ) and a singlino-dominate LSP (χ 0 1 ) accompanied by an on-shell gauge or a Higgs boson. Here also, χ 0 3,4 's decays to the singlino-dominated state (the LSP in this case) are favored as 'λ' is on the larger side. On the other hand, the bino-dominated χ 0 2 now undergoes a dominant decay to the singlet-like CP -even Higgs boson h S and the singlino-dominated χ 0 1 . This can play a crucial role in relaxing the relevant collider bounds [117,120] when the heavier higgsinos first decay to χ 0 2 which all have branching fractions around 20% for the present benchmark point (see table 4).
Among the implemented analyses in the CheckMATE package an older CMS one (with 35.9 fb −1 of data) [170] and another from the ATLAS (with 139 fb −1 of data) [183] show maximal sensitivities in the final states with 3 + / E T and 1 + 2γ + / E T , respectively. The corresponding 'r' values are found to be 0.55 and 0.53 which signify that the BP-A4 is still allowed by a wide margin by the electroweakino searches at the LHC. This may not be unexpected given that the heavier higgsinos do not always undergo one-step decays to the LSP as is assumed by the experimental collaborations. 184], are expected to have heightened sensitivities to the present benchmark scenario but are yet to be implemented in the recast packages. However, we have managed to check the constraints from the ATLAS analysis [121] for this benchmark scenario and we find [196] BP-A4 to be still allowed. The scenario is expected to get probed at the future LHC runs in the electroweakino searches first rather than in the searches for the heavy Higgs bosons. This is since these doublet-like heavy Higgs bosons are heavier in the present case (∼ 1.3 TeV).
As we have just discussed, benchmark points BP-A1, BP-A2, and BP-A3 have similar phase transition patterns while BP-A4 has one of a different kind. Hence, in figure 2, we choose to show the relevant phase diagrams for only BP-A2 (having an SFOEWPT in two steps with a palpable split between T c and T n ) and BP-A4 (one-step phase transition with a reasonably large T c and T n for the SFO) which may serve as the representative scenarios for BP   the purpose.

Prospects of GW detection
Values of various key parameters (T n , α, β/H n ) pertaining to the GW spectra arising from the FOPTs for the benchmark scenarios BP-A1 to BP-A4 are shown in table 6. The corresponding GW (frequency) spectra are calculated using equations 3.8-3.19 and are shown in figure 3.
These are further compared with the sensitivity of some space-and ground-based gravitational wave detectors, viz., LISA [197], Taiji [198], TianQin [199], aLigo+ [200], Big Bang Observer (BBO) [201] and Ultimate(U)-DECIGO [202]. Note that for each of the benchmark points BP-A1, BP-A2 and BP-A3 we observe a two-step first order phase transition. The phase transitions along the singlet-direction for these benchmark points happen at relatively larger temperatures (happens to be at larger β/H n and for smaller α). Thus, the contributions of the first FOPT (along the singlet direction) to the GW spectrum are relatively much smaller compared to the ones from the second FOPT which occur along the SU (2) field directions. This is why a spectral peak due to the first FOPT along the singlet direction for these benchmark points does not appear in figure 3. 6 In all four plots, the individual contributions from sound waves (within the bag model) and turbulence are shown with broken lines, in green and red colors, respectively. The total GW spectra, within the bag model, are denoted by solid, black lines. On the other hand, the same beyond the bag model (taking into account equations 3.15-3.17) are indicated in these plots by broken black lines. As can be found from these plots, the peak region of the GW spectrum obtained in the bag model for BP-A1 and BP-A3 lie only within the sensitivity of U-DECIGO. However, beyond the bag model, calculations suggest that the peak GW intensities for these two points fall short of the sensitivity of U-DECIGO. As for BP-A2, the peak of GW spectrum in the bag model lies within the sensitivity of ALIA, BBO and U-DECIGO while beyond the bag model it only falls within the sensitivity of U-DECIGO. Also, note that the peaks of these spectra are higher for BP-A2 when compared to BP-A1. Given that these two benchmarks have otherwise very similar types of phase transitions, the reason behind this can be traced back to the fact that in BP-A2 there occurs a stronger FOPT in the SU (2) field directions due to a relatively late-time nucleation and the associated values taken by 'α' and β/H n are such that the GW peak intensities shoot up thus increasing the prospects of observing the same. In contrast, for BP-A4, an extended section of the GW spectrum, around its peak, in the bag model, falls within the reaches of multiple experiments like LISA, Taiji, ALIA, BBO and U-DECIGO and, beyond the bag model, the same falls within the sensitivities of ALIA, BBO and U-DECIGO.
The quantity signal-to-noise ratio (SNR) is used to measure the detectability of the GW signal at the experiments. SNR is defined as [123] where T is the duration of the experimental mission in years, δ stands for the number of independent channels employed by an experiment to exploit cross-correlations (required to pin down the stochastic origin of the GW) and Ω exp (f ) h 2 denotes the effective power spectral density of strain noise of the experiment. Here, we consider δ = 2 for BBO and U-DECIGO and δ = 1 for LISA while for all of them we take T = 5. The SNR values for the benchmark points are found to be way below 1, with the exception of BP-A4, for which it is comparatively large but still < 1. These are to be compared with the reference minimum threshold value of SNR for the detection of GW which is taken to be 10 [123]. Thus, none of the benchmark scenarios meets this detectability criterion. It is, however, expected that there is a region of parameter space in the NMSSM which might give rise to stronger FOPTs (larger α) and the corresponding GW spectra lie deeper within the regions of experimental sensitivity (depending upon β/H n ) thus yielding an SNR value larger than 10. 7

Summary and outlook
Inspired by the prospects the Z 3 -NMSSM scenario holds in explaining the baryon asymmetry of the Universe via EWBG, which in turn requires SFOEWPT, we have sought to figure out how accommodating the scenario appears in the face of recent LHC results, in particular, the ones pertaining to the searches of the lighter electroweakinos which might happen to be higgsino-like and are favored by SFOEWPT. Various pertinent theoretical requirements, constraints on the Higgs sector (including the observed properties of the SM-like Higgs boson) from the LHC, flavor-constraints and bounds on various DM observables obtained from a host of dedicated experiments do already play their parts in delineating the allowed region of the NMSSM parameter space that still remains compatible with SFOEWPT. We further look into the prospects of detecting the (stochastic) GW arising from an SFOEWPT at various future experiments. The backdrop of our present study has been the looming tension between the physics of SFOEWPT and the recent LHC results from the electroweakino searches. While the former prefers µ eff in the range of a few hundreds of a GeV, the latter are tending to push µ eff steadily above such a ballpark. The general goal of such a study could then be to check if there is a meeting ground somewhere in the middle where both constraints are simultaneously complied with. A further observation is that EWPT is somewhat stubborn in its need for relatively small µ eff . A middle ground can thus only be found if the reported constraints from the LHC could be evaded under circumstances that have not been considered explicitly 7 The exploration, however, demands a more detailed study which we will take up in a future work.
by the LHC experiments. In this work, we exploit such caveats to our advantage via recasts of the relevant LHC analyses using popular packages like CheckMATE and SModelS.
Thus, the region of the Z 3 -NMSSM parameter space that concerns us in this work is characterized by reasonably small µ eff that yields relatively light higgsino-like states. Also, SFOEWPT prefers a relatively light CP -even singlet-like scalar, h S , thus requiring 'κ' to be small. This leads to a relatively light singlino in the spectrum which can be the DM particle while the light singlet scalars play crucial roles in the DM phenomenology. Also, for our benchmark scenarios, we choose relatively large values of 'λ' ( 0.5) and smaller values of tan β ( 5). Together, these yields m h SM in the right ballpark (∼ 125 GeV) without requiring too large SUSY radiative corrections. Furthermore, a smaller tan β could allow the heavier doublet Higgs bosons to remain relatively light (reminiscent of the "alignment without decoupling" scenario) which might aid SFOEWPT. With 'λ' on the larger side, the mixing among the higgsinos and the singlino, all of which can be relatively light thus being the candidates for the LSP DM and the NLSP, can be sizable. Hence, with such choices of theory parameters, the physics of the EWPT (and hence EWBG) becomes intricately connected to the DM and collider (LHC) phenomenologies.
Results of a scan over the parameter space are presented depicting first how relatively small µ eff ( 500 GeV) fairs against various theoretical and basic experimental bounds including those from the observed Higgs sector and the ones from the DM experiments. Two sets of benchmark scenarios are then presented to demonstrate the SFOEWPT-DM-LHC connection. These scenarios are checked to give rise to SFOEWPT by using the package CosmoTransitions in which we implemented the framework of Z 3 -NMSSM, matched to THDMS, as has been a pretty standard practice for the purpose.
With one set of benchmark scenarios, we have sought to find out up to what a ballpark maximum value of not so large a µ eff can still be ruled out by recent LHC analyses, in particular, when one departs from the simplified assumptions on the decays and branching fractions of the cascading electroweakinos which is expected to relax the reported bounds on the electroweakino masses. Subjecting this set of otherwise highly motivated scenarios to thorough recasts of some pertinent LHC analyses (with both 36 fb −1 and 139 fb −1 of data) with the help of CheckMATE and SModelS reveals that µ eff 300 GeV, with low values of 'κ' ( 0.1) and larger 'λ' ( 0.5), is mostly ruled out. It should be noted that smaller values of µ eff already attract severe constraints from the DM direct detection experiments. Thus, in this regard, the LHC searches might not always yield a robust improvement over the DM bounds. An even smaller µ eff could, however, survive the LHC bounds for a compressed electroweakino spectrum. In this work, being conservative, we do not consider this possibility.
With the other set of benchmark scenarios, we have demonstrated how low a µ eff could still be allowed instead. A similar exercise shows that, under favorable circumstances, upwards of µ eff ∼ 335 GeV could survive the LHC onslaught. This is rather encouraging since we find that EWPT could still remain to be of strong, first-order type even for µ eff as large as ∼ 425 GeV which is the case for a couple of benchmark scenarios that we have presented. These also show that a viable LSP DM can be bino-or singlino-like or even a mixture of bino, singlino and higgsino states. We have thoroughly studied the properties of EWPT in these scenarios with the help of CosmoTransitions and have found that for µ eff on the larger side, a two-step phase transition is a more likely phenomenon with the first transition taking place in the singlet field direction followed by the other in the SU (2) field directions.
For these latter set of scenarios, we have thoroughly studied the stochastic GW (background) spectra that might carry the imprints of FOPT from new physics beyond the SM. We find that the signal intensities lie inside the sensitivity limits of one or more of the future/proposed experiments like the LISA, BBO, UDECIGO, Taiji, Alia, etc. However, the SNR values, as such, are not found to be healthy enough to guarantee a positive detection.
In summary, the present work corroborates the basic findings reported in the literature pertaining to SFOEWPT, in particular, and EWPT, in general, in the framework of the Z 3 -NMSSM. We broadly concur with various reported patterns and features of EWPT in such a scenario and the different conditions under which those manifest. We then go beyond to shed light on what the recent searches of the electroweakinos at the LHC have to say about the viability of SFOEWPT in the current framework while compatibility with the constraints from various pertinent theoretical and experimental sectors including the DM sector is ensured all through. Furthermore, it appears that the GW signals resulting from the strong FOPTs in these scenarios are likely to remain too weak to be detected at future dedicated experiments.
As for an outlook, new LHC studies with data from the recently terminated LHC Run 2 and those that would arrive soon from high luminosity LHC (HL-LHC) are likely to shed a more unambiguous light on the broad viability of EWBG within the Z 3 -NMSSM while these continue to explore electroweakinos at larger masses and in difficult scenarios like the compressed ones. Furthermore, improvements are possible in the theoretical calculations of several key EWBG objects, viz., the bubble wall profile, wall velocity and CP -violation and in the dealing of the transport equations which could lend a more accurate estimate of the relation between the NMSSM parameters and EWBG. With these, a reassessment of the detectability of such GW signals may be warranted which might prove the latter's role as complementary to the LHC searches. A synergy like this between LHC and GW physics is likely to be rather intriguing and we reserve such a study for a future work.

Acknowledgments
AC thanks Harish-Chandra Research Institute (HRI) for hosting him during the course of this collaborative work. AC also acknowledges partial support from the Department of Science and Technology, India, through INSPIRE faculty fellowship (grant no: IFA 15 PH-130, DST/INSPIRE/04/2015/000110). SR is supported by the funding available from the Department of Atomic Energy (DAE), Government of India for the Regional Centre for Acceleratorbased Particle Physics (RECAPP) at HRI. SR would like to thank Peter Athron, Sebastian Baum, Junji Cao, Ulrich Ellwanger, Andrew Fowlie, Tathagata Ghosh, Thomas Konstandin, Sabine Kraml, Indrani Pal, Avik Paul, Krzysztof Rolbiecki, Tim Stefaniak, Di Zhang, Yang Zhang for helpful discussions and communications. SR also acknowledges the use of cluster computing available at the High Performance Scientific Computing facility at HRI and thanks Rajiv Kumar for his technical help at this facility.

Appendices
A Matching the NMSSM parameters to those in the THDMS potential In terms of the CP -even Higgs fields h d , h u and s, the tree-level Z 3 -NMSSM potential of equation 2.3, which is relevant for the study of phase transitions, can be written as [57] V NMSSM tree (h d , h u , s) = 1 32 On the other hand, the tree-level Z 3 -symmetric THDMS potential is given by [43,[100][101][102] V THDMS All parameters in equation A.1 and A.2 are taken to be real as we do not consider any CP -violation in the Higgs sector in this work. We use NMSSMTools to obtain the particle spectrum of the Z 3 -NMSSM at the scale M SUSY . Except for the electroweakinos with masses around a few hundred GeV, we consider all other SUSY excitations to be much heavier such that those may be considered effectively decoupled from the physics of phase transitions. However, to avoid large logarithmic corrections from appearing in V NMSSM At one-loop, the only relevant threshold correction that arises (as the top squarks are integrated out at the scale M SUSY ) is to λ 2 and is given by [203][204][205][206] where A t is the soft-SUSY-breaking top squark-Higgs trilinear coupling in the scalar potential and y t is the top quark Yukawa coupling, both defined at the scale M SUSY . Note that all the NMSSM parameters are also provided at the scale M SUSY and in the DR scheme. Thus, after matching, all the THDMS parameters also get defined at the same scale and in the same renormalization scheme.

B RGEs in the THDMS
We borrow the set of relevant RGEs from reference [38] which we use to run the THDMS model parameters from the scale M SUSY to the scale m t . Contributions to the β-functions from the SM gauge bosons, the Higgs bosons, the top quark, the higgsinos and the singlino are included. As for the gauginos, the contribution from the bino is known to be small (even when M 1 does not get to be too large, as is the case in our present analysis) while M 2 is set at a rather large value. Hence we ignore their effects. With these, the one-loop β-functions for the model parameters q i are given by where 'Λ' is the energy scale. The one-loop RGEs for the quartic couplings λ i , i ∈ {1, 2, 3, ..., 8}, the mass parameters (m 4,5 ) and v s appearing in V THDMS

C Field-dependent masses and the daisy corrections
Here we present the field-dependent mass-squared matrices for the scalar sector which are derived from equation A.2 [43,100]. Note that in the tree-level field-dependent mass-squared matrices (see equations C.1, C.2 and C.3) the values of m 2 1 , m 2 2 and m 2 3 are without this modification since the latter are solutions of the corresponding tree-level tadpole equations as presented in equation C.4.
In the fermionic sector, we consider the top quark, the bottom quark, the tau lepton along with the four neutralinos (χ 0 1,2,3,4 ) and the one chargino (χ ± 1 ), since the wino-like states are taken to be much heavier and hence are decoupled from the physics of phase transitions. The field-dependent masses of the top quark, the bottom quark and the tau lepton are given by [57] m t = 1 The field-dependent 4 × 4, symmetric neutralino mass matrix in the basis { B, H 0 d , H 0 u , S} is given by On the other hand, the field-dependent mass of the higgsino-like chargino is given approximately by m χ ± 1 λs √ 2 . Note that the masses of these electroweakinos are given in terms of the NMSSM model parameters which are defined at the scale M SUSY . This is acceptable for our purpose since these masses do not appear in the tree-level potential. 9 The field-dependent masses of the gauge bosons, W ± and Z, are given by Note that in the daisy potential of equation 3.6, M 2 h and M 2 V are the eigenvalues of the thermally improved (i.e., Debye-corrected) mass-squared matrices for the Higgs and the gauge bosons, respectively, i.e., generically, M 2 = eigenvalues[M 2 + ∆(T 2 )] where ∆(T 2 ) = c ij T 2 and c ij 's are the so-called daisy coefficients. From the high temperature expansion of the thermal one-loop potential V T (of equation 3.3) using equation 3.5, the daisy coefficients can be found from the following relation: For an FOPT, the daisy correction is especially important since it has an impact on the all very critical cubic term of the potential at a finite temperature. Effects of only the scalars and the longitudinal modes of the vectors are included in this contribution. Thermal contributions to the transverse modes are suppressed due to gauge symmetry [207]. With these in mind, the various daisy coefficients (neglecting the electroweakino contributions) are as follows [43,90,208,209]: where the subscripts {1, 2, 3} refer to the fields {h d , h u , s}. Note that the gauge symmetries plus the discrete Z 3 symmetry of the model set the off-diagonal terms of the ∆(T 2 ) matrix to zero (i.e. ∆(T 2 ) is a diagonal matrix). The longitudinal components of the gauge bosons receive thermal corrections. For W ± bosons the correction is c W ± L = 2g 2 2 T 2 . Thus, the thermally improved mass of the longitudinally polarized W ± bosons is given by Similarly, the longitudinal components of the Z-boson and the photon (A) fields also receive thermal corrections. Their masses can be determined by diagonalizing the following matrix: (C. 12) The thermally improved masses of the longitudinally polarized Z-boson and the photon are given by (C.14)