Nonperturbative Analysis of the Electroweak Phase Transition in the Two Higgs Doublet Model

We perform a nonperturbative study of the electroweak phase transition (EWPT) in the two Higgs doublet model (2HDM) by deriving a dimensionally reduced high-temperature effective theory for the model, and matching to known results for the phase diagram of the effective theory. We find regions of the parameter space where the theory exhibits a first-order phase transition. In particular, our findings are consistent with previous perturbative results suggesting that the primary signature of a first-order EWPT in the 2HDM is $m_{A_0}>m_{H_0} + m_Z$.


I. INTRODUCTION
Accounting for the observed baryon asymmetry in the present universe is a major unsolved problem in cosmology. One of the leading candidates for a viable mechanism, electroweak baryogenesis (EWBG) [1], suggests that the asymmetry originates from the electroweak phase transition (EWPT) in the early universe. According to the Sakharov conditions [2] the transition would have to be of first order, accompanied by a sizable violation of CP-symmetry. Unfortunately, these conditions immediately rule out EWBG within the minimal Standard Model (SM), as it was shown in the 1990's that the SM EWPT is a crossover [3][4][5][6]. Furthermore, CPviolating effects within the SM are heavily suppressed already at temperatures of order 1 GeV, let alone at the EW transition temperature of O(150 GeV) [7][8][9].
Independently of the question of the baryon asymmetry of the universe, a host of beyond the Standard Model (BSM) theories have been proposed to solve open problems in physics. Determining whether BSM theories can produce a first order EWPT and thus facilitate EWBG is nontrivial: quantitatively reliable conclusions about the nature of the phase transition typically require a nonperturbative-in practice lattice field theory-approach, which has been deemed unmanageable for large parameter spaces. Because of this difficulty, analyses based on the finite-temperature effective potential have become the standard approach [10][11][12][13][14][15][16]. Such studies can, however, have considerable uncertainties: in one study [17], an error in excess of 10% in the critical temperature was seen comparing the 1-loop effective potential to nonperturbative results. The order of the phase transition can also differ.
One fruitful approach is the derivation of threedimensional high-temperature effective theories via dimensional reduction (DR). This method was originally applied to the SM in Ref. [18], and non-perturbative simulations were carried out in [3][4][5]; unfortunately these did not reveal the potential for a first order EWPT (see however [19]). Recently, the same method was applied to the SM accompanied by a real singlet [20].
In this paper, we turn our attention to one of the most widely studied BSM models, namely the two Higgs doublet model (2HDM), where the SM is augmented with an additional Higgs doublet (see Ref. [21] for a review, and Refs. [22][23][24] for earlier work on DR in the 2HDM). We derive a three-dimensional high-T effective theory, studying a region of parameter space where the effective theory has the same form as that of the Standard Model, similar to Ref. [25]. This reduces the problem of determining the phase diagram of the theory to mapping its parameter space to that of the SM effective theory. Equipped with the analysis of [3][4][5], we can discover interesting and phenomenologically viable regions of parameter space where the EWPT is first order. Our work corroborates key findings of perturbative studies of EWBG in the 2HDM.

II. DIMENSIONAL REDUCTION OF THE 2HDM
The 2HDM is an extension of the SM featuring an additional Higgs doublet. Our four-dimensional original theory can therefore be described by the schematic action suppressing the presence of counterterm and ghost contributions. The field content of the theory includes SU(3) c , SU(2) L and U(1) Y gauge fields, two scalar doublets φ 1 and φ 2 , as well as all fermions present in the SM. In our present treatment, we will consider only one quark flavor in the Yukawa sector, namely the top, since it has the largest coupling to the Higgs field. The top quark couples to one doublet only (by convention φ 2 ), and we have not yet committed to a specific type of 2HDM (Type I or II).
We will content ourselves with explicitly writing down only the BSM part of our model, i.e. the extended scalar sector. This has Lagrangian with usual covariant derivative D µ . The potential has the form In general, C(P) symmetry is broken when λ 5 or µ 2 12 are complex; we have discarded so-called hard CP-breaking terms, often parametrised by λ 6,7 . For a detailed discussion see Refs. [21,26].
The first three-dimensional effective theory, obtained by integrating out the hard scale πT from the theory (see e.g. Ref. [20] for details of the general procedure), has schematic form where we have again suppressed the ghost and counterterm contributions. The field content is now SU(2) L and U(1) Y gauge fields; two Higgs doublets; and temporal scalar fields A a 0 , B 0 , C α 0 . The fermions are integrated out and the SU(3) c gauge fields can be neglected [20]. The fundamental scalar sector remains similar to the full theory, where r = 1, 2, 3 is summed over. In the second step of DR, the heavy temporal scalar fields are integrated out. The final step is to notice that φ 1 and φ 2 mix when µ 2 12 = 0, and near the phase transition one mass eigenvalue will generically be small while the other is large. This observation-specific to the 2HDM-allows us to integrate out the heavy mode and study the phase transition with only one scalar field coupled to the gauge fields. Our final effective theory therefore becomes Now φ is the remaining light φ 1 -φ 2 mixing mode, and the parameters of the theory areμ 2 3 ,λ 3 and the 3-d gauge couplingsĝ 3 andĝ 3 for the U(1) Y and SU(2) L interactions, respectively. As in the original analysis of Refs. [3,18], we omit all non-perturbative effects related to the U(1) Y gauge field.
The main task of the DR process is to perturbatively match the parameters of the original 4-d theory, Eq. (2), to those of the final effective theory, Eq. (7). This is accomplished by demanding that the effective theory reproduces the static Green's functions of the original theory at large distances R 1/T . This, in turn, results in a number of matching relations from which the effective theory parameters are solved. This procedure will be presented in Ref. [26]; the results are summarised in the Appendix.
As discussed above, the effective theory of Eq. (7) has the same form as that of the SM, studied in Refs. [3][4][5]. This allows us to adopt existing numerical results for the strength of the phase transition, and determine the form of the phase diagram through our matching procedure alone.
The validity of the DR can be quantified by estimating the effect of neglected dimension-6 operators. While it is difficult to comprehensively evaluate their effect, one can evaluate the change in the vacuum expectation value (vev) of the Higgs field in the effective theory caused by the (φ † φ) 3 operators. In Eq. (201) of Ref. [18], it was shown that in the SM the dominant neglected contribution comes from the top quark; its effect is about one percent. We therefore estimate the effect of new BSM contributions to these operators by comparing their magnitude to the contribution from the top quark [see Eqs. (A34,A35) in the Appendix]. This is only a rough estimate: we consider only contributions coming from the first step of DR when the superheavy modes are integrated out. Furthermore, we do not consider other dimension-6 operators or the final (φ † φ) 3 operators appearing in the effective theory at the light scale after integrating out the heavy mode in the diagonal basis.
Finally, although the parameter matching is perturbative, the study of the 3-d phase diagram is non-perturbative and-within the limitations of lattice methods-exact. Hence, the existence and strength of a phase transition is assured. We emphasise that the main advantage of our approach lies in a proper handling of the infrared physics, which causes trouble in traditional perturbative studies of the EWPT. Resummations are performed when the superheavy and heavy scales are integrated out perturbatively, and the problematic light modes are treated non-perturbatively on the lattice. However, the mapping to precise values of the 4-d parameters, where this phase transition occurs in the 2HDM, is limited by the accuracy of the perturbative truncation. We organise the expansion in terms of the gauge coupling g, and perform the DR up to O(g 4 ).
Thus the calculation is carried out at the one-loop level for quartic couplings, and the two-loop level for mass parameters. This exceeds the accuracy used in the perturbative calculations of (for instance) Ref. [27] (but see Ref. [28] for a recent two-loop perturbative calculation in the inert doublet model).

III. SCANNING THE PARAMETER SPACE
The phase diagram of the dimensionally-reduced theory can be mapped using the dimensionless parameters x ≡λ 3 /ĝ 2 3 , y ≡μ 2 3 /ĝ 4 3 . It is known that within this theory the EWPT occurs near y 0; the Higgs mass parameter becomes negative at this point. In Refs. [3][4][5], it was found that the transition is first order for x 0.11, and strongly so for x 0.04.
We therefore search for areas of 2HDM parameter space that map onto regions of the 3-d effective theory with x < 0.11 and y 0. Since there are ten real parameters in the 4-d theory and only three in the 3-d one, inverting the mapping process is not unique. We need to perform scans of the 2HDM parameter space, where we will be guided by the results of Ref. [29] that combine phenomenological constraints with a one-loop resummed perturbative determination of the effective potential. Another recent treatment is found in Refs. [30,31].
A uniform scan through a 10-dimensional parameter space is computationally expensive; we must therefore make some simplifying assumptions. We take all parameters of the 2HDM to be real, setting Im(λ 5 ) = 0, Im(µ 2 12 ) = 0. This eliminates extra CP violation in the model, which would be crucial for baryogenesis. However, the effect of these imaginary parts on the order and strength of the transition is expected to be negligible; the CP-violating phase must necessarily be small due to EDM constraints [32][33][34].
We then reparametrise the model following Ref. [29], applying tree-level relations between the MS parameters and physical quantities. The masses of the CP-even scalars are denoted by m h = 125 GeV and m H0 ; that of the CP-odd scalar by m A0 ; and that of the charged scalar by m H ± . We also employ two angles α and β: α parametrises mixing between the CP-even states, while β is related to the ratio of the vevs tan(β) ≡ ν2 ν1 . Here, ν 1 and ν 2 are the vevs for φ 1 and φ 2 respectively with ν 2 1 + ν 2 2 = ν 2 and ν = 246 GeV. Finally, there is the squared mass scale M 2 ≡ µ 2 (tan(β) + 1/ tan(β)), where we treat µ 2 ≡ −Re µ 2 12 as an input parameter. The relations between the physical states and gauge eigenstates can be obtained from Ref. [29].
We also fix m H ± = m A0 , since EW precision tests require the mass of the charged Higgs to be roughly degenerate with either H 0 or A 0 [35,36]. Furthermore, we work in the alignment limit by setting cos(β − α) = 0. In this limit, the CP-even scalar h couples to SM particles exactly like the SM Higgs. We investigate relatively few values for tan(β), whereas we perform a more exhaustive scan in a three-dimensional parameter space spanned by m H0 , m A0 , and µ 2 . At each point, we require that treelevel stability and unitary constraints be satisfied; for details, see Ref. [26]. Note that Ref. [27] imposes stronger conditions than these, so we will overestimate the num-ber of physical points. Furthermore, for the DR to be valid, the tree-level mass parameters µ 11 , µ 22 and µ 12 should be comparable to the Debye mass m D ∼ gT near the phase transition. This sets an upper bound for the input parameter µ 200 GeV. Finally, we verify that in the effective theory the other doublet really is heavy near the phase transition, so it is justified to integrate it out.

IV. RESULTS
Following our scanning protocol outlined above, we fix tan(β) and scan in the two scalar masses m H0 and m A0 between 137.5 and 562.5 GeV at spacings of 6.25 GeV; a total of 4624 points. A dense scan in µ is then carried out for each pair, from 10 to 150 GeV at intervals of 2.5 GeV for a total of 56 values. In all, each of our fixed-tan(β) plots results from scanning approximately 260 000 points. The upper limit on µ is chosen to ensure that the DR is valid, as explained above.
We first check whether each point is physical, according to our criteria. If so, we then perform the DR for evenly-spaced temperatures between 80 and 200 GeV, at intervals of 20 GeV. This allows us to find the value of x when y = 0-on the critical line-by interpolation. We then use x to characterise the phase transition.
We take 0.0 < x < 0.11 as the indicator of a firstorder EWPT. The upper limit comes from previous nonperturbative studies. At small but positive x, the DR no longer works well, and at negative x it fails as the potential of the dimensionally reduced system is no longer bounded from below; one must include dimension-6 operators in the three-dimensional theory.
Combining different values of µ, we indicate the relative number of points with a first-order phase transition as a heat map in Fig. 1, for three separate values of tan(β). The majority of our points reside in the region m A0 > m H0 + m Z . This was claimed in Refs. [27,29] as the key signature of a first order phase transition (but see Refs. [30,31]). However, at small tan(β) we also see a considerable number of points in regions where this mass hierarchy does not hold.
In Fig. 2 we show a breakdown of the heatmap plot with fixed tan(β) = 2.0 for two values of µ. An estimation of the effect of the neglected 6-dimensional operators is included. Generally, decreasing values of x correspond to increasing importance of dimension-6 terms; when the effect of these terms becomes large, the DR also breaks down. These plots also show how the lower first-order region disappears as µ increases. The lower region extends further into the domain for which the dimension-6 operators become important. We therefore have greater confidence in our results for points in the upper region.
Experimental constraints on the 2HDM parameter space depend strongly on the way in which fermions couple to the Higgs doublets. With the exception of the top quark, other Yukawa couplings have little effect on our  EWPT analysis, and we have not yet had to indicate whether we are considering Type I (all quarks couple to φ 2 ) or Type II (up-type quarks couple to φ 2 , down-type to φ 1 ). The most stringent constraints come from flavour physics, where B-decays set the bound m H ± 580 GeV for the charged Higgs mass in the Type II 2HDM [38].
Assuming that m ± H is degenerate with m A0 in accordance with EW precision tests, this rules out our regions of firstorder EWPT in Type II, but no such lower bound exists in Type I for tan β ≥ 2 [38,39].
Additional restrictions come from direct searches for neutral Higgses at the LHC [40]. For Type I, the H 0 → τ τ cross section is suppressed by cot 2 β, and our choices of tan β are within current experimental bounds. Finally, we have verified that the mass range we scan in is allowed by measurements of the h → γγ decay [41], as well as the relatively recent search for A 0 → Zh processes [42].
Having not scanned in the hidden-Higgs region where constraints from charged-scalar searches become important [43], we conclude that our first-order EWPT regions are currently not ruled out by experiments if a Type I 2HDM is assumed.

V. DISCUSSION
It is a shortcoming of present-day particle cosmology that it is still impossible to reliably determine the nature and strength of the EWPT for a given BSM scenario. This information would be valuable not only for EWBG, but also in the context of gravitational wave physics, as a first-order EWPT would leave an imprint in the sensitivity range of the future LISA mission and other proposed space-based gravitational-wave detectors [37].
We have taken a step towards tackling this challenge, studying the mapping of the phase diagram of one viable BSM theory, the 2HDM. Our results concern the EWPT in the alignment limit cos(β − α) = 0. Our work so far supports the idea that the primary signature of a strong first order transition in this theory is indeed m A0 > m H0 + m Z , as suggested by Refs. [27,29].
We perform a perturbative construction of an effective 3D theory, going beyond previous DR studies in the 2HDM. We then apply non-perturbative results in our scan of the parameter space. Our procedure avoids the IR problems that usually plague perturbative expansions at high T .
The techniques discussed in this paper can be applied, with suitable modifications, to a host of other models where a substantial region of parameter space can be mapped onto the three-dimensional theory of the minimal Standard Model. In the 2HDM, nearly the entire parameter space is accessible in this way. For other models, further simulations will be required to fully explore the phenomenologically interesting parameter space.
In the future, our aim is to perform a thorough comparison of perturbative and non-perturbative results in this model. Similar projects to study the EWPT and benchmark the accuracy of perturbation theory are already underway in the Standard Model augmented by a real singlet [44] or triplet field [45]; the EWPT has been perturbatively analysed for the former in Refs. [46,47], and for the latter in Ref. [48]. In this Appendix we collect the matching relations between the full four-dimensional theory and effective theories. A detailed derivation can be found in Ref. [26].

Three-dimensional effective theories
We denote the fields of the effective theories with the same symbols as those of the four-dimensional theory, even though their normalisation is different and will affect the mapping between full and effective theories. These normalisations between 4-d and 3-d fields have been listed below.
The schematic form of classical Lagrangian density of the effective theory was given in Eq. (4) of the main paper. The temporal part reads Here the covariant derivative of an isospin triplet reads D r A a 0 = ∂ r A a 0 + g 3 abc A b r A c 0 , and for the temporal gluon C α 0 ordinary derivative is used instead of covariant derivative as gluons are discarded for only contributing at a higher order [20]. After the heavy temporal scalars have been integrated out, their effects are encoded by the parameters and fields of a new theory where the parameters are denoted with a bar asḡ 3 ,ḡ 3 ,μ 2 11,3 etc. In this theory, the phase transition takes place close to a point where the mass matrix has zero eigenvalue, and then generically in the diagonal basis the other mass parameter is heavy. By performing a unitary transformation, one can diagonalise the scalar potential.
Denoting Ω ≡ (μ 2 11,3 −μ 2 22,3 ) 2 + 4μ 2 * 12,3μ 2 12,3 , this transformation reads . (A3) The mass parameters in the diagonal basis read Generally µ 2 θ is heavy when µ 2 φ is light, and therefore the field θ can be integrated out. The scalar self-couplings in the diagonal basis are given by The scalar potential in the diagonal basis reads where φ and θ are light and heavy fields, respectively. When the heavy doublet θ has been integrated out, the final effective theory is same as in that of the SM, as given in Eq. (6) of the main paper.

Matching relations and normalisations of fields
Our calculations are carried out in the MS scheme. We use the following notation: where Λ is the renormalisation scale of the 4-d theory and γ is the Euler-Mascheroni constant. The normalisations relating three-and four-dimensional fields read For dimensional reduction, the required ingredients include matching relations between the 4-d and 3-d theories, oneloop β-functions (to make the matching relations renormalisation scale independent) and finally the relations between MS-parameters and physical quantities. We use the tree-level relations, despite the fact that for consistent O(g 4 ) accuracy one should use the one-loop corrected relations. This would require performing one-loop vacuum renormalisation of the physical quantities. This is a non-trivial task, and is left for the future. In the special case of the inert doublet model, the one-loop vacuum renormalisation can be found in Ref. [28]. Below we list all needed matching relations, while β-functions and relations of MS-parameters and physical quantities can be found in Ref. [26] with detailed derivations and explicit, step-by-step intermediate results.
a. Integration over superheavy scale A full O(g 4 )-accurate dimensional reduction requires the evaluation of the mass parameters at two-loop and couplings at one-loop order. The results are listed below. and We emphasise that these relations are independent of the four dimensional theory renormalisation scale Λ, which can be seen by applying the β-functions to the tree-level terms. This serves as a crosscheck of the correctness of our calculation. These relations have been calculated previously in Ref. [22], with the restriction of λ 5 being real rather than complex. We have corrected two minor errors in the expressions for λ 4,3 and λ 5,3 for the terms involving λ 5 .
We estimate the validity of the dimensional reduction by evaluating to order O(g 6 ) the terms Λ 6,1 (φ † 1 φ 1 ) 3 3d and Λ 6,2 (φ † 2 φ 2 ) 3 3d , where the coefficients are given by where λ ± ≡ 1 2 (λ 3 + λ 4 ± λ 5 ). The last term in the above equation is the contribution from the top quark. Comparing the magnitude of the other terms to that term yields a rough estimate of the importance of neglected dimension-6 operators.

b. Integration over heavy scale
Here we list the matching results for the parameters in the simplified 3-d effective theories, where the heavy temporal scalars and the heavy second doublet have been integrated out. The two-loop contributions to mass parameters are highlighted with subscripts. One could give the coefficients of logarithmic contributions in terms of the running in the final theory here as well, but we have omitted this for simplicity.