Electroweak Symmetry Non-Restoration from Dark Matter

Restoration of the electroweak symmetry at temperatures around the Higgs mass is linked to tight phenomenological constraints on many baryogenesis scenarios. A potential remedy can be found in mechanisms of electroweak symmetry non-restoration (SNR), in which symmetry breaking is extended to higher temperatures due to new states with couplings to the Standard Model. Here we show that, in the presence of a second Higgs doublet, SNR can be realized with only a handful of new fermions which can be identified as viable dark matter candidates consistent with all current observational constraints. The competing requirements on this class of models allow for SNR at temperatures up to $\sim$TeV, and imply the presence of sub-TeV new physics with sizable interactions with the Standard Model. As a result this scenario is highly testable with signals in reach of next-generation collider and dark matter direct detection experiments.


Introduction
Given the form of the Higgs potential and the particle contents of the Standard Model (SM), it is a common expectation that at temperatures around 160 GeV the electroweak symmetry should be restored to an unbroken gauge symmetry. However, if the Higgs potential is appropriately modified (or some other source of electroweak symmetry breaking is added) then the electroweak symmetry can remain unbroken to higher temperatures. The general concept of symmetry non-restoration (SNR) was first explored by Weinberg [1] and then expanded upon in [2][3][4][5][6][7][8][9][10][11][12][13][14]. Recently, scenarios of SNR for the electroweak sector, and the associated phenomenological observables, have become a topic of investigation [15][16][17][18][19][20][21][22]. Electroweak SNR is of special interest in the context of certain baryogenesis scenarios [23]. In particular, electroweak SNR can change the minimal temperature at which sphalerons are active if the scale of electroweak restoration is altered [16,17]. For instance, this allows the realization of electroweak baryogenesis at higher temperatures, with a corresponding increase of the new physics scales involved, and a suppression of certain unwanted new physics effects such as contributions to electron EDM. Moreover, if the temperature of the universe following inflationary reheating is lower than the scale of electroweak symmetry restoration, then sphaleron effects can be entirely eliminated thus prohibiting sphalerons from transferring or washing-out of asymmetries generated in the early universe.
One of the defining features of electroweak SNR is the presence of new states interacting with the SM, which are abundant in the early universe and can a priori be stable on cosmological time scales. However, within the previously proposed scenarios, these new states cannot serve as viable dark matter (DM) candidates. We show that this issue can be resolved by extending the scenario with a second Higgs doublet 1 . Besides linking SNR with DM, such an extension allows one to significantly decrease the number of new degrees of freedom required for SNR [17,21]; notably, the original one-Higgs SNR mechanisms required at least O(10) new Dirac fermions, or at least O(100) new singlet scalars.
In this work we will concentrate on SNR induced by new singlet fermions [18]. Compared to the case of singlet scalars [15][16][17], this class of SNR allows for a significant reduction in the overall number of new degrees of freedom, at the price of lower maximal temperatures where it is effective. The number of required new fermions grows with the square of the maximal temperature and starts being a disadvantage compared to the scalar case for temperatures in excess of a few TeV. This limitation however stops being relevant for instance when electroweak baryogenesis is considered in the context of such gauge hierarchy problem motivated scenarios as composite Higgs models [24][25][26][27] or the relaxion mechanism [28]. In the former, the Higgs field "dissolves" at temperatures above ∼TeV, while in the latter the reheat temperatures after inflation are often constrained to be similarly low. Therefore electroweak baryogenesis in both cases has to operate at T 1 TeV even in the presence of SNR. Moreover, it was shown in [19] that in the specific case of composite Twin Higgs models the new states triggering SNR can be naturally present in the spectrum of the model. It might also be highlighted that, generally speaking, embedding of scalar-driven SNR in a naturalness-motivated scenario would at least require a reconsideration due to the presence of new physics needed to explain the lightness of the Higgs and the multitude of SNR scalars.
Notably, here we show that with the inclusion of the second Higgs doublet the number of new states required for SNR can be as little as a single pair of new fermions. We also find that the competing requirements of SNR models with DM tend to force the viable parameter space to be relatively constrained, making these scenarios eminently testable. Moreover, the model constraints imply that a relatively low mass for the second Higgs is a general feature, with a sizable coupling between the two Higgs doublets, leading to an abundance of complementary phenomenological signals at next-generation collider and direct detection experiments. This paper is structured as follows. We start by exposing the difficulties in connecting the SNR states with DM for the case of models with a single (SM) Higgs doublet in Section 2. Subsequent, in Section 3 we introduce SNR in the context of two Higgs doublet models. Identifying the new fermion states which induce SNR with DM, in Section 4 we examine the compatibility of this class of models with the current experimental bounds. Section 5 contains an overview of our results and some concluding remarks.

SNR and DM
We shall start by making some comments regarding the need for a second Higgs doublet in scenarios in which SNR is linked to the DM relic density. Specifically, we shall highlight the issues that arise in models in which the SM Higgs alone is the only significant source of electroweak symmetry breaking.
We will consider scenarios where the electroweak symmetry is broken starting from temperatures of the order 1 TeV and remains broken all the way down to T = 0, while maintaining h/T > 1, which is a critical condition to ensure that the baryon asymmetry is not washed out by electroweak sphalerons [29]. For T 130 GeV a sufficient amount of electroweak symmetry breaking is produced by the SM zero-temperature Higgs potential. However, to ensure that h/T stays greater than unity at temperatures above 130 GeV one needs to compensate the positive Higgs thermal mass induced by the plasma of SM particles which acts to drive the Higgs vacuum expectation value (VEV) to zero (we collect the relevant expressions for the thermal corrections in Appendix A) where λ t is the top Yukawa coupling, λ is the Higgs quartic, g and g are the electroweak couplings, and h is the physical Higgs boson. We will focus on the case in which hightemperature SNR is achieved by adding n χ new fermions 2 with the Higgs-dependent mass term which induces the following thermal correction to the Higgs potential T , at all the relevant temperatures. The lowest such a temperature is T 130 GeV, below which the SM zero-temperature potential is sufficient to ensure h/T > 1. For the estimates in this section we will therefore assume m χ ∼ m χ0 ∼ 150 GeV, although, as we will show later on, the optimal χ mass for SNR is in fact somewhat higher which however would not affect qualitatively the conclusions derived below. It follows from Eq. (2.6) that at least 10 − 40 new fermions χ are required to ensure SNR up to temperatures of T = 0.5 − 1 TeV (which is the typical scale of maximal SNR temperatures that we will be interested in).
Let us now estimate the density of χ states today, where we will assume that the χ are stable (this can be achieved via a Z 2 symmetry χ → −χ) and their relic abundance is set by the standard freeze-out mechanism controlled by the χχ|H| 2 interactions. The χ states dominant annihilation channels are χχ → hh, W + W − , ZZ, the cross-section for which is parametrically and it follows that the χ relic density is given by More precise analytic forms for the cross section and relic abundance are given in Appendix B. Observe from Eq. (2.8) that even with just n χ = 10 new fermions, the correct DM relic density can only be reproduced at the price of an unacceptably low cutoff: Λ 350 GeV. Moreover, perturbativity is limited to the temperatures which are even lower, T Λ/ √ n χ , see Eq. (2.6).
One of the reasons for this outcome is the p-wave suppression of the annihilation cross-section. This can be remedied by adding a pseudo-scalar interaction χiγ 5 χh 2 /Λ which provides an s-wave annihilation channel, permitting the correct χ relic density to be reproduced with higher energy cutoff scales: (2.9) For n χ = 10 the scale of pseudoscalar interactions is thus fixed to a more moderate value of Λ 2 TeV. Therefore, henceforth we will consider models in which χ has both a scalar and psuedoscalar coupling to the Higgs and whileΛ is fixed by the DM relic abundance condition, the scalar interaction scale Λ is unconstrained at this stage. The second major concern with linking χ to the DM are the experimental constraints from DM direct detection. The spin-independent DM-nucleon cross-section induced by the scalar interaction (see Appendix C) is currently bounded as follows [30] σ SI 2.8 × 10 −43 cm 2 (1 TeV/Λ) 2   Such a large suppression of Higgs-χ interactions then has to be compensated by a large number of new fermions in order to have SNR, see Eq. (2.4), n χ > Λ/m χ0 300. This number is in turn incompatible with the relic density requirement of Eq. (2.9) unlessΛ 350 GeV. At the same time the maximal temperature at which the effective field theory (EFT) is valid in the presence of the non-renormalizable pseudo-scalar operator is constrained to be T <Λ/ √ n χ , analogously to Eq. (2.5), which again confines the validity of the theory to unacceptably small temperatures. Moreover, we find that even the special case of resonant DM annihilation (with m h ≈ 2m χ ) is excluded due to the limits on Higgs physics. We conclude that in the minimal models with a single SM Higgs boson presented here, it is not possible to identify the SNR fermions as viable DM candidates. The tension between direct detection and the relic abundance was also found in the models with SNR induced by singlet scalars [16,17], which do not rely on non-renormalizable interactions. The authors of [17] identified some variants in which DM could be included in the model, however the DM candidate was largely decoupled from the SNR mechanism. In this work we specifically restrict our attention to the case that the fields which induce electroweak SNR are at the same time viable DM candidates. In order to allow for that we will separate the SNR sector from the SM quark sector by introducing the second Higgs doublet, as we will discuss now.

SNR with a Second Higgs Doublet
We now consider a modified scenario in which electroweak symmetry breaking at high temperature is due to a second Higgs doublet. In Figure 1 we show an example of the thermal evolution of the Higgs VEVs, which we will analyze in detail in the remainder of this section. The parts of the Lagrangian relevant for electroweak symmetry breaking are (3.1) In the above Lagrangian we assume no interactions of the following types (we omit h.c. coun- , these can be forbidden by H 2 transforming under a combination of Z 2 and global U (1) symmetries. 3 A small coupling between H 1 and χ is induced at one-loop level, which we will account for when analyzing the DM direct detection constraints. All the parameters of Eq. (3.1) are assumed to be real.
The crucial assumptions which allow this scenario to circumvent the issues encountered in Section 2 are that we insist that second Higgs doublet has negligible couplings to SM quarks 4 and modest couplings to the new SNR inducing fermions χ. The first assumption decreases the size of the positive thermal mass, facilitating SNR, and at the same time suppresses the direct detection cross-section once the χ states are identified as DM. Furthermore, in order to comply with the current Higgs couplings measurements and eliminate tree-level direct detection processes we require that at zero temperature the second Higgs doublet has a vanishing VEV, such that all the electroweak symmetry breaking is induced by the first, SMlike Higgs doublet. While postulating suppressed couplings to quarks, we will still assume the presence of small Yukawa couplings between the second Higgs and the SM fermions, allowing the former to decay, and hence not to contribute to the DM density.
Having in mind models of electroweak baryogenesis as the main motivation, we will now analyze how to ensure the condition for all temperatures below some T ∼ 0.5...1 TeV. Provided that the condition Eq. (3.2) is satisfied this will preserve any baryon asymmetry from sphaleron washout.

Vacuum Structure
We will first analyze the global structure of the zero-temperature potential. The physical masses of the Higgs boson and four massive components of H 2 in the desired global minimum at {h 1 , h 2 } = {v SM , 0} are given, at tree level, respectively by Besides taking λ 1 , λ 2 > 0, we will also assume λ 12 > 0, which will be required for SNR as we will see shortly. These assumptions ensure the absence of run-away directions at tree level. For m 2 2 > 0 the potential has only one minimum, which is the correct Standard Model electroweak symmetry breaking vacuum: {h 1 , h 2 } = {v SM , 0}. However for m 2 2 < 0 the scalar potential can feature a second minimum, which can be deeper or shallower. 3 Specifically, H † 1 H2, H † 1 H2|H1| 2 , H † 1 H2|H2| 2 is forbidden by imposing a Z2 symmetry acting on H2. (H † 1 H2) 2 is forbidden by U (1)H2 acting on H2 [21]. H † 1 H2H † 2 H1 is related to the previous one by a custodial symmetry (weakly broken at one-loop level), the term hence vanishes if the custodial symmetry is imposed together with U (1)H2 [21,31]. 4 This coupling structure also avoids dangerous tree-level flavor changing neutral currents, similar to the classic Type I two Higgs doublet model [32].
The case of a deeper minimum would make the Standard Model a valse vacuum, which can be phenomenologically problematic [33]. At tree level, a second deeper minimum develops if m 2 2 < 0, and m 4 2 /λ 2 > m 4 1 /λ 1 , thus to avoid such situations we require that correct global minimum: On the other hand, there is a phenomenologically viable possibility of having a second minimum with h 1 = 0, h 2 = 0 which is shallower than the global minimum at Existence of the shallower second minimum requires, besides satisfying m 2 2 < 0 and Eq. (3.5), also the condition V h 1 > 0 at h 2 2 = |m 2 2 |/λ 2 , which leads to the constraint second local minimum: In summary, condition (3.5) ensures that the global minimum reproduces the Standard Model, but contains no restriction on the presence of possible additional metastable minima. While condition (3.6) corresponds to the existence of a second (global or false) minimum at h 1 = 0 and h 2 = 0. Therefore, these two conditions combined ensure that the second minimum at h 1 = 0 and h 2 = 0 is only a local minimum with the global minimum at {h 1 , h 2 } = {v SM , 0}. Moreover, taking the correct global minimum condition (3.5) along with the converse of condition (3.6) -by flipping the inequality-identifies parameter points with only a single minimum at {h 1 , h 2 } = {v SM , 0} for m 2 2 < 0. These tree-level conditions will be used below to derive analytic estimates for SNR, while in our parameter space scans we will use one-loop zero-temperature potential (see Appendix D) and check the presence of the extra minima numerically.

H 1,2 Thermal Potential
In the early universe the evolution of the Higgs fields is governed by the thermal potential whose approximate form in the high-temperature expansion is where incorporate the fixed contributions proportional to electroweak gauge couplings, top Yukawa coupling and H 1 self-quartic. Here h 1,2 denote the H 1,2 components taking a VEV, whose directions in SU (2) L we assume to be aligned (though in most cases this is unimportant).
Notice the following feature of the contribution of new fermions to Eq. (3.7): while the overall correction to the Higgs mass in the high-T expansion is δm 2 h2 n χ m χ0 T 2 /3Λ and formally grows with m χ0 , this only holds as long as m χ0 T , where the high-T expansion is valid. Otherwise the correction becomes exponentially suppressed, which reflects the depletion of the χ density in plasma. The optimal balance between the linear growth and the suppression is achieved for (m χ0 /T ) 2 ∼ 2, where the thermal correction to the mass reaches its maximal absolute value for a given T . The size of the correction at this point is about a half of what would be predicted from the high-T expansion for given values of m χ0 and T (see Appendix A). Note that we are only using the high-T expansion for the sake of gaining analytic understanding, and a full one-loop thermal potential (see Appendix A) is used in the scans presented in Section 3.4.

Conditions for SNR
Using the analytic approximation Eq. (3.7) we will now discuss the conditions necessary to maintain h Σ /T > 1 from temperatures of at least T ∼ 0.5 TeV down to T = 0, and then move to the numerical results in Section 3.4. At very high temperatures the evolution of the first Higgs doublet is dominated by the SM thermal corrections ∝ c T 1 , which set h 1 = 0 analogously to what happens in the SM. The VEV of the second Higgs at high temperature has the following form where m 2 h2 (T ) and λ 2 (T ) are effective mass and quartic coupling of h 2 which can be read from the second line of Eq. (3.7) and are given by It is useful to write the requirement h Σ = h 2 > T , using the expression for the h 2 VEV in Eq. (3.10), as a bound on the combination of the number of new fermions n χ , their mass m χ0 and coupling strength 1/Λ. Thus we rewrite h 2 > T as follows At very high temperatures, the bound on α SNR is mostly determined by the last terms of Eq. (3.13) which are proportional to T 2 and originate from the thermal correction to the quartic term of the h 2 potential. From inspection of Eq. (3.13) we can derive the maximal T satisfying Eq. (3.13), to obtain (3.14) Another upper bound on the SNR temperature comes from the validity of perturbative expansion in our finite-temperature EFT [18] T In our numerical study we find that both upper bounds on the temperature are important, depending on n χ . At lower temperatures, the right-hand side of Eq. (3.13) is dominated first by the constant terms and then by the terms ∝ 1/T 2 . Notice the presence of the negative contribution ∝ −λ 12 /T 2 , this can be understood from the fact that λ 12 controls the change of the second Higgs mass between h 1 = v SM minimum and h 1 = 0, see Eq. (3.4). For a fixed physical mass of the second Higgs at h 1 = v SM , its mass at h 1 = 0 decreases if λ 12 is positive, which facilitates SNR at low temperatures.
Let us now derive an analytic lower bound on α SNR , and hence on the number of new fermions. To minimize the right-hand side of Eq. (3.13) we will set the second Higgs mass close to the minimal phenomenologically allowed value m h2 = 63 GeV and, at first, also assume negligible quartic coupling λ 2 . The corresponding dependence of the minimal α SNR on λ 12 is shown in Figure 2 as the green line. This bound is obtained by imposing Eq. (3.13) at T = 130 GeV and 500 GeV, the lower-T part is dominated by 1/T 2 terms, while the higher-T part is dominated by constant terms. We omit the ∝ T 2 contribution of Eq. (3.13) whose effect was already estimated in the n χ -independent bound of Eq. (3.14). Also, at T = 130 GeV we assume the bound to be twice what is actually obtained in high-T expansion Eq. (3.13) in order to account for the suppression of the fermionic contribution to the Higgs potential at low T , as was discussed after Eq. (3.7).
A more refined bound on α SNR can be obtained by noticing that at negative m 2 2 the h 2 quartic has to be positive to ensure stability of scalar potential. Imposing a lower bound on λ 2 derived from the correct global minimum condition Eq. (3.5) and from the absence of false vacuum Eq. (3.6), we obtain respectively the orange and blue curves in Figure 2. Gray dots of different shades show the results of numerical scans for n χ = 2, 3, 6 (more details are given in Section 3.4) with a requirement to have no false vacua. Using the analytic results we obtain The fact that the numerically obtained values of α SNR somewhat exceed the analytically derived lower bound is explained by the effect of the terms ∝ T 2 in Eq. (3.13) that we neglected. Furthermore, the bound on α SNR can be recast as a bound on n χ where we have assumed that Λ ≥ 1 TeV and m χ0 ≤ 200 GeV. As we will show in the next subsection, SNR can be achieved for even higher m χ0 , but an increase of m χ0 would need to be accounted for by a further suppression of the fermionic contribution to the Higgs mass at low temperature compared to the high-T expansion result. The estimate in Eq. (3.17) is confirmed by our numerical parameter space scan.
As the temperature drops further, the Higgs fields will eventually transit to the h 1 = 0, h 2 = 0 vacuum. Once this happens, the electroweak symmetry breaking will be supported by the first Higgs doublet which alone provides h/T > 1 for T 130 GeV. In the presence of the shallower minimum at h 1 = 0 the system can remain in false vacuum for some time and then tunnel. In some cases the transition can be smooth and go through a global minimum with both Higgses having non-vanishing VEVs.

Numerical Scan
In this section we present the results of the numerical scans of the model parameter space in Figure 3 and  [1,3] TeV. The scaleΛ is fixed to reproduce the DM relic abundance withΛ 6 TeV, for n χ = 2,Λ 4.5 TeV, for n χ = 3 andΛ 2.7 TeV for n χ = 6 (see Appendix B). Only points with a single minimum of the scalar potential within h 1,2 < Λ,Λ are shown. 5 Figure 3 shows the dependence of the maximal SNR temperature on m h2 and λ 12 , which are the variables most relevant when considering the collider tests of our scenario. Figure 4 5 Models with two minima are numerically more complicated, as one needs to track the metastable vacua.
The tunneling action at each temperature must be computed to ascertain the lifetime of the metastable vacuum to check if it is phenomenologically acceptable. To simplify the analysis we rejected all points with metastable minima, thus our scans are conservative and potentially the viable parameter space could be somewhat larger.   shows the maximal SNR temperature as function of m χ and λ 12 /Λ, which are the parameters relevant for DM direct detection experiments. In Section 4 we will discuss the experimental tests in more detail. As for the remaining parameters not shown in the plots, successful SNR requires Λ 1.5 TeV (with SNR temperatures maximized at maximal Λ) and λ 2 0.1 − 0.2.
We first note that SNR with n χ = 1 requires a rather low scale Λ, around 0.6 TeV, and therefore does not allow to reproduce the DM density. Having performed a separate scan with correspondingly low Λ, we found the maximal SNR temperatures in this case to be ∼ 450 GeV at λ 12 0.18 and m h2 63 GeV. We comment further on this case in Section 5.
For n χ > 6 the maximal SNR temperatures drop (this is the effect of decreasingΛ as we discuss below), while α SNR increases leading to the overall parameter space region with SNR growing. The estimate of minimal α SNR of Eq. (3.16) together with the perturbativity bound Eq. (3.15) predicts the growth of the maximal SNR temperature with n χ as follows however this deviates from what we find numerically. The reason for this is that once we fix the DM density via Eq. (2.9), the scale of the pseudo-scalar interaction becomes tied to the number of new fermionsΛ 6 TeV/ √ n χ . At larger n χ the scaleΛ decreases, leading to a larger one-loop contribution to the Higgs potential, which can eventually produce a barrier within h 2 < Λ, followed by a run-away region. Parameter space points with correspondingly large Λ, which cannot be reliably described within our EFT, are therefore not considered in our scan. Discarding large Λ in turn implies lower SNR temperatures, subject to the perturbativity bound T Λ/ √ n χ . Specifically, we find that the SNR temperatures of at least 0.5 TeV can only be achieved with n χ ≤ 9. While there might be ways to deal with the run-away problem (e.g. arranging for appropriate higher-dimensional operators or using some specific UV-completion), another perturbativity bound T Λ / √ n χ combined with the condition for the relic densityΛ 6 TeV/ √ n χ also eventually limits the number of new fermions to at most 10 for SNR temperatures in excess of 500 GeV. Importantly, the results derived for n χ ≤ 3 are not affected by the mentioned perturbativity bounds or the runaway behavior of the Higgs potential, which allows to consider them as a more robust prediction of our EFT.

Experimental bounds
Having identified the parameter space for which electroweak SNR can be achieved while simultaneously reproducing the DM relic abundance, we next examine the experimental constraints on these models, as well as the future detection prospects. Specifically, we examine the constraints coming from DM direct detection and collider bounds on exotic Higgs decays.
We do not explore indirect detection bounds as these are expected to be model dependent and subleading to those from direct detection over the parameter space of interest. The dominant annihilate channel is χχ → h 2 h 2 and with the h 2 decaying to fermions via their small Yukawa couplings (which could differ from the SM Higgs). This introduces a degree of model dependence in the indirect detection bounds. Assuming, the h 2 have modest couplings to bquarks such that bb-pairs are the dominant annihilation product, we find that the annihilation cross-section is about 3 orders of magnitude below the current sensitivity of Fermi-LAT [34].

DM Direct Detection
The bounds from DM direct detection experiments were fatal to the viability of SNR fermions as DM candidates in the original (single Higgs boson) SNR scenarios. In the model extended with the second Higgs doublet these bounds become significantly relaxed: in the absence of a H 2 VEV today and vanishing H 2 couplings to SM quarks, the direct detection cross-section becomes one-loop suppressed (see Figure 5 for the diagram of the dominant contribution). Moreover, the one-loop interaction between nucleons and χ has to proceed via the crossquartic λ 12 which is also typically less than 1. The leading one loop diagram generates the following effective interactions between quarks and χ particles where the parametric estimate (see Appendix C) of the effective couplings a q reads While the overall quark-χ coupling is also expected to receive UV contributions from the cutoff physics, Eq. (4.2) can be used as a lower bound (assuming no cancellations happen between different contributions). The resulting DM-nucleon cross section is [35,36] where µ χN is the reduced mass of DM and a nucleon, Z and A are the atomic and mass numbers of the target nuclei, and we have neglected the subleading contribution due to the pseudo-scalar interactions. The effective DM-nucleon couplings g N are defined as with the matrix elements of quark currents in a nucleon N = {p, n} given by [37,38] .

(4.5)
The values of f (n) T q can be computed in chiral perturbation theory using measurements of pion-nucleon sigma term [39][40][41][42], where we use a representative set of values adopted in the DarkSUSY package [43] f p T u = 0.023, f p T d = 0.034, f p T s = 0.14, f n T u = 0.019, f n T d = 0.041, f n T s = 0.14.
Darks ide-2 0k ⋁ floor Figure 6. Current (XENON1T [30], LUX [44], PandaX-II [45], PandaX-4T-Commissioning Run [46]) and projected (XENONnT [47], DARWIN [48], LZ [49], PandaX-4T [50], Darkside-20k [51]) experimental bounds on spin-independent nucleon-DM cross-sections as a function of DM mass. The green regions indicate the cross-sections corresponding to SNR with T ≥ 0.5 TeV for n χ = 2, 3, 6 in the parameter range satisfying the Higgs physics bounds, for the χ-quark couplings given in Eq. (4. 2) The gray areas are obtained by increasing or decreasing the χ-quark couplings by a factor of 3, for n χ = 2. The neutrino floor is also displayed. Figure 6 shows that the DM-nucleon cross sections favored by SNR lay in a region between the excluded area and the reach of currently constructed experiments. Note that even a factor of 3 enhancement (or suppression) of a q , which in fact are free parameters in our EFT, can bring a large fraction of the SNR-favored parameter space into the excluded region (or outside the region that will be probed in near future experiments). Notably, in case of such a coupling suppression most of the green region will remain above the neutrino floor.

Collider Bounds
If m h1 > 2 m h2 the Higgs boson can decay into a pair of H 2 components, with the total decay width given by which, depending on the H 2 decay mechanism can be classified as invisible or untagged h 1 decays [52]. Such exotic Higgs decays into H 2 are bounded by BR BSM < 0.2 [53]. Neglecting the phase space suppression factor, this results in a bound λ 12 0.007.
In the opposite case m h1 < 2 m h2 the strongest bound on λ 12 was derived from a modification of the 1-loop h 1 -photon coupling [21] (assuming m h2 = m h1 /2) and reads λ 12 0.4. This bound will potential be improved by about one order of magnitude at HL-LHC [54]. Both discussed current bounds were imposed in the plots presented in Figures 3 and 4.

Concluding Remarks
This paper has highlighted that DM can extend electroweak symmetry breaking to higher temperatures. The main ingredient allowing for the link between DM and SNR is the addition of a second Higgs doublet, which both facilitates the symmetry breaking and suppresses direct detection signals. Unlike scenarios in which SNR is not linked to DM, our model has an upper bound on the number of SNR states: n χ 10. High-temperature symmetry breaking at this low number of new states relies on the non-vanishing positive cross-quartic between the two Higgs doublets and requires the second Higgs to be relatively light. These facts imply that the model will be tested at the forthcoming collider and, potentially, direct detection experiments. The latter can, however, be evaded with a mild tuning of the DM-quark interactions (with the exception of DARWIN). While a positive identification at a direct detection experiment would not be sufficient to confirm this model, however it is notable that the model has implications for collider signatures which would allow for correlated signals. In particular, one might be optimistic regarding probing the second Higgs state H 2 through precision Higgs studies.
The model studied here features an explicit cutoff Λ, which calls for a UV completion. The desired χχ|H 2 | 2 interactions can be produced in multiple ways, for instance, by integrating out a heavy scalar with a Yukawa coupling to χ fermions, integrating out new SU (2) L doublet fermions, or can be produced by some new strong dynamics at the scale Λ [18]. Since the cutoff is necessarily TeV-scale this implies that the UV physics can be potentially probed in the near future collider experiments, and may be linked to the hierarchy problem (as in [19]).
Notably, the fact that SNR is achieved with a small number of additional states is very elegant (as opposed to O(100) in the original bosonic models [15][16][17]). Indeed, it may be tempting to draw a connection between the three generations of the Standard Model and the n χ = 3 SNR model, which is seen to work well in Figure 6. The first detailed analysis of the effect of the second Higgs doublet on SNR was recently performed in [21] for the case of SNR with new singlet scalars, this work also predicting a lower number of SNR scalars (order-20) compared to the original single Higgs models [15][16][17]. Moreover, these authors presented a detailed phenomenological analysis regarding detection prospects of the second Higgs doublet, which is highly relevant to the models outlined in this work.
Let us also comment further here on whether the case of n χ = 1 can be accommodated within our scenario. Setting n χ = 1 implies that the SNR condition Eq. (3.16) leads to an upper bound: Λ 0.7 TeV (close to the value obtained in our numerical scan). While such values of the cutoff appear dangerously close to the temperatures at which we would like SNR to operate, we now also encounter a problem of DM underproduction. With such low n χ and Λ, the scalar interaction alone leads to the χ relic density Ω χ 0.05. Thus, while n χ = 1 is interesting for SNR (provided that the low cutoff does not interfere with SNR, which can be checked in a specific UV completion), it cannot explain the observed value Ω DM 0.12.
Small changes to our model give alternative scenarios which may also be of interest. For instance, new SNR fermions might be involved in the generation of neutrino masses, however, in this case it would seem unlikely that they could also play the role of DM. Alternatively, the neutral component of the second Higgs doublet could also be stable, while the sizable cross-quartic required by SNR with low n χ implies that this state cannot account for all of the DM without tension with direct detection limits [55], giving a motivated route to a two-component DM model. We leave these variations for future studies.
As alluded to in the introduction, the principle motivation for delaying electroweak symmetry restoration until high temperatures is the realisation of electroweak baryogenesis at higher temperatures. Since the analysis of the strength of the phase transition and baryon number violation are distinct sets of calculations, which also require additional model building, we have chosen to study the connections between electroweak symmetry non-restoration and dark matter in isolation from these other issues. Implications of such high temperature electroweak baryogenesis have been explored in [16][17][18]. However, the scenario that we have studied here is distinct from these previous models since in order to reduce the number of χ states in our theory and satisfy direct detection bounds, a second Higgs field H 2 was introduced which significantly alters the vacuum structure of the theory.
Notably, at temperatures above ∼100 GeV it is H 2 which is the primary source of electroweak symmetry breaking. Thus the universe essentially undergoes a two-step phase transition evolving from the unbroken phase, through a phase in which the H 2 VEV breaks the electroweak symmetry, and then finally transitioning to the observed low energy vacuum state with the Standard Model Higgs being the main source of electroweak symmetry breaking. Successful mechanisms of baryogenesis via (low-temperature) two-step electroweak phase transitions have been previously explored in the literature [65][66][67], and the EW phase transition can be made first order with appropriate parameter choices. Building on the results of this work, we plan to report on the prospects for realisation of electroweak baryogenesis in the context of this model in a future publication.

A Thermal Corrections
At finite temperature T the scalar potential receives the following corrections respectively for one thermalized bosonic degree of freedom and one Dirac fermion with mass m. Their interactions with the Higgs field are encoded in the Higgs-dependent masses m.
The thermal loop functions are defined as In the high-temperature limit m 2 /T 2 1 they simplify to Let us consider the correction induced by the new fermions with mass m = m χ0 − h 2 /Λ. Applying the high-T expansion, we find the Higgs mass correction grows with m χ0 while the true quadratic correction is given by and its absolute size is maximized at m 0 χ /T 1.4 where it reaches ∼ 1/2 of the naively expected value in Eq. (A.4) The overall thermal correction to the Higgs potential is given by

B Relic Abundance
In this appendix we present the annihilation cross-sections and DM relic density for n χ Dirac DM fermions χ which annihilate into n φ scalars φ, with the Lagrangian The annihilation cross-section is given by where the physical χ mass is The annihilation cross-section can be decomposed into s and p wave components as where v is the relative velocity between two χ particles in the centre of mass frame and It follows that the relic abundance of χ is given by [57] Ωh 2 1.07 × 10 9 where M Pl 1.22 × 10 19 GeV is the Planck mass. In our numerical calculations we take g * 86 and x F 25.

C Direct Detection Cross-Sections
In this appendix we provide details on the calculation of the direct detection cross-sections used in Section 4.

C.1 Tree level
We parametrize the effective χ-quark interaction obtained after integrating out heavy mediators in the following way and the corresponding scattering cross-section is given in Eq. (4.3). The following tree-level contribution to a q would arise in the SNR scenarios with a single Higgs doublet (C.2) Figure 7. DM-quark interaction at one-loop level (produced using TikZ-Feynman [58]).

C.2 One loop
The leading contribution to χ-quark interactions in the discussed two Higgs doublet model with SNR arises at one loop level. The corresponding Feynman diagrams are depicted in Figure 7. The h 1 -mediated scattering (1A) leads to the following correction to the effective four-fermion interaction of Eq. (C.4) where we used dimensional regularization with the M S scheme. In the initial EFT containing h 1,2 the divergence appearing in the computation of a q would be absorbed in the renormalization of the operator h 2 1 χχ, which contributes to the quark-χ interactions as where c h 1 χ is a free parameter of the EFT. Not aiming at matching our theory to some specific UV completion, we simply estimate the resulting χ − q interactions by substituting µ = Λ in Eq. (C.3) and omitting δa (contact) q . Note that the one-loop diagram involving gauge bosons (1B) instead vanishes, as the corresponding amplitude requires to couple (pseudo)scalar and vector currents [59].

C.3 Two loops
The parametric estimate of the size of the two loop contributions (see Figure 8) to the fourfermion interaction is which represents a small correction to the 1-loop contribution with a relative size g 4 (4π) 2 /λ 12 10 −3 /λ 12 . (C.6) These short-range corrections do not include the 2-loop contribution with two photon lines which has to be treated separately. Given that the latter has a different energy dependence, we will perform the comparison at the level of differential event rates. We will calculate the photon-mediated two-loop scattering in two steps by first integrating out the second Higgs doublet and then using the resulting effective operator to compute the cross-section.
In the first step we obtain (see Figure 9) which results in the following DM-nucleus scattering amplitude [59]   where E R = − q 2 2m i is the recoil energy of a nucleus, q is the four-momentum transfer, m i is the target nucleus mass, Z is its atomic number andF (q) is the corresponding two-loop nuclear form factor. Correspondingly, the DM-target nucleus i differential cross section reads where v χ is the DM velocity in the lab frame and µ χi is the reduced mass of χ and a nucleus.
The minimum value for v χ to induce certain nuclear recoil energy E R is v min χ = m i E R /2µ 2 χi , which explains the presence of the step function.

C.4 Event rates
The differential event rate per unit target mass and per unit time is given by [60,61] where N i is number of target nuclei per unit target mass, f ( v χ ) is the normalized DM velocity distribution in the lab frame, and ρ χ = 0.3 GeV/cm 3 is the local DM density. The differential cross-section dσ χi /dE R for the photon-mediated scattering was given in Eq. (C.9) whereas for the one-loop Eq. (C.3) and the remaining two-loop Eq. (C.5) contribu-tions it reads where σ χN was given in Eq. (4.3) as a function of a q and A is the atomic mass number.
To compare differential rates we take Helm form factor [62] for both F (|q|) andF (|q|) (in the latter case the exact form factor is not expected to be important as this contribution is very suppressed), and we use the conventional truncated Maxwellian distribution [60] for the local DM velocity distribution. The comparison of differential rates for the 1-loop contribution, 2-loop contribution without photons, and the 2-loop contribution with photons is presented in Figure 10 for Xenon targets with A = 131. For a typical SNR-favored parameter choice with m χ = 0.5 TeV, m h2 = 70 GeV, λ 12 = 0.2, and Λ = 1 TeV the 1-loop contribution completely dominates over the others, which can therefore be safely neglected in our analysis.

D The h 1,2 Effective Potential
The tree-level zero-temperature potential is fixed by Eq. (3.1). The one loop zero-temperature correction to the h 1 and h 2 potential in the Landau gauge is given by with F = 0 for bosons and F = 1 for fermions, g {h 1 ,G 1 ,h 2 ,G 2 ,W,Z,t,χ} = {1, 3, 1, 3, 6, 3, 12, 4n χ } and m i are corresponding h 1,2 -dependent tree-level masses. We split H 1,2 into h 1,2 (which obtain a VEV at some T ) and three G 1,2 (with vanishing VEVs). When presenting the phenomenological bounds in the main text we assumed that in the current minimum all H 2 components have the same mass and same coupling to h 1 , thereby neglecting a small splitting induced by gauge bosons loops [63]. We fix the counterterms ∆X by the following conditions where λ 12 is defined from the h 1 h 2 2 coupling at q 2 = −(0.5 GeV) 2 −(10 −3 m χ ) 2 relevant for direct detection experiments, and m 2 h1 and m 2 h2 are pole masses. In ∆Π h 1 and ∆Γ h 1 h 2 2 we account for the leading contributions to the difference between the zero momentum transfer quantities (derivatives of the effective potential) and the quantities at relevant finite momentum [56]. The top quark and IR-divergent H 1 contributions to the h 1 self-energy difference reads [64]  Note that m 2 G vanishes at the minimum of the effective potential, but the corresponding IRdivergent pieces from ∆Π h 1 , ∆Γ h 1 h 2 2 and derivatives of the effective action cancel out from the renormalization conditions Eq. (D.2). We do not include ∆Π h 2 in Eq. (D.2) since it does not receive contributions from the top quark or with IR-divergences.
We do not impose any renormalization condition on λ 2 since in our study it is a mute parameter which is scanned over, and we set ∆λ 2 = 0. The non-renormalizable interaction h 2 2 χ 2 induces UV-divergent corrections to the operators ∝ h 6 2 , h 8 2 . We set corresponding counterterms ∆c 6,8 to zero as well. As a result of these prescriptions, the one-loop effective action is explicitly µ-dependent after renormalization, and we fix µ = 1 TeV.