Gravitational Waves, baryon asymmetry of the universe and electric dipole moment in the CP-violating NMSSM

In this work, we make the first study of electroweak baryogenesis (EWBG) based on the LHC data in the CP-violating next-to-minimal supersymmetric model (NMSSM) where a strongly first order electroweak phase transition (EWPT) is obtained in the general complex Higgs potential. With representative benchmark points which pass the current LEP and LHC constraints, we demonstrate the structure of EWPT for those points and how a strongly first order EWPT is obtained in the complex NMSSM where the resulting gravitational wave production properties are found to be within the reaches of future space-based interferometers like BBO and Ultimate-DECIGO. We further calculate the generated baryon asymmetries where the CP violating sources are (1): higgsino-singlino dominated, (2): higgsino-gaugino dominated or (3): from both sources. It is shown that all three representing scenarios could evade the strong constraints set by various electric dipole moments (EDM) searches where cancellations among the EDM contributions occur at the tree level (higgsino-singlino dominated) or loop level (higgsino-gaugino dominated). The 125 GeV SM like Higgs can be either the second lightest neutral Higgs H2 or the third lightest neutral Higgs H3. Finally, we comment on the future direct and indirect probe of CPV in the Higgs sector from the collider and EDM experiments.


Introduction
The standard model (SM) of particle physics provides a satisfactory description of particle physics phenomena over the past few decades together with the discovery of the Higgs boson [1,2], a key element to the electroweak symmetry breaking mechanism. However the theory is not yet perfect since the SM can not provide a dark matter candidate and furthermore can not explain the baryon asymmetry of the universe (BAU). The BAU is traditionally characterized by the baryon to entropy density ratio and has been measured to a high precision by Planck [3], Among the various baryogenesis mechanisms, the electroweak baryogenesis (EWBG) [4] is among the most theoretically well motivated and experimentally testable scenarios since it connects BAU generations to the details of the electroweak symmetry breaking (see [5] for a recent review). According to Sakharov [6], generation of a nonvanishing baryon asymmetry requires three ingredients in the particle physics of the early universe: baryon number violation, C and CP-violation (CPV) and out of equi-librium conditions. In the framework of of EWBG, the non-equilibrium environment for baryon generation is provided through a first order electroweak phase transition(EWPT). Even though the SM provides baryon number violation through the electroweak sphaleron process, it fails to provide a first order EWPT with a 125 GeV Higgs as well as a large enough CPV for sufficient baryon generation. Physics scenarios beyond the SM, with the capability of providing a first order EWPT and with new sources of CPV, are therefore resorted to for a successful baryon number generation. Among the various new physics scenarios, the supersymmetric theory is among one of the most popular and intensively studied models due to its theoretical attractions that it can solve the gauge hierarchy problem, has the dark matter candidate, leads to strong-electroweak gauge unification and also can potentially explain the origin of baryon asymmetry in the universe. The minimal supersymmetric standard model (MSSM) accommodates two Higgs doublets, where the lightest of neutralinos can serve as a dark matter candidate and it can also provide new sources of CPV needed for baryon asymmetry. However the MSSM suffers from the µ-problem [7-9] and it is also hard to obtain a strongly first order EWPT (SFOEWPT) in this model [10][11][12][13]. Both of these problems can be solved in a minimal way by extending the MSSM with an electroweak singlet chiral superfied. This is known in the literature as the next-to-minimal supersymmetric standard model (NMSSM) [14,15]. With this addition of the extra field content, the particle spectrum of NMSSM now accommodates two extra scalars giving then a total of five neutral scalars and a pair of charginos. This makes it quite easier to obtain a SFOEWPT in this model which has been studied in the CP conserving case [16][17][18][19][20][21][22][23][24][25]. Besides, there can be CPV at tree level from the new interactions in contrast to the case in MSSM where CPV occurs at loop level. This new CPV in NMSSM is manifested in terms of a tree level mixing between the set of three CP-even scalars and the set of two CP-odd scalars, similar to what happens in the two Higgs doublet model [26][27][28][29]. With the presence of these two kinds of CPV sources, there appears the chance in the parameter space that their CPV effects cancel [30] in contributions to the electric dipole moments (EDM) of electron, neutrons and atoms, etc,. and therefore can evade the stringent constraints from null search results of the EDMs. Moreover, even though the studies of BAU in the context of MSSM and in the NMSSM have been performed in the literatures [31,32], a joint analysis with these two different origins of CPV is still lacking. We therefore make an updated study of the EDM and BAU phenomenologies with these two different CPV sources taken into account. Furthermore, with the addition of the gauge singlets as well as the inclusion of the CPV in the Higgs potential, the effective potential is now of relatively high dimension and this can potentially leads to a rich structure of phase transition patterns. In particular, the phase transition does not necessary proceed through one step to the electroweak vacuum and generally need multiple steps [19] for this to happen which therefore deserves a detailed scrutiny. In addition, after the discovery of the gravitational waves from merging black holes detected by LIGO [33], there has been increasing interests in the discussions of the gravitational waves from the EWPT in various new physics models (see Ref. [34] for a recent review). Therefore it deserves a similar study in the NMSSM to augment the analysis of the baryogenesis and EDMs.
The rest of this paper is organized as follows. We first introduce our conventions for the complex NMSSM in Sec. 2. We explore the phase structures of this model in Sec. 3 where the impact of CPV on the EWPT is studied with representative benchmark parameter space points. We then present the properties of gravitational wave signals generated during the EWPT for these benchmark points in the following Sec. 4. We perform a joint analysis of the baryon asymmetry in the framework of EWBG and the constraints from null searches of EDM in Sec. 5 and briefly comment on the collider probes of CPV in Sec. 6 after which we give the summary in Sec. 7.

The CP-violating Higgs sector
We consider the Z 3 -invariant complex NMSSM in which the Higgs superpotential and the soft supersymmetry breaking terms are given respectively by [14] with here the fields with a hat being chiral superfields and those without a hat being the corresponding scalar components. Here the limit κ → 0 corresponds to the Pecci-Quinn limit. Collecting the F-terms from above superpotential W Higgs , the Higgs interactions from ∆L soft and the D-terms from the Higgs gauge interactions, one obtains the tree-level Higgs potential: and each of above contributions is given explicitly by,

093106-2
We take here the parameters λ, A λ , κ and A κ to be generically complex as opposed to the CP-conserving case where these parameters are real. Thus these complex parameters provide additional sources of CPV aside from those originally appearing in the MSSM case. After the electroweak symmetry breaking, two further phases θ 0 u and θ 0 S appear in the expansion of the Higgs fields about the vacuum expectation values(VEVs) of the neutral components of H u , H d and S which are defined with the following convention [19], Here v u and v d are the VEVs of H u and H d respectively, h u , h d , h S are the CP-even scalars, a u , a d , a S are the CPodd scalars and we have removed the phase of H d by the gauge transformations. Then, the effective higgsino mixing parameter in NMSSM is given by µ = |λ|v s e iφµ / √ 2 with φ µ =θ 0 s +φ λ . For translating into physical parameters, three of the minimization conditions about the VEVs allow us to replace the soft mass parameters m 2 Hu , m 2 H d and m 2 S by v u , v d and v S while the remaining three conditions yield relations among the complex parameters indicating that the CP phases are not all independent. Considering this, we follow Ref. [19] and introduce the following notations to make our discussions independent of the phase conventions, with here φ ′ λ ≡φ λ +θ 0 u +θ 0 S ,φ ′ κ ≡φ κ +3θ 0 S and φ ′ λ , φ A λ , φ ′ κ , φ Aκ are the phases of λ, A λ , κ, A κ respectively. Here I λ and I κ denote the CP phases of the first and second terms in the soft SUSY breaking Lagrangian and due to the minimization conditions, they are related to the single CPV source I by So there is only one physical CP phase at the tree level: κ +A κ can be determined from this up to a two-fold ambiguity. However above con-ditions face modificatioins by loop corrections and this will be discussed in the one-loop effective potential in the following.
With the presence of the Higgs sector CPV corresponding to φ ′ λ −φ ′ κ = 0, the neutral Higgs bosons need not carry any definite CP parities already at the tree level and their mixing can be described by an real orthogonal 5×5 matrix O as with here H 1(5) defined as the lightest (heaviest) Higgs mass eigenstate, a ≡ a d sinβ +a u cosβ and the massless Goldstone boson G 0 ≡ a u sinβ −a d cosβ decouples from above scalars in the quadratic terms. This mixing constitutes the most significant difference compared with the MSSM case, in which the CPV can only be induced at one-loop level, mostly through loop corrections to the Higgs self-energy, from the phase of the soft SUSY breaking trilinear couplings A u ,A d ,A e of the up-type, downtype and charged lepton-type sfermions, as well as the soft SUSY breaking mass parameters of the gauginos M 1 ,M 2 and M 3 . To capture the main features of EDMs and EWBG, we study the scenarios with the tree-level CPV phase coming from φ κ and the loop-level induced CPV phase being φ ′ µ ≡ φ µ −φ M 2 . To describe the EDMs properties clear, the three phases of φ κ,µ,M 2 will be used as input parameters in the numerical analysis of Sec. 5.

Electroweak phase transition
The phase structure of the model is governed by the finite temperature effective potential V eff , composed of the following contributions, Firstly, the tree level effective potential V Tree is given by Eq. (4), where ϕ d ,ϕ u ,ϕ S ,θ u and θ S are the background fields of the Higgs scalars. For convenience of the calculations of EWPT, we construct the 5-dimensional (5D) order parameters of EWPT [19], with here ∆θ u = θ u −θ 0 u and ∆θ S = θ S −θ 0 S . Using these parameters and with the definitions in Eq. (7), the the tree-level contribution can be rewritten as The second part in the effective potential is the oneloop Coleman-Weinberg potential V CW [35], given by where i runs through all particles in the model, with degrees of freedom n i , field-dependent mass m i ( ϕ) and spin s i . The above result is given in the Landau gauge with renormalization scale Q and the constant c i is 5/6 for gauge bosons and 3/2 for the others. With above Coleman-Weinberg potential included, the tadpole conditions are modified by these one-loop corrections. Hence, the counterterms V CT are introduced to preserve the tree-level relations for VEVs [36], i.e.
where it is sufficient to include the following terms in V CT : The last piece in the effective potential is the thermal corrections V T at the one-loop level [37], (17) Note these lowest order thermal corrections need to be improved by daisy resummations [38,39] at high tem-peratures when the perturbative expansion of the potential fails. This in practice can be done by insertions of thermal mass corrections, [36]. Here, we note that the Goldstone contributions may leads to a negative m 2 i which has been treated as in Ref. [40].
With the finite temperature effective potential given above, the phase structure of the model can be readily obtained by searching for its minima at each step as temperature drops. The basic picture of the phase evolution is as follows. At sufficiently high temperature, the universe sits at the origin of the 5D field space where the electroweak symmetry is manifest.
As temperature drops, other minima will develop in this space and the universe will transit to the lower minimum. Depending on the presence or absence of a barrier between these two adjacent minima, this phase transition is classified as first order or second order respectively. The phase transition occurs through nucleations of the bubble within which the lower minimum is developed [41][42][43]. The bubbles expand, collide and coalesce with each other leaving eventually the universe in the vacuum of the lower minimum. During this process, the temperature at which these two minima are degenerate is the critical temperature T C . For baryon asymmetry to be generated within the non-equilibrium phase boundaries, the phase transition that is connected to the electroweak symmetry breaking needs to be strongly first order to sufficiently quench the weak sphaleron process inside the bubble and thus to preserve the generated baryons. This generally translates into the widely adopted criteria φ EW (T C )/T C 1 [44][45][46] where φ EW ≡ ϕ 2 u +ϕ 2 d 1) . Due to the addition of the extra scalar singlet and the inclusion of CPV in the Higgs potential, the field space is of relatively high dimension and thus the phase transition history can be of quite rich structures [19] and typically have a two step EWPT [25,[48][49][50][51][52][53][54][55][56] 2) . Due to the relatively large parameter space of this model, we seek here only representative benchmark points to illustrate the phase structures while a more comprehensive analysis of the model parameter space is deferred to future works. In order to show the effect of CPV to the EWPT, here we present and compare two representative benchmark points corresponding to the CP-conserving and CP-violating cases respectively, 1) Both the φ EW (T C ) and T C involve gauge dependence issues [47]. One should keep in mind that we analyze the EWPT in Landau gauge when interpreting our results.
2) Notice that those models can be probed through the resonance di-Higgs searches [57][58][59].  Table 1 (top). The bottom plots are for the CP-violating case obtained as described in the caption of Table 1. For the right plots, also plotted are the critical temperature TC in brown and the nucleation temperature Tn in green.
with the numerical results here obtained using Cosmo-Transitions [60]. The resulting phase structures for the CP-conserving case is shown in the top panel of Fig. 1 corresponding to the benchmark point in Table 1.
The CP-violating case, which is obtained as described in the caption Table 1, is shown in the bottom panel.
For both cases, the right plots show the evolution of the electroweak background field amplitude φ EW (T ) as temperature drops from right to left. Here blue lines denotes the high temperature phases which are all observed to leave the origin continuously, thus signifying its second order nature. As temperature decreases, the phases corresponding to the red lines appear and denote the minima which eventually evolve to the electroweak minima. During the temperature ranges where these two phases overlap, the critical temperatures are shown with brown dashed vertical lines and is 157 GeV for the CPconserving case while a slightly larger value of 158 GeV is found when the CPV is turned on. So the inclusion of CPV tends to drive the onset of EWPT earlier while the duration of the high phase is seen to be shortened. Also shown in these plots are the nucleation temperature to be defined in the following section. For both cases, the left plots show the evolution of the difference of V eff at the high and low phases during the overlapped temperature ranges, where the critical temperatures correspond to ∆V = 0. As was observed in these plots, the EWPT for these two benchmark points proceeds through two steps with the first step being second order and the subsequent one being first order with its strength being φ EW (T C )/T C =1.039, consistent with the SFOEWPT criterion 1) .
1) For the second step EWPT, the weak sphaleron process outside the electroweak bubble is suppressed compared with that in the symmetric phase since electroweak symmetry is already broken outside. A benchmark point with a better EWPT pattern could generally be found from a dedicated scan over the NMSSM parameter space which however is beyond the scope of this work.

Gravitational waves
During the first order EWPT, there can be gravitational waves generated, mainly coming from three processes: bubble collisions, sound waves in the plasma and Magnetohydrodynamic turbulence [61] (see [34] for a recent review). The total energy density of the resulting gravitational waves is approximately a sum of these three contributions, where we adhere to the Hubble constant definition H = 100hkm s −1 Mpc −1 . The energy spectrums depends on the electroweak bubble profiles during the EWPT which corresponds to the bounce solutions minimizing the action, leading to the equation of motion for solving the bubble profile φ b , In addition, for a bubble to be effectively growing, triumphing the force from surface tension, the following condition needs to be met [62] which serves as the definition of the nucleation temperature T n . This then translates into finding the nucleation temperature T n such that S 3 (T * )/T * | T * =Tn ≈ 140 [62] is satisfied. From the bubble profiles and the nucleation temperature, the following two key parameters α and β, relevant for gravitational wave calculations, can be obtained α= 30∆ρ where in the evaluation of α, ∆ρ is the difference of energy density between the false and true vacua [22], and H n is the Hubble constant evaluated at the nucleation temperature T n . With these parameters solved, one can obtain the energy spectrum of the gravitational waves.
Firstly for the gravitational waves from the bubble collision, its contribution can be estimated using the envelop approximation [63][64][65] and results in the following spectrum [66], with here v w being the bubble wall velocity and κ characterizing the fraction of latent heat deposited in a thin shell, both of which are functions of the previously defined parameter α [67], Moreover f env is the peak frequency at present time and is approximately given by Hz.
Secondly, for the contribution from sound waves, it is given by Here κ v is the fraction of latent heat transformed into the bulk motion of the fluid and is approximately given by κ v ≈α(0.73+0.083 √ α+α) −1 in the case of v w ≈1 while for the peak frequency f sw , we have Hz. (27) Finally the MHD turbulence contributes with here κ turb ≈ 0.1κ v and in this case, the peak fre- 10 -4 f(Hz)  quency f turb is Hz. (29) After the discovery of the SM-like Higgs at LHC, the gravitational waves produced by bubble collision in the CP-conserving NMSSM has been studied in Ref. [22]. The gravitational waves produced by the bubble collision and sound waves has been studied later [25]. We explore the prospects of GW signals for the CP-conserving and CP-violating benchmarks presented in Table 1. Specifically we use the package CosmoTransitions [60] to solve the bounce solutions and to find the nucleation temperatures T n with which we calcualte the parameters α, β and finally the gravitational wave spectrums. In this process, we observe that the benchmarks both satisfy the nucleation criterior S 3 (T )/T <140 below certain nucleation temperatures T n . The calculated nucleation temperatures for these two cases are shown in the respective plots in Fig. 1 with green dashed vertical lines. We further show the resulting gravitational wave spectrums for these two cases: Fig. 2 for the CP-conserving case and Fig. 3 for the CP-violating case. In the left plot of Fig. 2, we also show the profile of S 3 (T )/T in the neighborhood of its nucleation temperature T n =134.6 GeV for illustration. In these gravitational wave spectrum plots, the three individual contributions from bubble collision, sound waves and turbulence are ploted using blue dot-ted, brown dotdashed and yellow dashed lines with their sum denoted by a solid cyan line. In each of these cases, the sound wave contribution dominates and is almost indistinguishable with the cyan line. The color shaded regions in these plots are the experimental sensitivity  Table 1 and the bottom plots in Fig. 1. The plotting conventions are the same as that in the right plot of Fig. 2.

093106-7
regions for several proposed space based interferometers: LISA [68] with two design configurations in notation NiAjMkLl [61,69], BBO, DECIGO (Ultimate-DECIGO) [70] and ALIA [71] For the CP-conserving case corresponding to Fig. 2, the gravitational wave energy spectrum falls within the sensitivies of BBO and Ultimate-DECIGO while unreachable by the others. Turning on CPV as plotted in Fig. 3, the magnitude of the energy spectrum decreases a little bit and can barely be detected by BBO. We note again that, these benchmark points presented here only constitute a tiny fraction of the whole parameter space of the NMSSM and therefore there might be other parameter space points which can give significantly enhanced GW spectrum. Exploration of this vast parameter space however requires a significantly improved calculations of the bounce solutions, a task which we will defer to future works.

Anatomy of BAU and EDM
In this section, we study the explanation of BAU in the framework of EWBG under the constraints of EDMs measurements. Typically, we focus on the CPV phases coming from tree level and loop level together.

The Electroweak Baryognesis in the CPV NMSSM
With the tree level CP-phases taken into account, there are in general three different CPV source terms in the quantum transport equations governing the dynamics of the particle densities in the framework of the closed-time-path-formalism (CTP) (see Ref. [72] for pedagogical discussions). Under assumptions of chemical equilibrium of Yukawa interactions, strong sphaleron and gaugino interactions, the set of coupled transport equations can be reduced to a single equation of H [72][73][74], where a one dimensional picture of the expanding bubble wall assumption is chosen as usually adopted in the literature to simplify the calculations, v w is the bubble wall velocity andz is the spatial coordinate perpendicular to the wall in the frame where the wall is at rest with negative values ofz corresponding to the symmetric electroweak phase (bubble exterior). It should be note that the wall velocity for the gravitional wave and baryogenesis studies should be different due to the the fist one require a relatively larger v w in comparison with the latter one. The v w for the calculation of gravational wave and baryogenesis can be different when taking into account the hydrodynamics of the bubble [75,76], the consistent calculation can be performed with the one for gravational wave being the speed of the wall and the one for the calculation of baryogenesis being the relative ve-locity between the wall and the plasma.
In the above equation, the quantitiesD,Γ andS are the effective diffusion constant, effective relaxation rate and effective source term for H(z) respectively and the explicit formulae are [73,74], Here the Higgs rate Γ h is and the relaxation rates for the stop, sbottom, and stau cases could be computed following Ref. [72-74, 77, 78]. The contributions from the Higgs sector CPV are summarized in the quantity S CP H given by Here the first term S CP H ± is the gaugino-higgsino driven CP-violating source term and can be computed in the vev insertion approximation from a tree-level wino W ± mediation. The result is given by [72] S CP H ± (x) = whereβ(x) = dβ/dt and can be approximated byβ ≈ v w ∆β/L w with here v w being the bubble wall velocity, L w the wall width and ∆β the difference of the value β(x) outside and inside the bubble. So above source term grows linearly with ∆β and our choice is optimal. This dependence of ∆β is to be regarded as a theoretical uncertainty and the precise determination of ∆β in the NMSSM is however beyond the scope of this work. For the more precise calculation, the profiles of the bubble wall and the profile-dependent masses and widths should be taken into account when solving the quantum transport equations during the EWPT. The second term S CP H 0 (x) in Eq. (33) denotes the contribution from the neutral higgsino mixing terms and the corresponding result for W and B intermediate states can be obtained from Eq. (34) by the replacements: Since these above two source terms depend on sinφ ′ µ and therefore they vanish when the phase φ ′ µ is set to zero. For further notations about the quantities appearing here, we refer the readers to [72]. We note that in the NMSSM, these gaugino-higgsino sources driven EWBG have been studied in [32], wherein tree level CP phases are all imposed to be zero [73]. The third term in Eq. (33) is the higgsino-singlino driven CPV source and has been studied in Ref. [31], including the singlino thermal mass term. Here we have assumed that there is no spontaneous CPV and the functional form of the Fermionic source function I f SH 0 can be found in Ref. [31]. It is obvious from Eq. (36) that S CP S H 0 vanishes when the phase combination vanishes, that is, when sin(φ λ −φ κ )=0.
Aside from the Higgs sector CP-violating source terms in Eq. (33), we also included the source terms from sfermion sector: Since the source terms are negligibly small forz <−L w /2, and that the relaxation terms have the approximate form Γ (z)= θ(z)Γ , the solution of H in the symmetric phase accepts the following solution from Eq. (30), with the prefactor given by where the usually encountered factor γ ± is, With the profile for H solved under previous assumptions, the profiles for the other particle densities can be obtained. In particular the left-handed Fermionic charge density, which serves as the source for the weak sphaleron process for generating the baryons, can be obtained, The left handed density gets converted into a baryon density n B through weak sphaleron transitions. The obtained baryon number density in the electroweak broken phase is a constant given by [79], Γws vw z .
In above formulation, we have used the approximation that due to the much smaller weak sphaleron rate Γ ws and thus the rates for both the creation of n left and its diffusion ahead of the bubble wall, above Eq. (43) is usually decoupled from the set of diffusion equations.

CPV Phases v.s. EDM in CPV NMSSM
In the MSSM, the CPV could only occur at loop level. Assuming mass hierarchy between the first two generations of sfermions, sleptons and the third generation, the EDMs could be dominated by the Barr-Zee diagrams. Details on this case can be found in our previous work [30]. When the phase combination φ ′ λ −φ ′ κ =0 and the phases of A u ,A d ,A e ,M 1 ,M 2 ,M 3 are nonzero, the CPV in NMSSM occurs at loop level just like the CPV MSSM case, and the constraints from the EDMs measurements are supposed to be smaller than the case that CPV occurs at tree level. In this situation, EWBG could be dominated by the higgsino-wino mixing CP-violating source, i.e., Eq. (34,35). For the case of phase combination φ ′ λ −φ ′ κ = 0 and the phases of A u ,A d ,A e ,M 1 ,M 2 ,M 3 being zero. The CPV of NMSSM could occur at tree level and the EWBG is driven by the higgsino-siglino CPV source, i.e., Eq. (36). It should be mentioned that, the cancellations of the theoretical predictions of electron EDM (eEDM) between the CP phase at tree level and the one at loop level may help evade a lot of the parameter spaces from the stringent ACME2013 constraint [80], as will be explored with the CP phases of the Higgs sectors (φ κ , φ µ ) and chargino sector (φ M 2 ). φ M 2 characterizes the couplings of H−χ ± −χ ± and thus determines the magnitude of W-loop contribution to Barr-Zee diagram of eEDM. It should be mentioned that, the phase φ M 2 also plays an important role in the coupling of W ± −χ ± −χ 0 ,  Table 1.
thus a larger magnitude of φ M 2 might make a larger contribution to eEDM through the Barr-Zee diagrams with γW ± W ∓ coupling [81].

Numerical analysis of EWBG and EDM
In this section, we calculate the Higgs spectrum using NMSSMCALC [82], implement the EDM analysis with the similar approach as Ref. [30], and show the combined results of EWBG and EDMs for two typical scenarios after imposing the constraints from Higgs-Bounds [83]: the CPV NMSSM with a small tanβ and the CPV NMSSM with a moderate tanβ (the semi-PQ limit case with κ ≈ 0). In both scenarios, we vary the imaginary parts of M 2 , µ and κ to study the different CP phases of φ M 2 , φ µ and φ κ , thereby taking into account both loop and tree level CPV effects. In both cases the neutron EDM [84] and Mercury EDM experimental constraints [85] can be satisfied for parameter spaces shown in the figures 1) .
For the eEDM calculations in above two benchmark scenarios, both cases embrace the same property, that is, the top, W-and chargino-loop Barr-Zee diagrams dominate the eEDM contributions and the cancelation among these makes the magnitude of eEDM falling into the allowed regions set by the experiment ACME 2013. We perform the numerical analysis of eEDM in the scenarios that m H 2 or m H 3 is the SM-like Higgs respectively. For the scenario with m H 2 being the SM-like Higgs, we explore the case with a small tanβ = 1.5. At last, the PQ-limit case with a relatively moderate tanβ and with λ≫κ is studied in the parameter space of φ M 2 and φ κ .

H 2 as the SM-like Higgs with a small tanβ
In this section, we study the scenario of H 2 being the SM-like Higgs and work with a small tanβ = 1.5. We investigate CPV physics in the parameter space of (φ κ ,φ M 2 ) and (φ κ ,φ µ ). Different from the case of Ref. [30], here the eEDM is mostly dominated by the H 2 and H 3 mediated Barr-Zee diagrams contributions. In this case, the pseudo-scalar a s constitutes the main proportion of H 3 , and the H 2 is dominated by h u with a s being another leading mixture component.
The left plot of Fig. 4 depicts that there is an overall larger parameter space in the plane of (φ κ ,φ µ ) where all the physical constraints are satisfied. The requirement of a H 2 mass close to 125 GeV excludes the region with relatively large φ µ . As for the BAU allowed regions, we can see it is characterized mostly by the horizontal φ κ and is not so sensitive to φ µ . This is because the higgsino-singlino CPV source drives the generation of BAU in this case, i.e., the CPV source given by Eq. 36, and the same for the case of the middle plot of Fig. 4. For the eEDM constraints, the exclusions locate at the lower-left and upper-right corners leaving the band between them as viable parameter space. The appearance of this viable band between two excluded regions is again due to a cancellation among the three contributing parts (top quark, charginos and W boson loop contributions of the γH Barr-Zee diagrams), which are plotted as contours to help understand the behavior of the eEDM constraints. More explicitly, the magnitude of both the top and W contributions decrease as φ µ increases or as |φ κ | decreases. However the W contribution comes with a minus sign and therefore these is a partial cancella-1) If we decrease the current neutron or Mercury EDM upper bound by an order of magnitude, in both cases all parameter spaces shown in the paper are ruled out, which indicates that those region will be probed in the near future by EDM experiments. tion between these two parts. The third part from the chargino contribution increases from a negative value to a large positive value as one goes from the lower-left corner to the upper-right corner. The net effect from all these three contributions gives a negative eEDM at the bottom-left corner of this plot and is excluded by experiment. This negative EDM increases gradually to zero and increases to a larger positive value as φ µ increases and as |φ κ | decreases. Therefore a eEDM allowed region is left at the middle of the plot which turns out having relatively large overlap with the BAU and the mass allowed regions.
If we change the signs of φ κ and φ µ , we then obtain the plot in the middle panel. In this case, the shape of the excluded region from the mass requirement on m H 2 is changed while the BAU allowed region shows similar behavior as previous case and is determined mostly by φ κ . For the eEDM constraints, once again, we observe cancellations among the three main contributing parts and the behavior of which can be understood with the help of the plotted contours of individual parts. Putting all these constraints together, we can see there is sizable parameter space in this plane where all physical constraints can be satisfied.
The right plot shows the constraints in the plane (φ κ ,φ M 2 ) where we observe no 125 GeV Higgs mass exclusion. As for the eEDM experiment exclusions, it is only the upper-left corner that is excluded while the majority of the parameter space gives an eEDM prediction that is compatible with the experimental limits. This behavior once again results from a cancellation among the three contributing parts. We can see there is a large parameter space that is compatible with all physical constraints in this case. Different from the the left and middle plots cases, here the higgsino-singlino CPV source Eq. (36) and higgsino-gaugino CPV source Eq. (34) and Eq. (35) together determines the BAU allowed regions.

H 3 as the SM-like Higgs for a moderate tanβ
In this section, we study the CPV scenario of NMSSM in the phase plane of (φ κ ,φ M 2 ) with the H 3 being identified as the SM-like Higgs in the approximate PQ limit (κ ≈ 0). In this scenario, H 1 is dominated by h s , H 2 is dominated by a s , and H 3 is dominated by h u 1) . The H 1,2,3 mediated Barr-Zee diagrams dominate the eEDM, all three Higgs are mixture of CP-even Higgs and the CP-odd a s .
In this case, the mass of the SM-like Higgs m H 3 is characterized by φ κ and increases from 126 GeV to 127.5 GeV as φ κ varies from the right side to the left in the region of Fig. 5. Different from the scenarios of the last section, the BAU is determined mostly by φ M 2 and is not very sensitive to the variation of φ κ . This is because the BAU is driven mostly by the higgsino-gaugino CPV sources Eq. (34) and Eq. (35) while the contribution from the higgsino-singlino CPV source term Eq. (36) is negligible since it is proportional to φ κ , which is small as a result of the definition of the PQ limit. Here we present a benchmark point in Table 2 for this moderate tanβ scenario. It should be noted that the suppressed weak sphaleron might make the generation of the BAU number harder and the PTS might not be strong enough to prevent the generated BAU being washed out. As for the eEDM constraints, we can understand the behavior of the eEDM exclusion regions from the contours of individual eEDM contributions. Firstly, the t/W Barr-Zee diagrams contributions are presented by vertical contours in this plots and is characterized solely by φ κ . Furthermore, they both give contributions that decreases as the magnitude of φ κ decreases. On the other hand, the charginos contribution to the Barr-Zee diagram of eEDM depends on both CP-phases. Its magnitude decreases first and then increase as |φ M 2 | increases or as |φ κ | decreases.
The BAU allowed horizontal band lies at the top of this figure which has sizable overlap with the allowed region from EDM constraints. This viable parameter space sits at the right-top corner of this plot corresponding to relatively small values of the two phases |φ κ | and |φ M 2 |. The combined plots of the eEDM and BAU in the plane (φκ,φM 2 ) for the scenario with H3 being the SM-like Higgs in the moderate small tanβ =3.05 of NMSSM. The color-codes are chosen to be the same as in Fig. 4. Other parameters are set as Table 2. 1) Here, the H 3 will not decay to H 1,2 due to its kinetically forbidden.  In the CPV NMSSM, there are two types of CPV sources as in aforementioned arguments, i.e., tree level and loop level. The type of the CPV in the Higgs sector is tightly related with the CPV in the Higgs mass matrix. When the the CPV occurs at the tree level, the Higgs mass matrix of the complex NMSSM are dominated by the mixing of the SM-like Higgs and heavy CP-odd Higgs like the CPV 2HDM case [29,30]. This type of CPV might be able to be detected through the associated production of Higgs with att pair [86,87] or a single t(t) [88], as well as the Higgs decay into top pairs [89][90][91] with top pair polarization being implemented at LHC [92] and linear collider [93], or the τ −lepton decays studies at LHC13 [94,95]. In the meantime, the CPV in the coupling of HZγ and HZZ can be detected through forward-backward asymmetry of the charged leptons [96] and the azimuthal angular distribution of Z boson decay states [97].

Summary and outlook
With the LHC data accumulating at the 13 TeV energy scale, we would have more access to the Higgs CP properties, which is essential to the cosmic baryon asymmetry generation at the electroweak scale together with a strongly first order phase transition [98][99][100]. In this work, we analyse the strength of the EWPT, gravitational wave production and the implementation of EWBG in the CPV NMSSM with large cancellation in the electron EDM allowed by the current eEDM measurements. We first explore the possibility that this model can provide a strongly first order EWPT required to explain the observed baryon asymmetry in the universe at the electroweak scale. Then we investigate the gravitational wave production during the strongly first order EWPT for the benchmark points and show that there exists the possibility for such gravitational waves to be detected in the future space-based interferometers like BBO and Ultimate-DECIGO. We further calculate the baryon asymmetry through the CP violating sources from either higgsino-singlino or higgsino-gaugino and analyze the constraints from the current search limit of the electron EDM. We find that the right amount of baryon asymmetry can be generated without contradicting the EDM limits for small tanβ where the 125 GeV SM like Higgs is the second lightest neutral Higgs H 2 or moderate tanβ where the 125 GeV SM like Higgs is the third lightest neutral Higgs H 3 in the CP-violating NMSSM. Such scenarios can be searched in the future for the Higgs CP properties at the LHC or some future neutron/Mercury EDM experiments (footnote 4). We note that the BAU evaluation method we adopt in this work(in the Sec.V.A) implicitly assumes the electroweak symmetry being broken inside the bubble and unbroken outside the bubble so that the weak sphaleron process is active outside the bubble and quenched inside the bubble. While, in the two benchmark models being chosen in this work the electroweak symmetry is already broken in the first step of the phase transition (which is a second order phase transition), prior to the second step first order phase transition. Therefore, the weak sphaleron process outside the bubble might be already exponentially suppressed by a factor of exp(−E sph (T )/T ), with E sph being the sphaleron energy, and thus highly suppressed the magnitude of the BAU being generated. The successful explanation of BAU requires future comprehensive studies of phase transition in the model with the electroweak symmetry being preserved outside the bubble. Furthermore, one should be aware the cosmologically domain walls problems which is introduced when the Z 3 is broken by the singlet of the model.