The CP-Violating 2HDM in Light of a Strong First Order Electroweak Phase Transition and Implications for Higgs Pair Production

We investigate the strength of the electroweak phase transition (EWPT) within the CP-violating 2-Higgs-Doublet Model (C2HDM). By applying a renormalisation scheme which allows efficient scans of the C2HDM parameter space, we analyse the possibility of a strong first order EWPT required for baryogenesis and study its phenomenological implications for the LHC. Like in the CP-conserving (real) 2HDM (R2HDM) we find that a strong EWPT favours mass gaps between the non-SM-like Higgs bosons. These lead to prominent final states comprised of gauge+Higgs bosons or pairs of Higgs bosons. In contrast to the R2HDM, the CP-mixing of the C2HDM also favours approximately mass degenerate spectra with dominant decays into SM particles. The requirement of a strong EWPT further allows us to distinguish the C2HDM from the R2HDM using the signal strengths of the SM-like Higgs boson. We additionally find that a strong EWPT requires an enhancement of the SM-like trilinear Higgs coupling at next-to-leading order (NLO) by up to a factor of 2.4 compared to the NLO SM coupling, establishing another link between cosmology and collider phenomenology. We provide several C2HDM benchmark scenarios compatible with a strong EWPT and all experimental and theoretical constraints. We include the dominant branching ratios of the non-SM-like Higgs bosons as well as the Higgs pair production cross section of the SM-like Higgs boson for every benchmark point. The pair production cross sections can be substantially enhanced compared to the SM and could be observable at the high-luminosity LHC, allowing access to the trilinear Higgs couplings.


Introduction
The discovery of the Higgs boson by the LHC experiments ATLAS [1] and CMS [2] has marked a milestone for particle physics while, at the same time, leaving many open questions. Despite the Standard Model (SM) nature of the Higgs boson [3][4][5][6] new physics (NP) beyond the SM is called for in order to solve the puzzles within the standard theory. The observed baryon asymmetry of the Universe (BAU) [7] is one example that requires NP extensions. It can be generated dynamically in the early Universe during a first order electroweak phase transition (EWPT) through the mechanism of electroweak baryogenesis (EWBG) [8][9][10][11][12][13][14][15][16] provided that all three Sakharov conditions [17] are fulfilled. These are baryon number violation, C and CP violation and departure from the thermal equilibrium. A strong first order phase transition (PT) [14,16] proceeds through bubble nucleation and can generate a net baryon number near the bubble wall. Diffusion through the bubble preserves this asymmetry for the later evolution of the universe by suppressing the sphaleron transitions in the false vacuum [18,19]. A first order EWPT takes places when bubbles of the broken phase nucleate in the surrounding plasma of the symmetric phase. The bubble walls interact with the various fermion species in the plasma. If there is CP violation in the bubble wall or some CP violation in the hot i.e. the symmetric phase that is disturbed by the wall, particles with opposite chirality interact differently with the wall and CP and C asymmetries in particle number densities can be generated in front of the bubble wall. These asymmetries diffuse into the hot plasma ahead of the bubble wall biasing baryon number violating electroweak (EW) sphaleron transitions to generate a baryon asymmetry. The latter is transferred through the expanding wall into the broken phase. The rate of sphaleron transitions is strongly suppressed in this phase such that a washout of the baryons generated before is avoided. The required departure from thermodynamic equilibrium is guaranteed by the passage of the bubble walls that rapidly expand through the cosmological plasma. Additionally, gravitational waves produced by a strong first order EWPT [20] are potentially observable by the future space-based gravitational wave interferometer eLISA [21]. The interplay between gravitational waves and a strong first order PT and/or collider phenomenology has recently been studied in [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].
In the SM in principle all three Sakharov conditions could be realized. However, the EWPT is not of strong first order [41], since this would require a SM Higgs boson mass of around 70-80 GeV [42]. Additionally, the CP violation of the SM arising from the Cabibbo-Kobayashi-Maskawa (CKM) matrix is too small [16,[43][44][45]. Extensions beyond the SM provide additional sources of CP violation as well as further scalar states triggering a first order EWPT also for a SM-like Higgs boson with a mass of 125 GeV. This is the case for the 2-Higgs Doublet Model (2HDM) [46,47] which belongs to the simplest NP extensions that are compatible with present experimental constraints. Previous studies have shown that 2HDMs provide a good framework for successful baryogenesis, both in the CP-conserving [48][49][50][51][52] and in the CP-violating case [29,[53][54][55][56][57][58][59]. 2HDMs feature five physical Higgs bosons, and the tree-level 2HDM Higgs sector provides sources for explicit CP violation. Allowing for CP violation, there could in principle be a complex phase between the vacuum expectation values (VEVs) of the two Higgs doublets. This phase can, however, always be removed by a change of basis [60] so that without loss of generality it can be set to zero. This may not be the case any more at finite temperature. CP violation might be generated spontaneously only in the bubble wall around the critical temperature [54] and provide a source for the generation of the matter-antimatter asymmetry through EWBG. The CP-violating phase which at zero temperature is just another parameter of the theory, during EWPT becomes a spatially varying field, and its value depends on the position relative to the bubble walls. In order to study the effect of CP violation on the generation of the baryon asymmetry the detailed form of the spatially varying field has to be determined.
CP-violating Higgs sectors are strongly constrained by the electric dipole moments (EDMs). The strongest constraint [61] is imposed by the limit on the electron EDM provided by the ACME collaboration [62]. The possibility of spontaneous CP violation generated at the EWPT which vanishes at zero temperature may provide an attractive scenario to lift the possible tension between the restrictions imposed by the EDMs and the requirement of a substantial amount of CP violation by the EWBG.
In [51] we investigated the implications of a strong first order PT in the CP-conserving or real 2HDM (R2HDM) on the LHC Higgs phenomenology. We found a strong interplay between the requirement of successful baryogenesis and LHC Higgs phenomenology. In this work, we extend our analysis to the CP-violating 2HDM (C2HDM). 1 The computation of the equation of motion for the CP-violating phase between the two Higgs doublets and the computation of the actual baryonantibaryon asymmetry generated through EWBG within the framework of the C2HDM is beyond the scope of this paper. We focus instead on the interplay between the requirement of a strong first order phase transition and LHC phenomenology in the presence of explicit CP violation in the tree-level 2HDM Higgs sector. We investigate the possible spontaneous generation of a CP-violating phase at the EWPT. We furthermore analyse in detail the effect of higher order corrections on the trilinear Higgs self-couplings extracted from the one-loop corrected effective potential. 2 We discuss the impact of the requirement of a strong phase transition on their size and the resulting implications for LHC phenomenology, namely Higgs pair production. We present several benchmark scenarios, emphasizing the specific features of C2HDM parameter points compatible with all constraints and a strong phase transition.
For the purpose of this paper we compute the one-loop corrected effective potential at finite temperature [69][70][71] including daisy resummations for the bosonic masses [72]. For the numerical analysis the parameter space of the C2HDM is scanned and tested for compatibility of the model with the theoretical and experimental constraints. Subsequently, the implication of a strong first order PT on the surviving parameter sets is determined and interpreted with respect to collider phenomenology. The former necessitates the minimisation of the loop-corrected Higgs potential at increasing temperature in order to find the vacuum expectation value v c at the critical temperature T c , which is defined as the temperature where two degenerate global minima exist. A value of v c /T c larger than one is indicative of a strong first order PT [11,73]. 3 In order to be able to perform an efficient scan, like in [51], we renormalise the loop-corrected potential in such a way that not only the VEV and all physical Higgs boson masses, but also all mixing matrix elements remain at their tree-level values. In our analysis, we will focus on the C2HDM with type I and type II couplings of the Higgs doublets to the fermions. We will discard parameter points inducing a 2-stage PT [77,78]. Our analysis reveals a strong link between the demand for a strong first order PT and testable implications at the collider experiments.
The outline of the paper is as follows: In section 2 we set our notation and present the loopcorrected effective potential of the C2HDM at finite temperature. Our renormalisation procedure is described in section 3. Section 4 is dedicated to the description of the numerical analysis. It includes the outline of the minimisation procedure of the effective potential and the details of the scan in the C2HDM parameter space as well as of the applied theoretical and experimental constraints. Sections 5-10 contain our results. In Section 6 we analyse the type I C2HDM with the lightest Higgs boson being the SM-like Higgs state. We first investigate the spontaneous generation of a CP-violating phase and its relation to explicit CP violation in the tree-level potential. We then present the parameter regions compatible with the applied constraints and a strong first order PT and analyse the implications for collider phenomenology. In Section 7 we discuss in detail the role of the trilinear Higgs self-couplings in the EWPT, the impact of the next-to-leading order (NLO) corrections derived from the effective potential and the implications of the requirement of a strong phase transition on the Higgs self-couplings and their corrections. In Section 8 we briefly summarise the results for the type I C2HDM with the next-to-lightest Higgs boson representing the SM-like Higgs scalar. Sections 9 and 10 are dedicated to the type II C2HDM with the lightest Higgs boson being SM-like, and we present our results in analogy to the type I case. Section 11 contains our conclusions.

The effective potential in the C2HDM
In this section we provide the loop-corrected effective potential at finite temperature for the CPviolating 2HDM. We start by setting our notation.

The CP-violating 2-Higgs-Doublet Model
The tree-level potential of the C2HDM for the two SU (2) L scalar doublets reads It incorporates a softly broken Z 2 symmetry, under which the doublets transform as This ensures the absence of tree-level Flavor Changing Neutral Currents (FCNC). The hermiticity of the potential V tree forces all parameters to be real apart from the soft Z 2 breaking mass parameter m 2 12 and the quartic coupling λ 5 . For arg(m 2 12 ) = arg(λ 5 ) the complex phases of these two parameters can be absorbed by a basis transformation. If furthermore the VEVs of the two doublets are assumed to be real, we are in the real or CP-conserving 2HDM. Otherwise, we are in the C2HDM, for which we will adopt the conventions of [79] in the following.
After EW symmetry breaking the two Higgs doublets acquire VEVsω i ∈ R, about which the Higgs fields can be expanded in terms of the charged CP-even and CP-odd fields ρ i and η i , and the neutral CP-even and CP-odd fields ζ i and ψ i , i = 1, 2. In the general 2HDM, there are three different types of minima, given by the normal EW-breaking one, a CP-breaking minimum, and a charge-breaking (CB) vacuum. In Refs. [80][81][82] it has been shown that at tree level minima that break different symmetries cannot coexist. If a normal minimum exists, all CP or CB stationary points are proven to be saddle points. These statements may not be true any more at higher orders, as recent studies have shown for the Inert 2HDM at one-loop level in the effective potential approach [83]. Consequently, we allow for the possibility of a CP-breaking vacuum as well as a charge-breaking one. Through the VEVsω CP andω CB we include the possibility of generating at one-loop and/or non-zero temperature a global minimum that is CP-violating and/or charge breaking. As a charge-breaking VEV breaks electrical charge conservation inducing a massive photon, this unphysical configuration of the vacuum will not further be discussed in the numerical analysis. Denoting the VEVs of the normal vacuum byω 1,2 and the CP-and charge-breaking VEVs byω CP andω CB , respectively, the expansion of the two Higgs doublets about the VEVs is given by where, without loss of generality, the complex part of the VEVs and the charge-breaking VEV have been rotated to the second doublet exclusively. The VEVs of our present vacuum at zero temperature 4 are denoted as with The VEVs of the normal vacuum, v 1,2 , are related to the SM VEV v ≈ 246 GeV by The angle β defines the ratio of v 1 and v 2 , so that The minimum conditions of the potential Eq. (2), where the brackets denote the Higgs field values in the minimum, i.e.
where we have introduced the abbreviation Equations (11a) and (11b) can be used to trade the parameters m 2 11 and m 2 22 for v 1 and v 2 , while Eq. (11c) leads to a relation between the two sources of CP violation in the scalar potential so that one of the ten parameters of the C2HDM is fixed. Introducing the neutral mass eigenstates H i (i = 1, 2, 3) are obtained from the C2HDM basis ζ 1 , ζ 2 and ζ 3 through the rotation The corresponding Higgs masses are obtained from the mass matrix through the diagonalisation with the matrix R, The Higgs bosons are ordered by ascending mass as the mixing matrix R can be parametrised as Type I Type II Lepton-Specific Flipped Exploiting the minimum conditions of the potential at zero temperature, we use the following set of 9 independent parameters of the C2HDM [84], The m H i and m H j denote any two among the three neutral Higgs boson masses, and the mass of the third Higgs boson is obtained from the other parameters [84]. For the analytic relations between the above parameter set and the coupling parameters λ i of the 2HDM Higgs potential, see [79].
The limit of the CP-conserving 2HDM is obtained for α 2 = α 3 = 0 and α 1 = α + π/2 [85]. The mass matrix Eq. (15) becomes block diagonal in this case leading to the pure pseudoscalar A, which is identified with H 3 , and the CP-even mass eigenstates h and H, which are obtained from the gauge eigenstates ζ 1,2 through the rotation with the angle α, The imposed Z 2 symmetry ensures that each of the up-type quarks, down-type quarks and charged leptons can only couple to one of the Higgs doublets so that FCNCs at tree level are avoided. Table 1 lists the possible different 2HDM types given by type I, type II, lepton-specific and flipped. For the exact form of the C2HDM Higgs couplings to the SM particles in terms of the input parameters, we refer to Refs. [86][87][88]. In this work we focus on C2HDMs of type I and type II.

One-loop effective potential at finite temperature
The form of the one-loop effective potential at finite temperature for the C2HDM case does not change with respect to the one introduced for the CP-conserving 2HDM in Ref. [51]. For convenience of the reader, we briefly repeat the main ingredients, also in order to set our notation.
The one-loop contribution V 1 to the effective potential consists of the Coleman-Weinberg (CW) contribution V CW [69] already present at zero temperature, and the contribution V T for the thermal corrections at finite temperature T . The one-loop corrected effective potential reads with the tree-level potential given in Eq. (2) after replacing the doublets Φ 1,2 with their classical constant field configuration In the MS scheme the Coleman-Weinberg potential for a particle i reads [71] The sum extends over the Higgs and Goldstone bosons, the massive gauge bosons, the longitudinal photon and the fermions f , with the exception of the neutrinos, which we assume to be massless, i = h, H, A, H ± , G 0 , G ± , W ± , Z, γ L , f . Here, m 2 i denotes the respective eigenvalue for the particle i of the mass matrix squared expressed through the tree-level relations in terms of ω i (i = 1, 2, CP, CB). The sum also extends over the Goldstone bosons and the photon. Although in the Landau gauge applied here the Goldstone bosons are massless at T = 0, they can become massive for field configurations different from the tree-level VEVs at T = 0, which are required in the minimisation procedure. This is also the case for the photon, as we allow for non-physical vacuum configurations with a non-zero charge breaking VEV. Furthermore, the Goldstone bosons and the longitudinal photon can become massive due to the temperature corrections discussed below. Because of the Landau gauge we need not consider any ghost contributions. The spin of the particle is denoted by s i and the number of degrees of freedom by n i . For the scalars Φ = H i , G 0 , H + , H − , G + , G − , the charged leptons 5 l + and l − , the quarks and antiquarks q andq, the longitudinal and transversal gauge bosons V L = Z L , W + L , W − L , γ L and V T = Z T , W + T , W − T , γ T , they are In the MS scheme employed here the constants c i read The renormalisation scale µ is fixed to µ = v = 246.22 GeV.
The thermal corrections V T comprise the daisy resummation [72] of the n = 0 Matsubara modes of the longitudinal components of the gauge bosons W + L , W − L , Z L , γ L and the bosons Φ, so that their masses receive Debye corrections at non-zero temperature. The potential V T can be cast into the form [70,71] Since the Goldstone bosons and the photon acquire a mass at finite temperature, they have to be included in the sum. Denoting the 5 Because of the CB-breaking VEV we have to take into account different masses for the charge conjugated particles. mass eigenvalue including the thermal corrections for the particle k by m k , we have for J (k) ± (see e.g. [89]) with the thermal integrals where J + (J − ) applies for k being a fermion (boson). The masses m i depend implicitly on the temperature T , since for each T we determine the VEVs, respectively the field configurations, ω i =ω i (T ), that minimise the loop-corrected potential V , Eq. (21). These field configurations enter the tree-level mass matrices. The m k in addition depend explicitly on T through the thermal corrections. With the definition of J (k) ± Eq. (27) we follow the 'Arnold-Espinosa' approach of Ref. [90]. A different approach has been proposed in [91], to which we refer as 'Parwani' method.
Here the Debye corrections are included for all the bosonic thermal loop contributions and the Debye corrected masses are also used in the CW potential. Since the 'Parwani' method admixes higher-order contributions, possibly leading to dangerous artifacts at one-loop level, we apply the 'Arnold-Espinosa' method. For a discussion and comparison of the two methods, see also [56,58].
with Γ(x) denoting the Euler Gamma function. For the interpolation between the two approximations the point is determined where the derivatives of the low-and high-temperature expansions can be connected continuously. At this point a small finite shift to the small x 2 expansion is added such that also the two expansions themselves are connected continuously. Denoting the values of x 2 where the connection is performed, by x 2 + and x 2 − and the corresponding shifts by δ ± for the fermionic and bosonic contributions, respectively, they are given by For small x 2 the exact result is well approximated by including terms of up to order n = 4 in the expansion J +,s for fermions, for bosons this is the case for n = 3 in J −,s . For large x 2 , the integral is well approximated by n = 3 in both the fermion and the boson case, J ±,l . The deviation of the approximate results from the numerical evaluation of the integrals is less than two percent.  28), is set to the real part of its numerical evaluation which is the relevant contribution for the extraction of the global minimum [92]. In practice, the integral is evaluated numerically at several equidistant points in m 2 /T 2 < 0, and in the minimisation procedure the result obtained from the linear interpolation between these points is used. This allows for a significant speed-up. We explicitly verified that the difference between the exact and the interpolated result is negligible for a sufficiently large range of m 2 /T 2 .

Renormalisation
The masses and mixing angles extracted from the loop-corrected potential differ from those extracted from the tree-level potential. In the tests for the compatibility of the model with the experimental constraints these corrections have to be taken into account. In order to ensure an efficient scan over the parameter space of the model in terms of the input parameters Eq. (19), it is more convenient to directly use the loop-corrected masses and angles as inputs. We therefore modify the MS renormalisation applied in the Coleman-Weinberg potential Eq. (23) and choose a renormalisation prescription by which we enforce the one-loop corrected masses and mixing matrix elements to be equal to the tree-level ones. This follows the approach chosen in our analysis of the CP-conserving 2HDM [51]. The counterterm potential V CT , which is added to the one-loop effective potential Eq. (21), reads In the last line we explicitly included the tadpole counterterms δT for the directions in field space in which we allow for the development of a vacuum. 6 Since we check for the compatibility with the experimental constraints at T = 0 we apply our renormalisation conditions at this temperature. They are given by (i, j = 1, ..., 8) with and φ c T =0 denoting the field configuration in the minimum at T = 0, The first set of conditions, Eq. (36), ensures that at T = 0 the tree-level position of the minimum yields a local minimum. We check numerically if it is also the global one. The second set of conditions, Eq. (37), guarantees that at T = 0 both the masses and the mixing angles remain at their tree-level values. In Ref. [93] formulae for both the first and the second derivatives of the CW potential have been derived in the Landau gauge basis. We employ these formulae to calculate the required derivatives. Since the system of equations resulting from the conditions Eqs. (36) and (37) is not sufficient to fix all renormalisation constants, one of them is left free. In analogy to our previous paper [51], we choose to set This finally yields the counterterms in terms of the derivatives of the CW potential, Since the physical vacuum is required to be neutral, there exist no charge breaking diagrams contributing to the one-loop effective potential. with We note that the second derivative of the CW potential, required for our renormalisation procedure, leads to the well-known problem of infrared divergences for the Goldstone bosons in the Landau gauge [49,56,58,[93][94][95][96]. For the procedure on how to treat this problem, we refer to the investigation within the CP-conserving 2HDM in Ref. [51]. We checked that by applying these formulae in the limit of the real 2HDM we reproduce our results of the R2HDM.

Minimisation of the Effective Potential
The electroweak PT is of strong first order if the ratio between the VEV v c acquired at the critical temperature T c , and the critical temperature T c is larger than one [11,73], The value v at a given temperature T is given by whereω i are the field configurations that minimise the loop-corrected effective potential at non-zero temperature. The critical temperature T c is defined as the temperature where the potential has two degenerate minima. In order to obtain T c , the complete loop-corrected effective potential Eq. (34), is minimised numerically for a given temperature T . If the PT is of strong first order, the VEV jumps from v = v c at the temperature T c to v = 0 for T > T c . For the determination of T c we employ a bisection method in the temperature T , starting with the determination of the minimum at the temperatures T S = 0 GeV and ending at T E = 300 GeV. The minimisation procedure is terminated when the interval containing T c is smaller than 10 −2 GeV, and the temperature T c is then set to the lower bound of the final interval. Parameter points that do not satisfy |v(T = 0) − 246.22 GeV| ≤ 2 GeV are excluded as well as parameter points where no PT is found for T ≤ 300 GeV. Adding a small safety margin, we ensure by the latter condition that possible strong first order PTs are obtained for VEVs below 246 GeV. We furthermore only retain parameter points with T c > 10 GeV.

Constraints and Parameter Scan
The points, for which the value of ξ c is determined, have to satisfy theoretical and experimental constraints. We use a pre-release version of ScannerS [97,98] to perform scans in the C2HDM parameter space in order to obtain viable data sets. In these extensive scans we check for compatibility with the following constraints. The potential is required to be bounded from below and the treelevel discriminant of Ref. [99] is used to enforce that the electroweak vacuum is the global minimum of the tree-level potential at zero temperature. We further require perturbative unitarity to hold at tree level. We take into account the flavour constraints on R b [100,101] and B → X s γ [101][102][103][104][105]. They can be generalized from the CP-conserving 2HDM to the C2HDM as they only depend on the charged Higgs boson. These constraints are checked as 2σ exclusion bounds on the m H ± − t β plane. According to the latest calculation of Ref. [105] the charged Higgs boson mass is required to be rather heavy, in the type II and flipped 2HDM. In the type I and lepton-specific model this bound is much weaker and depends more strongly on tan β. Agreement with the electroweak precision measurements is verified using the oblique parameters S, T and U . The formulae for their computation in the general 2HDM can be found in [47]. For the computed S, T and U values 2σ compatibility with the SM fit [106] is demanded including the full correlation among the three parameters. One of the Higgs bosons, called h in the following, is required to have a mass of [107] m h = 125.09 GeV .
Compatibility with the Higgs data is checked by using HiggsBounds [108] and the individual signal strength fits of Ref. [109] for the h. The required decay widths and branching ratios are obtained from a private implementation of the C2HDM into HDECAY v6.51 [110,111], which will be released in a future publication. Additionally, the Higgs boson production cross sections normalized to the SM are needed, including the most important state-of-the-art higher order corrections. Where available, we include the QCD corrections which can be taken over from the SM and Minimal Supersymmetric Extension of the SM (MSSM). Electroweak corrections are consistently neglected in both the production and decay channels, as they cannot be taken over and are not available yet for the C2HDM. Details on how the production cross sections are determined can be found in [87]. This information is passed via the ScannerS interface to HiggsBounds which checks for agreement with all 2σ exclusion limits from LEP, Tevatron and LHC Higgs searches. As mentioned above, the properties of the h are checked against the fitted values of the signal strengths given in [109]. For details, we again refer to [87]. We use this method for simplicity. Note that performing a fit to current Higgs data is likely to give a stronger bound than this approach. As we include CP violation in the Higgs sector we also have to check for compatibility with the measurements of electric dipole moments (EDM), where the strongest constraint originates from the EDM of the electron [61]. The experimental limit has been given by the ACME collaboration [62]. For the check, we have implemented the calculation of the dominant Barr-Zee contributions by [112] and require compatibility with the bound given in [62] at 90% C.L.
For the scan, the SM VEV is fixed to The ranges chosen for the remaining input parameters of Eq. (19) are as follows. The mixing angle t β has been varied as The angles parametrising the mixing matrix Eq. (18) are chosen in the intervals For Re(m 2 12 ) we use the range 0 GeV 2 ≤ Re(m 2 12 ) < 500 000 GeV 2 .
Note that, although possible, physical parameter points with Re(m 2 12 ) < 0 are extremely rare so that we neglect them in our study. This is mainly a result of requiring absolute stability at tree level. One of the neutral Higgs bosons H i is identified with h. In type II, the charged Higgs mass is chosen in the range and in type I in the range The electroweak precision constraints combined with perturbative unitarity require at least one neutral Higgs boson to be close in mass to m H ± . For an increased scan efficiency we therefore choose the second neutral Higgs mass m H i =h in the interval 500 GeV ≤ m H i < 1 TeV (54) in type II and in type I. In the C2HDM the third neutral Higgs boson m H j =H i ,h is not an independent input parameter and is calculated by ScannerS. It is, however, required to lie in the interval We further impose that the m H i =h deviate by at least 5 GeV from 125.09 GeV to avoid degenerate Higgs signals. To improve the coverage of the CP-conserving limit we have performed dedicated scans in the CP-conserving 2HDM and merged the resulting CP-violating and CP-conserving samples. The scans in the CP-conserving 2HDM were also performed using ScannerS with the same constraints 7 and parameter ranges. The final samples are composed of more than 2 · 10 5 valid parameter points for each Yukawa type.
For the SM parameters we have chosen the following values: Apart from the computation of the oblique parameters, where we use the fine structure constant at zero momentum transfer, 7 With the exception of the EDM constraint which is trivially satisfied if CP is conserved.
the fine structure constant is taken at the Z boson mass scale [113], The massive gauge boson masses are set to [113,114] M W = 80.385 GeV and the lepton masses to [113,114] and the light quark masses to following [115]. For consistency with the ATLAS and CMS analyses the on-shell top quark mass m t = 172.5 GeV (62) has been taken, as recommended by the LHC Higgs Cross Section Working Group (HXSWG) [114,116]. The charm and bottom quark on-shell masses are [114] m c = 1.51 GeV and We take the CKM matrix to be real, with the CKM matrix elements given by [113] 8

Results
In our analysis we investigate the question to which extent the allowed parameter space of the C2HDM is constrained by the requirement of a first order phase transition and what are the consequences for LHC phenomenology. We compare with the case of the CP-conserving 2HDM which has been analysed in [51]. We investigate the impact on the trilinear Higgs self-couplings and Higgs pair production. 9 We analyse the size of the electroweak corrections to the Higgs selfcouplings derived from the effective potential. We furthermore study the possible spontaneous generation and size of a CP-violating phase at the electroweak phase transition. We will show results both for the type I and the type II C2HDM and for the cases where the lightest of the neutral Higgs bosons is identified with the discovered Higgs boson, i.e. H 1 ≡ h, and where the next heavier one is the SM-like Higgs boson, H 2 ≡ h.
For the interpretation of the results, we note that the strength of the phase transition increases with the size of the couplings of the light bosonic particles to the SM-like Higgs boson and decreases 8 In the computation of the loop-corrected effective potential we choose VCKM = 1 for simplicity. The impact of this choice on the counterterms and thereby on the potential and its minimisation is negligible. 9 For previous studies on Higgs pair production in the real 2HDM, see [117,118].
with the Higgs boson mass [89]. Since in the C2HDM all non-SM-like neutral Higgs bosons receive a VEV through mixing and hence contribute to the PT, a strong electroweak PT requires the participating Higgs bosons either to be light or to have a VEV close to zero. In the latter case we are in the alignment limit where only one of the physical Higgs bosons has a VEV [119]. In the type II C2HDM the requirement of a light Higgs spectrum puts the model under tension, as EW precision tests combined with perturbative unitarity enforce one of the neutral Higgs bosons to be close to m H ± . Charged Higgs masses below 580 GeV are already excluded by B → X s γ, however.
Note, that in contrast to the analysis of the CP-conserving 2HDM in [51], in the C2HDM we have a larger number of parameters to be scanned over. Also the minimisation procedure requires more computing power due to a possible CP-violating VEV. The result is, that the overall density of parameter points compatible with our applied constraints is smaller than in the real 2HDM. Consequently, we found for type I only very few parameter points where H 2 = h and where we have both a strong phase transition and CP violation in the Higgs sector. A considerably enlarged parameter scan might lead to more points fulfilling these criteria. Here we content ourselves to demonstrate that such configurations are possible in principle. In the type II C2HDM, no parameter sets were found where the SM-like Higgs boson is given by the heavier neutral Higgs bosons, due to the constraint m H ± ≥ 580 GeV [105]. 10 Validity of the global minimum and of the unitarity constraint at NLO Since we compute the global minima of the effective potential at NLO, an interesting question to ask is how the inclusion of NLO effects influences the absolute stability of the EW vacuum. In the scan of the type I C2HDM with H 1 ≡ h we found that the inclusion of the NLO computation eliminated about 6% of the parameter points due to the EW vacuum no longer being the global minimum at NLO. In case H 2 ≡ h, 26% of the tree-level points did not lead to a global minimum any more. For the type II scan we found that the requirement of a global minimum at NLO eliminated 9% of the points with a valid tree-level global minimum in case H 1 ≡ h. If, however, the heavier H 2 is the SM-like Higgs boson, there are practically no scenarios that represent a valid minimum of the potential at NLO. It turns out that the request of a global NLO minimum at 246 GeV implies reduced mass differences between the different Higgs bosons. If, however, the H 2 is SM-like then the mass difference between the lightest Higgs boson H 1 with m H 1 < 125 GeV and the charged Higgs boson with a required mass above 580 GeV is too large to allow for an NLO minimum at 246 GeV. Consequently, there are also no scenarios with a strong first order PT in this case. In fact, as can be read off the formula Eq. (23) for the Coleman-Weinberg potential large masses imply large one-loop corrections. Through our renormalisation procedure we move these large corrections into the quartic couplings, so that they may become too large to guarantee a stable vacuum.
This discussion shows that the inclusion of the NLO effects is important in order to correctly define the parameter regions that are compatible with the requirement that the EW minimum 10 This mass configuration corresponds in the limit of the real 2HDM to the cases where the lighter neutral Higgs boson h corresponds to the 125 GeV Higgs boson and mA < m h = 125 GeV or the heavier one, H, represents the discovered Higgs boson and m h < mH = 125 GeV. Already in the R2HDM where we applied the older constraint of m H ± ≥ 480 GeV, we found very few scenarios in this case that are compatible with a strong PT, cf. Figs. 5 and 11 in [51]. With the stricter lower limit of 580 GeV on the charged Higgs mass we do not find any allowed scenarios with a strong PT any more. In the CP-violating case where in general we find less scenarios compatible with a strong PT the situation becomes even more severe.
represents the global minimum. In our analysis we only keep points compatible with a global NLO minimum at 246 GeV. Additionally, we have to make sure that the possibly large corrections to the Higgs self-couplings due to our renormalisation procedure do not spoil the unitarity constraint. In order to check this, in a first rough approximation we insert the renormalised λ i , derived from the loop-corrected effective potential, into the tree-level formulae for the unitarity constraints [47]. We keep only those points for which the bound of 8π is not violated. The unitarity check further reduces the sample of points fulfilling the experimental and global minimum constraints by 11% and 18% in the type I C2HDM with H 1 ≡ h and H 2 ≡ h, respectively, and by 9% in the type II C2HDM with H 1 ≡ h.

The CP-violating phase
We start by investigating the size of a possible CP-violating phase that is spontaneously generated at the EWPT, and its relation to the explicit CP violation through a complex phase of m 2 12 . In Fig. 1, we depict the tangent of the CP-violating phase tan ϕ spont =ω 3 (T c )/ω 2 (T c ) at the critical temperature T c as a function of the tangent of the CP-violating phase at zero temperature 11 tan ϕ explicit = Im(m 2 12 )/Re(m 2 12 ) for the points that fulfill all constraints 12 and are compatible with a strong PT. The colour code indicates the value of ξ c . Note that in Fig. 1 we only plot C2HDM parameter points, i.e. Im(m 2 12 )/Re(m 2 12 ) = 0, although this ratio can become very small. As can be inferred from the figure, the phase ϕ spont of spontaneous CP violation and the one of explicit CP violation, ϕ explicit , are correlated. A CP-violating phase at T c is generated spontaneously only if already in the zero-temperature potential there is non-vanishing CP violation. As we set the CKM matrix to unity in the computation of the effective potential no CP violation can be generated through loop effects if CP is conserved at T = 0. The size of ξ c is not correlated to the size of the CP-violating phase. We observe, however, that the maximum obtained value of ξ c , which quantifies the strength of the PT, is ξ c = 1.89. This is below the value found in the CP-conserving 2HDM, cf. Ref. [51] and the following section. (65) 11 In the C2HDM, the only CP-violating source in the Higgs potential is given by a complex m 2 12 or alternatively a complex λ5 which is related to m 2 12 through the tadpole conditions. Any other CP-violating phase, namely the phase of the VEV, can be absorbed by a redefinition of the fermion fields. 12 Here and in all following plots this means that they fulfill the experimental constraints and also the unitarity and global minimum constraints at NLO. We observe, in the upper left plot the generation of the VEV v = ω 2 1 +ω 2 2 +ω 2 3 (ω 3 ≡ω CP ) at T c = 136.011 GeV with a value of v c = 156.245 GeV and hence ξ c = 1.149. With decreasing temperature the VEV increases to the value 246 GeV at T = 0. Note, that the VEV also includes ω CB , which, however, always turns out to be zero, so that we did not write it explicitly here. The upper right and lower left plots show the development of the absolute values of the individual VEVs, i.e. the ones ofω 1 andω 2 , the CP-conserving VEVs coinciding with v 1 and v 2 at zero temperature, and of the CP-violating VEVω 3 . We also show, in the lower right plot, the development ofω 3 /ω 2 . Because of the value of tan β above 1, the value |ω 2 | is larger than |ω 1 |. The spontaneously generated CP-violating VEVω 3 at the PT amounts with |ω 3 (T c )| = 17.06 to 11% of the absolute value ofω 2 and decreases monotonously to zero with decreasing temperature, whileω 1 andω 2 monotonously increase to reach ω 2 1 +ω 2 2 (T = 0) = 246 GeV. Figure 3 shows the total CP-violating angle at the PT as a function of the CP-violating angle at T = 0. The former varies between about -34 • and 49 • for CP-violating angles of -35 • to 51 • allowed at T = 0. We observe almost maximal CP violation which should be enough for successful baryogenesis [10,43]. Figure 4 shows the mass of the heavier non-SM-like Higgs boson, m ↑ , versus the lighter one, m ↓ , where the grey points pass the applied constraints and include both CP-conserving and CPviolating points in the left plot, but only CP-violating points in the right plot. The coloured points additionally feature a strong first order PT. The maximum possible value of ξ c is found to be ξ c = 5.7 for all 2HDM points (left plot). 13 If we only consider CP-violating points (right plot) the 13 Barring the few CP-violating points, the left plot can be compared to the results of Ref. [51]. There, we found ξ max c = 4.5. The difference to ξ max c found here (and also the difference in the shape of this plot and all following ones, that can be used for a comparison) arises from a different constraint on tan β due to different applied flavour   of a strong PT overall prefers somewhat lighter non-SM-like Higgs bosons. The allowed maximum masses are further reduced in case of CP violation where m ↑ remains below about 753 GeV and m ↓ does not exceed 636 GeV. Through the CP mixing all Higgs bosons participate in the PT, also the heavier ones. Additionally, a CP-violating VEVω 3 is generated at T c and feeds into the VEV of also the heavier Higgs bosons. With heavier Higgs bosons participating in the PT, the PT is weakened inducing smaller ξ c values. In order to counterbalance these effects, the Higgs masses overall become lighter and/or they move closer together (cf. the coloured points on the diagonal axis), thus also distributing large portions of the VEV to the lighter among the Higgs bosons.

Implications for LHC phenomenology and benchmark scenarios
In Fig. 5 we investigate how a strong PT affects LHC phenomenology. The left plot shows in the plane (m ↑ − m H ± ) versus (m ↓ − m H ± ) the frequency of the points that pass the constraints, the middle and right plot display the frequency of the points when additionally a strong EWPT is required and only the CP-conserving points (middle) or the CP-violating points (right) are taken into account. The LHC constraints favour degenerate non-SM-like neutral Higgs bosons that are lighter (by at most 70 GeV) or equal to the charged Higgs boson mass, cf. yellow points in Fig. 5 (left). In the CP-conserving case, a strong PT, however, favours a hierarchy between the neutral masses, see Fig. 5 (middle): While the heavier H ↑ is approximately mass degenerate with the charged Higgs, H ↓ is lighter by about 200 GeV. A slight preference is also found for mass degenerate H ↓ and H ± with H ↑ being heavier by 250-350 GeV.
As can be inferred from the comparison of the middle and right plot, most of the points with a strong PT are found in the CP-conserving limit, and we find that the strongest PT is obtained in the alignment limit where H 1 is totally SM-like. In line with Ref. [51] we hence conclude, that a strong PT favours scenarios where decays of the heavier neutral Higgs boson into the lighter one together with a Z boson are kinematically allowed and may have a considerable branching ratio due to the involved coupling being large in the alignment limit. We also find that the mass configuration with a mass gap of 290-330 GeV between H ↓ and H ↑ induces the largest ξ c values (not shown in the plots). When keeping only the CP-violating points (right plot), we find three different types of mass configurations that are compatible with a strong PT. (i) We have the parameter sets where H ↑ and H ± are approximately mass degenerate with H ↓ being lighter by 180 to 320 GeV. (ii) We find two points with mass degenerate H ↓ and H ± and H ↑ begin heavier by 280 GeV and 330 GeV, respectively. (iii) Finally, H ↓ and H ↑ are approximately mass degenerate and H ± is either lighter, heavier or has the same mass. We investigate the features of these points closer by identifying benchmark points for each of these regions. We denote these points by BPi1-3, BPii1-2 and BPiii1-3 for three mass configuration types (i), (ii) and (iii), respectively. In Tables 2, 4 and 6 we list the input parameters of the 3 sets of benchmarks. We also give the derived third neutral Higgs boson mass, the strength of the PT ξ c and the CP admixtures of the Higgs bosons. These are quantified by the mixing matrix element squared relating to the CP-odd neutral component of the Higgs doublets ζ 3 , namely R 2 i3 (i = 1, 2, 3). Finally, we give for all benchmark points the result for the production of a SM-like Higgs pair h through gluon fusion at a c.m. energy of √ s = 14 TeV including the NLO QCD corrections in the heavy top mass limit [121]. We come back to the discussion of Higgs pair production in Section 7. In Tables 3, 5 and 7 we summarise the dominant branching ratios of the benchmark points of the three sets, which determine their phenomenology. For H 1 we have SM-like branching ratios and do not give them separately here.
We start our discussion with the benchmark set BPi1-3. Denoting by φ the heavier mass degenerate Higgs bosons H ↑ = H 3 and H ± , the three benchmarks differ in the their mass difference m φ − m H ↓ =H 2 which is about 200, 250 and 300 GeV for BPi1, BPi2 and BPi3. Since the masses of the heavier Higgs bosons φ are between about 490 and 550 GeV in all three scenarios, this means that m ↓ becomes successively smaller with increasing mass difference. The interplay of the kinematically available phase space and the CP nature of the Higgs bosons determines their branching ratios. In particular, the decays of H ↓ ≡ H 2 turn out to be interesting. In BPi1 and BPi2, where H 2 is mainly CP-even 14     idea has been proposed and discussed before in [122,123]. 15 The reason is that the former decay requires H 2 to be CP-even, whereas the latter requires it to be CP-odd in a purely CP-conserving theory. In BPi3, H 2 has a large CP-odd admixture. Due to its small mass, however, the off-shell decay into ZH 1 is less important than the on-shell decays into massive gauge bosons. Here we make the interesting observation that H 2 , despite its rather CP-odd nature, mainly decays into massive gauge bosons as a consequence of the available phase space. These decays are only possible because H 2 also has a CP-even admixture. The heavier non-SM-like Higgs H ↑ ≡ H 3 in all three scenarios mainly decays into ZH 2 . In BPi1 and BPi2, where H 3 is mainly CP-odd, the next important decay is the one into top-quark pairs. In BPi3 H 3 is more CP-even, so that the second important decay becomes the one into W W . Note, that the decay into ZH 1 is less important than into ZH 2 because of a much smaller involved coupling. In our scenarios, the coupling of H ± to W H 1 is smaller than the one to W H 2 , so that the charged Higgs decays mainly into W H 2 followed by the decay into tb in BPi1 and BPi2 and by W H 1 in BPi3.
The benchmark points BPii1 and BPii2 are rather similar. They differ in the mass gap between the heavier H ↑ ≡ H 3 and the lighter H ↓ ≡ H 2 , which is now almost mass degenerate with H ± . Denoting the latter two by ϕ, we have m H ↑ − m ϕ ≈ 260 GeV in BPii1 and around 310 GeV in BPii2, where H ↓ and H ± are lighter than in BPii1. In both scenarios H 1 ≡ h and H ↓ are mostly CP-even and H ↑ is mostly CP-odd. Therefore, H 2 mainly decays into the massive gauge bosons, cf. Table 5. In BPii1 also the decay into H 1 H 1 has a substantial branching ratio of 0.34. In BPii2, this decay is kinematically closed. 16    are into a gauge+Higgs boson final state, namely into W ± H ∓ and ZH 2 with branching ratios of about 0.7 and 0.3 in both scenarios. In contrast to the parameter set (i), the charged Higgs boson is considerably lighter, so that it mainly decays into tb and the decay into a gauge+Higgs final state is very small. Again, we find non-SM-like decays in these scenarios, where the possible final states reflect the mass hierarchy among the Higgs bosons. Note, finally, that both benchmark points have a ξ c just slightly above the strong PT value of ξ c = 1 and are almost ruled out.
We find four scenarios for which all non-SM-like Higgs bosons are approximately mass degenerate. Denoting by Φ generically H ↓ , H ↑ and H ± , we have for their average mass m Φ ≈ 498 GeV, 547 GeV, 555 GeV and 563 GeV, respectively, for these scenarios. All four scenarios feature the decays into Higgs pairs of the C2HDM is deferred to a future project.  same dominant branching ratios. We exemplary give as benchmark point BPiii1 the scenario with the lightest mass spectrum. The two other benchmark points feature m ↓ ≈ m ↑ , where the charged Higgs boson is heavier (BPiii2) or lighter (BPiii3). In BPiii1, there is no mass gap between the non-SM-like Higgs bosons, and in BPiii2 and BPiii3 the largest mass gap is 80 and 77 GeV, respectively. Furthermore, in all these scenarios the couplings between H 2,3 , Z and h are small. Therefore the non-SM-like Higgs bosons decay dominantly into SM final states, which due to their mass values are tt for the neutral Higgs bosons and tb for the charged Higgs boson. For H 2 which has a significant CP-odd admixture but still is dominantly CP-even, the next important decay channels are those into W W and ZZ.
We can therefore summarise that the requirement of a strong PT induces Higgs spectra with mass gaps that are characterized by large Higgs branching ratios of the non-SM-like Higgs bosons into gauge+Higgs final states or into Higgs pairs which should be testable at the LHC. Some of the decays of our scenarios would be forbidden in a purely CP-conserving 2HDM. In contrast to the CP-conserving case, due to the CP-mixing of all Higgs bosons, also scenarios with the non-SM-like Higgs bosons being close in mass or even mass degenerate are preferred. In these cases the dominant decays are those into SM final states. Figure 6 shows in grey the distribution of the Higgs signal strengths for the scenarios passing the constraints and in colour the ones that are additionally compatible with a strong PT. The fermion initiated cross section (gluon fusion and associated production with a heavy quark pair) of the SM-like Higgs boson h normalised to the SM is denoted by µ F , and the normalised production cross section through massive gauge bosons (gauge boson fusion and associated production with a vector boson) by µ V . The value µ xx is defined as where H SM is the SM Higgs boson with mass 125 GeV. In the right plot we retained only the points with explicit CP violation. Photonic rates of up to 1.15 are still allowed. When imposing a strong PT, however, this reduces to 1.02 for parameter points with CP violation and below 1 for CP-conserving scenarios. Note that we did not find any points with both reduced µ V /µ F and µ γγ by more than 10% that lead to ξ c ≥ 1 in the CP-violating case.
In Fig. 7 we see in the µ τ τ − µ V V plane the distribution of points passing the constraints and how this compares to the points when we require a strong PT (left) and only keep those points that are CP-violating and feature ξ c ≥ 1 (right). Clearly, in the CP-violating case the points with a strong PT are much more sparse and disfavour points with µ τ τ above the SM value of 1 and reduced µ V V . The differences in the rates that we found can be exploited to distinguish the C2HDM from the R2HDM or to exclude the C2HDM. 7 Analysis of the trilinear Higgs self-couplings and Higgs pair production in the C2HDM Type I Having computed the loop-corrected effective potential at finite temperature, we now investigate the effects of the NLO corrections on the trilinear Higgs self-coupling as well as the interplay between a strong PT and the trilinear Higgs self-couplings. The loop-corrected trilinear Higgs self-couplings are obtained from the loop-corrected effective potential by performing the third derivative of the Higgs potential 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 [93].

The Higgs self-coupling between three SM-like Higgs bosons
In Fig. 8  ). Still, a strong PT clearly requires a large enough Higgs self-coupling, larger than the one realized in the SM. Along these lines, we see that the largest ξ c values are obtained for the largest enhancement factors in this range. Including only points with explicit CP violation the maximum enhancement factor is found to be ∼ 3.6 for all points passing the constraints. The scenarios with ξ c > 1 reduce the maximum enhancement factor somewhat with λ C2HDM,NLO  The larger the fraction of the VEV carried by the Higgs boson, the stronger is its participation in the EWPT, so that it should not be too heavy in order not to spoil the strong PT. The largest fraction of the VEV is carried by h, with < h >= (0.87...1) v. The largest ξ c values are obtained in the alignment limit, which is preferred by a strong PT, where the SM-like Higgs boson carries the entire VEV. In this case the remaining two neutral Higgs bosons do not take part in the PT, so that a strong first order PT requires a substantial trilinear Higgs self-coupling for h. A conservative estimate of the prospects of the high-luminosity LHC to measure the trilinear Higgs self-coupling of the SM, concludes that an accuracy of about 50% on its value might be feasible [126]. This allows to distinguish some of the C2HDM scenarios compatible with a strong PT from the SM case.
The impact of the NLO corrections on the trilinear Higgs self-coupling between three SM-like Higgs bosons in the C2HDM is shown by Fig. 9. The left plot shows the NLO coupling as function of the leading order (LO) coupling. The NLO corrections can both suppress and enhance the treelevel coupling. The corrections can be substantial. For the points with a strong PT, the increase can be by up to a factor 8.3, while it is 2.5 for the parameter point with the largest ξ c . As outlined in [125], large corrections, beyond the top loop contribution also present in the SM, arise from Higgs loop contributions in the 2HDM. They increase with the fourth power of the Higgs boson mass, m 4 Φ , where Φ stands generically for the 2HDM Higgs bosons H, A, H ± . And they are suppressed by with M 2 ≡ m 2 12 /(sin β cos β). The masses of the heavy Higgs bosons schematically take the form where f (λ i ) denotes a linear combination of the quartic couplings of the Higgs potential. In case of small M 2 large Higgs masses are generated through large values of λ i . In this case the loop contributions to the Higgs bosons do not decouple and increase with m Φ , see also [127]. In case M 2 f (λ i )v 2 the loop contributions decouple in the limit m 2 Φ ≈ M 2 → ∞ due to the suppression factor Eq. (67). For a strong PT we need large Higgs boson self-couplings, respectively large couplings λ i , inducing the observed large corrections in the non-decoupling regime. On the other hand, the masses of the Higgs bosons participating in the EWPT should not become too heavy, thus restricting the size of the quartic couplings and hence also of the enhancement through the NLO corrections. This explains why in the regime compatible with a strong PT the enhancement of λ hhh remains below the maximum enhancement factor compatible with the applied constraints. The right plot of Fig. 9 displays the ratio between the C2HDM coupling and the SM counterpart at NLO versus this ratio at LO. The NLO ratio deviates substantially from the LO ratio for a large fraction of the parameter points, showing that the NLO effects in the two models can be quite different. This was to be expected as in the C2HDM further Higgs bosons contribute to the loop corrections and their impact can be quite sizeable due to large Higgs self-couplings and/or light masses. Inspecting only pure CP-violating scenarios, we find overall less scenarios compatible with the constraints and a strong PT and the size of the NLO corrections is somewhat reduced.

The benchmark scenario BP3HSM
To quantify the impact of the PT on the Higgs self-coupling and on Higgs pair production, we exemplary give one benchmark point in the CP-violating case, BP3HSM, with the input parameters (M H 3 is derived from the input parameters) with the NLO to LO C2HDM coupling ratio being showing the importance of the loop corrections.
We now discuss the effect of these corrections on continuum Higgs pair production. At the LHC the dominant process is given by gluon fusion [121,126,128,129]. The contributing diagrams are the triangle diagrams with the production of a Higgs or Z boson that subsequently decays into a Higgs pair, and the box diagrams [88]. For the NLO cross section of gluon fusion into the SM-like Higgs pair hh, computed with a private version of HPAIR [130] adapted for the C2HDM [88], we find at a c.m. energy of 14 TeV The QCD corrections computed in the heavy top mass limit yield a K-factor, i.e. a ratio of NLO to LO cross section (the latter calculated with LO strong coupling constant and parton distribution functions), of showing the importance of the QCD corrections. The NLO cross section for SM Higgs pair production computed with full top quark mass dependence amounts to [131][132][133] σ NLO (H SM H SM ) = 32.91 fb .
The C2HDM cross section is hence by a factor of 3.8 larger than in the SM. This cross section does not include any EW corrections, and in particular not the ones given in Eq. (75). The enhancement of 3.8 is mostly due to the resonant production of H ↓ that subsequently decays into an h pair. Without this resonant enhancement the cross section amounts to The quantification of the effect of the EW corrections requires the complete calculation of the Higgs pair production process at NLO EW, which is clearly beyond the scope of this paper. The computation of the loop-corrected effective trilinear Higgs self-couplings gives a flavour of the importance of the EW corrections and the impact of the EWPT on this value In particular, we note that the increase of the trilinear Higgs self-coupling may also decrease the total size of the cross-section due to the destructive interference between triangle and box diagrams. Electroweak baryogenesis which requires a certain size of the trilinear Higgs self-coupling between the SM-like Higgs bosons in order to be successful hence has a direct influence on the size of resonant and continuum Higgs pair production that is significant enough to be tested at the LHC (and future colliders).
We finish this section by commenting on the size of the Higgs pair production cross sections of the benchmark points given in Tables 2, 4 and 6. As can be inferred from the values given in the tables, Higgs pair production in the C2HDM is significantly enhanced compared to the SM value for scenarios where resonant heavy Higgs production with subsequent decay into hh is kinematically possible. In the scenarios we looked at, the H ↑ hh couplings is usually very small, so that the main resonant contribution comes from H ↓ production. In case this is kinematically not allowed or the H ↓ hh coupling is also small, the cross section value compares to the one of the SM.

Further Higgs self-couplings
The inspection of the other trilinear Higgs couplings not involving only the SM-like Higgs boson shows the following: The trilinear C2HDM Higgs self-couplings can be suppressed but also be substantially enhanced compared to the SM trilinear coupling and still be compatible with all constraints. The enhancement factor is less important for scenarios that additionally feature a strong PT. However, it can still be considerable, depending on the self-coupling and the scenario. Also the NLO corrections can be important. Barring the case where the LO coupling is close to zero and hence the coupling becomes effectively loop-induced, the maximum enhancement factor for CP-violating scenarios with ξ c > 1 is 5.7. Due to the large amount of possible trilinear couplings, we exemplary show in Fig. 10 the coupling between the SM-like Higgs boson and H ↓ and H ↑ , λ hH ↓ H ↑ . For this coupling, the enhancement in the C2HDM compared to the SM coupling can be up to a factor 1.5 at NLO in accordance with all constraints, cf. Fig. 10 (left). When additionally a strong PT is demanded, the ratio drops to values between -0.34 and 0.47 the SM coupling. The largest ξ c value is obtained close to the alignment limit where h ≈ v. The right plot shows the impact of the NLO correction. For the sample compatible with all constraints the ratio of the NLO to the LO coupling can be quite large. The demand of a strong PT has a considerable impact, as in this case the ratio becomes much smaller, as can be inferred from the plot.

Type I: Parameter sets with H 2 ≡ h
In this mass configuration our scan resulted in only three scenarios compatible with all constraints that both allow for a strong PT and include CP violation. The results are therefore those of a real 2HDM with the heavier of the two CP-even Higgs bosons being SM-like with a mass of 125 GeV. We reproduced the results of our previous publication on the PT in the CP-conserving 2HDM given in [51], which we briefly summarise here: The scenarios compatible with ξ c > 1 require (neglecting a few outliers) a mass hierarchy where H ↑ , i.e. the pseudoscalar A in the R2HDM, is mass degenerate with H ± and lies in the mass range ∼130 to ∼490 GeV, so that there is a mass gap between m ↑ and m ↓ < 125 GeV. The reason is that due to the required small m ↓ , coinciding with the mass of the lighter of the two CP-even R2HDM Higgs bosons h, the quartic coupling λ 2 and m 2 12 have to be small. Hence we are left with λ 4 and λ 5 that drive the phase transition, implying large H ↑ and H ± masses and the mass gap to H ↓ . Keeping in mind that H 2 , i.e. the heavier of the 2 CP-even Higgs bosons, is SM-like induces sin(β − α) = 0 in the limit of the R2HDM, so that this mass hierarchy allows for A → hZ decays, involving the coupling g AhZ ∼ cos(β − α), and can be probed at the LHC. The upper bound on the masses of the heavy Higgs bosons is given by the fact that the Higgs bosons participating in the PT must not be too heavy.
In the CP-violating 2HDM we barely find any points compatible with ξ c > 1. We have seen that explicit CP violation comes along with spontaneous CP violation at the PT. The thus generated CP-violating VEV at the EWPT feeds into all Higgs bosons as they are all mixing in the CPviolating 2HDM. As the SM-like Higgs boson must not receive a large CP admixture it is either H ↓ or H ↑ that develop a non-negligible CP-violating VEV. Due to the above described mass hierarchy with a heavy H ↑ its VEV should not become too large, however, in order not to weaken the PT. This favours the lighter H ↓ to receive a more important fraction of the VEV or else a hierarchy where all neutral Higgs bosons are rather light and hence close in mass. Already in the R2HDM we see that such hierarchies together with a strong PT are very rare, so that the scenarios that can be found are very sparse.
The phenomenological implications of the R2HDM are the same as found in [51] with the main feature that there are only very few scenarios with a strong PT that yield photonic rates µ γγ beyond 0.9, although values of up to 1.45 would still be compatible with the applied constraints. In contrast, however, to [51] the rate into τ τ can go up to the maximum allowed experimental value of 1.4 also for ξ c > 1, which is due to different, i.e. newer, limits on tan β applied in this work. The three explicitly CP-violating scenarios lie in the same boundaries for the µ-rates as the ones of the R2HDM. As for the trilinear Higgs self-couplings, the overall picture is the same as in the H 1 ≡ h case and we content ourselves to summarise the main features in the conclusions.

Features of the CP-violating scenarios with strong PT
The closer inspection of the 3 CP-violating scenarios reveals that they all feature Higgs spectra with rather close mass values. The largest difference between the heaviest and the lightest Higgs boson mass is 256 GeV. In Table 8 we list the input parameters of the 3 benchmark point scenarios, denoted by BPCPV1-3, featuring a strong PT in the CP-violating case. We additionally give the derived third neutral Higgs boson mass, ξ c , the CP admixtures of the Higgs bosons and the SM-like Higgs pair production cross section through gluon fusion. All three scenarios have ξ c values rather close to 1 underlining the difficulty in finding parameter sets inducing a strong PT in case H 2 is SM-like. In BPCPV1, H 3 has the largest mass of all three benchmark points with 376 GeV. A strong PT is possible as H 3 receives a smaller fraction of the  and 3, all masses are rather close with the mass of H 3 being below 160 GeV, which now also carries a larger fraction of the VEV.
The phenomenological features of the three benchmarks are summarised in Table 9 where we depict the dominant branching ratios of the various Higgs bosons. For H 2 we have SM-like branching ratios and do not give them separately here. The branching ratios are in accordance with the nature of the Higgs bosons. As can be inferred from the R 2 i3 given in Table 8 the lightest Higgs boson of BPCPV1 is mostly CP-even whereas the heaviest one is mostly CP-odd, so that its branching ratios into massive gauge bosons are suppressed. Instead, the large mass gap of more than 250  In type II, all scenarios compatible with a strong PT, where the SM-like Higgs boson is given by the heavier neutral Higgs states are excluded due to the lower bound on the charged Higgs mass, so that necessarily H 1 ≡ h. We start our discussion of the EWPT in the C2HDM with the investigation of the CP-violating phase. Figure 11 shows the relation between tan ϕ spont =ω 3 (T c )/ω 2 (T c ) at the critical temperature T c and tan ϕ explicit = Im(m 2 12 )/Re(m 2 12 ) at vanishing temperature. All points comply with our constraints and feature a strong PT. Furthermore, only points with explicit CP violation are included, i.e. Im(m 2 12 )/Re(m 2 12 ) = 0 although we allow it to be very small. The plot shows the correlation between the two types of CP-violating phases, with the absolute value of the spontaneously generated phase decreasing with increasing absolute value of the explicitly CP-violating phase. Like in the type I, we observe that the spontaneous generation of a CP-violating phase at the critical temperature only appears for scenarios with explicit CP violation at zero temperature. The color bar visualizes the size of ξ c at T c . The maximum value attained in case of CP violation is somewhat smaller than in the type I case, with ξ max c = 1.3. The smaller size of ξ c is to be attributed to the overall heavier mass spectrum in the type II model due to the lower bound on the charged Higgs mass. In the CP-violating case all Higgs bosons mix and hence the heavier Higgs bosons receive contributions from all VEVs, increasing their participation in the PT. Figure 12 displays the total phase at T c as a function of the CP-violating angle at T = 0. We see that the total phase varies between -7.6 • and +7.6 • and is hence much smaller than in type I. This is closely related to the fact that at zero temperature the applied constraints, mainly the EDM constraint, restrict a possible CP-violating phase at T = 0 more strongly than in type I [87], namely to values between about -9 • and 8.6 • . Due to the correlation between ϕ explicit and ϕ spont the total phase at T c is then smaller, so that in the type II C2HDM smaller CP-violating effects are to be expected from this mechanism.  Figure 13 shows the mass values of the neutral non-SM-like Higgs bosons compatible with all constraints (grey) that are additionally compatible with a strong PT (color), taking all points of the CP-conserving and CP-violating scenarios (left) and restricting to purely CP-violating scenarios (right). The results of Fig. 13 (left) basically agree with the results found in [51], taking into account the fact that the lower bound on m H ± has moved up to 580 GeV. Furthermore, we do not find valid points with m ↓ < 250 GeV, which would come along with large mass gap of the heavier Higgs boson masses to m ↓ . Overall, grey points with mass gaps between H ↓ and H ↑ above 332 GeV are not allowed any more. This exclusion results from the unitarity check with the NLO Higgs self-couplings. The plot confirms that for the CP-conserving parameter points with ξ max c = 6.5 larger ξ c values can be obtained than in the CP-violating case where ξ max c = 1.3. The right plot shows that the inclusion of CP violation implies mass spectra where overall the non-SM-like Higgs masses move closer, cf. also Fig. 14. In particular most of the points with a strong PT feature H ↓ and H ↑ which are close in mass, and the largest ξ c values are found for the lightest possible values that they can have within the given constraints. Again this can easily be understood by reminding that due to CP violation all Higgs bosons mix and have a non-negligible VEV and that additionally in type II with H 1 ≡ h the non-SM-like neutral Higgs bosons are rather heavy. In order not to weaken the PT too much, in case of H ↑ having a large portion of the VEV, it should be as light as possible and hence be mass degenerate with H ↓ , or in case it is not mass degenerate, the lighter of the two should acquire most of the VEV.

Implications for LHC phenomenology
The implications for LHC phenomenology can be read off Fig. 14, showing the mass differences of H ↑ and H ± versus the ones of H ↓ and H ± . The colour code shows the relative frequency for all points passing the constraints (left), with additionally a strong PT and only CP-conserving points (middle) and for only CP-violating points with ξ c ≥ 1 (right). While the application of the constraints favours scenarios with degenerate neutral non-SM-like Higgs masses, the requirement of a strong PT favours a mass hierarchy between m ↑ ≈ m H ± and m ↓ with a mass gap of about 200 GeV. The comparison with the right plot shows that these are mostly CP-conserving scenarios. This allows for decays of H ↑ → ZH ↓ that can be searched for at the LHC. In contrast the right plot shows that the CP-violating case favors more degenerate masses of the non-SM-like Higgs bosons. This has consequences for LHC phenomenology. In order to quantify this we have chosen four exemplary benchmark scenarios featuring different mass patterns. We have one benchmark point, BP1T2, with almost degenerate H ↑ and H ± and a mass gap to the lighter H ↓ . In BP2T2, H ↓ and H ± are closer in mass and H ↑ is the heaviest Higgs boson. The benchmarks BP3T2 and BP4T2 feature nearly mass degenerate H ↓ and H ↑ with a heavier charged Higgs boson in the former and a lighter one in the latter case. The input parameters are given in Table 10. The overall mass spectrum is heavier than in the type I C2HDM as expected from the lower bound on the charged Higgs mass in type II. The dominant branching are summarised in Table 11. These are determined by the mass pattern together with the fact that the H 3 V V (V = W, Z) coupling is very small despite R 2 33 being small in BP1T2-BP3T2. Besides the dominant decay into tt, H 3 has a substantial branching ratio into ZH 2 in BP1T2 and in into W ± H ∓ in BP2T2. In BP3T2 and BP4T2, the mass pattern forces H 3 to mainly decay into tt. This is the dominant decay channel for H 2 in all four scenarios. The mass ordering of BP1T2 allows H ± to decay with a significant branching ratio into W ± H 2 besides the dominant decay into tb. For all the other scenarios H ± almost exclusively decays into tb. Again we find that the fact that in the C2HDM not only mass hierarchies with large mass gaps are preferred by the strong PT, induces also decay patterns with SM particles in the final state.
In Fig. 15 we depict µ V /µ F versus the photonic rate µ γγ in grey for all points passing the constraints and in colour for those with ξ c > 1. The left plot comprises all 2HDM points, while in the right plot only the CP-violating points are retained. The requirement of a strong PT restricts the region of allowed µ values, and the results of the left plot reconfirm our findings of [51] in the real 2HDM. Due to a more efficient scan procedure applied here, we have now more scenarios compatible with ξ c ≥ 1 in the wrong-sign limit region [98,134,135]    in the lower left corner) where the h coupling to the massive gauge bosons has an opposite sign with respect to its b-quark Yukawa coupling. The inclusion of CP violation restricts the µ-values further. We see a strong correlation of µ γγ and µ V /µ F , with the latter decreasing with increasing µ γγ . The wrong-sign limit is completely excluded which, however, is already almost excluded due to the applied constraints and not because of a strong PT. The CP-violating scenarios with ξ c ≥ 1 preclude µ V /µ F above 1.17 and below 0.86 and µ γγ above 1.12 and below 0.76. Any values outside these ranges point to the CP-conserving 2HDM. Note also, that the SM point (red triangle) is not compatible with a strong PT in the C2HDM.
As can be inferred from Fig. 16, which is the same as Fig. 15 but for µ τ τ versus µ V V , in the CP-violating case there is also a strong correlation between µ τ τ and µ V V . The former increases with µ V V . Furthermore, no scenarios in the wrong-sign regime (corresponding to the points in the upper left corner of Fig. 16 (left)) are realized, which, again, is mostly due to the applied constraints. The value of µ τ τ is restricted to values between 0.9 and 1.10, as well as µ V V to the range between 0.8 and 1.22. 10 Analysis of the trilinear Higgs self-couplings and Higgs pair production in the C2HDM Type II We conclude our investigation with the discussion of the interplay between the requirement of a strong PT and the trilinear Higgs self-coupling among the SM-like Higgs bosons. Figure 17 (  with respect to the SM, with a maximum enhancement factor of 5.1. If additionally the PT is required to be of strong first order, the trilinear coupling has to be somewhat larger than in the SM, but must not be above a factor of about 2.9 the SM value, in order not to weaken the PT again due to too large Higgs masses. If only CP-violating points are taken into account this ratio is slightly lowered to 2.4. The right plot shows the NLO coupling versus the LO one and underlines the importance of the loop corrections. The wide spread points in the upper half of the plot and the two isolated points in the lower right part are due to the wrong-sign limit. In the regions with a strong PT the maximum corrections amount to about a factor 3.3, and the largest values of ξ c are found for the largest correction factor. Taking into account only the CP-violating scenarios, the maximum enhancement is lowered to 2.2. Note, that the sharp cut of the tree-level couplings at large values are due to tree-level unitarity. Interestingly, the plot shows that a vanishing trilinear Higgs self-coupling between the SM-like Higgs bosons is not compatible with the constraints any more, while in type I this was not excluded. Concerning the trilinear Higgs self-couplings also involving non-SM-like Higgs bosons we find the same overall features as in the type I case.
In Table 10 we give the NLO QCD hh production cross sections through gluon fusion at √ s = 14 TeV for the four benchmark points BP1-4T2. The size of the non-SM-like Higgs boson masses allows for resonant enhancements through their on-shell production with subsequent decay into hh. Depending on the size of the Higgs self-couplings between them and hh, the C2HDM hh cross section is more or less enhanced compared to the SM.

Conclusions
We have found that the NLO effects derived from the effective potential have a non-negligible influence on the global minimum and perturbativity at NLO. The requirement of the EW minimum to be the global minimum also at NLO excludes O(5 − 25%) of the generated scenarios compatible with the applied constraints in the type I C2HDM with H 1 or H 2 ≡ h, denoted by C2HDM(I H 1 ), and C2HDM(I H 2 ) and the type II C2HDM with H 1 ≡ h, denoted by C2HDM(II H 1 ). In the type II C2HDM with H 2 ≡ h, the mass differences enforced by the experimental constraints become large and, due to our renormalisation procedure, induce large corrections to the Higgs self-couplings so that a stable vacuum cannot be guaranteed any more so that we did not investigate this version of the C2HDM further in this paper. Along the same lines accordance with unitarity at NLO eliminates O(9 − 18%) of the remaining points.
We showed that the presence of explicit CP violation at zero temperature induces spontaneous CP violation at the EWPT. The size of the induced phase correlates with the amount of explicit CP violation at T = 0. In type I a larger CP-violating phase at zero temperature is compatible with the experimental constraints than in type II. Consequently, the total CP-violating angle, including the spontaneously generated CP-violating phase, at T c amounts up to 49 • in the C2HDM(I H 1 ) and up to about 8 • in the C2HDM(II H 1 ). In particular, in type I this should be large enough to ensure successful baryogenesis.
We overall re-confirm our findings of our previous investigation of the PT in the CP-conserving 2HDM for the scenarios in the CP-conserving limit. Deviations in type I occur due to an updated limit on tan β and in type II due to a more efficient scan procedure. A strong phase transition is found to be possible in all three different remaining set-ups including CP violation, i.e. the C2HDM(I H 1,2 ) and the C2HDM(II H 1 ). Overall, the strength of the phase transition is smaller in the CP-violating case, with ξ max c = 1.9 in C2HDM(I H 1 ) and ξ max c = 1.3 in C2HDM(II H 1 ). In C2HDM(I H 2 ) we find only three points (BPCPV1-3) that are compatible with a strong PT, with ξ max c = 1.48. Although the SM-like Higgs boson has the largest fraction of the VEV at the PT and hence mainly drives its strength, through CP mixing all three neutral Higgs bosons, also the non-SM-ones, receive a VEV. Scenarios where heavy non-SM-like Higgs bosons receive a large portion of the VEV at the EWPT weaken its strength. Therefore either the lighter of the non-SM-like Higgs bosons receives the larger portion of the VEV among the two, or the overall spectrum is as light as possible.
Concerning the mass hierarchies and hence the implications for phenomenology, we confirm for scenarios compatible with a strong PT in the CP-conserving limit the preference for a mass hierarchy with a gap between the heavier and the lighter neutral non-SM-like Higgs bosons. The heavier Higgs boson therefore mainly decays into gauge+Higgs final states. Depending on the value of the charged Higgs mass, which due to electroweak precision physics constraints has to be mass degenerate with one of the neutral Higgs bosons, these final state particles are neutral (ZH i ) or charged (W ± H ∓ ). Such decays are a clear indication of non-SM physics arising from an enlarged Higgs sector, which can be searched for at the LHC. For CP-violating scenarios, however, no such preference is found: The requirement of a strong phase transition combined with a CP-violating setup where all Higgs bosons receive a VEV, besides this mass pattern, also favours scenarios where all three neutral Higgs bosons are rather close in mass. There are scenarios with one or two of the non-SM-like neutral Higgs bosons degenerate with the charged Higgs boson. Without mass hierarchies between the non-SM-like Higgs bosons (and small gauge-Higgs pair couplings involving the h) they all decay into SM final states, with the specific nature of the final state particles determined by the mass of the decaying Higgs boson. Among the scenarios with a mass hierarchy we also have Higgs-to-Higgs decays and decay patterns that allow to identify the CP-violating nature of the Higgs boson through the combination of its decay channels. We have provided several benchmark mass pattern excluded µ V /µ F − µ γγ µ τ τ − µ V V C2HDM(I H 1 ) µ V /µ F < 0.9 and µ γγ < 0.9 µ τ τ > 1 and µ V V < 1 C2HDM(II H 1 ) µ γγ > 1.12 µ V V > 1.22 µ γγ < 1 and µ V /µ F < 0.82 (wsl) µ V V < 1.13 and µ τ τ > 1.07 (wsl) µ V /µ F = µ γγ = 1 (SM point) Table 13: Distinguishing features between the C2HDM and 2HDM rates for scenarios with ξ c > 1: The C2HDM excluded regions given here are still allowed by the 2HDM. Note that the wrong-sign limit (wsl) regions are excluded due to the applied constraints and not due to ξ c > 1.
in Table 13. The regions are excluded by the C2HDM but not by the 2HDM. Any measured value outside these regions rules out the CP-violating set-up for successful baryogenesis. Note, however, that the wrong-sign-limit is mainly disfavoured by the applied constraints and not by the requirement of ξ c > 1 in contrast to the other regions given in the Table. From the loop-corrected effective potential we derive the NLO-corrected effective trilinear Higgs couplings. The C2HDM couplings can be enhanced or suppressed relative to the SM case for all scenarios compatible with the experimental and theoretical constraints. (Note, that the maximum allowed value is constrained by the requirement of unitarity.) The scenarios with a strong PT, however, require enhanced trilinear Higgs self-couplings between the SM-like Higgs bosons, which is the Higgs boson with the largest VEV and hence mainly drives the phase transition. Including also the points in the CP-conserving limit, the two enhancement regions are λ C2HDM,NLO Higgs bosons shows that the C2HDM values can be enhanced or suppressed compared to the SM, that the enhancement is less important for scenarios with a strong PT and that the NLO corrections are significant. We summarise that the requirement of a strong PT induces enhanced trilinear Higgs self-couplings among the SM-like Higgs bosons that deviate from the SM value significantly enough to be measurable at the LHC. The experimental confirmation of large self-couplings establishes a clear connection between a strong PT, i.e. cosmology, and collider physics. The self-coupling is extracted from Higgs pair production at the LHC, with the main channel given by gluon fusion. We calculated for our benchmark scenarios the cross sections for hh production including the NLO QCD corrections in the large top quark mass limit. Due to the possible resonant production of heavy Higgs bosons with subsequent decay into hh they are all enhanced compared to the SM case, so that the prospects for measuring them are promising. These cross sections do not include EW corrections, as they are not available. The derived sizes of the EW-corrected effective Higgs self-couplings indicate that these corrections can be sizeable. They could lower the Higgs pair production cross section.
In conclusion, the requirement of a strong PT with successful baryogenesis demands an enlarged Higgs sector and has measurable consequences for the C2HDM. In contrast to the real 2HDM not only mass gaps between the Higgs bosons, but also degenerate scenarios are favoured. The C2HDM with a strong PT predicts strong correlations among the signal strengths of the SM-like Higgs boson. Finally, the trilinear hhh coupling must be enhanced compared to the SM, and the additional Higgs bosons induce hh production cross sections that are larger than in the SM case and can be measured at the LHC. The combination of successful baryogenesis with collider phenomenology is a powerful tool to further restrict the underlying model and to identify its true nature.