The Electroweak Phase Transition in the Inert Doublet Model

We study the strength of a first-order electroweak phase transition in the Inert Doublet Model (IDM), where particle dark matter (DM) is comprised of the lightest neutral inert Higgs boson. We improve over previous studies in the description and treatment of the finite-temperature effective potential and of the electroweak phase transition. We focus on a set of benchmark models inspired by the key mechanisms in the IDM leading to a viable dark matter particle candidate, and illustrate how to enhance the strength of the electroweak phase transition by adjusting the masses of the yet undiscovered IDM Higgs states. We argue that across a variety of DM masses, obtaining a strong enough first-order phase transition is a generic possibility in the IDM. We find that due to direct dark matter searches and collider constraints, a sufficiently strong transition and a thermal relic density matching the universal DM abundance is possible only in the Higgs funnel regime.


Introduction
The simplest extension of the Standard Model (SM) that includes two SU (2) Higgs doublets is known as the inert Higgs Doublet Model (IDM). In the IDM, the extra doublet has no coupling to SM fermions and is odd under a postulated new Z 2 discrete symmetry, whereas all SM fields are Z 2 -even. Such symmetry makes the lightest Z 2 -odd particle (LOP) from the extra doublet stable and, thus, a potential weakly interacting massive particle (WIMP) dark matter candidate. The symmetry also eliminates numerous terms in the interaction Lagrangian of the model containing an odd number of extra "inert" scalars.
An additional early motivation to consider the IDM as an appealing augmentation of the SM scalar structure was to allow for a relatively heavy SM-like Higgs while remaining compatible with constraints from electroweak precision observables, and without large fine tuning [3,10]. Although this motivation has somewhat faded after the discovery of a SM-like Higgs boson at the LHC with a mass of ∼ 125 GeV [11,12], this important discovery decreases the number of free parameters in the theory by one, and places interesting and stringent constraints on the IDM phenomenology [13].
In the present study we are concerned with the nature of the electroweak phase transition (EWPT) in the IDM, and, specifically, with determining which physical parameters drive the strength of the phase transition, making it more or less strongly first-order, or second-order. This question is intimately related with the possibility to produce the observed baryon-antibaryon asymmetry in the universe at the electroweak phase transition: A strongly first-order phase transition (in a quantitative sense we shall make clear below) is a necessary ingredient to (i) achieve the necessary out-of-equilibrium conditions, occurring on the boundary of broken and unbroken electroweak phase, and to (ii) shield a baryon asymmetry captured in the broken electroweak phase region from sphaleron wash-out.
While necessary, a strongly first-order phase transition is not a sufficient condition. The CP violating sources of the SM are known to be insufficient to generate the necessary asymmetry in the number density of baryons compared to antibaryons during the electroweak phase transition. The unbroken Z 2 symmetry in the IDM precludes any new source of CP violation, and thus this model per se cannot accommodate successful electroweak baryogenesis (EWBG). However, the IDM might be in effect a good approximation at low energy of a broader construction that includes such additional CP violating sources at higher energies. Several suggestions of plausible effective higherdimensional operators have been made in the literature [14,15]. We will not discuss this aspect any more, as it falls outside the scope of this study.
The nature and strength of the electroweak phase transition in the IDM has been subject of several studies, with increasingly refined treatment of the effective potential [16][17][18][19][20]. For example, Ref. [16] utilized only the high-temperature form of the effective potential without including the zero-temperature Coleman-Weinberg terms. These were then shown to be quantitatively important for the phase transition strength in Ref. [18], where the full one-loop effective potential was used. Alternative SU (2) L representations of the inert scalar were considered in Ref. [19] where it was argued that in general, higher representations are less successful in satisfying experimental and theoretical constraints, thereby further motivating the study of the doublet case.
With the exception of Ref. [20], the primary focus has been on the Higgs funnel regime (described in more detail in Section 3). Indeed, we will confirm the findings of Refs. [17,18] that this is the only region of parameter space that can successfully saturate the DM abundance and provide a strong-enough first-order EWPT. In Ref. [20] it was emphasized that the IDM can be useful for the EWPT even if the LOP provides only a sub-leading component of DM.
In this work we go beyond previous studies by utilizing a state-of-the-art treatment of the finite-temperature and zero-temperature effective potential including renormalization group, daisy resummation improvements and one-loop model parameter determination. As we discuss in great detail in what follows, strongly first-order EWPT in the IDM requires sizable quartic couplings that enhance quantum corrections to masses. This is important in the context of DM phenomenology since DM particle production in the early universe often relies on resonance and threshold effects [13]. In addition, we also ensure that the phase transition completes by evaluating bubble nucleation rates.
Unlike previous studies which primarily utilized large numerical scans of the parameter space, here we take an orthogonal approach: we restrict our attention to a few benchmark models, motivated by the requirement of having a viable dark matter particle candidate and representing different features in the DM phenomenology. Based on these benchmarks, we then discuss how the EWPT depends on the physical model parameters. We identify the key physical inputs that drive the phase transition to the interesting regime where it is strongly enough first-order to accommodate successful electroweak baryogenesis. We will see that in all but one case the demand for a strongly first-order EWPT is in tension with either the relic abundance requirement or with experimental probes.
Our central finding is that the main driver of the strength of the phase transition is the mass difference between the lightest inert scalar and the heavier scalars. Thus, we extend the results of Refs. [16,18] to other regions of IDM parameter space. For large enough mass splittings, but for light enough heavy scalars, we find a phase transition strength (as measured by the ratio v c /T c , as we discuss in detail below) which increases with the mass splitting.
The remainder of this paper is organized as follows: In Section 2 we give a brief introduction to the IDM, thereby clarifying our conventions, discuss quantum and finitetemperature corrections to the effective potential, and outline the computation of the phase transition strength. The essential features of DM phenomenology are reviewed in Section 3. In Section 4 we study the electroweak phase transition in several benchmark models motivated by the various DM scenarios available in the IDM. We conclude in Section 5.
2 Phase Transitions in the Inert Doublet Model (IDM)

IDM at Tree-Level
The IDM is a particular realization of the general type I Two Higgs Doublet Model (2HDM) (see, e.g.,, Ref. [21] for a review) which features an additional Z 2 symmetry. The SM doublet H is even under Z 2 , while the new (inert) doublet Φ is odd. If we take Φ to have hypercharge +1/2, the most general renormalizable potential consistent with these symmetries is then given by [13]: (2.1) Conventionally, within CP -conserving Higgs sectors, the physical states are decomposed into CP -even and CP -odd scalars. One should keep in mind, however, that in the IDM there is no observable that can actually distinguish between the CP -even or CP -odd character of the inert Higgs bosons. In the absence of a vacuum expectation value (VEV) for Φ, the doublets decompose as .

(2.2)
Below we consider the thermal evolution of the effective scalar potential in the early universe. In general, spontaneous breaking of Z 2 can occur, in which case we must also include a VEV for the neutral component of Φ which we will indicate with φ.
The lightest Z 2 -odd particle is stable, and potentially provides a viable particle dark matter candidate. The Z 2 symmetry also forbids Yukawa couplings of Φ to SM fermions (assumed to be even under Z 2 ), which eliminates tree-level flavor-changing neutral currents. Either H or A can be the LOP, and since gauge interactions with SM states do not distinguish between the two, they are effectively equivalent from the standpoint of phenomenology. Below we will indicate the LOP as H, but all statements made with respect to DM phenomenology and the electroweak phase transition remain true after the replacements H → A and λ L → λ S (the latter determines the coupling of the LOP to the SM Higgs) [17].
In the electroweak vacuum, the tree-level masses of the new states are given by where λ L = (λ 3 + λ 4 + λ 5 )/2 and λ S = (λ 3 + λ 4 − λ 5 )/2. In what follows we employ the one-loop effective potential defined in the next section to fix µ 2 1 and λ 1 from the physical values v = 246.22 GeV and m h ≈ 125 GeV. The remaining parameters of the model are specified using the three physical masses m H , m A and m H ± , along with λ L and λ 2 . The masses are related to potential parameters using the one-loop relations from Ref. [13], while λ L and λ 2 are given at scale M Z .

Finite-Temperature Corrections
The effective potential at finite temperature T can be written as where V 0 , V 1 and V T are tree-level, one-loop temperature-independent and -dependent pieces, respectively. The tree-level potential V 0 has been given in Eq. (2.1). The temperature-independent one-loop correction has the Coleman-Weinberg form [22][23][24]: The sum is over all particle species coupling to the doublets; n i is the number of degrees of freedom (positive for bosons and negative for fermions), C i are renormalizationscheme-dependent constants (C i = 1/2 for transverse gauge bosons and 3/2 for everything else in the MS scheme); m 2 i (v, φ) is the field-dependent squared mass for each species. In writing the above, we have implicitly absorbed the counterterms into V 1 ; the temperature-dependent part is ultraviolet finite.
The field dependent masses in the IDM for the SM vector bosons and fermions are, respectively, and with the corresponding bosonic degrees of freedom n i = 6, 3, 2 for i = W, Z, A, and fermionic degrees of freedom The field-dependent neutral CP -even, CP -odd and charged scalar mass eigenstates are obtained by diagonalizing (2.10) Notice that for φ = 0 the (22) components reduce to the expressions in Eq. (2.3). The leading order quantum corrections give rise to a renormalization scale-dependent potential. One can choose the renormalization scale Q to minimize the size of higher order k-loop corrections which scale with (ln m 2 /Q 2 ) k . The scale choice can be important when a parameter in the potential is very different from the electroweak VEV ∼ 246 GeV. We thus choose to use the renormalization group (RG) improved effective potential to minimize the scale dependence. The potential parameters are replaced by their running values, evaluated at the scale Q. The relevant one-loop β functions are given in Appendix A.
The leading order temperature-dependent corrections to the effective potential take the form [24] where the J functions are defined as These functions admit useful high-temperature expansions which allow us to study the phase structure as a function of T analytically (as long as the expansion is justified): The T 2 terms in the expressions above illustrate symmetry restoration at high temperatures. The non-analytic m 3 term in Eq. (2.14) can be responsible for the barrier between the high T phase (at the field origin) and low T phase that breaks Note that symmetry restoration signals the breakdown of perturbation theoryhigher order diagrams become important. This can be accounted for by performing a resummation of daisy diagrams [25][26][27]. The resummation is performed by adding finite-temperature corrections to the boson masses in Eq. (2.12): where c is computed from the infrared limit of the corresponding two-point function.
For the SM Higgs doublet we find The various components of the inert doublet receive similar contributions (but without contributions from the fermions): These expressions are in agreement with those in Refs. [18,28,29]. We implement these corrections by replacing µ 2 i → µ 2 i + c i T 2 in the scalar mass matrices, Eqs. (2.8, 2.9, 2.10). 1 The thermal masses of the gauge bosons are more complicated. Only the longitudinal components receive corrections. The expressions for these in the SM can be found in Ref. [28], but it is easy to modify them to include the contribution of an extra Higgs doublet. For the longitudinally polarized W boson, the result is This includes contributions from gauge boson self-interactions, two Higgs doublets and all three fermion families. The masses of the longitudinal Z and A are determined by diagonalizing the matrix (2.20) The eigenvalues can be written as

Electroweak Phase Transition (EWPT)
Armed with the finite-temperature effective potential, we now proceed to study the structure of the EWPT. The key property we intend to investigate is the transition strength, which sets the baryon number wash-out rate inside a bubble of broken phase (for a recent review of electroweak baryogenesis, see, e.g., Ref. [30]). In order to suppress sphaleron wash-out in the regions of broken electroweak phase, the relevant condition is typically quantified by requiring that [31] v c T c 1, where v c is the Higgs VEV at the critical temperature T c , defined as the temperature at which the origin is degenerate with the electroweak-breaking vacuum. Note that it has been shown that this baryon number preservation condition (BNPC) is a quantity which is manifestly not gauge-invariant [32]. A gauge invariant BNPC can be however derived from the high-T expansion of the dimensionally reduced effective action and the critical temperature T c must be obtained using the gauge invariant prescription of Ref. [32], which employs expansions in powers of of the potential and VEV. However, near the critical temperature, O( ) contributions to the potential are as important as the tree-level terms, so the expansion fails. This is also why an all-orders ring diagram resummation discussed in Sec. 2.2 is needed. For these reasons, below we employ the standard BNPC of Eq. (2.23) and use the full one-loop effective potential to study IDM phases. We will argue that our results do not depend strongly on the issues of gauge invariance.
Finally let us note that the physical phase transition does not begin at T c , but rather at a lower nucleation temperature T n , at which the bubble formation rate exceeds the Hubble expansion rate. If this rate is too slow, the false vacuum is metastable and the transition does not complete. We evaluate the nucleation temperature for a given model with a first-order phase transition using the CosmoTransitions package [33].

Dark Matter
The requirement of a thermal relic abundance for the LOP matching the observed DM density in the Universe of Ω cdm = 0.1199 ± 0.0022 [34], or at least of not over-producing such density via thermal production ("subdominant IDM", see, e.g., Ref. [20]) naturally selects four distinct sectors of the model's parameter space: 4. a heavy mass regime, with a LOP mass between 500 GeV a few TeV.
In the first case, the low mass regime, the DM pair-annihilation predominantly proceeds via the pair production of the heaviest kinematically accessible fermion (τ leptons, b quarks) through h exchange. The lower the LOP mass, the larger the λ L,S couplings need to be in order to produce a large enough pair-annihilation cross section. The allowed mass values range down to values close to the classical Lee-Weinberg lower mass limit for WIMPs [35], for this class of models somewhere in between 3 and 4 GeV. Direct detection limits from XENON10 [36] probe such combinations of masses and couplings quite tightly, such that only a small mass window below 5 − 7 GeV remains. 2 As the mass of the LOP approaches the resonant condition m H ∼ m h /2, the resonant Higgs exchange allows for much smaller values of the λ L,S couplings, and direct detection constraints can be readily evaded. The relevant mass window left unconstrained by XENON100 [38] and LUX [39] has a width of approximately 10 − 15 GeV centered around m h /2 [40].
The mass regions above and below the resonance m H = m h /2 are actually slightly different from each other: Above the resonance, the pair production of W W * in the final state of DM pair annihilation processes becomes increasingly important, even if λ L,S = 0, because the four-point interaction through gauge couplings, independent of λ L,S , starts contributing significantly. As a result the values of λ L,S giving the "correct" relic density are pushed to increasingly (with LOP mass) large, negative values.
For larger and larger LOP masses, the cross section for LOP pair annihilation to gauge bosons becomes very large, such that the thermal relic density is systematically below the universal dark matter density for any combination of model parameters. Barring non-thermal production mechanisms, in this intermediate mass region the LOP cannot be the dominant dark matter constituent [13,20].
Finally, at about m H 500 GeV, for λ L,S 0 cancellations between scalar tand u-channel exchange diagrams and the four-point interaction diagram alluded to above allow, again, for a sufficiently large thermal relic density. Such cancellations are suppressed by driving λ L,S away from zero. Thus, tuning λ L,S for increasing values of m H generally allows one to achieve the correct relic density for mass values from m H 500 GeV up into the multi-TeV range.  Table 1. Input parameters for the three benchmark scenarios discussed in the text along with critical and nucleation temperatures, the transition strength and the signal strength for h → γγ. The masses, given in GeV, are pole masses and the couplings λ i are specified at Q = M Z . Temperatures are also given in GeV.

Benchmark Models
The benchmark models specified in Ref. [13] demonstrated various aspects of DM phenomenology and the possibility for the IDM to influence the h → γγ rate. Unfortunately, none of the suggested scenarios exhibits a strongly first-order EW phase transition.
In this Section, we identify alternate benchmark models which can potentially yield a strongly first-order EW phase transition, while having disparate properties for the lightest Z 2 -odd particle. All our benchmark models are compatible with constraints from Higgs collider bounds and rate measurements, which has been explicitly checked using the tools HiggsBounds [41][42][43] and HiggsSignals [44], where the model predictions have been calculated using a SARAH-generated SPheno version [45][46][47][48]. In the following discussion we mostly focus on the interplay between the dark matter phenomenology and the strengths of the EWPT.
Our key finding is that the requirement of a strongly first-order phase transition generally leads to a large mass splitting between the LOP and the other scalars in the IDM. Our benchmark models are summarized in Tab. 1, along with the corresponding critical and nucleation temperatures, as well as phase transition strengths, as parametrized by the ratio v c /T c . In each case the masses of the A and H ± are chosen to ensure a strongly first-order phase transition. In Fig. 1 we show the dependence of the transition strength on these parameters. The lines corresponding to BM1 and BM3 terminate where the potential develops a non-inert (φ = 0) vacuum first during thermal evolution. This vacuum can then either continuously evolve into the SM/inert (φ = 0) vacuum at T = 0 or it can persist to low temperatures. In the latter case, the EWPT can occur in two steps. Such models are viable if the inert vacuum is deeper than the new one at T = 0 and both transitions complete (i.e., nucleation rate(s) are large enough). Two step electroweak phase transitions have been investigated in detail in Refs. [49,50]. In this work we consider only simple one step transitions, hence the truncation. Notice that in Ref. [18] the strength of the EW phase transition in models with multiple phase transition steps was always weaker, see Fig. 3 and 4 in Ref. [18].
First, let us consider model BM1, where LOP production in the early universe is predominantly set by near-resonant s-channel Higgs exchange. This scenario has been recently examined in the context of phase transitions in Ref. [18]. Even more recently, it has also been suggested as a possible explanation [51] of the Fermi-LAT gammaray excess (see Ref. [52] and references therein). As discussed above, the on-resonance requirement forces m H ∼ m h /2, but allows λ L to be small enough to be consistent with direct detection constraints. Here DM production does not rely on interactions with A or H ± , so their masses can be essentially chosen freely, as long as the resulting quartic couplings λ i (through Eq. (2.3)) satisfy perturbativity and constraints from electroweak precision observables (EWPO), which we check with 2HDMC [53]. In order to satisfy the BNPC of Eq. (2.23) one needs to increase the coupling of the new scalars to h, which, in turn increases the splitting of A and H ± relative to H. For this and the following models we choose m A = m H ± to minimize the impact of splitting these states from H on the Peskin-Takeuchi T function [3] and to reduce the number of parameters. This assumption can be easily relaxed, but the results are qualitatively similar. This benchmark represents the only class of scenarios where the thermal LOP relic density (which we calculated with the micrOMEGAs code [54]) matches the observed DM universal density, and where the EW phase transition is strong-enough first order.
When m A , m H ± 340 GeV, the corresponding loop corrections to m H are large and require µ 2 2 < 0. This causes a second minimum to appear in the potential at T = 0. As m A = m H ± is increased further, this minimum quickly becomes deeper than the SM one, corresponding to the termination of the blue curve in Fig. 1 at m A ∼ 350 GeV. This behaviour was also observed in Ref. [18].
In this scenario, the LOP mass m H has been chosen slightly above the kinematic threshold of the decay h → HH in order to evade constraints from direct searches for invisible Higgs decays and Higgs rate measurements. These however become important for our benchmark scenario BM3 (see below).
The second benchmark BM2 in Tab. 1 represents the intermediate mass regime. Here annihilation into gauge bosons is efficient and DM is generally underabundant, unless there is a cancellation among different amplitudes [13]. The cancellation depends, as indicated above, on how close λ L,S → 0, i.e., on how degenerate the IDM Higgs sector is. In our benchmark, such a cancellation requires m H ≈ m A ≈ m H ± with a maximum splitting of ∼ 10 GeV. These small splittings lead to small couplings of the new states to h and therefore an insufficiently strong phase transition. Thus the phase transition requirement forces thermal relic DM to be underabundant. The observed DM density can be explained here, however, by invoking non-thermal production mechanisms (e.g., the decay of a heavy particle) or with the existence of additional DM particles (e.g., axions). The multitude of "non-standard" production mechanisms has been recently reviewed in Ref. [55].
The final benchmark model, BM3, belongs to the light-mass regime, and is another example that requires further ingredients to be fully consistent with the phenomenology of the DM sector. For m H < m h /2, decays of the SM Higgs to invisible final states become possible, with a decay rate [3] Requiring consistency with the observed 95% C.L. upper limit on the branching fraction, BR(h → HH) ≤ 17% [56], provides a strong constraint on the coupling λ L of  Table 1, which are shown by black dots. The lines for BM1 and BM3 terminate where the inert doublet develops a non-zero vev, φ = 0, as described in the text.
|λ L | 0.007, while a large value |λ L | 0.4 is required to sufficiently deplete the DM abundance [13]. These problems can be remedied by softly breaking the Z 2 , which would allow H to decay [57,58]. As in the previous example, another explanation for DM is then needed. An alternative possibility is to provide the LOP with new annihilation modes, e.g., to new light vector bosons [59], or a mechanism to dilute the thermal relic density, such as an episode of late entropy injection [60,61].
DM phenomenology aside, it is again easy to get a strongly first-order phase transition with a large mass splitting between H and A, H ± . We note that this scenario requires a significant tuning of parameters, because a small LOP mass requires near cancellation of tree-level and loop contributions. For λ L > 0, this leads to negative values of µ 2 2 which can result in the appearance of a new φ = 0 minimum. In all three cases, the first-order transition is driven by the non-analytic (m 2 ) 3/2 terms (see Eq. 2.14) due to A and H ± , while the gauge boson contributions are not as important. This explains the common feature of large splittings between H and A, H ± among the benchmark scenarios. These lead to large couplings between h and the new states, enhancing the size of thermal corrections. This appears to be a generic requirement for increasing the strength of the phase transition in the IDM. It is also important to emphasize that thermal corrections to the crucial (m 2 ) 3/2 terms from A and H ± are not subject to gauge invariance issues that affect the gauge sector contributions.
As a result, we expect these arguments to remain valid in the context of a fully gauge invariant treatment.
The high-T expansion of the effective potential also provides a simple explanation for the shape of the curves in Fig. 1. In this limit the transition strength is proportional to the coefficient of the v 3 term [24]. For the IDM scalars such terms arise from the non-analytic contributions proportional to (µ 2 2 + λ S v 2 c ) 3/2 (assuming m A = m H ± , as above, and ignoring daisy contributions for simplicity), which behaves as v 3 only when λ S v 2 c µ 2 2 . Thus when λ S v 2 c µ 2 2 , the transition strength is independent of IDM parameters, corresponding to the plateau of the green curve in Fig. 1. In the opposite limit, the IDM gives an additional contribution to the cubic coefficient, so the transition strength scales as v c /T c ∼ λ A , as illustrated by the monotonically increasing sections of the curves in Fig. 1. 3 For heavy masses m 2 /T 2 1 (with T ∼ 100 GeV), the IDM states thermally decouple, but this does not mean that they have no impact on the phase transition. When µ 2 2 |µ 1 | 2 , the heavy doublet can be integrated out to yield a SM effective theory with the potential where the dots stand for higher mass dimension operators. The parameters µ 2 , λ and κ can be related to those in the fundamental IDM by equating the effective potentials for the two models at a matching scale Q ∼ µ 2 . For example, one-loop matching yields while µ 2 and λ are determined below the matching scale by fixing the VEV and the Higgs mass. With the presence of a dimension-six term in the potential, the barrier required for a strongly first-order transition can be generated if λ 0 and µ 2 + cT 2 > 0 for T ∼ T c , where c encodes thermal corrections from SM states only [28]. Scenarios of this type have been considered, e.g., in Refs. [62][63][64][65]. One immediate difficulty is that in the IDM κ is generated only at one-loop, so in order for this operator to be significant for field values of around the electroweak VEV, one must overcome the loop suppression, suggesting that the combination λ 3 L + λ 3 S + λ 3 3 /4 cannot be too small. This again forces a large splitting between the IDM states, meaning that the heavy DM scenario described in Section 3 cannot be realized together with a strongly first-order phase transition. Such large couplings can run into perturbativity problems and invalidate the expansion used to generate the effective field theory.
We briefly comment on the discovery prospects of the new IDM states at the LHC. Due to the Z 2 symmetry, the IDM states can only be produced pairwise at colliders. Successive decays of the heavier IDM states A and H ± into the LOP and a Z or W boson, respectively, can give rise to multilepton signatures [10,66,67]. In a recent analysis [67] of LHC searches for supersymmetric particles with two leptons plus missing transverse energy in the final state in the context of the IDM, mass limits of up to m A 140 GeV for LOP masses m H 55 GeV and charged Higgs masses around 85 − 150 GeV have been derived. While these limits partly exceed previous limits from the LEP collider, they are not yet sensitive to the parameter regions that yield a strongly first-order phase transition as required for successful electroweak baryogenesis, see Fig. 1.
The new IDM states can also have an indirect effect on precision Higgs measurements. In particular, the new charged state H ± provides an additional contribution to the loop-induced h → γγ and γZ rates. These effects have been recently studied in Refs. [68,69] in the context of the 125 GeV Higgs boson. Modifications of these branching fractions by O(10%) are a generic feature of our benchmark scenarios, as we show below. The h → γγ rate has the form [21,[68][69][70][71][72] where the leading contributions to the SM amplitude A SM ≈ −6.56 + 0.08i come from W bosons and top quarks. The second term is the new contribution from H ± , where A 0 is a loop function with the property lim x→0 A 0 (x) = 1/3 [72]. For our benchmarks we have λ 3 > 0, so one expects a suppression of h → γγ relative to the SM. 4 Note that for fixed µ 2 2 , the amplitude for the H ± contribution tends to a constant value 1/3 as λ 3 is increased. This means that in the limit of a large mass splitting between H and H ± , which is required for a strongly first order phase transition, the branching fraction is reduced by ∼ 10%. This effect was also noticed in Refs. [17,20] for models similar to our BM1 and BM2, respectively. The deviations to BR(h → γγ) induced by H ± are shown in Tab. 1 in terms of the SM normalized signal strength µ γγ = (σ BR)/(σ BR) SM . While they are still consistent with the present measurements from ATLAS [73] and CMS [74], the LHC should reach a precision of 4-8% for µ γγ [56,75], thereby definitively testing the benchmark scenarios in Tab. 1.
While our benchmarks feature sizable deviations of BR(h → γγ) from the SM expectation, we note that it is possible to avoid this by taking H ± to be nearly degenerate with H, and using A alone to drive the phase transition to be strongly first order. However, in this case, efficient coannihilation of H with H ± during freeze-out generally results in a very small relic abundance [76]. The near degeneracy is also required by constraints on the oblique T parameter when m A m H ± [3]. For example, taking m H ± = 70 GeV, m A = 370 GeV and other parameters as in BM1 results in a strongly first order phase transition, an order of magnitude smaller relic abundance and only a ∼ 3% depression of µ γγ relative to the SM.

Discussion and Conclusions
We studied the structure of the electroweak phase transition in the inert Higgs doublet model, utilizing a set of three benchmark scenarios that feature a potentially viable dark matter particle. Our choices for the three benchmark models essentially exhaust all possible prototypical setups for particle dark matter in the inert doublet model. While only one of the benchmarks has a dark matter particle with a thermal relic density matching the observed dark matter density, the other two (under-and over-abundant) can be made viable by invoking additional production mechanisms or a scenario where the thermal relic density is diluted away, respectively.
The key finding of our study is that in all cases where the model possesses a reasonable particle dark matter candidate, the inert scalar spectrum can be arranged in such a way so as to produce a strongly first-order electroweak phase transition. Central to achieving such a phase transition is to postulate a large enough splitting between the dark matter candidate and the heavier inert scalars. The physics driving this result is simple: Large mass splittings generically correspond to large couplings between the inert scalars and the Standard Model-like Higgs; These, in turn, increase the magnitude of non-analytic ∼ (m 2 ) 3/2 terms in the temperature-dependent effective potential and thus the potential barrier between the field origin and the SU (2) × U (1)-breaking phase.
The mass splitting under consideration cannot be arbitrarily large. For large enough values, for example, the phase structure of the model becomes more complicated, with possible non-zero vacuum expectation values for the inert doublet and multiple-step phase transitions. While, based upon the results of Ref. [18] the latter possibility is not expected to yield stronger electroweak phase transitions than in the single-step case, this is an interesting possibility which we leave for future studies.
The question of how to embed large-enough CP violating sources in detail was also left unanswered here. It will be interesting to study whether such a source (for example an additional gauge-singlet complex scalar, see Ref. [77]) significantly impacts the electroweak phase transition and dark matter phenomenology, and, with this, the conclusions reached in the present study.

A Renormalization Group Equations
Here we list the beta functions for the inert 2HDM at one-loop order. The general form of the RG equations is dλ dt = 1 16π 2 β λ , (A.1) where t = ln Q/Q 0 and Q 0 is a reference scale. We take Q 0 = M Z . The U (1) Y gauge coupling has the GUT normalization: g 1 = 5/3g . The beta functions below have been checked with SARAH [48]. Partial one loop results can be found in Refs. [13,21,78] for dimensionless parameters only. These agree with the formulae below. The gauge coupling evolution is determined by Next we consider the scalar potential parameters. The evolution of the dimensionless quartic couplings λ i is governed by