Electroweak Phase Transition in Non-Minimal Higgs Sectors

Higgs sector extensions beyond the Standard Model (BSM) provide additional sources of CP violation and further scalar states that help to trigger a strong first order electroweak phase transition (SFOEWPT) required to generate the observed baryon asymmetry of the Universe through electroweak baryogenesis. We investigate the CP-violating 2-Higgs-Doublet Model (C2HDM) and the Next-to-Minimal 2-Higgs-Doublet Model (N2HDM) with respect to their potential to generate an SFOEWPT while being compatible with all relevant and recent theoretical and experimental constraints. The implications of an SFOEWPT on the collider phenomenology of the two models are analysed in detail in particular with respect to Higgs pair production. We provide benchmark points for parameter points that are compatible with an SFOEWPT and provide distinct di-Higgs signatures.


Introduction
The discovery of the Higgs boson by the LHC experiments ATLAS [1] and CMS [2] has been a great success for particle physics. While the Higgs boson behaves very Standard Model (SM)like [3][4][5][6] it is clear that there must be physics beyond the SM (BSM) in order to solve the remaining puzzles that cannot be explained within the SM. Thus the observed baryon asymmetry of the Universe [7] calls for new physics extensions. Electroweak baryogenesis [8][9][10][11][12][13][14][15][16] provides a mechanism for its generation if all three Sakharov conditions [17] are fulfilled. These are baryon number violation, C and CP violation and departure from the thermal equilibrium. The asymmetry can be generated if the electroweak phase transition (EWPT) is of strong first order [14,16]. The strong first order phase transition proceeds through bubble formation and suppresses the baryon number violating sphaleron transitions in the false vacuum [18,19]. While the SM in principle fulfills all three Sakharov conditions the EWPT is not of strong first order [20]. This would require an SM Higgs boson mass of around 70-80 GeV [21] in contrast to the mass value of 125 GeV of the discovered Higgs boson [22]. Moreover, the SM CP violation arising from the Cabibbo-Kobayashi-Maskawa (CKM) matrix is not large enough [16,23]. Extended Higgs sectors beyond the SM provide additional sources of CP violation and further scalar states that help to trigger a strong first order EWPT (SFOEWPT). Thus previous studies within the 2-Higgs-Doublet Model (2HDM) [24,25], which belongs to the simplest BSM extensions, have shown that it provides a framework where an SFOEWPT can be realized in accordance with the relevant theoretical and experimental constraints, both in the CP-conserving [26][27][28][29][30][31][32][33] and in the CP-violating case [34][35][36].
With the increasing precision in the Higgs property measurements new physics extensions get more and more constrained and also models being so far successful in generating the observed baryon-antibaryon asymmetry get more and more under tension. In view of this situation, we first of all revisit our results of Ref. [35], where we investigated the CP-violating 2HDM (C2HDM) [37] with respect to an SFOEWPT, by taking into account the newest collider constraints. We then move on to the Next-to-Minimal 2HDM (N2HDM) [38,39] in order to investigate the question if it can provide an SFOEWPT and how this connects with collider phenomenology. The N2HDM is based on the extension of the CP-conserving 2HDM by a real singlet scalar field, inducing a Higgs sector consisting of three scalar, one pseudoscalar and two charged Higgs boson fields. The C2HDM and N2HDM resemble each other in the sense that they both provide at least three neutral Higgs bosons. While in the C2HDM their phenomenology is determined by the amount of CP admixture to the mass eigenstates, in the N2HDM it is their singlet admixture that governs phenomenology. Moreover, the N2HDM features more input parameters that can be tuned to render the model compatible with all theoretical and experimental constraints, so that it provides more flexibility in finding parameter points that are both compatible with these constraints and provide an SFOEWPT. We compare both models, by taking into account all relevant and recent theoretical and experimental constraints, with respect to their ability to provide an SFOEWPT. We investigate the potential to distinguish between both models in view of these requirements. In particular, the size of the trilinear Higgs self-couplings that is necessary for an SFOEWPT is analysed and their implication for Higgs pair production. We provide benchmark points that connect the requirement of an SFOEWPT with specific features in the collider phenomenology, in particular Higgs pair production. Both models are interesting here also because their non-minimal Higgs sectors allow for the production of mixed Higgs pair final states. With our analysis we hope to provide a link between collider phenomenology and cosmology.
The paper is organized as follows. In Sec. 2 we introduce the computation of the effective potential in the C2HDM and N2HDM and set our notation. In Sec. 3, we describe the renormalisation of the N2HDM that is new. The one for the C2HDM has been provided previously in [35]. The basics of our numerical analysis are introduced in Sec. 4. Section 5 is devoted to the presentation of our results, where we first give in Subsec. 5.1 an update of the C2HDM by including the newest constraints. We then move on to the detailed presentation of the N2HDM phenomenology of the SFOEWPT in Subsec. 5.2, both discussing the related mass spectra and trilinear Higgs self-couplings and providing benchmark points. In Subsec. 5.3 we compare the C2HDM and N2HDM rates for Higgs pair production induced by the requirement of an SFOEWPT. Our conclusions are given in Sec. 6.

The Effective Potential
In order to decide if the EWPT is of strong first order, we have to determine the value v c of the VEV at the critical temperature T c . The VEV at the temperature T is given by the minimum of the one-loop corrected effective potential at non-zero temperature T . In this section we provide the loop-corrected effective potentials at finite temperature for our two models under investigation, the C2HDM and the N2HDM. We start with the tree-level potentials and thereby set our notation.

The Tree-Level C2HDM Potential
We briefly introduce the C2HDM Higgs sector and refer to [35,40] for a more detailed introduction. In 2HDMs [24,25] the SM Higgs potential is extended by an additional SU(2) L scalar doublet yielding the Higgs potential, with a softly broken discrete Z 2 symmetry under which Φ 1 → Φ 1 , Φ 2 → −Φ 2 , which ensures the absence of tree-level flavour-changing neutral currents (FCNC) upon extension to the Yukawa sector. The hermiticity of the Lagrangian requires all couplings to be real except for m 2 12 and λ 5 . If their complex phases are unrelated the model is CP-violating and called C2HDM [37]. In the following, we will adopt the conventions of Refs. [40,41]. We denote the VEVs of the EW minimum by ω i ∈ R (i = 1, 2) and write the two Higgs doublets as expansions around the VEVs in terms of the charged field components ρ i and η i and the neutral CP-even and CP-odd fields ζ i and ψ i . At tree level, the general vacuum structure of the 2HDM allows for three different possible vacua that are given by the normal EW-breaking vacuum, a CP-breaking and a chargebreaking (CB) vacuum. It has been shown that vacua breaking different symmetries cannot coexist at tree level in the 2HDM [42][43][44]. Since this statement might not hold at higher orders, or be broken by finite temperature effects, we allow for a more general vacuum structure. We therefore include the possibility of a CB-and CP-breaking VEV in the field expansion, denoted by ω CB and ω CP , respectively, , where the brackets around the doublets stand for their vacuum state. The VEVs of our present vacuum at zero temperature are denoted by A non-zero CB VEV would break the conservation of electric charge and introduce massive photons so that we will neglect in our analysis parameter points evolving such unphysical vacuum structures. The VEVs of the normal EW minimum are related to the SM VEV v ≈ 246 GeV The angle β is defined by the ratio of v 1 and v 2 , Using the minimum condition of the potential We use Eq. (2.9a) and Eq. (2.9b) to trade m 2 11 and m 2 22 for v 1 and v 2 . The two possible sources of CP violation are related to each other through Eq. (2.9c) so that one independent CP-violating phase remains in the C2HDM. Rotating the fields ψ 1 and ψ 2 with β yields where G 0 denotes the neutral Goldstone boson. The three neutral mass eigenstates H k (k = 1, 2, 3) of the C2HDM are obtained by the rotation of the gauge eigenstates ζ k (k = 1, 2, 3) to the mass basis, with the rotation matrix R parametrized in terms of the mixing angles α k as (cos (2.12) Without loss of generality, the rotation angles α 1,2,3 can be chosen in the interval (2.13) The mass eigenvalues are given by 14) The mass eigenstates are ordered by ascending masses as Using the minimum conditions and the rotation to the mass eigenstates, the following set of nine independent parameters of the C2HDM remains [45], v , tan β , α 1,2,3 , m H i , m H j , m H ± and Re(m 2 12 ) . (2.16) The m H i and m H j denote any of the masses of two among the three neutral Higgs bosons. The third neutral Higgs boson mass is determined through One of the three Higgs bosons is identified with the SM-like Higgs boson with the measured mass value of m h = 125.09 GeV [22]. As mentioned above, extending the Z 2 symmetry to the Yukawa sector ensures that each type of the up-and down-type quarks and charged leptons can couple to only one of the two Higgs doublets so that FCNCs are avoided at tree level. In Tab. 1 the different types of the 2HDM, type I, type II, lepton-specific and flipped, are listed. In this analysis we will concentrate on type I and type II.

The Tree-Level N2HDM Potential
We give a brief introduction in the N2HDM and refer for more details to [39]. The tree-level potential of the N2HDM consists of a CP-conserving 2HDM which is extended by a real singlet field Φ S . The potential is invariant under two discrete Z 2 symmetries. We call the one given by the generalisation of the 2HDM symmetry to avoid tree-level FCNCs Z 2 under which We call the second one Z 2 . It is defined as The most general tree-level potential invariant under these transformations apart from a term proportional to m 2 12 that softly breaks Z 2 , reads [39] where all parameters are real due to CP-conservation. After electroweak symmetry breaking the two Higgs doublets and the real singlet acquire VEVs ω i ∈ R. As a first analysis of the N2HDM vacuum structure performed in Ref. [39] has shown, the N2HDM exhibits a different vacuum structure than the 2HDM. The impact of the N2HDM vacuum structure has been further studied in [46] by applying the method of Ref. [47]. Additionally, loop corrections and finite temperature effects may change the vacuum texture, so that we expand the doublet and singlet fields including the most general vacuum structure, as 21) where we have expanded the Higgs fields in terms of the charged field components ρ i and η i (i = 1, 2), and the neutral CP-even and CP-odd fields ζ i and ψ i , respectively. As mentioned above, already the tree-level zero-temperature vacuum can exhibit a quite general structure as compared to the 2HDM. For simplicity we choose for our present vacuum at zero temperature a vacuum structure with the same properties as in the C2HDM, given by As in the (C)2HDM, the ratio of the electroweak VEVs v 1 and v 2 is defined by the angle β, Requiring the tree-level potential in Eq. (2.20) to be minimized at the electroweak vacuum Eq.
Equation (2.27) allows to trade the parameters m 2 11 , m 2 22 and m 2 S for the VEVs v 1,2,S resulting in the mass matrix for the CP-even neutral fields ζ 1,2,3 at zero temperature The matrix R rotates the mass matrix M 2 Scalar in the mass eigenstates H 1,2,3 with the mass eigenvalues We use the same convention as in the C2HDM where the mass eigenstates are ordered by ascending masses as The N2HDM has 12 real independent parameters where we chose as many parameters as possible with a physical meaning. We take Eq. (2.27) to trade the potential parameters m 2 11 , m 2 22 and m 2 S for the SM and singlet VEVs v and v S , and tan β. Moreover, we express the quartic couplings λ i in terms of the physical masses and mixing angles. The Z 2 breaking mass term m 2 12 is kept as independent parameter. This yields the following set of input parameters To avoid FCNC in the N2HDM the same types as in the C2HDM, given in Tab. 1, can be used since the singlet field Φ S does not couple to fermions. For more details and a phenomenological discussion of the N2HDM we refer to [39].

One-loop Effective Potential at Finite Temperature
In the following we briefly repeat the main ingredients for the one-loop effective potential at nonzero temperature. For details, we refer to [27,35,48] which can be generalized to the N2HDM.
The loop-corrected effective potential splits into the tree-level potential V tree , the one-loop correction at zero temperature given by the Coleman-Weinberg (CW) potential V CW , and the temperature-dependent part V T , so that it reads where V tree is given by Eq. . The CW potential [49] in the MS scheme is given by where n j are the degrees of freedom, s j the spin and m j the mass of the specific particle j. The sum extends over the Higgs and Goldstone bosons, the massive gauge bosons, the photon and the fermions, so that both for the C2HDM and the N2HDM we have The mass m j is the mass eigenvalue of particle j obtained from the tree-level mass matrix expressed in terms of the general VEV configuration ω i , with i = 1, 2, CP, CB for the C2HDM and i = 1, 2, S, CP, CB for the N2HDM. Applying the Landau gauge in the calculation of the CW potential allows us to drop the ghost contributions in the analysis, but we need to account for the possibility that the Goldstone bosons get massive. The Goldstone bosons as well the photon are massless at T = 0, but for field configurations different from the tree-level VEVs at T = 0, which is needed for the minimization procedure, they can acquire an effective mass term. Furthermore, we allow for unphysical vacuum structures with non-zero ω CB inducing additionally unphysical masses. We also account for the possibility of the generation of a CP-violating VEV ω CP . For the neutral scalars Φ 0 = H 1 , H 2 , H 3 , A, G 0 , the charged scalars Φ ± = H + , H − , G + , G − , the leptons l + ,l − , quarks q, q and longitudinal and transversal gauge bosons V L = Z L , W + L , W − L , γ L and V T = Z T , W + T , W − T , γ T the degrees of freedom in Eq. (2.35) read The renormalisation scale is chosen as µ = v ≈ 246.22 GeV and the MS renormalisation constants read The thermal corrections V T comprise the daisy resummation [50] of the n = 0 Matsubara modes of the longitudinal components of the gauge bosons W ± L , Z L , γ L and the bosons Φ 0 , Φ ± . This Debye correction adds to their masses at zero temperature. The thermal corrections V T can be cast in the form [51,52] The sum extends over k = W ± T,L , Z L,T , γ L,T , Φ 0 , Φ ± , l ± , q, q. Denoting the mass eigenvalues including the thermal corrections for the particles k by m k , we have for J ± (cf. e.g. [53] with the thermal integrals where J + (J − ) applies for k representing a fermion (boson). The general formulae for the thermal masses can be found in [48]. For the numerical evaluation of the effective potential at finite temperatures and the further minimization procedure we use the code BSMPT v1.1.2 [48] and we refer to [48] for a more detailed discussion of the used numerical approximations in the thermal integrals in Eq. (2.39).

Renormalisation
At one-loop level the masses and mixing angles differ from those extracted from the tree-level potential. In order to take into account the one-loop effects and at the same time enable an efficient scan in the parameter space of the model, we renormalize the loop-corrected masses and mixing angles such that they are equal to their tree-level values. This allows us to use them directly as input values for our scan. The scheme has been introduced in Ref. [27] where it was applied to the 2HDM. In [35] it was extended to the C2HDM. We therefore show here only the procedure for the N2HDM. The renormalised loop-corrected potentialV is obtained by adding the counterterm potential V CT with the counterterm potential given by where p i stands for the parameters of the tree-level potential. Furthermore, for each field direction φ k which is allowed to develop a VEV, a tadpole counterterm δT k is introduced. For the N2HDM these are δT 1 , δT 2 , δT S , δT CP and δT CB . The renormalisation conditions for the scheme described above read and where φ T =0 denotes the tree-level vacuum state at zero temperature. The conditions Eq. (3.43) and Eq. (3.44) ensure that at zero temperature the tree-level minimum remains a local minimum at one-loop level. We check numerically if it also the global one. Additionally, the second set of conditions in Eq. (3.44) ensures that the masses and mixing angles derived from the loopcorrected effective potential remain at their tree-level values. In general, the renormalisation conditions result in an overconstrained system of equations which can be solved by imposing additional assumptions. In the case of the N2HDM Eqs. (3.43) and (3.44) lead to a twodimensional solution space which can be fixed by imposing Solving then for the counterterm parameters yields and In Ref. [58] formulae for both the first and the second derivatives of the CW potential have been derived in the Landau gauge basis.
4 Numerical Analysis

Minimisation of the Effective Potential
The EWPT is of strong first order if the baryon-washout condition is met which requires that the ratio ξ c of the critical VEV v c and the critical temperature T c is larger than one [11,52,59], The VEV v at the temperature T is given by where the ω i are the field configurations that minimise the loop-corrected effective potential at finite temperature T . Note that we do not include the singlet VEV ω S (present in the N2HDM) in Eq. (4.50), but we take ω S into account for the minimisation procedure. Since the electroweak sphaleron couples only to particles charged under SU(2) L , the singlet VEV can be dropped in the calculation of the critical VEV in Eq. (4.50). The critical temperature T c is defined as the temperature where the potential has developed two degenerate minima. In order to compute the global electroweak minimum of the one-loop corrected effective C2HDM and N2HDM potentials of Eq. (3.41) we implemented both models in BSMPT v1.1.2 [48] which also calculates the strength ξ c of the phase transition.

Constraints and Parameter Scan
For simplicity, in the numerical analysis we discuss only Type I and II of the models. The parameter samples used for the numerical investigation in this paper have to satisfy theoretical and experimental constraints. We obtained them by performing scans in the parameter spaces of the C2HDM and N2HDM, respectively. For the scans we required that one of the neutral Higgs bosons, denoted by h in the following, behaves SM-like and has a mass of m h = 125.09 GeV [22]. The scan ranges for all the input parameters are given in Tab. 2 for Type I (T1) and in Tab. 3 for Type II (T2) of the N2HDM, for the C2HDM they are given in Tab. 4 for Type I and in Tab. 5 for Type II. We introduce here the notation m H ↓ /H ↑ for the masses of the lighter/heavier of the two non-SM-like neutral Higgs bosons.
As for the SM parameters of our analysis, we use the fine structure constant taken at the Z boson mass scale [60,61],    To be consistent with the CMS and ATLAS analyses, we take the on-shell top quark mass as [62] m t = 172.5 GeV The CKM matrix is taken to be real, with the CKM matrix elements given by [60] Finally, the electroweak VEV is set to To be consistent with recent flavour constraints, we test for the compatibility with R b [63,64] and B → X s γ [64][65][66][67][68] in the m H ± − tan β plane. For T2, this implies that the charged Higgs mass has to be above 580 GeV [68] whereas in T1 this bound is much weaker and is strongly in GeV in GeV 2 125.09 [30,1500] [30,1500] [30,1500] 10 −3 , 10 5  correlated with tan β. 2 Compatibility with the electroweak precision data is checked through the oblique parameters S, T and U [69] where we apply the general procedure for extended Higgs sectors given in [70,71]. Including the full correlations, we demand 2σ compatibility with the SM fit [72].
In order to avoid degenerate Higgs signals we impose a mass window between the non-SM-and SM-like Higgs bosons so that masses |m H i =h − m h | < 5 GeV are excluded from the analysis. We use the program ScannerS [39,[73][74][75] to search for valid parameter points. The program allows to check for boundedness from below of the tree-level potential and, for the C2HDM, uses the tree-level discriminant of [76] to ensure the electroweak vacuum to be the global minimum at tree level. For the N2HDM, all tree-level minima have been implemented in ScannerS and are compared numerically to find the global minimum. To check for consistency with the Higgs exclusion limits from LEP, Tevatron and LHC HiggsBounds5.5.0 [77][78][79] is used, and the SM-like Higgs rates are required to be within the 2σ range of the SM which is checked by HiggsSignals2.3.0 [80]. The required decay widths and branching ratios are obtained from C2HDM HDECAY [40] for the C2HDM and from N2HDECAY [39,81] for the N2HDM. Both codes are based on the implementation of the C2HDM and the N2HDM each in the existing code HDECAY [82,83]. Due to the CP violation, in the C2HDM also compatibility with the electric dipole moments [84] has to be checked.The most stringent limit is provided by the ACME collaboration [85]. Moreover, we take into account for both models the impact of the recent di-Higgs searches in the final states 4b [86,87], (2b)(2τ ) [88,89] and (2b)(2γ) [90] on the viable parameter space. These searches are also implemented in HiggsBounds and HiggsSignals.
In addition to these checks, we impose the requirement that the tree-level minimum of the potential is still the global electroweak minimum at one-loop level. The one-loop minimum is determined by numerical minimisation of the one-loop potential at zero temperature and checked 2 Many of the experimental constraints applied on the 2HDM also hold for the N2HDM, since these constraints are only sensitive to the charged Higgs boson so that the calculation of the 2HDM can be taken over to the N2HDM [39]. In this way we are able to use the same constraints for the C2HDM and N2HDM.
against the tree-level value. In the following, all parameter points providing an SFOEWPT also have an NLO stable global minimum at zero temperature. Furthermore, we demand an approximated NLO unitarity. The tree-level perturbative unitarity relations for the N2HDM read [39,48] and the eigenvalues a 1,2,3 that are the real roots of the cubic polynomial By replacing λ i → λ i + δλ i in Eq. (4.60), the NLO effects on the unitarity of the S-matrix can be approximated and checked in the parameter scan. The counterterms δλ i in the N2HDM are given in Eq. (3.46) and are calculated numerically during the minimisation procedure. The corresponding relations for the C2HDM counterterms and perturbative unitarity relations are given in [25,40,48]. By imposing these additional constraints at NLO, it can happen that the parameter sample is significantly reduced. In Tab. 6 a list of the remaining parameter sample is given after imposing the NLO constraints and the requirement for an SFOEWPT. These numbers have to be taken with a grain of salt, however, as dedicated scans adapted to specific requirements like e.g. not loosing points due to NLO unitarity, would change this picture. The numbers are simply meant to show that NLO constraints and/or the requirement of an SFOEWPT have an effect on the parameter sample.
The Higgs spectra of the C2HDM and N2HDM consist of at least three neutral Higgs bosons (in the N2HDM, we additionally have the pseudoscalar). Denoting the lighter and heavier non-SM-like Higgs bosons with masses m H ↓ and m H ↑ by H ↓ and H ↑ , we can have three different mass configurations that we will refer to in the following as heavy mass hierarchy with the mass hierarchy  semi-inverted mass hierarchy with and inverted mass hierarchy with the hierarchy

Results
In the following analysis we investigate to which extent the viable parameter spaces of the C2HDM and N2HDM are constrained by the additional requirement of a strong first order EWPT. This also allows us to investigate the differences that arise due to CP violation on the one hand and singlet admixture on the other hand. Since we discussed the implications for the phenomenology of the C2HDM already in detail in Ref. [35], we start by providing a rather short update of our analysis of the C2HDM by taking into account the newest results for the Higgs data. Subsequently, we discuss in detail the interplay between a strong first order phase transition and the collider phenomenology of the N2HDM. In particular, we study the impact of ξ c ≥ 1 on the size of the trilinear Higgs self-couplings and the overall mass spectrum. We will provide benchmark scenarios that highlight the connection between a successful SFOEWPT and collider phenomenology with emphasis on Higgs pair production at the LHC. In the end a comparison between the two models will be drawn and the characteristic differences in their phenomenology with respect to a successful SFOEWPT will be discussed. For our analysis we produced for each model of each type about half a million parameter points that respect the theoretical and experimental constraints listed in Sec. 4.2. These parameter points are checked for a successful SFOEWPT by using the program BSMPT v1.1.2 [48].

C2HDM -Update
We start with an update of our analysis in Ref. [35] by taking into account the new collider constraints that have been implemented in HiggsSignals 2.3.0 and HiggsBounds 5.5.0. In addition, we increased the scan ranges of all scalar masses from 1 TeV up to 1.5 TeV with the aim to find new valid parameter points featuring a heavy scalar spectrum that provides an SFOEWPT. The subsequent discussion will show that no additional heavy parameter points were found fulfilling the requirement of an SFOEWPT and compatibility with recent collider

constraints.
We start the discussion with Fig. 1, where the charged Higgs boson mass m H ± is shown as a function of tan β for the C2HDM T1 (left) and T2 (right). The grey points are all parameter points compatible with the theoretical and recent experimental constraints as described in Sec. 4.2. The brown points additionally provide an NLO stable vacuum and fulfill NLO perturbative unitarity (see also Sec. 4.2). The color code indicates parameter points with values of ξ c ≥ 1 and thereby all points with an SFOEWPT. In the C2HDM T1 two distinct possible scenarios for parameter points providing an SFOEWPT can be observed. The first region has charged Higgs boson masses of ∼ 450 GeV up to ∼ 690 GeV and quite small tan β values around 1. Only one point provides an SFOEWPT with a medium charged Higgs boson mass and a tan β value around 5. All the other points of the second region with larger values of tan β have a charged Higgs boson mass below 200 GeV. Compared to our previous analysis [35], parameter points with medium charged Higgs boson masses and large tan β values could not be found any more, so that we have this strict separation of small masses in combination with large tan β values and medium masses in combination with small tan β. The maximum strength of the phase transition that we found for the C2HDM T1 is ξ c ≈ 1.7, and for the C2HDM T2 it is ξ c ≈ 1.18, which is compatible with our findings in [35].
In the C2HDM T2, the flavor constraints B → X s γ [64][65][66][67][68] require the charged Higgs boson mass to be above 580 GeV, which is reflected in Fig. 1 (right). Compatibility with the electroweak precision observables enforces the degeneracy of two Higgs boson masses, so either m H ± ≈ m H ↓ ,H ↑ ,h or a pair of the neutral Higgs bosons are mass degenerate. In order to quantify this effect, we first of all look for that pair of Higgs masses (m i , m j ) that has the minimum mass difference out of all possible Higgs mass pairings, i.e. we define the mass gap of the (almost) degenerate pair (m H i , m H j ) as The requirement of an SFOEWPT tightens this mass gap to even smaller values Type II : max In both types of the C2HDM the mass degeneracy that is realized in most of the points is the one of H ↓ and H ↑ , as shown in Fig. 2. It shows the non-SM-like Higgs boson mass plane m H ↑ − m H ↓ with the color code denoting the values of ξ c above 1, hence an SFOEWPT. In the other cases the mass degeneracy occurs with the charged Higgs boson, i.e. m H ± ≈ m H ↓ or m H ± ≈ m H ↑ .
In the C2HDM T1 in Fig. 2(a) two mass hierarchies are possible, the heavy mass hierarchy m h < m H ↓ < m H ↑ and the fully inverted mass hierarchy m H ↓ < m H ↑ < m h . As we will see later, in the N2HDM the semi-inverted and the heavy mass hierarchies are possible in T1 and T2.
As noted in [27] the SFOEWPT favors a light scalar spectrum, where additional heavy degrees of freedom, that do not obtain an SU(2) L VEV, can help to strengthen the phase transition. At the same time, the mass degeneracy of H ↑ and H ↓ enforces the overall scalar spectrum to be in a medium range of O(m H ↓ , m H ↑ ) ∼ 500 GeV. The parameter region with could also produce an SFOEWPT as remarked in [27], since this region has a light H ↓ enabling the SFOEWPT and in the meantime the heavy degree of freedom H ↑ could build up a deep potential barrier between the symmetric and broken minimum. Since this parameter region is restricted by the di-Higgs searches through Higgs decays H ↑ → H ↓ H ↓ /hh, the parameter space in this region is already sparse due to the collider constraints. Consequently, with the updated HiggsSignals 2.3.0 and HiggsBound 5.5.0 versions taking into account the recent di-Higgs searches, it is more involved to find parameter points compatible with the collider constraints and an SFOEWPT. As we will see later, this restriction can be circumvented in the N2HDM due to the singlet admixture.
In the C2HDM T2 the overall scalar spectrum in Fig. 2(b) is heavier compared to T1. This is due to the already required heavy charged Higgs boson mass and the small mass gaps. As the SFOEWPT still favors a light scalar spectrum, the only parameter points providing an SFOEWPT are found in the edge of smallest masses. The overall order of neutral non-SM like Higgs masses providing an SFOEWPT is also O(m h i ) ∼ 500 GeV as in the C2HDM T1.
To conclude the C2HDM update, we find that the mass spectrum compatible with the recent collider and theoretical constraints, for T1 is mainly constrained by the recent di-Higgs measurements whereas for T2 the flavor constraints are the most restrictive ones. Future analyses with increasing constraining power in these mass regions might exclude significant regions of the C2HDM parameter space providing an SFOEWPT.

N2HDM -Phenomenology of the SFOEWPT
In the following we discuss the implications of an SFOEWPT on the phenomenology of the N2HDM. The N2HDM has one more degree of freedom compared to the C2HDM due to the additional singlet in the Higgs sector. The larger number of free parameters in this model reduces the influence of the constraints on the parameter space so that heavier Higgs spectra compatible with an SFOEWPT are still possible. In this kind of scenarios, however, one of the non-SM like Higgs bosons is almost completely singlet-like. The maximum values of ξ c that we find are ξ N2HDM T1 c = 2.04 and ξ N2HDM,T2 c = 1.43 for the N2HDM T1 and T2, respectively. We start with the investigation of the mass spectrum followed by the discussion of the trilinear Higgs self-couplings. Afterwards, we will present several benchmark points providing interesting scenarios with different features.

Mass Spectrum of the N2HDM T1
For the following discussion of the N2HDM we introduce the singlet admixture Σ H i of the respective CP-even Higgs boson H i as It describes the amount of admixture of the singlet field ρ S to the corresponding mass eigenstate H i . In Fig. 3   As these mass windows show, in the N2HDM T1 the heavy mass hierarchy (m h < m H ↓ < m H ↑ ) and the semi-inverted mass hierarchy (m H ↓ < m h < m H ↑ ) are possible, whereas the inverted hierarchy (m H ↓ < m H ↑ < m h ) is not realised. We will provide benchmark scenarios for all possible cases. Like in the C2HDM, an SFOEWPT favors light Higgs mass spectra below ∼ 500 GeV except for spectra with singlet-like H ↑ . If the heaviest CP-even non-SM like Higgs boson H ↑ has a singlet admixture above Σ H ↑ 80%, the SFOEWPT opens the window for larger masses m H ↑ . In case of singlet admixtures below about 80%, Σ H ↑ 80%, on the other hand the same preference for intermediate Higgs mass spectra as in the C2HDM can be observed for ξ c ≥ 1. The possibility of large neutral Higgs boson masses allows for heavy Higgs decays into pairs of lighter Higgs bosons, H ↑ → H i H j . Since almost all of these heavy states are singlet-like, their couplings to SM particles are suppressed and the Higgs-to-Higgs decay channel becomes an important discovery channel. We will provide benchmark scenarios where these channels may become accessible at the LHC.
As already mentioned, compatibility with the EW precision data checked through the S, T and U parameters requires a small mass gap between the charged Higgs boson H ± and one of the neutral Higgs bosons or between two neutral Higgs bosons, and the SFOEWPT enforces even more the mass degeneracy between at least one non-SM like Higgs boson pair. In Fig. 4(a) we show the minimum mass gap out of all possible neutral Higgs pairings, ∆ m H i ,m H j , defined analogously to Eq. (5.65), versus the maximum mass in the spectrum, max(m H i ). The color code is the same as in Fig. 3. exceptions of up to 100 GeV. Scenarios with degenerate neutral Higgs boson masses are rather rare so that in general the charged Higgs boson mass lies in the same region as the neutral Higgs boson masses. Consequently, the mass spectrum of the charged Higgs boson is also reduced in case of an SFOEWPT. This is reflected in Fig. 4(b) where the charged mass m H ± is depicted versus tan β, with the same color coding as in Fig. 4(a).
The last phenomenological effect on the N2HDM mass spectrum induced by the SFOEWPT we want to discuss, is shown in Fig. 5. It displays the mass difference m A − m H ± versus tan β. This plots shows that an SFOEWPT does not only allow for scenarios with m A ≈ m H ± , but also mass spectra with large gaps between m A and m H ± are allowed. This possibility should also be taken into account when the EWPT in investigated. For parameter points in M deg , stronger EWPTs with ξ c up to ξ c 2.1 can be observed compared to the points in M sep with ξ c 1.93. To discuss the slight tendency of a stronger EWPT in M sep for increasing mass gaps m A − m H ± , we first note the observations made in Refs. [35,53]: The strength of the phase transition ξ c increases with the size of the couplings of the light bosonic particles to the SM-like Higgs boson and decreases with the Higgs boson mass. Additionally, particles that contribute to the EWPT necessarily have a non-vanishing electroweak VEV. All non-SM-like neutral Higgs bosons H ↓ , H ↑ receive an electroweak VEV through mixing and therefore, for an SFOEWPT, their masses m H ↓ , m H ↑ have to be either light or their VEV has to be small. We note here again explicitly that we take the singlet VEV v S into account for the minimisation of the effective potential, but we do not include v S in the calculation of the electroweak VEV in Eq. (4.50). So the EWPT is not directly affected by the singlet VEV, just indirectly through the minimisation. Thus particles that do not obtain an electroweak VEV (A, H ± ) and singlet-like CP-even neutral Higgs bosons (Σ i ≈ 1) are still allowed to be heavy without decreasing the strength of the EWPT. The heavy degrees of freedom even help to strengthen the EWPT by enabling a deeper potential barrier between the broken and symmetric phase [29]. For increasing mass gaps m A − m H ± and therefore increasing m A masses, we have an additional heavy degree of freedom in the spectrum, whereas the bosonic degrees of freedom that obtain an electroweak VEV remain light. Consequently, the strength of the EWPT will increase for the parameter points in M sep with increasing mass gap. Since there are additional interplays between ξ c and the mass spectrum or other effects, the effect on the size of ξ c is not significantly enhanced, however.

N2HDM T1 -Trilinear Higgs Self-Couplings
The one-loop-corrected trilinear Higgs self-couplings are obtained from the one-loop effective potential by performing the third derivative with respect to the Higgs fields. The problem of infrared divergences related to the Goldstone bosons in the Landau gauge is treated analogously to the extraction of the masses from the second derivative of the potential, cf Ref. [58] for details. In Fig. 6 Figure 6: N2HDM T1: The leading-order trilinear self-coupling between three SM-like Higgs bosons normalised to the SM reference value versus the NLO trilinear self-coupling between three SM-like Higgs bosons normalised to the SM reference value. The color code is the same as in Fig. 3. The di-Higgs search constraints are included.
As observed in [35] the SFOEWPT favors large trilinear Higgs self-couplings which is also observed here, since the strongest EWPTs are located at the maximum values of the trilinear Higgs self-couplings. On the other hand, the upper bound of the trilinear Higgs self-coupling is significantly reduced by the SFOEWPT which can be explained by the interplay between the quartic coupling and the masses of the Higgs bosons participating in the EWPT. Note, that besides the dominant top-loop contributions to the NLO coupling, also the Higgs-loop contributions present in the C2HDM and N2HDM can be large. The masses of the heavy N2HDM Higgs bosons Φ can be cast in the following schematic form (see [91] for the 2HDM) with M 2 denoting the mass scale independent of the VEVs and f v,m,s (λ i ) a linear combination of the quartic couplings of the Higgs potential. For an SFOEWPT large couplings λ i are required. On the other hand, the masses should not become too heavy, which we observed in the previous discussion of the general mass spectrum, thus limiting the maximum values for the quartic coupling due to Eq. (5.73). This explains why we observe the strongest EWPT for the largest trilinear coupling, but the maximum enhancement of the trilinear Higgs self-coupling for a strong EWPT remains below the value compatible with the applied constraints.
Since the N2HDM, in contrast to the C2HDM, has a singlet VEV which does not contribute to the EWPT, heavy degrees of freedom are compatible with an SFOEWPT if they are singletlike. In Fig. 7(a) we display the NLO trilinear Higgs self-coupling between three heavy Higgs bosons H ↑ normalised to the NLO trilinear Higgs self-coupling of the SM versus one minus the singlet admixture Σ H ↑ of H ↑ (which corresponds to the doublet admixture of H ↑ ). The meaning for the grey and brown points is the same as in the previous plots. The colored points now, however, denote the mass value m H ↑ for those parameter points that are compatible with an SFOEWPT. We observe that all singlet admixtures from 0% up to 100% are possible for intermediate mass ranges, but only for the singlet-like heavy states we observe large masses m H ↑ of up to 1490 GeV. For these heavy masses the trilinear self-coupling λ NLO H ↑ H ↑ H ↑ is enhanced compared to the parameter points with intermediate masses. Simultaneously these large masses enhance the trilinear coupling λ hhH ↑ opening a window for heavy Higgs decay chains like H ↑ → hh → XX. The enhancement can be observed in Fig. 7(b) where the NLO trilinear self-coupling λ hhH ↑ normalised to the NLO trilinear Higgs self-coupling of the SM is plotted against one minus the singlet admixture Σ H ↑ , with the same color code as in Fig. 7(a). Like for the self-coupling λ H ↑ H ↑ H ↑ the largest enhancements of λ hhH ↑ compatible with an SFOEWPT are observed for large singlet admixtures. These heavy Higgs decays allow to distinguish between the C2HDM and the N2HDM as we will discuss later.

N2HDM T1 -Benchmark Points
In the following we will present benchmark points that provide an SFOEWPT and have interesting phenomenological features. In total we generated 271743 parameter points with ScannerS, fulfilling theoretical and experimental constraints. After applying the NLO constraints, NLO vacuum stability and NLO unitarity, and demanding an SFOEWPT 920 parameter points are left in the sample. 836 of these points feature the heavy mass hierarchy while 84 have the semi-inverted.

Semi-inverted Mass Hierarchy
In Tab. 7 two benchmark points are listed with a semi-inverted mass hierarchy where   for a c.m. energy of √ s = 13 TeV. We then have for BPii2 the following signal rates for SM-like di-Higgs production in the 4b, (bb)(τ τ ), (bb)(W W ), 4W and (bb)(γγ) final states,

Mass Spectrum of the N2HDM T2
We start our discussion of the N2HDM T2 mass spectrum with the minimum mass gap ∆ m H i ,m H j among all possible pairings of the N2HDM Higgs bosons. Figure 8 Higgs bosons with masses above ∼ 800 GeV are again singlet-like states. In Fig. 8(b) the mass gap between A and H ± versus the mass gap between H ↓ and H ↑ is shown. Unlike the results in the N2HDM T1, we do not find the separation into the two distinct parameter regions M deg and M sep . In the N2HDM T2, only parameter points with nearly mass degenerate A and H ± fulfill the experimental and theoretical constraints and simultaneously provide an SFOEWPT.
In Fig. 9 we

N2HDM T2 -Trilinear Higgs Self-Couplings
In Fig. 11  As already mentioned in the discussion of the N2HDM T1, the interplay of the requirement of large quartic couplings and small to medium Higgs boson masses, pushes the trilinear couplings to enhanced values, but their upper limit remains under the allowed possibilities concerning theoretical and experimental constraints. In contrast to the N2HDM T1, in the N2HDM T2 trilinear couplings of at least 1.52 times the SM value are required for an SFOEWPT, the SM size is not sufficient.

Comparison of the C2HDM and the N2HDM Di-Higgs Rates
So far we discussed both models independently with respect to their phenomenology of the mass spectra and the trilinear Higgs self-couplings. In the following we will compare the two models with respect to their expected signal rates in the di-Higgs final states. We start the comparison for the (4b) final state of the T1 models. In Fig. 12 we show the signal rates for the process pp → H ↑ → H ↓ → bb H ↓ → bb versus the mass of the heaviest neutral Higgs boson H ↑ for the C2HDM (left) and the N2HDM (right). Figure 13 displays the signal rates pp → H ↑ → h → bb h → bb versus m H ↑ . The rates decrease with increasing mass m H ↑ . In Fig. 12 (left) a strong reduction of the C2HDM rates can be observed at the threshold for the decay into a top-quark pair, m H ↑ = 350 GeV, down from O(1 pb) by more than three orders of magnitude. 4 In the N2HDM, no such reduction appears as the addition of the singlet to the Higgs sector allows for singlet-like H ↑ granting Higgs-to-Higgs decays being of the order O (1 pb). Reaching m H ↑ ≈ 700 GeV also H ↓ can be heavier than 350 GeV so that the decays H ↓ → tt can become possible. In the C2HDM, this possibility reduces the overall cross section so effectively that no Higgs-to-Higgs decays above O 10 −10 pb are realized in the C2HDM for masses above 700 GeV. With respect to a successful SFOEWPT the differences of both models are significant. An analogous behaviour of the maximum rates is observed in Fig. 13 in the decay channel In the SM-like decay chain H ↑ → hh → XX (and also in the decay chain H ↑ → H ↓ H ↓ → XX) increasing sensitivity in the di-Higgs searches for intermediate Higgs boson masses (up to 500 GeV) could allow for probing valid N2HDM candidates with a successful strong first order EWPT, in contrast to the C2HDM. In the C2HDM the cross section for the SFOEWPT parameter points is orders of magnitudes below the LHC sensitivity. Since in C2HDM and N2HDM T1 the couplings of the Higgs bosons to the down-type quarks and leptons are the same [39,40] the rates in the (bb)(τ τ ) final state are simply reduced by about a factor 10 without changing the qualitative behaviour so that we do not show separate plots for this final state.

bb (γγ) Final State
The rates for the (bb)(γγ) final states from H ↑ production with subsequent decay into H ↓ H ↓ are shown in Fig. 14 and those from H ↑ → hh in Fig. 15. We observe the same behaviour as in the final state bb bb with the difference that the overall signal strength is smaller due to the smaller branching ratio for the Higgs decays into photons. Again, in the C2HDM the signal rate is significantly reduced by requiring an SFOEWPT whereas this is not the case in the N2HDM. The top threshold reduces the overall signal in the C2HDM while again the N2HDM is not affected due to the possibility of the singlet admixture.

C2HDM and N2HDM T2
We now turn to the comparison of the C2HDM and N2HDM T2 models. We display the signal rates σ pp → H ↑ → h → bb h → bb versus m H ↑ in Fig. 16. In the left plot, the results for the C2HDM are shown and in the right plot those for the N2HDM. The signal rates obtained in the N2HDM and C2HDM T2 are smaller than in T1 of the models. The maximal signal rates  Although the N2HDM still provides larger signal rates than the C2HDM the differences are not as significant as in T1. Since we do not have parameter points with the mass relation 2m H ↓ m H ↑ in the C2HDM T2, we cannot discuss the comparison of the decay channel H ↑ → H ↓ H ↓ in the T2. Furthermore, the rates for the decay channel H ↑ → H ↓ h are not sufficiently enhanced to be interesting for di-Higgs searches so that we do not discuss these final states further here. The rates for the (bb)(γγ) final states are below 1 fb in both models so that we also do not consider these channels here.
We conclude by remarking that the signal σ pp → h → H ↓ → bb H ↓ → bb in the N2HDM T2 can reach max σ pp → h → H ↓ → bb H ↓ → bb N2HDM = 704 fb . (5.99) To enable such decays a semi-inverted mass hierarchy is required where m h 2m H ↓ . (The inverted mass hierarchy is not realised in the N2HDM T2.) While the signal rate is large, it has to be kept in mind, that this mass scenario is already under pressure to fulfill all required constraints.

Conclusions
In this paper, we investigated the possibility of an SFOEWPT in models with non-minimal Higgs sectors, the C2HDM and the N2HDM. For the C2HDM we updated our analysis [27] by allowing for heavier neutral scalar masses up to 1.5 TeV and including the most recent collider constraints. Still similar mass regions compatible with theoretical and experimental constraints and the requirement of an SFOEWPT were found. The inclusion of the new constraints from the Higgs data significantly reduced, however, the scenarios compatible with an SFOEWPT. The strength of the phase transition ξ c in both types of the C2HDM is rather small (max ξ C2HDM,T1 c = 1.7 , max ξ C2HDM,T2 c = 1.18). In the N2HDM the inclusion of the additional real scalar singlet slightly enhances the strength of the phase transition (max ξ N2HDM,T1 c = 2.1 , max ξ N2HDM,T2 c = 1.5) compared to the C2HDM. The compatibility with the EW precision constraints requires two of the Higgs boson masses in the spectrum to be close. The additional requirement of an SFOEWPT reduces this mass gap further. For the N2HDM, we found that not all mass hierarchies were compatible both with the theoretical and experimental constraints and an SFOEWPT. Thus, the inverted mass hierarchy (m H ↓ < m H ↑ < m h ) was not found in the N2HDM T1 and T2. In the N2HDM T2 the interplay of the requirement of small mass gaps and a heavy charged Higgs boson mass yields an overall heavier Higgs boson mass spectrum than in the N2HDM T1. We showed that the SFOEWPT favors two different mass configurations in the N2HDM T1. The first region M deg features parameter points with m A ≈ m H ± , in the second one, M sep , the pseudoscalar Higgs boson A is significantly heavier than the charged Higgs boson H ± . In particular the mass region M sep is often neglected in the context of electroweak baryogenesis. For simplicity, it is often assumed that m A ≈ m H ± . In the N2HDM T2 only points in M deg were found to be compatible with an SFOEWPT.
We furthermore investigated in both models the trilinear Higgs self-couplings between neutral Higgs bosons. An SFOEWPT favors enhanced trilinear self-couplings. Their values remain, however, significantly below the values allowed by the experimental and theoretical constraints. For the Type I of both models they are significantly constrained by the recent updates in the di-Higgs searches. Since the SFOEWPT favors Higgs boson masses in the range ∼ 400 − 600 GeV, which is also the most sensitive region in the di-Higgs searches, a large part of our previous sample of parameter points providing an SFOEWPT is excluded. Future updates in the di-Higgs searches will therefore have a significant impact on the valid parameter space of the N2HDM and C2HDM with respect to an SFOEWPT. On the other hand, Type II of both models is strongly constrained due to the results for B → X s γ [64][65][66][67][68]. We observed that the charged Higgs mass tends to smaller mass regions if an SFOEWPT is required. This cannot be realized in Type II as the charged Higgs mass must be above 580 GeV because of the B → X s γ constraints. Most of the valid parameter points in the C2HDM T2 in Fig. 1 and in the N2HDM T2 in Fig. 10 are located on the left edge in the m H ± − tan β plane. Exactly this region is sensitive to the B → X s γ constraints. We therefore expect the strongest constraints for the Type II to come from future updates on the charged Higgs mass constraint in the flavor sector. We finally investigated the phenomenology of both models with respect to the expected di-Higgs signals at the LHC. Requiring an SFOEWPT, the expected di-Higgs signals in the C2HDM are significantly suppressed compared to the N2HDM.