A New Approach to Electroweak Symmetry Non-Restoration

Electroweak symmetry non-restoration up to high temperatures well above the electroweak scale offers new alternatives for baryogenesis. We propose a new approach for electroweak symmetry non-restoration via an inert Higgs sector that couples to the Standard Model Higgs as well as an extended scalar singlet sector. We implement renormalization group improvements and thermal resummation, necessary to evaluate the effective potential spanning over a broad range of energy scales and temperatures. We present examples of benchmark scenarios that allow for electroweak symmetry non-restoration all the way up to hundreds of TeV temperatures, and also feature suppressed sphaleron washout factors down to the electroweak scale. Our method for transmitting the Standard Model broken electroweak symmetry to an inert Higgs sector has several intriguing implications for (electroweak) baryogenesis, early universe thermal histories, and can be scrutinized through Higgs physics phenomenology and electroweak precision measurements at the HL-LHC.


I. INTRODUCTION
The Standard Model (SM) of particle physics accurately describes the behavior of the particles making up the ordinary matter, but it fails to provide an explanation of how they came to be. Under the assumption that particles and anti-particles are produced in equal numbers in the early Universe, the SM predicts that they would have long annihilated each other without leaving any remnant matter today. Sakharov [1] enunciated that producing a Baryon Asymmetry (BA), i.e., more matter than anti-matter, requires baryon number violation, C and CP violation, and out-of-equilibrium processes to all occur at the same time. Although the SM provides sources of C, CP, and baryon number violation through the electroweak interactions and sphalerons, respectively, it fails to explain the observed BA. Indeed, the SM Electroweak Phase Transition (EWPT) is a smooth crossover and, thus, is not giving rise to sufficient deviations from thermal equilibrium [2]. In addition, the amount of C and CP violation in the SM is insufficient to generate the observed baryon asymmetry [3]. In order to generate the observed baryon asymmetry, sources of CP violation and out-of-equilibrium processes beyond those found in the SM must be realized in nature.
There are many mechanisms proposed in the literature to explain the generation of a net Baryon number B, and in most cases, sphaleron processes that are capable of violating B+L, but conserve B-L, play a relevant role (with L the lepton number). One interesting possibility to achieve sphaleron-induced B number generation is via a Strong First Order Electroweak Phase Transition (SFOEWPT), yielding promising conditions for electroweak baryogenesis [4]. Accommodating a SFOEWPT demands modifications of the Higgs potential. Such modifications may be induced predominantly by thermal effects, as it happens e.g., in the Minimal Supersymmetric extension of the Standard Model (MSSM) [5][6][7][8][9][10][11], or by zero-temperature effects that have a lasting consequence after thermal effects are taken into account. The latter situation naturally occurs in models of new physics containing additional light scalar particles with sizable couplings to the Higgs.
In this study, we are interested in models in which the electroweak (EW) symmetry is broken at temperatures well above the EW scale. Taking a bottom-up approach, we called these scenarios: i) delayed restoration, if the electroweak symmetry is restored at very high temperatures, or ii) non-restoration if the electroweak symmetry remains broken all the way up to some high energy scale Λ of validity of the theory. Electroweak non-restoration or delayed restoration scenarios have advantages in modeling mechanisms for baryogenesis. For example, in the case of electroweak baryogenesis (EWBG), one important advantage is that the additional, required sources of CP violation will only be effective at high energies and, therefore, will avoid current electric dipole moment experimental bounds.
Symmetry non-restoration at high temperatures has been first studied a long time ago [12][13][14][15][16][17][18] and recently [19][20][21][22][23][24][25][26][27]. In particular, new ideas of electroweak symmetry non-restoration or delayed restoration have been discussed [22][23][24]26] by extending the SM Higgs sector with additional singlet scalars that couple to the SM Higgs and provide it with a negative thermal mass at very high temperatures. Such models typically require several hundreds of new scalar fields. On top of the new scalar sectors, for models with delayed restoration, the Ultraviolet (UV) completions typically require additional scalar and/or fermion fields that couple with the EW sector and yield electroweak symmetry restoration, as well as a strong first-order phase transition, at very high temperatures [23,24,26].
In this work, we explore the EW non-restoration or delayed restoration with an extended-Inert Two Higgs Doublet Model (I2HDM) [28,29], where, instead of the SM Higgs, it is the inert Higgs who acquires a non-zero vacuum expectation value (vev) up to very high temperature by coupling to an additional scalar sector. Such a model requires minimal couplings between the new scalar fields and the SM Higgs boson, and opens the window to different realizations for baryogenesis at very high energy scales. Due to the lack of large Yukawa couplings to the inert sector, the number of scalars required to achieve negative thermal masses is somewhat reduced. Specific new physics models for high scale EW baryogenesis in the context of extensions of the I2HDM will be the topic of a forthcoming publication.
The EW non-restoration sets the boundary conditions at high temperatures ( T UV ), while the observed EW vacuum defines them at zero temperature; see the top row of Figure 1 for a schematic view. For intermediate temperatures, the I2HDM allows different phase histories that we depict in the bottom row of Figure 1. There could either be a temperature range (between T c H and T c Φ -to be precisely defined below) for which the global minumum is given by non-vanishing vevs of both the Higgs and the inert fields (left plot), or there could be a discrete jump between the Higgs and the inert phases at a critical Temperature T c (central plot). A third option is given by a scenario in which the Higgs vev goes to zero at a restoration temperature T r H lower than the temperature T r Φ above which the inert vev starts to grow (right plot). In the temperature range between T r H and T r Φ the system is in a EW preserving vacuum.
In this work, we utilize the perturbative effective potential (EP) method to calculate the finite-temperature phase structure and quantities relevant to the baryon asymmetry. However, unlike for typical EWPT calculations for which the electroweak symmetry breaking takes place close to the EW scale, here we need to take into account important effects due to the large scale separation between the high temperatures ( O(1 − 100) TeV)high field values and the EW scale, which requires careful treatment and improvement of the perturbative calculation. For this purpose, we will implement a Renormalization Group (RG) improvement and daisy resummation of the EP to ameliorate the perturbative convergence.
This paper is organized as follows: in section II, we introduce our model and discuss its zero temperature constraints. In section III, we investigate the validity of the radiatively corrected, finite temperature effective potential in calculating the phase structures, introducing the RG improvement and daisy resummation, and we set up schemes for the improved perturbative calculation. In section IV, we present an analytical study of the possible thermal histories based on a mean-field approach. In section V, we present the full numerical computation of the finite-temperature phase structure for two benchmark (BM) scenarios. In section VI, we discuss the baryon washout conditions and consider them in light of the thermal history results for the two BM scenarios presented in the previous section. We also discuss the impact of future model building on high-temperature baryogenesis. In section VII, we discuss phenomenological constraints in this type of model. Finally, we present our conclusion in section VIII. We collect various technical aspects in appendices.

II. THE MODEL
A. The effective potential at tree level We consider an extension of the SM Higgs sector that includes an Inert Higgs Doublet with additional singlet scalars. In such case, the most general Z 2 -symmetric potential reads 1 , where the two Higgs doublets are written as and the fields χ i represent N real, singlet scalars. Assuming that extra sources of CP violation will come from a new sector, once we study the complete UV theory, we impose CP invariance in the Higgs sector and define all model parameters to be real. The assumed Z 2 -symmetry forbids couplings of the type ert, which forbids additional terms that we omitted here in the potential. µ 2 12 (H † Φ), λ 6 (H † ΦH † H), and λ 7 (H † ΦΦ † Φ). Portal couplings of the form (Φ † H)(Φ † H) and (H † Φ)(H † Φ) are allowed by the Z 2 symmetry and are related to the operator (H † Φ)(Φ † H) by custodial symmetry [30]. However, assuming a U (1)-symmetry on (one of the) doublets forbids these additional portal couplings and simplifies the potential. Given the custodial symmetry and the additional U (1)-symmetry, we can set the coupling λ HΦ to 0 as well. However, this is not stable under RG-running, as the hypercharge gauge coupling breaks custodial symmetry. We therefore keep track of the operator with the coefficient λ HΦ for future RG improvement of the EP; see discussion below in section III. In addition, to better accommodate phenomenological constraints, we set λ Hχ = 0, although, similarly to λ HΦ this coupling will also be induced by the renormalization group evolution (RGE), and we will keep track of its effects. Finally, observe that λ χ = 0 is protected by an SO(N ) symmetry of the singlet sector, and we shall impose such symmetry. In the case of a potential with generic values of λ χ , the singlet sector exhibits a discrete Z N symmetry.
To summarize, parameters in the above potential can be separated as follows: • free parameters set to zero: { λ HΦ , λ Hχ , λ χ }, where the two fixed parameters are given by the current observation of the EW vacuum expectation value (vev) v 0 = 246 GeV and the SM Higgs mass m h = 125 GeV.
In general, there could be charge breaking and CP breaking minima in two Higgs doublet models. However, [31,32] showed that at tree level, if an EW breaking minimum exists, any possibly existing charge breaking or CP breaking extremum is necessarily a saddle point above the EW breaking minimum. Although the validity of this result may not hold after the inclusion of radiative corrections, and its validation requires a more detailed analysis beyond the scope of this work, we shall only allow for the neutral CP even components to develop non-zero vacuum expectation values at any temperature. Therefore, from now on, we focus on analyzing the effective potential of the CP-even components of the two Higgs doublets and the singlet sector. The tree-level CP even potential reads, The particles in the plasma include bosons {h, G 0 , G ± , ϕ, φ 0 , φ ± , χ, γ, W ± , Z} with corresponding particle degrees of freedom (d.o.f.) n bos = {1, 1, 2, 1, 1, 2, N, 3, 6, 3}, and fermions, {t} with corresponding particle d.o.f. n f erm = {12} that couple (self-couple) to the dynamical fields. Notice that we work in the Landau gauge so there are no ghost d.o.f. We collect the effective, field-dependent masses of these particles in appendix section A.

B. Zero temperature constraints
In this section, we present the tree-level, zero temperature constraints on our model, including the bounded from below (BFB) conditions, and the correct vacuum structure of the tree-level potential. This study provides guidance, later on, in defining the viable parameter space for which we shall perform numerical calculations to constrain the model after the inclusion of radiative corrections.

Bounded From Below Conditions
The bounded from below (BFB) conditions, which need to be satisfied simultaneously, for the generic tree level potential given in eq. (1) are where for simplicity we define the effective couplings Λ χ,n ≡ 1 n λ χ + λ χ and Λ HΦ ≡ λ HΦ + λ HΦ ρ 2 . (6) There are two variables in these conditions, n ∈ {1, . . . , N } and ρ 2 ∈ [0, 1], see appendix section B for details. The conditions (5) have to hold for all values of n and ρ. Notice that they only enter the conditions through Λ χ,n and Λ HΦ . If λ χ > 0, Λ χ,n is the smallest when n = N , while if λ χ < 0, the smallest Λ χ,n is found for n = 1. Similar considerations apply to Λ HΦ and ρ. A detailed derivation of these conditions can be found in appendix section B.

Vacuum Structure
In order to be consistent with the current Higgs and EW precision measurements, as the inert doublet is charged under the EW gauge group, we consider the case that at zero temperature, both the inert Higgs and the singlets have zero vev, say the physical vacuum is where v 0 = 246 GeV, and we require such vacuum state to be the global minimum of the zero temperature potential. Firstly, for the physical vacuum to be a minimum, one needs to avoid tachyonic solutions, which give constraints on the bare mass parameters of the potential (at tree level) Equation (8) does not involve the RG-generated parameters, λ HΦ and λ Hχ since it refers to the couplings at the physical minimum. As stated above, at tree level, any possibly existing CP or charge breaking extrema are saddle points lying above the EW vacuum, which, therefore, do not put any further constraints on the viable parameter space. To secure that the EW vacuum is the global minimum of the tree-level potential in the subfield space of the two CP even components and the singlet degrees of freedom, we find all possible extrema of the polynomial potential (see all possible extrema in appendix section B at tree level) and we numerically impose the necessary conditions to establish that for each extremum either it cannot exist, or it is above the physical one.

III. RADIATIVE CORRECTED, FINITE TEMPERATURE POTENTIAL
In the perturbative effective potential calculation, two types of radiative corrections to the tree-level potential need to be considered, i.e. the zero temperature loop corrections and the finite temperature radiative corrections.
At one-loop order, the zero temperature loop correction can be taken into account through the Coleman-Weinberg (CW) potential [33,34] under the MS-renormalization scheme, and where S i = 1 or 1/2 for i =B or F, respectively. The short-handed notation has been introduced for the dynamical fieldŝ Φ ≡ {h, ϕ, χ 1 , χ 2 , · · · , χ N }. The specie i is summed over all degrees of freedoms in the plasma. The constant a i has a value of 3 2 for scalars, longitudinal gauge bosons and fermions while 1 2 for transverse gauge bosons. µ R is the renormalization scale, and finally M 2 i (Φ) is the fielddependent mass eigenvalue. The field dependent masses of all degrees of freedom in the plasma for our model are given in appendix section A. We work in the Landau gauge [34], which introduces a gauge-dependence of the EP [35][36][37][38][39][40][41][42][43] The CW potential changes the shape of the zero temperature potential, introducing deviations from the tree level constraints at zero temperature that we discussed in the last section. Specifically, to accommodate the Higgs vev of 246 GeV and a 125 GeV mass eigenstate the parameters µ 2 H and λ H have to be adjusted to recover the two physical conditions at T = 0. Other zero temperature constraints, including the BFB and correct vacuum structure, also need to be adjusted numerically, as necessary, so that they remain robust after the inclusion of loop corrections.
The leading temperature dependence is given by the thermal one-loop effective potential (see reviews e.g. [44]) 2 Given that we observe that the high-temperature expansion approximation is in good qualitative agreement with the full treatment of the temperature effects when considering the electroweak symmetry non-restoration analysis, we argue that the main results of this work will not be qualitatively changed by effects of gauge dependence. Indeed, the EW non-restoration at high temperatures relies on a negative thermal mass for the inert Higgs that is governed by the leading order term in the hightemperature expansion, which in turn does not exhibit gauge dependence. Indeed Refs. [40,41] show that the gauge dependence appears only in the sub-leading temperature-dependent terms in the high-temperature expansion. A dedicated study of the gauge dependence considering a numerical analysis of the full temperature-dependent EP would be necessary to fully understand the relevance of gauge-dependent effects in the analysis of EW non-restoration, which is beyond the scope of this work.
where relevant notation has been introduced above, and where in the thermal potential the argument y = M 2 i (Φ)/T 2 . The J B/F functions can be evaluated numerically, see appendix section F for the details of our implementation. To gain analytical understanding of the thermal history, a high-temperature expansion can be used to obtain an analytical expression for the thermal potential: where and γ E is the Euler constant. The high-temperature expansion in eq. (12) guarantees a good convergence for values of the argument of the J B/F functions up to 2 − 5, while values are constrained to be below/about 1 without inclusion of the logarithmic terms. At very high temperatures and very large field values, which are the relevant scales for the electroweak symmetry non-restoration or delayed restoration scenarios, perturbative convergence of the fixed order calculation becomes compromised, for both the CW and the oneloop thermal potential. In the following, we discuss the improvements that we will implement to deal with both shortcomings.
As it is well understood in the literature, at finite temperature, the self energy of a particle receives higher loop corrections from daisy diagrams, e.g. see [45,Fig. 3a]. Such corrections at N -loop order contain powers of a field-and temperature-dependent parameter α (up to a normalization factor) [12,[45][46][47][48], where λ i is the coupling corresponding to M i (Φ) in the theory. At large temperatures, such contributions exhibit severe IR divergence for some field values such that M i (Φ) T , for example around the origin, where higher loop contributions dominate and the fixed order calculation becomes problematic. Various treatments have been proposed to resum higher loop thermal contributions and solve the associated IR problem [12,[45][46][47][48][49][50][51][52]. A full dressing daisy resummation involves adding thermal corrections to the tree level effective masses in the effective potential. For the one loop EP it follows, where Π 2 i is the squared thermal mass for the specie "i". Such a procedure effectively resums higher order corrections from daisy diagrams 3 .
The squared thermal masses Π 2 i are in general fieldand temperature-dependent and can be solved by gap equations. At one-loop level the gap equations read where the degree of freedom i appears as a background field in the EP. A truncated treatment involves doing an expansion of the right hand side of the gap equation with respect to Π 2 k and truncate to a given order. To the leading order, the truncated squared thermal mass reads If the thermal potential is evaluated to leading order in high-temperature expansion, one obtains the well known field-independent form of the squared thermal masses The c i are constant coefficients dependent on couplings determined by the theory, and we collect the thermal mass coefficients c i for all degrees of freedom in our model in appendix section C. We implement the high-T thermal masses in eq. (17), the truncated thermal masses in eq. (16), and the gap thermal masses in eq. (15) in comparison, to effectively resum higher-order daisy diagrams. In Figure 2, we show the squared thermal mass of the scalars as a function of the inert field values at a temperature of T = 5000 GeV, computed with the different levels of accuracy described above, for the BM scenario B to be defined in Table I. The high-T thermal masses should be independent of the inert field value φ, however, Figure 2 shows a small variation with respect to the field value due to the RG improvement implementation to be discussed below. The truncated thermal masses have an enhanced dependence of the inert field value, especially for the inert thermal mass itself, but a more sizable variation occurs for the gap thermal masses. The differences among thermal masses for different implementations as shown in Figure 2 will end up, however, having a very small impact on the results relevant for the phase structure of the EW non-restoration BMs. The fixed order EP at finite temperature, including both the zero temperature and thermal contributions, depends on the scale µ R at which the theory is renormalized. For example, at one loop order, using the hightemperature expansion in eq. (12), the potential has a logarithmic dependence on the renormalization scale as where the log(M 2 i (Φ)+Π 2 i ) piece is cancelled between the CW and logarithmic term in the high-temperature expansion of the thermal potential contribution. By implementing RG improvement, where the parameters, fields and vacuum energy of the potential are evaluated at the scale µ R , one would cancel the scale dependence to the order of the calculation. As we only calculate the effective potential and the RG improvement at one-loop order, the scale µ R needs to be chosen wisely to avoid un-resummed large logarithms from higher-order loop effects. Formally, at T = 0 with a convenient choice of the renormalization scale, the L-loop effective potential with an RG improvement at (L + 1)-loop order, is exact up to L-th-to-leading log order [54][55][56]. At finite temperature, the choice of the renormalization scale, should vanish or minimize the un-resummed logarithms such as [50,53]. Our model at hand involves multiple degrees of freedoms, therefore, there is no single choice of the scale to make all the logarithms negligible. In this work, we choose where i runs over all degrees of freedom (mass eigenstates) in the plasma. This is a convenient choice as long as there is no large separation between scales of the particles' masses, including the thermal mass contribution, as well as between the particle masses and the temperature, as it is the case in our study. The CW potential further includes polynomial contributions of the radiative corrections. It also partially accounts for multi-scale particle threshold effects beyond the one single scale threshold taken into account through the RG improvement. We collect the one-loop beta functions and wave function renormalization factors for our model in appendix section D and implement the RG improvement for all numerical calculations 4 .

IV. MEAN FIELD ANALYSIS FOR THE THERMAL HISTORY
This section provides an analytical understanding of the model parameter space compatible with the desired thermal history -the electroweak symmetry stays nonrestored in the inert sector up to temperatures much higher than the EW scale, whereas the agent of the electroweak symmetry breaking changes at temperatures around the EW scale from the inert Higgs sector to the SM one. In this work, we do not explicitly discuss the UV scale physics completion that may lead to electroweak symmetry restoration at even higher energies and hence would allow for the possibility of EWBG. However, we will study the conditions necessary for the suppression of the sphaleron rate as a function of the model parameter space through the whole temperature regime for which the electroweak symmetry is broken. More specifically, we will explore the constraints on the ratio between the electroweak symmetry breaking vevs to the temperature that may allow for such a suppressed sphaleron rate. This will provide a framework for future EWBG model building. If, instead, the new physics UV completion would directly provide a source of baryon asymmetry at the high scale, such as, for example, in the case of Leptogenesis, GUT-baryogenesis or Affleck-Dine baryogenesis [60], then the requirement on the sphaleron rate could be ignored. A discussion of possibilities for baryogenesis as well as specific details on the sphaleron rate relevant for our model will be presented in section VI.
We summarize the above desired thermal history with three conditions as follows • C1: Non-restoration of the electroweak symmetry This is realized up to very high temperatures by having a non-trivial inert phase: ϕ highT = 0; • C2: Phase transitions from the inert Higgs phase to the SM Higgs phase This condition secures that the Universe is at the SM vacuum at zero temperature, while being compatible with C1.
• C3 5 : Sufficiently suppressed sphaleron rate after EWSB This would allow preserving any baryon number density that may be generated through an EWBG mechanism at the ultraviolet.
To gain an analytical understanding of the model parameter space compatible with the above conditions, we use a mean-field approximation of the finite temperature effective potential, where the thermal potential is evaluated up to leading order of the high-temperature expansion where Such a mean-field potential is a reliable approximation before considering RG improvement and daisy resummation, especially at high temperatures. We shall include resummations in the next section for a full numerical study at high field values and temperatures. Here we provide an analytical study based on the mean-field potential to obtain a coarse understanding of how the desired thermal history is achieved within our model. Let us first study the SM and inert Higgs sector phases of the potential in eq. (20) . An inert phase P Φ , where only the inert Higgs field has a non-zero field value, reads with At very high temperatures, T 2 µ 2 Φ , one can approximate Given the BFB condition that λ Φ > 0, a negative thermal mass coefficient c ϕ , 5 As discussed above, this condition is optional.
generates a non-zero inert phase at very high temperatures, which is the key to achieve electroweak symmetry non-restoration (or delayed restoration) in the inert sector in our model. This provides for condition C1 in the mean field approximation as The main driver of a negative c ϕ is a negative cross quartic between the inert and the singlet sector λ Φχ , whose negative contribution is magnified by the number of singlets N . If the inert mass parameter µ 2 Φ ≥ 0, such a phase where only the inert field has a non-zero vev would disappear at a temperature T r Φ (either as a global or local minimum), where A low restoration temperature T r Φ facilitates the existence of phase transitions between the inert and SM Higgs phases as well as the associated condition for a suppressed sphaleron rate, which will be discussed in more detail below. Instead, if µ 2 Φ < 0, this inert phase exists at zero temperature, which puts a constraint for it to be above the EW vacuum at T = 0, i.e. V 0 (0, w(0), 0, · · · 0) > V 0 (v 0 , 0, 0, · · · 0), in addition to condition in eq. (8).
A SM Higgs phase P H of the potential, where only the SM Higgs has a non-zero field value, reads where at zero temperature it becomes the EW vacuum with v(0) = v 0 . Such a phase appears at a temperature Another phase that possibly exists during the thermal history is when both the SM Higgs and the inert Higgs fields acquire simultaneously non-zero values implying that this phase is governed by the Higgs-Inert mixing coupling Λ HΦ defined in eq. (6). An important feature of this phase is that given the potential in eq. (20), the potential difference reads where the proportionality coefficients are always positive independent of the temperature. Thus, if 4λ Φ λ H −λ 2 HΦ ≤ 0, the Higgs-inert phase P HΦ is irrelevant as it is always shallower than either the SM or inert Higgs phases. On the contrary, if 4λ Φ λ H − λ 2 HΦ ≥ 0, as long as such a Higgs-inert phase exits, it is deeper than both the SM or inert Higgs phases, thus becoming the global minimum.
Concentrating on the case where P HΦ is the global minimum at a given temperature, notice that the situation 4λ Φ λ H − λ 2 HΦ ≥ 0 coincides with the BFB condition if λ HΦ ≤ 0, hence for negative/zero cross quartic, the Higgs-inert phase will be the global minimum at finite temperature. Moreover, at zero temperature, the non-tachyonic condition enforced in eq. (8) implies (32), as expected since the non-tachyonic solution was derived under the assumption that the P H at T = 0 being the physical vacuum. In addition, let's recall that at very high temperatures we have restricted our case to the inert phase P Φ being the global minimum (no electroweak symmetry breaking in the SM Higgs sector), hence eq. (32) implies that we voluntarily enforced whenever T 2 µ 2 H(Φ) . Given the above constraints (P H and P Φ are the global minimum at T = 0 and high temperatures, respectively), if the phase P HΦ ever appears, in a temperature regime T 2 ∼ µ 2 H(Φ) , it develops at a temperature Max{ T r H , T r Φ } and must disappear at a lower temperature Min{ T r H , T r Φ }. These two characteristic restoration temperatures are defined from eq. (32) demanding that either v( T r H ) = 0 or w( T r Φ ) = 0, respectively. Observe that, within the mean field approximation we are considering, from eqs. (22), (29) and (32) , which implies that, the critical temperature defining the transition between the P HΦ and P Φ (P H ) phases is given by indicating that these transitions are second order within the mean field approximation. These equalities imply that, and their existence demands In the numerical study where we consider the full thermal potential as well as daisy contributions, such phase transitions could be affected and become first order. However, they would hardly be strongly first order in the absence of large thermal or tree level barriers.
Other possible phases associated with the finite temperature potential (20) include the trivial point, which, as long as any of the above phases exist, yields a shallower value of the potential, as well as phases involving singlets with non-zero field values. The latter will not be further considered in this section as they are unlikely to participate in the thermal history. When evaluating the thermal history in the numerical section, however, all possible phases will be taken into account.
After having considered the existence of all possible phases and some of their properties, let us now concentrate on the specifics of the phase transitions from the inert sector to the SM Higgs sector.
First, we discuss the simpler case where the phase P HΦ either never appears or is irrelevant. In such a case, there should be a phase transition from P Φ to P H , as illustrated on the middle penal of the second row in Figure 1. Given the potential in eq. (20), such a transition happens at a critical temperature at which the potential becomes degenerate The condition for such a T c to exist reads (with help from the zero temperature constraint eq. (8)) and this will have a relevant impact on the allowed values of the inert Higgs boson mass, as will be discussed later on.
As mentioned in C3, to allow for the possibility of a EWBG after UV completion, we will look at the conditions on the sphaleron rate at finite temperatures. The dilution of the baryon number density after the onset of a UV induced EWPT responsible for the EWBG will be double exponentially suppressed by the ratio of the sphaleron energy to temperature, see discussion in section VI. Hence successful EWBG in the complete model will require (see e.g. [44]) where ϕ and h are the inert and SM Higgs fields charged under the EW gauge group. This condition should be satisfied at any temperatures throughout the thermal history from the creation of baryon asymmetry up to present times. It can be shown that such a condition can be satisfied if the phase transition P Φ to P H fulfills This follows from the fact that as long as µ 2 Φ ≥ 0, as will be implemented in our BM scenarios, Next, we discuss the case where the phase P HΦ is relevant and appears as a global minimum in the thermal history, as illustrated on the left panel of the second row in Figure 1. To have a two step phase transition near the EW scale one needs which corresponds to the condition for these two temperatures to exist given by eq. (37). Analogous to the previous case, the condition to avoid baryon asymmetry washout in the context of a EWBG in an UV completed theory, would require Another thing to notice in this case is the role played by the mixing quartic λ HΦ , which controls the deviation from T c H to T r H and from T c Φ to T r Φ . The smaller the mixing quartic, which is the region that we are mainly interested in, the smaller the deviations are. Moreover, in the region of small λ HΦ , the phase transition pattern P Φ → P HΦ → P H is most likely to happen due to the decoupled contributions from the inert and SM Higgs minima to render the P HΦ minimum in the intermediate temperature range. This is apparent in Figure 3 to be discussed below.
It is also possible to have a temporary electroweak symmetry restoration at temperatures between those supporting the two EW breaking phase structures P Φ and P H . This is the case when T r Φ is higher than T r H , as illustrated on the right penal of the second row in Figure 1.
Since in the temperature range between T r H and T r Φ the system is in a EW restoring phase, this scenario would allow for the EW sphaleron to be active in this regime. The sphaleron will wash out any baryon asymmetry that could have been generated by high scale EWBG. At this moment, we will mainly focus on the previous cases that are compatible with an UV EWBG mechanism.
Another possible case is a more fine-tuned four-step phase transition when T c H ≤ T c Φ and T r H ≥ T r Φ . This case will require large mixing quartic and significant finetuning of the parameter space. We do not further concentrate on this case.
In Figure 3, we show the parameter space spanned by N λ Φχ − λ HΦ considering the zero temperature constraints discussed in section II and the different thermal history possibilities discussed above. The region violating condition C1 MF is shaded gray, while the regions satisfying the thermal history patterns and the nonwashout conditions are highlighted with light and dark orange (light and dark maroon) for the transition pattern P Φ → P H (P Φ → P HΦ → P H ), respectively. There is no region that satisfies the rare four-step phase transition. The conditions for the correct zero temperature vacuum structure are satisfied on the whole parameter space if we impose µ 2 Φ , µ 2 χ ≥ 0. The region giving a treelevel BFB potential, calculated from conditions (5), is at the right side of the black solid lines for different number of singlet scalars N . Notice that, within the meanfield approximation, the thermal history patterns, as well as the non-washout requirements, are independent on N as long as the value of N λ Φχ is kept a constant. Since both the thermal histories and non-washout conditions are strongly correlated to the inert mass parameter, the mass of the inert Higgs boson is in turn also constrained. In Figure 3, we show solid blue lines that determine the maximal value of the inert Higgs boson mass compatible with the corresponding phase transition patterns for a given value of N λ Φχ and λ HΦ . Higher values of the inert Higgs boson mass can be achieved to the left of the lines. Similar lines for the suppressed sphaleron rate conditions are shown by the dotted blue lines. Other parameters have been fixed in Figure 3 to be λ Φ = 0.1, N Λ χ,n = 2.5 and λ Hχ = 0, λ HΦ = 0. The SM Higgs sector parameters are fixed to satisfy the Higgs vev and mass at the tree level. We constrain the discussion to the case µ 2 Φ ≥ 0, which makes conditions (41) and (47) sufficient to secure a suppressed sphaleron rate within the mean field approximation as discussed above. In addition, in Figure 3, we also show the two benchmark points A and B 6 , which will be discussed in the full numerical study in the next section.
From Figure 3, one notices that the region where the cross quartic coupling between the inert and the SM Higgs sectors almost vanishes, i.e. λ HΦ ∼ 0 and hence the SM Higgs sector is minimally perturbed, can be compatible with the desired thermal history. Main constraints on the parameter space come from the tension between the BFB and desired thermal history: the more negative the cross quartic N λ Φχ , the easier the non-restoration and the lower the critical temperatures which yield larger EW vev to temperature ratios ξ(T ). A more negative cross quartic coupling N λ Φχ makes it harder for the potential to be BFB, as shown in eq. (5). Moreover, a larger number of singlets in turn helps to relax the BFB condition on N λ Φχ by relaxing its lower bound while increasing the singlet effective quartic coupling N λ χ . As mentioned above, another constraint is on the mass of the inert Higgs boson. The restriction on the parameter space is alleviated for a lighter inert Higgs boson mass, especially in the region where the cross quartic λ HΦ is small. This can be easily understood, for example in the P Φ → P H phase transition pattern, since a smaller inert mass parameter µ 2 Φ yields a lower critical temperature T c , as is shown in eq. (38). A similar argument, although more involved, applies to the two-step phase transition. The direct correlation between a smaller inert Higgs boson mass and a smaller inert mass parameter µ 2 Φ especially holds in the region of small λ HΦ , as the one considered here. Observe however, an inert Higgs boson mass above half of the SM Higgs mass can be achieved, even with λ HΦ ∼ 0, as long as the number of singlets is sufficient to be in the BFB allowed region.
The analysis in this section is based on the meanfield approach, where we consider the leading order hightemperature expansion of the thermal potential. For temperatures well above the EW scale however, including the RG improvement and the daisy resummation becomes necessary. In the next section, we perform a full numerical study for two benchmark points and present the results for different approximations.

V. NUMERICAL RESULTS ON BENCHMARK POINTS
In this section, we explore the thermal histories of two model benchmark points based on numerical calculation of the finite temperature effective potential prescriptions as described in section III. Appendix section F contains the details of our python implementation. The main result of the algorithm is the value of the global minimum at a given temperature. The set of all global minima at a given set of temperatures defines the phase history we consider. A phase transition is observed when there is a change of phase pattern (e.g. from an inert-only-phase to an inert-SM Higgs phase) at a certain temperature 7 . In  this section, we also explore the value of the EW vev to temperature ratio ξ(T ), which is relevant for obtaining information on the sphaleron rate. We define two characteristic benchmark points for our model -benchmark A that has inert Higgs eigenstates with masses slightly above half of the Z boson mass, and a benchmark B that has inert Higgs eigenstates with masses slightly above half of the SM Higgs boson mass. Inert mass eigenstates with masses above 100 GeV can be achieved, but they would either lead to restoration of the electroweak symmetry at intermediate temperature scales or would require a number of singlet scalars of order O(1000) or more. The specific values of the model parameters and masses are given in Table I.
We implement the RG improvement on the BFB conditions of eq. (5) and find that, at scales of the order of 10 5 GeV, these conditions are violated for both BMs 8 . Such an energy scale is of the order of the scale at which the SM Higgs quartic coupling becomes negative through its SM one loop RGE. This is expected since we consider that the SM Higgs only interacts with the extended scalar sector through a tiny inert-Higgs coupling, and therefore its quartic coupling evolution should be minimally perturbed compared to its SM behavior. The scale above gives a rough estimate of the energy scales up to which our results can be trusted. By minimizing the finite temperature potential numerically, we have checked that the potential remains stable up to high energy scales shown below for both BMs. We also ran the RGE of the model (see eqs.(D.1)) for BMs A and B and found that Landau poles appear at energies around 2.5 · 10 14 GeV and 10 15 GeV, respectively -well above the scale of validity of the theory at the one-loop RGE level.
In Figure 4 and Figure 5, we show the phase structure (upper panel) and EW vev-temperature ratio (lower panel) for BMs A and B, respectively, and for different implementations of the finite temperature effective potential as introduced in section III. In the phase structure plot, we are showing as a function of the temperature the field values of the SM Higgs (red), inert Higgs (blue), and singlet (green) 9 at the global minimum. In the vev- 8 There is a small dependence on the CW treatment that somewhat perturbs the SM Higgs quartic coupling as is explained in section III. 9 We assume all singlets have the same vev -it's either all or temperature ratio plot, we show the value of ξ(T ), as defined in eq. (40), as a function of the temperature. To showcase the uncertainties associated with different finite temperature implementations, we show results obtained with no daisy resummation (solid lines), daisy resummation with high-T thermal masses, as in eq. (17), (dashed lines), and daisy resummation with truncated thermal masses, as in eq. (16), (dashed-dotted lines). In addition, we have included the RG improvement for all calculations, and consider the uncertainties related to the CW potential, which takes care of multi-scale issues beyond the RG improvement. In the figures, we use the same type of lines to represent a given finite temperature approximation with or without the CW contribution. Hence the space in between the lines shows the uncertainty related to the CW effects. It is apparent from the figures that this accounts for a small effect, and we will not discuss it any further.
In Figure 4, for BM A, one observes that the major uncertainty is caused by the effects of daisy resummation and the impact of different thermal mass treatments within the daisy resummation. However, the most important feature of these results is that the qualitative behavior of the phase structure, in Figure 4 upper panel, and the EW vev-temperature ratio affecting the sphaleron rate, in Figure 4 lower panel, is not significantly modified by the different finite temperature treatments. In fact, BM A exhibits both the feature of EW non-restoration until high energies and ξ(T ) > 1. The plots of BM A are shown up to the temperature of 10 5 GeV, after which the potential becomes unbounded from below. Observe that there is a kink below/about 200 GeV, which is due to the phase transition pattern from the P Φ phase to the P HΦ phase, and it is a physical effect. In addition, there is a spike at T ∼ 110 GeV for the daisy resummation with truncated thermal masses, which is, however, a defect of this finite temperature implementation. We expect this effect to be smoothed out when implementing an improved treatment of the thermal masses 10 .
In Figure 5, we show similar results as for Figure 4, but for a heavier inert Higgs boson mass of the order of m h /2 that will allow for different phenomenology. Same as BM A, BM B exhibits both the feature of EW non-restoration until high energies and ξ(T ) > 1. The plots are shown up to the temperature of T = 4 · 10 4 GeV, after which the none. Given that λχ = 0, which we chose at tree level and is protected against RGE, we have the SO(N ) symmetry that we can use to rotate in that form.

FIG. 4.
Phase structure (upper panel) and EW vevtemperature ratio (lower panel) as a function of temperature, for different finite temperature implementations, for BM point A as defined in Table I. The singlets never obtain a non-zero vev in the entire temperature range. The black dotted line separates linear from logarithmic axes scales. The red dotted line in the lower panel marks ξ = 1. potential becomes unbounded from below. For the BM B, there is a kink above/about 100 GeV, which is due to the phase transition from the P HΦ phase to the P H phase. In addition, analogs to BM A, there is a spike at around 300 − 400 GeV for the daisy resummation with truncated thermal masses, which we understand is the same type of artifact as discussed above and will be cured by implementing an improved treatment of the thermal masses.
As described above, using the gap equation, eq. (15),  Table I. The singlets never obtain a non-zero vev in the entire temperature range. The black and red dotted lines are as in Figure 4.
to derive the thermal masses is the most robust procedure. However, solving the gap equation at every step in the minimization of the potential is computationally extremely expensive and is beyond the scope of this work. However, in order to secure that the non-restoration behavior at high temperatures survives the most precise treatment of the thermal masses through the gap equation, we checked for several high-temperature values all the way down close to the EW scale, that the nonrestoration behavior and ξ(T ) > 1 survive for both BM scenarios.

VI. BARYOGENESIS AND SPHALERON RATE SUPRESSION
In this section, we briefly discuss the possibilities of high-scale baryogenesis scenarios based on 1) EWsymmetry non-restoration up to scales as high as the GUT/Planck scale or 2) electroweak symmetry restoration around a UV scale of the order of validity of our model at which a new UV theory is in place. In the latter case, we expect to build a UV theory that allows for EWBG. Hence, in this case, we would like to explore in more detail the sphaleron washout constraints in our BM scenarios to preserve the created asymmetry down to zero temperature.
If the EW symmetry, through a specific UV completion, were to remain broken well above the scale of validity of our current model 11 , possibly up to the GUT or Planck scale, this would enable baryogenesis mechanisms with little dependence on how the EWSB is triggered. In such case, the baryon asymmetry can be generated by a mechanism that creates a source of B-L = 0, such as, for example, GUT-genesis, leptogenesis, or Afflect-Dine baryogenesis ( [61] and references therein). Recall that, sphaleron processes preserve B-L, and hence an asymmetry will subsist once generated. However, they tend to wash out B+L as long as they remain active, thereby enabling conversion of Baryon (anti-Baryon) number into anti-Lepton (Lepton) number. For any specific B-L = 0 mechanism, there will be additional model-building considerations for successful baryogenesis, including specifics of the new sources of CP violation and out of equilibrium conditions. It is important to notice that the two BMs we consider in this work imply that the sphaleron rate is suppressed during the broken-electroweak symmetry epoch, hence a mechanism such as Leptogenesis, that requires active sphalerons to convert Leptons into antiBaryons will not work. Other BMs could be studied that allow for sphalerons to become active at some intermediate energy scale during the temporary restoration of the electroweak symmetry, as in the lower right panel of Figure 1. Exploring these new ideas for baryogenesis will be the subject of future work.
In the case of a UV completion that induces a restoration of the electroweak symmetry at high energy scales of the order of the validity of our model, one can also require that such UV theory induces a strong first order phase transition and enables EWBG. Although building such UV theory will remain a topic of future work, let us briefly comment on the various ways that we can picture such a scenario. In our minimal model, the restoration can occur through the RGE of the quartic couplings. Under the 11 Our study is only including one-loop RGEs, but in analogy to the SM, we expect the validity of our model to be extended to higher energies by considering higher-order loop RGEs.
high-temperature expansion, one can visualize this possibility through the thermal coefficient c ϕ given in eq. (24). If c ϕ , which at lower temperatures has a negative value, were to become positive at a given high scale through the RGEs, this will render EW restoration at high temperatures. For simplicity, let's consider the limit where in the IR the mixing quartics λ HΦ , λ HΦ , and λ Hχ are zero, and neglect the leading log impact of these mixing quartics. The running of the thermal coefficient is then determined by the running of the linear combination of λ Φ /2 + (3g 2 + g 2 )/16 + N λ Φχ /24. In our model the inert doublet self-coupling λ Φ generically becomes larger at higher scales, while the mixing quartic λ Φχ , whose initial value is negative, could also increase, depending on the specific region of parameter space. However, we checked that the latter is not fulfilled for our BMs, hence, additional effects will be needed to restore the electroweak symmetry in these cases. There are indeed different ways to change the running behavior of c ϕ , to allow for EW restoration. For instance, one can consider that the inert doublet is charged under some new spontaneously broken U(1) gauge group with coupling g . This will affect c ϕ directly by adding a g 2 /16 term after crossing the scale where the new U(1) is restored, rendering its gauge boson massless such that it starts contributing to the thermal mass of the inert doublet. Similarly, one can also introduce some heavy vector-like fermions (under SM gauge groups) that have Yukawa couplings to the inert doublet. When above the heavy fermion mass scale, this will add new positive contributions to the thermal coefficient c ϕ by y 2 F N F /12, where N F is the color-factor or the specifies of new heavy fermions.
Beyond directly changing the thermal coefficient c ϕ above some mass threshold scale, one can also modify the running of the couplings contributing to c ϕ , by adding new gauge and/or matter content. The minimal and simplest implementation would be to charge the χ field under some new SU (N ) gauge group. It directly contributes positively to the beta function of λ Φχ , which is the only source of negative quantities in the thermal coefficient c ϕ , helping restore the EW symmetry at a higher scale. On the contrary, matter fields interacting with the inert Higgs field seem to contribute negatively to the beta functions of the quartics contributing to c ϕ , although, as discussed below, they may be required to secure a strong first-order phase transition. An additional source of symmetry restoration could be to add scalar fields that directly couple to the inert field and acquire masses at high energies at which restoration would take place [23].
An important additional issue related to the high energy EWBG mechanism in the framework of delayed electroweak symmetry restoration, is that one needs to secure that a strong first-order phase transition takes place at the time of electroweak symmetry breaking. This is required by the out-of-equilibrium condition of Sakharov. Here, it is possible to exploit the existence of an inert fermion sector that suppresses the strength of the inert self-coupling and thereby enhances the strength of the phase transition.
As it is clear from the above discussion, a successful UV model of high-temperature EWBG will demand detailed model building, which we leave for future publication. In the following, we will concentrate on the EW nonrestoration case at hand, where the SM Higgs sector is minimally perturbed, to discuss details of the sphaleron rate.
Once the UV completion allows for the creation of the baryon asymmetry through an EWGB mechanism at high temperatures, one needs to evaluate the sphaleron washout factor to preserve the asymmetry all the way down to zero temperatures. Our model generically predicts a slowly varying ξ = v EW (T )/T up to high temperature, as well as a low scale phase transition between the inert doublet and the SM Higgs doublet phases near the weak scale. Following a high-scale SFOPT triggered by a UV completion of the model, the sphaleron will become inactive quite fast after that transition, but there will be some dependence on its rate on the model parameters. To properly compute the washout (dilution) of the baryon number density, one should integrate the effects of the sphaleron rate over a large range of temperatures (a large period of time), instead of the usual assumption that the washout factor is dominated near the vicinity of the phase transition and is treated as a constant.
Specifically, the amount of sphaleron induced washout is determined by two quantities: the product of prefactors entering in the sphaleron rate and the energy of the sphaleron that appears in one of the exponentials. The latter is straightforward to compute and largely depends on the gauge structure of the theory. We provide the necessary steps to get the sphaleron energy [62,63] in detail in appendix section E. The computation of the prefactors of the sphaleron rate are more model dependent. There are two different sources of deviations from the SM results. First, we have an extended scalar sector and the additional particles might contribute through indirect effects in the prefactors. Second, we focus on the inert doublet, and its quartic coupling is different than the quartic of the SM Higgs. We therefore discuss the specifications of the prefactors from the SM values [64,65] in detail in appendix section E.
The sphaleron rate can be written as The survival rate of the baryon number density at any given time t, after the onset of the transition at t = 0 is [40,44,66] where we consider present time, t = t now , and with n f the number of fermion families. In a radiation dominated Universe, changing the integration variable from time to temperature, the above equation reads where M P l is the Planck mass and g * is the number of relativistic degrees of freedom. In our case, it is g * = 106.75 + 4 + N for the range of temperatures under consideration.
Based on the calculation presented above, we can define the washout or dilution factor as f w.o. = 1− n B (tnow) n B (t high ) . In Table II we show the values of f w.o. for our two benchmarks. We see that BM A has a negligible washout factor for all choices of parameters, even for the most aggressive assumption for the fluctuation determinant κ, which is a factor 100 larger than the value suggested in [67]. BM B shows sub-percent or even negligible washout for the majority of approximations. Only for the most aggressive choice of κ, we observe values that can be as high as 70%, which can be compensated by producing an initial asymmetry about three times larger than the asymmetry we observe now. If we consider the central value for κ, we see a washout of at most 1.2%. We also note that the washout factor governed by eq. (50) is much less sensitive to T high than to effects at temperatures close to the EW scale. This is the case since at higher temperatures the double exponential in eq. (50) is larger than at EW temperatures. Indeed, at temperatures around the EW scale, there is an enhancement from the inverse of the Hubble expansion rate, as well as from the exponent proportional to exp [−ξ(T )], where ξ(T ) has its lowest values. That makes the double exponent in eq. (50) to take its smallest values for temperatures close to the EW scale. Hence at such temperatures is when the main effect of the exponential washout takes place. In other words, the relevant contribution to the washout factor is only at scales between the EW scale and around 500 GeV, while at high temperatures the exponential washout remains negligible. This holds as long as ξ(T ) does not fall fast below 1 at high temperatures, which is the case for our BMs. This ensures that high-temperature EWBG could build in through a proper UV completion of our model.
The possible existence of light scalars, Φ and χ i will open the possibilities of the SM Higgs decaying into invisible particles, via the generic portal couplings The generic Higgs decay width into new scalars via this portal coupling is (per scalar degree of freedom): where the coupling λ Hs can be one of the above quartics, λ HΦ , λ HΦ or λ Hχ , and m s can be the mass of the Φ or χ states, respectively. The current LHC 95% confidence level (C.L.) limit on Higgs invisible decays is 11% [68] and the HL-LHC projection is 5.6% [69]. When m s m h , the phase space suppression is negligible and this translates into an upper limit for the SM-new scalars mixing quartics. The current and future limits on the mixing quartics read N λ 2 Hχ + 2(λ HΦ + λ HΦ ) 2 + 2λ 2 HΦ 0.010 (0.007) (53) for LHC (HL-LHC). In the above, by including 2λ 2 HΦ , we also include the Higgs decays into a pair of the charged states from the inert doublet.
In the absence of other mass splitting generating interactions, e.g. λ HΦ being zero, one-loop SM effects generate mass splittings between the charged and neutral eigenstate of the inert doublet of about 360 MeV [70]. The charged state will decay back to the neutral state via a soft charged pion, or via the three-body decay mediated by an off-shell W boson. The typical lifetime is independent of the inert doublet mass and is a few mm.
Hence, this charged state can also be treated as invisible at colliders. In fact, precision Z boson measurements of its invisible decays exclude all inert masses below 45 GeV, and hence we shall only consider inert masses beyond the 45 GeV value [71].
Still, one can attempt to look for signals beyond the missing energy at colliders. At high energy colliders, such as the LHC, although challenging, one can look for the disappearing track signatures from the charged eigenstate of the inert doublet. However, it is well-known that this channel is difficult for Higgsinos, to which our inert doublet model signature resembles most. The current sensitivity from LHC disappearing track searches can exclude pure Higgsinos up to 78 GeV [72]. The inert doublet production rate from the Drell-Yan process is roughly a factor of four lower than that of Higgsino production, due to the inert charged Higgs being a scalar rather than a fermion. Furthermore, for small mixing quartics such as λ HΦ , one can arrange additional contributions to the mass splitting between the neutral and charged inert doublet states. This will make the charged state decay promptly and therefore the disappearing track searches will no longer apply. Given the above, we are entitled to ignore the disappearing track search limits and only comply with the LEP Z invisible bounds for our benchmark scenarios. Future tests on disappearing tracks could still shed light on our model. Summarizing, considering direct search constraints for our electroweak symmetry non-restoring model, we observe that the mixing quartics λ HΦ and λ Hχ are bounded by constraints on invisible SM Higgs decay rates. This can give a strong handle for testing possible benchmarks, but at the same time there are models, like our BMs, in which they happen to have neglibible values. In this sense, the more direct and inevitable probe for our model at colliders are through the invisible Z decays, relying only on the gauge coupling structure. Disappearing charged track searches open a new window of opportunity, if not undermined by parameter choices of the various mixing quartic couplings.
There are additional U(1) global symmetries in the inert sector Φ as well as Z 2 symmetries under which the singlet fields χ i are odd, that prevent direct mixings between these states with our SM Higgs doublet. There are, however, loop-induced corrections to the SM that can be probed through precision observables. The leading contribution to the electroweak precision observables (EWPO) is from the custodial symmetry breaking term λ HΦ , inducing an operator contributing to the Tparameter [73] For an inert doublet mass scale µ Φ around half the Higgs mass, the EW precision measurement constrains the Tparameter with uncertainty 0.07 [74,75], constraining | λ HΦ | < 0.36 at 95% C.L. Although this estimation is subject to sizable corrections due to the fact that µ Φ is of the order the Higgs mass, this gives an estimate of the bounds on λ HΦ not being very stringent coming from one-loop suppressed effects. For our benchmarks, we simply set λ HΦ to zero at tree level. The next set of constraints comes from the Higgs boson coupling precision measurements, through the coefficient of the operator, This results in an overall reduction of the Higgs couplings by 1/2c H v 2 0 . We note here that this EFT matching is subject to large corrections and higher-order terms since the scales µ 2 Φ and µ 2 χ are not far from the Higgs mass squared. On the other hand, our non-restoration mechanism has limited dependence on these parameters. In particular, we have set λ HΦ and λ Hχ to be zero in our BM scenarios, leaving only a shift of the Higgs couplings of about −1/2c H = −λ 2 HΦ v 2 0 /(96π 2 µ 2 Φ ). For an inert doublet with µ Φ around half the Higgs mass, it yields a global shift in the Higgs couplings of around −λ 2 HΦ /(6π 2 ), bounding |λ HΦ | < 1.1 (at 95% C.L.) if we were to achieve 1% Higgs coupling precision at the HL-LHC [76]. This constraint is much weaker when we compare it to bounds from direct invisible Higgs decay searches discussed earlier in this section. It could however be relevant for scenarios with heavy inert masses, since the Higgs invisible decay bound no longer applies. However, in such case, as we shall see next, the precision measurements on Higgs to diphoton coupling provide a stronger constraint than the one derived from eq. (55).
The EW charged inert doublet also radiatively modifies Higgs couplings to EW gauge bosons, through Here τ a are the SU (2) generators. Consequently, the Higgs diphoton coupling is modified by where κ γγ ≡ g hγγ /g SM hγγ . Due to the fact that the SM Higgs to diphoton coupling is loop-induced, this provides a strong constraint on |λ HΦ | to be smaller than 0.04 (at 95% C.L.) for a 1.9% precision [76] on the Higgs to diphoton coupling at HL-LHC. The current Higgs precision uncertainty of 17% [77] translates to a constraint on |λ HΦ | < 0.4 (at 95% C.L.). Again, in deriving this limit, we assume that λ HΦ = 0, µ Φ being half the Higgs mass, and ignore the deviation of the form factor from unity from the inert doublet running in the loop.
Beyond the above, the model also generate less constraining effects on EWPO (W and Y parameter) and Higgs self-coupling [73,78], whose current and future perspective sensitivities can be found in Refs. [78][79][80]. This may provide, in the future, further complementary information about the model.

VIII. CONCLUSION
The exploration of electroweak phase transition patterns leading to electroweak symmetry breaking allows us to envision plausible paths for EWBG, as well as details of the cosmological history of our universe. In particular, the possibility of electroweak symmetry non-restoration up to high energy scales, conceivably up to the GUT or Planck scale, or the opportunity for delayed electroweak symmetry restoration up to scales of the order of 100 TeVs, opens new windows for baryogenesis mechanisms. In this paper, we propose a novel approach to realize new thermal histories, by enabling the agent of EWSB to be an inert doublet that yields electroweak symmetry nonrestoration up to high temperatures. These possibilities allow for diverse thermal histories with multi phase transition patterns, involving the SM Higgs, the inert Higgs and the SM-inert Higgs mixing phases at finite temperatures.
Our new approach for electroweak symmetry nonrestoration at high energies has interesting computational requirements. Since the thermal history of our model, as defined in section II, spans over large scale separations from the EW scale to high temperatures, in our study we carefully implement the effects of RGE and thermal resummation, as detailed in section III. When considering daisy resummation, we compute thermal masses with different approximations and observe that they lead to similar quantitive results. In section IV we perform an analytical study at leading order in the high-temperature (mean field) approximation that helps us zoom in into the promising region of parameter space for our numerical study. In section V, we report our numerical calculations for two benchmark points, and show that our results are robust under various treatments of thermal resummation while including RGE effects. Most importantly, the non-restoration patterns can hold at least up to high scales of the order of 10 5 GeV, within the one loop RG resumed effective potential. An UV completion of our model could take place at higher energy scales. In section VI, we present a detailed study of the sphaleron washout effects over a broad range of temperatures, and show that for our two benchmark scenarios, the washout rates are such that high temperature EWBG could be realized after a proper UV completion. Observe that the crucial ingredient of our BM scenarios is that the EW symmetry is non-restored from high temperatures all the way down to the EW scale.
Most importantly, our mechanism for transmitting broken electroweak symmetry from the SM sector to an inert sector has a specific interesting feature: It can work even if one decouples the two Higgs sectors in the tree level scalar potential, implying that the effect of the new doublet enters our zero-temperature particle physics tests at the electroweak-loop level. This enables the existence of large model parameter space compatible with experimental constraints and at the same time calls for new precision tests of the SM. As discussed in section VII, our model will find scrutiny at the HL-LHC through electroweak and Higgs precision tests, invisible decays and searches for disappearing tracks.
At high temperatures, our model opens up to possible UV completions that would enable various baryogenesis mechanisms. If we go through EWBG, where a strong first-order electroweak phase transition is necessary, it would give rise to gravitational wave signals. The peak frequency, instead of populating around the LISA band (mHZ), will increase to higher frequencies, at reach of facilities [81,82] such as BBO, DECIGO, and even aLIGO. Moreover, the additional singlets χ in our study can themselves go through phase transitions, further enriching the possible thermal histories of our universe. Beyond all the above, one can also explore such relay of the EW-broken phase between the SM Higgs and scalars under other EW representations.
Note added: During the completion of this work, [89] appeared and considered a specific realization of symmetry non-restoration in a scenario with a 2HDM and one singlet scalar. We note that the main cause for us to require more singlet scalars is to strictly forbid the EW restoration at low temperatures, and thereby avoid the situation as depicted in the bottom right panel of Figure 1. Furthermore, additional number of scalars are needed for the theory to obey tree-level perturbative unitarity up to the high temperature scales of nonrestoration. In particular, our BM scenarios satisfy unitarity up to 10 14 GeV. To the best of our understanding, in [89], only benchmark E1, plus gray points in Fig. 7, survive our requirement of no temporary restoration at low temperatures, but due to perturbativity, the validity of the model appears to be limited to scales not far above the EW scale 12 .
Lastly, the principle submatrices in the form of The corresponding eigenvectors read {{0, 0, −1, 1, 0, . . . , 0}, · · · , {0, 0, −1, 0, . . . , 0, 1}, {x 1 , y 1 , 1, · · · , 1}, {x 2 , y 2 , 1, · · · , 1}, {x 3 , y 3 , 1, · · · , 1}} , (B.23) where x i and y i are given by the solution of Notice that the eigenvectors of the n − 1 eigenvalues λ χ are non-positive, and accordingly they do not give any constraints on the copositivity of the matrix. For the last three eigenvalues e i with i = 1, 2, 3, the last entries of their eigenvectors is unity, and accordingly, the positivity of the eigenvectors is determined by the sign of x i and y i . Thus, such (n + 2) × (n + 2) matrices are copositive if and only if there is no e i < 0 with corresponding x i , y i > 0. Now let's use the fact that the roots {e 1 , e 2 , e 3 } are also eigenvalues of a 3 × 3 matrix with the corresponding eigenvectors being where x i , y i are solved by the same conditions given in eq. (B.24). Since the last entry of the above eigenvectors is positive, one can immediately see that the condition of the copositivity of such a 3 × 3 matrix is identical to the conditions for the above (n + 2) × (n + 2) matrices: there is no e i < 0 with corresponding x i , y i > 0. The copositivity of a 3 × 3 matrix in terms of its entries has been discussed in the literature. The conditions are [92]: Now, we have derived the copositivity conditions for all type of principle submatrices of the quartic coupling matrix of our two doublets + N singlets tree-level potential. Notice that we work in a generic basis including all CP even, CP odd and charged components of the doublets. The potential will be bounded from below if all the conditions are satisfied. We define for convenience: with which the BFB conditions can be written as  For later convenience, we define constants c i = Π 2 i /T 2 ĥ =φ=0 . The thermal masses of W and Z are as given in [90,91], they only contribute to the longitudinal components: Π 0,W L = 2g 2 T 2 (C.4) Π 0,Z L ,A L = − g 2 + g 2 8 (ĥ 2 +φ 2 ) + (g 2 + g 2 )T 2 ± ∆ (C.5) ∆ 2 = g 2 + g 2 8 2 ĥ 2 +φ 2 + 8T 2 2 − g 2 g 2 T 2 ĥ 2 +φ 2 + 4T 2 . (C.6) FIG. 6. Input values for the sphaleron decay rate in eq. (48). Ntr and Nrot are taken from [64, Fig. 5], ω− is taken from [64, Fig. 6], and κ is taken from [67], including uncertainties as explained in the text, and the energy prefactor B is computed also as discussed in the text. All quantities are plot against λ/g 2 , the ratio of the corresponding quartic to the gauge coupling.
potential calculation, we use the values of λ/g 2 as shown in Figure 6, where g = g(µ R ) is the SU(2) gauge coupling, and we consider λ = λ H (µ R ) when we are in the phase P H , and λ = λ Φ (µ R ) when we are in the phases P Φ or P HΦ . This is justified by the fact that the sphaleron solution depends on the SU (2) structure of the theory, and our model mostly has either the SM Higgs or the inert Higgs taking a vev. The sphaleron energy is then given by E sph (T ) = E sph (T = 0) vEW(T ) vEW(T =0) = 4π g B v EW (T ), where the energy prefactor B can be obtained by performing the volume integral of the stress-energy tensor using the previously obtained profile functions. Our choice of B as a function of λ/g 2 is shown in Figure 6, which is consistent with [44].
N tr and N rot are the normalization of the zero-frequency translation and rotation modes [64]. They can be computed from small fluctuations around the sphaleron solution. The resulting formula depends again on the profile functions and can therefore be either computed numerically or read off from [64,Fig. 5]. We pursue the latter and show the values we use in Figure 6.
ω − is the frequency of the unstable mode [64,98]. It can be found as a negative eigenvalue of a system of equations that also depends on the profile functions. We use the values of [64, Fig. 6] directly and show them in Figure 6 (Note that this plot shows ω 2 − in units of (gv) 2 ). κ is the fluctuation determinant. A first numerical evaluation was given in [65], and later improved in [67,[99][100][101]. We use the values given in [67] and assume a rather large uncertainty of [0.01κ, 100κ] to also partially parametrize uncertainties in the other prefactors [97].
Finally, V rot = 8π 2 is the volume of the rotation group; α w = g 2 /4π 2 is the weak coupling constant; and α 3 = α w /(gξ(T )) is the weak coupling in the three-dimensional high-temperature effective theory.
Appendix F: Details on the Numerical Implementation Here we discuss our numerical implementation of the effective potential discussed and thermal calculation in section III. We wrote the main code in python, which is available at https://gitlab.com/claudius-krause/ew nr. We have a second, independent implementation of the code using Mathematica [83], which we extensively cross-checked against the python code.
When including the Coleman-Weinberg potential, we shift the numerical values of µ 2 H and λ H so that the full potential satisfies at T = 0. These values are then used throughout the computation, including inside the RGEs. The finite temperature potential (defined in eq. (10)) can either be evaluated numerically for each point y = M 2 /T 2 (which is slow), or a pre-computed look-up table and subsequent spline interpolation can be used. In the benchmark points discussed in the main text, we use a modified version of the spline implementation of CosmoTransitions [88]. Compared to the original implementation, we extended the pre-computed grid of exact evaluations of the J B/F (y)functions to include more points in the negative y direction and re-wrote the exact evaluation of J B/F (y) to reduce numerical noise. The corresponding files are also included in the GitLab repository.
Daisy corrections beyond the high-T approximation require the second derivative of the thermal potential. We use a numerical, finite-difference derivative based on 9 points chosen symmetrically around the desired functional argument, with a stepsize that increases with large field values or temperatures. We checked that this choice gives a stable value for the derivative for various temperatures and field configurations. Automatic differentiation (AD), as nowadays widely used in the machine learning community [102], would greatly improve the computation of the derivatives. However, implementation of AD in the computation of the effective potential as we do it here would be beyond the scope of this work. In practice, we include the following Daisy approximations in the truncated full dressing scheme of [51]: vanishing thermal masses, i.e., no Daisy correction; leading-order thermal masses in the high-T expansion, i.e. the formulas given in eqs. (C.1)-(C.3); field-dependent thermal masses in the definition of eq. (16); and the thermal masses as defined by the gap equation (15). Note that in the latter two approaches we do not include the Coleman-Weinberg contribution, to properly have the limit Π i (T = 0) = 0.
To reduce the dimension of field space that we have to scan, we assume that all χ i acquire a vev simultaneously. This is justified as long as λ χ = 0 because then the χ-sector exhibits an additional SO(N ) symmetry that allows us to rotate them freely into each other. Because of this enhanced symmetry, the condition λ χ = 0 is also conserved under the RGE.
We use the renormalization scale that we discussed in eq. (19), which is given by the largest square root of the absolute values of the eigenvalues of the bosonic mass matrix including thermal masses in the high-temperature approximation (eqs. (C.1)-(C.3)), or the EW scale v 0 = 246 GeV, whichever is larger. Since the mass matrix at a given point in field space itself also depends on the renormalization scale via the couplings, we have to solve eq. (19) numerically. We use Brent's method [103], as implemented in SciPy [85] for this purpose. Note that we do not include the wavefunction renormalization factors of eq. (D.16) at this point, as this would be numerically more complicated. Instead, we compute these factors at the end, after the minimization of the potential, and rescale the minima positions accordingly. We checked that even at large scales around 100 TeV the factors are at most around 1.05 for the Higgs and at least 0.92 for the inert scalar, so the feedback effect we neglect is in fact small.
A given potential is then numerically minimized. To ensure that we found the global minimum instead of a local one, we use 8 different initial guesses in field space, which cover all possible directions in the three-dimensional space spanned by the Higgs, the inert, and the singlets. As field value, we choose 1.5 (2.5) times the current temperature