Strong First Order Electroweak Phase Transition in the CP-Conserving 2HDM Revisited

The discovery of the Higgs boson by the LHC experiments ATLAS and CMS has marked a milestone for particle physics. Yet, there are still many open questions that cannot be answered within the Standard Model (SM). For example, the generation of the observed matter-antimatter asymmetry in the universe through baryogenesis can only be explained qualitatively in the SM. A simple extension of the SM compatible with the current theoretical and experimental constraints is given by the 2-Higgs-Doublet Model (2HDM) where a second Higgs doublet is added to the Higgs sector. We investigate the possibility of a strong first order electroweak phase transition in the CP-conserving 2HDM type I and type II where either of the CP-even Higgs bosons is identified with the SM-like Higgs boson. The renormalisation that we apply on the loop-corrected Higgs potential allows us to efficiently scan the 2HDM parameter space and simultaneously take into account all relevant theoretical and up-to-date experimental constraints. The 2HDM parameter regions found to be compatible with the applied constraints and a strong electroweak phase transition are analysed systematically. Our results show that there is a strong interplay between the requirement of a strong phase transition and collider phenomenology with testable implications for searches at the LHC.


Introduction
In 2012 the LHC experiments ATLAS and CMS announced the discovery of the long-sought Higgs boson [1,2]. Although it looks very SM-like [3][4][5][6] it is quite possible that it is the scalar particle of a Higgs sector beyond the SM (BSM). Despite the success of the SM, which has been tested to highest precision at previous and current colliders, there are still a lot of open questions that cannot be answered within the SM and call for new physics (NP) extensions. One of the unanswered problems is the origin of the observed matter-antimatter asymmetry of the universe [7]. Electroweak (EW) baryogenesis is an elegant mechanism to explain this asymmetry [8][9][10][11][12][13][14][15][16], which is related to physics at the weak scale, establishing a link between collider phenomenology and cosmology. The asymmetry can be generated provided the EW phase transition (PT) taking place in the early universe is of strong first order [14,16] and that all three Sakharov conditions [17] are fulfilled, namely baryon number violation, C and CP violation and departure from the thermal equilibrium. The strong first order PT, proceeding through bubble formation, suppresses the baryon number violating sphaleron transitions in the false vacuum [18,19]. CP-violating reflections of top quarks from the bubble wall produce a hypercharge asymmetry which is converted into a baryon asymmetry in the false vacuum. This asymmetry is transferred to the true vacuum when it passes the bubble wall [20], provided there is departure from the thermal equilibrium. Although in the SM all three Sakharov conditions are fulfilled, the electroweak PT is not of first order [21]. Not only the Higgs boson mass is too large [22], but in addition the CP violation of the SM from the Cabibbo-Kobayashi-Maskawa matrix is too small [16,20,23]. This calls for physics BSM. Among the plethora of NP extensions the 2HDM [24,25] belongs to the simplest models that are in accordance with present experimental constraints. Its Higgs sector features five physical Higgs bosons, three neutral and two charged ones. Their contributions to the effective Higgs potential can strengthen the PT and in addition introduce new sources of CP violation. Previous studies have shown that 2HDMs provide a good framework for successful baryogenesis [26][27][28] (see [29][30][31][32][33] for studies in the CP-violating 2HDM).
In this work we will investigate the implications of a strong first order PT required by baryogenesis on the LHC Higgs phenomenology in the framework of the CP-conserving 2HDM. For this purpose we compute the one-loop corrected effective potential at finite temperature [34][35][36] including daisy resummations for the bosonic masses [37] in two different approximations for the treatment of the thermal masses [38,39]. The renormalisation of the loop-corrected potential is chosen such that not only the vacuum expectation value (VEV) and all physical Higgs boson masses, but, for the first time, also all mixing matrix elements remain at their tree-level values. This allows to efficiently scan the whole 2HDM parameter space with the tree-level masses and mixing angles as input and at the same time test the compatibility of the model with the theoretical and experimental constraints. The points passing these tests will be investigated with respect to a first order PT. The loop-corrected Higgs potential will be minimised at increasing non-zero temperature to find the vacuum expectation value v c at the critical temperature T c , defined as the temperature where two degenerate global minima exist. A value of v c /T c 1 larger than one is indicative for a strong first order PT [11,43]. In our analysis we will discard points leading to a 2-stage PT [44,45]. We will perform a systematic and comprehensive investigation of the 2HDM in four configurations, given by the 2HDM type I and type II where either the lighter or the heavier of the two CP-even Higgs bosons is identified with the SM-like Higgs boson. We will test the compatibility of the model with both the experimental constraints and a strong EW phase transition. The thus delineated regions in the parameter space will be further investigated with respect to their implications for collider phenomenology. We find that the link between cosmology and high-energy collider constraints provides a powerful tool to further constrain the allowed parameter regions of the 2HDM. At the same time, the demand for a strong first order PT leads to testable consequences at the collider experiments.
The outline of the paper is as follows: In section 2 we introduce our notation and provide the loop-corrected effective potential at non-vanishing temperature. In the subsequent section 3 we describe in detail the renormalisation procedure, which is chosen such that at zero temperature the tree-level position of the minimum and the masses and mixing matrix elements of the scalar particles are preserved by the one-loop potential. Using the Higgs boson masses and mixing angles as input parameters, this simplifies the verification of the compatibility of the model with the Higgs data. In section 4 the basic elements of our numerical analysis are described, namely the minimisation procedure of the effective potential in 4.1 and, in 4.2, the details of the scan in the 2HDM parameter space together with the applied theoretical and experimental constraints. Section 5 is devoted to our results. We present the parameter regions compatible with the applied constraints and a strong first order PT, and we then analyse the implications for collider phenomenology. We end in section 6 with our conclusions. The paper is accompanied by an appendix containing the formulae for the masses of the relevant particles and, where appropriate, for the daisy resummed mass corrections.

The Effective Potential
In this section we provide the loop-corrected effective potential of the CP-conserving 2HDM for non-vanishing temperature. First, we set our notation by introducing the model under investigation.

The CP-conserving 2-Higgs-Doublet Model
In terms of the two SU (2) L Higgs doublets Φ 1 and Φ 2 , the tree-level potential of the 2HDM with a softly broken Z 2 symmetry, under which the doublets transform as The mass parameters m 2 11 and m 2 22 and the couplings λ 1..4 are real parameters of the model. The mass and coupling parameters m 2 12 and λ 5 can in general be complex, thereby offering new sources of explicit CP violation in the Higgs sector. We take them to be real as we work in the CPconserving 2HDM. After EW symmetry breaking the two Higgs doublets acquire VEVsω i ∈ R (i = 1, 2, 3), 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, where, without loss of generality, the complex part of the VEVs has been rotated to the second doublet exclusively. Denoting the VEVs of our present vacuum at zero temperature by 2 whereas the remaining two VEVs are related to the SM VEV by Introducing the angle β through we have Through the complex part of the VEV,ω 3 , we include the possibility of generating at one-loop and/or non-zero temperature a global minimum that is CP-violating. 3 The angle β coincides with the angle of the rotation matrix from the gauge to the mass eigenstates of the charged Higgs sector, and also of the neutral CP-odd sector. The physical states of the charged sector are given by the charged Higgs bosons H ± with mass m H ± and the charged Goldstone bosons G ± which are massless at zero temperature, 2 Strictly speaking, T = 2.7 K. Setting T = 0 does not make a difference numerically. 3 In the 2HDM we can have three different types of minima: the normal EW breaking one, a CP-breaking minimum, and a charge-breaking (CB) vacuum. It has been shown that, at tree level, minima which break different symmetries cannot coexist [46][47][48]. This means that, if a normal minimum exists, all CP or CB stationary points are proven to be saddle points. Recent studies of the Inert 2HDM at one-loop level [49], which apply the effective potential approach, indicate that these statements may not be true any more once higher order corrections are included. We therefore allow for the possibility of a CP-breaking vacuum. Including the possibility of a charge breaking Higgs VEV makes the present analysis considerably more complex.
For the neutral CP-odd fields ψ 1 and ψ 2 the same rotation yields the physical states A with mass m A and the neutral Goldstone boson G 0 , massless at zero temperature, Finally, in the neutral CP-even sector the rotation with the angle α transforms the fields ζ 1 and ζ 2 into the two physical CP-even Higgs bosons H and h with masses m H and m h , respectively, In the minimum of the potential Eq. (2) the following minimum conditions have to be fulfilled, with the brackets denoting the Higgs field values in the minimum, i.e. Φ i = (0, v i / √ 2) at T = 0. This results in two equations Exploiting the minimum conditions of the potential at zero temperature, we use the following set of independent parameters of the model, Due to the imposed Z 2 symmetry each of the up-type quarks, down-type quarks and charged leptons can only couple to one of the Higgs doublets so that flavour-changing neutral currents at tree level are avoided. The possible combinations of Yukawa couplings of the Higgs bosons to up-type quarks, down-type quarks or charged leptons are classified as type I, type II, lepton-specific and flipped and are defined in Table 1. The resulting couplings of the fermions normalised to the SM couplings can be found in [50]. In this work we focus on real 2HDMs of type I and type II.

One-Loop Effective Potential at Finite Temperature
The one-loop contributions V 1 to the effective potential consist of two parts: the Coleman-Weinberg (CW) contribution V CW [34] which is already present at zero temperature, and the contribution V T accounting for the thermal corrections at finite temperature T . The one-loop corrected effective potential then reads The tree-level potential is given in Eq. (2) with the doublet Φ 1 replaced by the classical constant field configuration Φ c 1 = (0, ω 1 / √ 2) and Φ 2 by Φ c 2 = (0, (ω 2 + iω 3 )/ √ 2). The Coleman-Weinberg potential in the MS scheme is given by [36] where the sum extends over the Higgs and Goldstone bosons, the massive gauge bosons, the longitudinal photon and the fermions, The m 2 i is 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, 3). The explicit formulae can be found in App. A. The sum also includes the Goldstone bosons. Although we work in the Landau gauge, where they are massless at T = 0, they can acquire a mass if the mass eigenvalues are determined at field configurations other than the tree-level VEVs at T = 0, which is required in the minimisation procedure. Moreover, due to temperature corrections specified below, the masses of the Goldstones and the longitudinal photon can be non-zero, which enforces also the inclusion of γ L in the sum. Note, that due to the choice of the Landau gauge there are no ghost contributions. The variable s i denotes the spin of the particle, n i represents the number of degrees of freedom. Also for later use, we define the degrees of freedom of all particles involved in the model. These are the neutral scalars Φ 0 ≡ h, H, A, G 0 , the charged scalars Φ ± ≡ H ± , G ± , the leptons l, the quarks q and the longitudinal and transversal gauge bosons, V L ≡ Z L , W L , γ L and V T ≡ Z T , W T , γ T , with the respective n i , In the MS scheme employed here, the constants c i read We fix the renormalisation scale µ by µ = v = 246.22 GeV.
In the thermal corrections V T we include the daisy resummation [37] of the n = 0 Matsubara modes of the longitudinal components of the gauge bosons W L , Z L , γ L and the bosons Φ 0 , Φ ± , which adds to their masses at non-zero temperature the Debye corrections given in App. A. The thermal contributions V T to the potential can be written as [35,36] The sum extends over k = W L , Z L , γ L , W T , Z T , Φ 0 , Φ ± , f . Note, that the Goldstone bosons and the longitudinal part of the photon, which are massless at T = 0, acquire a mass at finite temperature and are included in the sum. Denoting the mass eigenvalue including the thermal corrections for the particle k by m k , J (k) ± is given by (see e.g. [51]) with the thermal integrals where J + (J − ) applies for k being a fermion (boson). For each temperature T we determine the VEVsω i , i.e. the field configurations {ω} ≡ {ω 1 ,ω 2 ,ω 3 }, that minimise the loop-corrected potential V , Eq. (17). These enter the tree-level mass matrices such that the masses m i depend implicitly on the temperature T throughω i =ω i (T ). The m k furthermore depend explicitly on T through the thermal corrections. The definition of J (k) ± Eq. (22) is the approach chosen in [38]. A different prescription for implementing the thermal corrections is proposed by [39] where the Debye corrections are included for all the bosonic thermal loop contributions 5 , so that we have In this case, the Debye corrected masses are also used in the CW potential Eq. (18) [31]. We refer to the first approach, Eq.(22), as 'Arnold-Espinosa' and to the second one, i.e. Eq.(24) together with V CW including the thermal corrections in the bosonic masses, as 'Parwani' method. The two approaches differ in the organisation of the perturbative series and hence by higher order terms. The 'Arnold-Espinosa' method consistently implements the thermal masses at one-loop level in the high-temperature expansion, leading to Eq. (22). The 'Parwani' method admixes higher-order contributions, which at one-loop level could lead to dangerous artefacts. Therefore, in the discussion of our results we will apply the 'Arnold-Espinosa' method. The 'Parwani' method will be used only to make contact to previous results in the literature.
Since in the minimisation procedure the numerical evaluation of the integral Eq. (23) at each configuration in {ω} and T is very time consuming, the integrals J ± are approximated by a series in x 2 ≡ m 2 /T 2 . For small x 2 we use [29] J +,s (x 2 , n) = − 7π 4 360 + π 2 24 x 2 + 1 32 with c + = 3/2 + 2 log π − 2γ E and c − = c + + 2 log 4 , where γ E denotes the Euler-Mascheroni constant, ζ(x) the Riemann ζ-function and (x)!! the double factorial. For large x 2 the expansion for both fermions and bosons reads [29] J ±,l ( with Γ(x) denoting the Euler Gamma function. In order to interpolate between the two approximations, first the point is determined where the derivatives of the low-and high-temperature expansions can be connected continuously. We then add a small finite shift to the small x 2 expansion such that also the two expansions themselves are connected continuously. We denote 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 We find that for small x 2 the expansion J +,s for fermions approximates the exact result well by including terms of up to order n = 4, while 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 . This way, the deviation of the approximate results from the numerical evaluation of the integrals is less than two percent. The above approximations Eqs. (25)- (28) are only valid for m 2 ≥ 0. For bosons this is not necessarily the case as the eigenvalues of the mass matrix of the neutral Higgs bosons can become negative depending on the configuration {ω} and the temperature in the minimisation procedure. If this happens, the value of the integral J − , given by Eq. (23), is set to the real part of its numerical evaluation which is the relevant contribution when extracting the global minimum [52] 6 . In practice, we evaluated the integral numerically at several equidistant points in m 2 /T 2 < 0, and in the minimisation procedure we use the result obtained from the linear interpolation between these points, which leads to 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 Coleman-Weinberg potential, Eq. (18), in the one-loop effective potential Eq. (17) contributes already at T = 0, so that the masses and mixing angles obtained from the one-loop effective potential differ from those extracted from the tree-level potential Eq. (2). The loop-corrected masses obtained in this way correspond to the full one-loop corrected masses in the approximation of vanishing external momenta. When we test for compatibility of the model with the experimental constraints the loop-corrected masses and the loop-corrected mixing angles, which enter the couplings, have to be taken into account. For an efficient scan over the parameter space of the model in terms of the input parameters Eq. (16), however, it is more convenient to have the one-loop masses and angles directly as inputs, i.e. they should be the same as the tree-level ones. This can be achieved by an appropriate renormalisation prescription, which will be described in the following.
The Coleman-Weinberg potential given in Eq. (18) has already been renormalised in the MS scheme. We modify this scheme by including finite terms in the counterterm potential that ensure the one-loop corrected masses and, for the first time, also the mixing matrix elements to be equal to the tree-level ones. 7 Introducing counterterms for each of the parameters of the tree-level potential Eq. (2), the counterterm potential V CT added to the one-loop effective potential Eq. (17), reads The complete potential of Eq. (30) will be minimised to find the global minimum at a given temperature T . As stated above, the counterterms δp for the parameters p of the tree-level potential contain only the finite pieces, as the divergent ones have already been absorbed by the MS renormalised V CW . We renormalise the effective potential such that at T = 0 the tree-level position of the minimum yields a local minimum, which is checked to be the global one numerically. Furthermore, through our renormalisation the masses and mixing angles of the scalar particles are preserved at their tree-level values by the one-loop potential. The corresponding renormalisation conditions are imposed at T = 0, which is where we test for the compatibility with the experimental constraints. The position of the minimum is determined by the first derivative of the potential, whereas the masses and angles result from the second derivative, namely the mass matrix. Formulae for both the first and the second derivatives of the CW potential in the Landau gauge have been derived in [53]. We employ these formulae in the gauge basis to calculate the required derivatives. Consequently, for the renormalisation we also express the counterterm potential and the tree-level potential in the 7 Previous works included only the VEVs and (subsets of) the masses in the renormalisation conditions and required them to be equal to their tree-level values [27][28][29][30][31]. In models with extended Higgs sectors, the mixing angles, which enter all Higgs boson observables through the Higgs couplings, are crucial for the interpretation of the results. They are determined from the diagonalisation of the loop-corrected mass matrix. The renormalisation of the mixing matrix elements to their tree-level values guarantees that the relevant quantities and observables constraining the model can be tested with the tree-level input parameters.
gauge basis. The renormalisation conditions for the first derivatives are then given by (i = 1, ..., 8) and φ c T =0 denoting the field configuration in the minimum at T = 0, This results in two non-trivial conditions for the tree-level minimum at T = 0 to be a CP-conserving extremum also at the one-loop level. In order to ensure that both the masses and the mixing angles remain at their tree-level values the complete 8 × 8 mass matrix of the scalar sector should be preserved at its tree-level value by the renormalised one-loop potential. This is achieved by However, since we have only eight counterterms and after imposing Eq. (32) we are left with six to be set, the resulting system of equations is overconstrained and cannot in general be solved. This means that we cannot renormalise all masses and mixing angles to exactly match their tree-level values. We therefore pursue the following approach: Both the tree-level and the one-loop mass matrix are rotated to the mass basis with the tree-level rotation matrix. From the resulting 8 × 8 matrix we extract only the 2 × 2 submatrix, that corresponds to the physical charged Higgs bosons, and the 3 × 3 submatrix for the neutral Higgs bosons. In the CP-conserving case treated here, the latter decomposes into a 2 × 2 matrix for the CP-even Higgs bosons h and H, and the entry for the pseudoscalar A. On these submatrices the renormalisation conditions are imposed, so that we have and The subscript 'mass' indicates that the mass matrix in the gauge basis is rotated into the mass eigenbasis by means of the rotation matrix that diagonalises the tree-level mass matrix. The superscripts H ± and h, H, A indicate that from the resulting matrix only the 2 × 2 block for the physical charged Higgs bosons and the 3 × 3 block for the physical neutral Higgs bosons is considered, respectively. Equations (36) and (37) provide five independent non-trivial renormalisation conditions. 8 Together with the two renormalisation conditions from Eq. (32) we have altogether seven renormalisation conditions to fix eight renormalisation constants, cf. Eq. (31), so that one renormalisation constant is left for determination. Inspecting the counterterm potential Eq. (31), we observe that the counterterms δλ 3 and δλ 4 only appear as sum. Hence, we choose to use only one of them and set δλ 4 = 0. The remaining seven renormalisation constants are fixed by the conditions Eqs. (32), (36) and (37).
We find that these renormalisation conditions allow us to preserve the minimum, the masses and the mixing angles of the Higgs sector at their tree-level values up to a very good approximation. Taking into account numerical uncertainties, the minimum at one-loop remains at v ± 2 GeV, and all masses and mixing angles are preserved up to tiny numerical fluctuations.
Equations (36) and (37) require the second derivative of the CW potential. It is a well-known problem that this derivative leads to infrared divergences for the Goldstone bosons in the Landau gauge [27,29,31,[53][54][55]. In order to circumvent this problem, in [29] the logarithm is redefined to capture on-shell effects regularising the divergence while in [27,28,31] a non-vanishing infrared mass for the Goldstones is employed to regulate the divergence. In the effective potential approach itself, which is the approximation of the full theory at vanishing external momenta, it is not possible to cancel these divergences. Building up the complete self-energy of the Higgs bosons from the effective potential and the momentum-dependent parts obtained by a diagrammatic calculation, however, it becomes apparent that the infrared divergences from the Goldstone contributions cancel between the CW part and the momentum-dependent part [53,55,56]. Taking the limit of vanishing external momenta afterwards we arrive at a finite expression for the second derivative of the CW potential. This cancellation was checked explicitly using the results from the diagrammatic calculation performed in [57,58]. In practice, this result can be obtained directly from the effective potential approach by regularising the logarithmic divergence with a regulator mass and then discarding the terms proportional to this logarithm [53]. The obtained results are independent of the regulator mass and reflect the correct contributions present in the effective potential approach. 9 4 Numerical Analysis

Minimisation of the Effective Potential
The electroweak PT is considered to be strong if the ratio between the VEV v c at the critical temperature T c and the critical temperature T c is larger than one [11,43], The value v at a given temperature T is obtained as Remind thatω 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. For the determination of T c the effective potential together with the counterterm potential, Eq. (30), is minimised numerically at a given temperature T . In a first order electroweak PT the VEV jumps from v = v c at the temperature T c to v = 0 for T > T c . In order to double-check the results of the minimisation procedure, we apply two different minimisation algorithms. One is the active CMA-ES algorithm as implemented in libcmaes [59].
This algorithm finds the global minimum of a given function. As termination criterion we choose the relative tolerance of the value of the effective potential between two iterations to be smaller than 10 −5 . The other algorithm that has been used is the local Nelder-Mead-Simplex algorithm from the GNU Scientific Library [60] (gsl multimin fminimizer nmsimplex2), also with a tolerance of 10 −5 . For a given temperature, we start with 500 randomly distributed points in the interval ω 1,2,3 ∈ [−500, 500] GeV for which we compute the minimum of the potential. Note that we have includedω 3 in Eq. (4) for the sake of generality. The candidates for the global minimum obtained with the two algorithms are compared to each other and the one with the lower value of the effective potential is chosen as the global minimum. Although there may be local minima that are CP-violating we find that in the global minimumω 3 always vanishes up to numerical fluctuations at both T = 0 and T = T c . Hence we will not comment on it any further. In order to determine the critical temperature T c where the phase transition takes place, 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. The temperature T c is then set to the lower bound of the final interval. We exclude parameter points that do not satisfy |v(T = 0) − 246.22 GeV| ≤ 2 GeV, and parameter points where no PT is found for T ≤ 300 GeV 10 . Moreover, we only retain parameter points with T c > 10 GeV.
The complete calculation and implementation was checked against an independent calculation in Mathematica. Profiting from significant speed-up, the implementation above was used for the results presented in this work.

Constraints and Parameter Scan
We determine the value of ξ c only for those points that are compatible with theoretical and experimental constraints. In order to obtain viable data sets we use ScannerS [61,62] to perform extensive scans in the 2HDM parameter space and check for compatibility with the constraints. The program verifies if the tree-level potential is bounded from below by applying the conditions given in [63] and checks for tree-level perturbative unitarity as described in [64]. In the CP-conserving 2HDM investigated here, the requirement that the neutral CP-even tree-level minimum is the global one is tested through a simple condition [65]. The consistency with the EW precision constraints has been checked through the oblique parameters S, T and U [66] by applying the general procedure for extended Higgs sectors as described in [67,68] and demanding for compatibility with the SM fit [69] within 2σ, including correlations. Constraints applied to the charged sector of the 2HDM are based on results from the measurement of R b [70,71] and B → X s γ [71][72][73] including the recent calculation [74] that enforces m H± > 480 GeV (40) in type II models. In type I models the bound is much weaker and more strongly dependent on tan β. Note, that the results from LEP [75] and the LHC [76,77] 11 require the charged Higgs mass to be above O(100 GeV) depending on the model type. For the check of the compatibility with the Higgs data we need the Higgs production cross sections normalised to the corresponding SM values and the Higgs branching ratios. The latter have been computed with HDECAY version 6.51 [79][80][81]. This program includes the state-of-the-art higher order QCD corrections and off-shell decays. The Higgs production cross sections through gluon fusion and b-quark fusion at the LHC have been obtained at NNLO QCD from an interface with SusHi [81,82] and normalised to the corresponding SM value at NNLO QCD. The cross section ratio for associated production with a heavy quark pair has been taken at LO. In the ratio involving CP-even Higgs bosons the QCD corrections drop out. This is not the case for the pseudoscalar. For associated production with top quarks the cross section is very small. The associated production with bottom quarks can be important for large values of tan β. However, here the QCD corrections in the associated production of the pseudoscalar with the bottom quark pair almost cancel against those of the SM counterpart due to the nearly realised chiral limit for the small b-quark masses. The remaining processes through gauge boson fusion and Higgs radiation off a W ± or Z boson only apply for a CP-even Higgs boson so that here the QCD corrections drop out when normalised to the SM cross section. Since not all EW corrections have been provided for the 2HDM so far they are consistently neglected in all production and decay processes. Agreement with the exclusion bounds from LHC Higgs searches has been tested with HiggsBounds [83]. Compatibility with the observed signal of the 125 GeV Higgs boson has been verified by calculating the reduced signal strengths and checking against the two times one sigma bounds in the six parameter fit of [84]. Further details on the various checks can be found in [62]. 12 For the minimisation procedure we only use parameter points that are in agreement with the described theoretical and experimental constraints. In order to find viable parameter points we perform a scan in the 2HDM parameter space given by the input parameters Eq. (16). The SM VEV given by the Fermi constant G F through v = 1/ √ 2G F , has been fixed to v = 246. 22 GeV .
The mixing angle α is varied in the theoretically allowed region, i.e.
In all scans, one of the masses of the CP-even Higgs bosons has been fixed to [85] m h 125 = 125.09 GeV .
This is the Higgs boson we identify with the SM-like Higgs boson discovered at the LHC, and we denote it by h 125 . We performed two separate scans for the cases where the lighter or the heavier of the two CP-even Higgs bosons is identified with the SM-like Higgs, i.e. m h = m h 125 and m H = m h 125 , respectively. The scan ranges for the remaining parameters are given in Table 2 in case of type I and in Table 3 for type II. In our scans we required the neighboring non-SM-like    experiments. The parameter m 2 12 is constrained by the tree-level global minimum condition to be positive. The upper limits on tan β and m 2 12 have been set by choice, but as we observe later, most of the points compatible with the constraints and a strong PT are found for rather small tan β so that the chosen upper limit does not pose a strong constraint. Type-specific choices for the ranges are the lower bound on tan β in type I and the lower bound on m H ± in type II. They have been chosen such that they already leave out part of the parameter space that is excluded by the constraints from B → X s γ measurements. Moreover, in type II the lower bound on m A in the second set, where H ≡ h 125 , is motivated by the fact that fulfilling constraints on the oblique parameters requires one Higgs to be in vicinity of the charged Higgs boson. In the second set this can only be the pseudoscalar Higgs A.
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, the fine structure constant is taken at the Z boson mass scale [86], The massive gauge boson masses are chosen as [86,87] M W = 80.385 GeV and M Z = 91.1876 GeV , the lepton masses as [86,87] m e = 0.510998928 MeV , m µ = 105.6583715 MeV , m τ = 1.77682 GeV , For consistency with the ATLAS and CMS analyses the on-shell top quark mass m t = 172.5 GeV (49) has been taken as recommended by the LHC Higgs Cross Section Working Group (HXSWG) [87,89]. The charm and bottom quark on-shell masses are [87] m c = 1.51 GeV and m b = 4.92 GeV .
We take the CKM matrix to be real, with the CKM matrix elements given by [86] 13

Results
We now turn to the presentation of our results. We will discuss the specific features of the 2HDM parameter space that is compatible with the theoretical and experimental constraints and at the same time provides a strong first order PT. We will show results both for the type I and the type II 2HDM. For comparison with results in the literature, we show one plot where we have applied the 'Parwani' method in the treatment of the thermal masses, cf. subsection 2.2. In the remaining discussion, however, we apply the 'Arnold-Esinosa' method for reasons discussed in [38] and alluded to in section 2.2. We will discuss scenarios where the lighter of the CP-even Higgs bosons is identified with the discovered Higgs boson, i.e. h ≡ h 125 , and where H ≡ h 125 .
For the interpretation of our results some general considerations on first order PTs are in order. The value of ξ c is proportional to the couplings of the light bosonic particles to the SM-like Higgs boson, and it decreases with the Higgs boson mass [51]. The additional Higgs bosons in the 2HDM spectrum allow for large trilinear bosonic couplings, in contrast to the SM, where bosonic couplings are only due to the weak gauge couplings between the Higgs boson and the EW gauge bosons. In the 2HDM, the second CP-even Higgs boson with a non-vanishing VEV contributes to the PT and can reduce its strength if H is not light enough. A strong electroweak PT therefore requires H either to be light or to have a vanishing VEV. The latter corresponds to the alignment limit where only one of the physical Higgs bosons has a VEV [90]. Previous investigations suggest that a first-order PT prefers a scalar spectrum, which is not too heavy [27,28,31], or else a large mass splitting between the heavy scalars [28,32]. In the type II 2HDM the requirement of a light Higgs spectrum puts some tension on the model, as compatibility with the EW precision tests requires one of the non-SM-like neutral Higgs bosons to be close to m H ± . Charged Higgs masses below 480 GeV on the other hand are already excluded by B → X s γ.

Type I: Parameter Sets with h ≡ h 125
We start with the analysis of the results in the 2HDM type I. Figure 1 shows in the m A versus m H plane all parameter points that pass the applied constraints, for scenarios where h ≡ h 125 . The coloured points are those for which we obtain a strong first order PT, i.e. where ξ c ≥ 1. In the treatment of the thermal masses we have applied the 'Parwani' method (left plot) in order to compare to the results of [28], where the 'Parwani' method was applied. In the right plot we show the results for the 'Arnold-Espinosa' method, which we will use in the remainder of the   discussion. As can be inferred from the plots, in the 2HDM type I first order PTs are still possible taking into account the up-to-date LHC Higgs data and all theoretical constraints on the 2HDM Higgs potential. The comparison of the left and right plot, however, also shows that the results obtained for ξ c are significantly different when the two different approximations in the treatment of the thermal masses are applied. Overall, the regions in the parameter space compatible with ξ c ≥ 1 are smaller when the 'Arnold-Espinosa' method is applied. Furthermore, the maximum values of ξ c that can be obtained with the 'Parwani' method are by a factor five larger than those obtained with the 'Arnold-Espinosa' method. Working with a one-loop effective potential only, the 'Parwani' method cannot be applied consistently, which is reflected in the very different results for both methods. Note also that the unrealistically large values for ξ c obtained in the 'Parwani' method imply very low critical temperatures T c where the phase transition takes place. This again questions the way the thermal masses are implemented so that the results of the 'Parwani' method have to be taken with care. In the following, we will only show plots for the 'Arnold-Espinosa' method.
In order to examine how the requirement of a strong first order phase transition translates into LHC Higgs phenomenology we show in Fig. 2 the mass differences between the non-SM-like Higgs bosons. The left plot shows the frequency of the points that pass the constraints. The right plot displays the frequency of the points when additionally a strong EW phase transition is required. As can be inferred from the left plot, the EW precision tests, namely the measurement of the ρ parameter, force the mass differences between the charged Higgs boson and at least one of the non-SM-like Higgs bosons to be small and strongly favour mass spectra where all of the non-h 125 masses are close to each other. The requirement of a strong EW phase transition, however, favours scenarios where the pseudoscalar mass is close to m H ± with a larger mass gap relative to a lighter H. In Fig. 3 we display the relative frequencies in the m A versus m H plane for all point passing the constraints (left) and for those points which additionally fulfill ξ c ≥ 1 (right). The comparison of the two plots shows that the requirement of a strong PT favours a mass spectrum where the heaviest Higgs bosons A and H ± have masses around 400-500 GeV and m H ≈ 200 GeV. H, which acquires a VEV, should be light, so that the strength of the PT is not reduced by a heavy H. Consequently m 2 12 is small 14 , which means that the strength of the phase transition is governed by the quartic couplings λ 4 and λ 5 , cf. also [27]. The next important mass configuration is given by scenarios where again the mass gap between A and H is large, but now overall pushed to higher mass values, i.e. m H ≈ m H ± and m A − m H ≈ 350 GeV, cf. Fig. 2 (right). Since h ≡ h 125 and hence sin(β − α) ≈ 1, the coupling g ZAH ∼ sin(β − α) between A, Z and H is significant. The requirement of a strong PT prefers scenarios where the decay A → ZH is kinematically allowed so that this decay can become important. These scenarios can be searched for at the LHC, as has been found earlier in [28] and proposed by the authors as possible benchmark scenarios. Still, Fig. 2 demonstrates that also scenarios are compatible with ξ c ≥ 1 where all three non-SM-like Higgs bosons are close in mass or where the decay H → AZ is possible, i.e. where m H > m A and either m H − m H ± ≈ 0 or m A − m H ± ≈ 0. While our results confirm earlier results in the literature [28,32], our results also show that a decay A → ZH is not unique for a 2HDM type I featuring a strong first order PT. The majority of the scenarios we find is very close to the alignment limit, i.e. sin(β − α) ≈ 1 with tan β close to its smallest possible value of about 2.5. While this is a feature resulting already from the constraints applied, the requirement of a strong PT overall pushes the Higgs rates towards SM values, as can be inferred from Fig. 4. It shows in grey the distribution of the Higgs signal strengths for the scenarios passing the constraints and in colour the scenarios that are additionally compatible with a strong PT. The colour code indicates the strength of the PT. The left plot shows µ V /µ F versus µ γγ and the right one µ τ τ versus µ V V . Here µ F denotes the fermion initiated cross section (gluon fusion and associated production with a heavy quark pair) of the SM-like Higgs boson (h 125 ) normalised to the SM, and µ V the normalised production cross section through massive gauge bosons (gauge boson fusion and associated production with a vector boson). The value µ xx is defined as where H SM is the SM Higgs boson with mass 125 GeV. The left plot shows that for µ V /µ F close to 1, enhanced signal rates in the photonic final states with µ γγ of up to about 1.5 are still allowed. However, including the requirement for a strong first order PT the possible range of an enhanced µ γγ is strongly restricted down to µ γγ ≈ 1.1. On the other hand, the limits on the τ or gauge boson final states are not as significantly changed, as can be inferred from Fig. 4 (right).

Type II: Parameter Sets with h ≡ h 125
We now turn to the discussion of the compatibility of the 2HDM type II with the requirement of ξ c ≥ 1 for scenarios with h = h 125 . Figure 5, which displays the values of ξ c for all parameter points compatible with our constraints, shows that also in the 2HDM type II there are scenarios allowing for a strong first order PT. The constraints from B-physics observables and the EW precision tests raise the mass scale for m H ± and at least one of the non-SM-like Higgs bosons to higher values. For m A < ∼ 350 GeV we only find few scenarios compatible with the experimental constraints. The pseudoscalar with m A < ∼ 350 GeV has a significant branching ratio into Zh (up to 10%). This final state has been searched for by the LHC experiments. The resulting exclusion limits severely constrain this parameter region so that there the amount of points compatible with the experimental constraints is substantially smaller than above the top quark pair threshold where A dominantly decays into tt. 15 When additionally a strong first order PT is required, the mass region 130 GeV < ∼ m A < ∼ 340 GeV is completely excluded. As can be inferred from the plot, for these values of m A the heavy Higgs mass ranges between ∼ 450 and 700 GeV. In this range the occurrence of a strong first order PT is strongly limited by deviations from the exact alignment limit at the per mille level. The small portions of the VEV assigned to H by these tiny deviations already suppress the strength of the PT strongly due to the large m H . Once m H > 650 GeV, even for the parameter points extremely close to the alignment limit, the H mass is finally too heavy to allow for a strong PT. The restrictions for m A < ∼ 120 GeV on the other hand are less severe, as there are less experimental studies in this mass region so that we have more points allowed by the experimental constraints. This increases the chances of finding a strong first order PT and explains why we have some coloured points for m A < ∼ 120 GeV.  The implications of the requirement of a strong first order phase transition in the type II 2HDM for LHC phenomenology can be read off Fig. 6. Scenarios with all non-SM-like Higgs masses being While again A → ZH is a typical decay that is possibly realised for strong first order PTs, the non-discovery of such a decay does not exclude ξ c ≥ 1 as other scenarios can be realised as well. We find that scenarios with m A > ∼ 460 GeV are preferred and namely those scenarios that are located in the alignment limit with tan β ≈ 1. This is, however, not due to the first order PT but already found by only imposing the theoretical and experimental constraints.
In the type II 2HDM there are parameter regions compatible with the experimental constraints where the coupling of the h 125 to the massive gauge bosons is of opposite sign with respect to the coupling to down-type fermions. This wrong-sign regime [62,91,92] has interesting phenomenological implications like the non-decoupling of heavy particles [58,91]. Future precision measurements of the signal rates will allow to constrain or exclude this parameter region [62,91,93]. The question arises to which extent the requirement of a strong PT is able to restrict the wrong-sign regime. ure 7 (left) displays µ V /µ F versus µ γγ . Among the grey points, which show the scenarios passing all constraints, the outliers in the left bottom corner of the plot correspond to the wrong-sign regime. The coloured points fulfill ξ c ≥ 1 and show that a strong PT strongly disfavours the wrong-sign regime. This can also be observed in Fig. 7 (right) where the distribution of µ τ τ versus µ V V is displayed. The wrong-sign regime scenarios are given by the outliers in the upper left corner of the plot. This behaviour can be understood by the fact that the VEV H of the heavy CP-even Higgs normalised to the SM VEV for h ≡ h 125 is given by In the wrong-sign regime non-zero values of cos(β − α) are still compatible with the data. This means that H can take a significant fraction of the VEV and drive the PT. If H is not light enough, the PT is reduced to values ξ c < 1. We also observe that the maximum value of µ γγ is reduced from about 1.46 to about 1.38.

Type I: Parameter Sets with H ≡ h 125
We now investigate scenarios where the heavier of the two CP-even Higgs bosons is the SM-like Higgs boson, i.e. H ≡ h 125 . Figure 8 displays in the m A versus m h plane in grey all points passing the constraints and in colour all parameter points also compatible with ξ c ≥ 1. First, we observe that independent of the strength of the PT, there are only few scenarios with m h < ∼ 65 GeV. This is due to the fact that the decay H ≡ h 125 → hh can change the total width of h 125 such that its branching ratios into SM final states lead to signal rates not compatible with the LHC data any more. In this mass region there are hardly any points with ξ c ≥ 1. The requirement of ξ c ≥ 1 also restricts the mass of the pseudocscalar to the region 280 GeV < ∼ m A < ∼ 480 GeV, with the exception of a few outliers. The strongest PTs are reached for larger m A , close to 480 GeV. Figure 9 displays the distribution of the masses for A and H ± after applying all constraints (left) and when in addition ξ c ≥ 1 is demanded (right). With the exception of a few outliers, the strong PT restricts the mass region of the charged Higgs boson to 300 GeV < ∼ m H ± < ∼ 480 GeV. As we demand the heavier of the two CP-even Higgs bosons to be light, the mass scale m 2 12 , which determines its mass, cannot be too large. For the PT to be strong we need large quartic couplings. Since λ 2 , which enters m H , must not be large, we are left with λ 4 and λ 5 driving the PT, as can be inferred from the rather large mass values for A and H ± , namely the mass gap between m H and m A,H ± . When, on the other hand, the masses of the heavy Higgs bosons A and H ± become larger by increasing the involved quartic couplings the interplay of Higgs self-couplings and masses reduces ξ c again. We conclude that a strong PT in the 2HDM type I with two light CP-even Higgs boson excludes heavy Higgs bosons above about 500 GeV and enforces a mass gap between m A ≈ m H ± and m H . The decay A → HZ, however, is suppressed because of sin(β − α) ∼ 0 for H ≡ h 125 . The decay A → hZ on the other hand, is allowed. For pseudoscalar masses above the top pair threshold, it competes, however, with the decay A → tt.
The implications of a strong PT for the Higgs data are shown in Fig. 10. There are practically no points any more with values beyond 0.9 for the photonic rate, although rates of up to about 1.46 are still compatible with the Higgs data. Also the decays into τ final states cannot exceed 1.11 in case of ξ c ≥ 1.

Type II: Parameter Sets with H ≡ h 125
In the 2HDM type II with H ≡ h 125 the implications of a strong PT on the mass pattern are very pronounced, as can be inferred from Fig. 11. The requirement of ξ c ≥ 1 excludes a large portion of the parameter space, which is still compatible with the applied constraints. Scenarios with m A ≈ m H ± > ∼ 480 GeV are forbidden if ξ c ≥ 1. Furthermore, very light scalars with m h < ∼ 110 GeV are not compatible with a strong PT. The tension between the requirement of light scalar masses and the wish to have a strong PT makes a strong link between baryogenesis and collider phenomenology.
Further implications for LHC phenomenology are shown in Fig. 12 where the signal rates are displayed before (grey) and after (coloured) imposing a strong PT. All scenarios with ξ c ≥ 1 are located in the correct-sign regime (given by the triangle areas in the plots), whereas the wrong-sign regime (given by the outliers) is completely excluded by a strong PT. For the Higgs measurements, this means that the observation of µ V /µ F < 1 together with µ γγ < ∼ 0.9 is excluded, as well as the observation of µ τ τ > ∼ 1.04. Furthermore, the region where µ τ τ < ∼ 0.9 because of possible decays H → hh, is excluded by the demand of a strong PT.

Conclusions
In this paper we investigated the strength of the EW phase transition in the framework of the CPconserving 2HDM. For this purpose we computed the loop-corrected effective potential at non-zero temperature including the resummation of the daisy graphs for the bosonic masses following the 'Arnold-Espinosa' method. We applied a renormalisation scheme that preserves the position of the  minimum and where both the loop-corrected masses of the Higgs bosons and the mixing angles are renormalised to their tree-level values. This is in contrast to earlier works which focus solely on the Higgs boson masses. Our renormalisation allows us to efficiently scan the whole 2HDM parameter space and test the compatibility of the model with the theoretical and experimental constraints. This is possible since our renormalisation fixes not only the Higgs mass values but, through the mixing angles, also the Higgs couplings to their tree-level values.
We performed an extensive scan in the parameter space of the 2HDM and retained only those points that are compatible with the state-of-art theoretical and experimental constraints. For these parameter points we determined the value of ξ c . Subsequently, we performed a comprehensive and systematic analysis in four 2HDM configurations: For the 2HDM type I and II, with either h or H identified with the SM-like Higgs boson, we investigated the implications of the requirement of a strong PT, i.e. ξ c ≥ 1, for LHC phenomenology. Our results can be summarised as follows: Both the 2HDM type I and type II, with either of the CP-even Higgs bosons being the SM-like Higgs boson, are found to be compatible with the theoretical and experimental constraints on the model and a strong first order PT. The strong PT, however, strongly constrains the enhanced rates into photonic final states for the 2HDM type I, and to some extent also for type II. Furthermore, in the 2HDM type II with h = h 125 the wrong-sign regime is strongly restricted by the requirement of ξ c ≥ 1. In case H = h 125 the wrong-sign regime is even excluded for a strong PT. In more detail, our results for the four different realisations of the 2HDM are: • For the 2HDM type I with h ≡ h 125 , we confirm earlier results which find that a large mass splitting between the heavy scalars is favourable for a strong PT. The preferred scenarios are the ones with m A ≈ m H ± ≈ 400 − 500 GeV and m A − m H > ∼ 200 GeV. However, we also find that scenarios with different hierarchies among the heavy Higgs bosons (but at least one of H and A nearly mass degenerate with H ± ) or with degenerate heavy Higgs bosons H, A and H ± are allowed, though much less frequent. The maximally allowed photonic rate is reduced from 1.5 to 1.1 in case of ξ c ≥ 1.
• We find in the 2HDM type II with h = h 125 that scenarios with 130 GeV < ∼ m A < ∼ 340 GeV, which are already strongly constrained by LHC searches in A → Zh 125 , are completely excluded by the requirement of a strong PT. This requirement also restricts the wrong-sign regime considerably. The maximum value of µ γγ is reduced from about 1.46 to about 1.38.
• In the 2HDM type I with two light CP-even Higgs bosons, namely H ≡ h 125 , the heavy Higgs masses cannot exceed 480 GeV, although experimentally still allowed, if ξ c ≥ 1. Furthermore, this enforces a mass gap of about 155 GeV between the heavier and the lighter Higgs bosons. A strong PT is found to exclude almost completely scenarios with m h < ∼ 65 GeV. It strongly reduces µ γγ from 1.46 down to 0.9 (with very few exceptions) and limits the τ final state rates to values below 1.11.
• In the 2HDM type II with H ≡ h 125 the tension between a light CP-even Higgs mass spectrum and a strong PT excludes large portions of the parameter space. The observation of heavy Higgs bosons with masses above 480 GeV or of a light Higgs boson with a mass below 110 GeV is excluded by the requirement of a strong PT. Furthermore, simultaneously reduced values of µ V /µ F < 1 and µ γγ < ∼ 0.9 are not compatible with ξ c ≥ 1, nor values of µ τ τ > ∼ 1.04. The reason is that the wrong-sign regime is excluded by a strong PT. The requirement of a strong PT also excludes parameter regions with reduced µ τ τ < ∼ 0.9, resulting from Higgs-to-Higgs decays H → hh, which contribute to the total width of the SM-like Higgs boson.
Our results show that there is a strong interplay between the requirement of successful baryogenesis and LHC phenomenology. The realisation of a strong EW phase transition leads to testable consequences for collider phenomenology. The systematic investigations performed in this work serve as basis for further analyses of the LHC phenomenology of 2HDM models featuring a strong EW phase transition.
The fermion mass at T = 0 is given by the tree-level VEV v k of the doublet Φ c k as

A.2 Gauge Boson Masses
The longitudinal gauge bosons get a Debye correction to their mass matrix. The masses including the thermal corrections, denoted in section 2.2 by m, in terms of the field configurations ω k are given by where g and g denote the SU (2) L and U (1) Y gauge couplings, respectively, and Again, the physical masses are obtained for ω i ≡ω i , and at T = 0 we recover the well-known relations for the physical gauge boson masses (v 2 = v 2 1 + v 2 2 = i=1,2,3ω

A.3 Masses of the Higgs Bosons
The tree-level relations for the mass matrices of the Higgs bosons in the interaction basis in terms of the ω k are obtained by differentiating the tree-level Higgs potential V tree Eq. (2) twice with respect to the real interaction fields φ i ≡ {ρ 1 , η 1 , ρ 2 , η 2 , ζ 1 , ψ 1 , ζ 2 , ψ 2 } and replacing the fields with their classical constant field configurations φ c i ≡ {0, 0, 0, 0, ω 1 , 0, ω 2 , ω 3 } , leading to the mass matrix (M) ij = 1 2 The physical masses are given by the field values in the global minimum of the potential where ω k ≡ω k , which at T = 0 reduces toω 1,2 | T =0 = v 1,2 andω 3 | T =0 = 0. Because of charge conservation the mass matrix of Eq. (65) decomposes into a 4 × 4 matrix M C for the charged fields ρ 1 , η 1 , ρ 2 , η 2 and a 4×4 matrix M N for the neutral states ζ 1 , ψ 1 , ζ 2 , ψ 2 . In the CP-conserving 2HDM the neutral CP-even and CP-odd fields do not mix so that the latter matrix further decomposes into two 2 × 2 matrices, one for the CP-even Higgs states ζ 1,2 and one for the pseudoscalar states ψ 1,2 . We introduce the following definitions where we take for the top and bottom quark masses at zero temperature, m t,b (T = 0), the input values given in Eqs. (49) and (50). The masses of the charged Higgs boson and the charged Goldstone boson including the thermal corrections are then given by M N 12 = λ 5 ω 2 ω 3 (83) M N 13 = −m 2 12 + (λ 3 + λ 4 + λ 5 ) ω 1 ω 2 (84) M N 23 = λ 5 ω 1 ω 3 (86) M N 24 = −m 2 12 + λ 5 ω 1 ω 2 (87) The physical masses at T = 0 are recovered after replacing the ω k with the VEVs at T = 0. In particular, the Goldstone masses become zero in the Landau gauge, in which we are working.