Ultra-relativistic bubbles from the simplest Higgs portal and their cosmological consequences

We analyze the phase transitions in the minimal extension of the SM with a real singlet scalar field. The novelty of our study is that we identify and analyze in details the region of parameter space where the first order phase transition can occur and in particular when the bubbles with true vacuum can reach relativistic velocities. This region is interesting since it can lead to the new recently discussed baryogenesis and Dark Matter production mechanisms. We fully analyze different models for the production of Dark Matter and baryogenesis as well as the possibilities of discovery at the current and future experiments.


Introduction
The origin of Dark Matter (DM) and matter antimatter asymmetry are one of the most important unresolved puzzles of the early Universe cosmology. In this paper we will investigate in detail the novel mechanisms proposed in [1][2][3] for DM production and baryon asymmetry generation during the first order phase transitions (FOPT).
The FOPT are known to be very useful for constructing baryogenesis models since the process is out-of-equilibrium as required by one of the Sakharov's conditions [4]. Applications of this mechanism to the electroweak phase transition lead to the seminal scenarios of electroweak baryogenesis [5,6]. In the Standard Model (SM), electroweak phase transition is expected to be a smooth crossover [7,8]. However, this is not the case for various Beyond Standard Model (BSM) frameworks where successful scenarios of electroweak baryogenesis can be built (see for a recent review [9]). The connection between DM and phase transitions is less direct. However, there are class of models where these two phenomenon are strongly related (see for example Ref [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]).
The focus of this paper will be electroweak FOPT with very relativistic bubbles. In this context, collision between bubbles and particles in the plasma, as was shown recently [29], can lead to the production of new states with mass scale significantly larger than the scale of the phase transition 1 . These heavy states can serve as DM candidates [2] or their production and subsequent decay [1,3], if accompanied with C, CP and baryon number violating interactions can lead to baryon asymmetry generation. The successful realization of this scenario require bubble expansion with ultra-relativistic velocities. Such bubble wall motion is known to be possible and is expected within certain classes of the potentials ( [31][32][33][34][35]). However in the case of electroweak phase transition, the requirement of ultra relativistic velocities leads to very non-trivial constraints on the effective potential. In this paper we analyze in detail the simplest extension of the SM that can lead to a FOPT: real singlet scalar field. There have been numerous studies of phase transition in this type of scenario [36][37][38][39][40][41][42][43][44][45][46]. The Lagrangian can be restricted further by considering Z 2 [47] or Z 3 symmetries [48,49]. Interestingly such a real singlet scalar field can appear in composite Higgs models [50] (see also [51,52] for studies pertaining to phase transition). However, most of the above mentioned works did not analyze in detail the region of parameter space with relativistic bubbles since slow bubble wall velocities are required [5,6] for the usual electroweak baryogenesis. We fill this gap by analyzing in detail phase transitions in the real singlet extension of SM with an approximate Z 2 discrete symmetry, focusing on the parameter space which results in very fast bubble expansions. Anticipating the results of this paper we found that the most promising scenario for relativistic bubbles is the case where the phase transition occurs in two steps: the first transition is related to the Z 2 symmetry breaking while keeping EW symmetry unbroken and in the second one the Higgs field acquires its vacuum expectation value (vev) which comes together with Z 2 symmetry restoration.
This will permit us to evaluate the region of parameter space where Baryogenesis via relativistic bubble walls is viable and to put constraints on the mass scale of the dark sector responsible for CP and B violation. Implications for non-thermal DM production from bubble wall plasma particle collisions will also be considered. The article is organised as follows: in section 2, we review the singlet extension of the SM and write the two field potential with thermal corrections. In section 3, we remind the reader of the basics of FOPTs in early Universe and present the computation of the terminal velocity of the bubble wall in the ultra-relativistic regime, for the generic and for SM. In section 4, we qualitatively discuss the different patterns of phase transitions and argue that only a two-steps PT can yield ultra relativistic bubble wall motion. In section 5 we present a numerical study of the two-step PT analyzing in detail the region with fast bubble motion. In section 6, we draw the consequences of our previous results for the production of heavy DM and Baryogenesis and in section 7, we comment on the GW signal induced by such strong transitions. Finally, in section 8, we summarize and conclude.
2 Review of the singlet extension of the SM Let us start by reviewing the effective potential of the SM with a real scalar field (s). The scalar potential will be given by where H is the SM Higgs doublet and m h ≈ 125 GeV is the physical mass of the Higgs. For simplicity, we will impose Z 2 symmetry on the potential to avoid the terms with odd powers of s field. As a result, when s = 0, we avoid any scalar mixing terms which are constrained by the recent Higgs signal strength measurements [53]. The Higgs doublet can be decomposed as usual where h is the usual Higgs boson getting a vev which is given by v EW = m 2 h /2λ ≈ 246 GeV if s ≡ v s = 0.

Coleman-Weinberg potential
Next, we take into account 1-loop corrections which are encapsulated by the Coleman-Weinberg (CW) potential [54] V CW = i=Z,h,W,t In this expression, M i stand for the masses depending on the Higgs and singlet fields values, M i ≡ M i (h, s), and M i0 are the field values from the tree-level vev, M i0 ≡ M i (v EW , 0). Then one can easily check that this potential corresponds to the following renormalization conditions Eq.
(3) cannot be used directly for contributions of the Goldstone bosons, since their masses vanish in the true vacuum and the CW potential is IR divergent. However the solution to this issue was suggested in Ref. [55] which emphasized that the physical Higgs mass is defined at p 2 = m 2 h and the effective potential at p 2 = 0 (see for details Ref. [55]). Thus it would be better to use modified renormalization condition where the differences of self-energies are taking into account the running of selfenergy from p 2 = 0 to p 2 = m 2 h . In this case the IR divergences in Σ(0) due to the virtual Goldstone bosons are cancelled against the IR divergences of V eff . In practice, at the end of this renormalization procedure, the contribution of the Goldstone bosons is given by : The number of d.o.f, masses of various particles and the scalar mass matrix as a functions of h, s are given by: n W ± = 6, n Z = 3, n G ±,0 = 3, n t = −12 , As a side remark, in the region where the Higgs h → 0 and s ∼ O(v EW ), there will be two scales involved in the problem: the value of the singlet field s and the masses of the SM particles M W,Z,t (h → 0, s) → 0 s. This type of two scale potential has been studied in the past [56], by using two different renormalisation scales. It was concluded that resummation is needed when the log (M i (0, s)/s) is large enough to cancel the loop suppression. Although we have two largely separated scales, we have checked that for our region of the parameter space, such a resummation is not necessary.

Finite temperature potential
The temperature and the density effects can be taken into account by complementing the zero temperature potential with thermal corrections (see for example [57,58]), In Eq.(8), V T =0 eff (M i ) is the potential we computed in the previous subsection and the thermal potential V T (M i ) is given by However, to save computation time, we use the approximate expansion of the function J B/F (y 2 ) as given in Ref. [57]: where γ ≈ 0.577 is the Euler constant, and K 2 (z) are the second-kind Bessel function. To account for dangerously divergent higher loops due to the Daisy diagrams at finite temperature, we follow the so-called "Truncated-Full-Dressing" procedure [57]. Doing so, the full one-loop potential becomes where Π i (T ) are the thermal masses of various degrees of freedom. In the real singlet extension of the SM [57], the expressions of the thermal masses read Scalar: Gauge: where Π L g (T ) denote the thermal mass of the longitudinal mode of the gauge bosons, while transverse modes Π T g (T ) are protected by gauge invariance and thus do not receive a mass at leading order in perturbation theory. 5

First order phase transitions with relativistic bubbles
Starting with the effective potential derived in the previous section, we can proceed to the analysis of the phase transition. FOPT happens when the minima of the potential corresponding to different phases are separated by a potential barrier and the transition happens via bubble nucleation. The probability of nucleation of a bubble is given by [59][60][61] where S 3 , S 4 are O(3, 4) bounce actions and R 0 is the initial bubble radius. Bubble nucleation is characterised by a critical temperature T crit , defined as the point when the two phases of the system have vacua with the same energy. Below T crit phase transition becomes energetically possible. The probability to find a specific point of the Universe to be in the false vacuum is given by [34,62,63]: In Eq.(14) the strongest dependence on the temperature comes from the Γ(T ) ∝ exp (−S 3 /T ), so that the quantity I(T ) is mostly controlled by the ratio Γ (T ) /H 4 (T ) and an order one fraction of the volume of the Universe will be in the true vacuum when Γ[T ] ∼ H 4 [T ]. The temperature that satisfies this condition is coined as the nucleation temperature T nuc or, to phrase differently, the nucleation temperature T nuc corresponds to the appearance of roughly one bubble of true vacuum per Hubble volume. Assuming the scale of the phase transition to be O(100) GeV we get In case we want to be more precise about the temperature when the phase transition completes, we can define the so-called percolation temperature T per as the temperature when around ∼30 % of the space has been converted to the true phase If the condition in Eq. (18) is not fulfilled, I(T ) < 0.34, then the bubbles do not percolate. At last the energy released during the PT is traditionally quantified by the ratio between the energy stored in the vacuum and in the plasma at the moment of the transition 6

Velocity of the EW bubble
Let us proceed to the computation of the bubble wall velocity v w . The dynamics of the bubble wall motion is controlled by the driving force due to the potential differences between the false and true vacuum and the friction due to finite temperature effects. The calculation of the friction is generically a very complicated problem, however for the ultra-relativistic bubble motion at leading order (LO -tree level) very simple expressions have been obtained [64,65] for the pressure force from friction, where ∆M 2 i is the change in the mass of the particle i during the PT, c i = 1 (1/2) for bosons (fermions) and g i is the number of d.o.f of the incoming particle. Eq. (21) assumes that the masses of the particles outside of the bubble are less than the temperature, otherwise the friction will have an additional Boltzmann suppression ∝ exp[−M false /T ]. Interestingly the production of the heavy particles can also contribute to the friction at the same order [29] ∆P mixing where v is the vev of the Higgs field and we review the heavy particle production later in the section 6.1. One can see that the friction (pressure from plasma on the bubble wall) becomes velocity independent so that permanent accelerating (runaway) behavior of the bubble expansion becomes possible. However for theories where the gauge bosons receive a mass during the phase transition, this is not the case and the effect with multiple gauge boson emissions, leads to the additional contribution (NLO-Next to Leading order) to the pressure which scales as [20,29,66,67] At this point we can see that bubbles will keep accelerating till the moment when NLO friction becomes large enough to balance the driving force which will set the γ terminal w the bubble can reach, Note that we can be as well in the situation where the percolation (bubble collision) starts before the terminal γ terminal w is reached. 2 There is a claim that this friction scales as γ 2 T 4 [68] see [29,67] for criticism.

Friction forces during the electroweak phase transitions
Let us apply the discussion of the previous section to the EW phase transition (EWPT). Our main interest will be the possibility of heavy particle production which can later source baryogenesis and DM production following the ideas in [1][2][3]29]. The heavy particle production (see also the discussion in section 6.1) happens due to the collision between the plasma particles and the bubble wall. The typical center of mass energy for such process will be roughly ∼ √ γ w T v EW , thus γ w will be controlling the maximal mass of the new fields that can be produced. To calculate γ w in the context of the EWPT, we need to know the forces acting on the bubble wall.
The LO friction from the top, Z, and W takes the form: There is also a contribution to the LO friction from the singlet and Higgs scalars, however it depends on the masses of these fields in the false vacuum and we find it to be numerically subleading compared to the estimate in Eq.(26) 3 . In our numerical calculation we took this additional contribution into account, however, for the current discussion it is sufficient to use Eq. (26). This gives a rough condition on the nucleation temperature for the transition to become ultra-relativistic The computation of the NLO friction in the SM has been carried out in [67] and the following approximate expression has been derived: where ν a = 1(3/4) for a a boson (fermion), g a is the number of degrees of freedom of a and M g,i (v EW ) is the mass of the gauge boson inside the bubble. C abc represents the coupling which appears in the vertex (see appendix C), β c=Z 0 = 1, and β c=W ± = cos θ W = M W /M Z , α = e 2 (v EW )/4π ∼ 1/128 is the electromagnetic fine structure constant. The renormalization scale µ has to be understood as the lower cut-off for the integration over soft momentum, typically the thermal mass µ ∝ α 1/2 i T nuc [29]. The κ factor is introduced [67] to account for the contributions of the both reflected and transmitted bosons and is approximately equal to κ ∼ 4. The sum in the square brackets is approximately equal to abc ν a g a β c C abc ≈ 157 , Taking into account that (∆V − ∆P LO ) /(100 GeV) 4 O(1), we can see that the bubbles will become ultra-relativistic, i.e. γ terminal 1 only if T nuc is significantly lower than the scale of the phase transition ∼ 100 GeV.

Phase transition in the singlet extension
After the preparatory discussion in the previous sections 2 and 3, we can proceed to the analysis of phase transition in the model with the singlet field, Eq.(1). Our analysis will be focused on the region of parameter space with relativistic bubble expansion, for other studies of phase transition in SM plus Z 2 real singlet scalar, see Refs. [47,[69][70][71].
In the previous section 3, we have seen that the velocity of the bubble expansion is fixed by the balance between the friction from the plasma and the driving force. At low temperatures the friction is suppressed (see Eq.(26)- (28)) so that we expect the bubbles to become relativistic (large Lorentz γ terminal w factor). Let us check whether low nucleation temperatures are feasible in the singlet extension. We will assume that Z 2 remains unbroken in the true vacuum in order to avoid to constraints from the Higgs-scalar mixing (see for example [72]). Then in the model with Z 2 odd singlet, the phase transition can occur in two ways: one-step ( h = 0, s = 0) → (v EW , 0) and two-steps 4 [44,47,70], and each of these phase transitions can be first or the second order. We review both of these scenarios of phase transitions in order to understand in what case it is possible to obtain relativistic bubbles.

One-step phase transition:
This case has been largely studied in the literature [47,70] and we will not provide a complete description of it. In this scenario the singlet never gets a vev, and all of its effect reduce to the additional contributions to the Higgs potential from Coleman-Weinberg terms and thermal corrections. However it turns out that relativistic bubbles are very unlikely for such phase transitions (see also results in [47]). In the limit when |m s | T c we can show analytically that this is indeed the case. Near the origin h → 0, the potential in h direction is dominated by the ∝ h 2 terms The effective mass m 2 eff include the tree-level terms and the leading thermal contributions and is approximately equal to The FOPT can happen only if m 2 eff (T ) > 0. Then assuming perturbative values of the coupling λ hs we can estimate the lowest temperature where the FOPT can occur to be T nuc min 100 GeV. Comparing this value with the discussion in the section 3.2 we can see that the bubble wall velocities will always satisfy γ w 10. As mentioned before, we are interested in the expansions with much larger γ w factors, so that we do not discuss one step phase transition further.

Two step FOPT with relativistic bubbles
Two-steps realisations of the EWPT have already been studied in many works, see for example [44,47,70,[73][74][75][76][77]. The novelty of our study is that we will be focusing on the parameter space with relativistic bubbles which was previously ignored. We organize the discussion as follows: In section 4.2, we show qualitative results based on approximate treatment of the potential and then in the section 5 we present the exact numerical results obtained with our code. The two step phase transition can happen if the m 2 s parameter of Eq.(1) is positive. In this case it is convenient to parameterize the Lagrangian in the following way where v EW = m 2 h /2λ GeV and v s = m 2 s /2λ s correspond to the local minima at ( h = v EW , s = 0) and ( h = 0, s = v s ) respectively. The origin of two-step PT can be intuitively understood from the following considerations. For simplicity let us ignore the Coleman-Weinberg potential and restrict the discussion by considering only the thermal masses. Then the potential will be given by From this expression we can clearly see that the temperatures when the minima with non zero vevs appear for the Higgs and singlet fields can be different. Then it can happen that the Z 2 breaking phase transition occurs before the EW one. This means that there will be first a phase transition from (0, 0) → (0, v s ). After this phase transition the Universe keeps cooling down and the minimum with h = 0 will be generated. Choosing the appropriate values of masses and couplings we can make sure that the minimum with (v EW , 0) is the true minimum of the system. The transition (0, v s ) → (v EW , 0) will be of the first order if there is a potential barrier in the between two minima. One of the necessary condition in this case will be ∂ 2 V /∂h 2 | s=vs,h→0 > 0, which using Eq.(36) we get From this discussion we can see that there are qualitatively two different cases depending on whether the potential barrier between two minima remains or disappears at zero temperature, i.e. when m 2 h ≷ λ hs v 2 s . In the first case the phase transition will necessarily happen before the "No Barrier"(NB) temperature since the bounce action drops once the potential barrier between two minima reduces.
In the other case, when the barrier remains even for zero temperature potential, the phase transition is obviously of the first order. However it might happen that the tunneling rate is too slow and the system remains stuck in the false vacuum.
From Eq. (37) we can see that size of the potential barrier between two minima will be controlled by the coupling λ hs . Increasing this parameter will enhance the potential barrier and will lead to the reduction of T NB so that (v s , 0) remains a local minimum even at zero temperature. At some point we expect the potential barrier to become so large that the system will remain trapped in the false vacuum forever. From this discussion we can expect that the lowest temperatures will be achieved at the boundary of the region where no PT can occur. The lowest temperature will correspond to the (global) minimum of S 3 /T controlling the tunneling rate, as explained later in section 5. If this minimum happens at temperatures much lower than T crit , then there will be super-cooling (T nuc T crit ). Even though this discussion was made by considering only the thermal masses in the potential, numerically we find that for the effective potential with truncated full dressing [57] the qualitative behaviour does not change, and only explicit conditions for m 2 h and T nuc min are modified.
At last we always check whether the condition for the EWSB minimum to be the global minimum at zero temperature is satisfied: where above equation is valid up to small loop-level corrections.

Numerical results
After the qualitative discussion, let us proceed to the numerical calculations. The bounce action for (0, v s ) → (v EW , 0) was computed using our own dedicated code (cross checked against FindBounce [78]) and we relegate the details and methods of this calculation to Appendix A. We parameterize our model in terms of (m s , λ hs , v s ) parameters (see Eq. In the plane (see Fig.1) v s − λ hs , we identify four regions with different behaviours under the phase transition. The blue region shows second order transition when there is no barrier between the two separated vacuum. Next to it there is a light red region where the transition is of the first order. We indicate separately (dark red) the region, where the transition is of the first order and the bubbles are relativistic. In particular the boundary between the regions with relativistic and non-relativistic bubbles is defined by the criteria of Eq. (27), i.e. when the LO pressure for relativistically expanding bubbles is less than the driving force. At last, there is NO PT region, where the system remains stuck in the false vacuum since the tunnelling rate is too small. The structure of the diagram (on the Fig.1 top−left panel) can be easily understood from the qualitative discussion in the previous section. Indeed, keeping v s fixed, the size the potential barrier is controlled by the coupling λ hs . Moving from left to right the size of the potential barrier increases and we are gradually moving from the region of second order phase transition → FOPT → FOPT with relativistic bubbles → no PT region. Similarly moving up (increasing v s for fixed values of λ hs ) also corresponds to the increase of the potential barrier as we pass though the regions with different PT in the same order. On the upper axis we report physical mass of the singlet in the true vacuum M s (v EW , 0) and exclude the constrained region where h → ss is kinematically allowed (the gray meshed region).
Next we plot the values of T nuc and γ terminal w as a function of (λ hs , v s ) ( Fig.1). As discussed in the previous section, the region with the smallest values of the nucleation temperature (thus the fastest bubbles) is located near the "NO PT" region, i.e.
where the system remains trapped in the false vacuum. The blue dot on the T nuc plot indicates the last point we have found before the system enters the regime of no phase transition (NO PT). In the bottom right we plot M max = √ γ w T nuc v EW quantity which indicates the maximal energy in the c.o.m frame for the plasma particle−bubble wall collision, which corresponds to the largest mass of the heavy particles we can produce (see discussion in sec 6.1).
In order to better understand the dependencies of T nuc and γ w on the parameters of the model, it is useful to separate further the parameter space depending on whether the potential barrier disappears at zero temperature or not. This is important since the bounce action in such cases has very different dependencies on the temperature. See for example Fig.2, where we have plotted the S 3 /T for v s = 170, m s = 125 GeV for the various couplings λ hs . In the case when the potential barrier disappears at some temperature T NB the function S 3 /T drops to zero, but if the barrier remains even at zero temperature S 3 /T has a global minimum for T = 0 which will be controlling the lowest nucleation temperature possible.

No potential barrier at zero temperature
Let us start by defining the region where there is no potential barrier at zero temperature. On Fig.1, we demonstrate the curves where the barrier disappears for various values of the temperatures (T NB ). For T NB = 0 case, approximate curve can be obtained analytically by looking at the leading terms in the zero temperature CW potential  The agreement between this equation and exact T NB = 0 curve is at the level of a few permille discrepancies. To the right of T NB = 0 curve, the potential barrier between the two minima remains even at zero temperatures. For the values v s 200 GeV, we find that the line T NB = 0 approximately coincides with the boundary of no phase transition region (where system remains stuck in the false vacuum) but obviously the boundary of "NO PT" is always to the right of T NB = 0 curve. The size of this narrow strip is of the order 10 −4 in λ hs values. One can see it from the T nuc panel of Fig.1 where we have indicated the value of λ hs when T NB = 0 by vertical thin line and red dot (for intersection) and the position of the blue dot which is the last point where the transition is of the first order before we enter NO PT region. The boundaries of this region were obtained by numerical calculations where we have scanned λ hs parameter with a step 10 −6 . We postpone the discussion of the FOPT in this narrow region to the next section 5.2.
In this section we restrict our discussion only on the region to the left of T NB = 0 curve. Then the phase transition will be always completed before the universe cools down to T NB , i.e. T NB < T nuc , which provides a lower bound for the nucleation temperature. At the same time the velocity of the bubbles become largest for the smallest possible values for the nucleation temperature. So that the fastest bubbles will be near T NB = 0 curve. Looking at Fig.1 we can see that the largest γ w (Lorentz boost factor) and lowest nucleation temperatures happen for v s 200 GeV, where the T NB = 0 curve passes very close to the NO PT boundary. The shape of the lines in Fig.1 clearly indicate the necessity of tuning in order to obtain low nucleation temperatures (large γ w ). In particular for the values of v s 200 GeV we can see that the nucleation temperature drops by choosing λ hs close to the NB value (similarly γ w becomes maximal see Fig.1). We can estimate this tuning by looking at ∂ log λ hs /∂ log T nuc quantity (analogue of Giudice-Barbieri [79] measure of the tuning) as a function of T nuc . This result agrees with our expectation from the steepness of the curves in the Fig.3 and with the naive tuning expectation which scales as ∼ T nuc /m h 2 5 . The dashed green line represent the naive tuning ∼ (T nuc /m h ) 2 . We observe that this naive estimation for the tuning is rather precise at large nucleation temperature but can underestimate the tuning by one order of magnitude for very low nucleation temperature.
At last we would like to remind that the discussion in this section always assumed that the phase transition completes before the potential barrier disappears. We have checked numerically that this is always the case. Indeed the time of the phase transition is given approximately by the bubble radius at the moment of percolation [80,81] This radius is related to the β ≡ − d dt T parameter by an approximate relation [81] where we find R −1 ∼ β typical ∼ (10 − 10 4 )H. At this point the temperature drops during the bubble expansion will scale as Due to the large value of β/H we find numerically that this drop of the temperature is not enough for the barrier to disappear or in other words Such behaviour can be understood from the following consideration near T NB the bounce action drops very quickly and so that the tunneling becomes very efficient of temperature.
almost instantaneous and typical bubble radiuses are much smaller than the Hubble scale. This leads to another prediction that GW signal will be suppressed as well since it is controlled by the (β/H) quantity Eq. (42). As we will see, even with this suppression the GW signal is efficient enough to be detected in the future.

Tunneling with potential barrier at zero temperature
Let us proceed to the analysis of the case when the potential barrier does not disappear at zero temperature. The parameter space with the lowest nucleation temperatures (fastest bubbles) will be located again near the "NO PT" boundary. However in this case the nucleation temperature will be controlled by the local minima of the S 3 /T function. (see Fig.2). At least a minimum is expected since the potential at low temperature becomes fully temperature independent and S 3 (T → 0) → const., so that S 3 /T necessarily starts to grow for T → 0. Numerically (see Fig.1) for the value of m s = 125 GeV we find that for v s 180 GeV entire region with the fast bubbles has a potential barrier at zero temperature. On Fig.2, we present the euclidean action for v s = 170 GeV and m s = 125 GeV. Going back to Fig.1, we see that for those values the "NO PT" curve and the "T NB = 0" curve are largely separated. This is not a surprise since in this region of parameter space the bounce action S 3 /T ∼ O(10 2 ) is small enough to guarantee the successful tunneling even when the barrier remains at zero temperature. In the range of λ hs from 0.544 to 0.58, numerically we find that nucleation temperatures are between 65 − 44 GeV, with a clear saturation at T sat nuc ≈ 44 GeV, and the corresponding Lorentz factor for the velocities of the bubble expansion in the ranges of 10. Interestingly we find that in this case the bubble radius R ∼ (8π) 1/3 vw β are a little bit larger than the ones discussed in the section 5.1, corresponding to a bit smaller values of β/H parameter.
Super fine-tuned region We finally comment on the parameter space with v s 200 GeV (again we are fixing m s = 125 GeV), where the curves "NO PT" and "T NB = 0" almost superimpose (the region between red and blue dots on the T nuc panel of the Fig.1). There will be a very narrow strip between the curves "NO PT" and "T NB = 0" regions, where the tunnelling will happen even though the barrier remains at zero temperature. We find (see Fig.5) that the region corresponds to the variations of the λ hs parameter of the order δλ hs ∼ O(10 −4 ), i.e. two order of magnitude smaller than the full region with relativistic bubbles. In this very small region various additional effects can start playing a role. For example let us look at the Fig.5 we can see that the bounce action S 3 /T has a local maximum and a deeper (global) minimum with respect to the standard scenario. Such a behaviour of the action is coming from the cancellations of various terms in the effective potential. For simplicity let us look at T = 0 case. Then there is a region of parameter space where purely polynomial potential has no local minimum at (0, v s ), but the effects of the Mt(v EW ) terms in Coleman-Weinberg contribution lead to the appearance of the local minimum at (δv h , v s ). On Fig.4 we plot the contributions of the various terms in the effective potential leading to the appearance of this local minimum and the trajectory of the typical bounce solution in this case. As a result the distance in the fields space between the two minima decreases and the tunneling becomes faster, which leads to the appearance of the second (global) minimum in S 3 /T .

Benchmark points
In Table 1 and 2, we give typical values of the nucleation temperature, the Lorentz factor γ w , the β/H factor, and we indicate if the barrier remains at zero temperature, applying the criterion in Eq. (44). We can see that the largest bubble radius at the collision (smallest β) correspond to the case when S 3 /T is monotonic and very flat near the tunneling temperature (c.f. the right panel of Fig. 5).
6 Consequences for production of dark matter and Baryogenesis One of the motivation for the study of a model of EWPT with relativistic bubbles is the relation between relativistic expansion and the out-of-equilibrium production of heavy states presented for the first time in [29], when the field undergoing the PT (here the Higgs) is coupled to some heavy dark sector at typical mass M N . In   Table 2: Same as Table 1, but for Fig.5 and with v s = 205 GeV. We observe that the last two points display a displacement of the false minimum.
this section, we remind the principle of the production mechanism and we study the scenario of the production of Dark Matter [2] and Baryogenesis [1], that were previously agnostic about the EWPT realisation. First of all, the strong FOPT involves a supercooling represented by a dilution factor, with g s (T ) (g sym s (T )) being the number of relativistic degrees of freedom of the entropy in the broken (symmetric) phase. This means that with D 1, any type of dark matter production or Baryogenesis mechanism that happens much earlier than the PT should provide values denser than the conventional estimation by a factor of 1/D (see for example [19,82]). For instance, the WIMP cross section should be σ ∼ D 10 −3 10 −29 cm 3 /s to produce a correct dark matter abundance. This is the case when freeze-out happens at temperatures much higher than the reheating.

Production of heavy states during ultra-relativistic expansion
There are few mechanisms which can lead to heavy particle production during FOPT. This can happen if the incoming massless particle in the unbroken phase gets a very large mass from the Higgs vev [3,64] (mass gain), or due to the bubblebubble collision [25,30] or due to the plasma particle−bubble collision [29]. Our study will be focused on the later one. Let us assume that FOPT happens and the bubble expansion is indeed relativistic with γ terminal w 1. The simplest model where the production of heavy particles during plasma−bubble wall collision can be realized, is described by the following Lagrangian [29]: where q is a massless particle in the symmetric phase and N is a heavy field with large vev−independent mass M N v EW . h is the Higgs field undergoing a FOPT. With no loss of generality, we go to the basis where fermion masses are real. Before the strong phase transition starts, the abundance of heavy states N in the plasma is strongly Boltzmann suppressed and they would naively seem irrelevant for the dynamics of the transition. In an homogeneous vacuum, the transition from light to heavy state q → N is obviously forbidden by the conservation of momentum. However, in the presence of the bubble wall, the conservation of momentum along the z direction is broken (assuming a planar wall expanding in the x − y plane) and a computation using WKB phases for the q and N fields demonstrates that the probability P(q → N ) is non-vanishing [29] and is given by with L w ∼ 1/v EW the width of the wall. Behind the bubble wall a large abundance of N andN , n BE N is produced. Let us emphasize that this abundance is much larger than its equilibrium value.
Another possibility of the heavy particle production can be realized for the following interaction In this case φ is a heavy scalar field with mass M φ v EW , then in the vicinity of the wall, the process h → φφ has the probability [29] The results in Sec. 3.1 on the terminal velocity in the singlet extension of SM allow us to compute the maximal mass of the particles which can be produced during the electroweak FOPT in the singlet extension. Indeed saturating the step function in the above equation and assuming the L w ∼ 1/v EW we get approximately: Numerical results for the maximal mass M MAX are reported in Fig.1. We can see that the maximal mass we can produce is roughly ∼ 10 TeV scale. We would like to note that our results can be easily applied for the mass gain mechanism of the heavy state production [3]. Indeed in this case the maximal mass will be M mass gain γ w T , and can be read off from the bottom right plot of the Fig.1 by noting that it will scale as M mass gain ∼ M 2 MAX /v EW . Since the mass of the heavy field comes from the vev of the Higgs, it will additionally be bounded by the unitarity considerations to be below 2 TeV.

Dark Matter production
In this section we will apply the results for the velocity of the bubble expansion for DM model building.

Scalar DM coupled to the Higgs portal
We assume a heavy scalar φ coupled to the SM via the traditional Higgs portal The DM (φ) field is stabilized by some additional Z φ 2 (we use this subscript to differentiate it from Z 2 of the singlet potential). After the Higgs transition, the abundance of massive φ, n BE φ , behind the wall is given by We can see that DM production during the bubble expansion is strongly dependent on the density of the Higgs field available at the nucleation temperature f h (p, T nuc ). The relevant parameter for the discussion is the ratio where fv denotes the position of the false vacuum and d 2 V /dh 2 fv is the mass of the Higgs m False H in the false vacuum. As soon as this quantity becomes larger than 1, we expect exponential suppression of the Higgs abundance Here we define C eff to take into account the Boltzmann suppression. After redshifting to today, the stable produced abundance takes the form This expression has to be supplemented with the freeze-out(FO) contribution which is produced before the phase transition Note that the FO contribution is suppressed by the factor T nuc /T reh 3 due the brief stage of inflation during the phase transition. Obviously the prediction for relic density must match the experimental observations: Ω today φ,tot h 2 ≈ 0.1. We can see  Table 3. from Eqs. (55),(56) that for small values of the portal coupling λ hφ , DM production will be dominated by the freeze out mechanism while bubble expansion takes over for larger values of λ hφ .
Next we can check whether this mechanism for DM production can lead to viable phenomenology, given the results on bubble dynamics in section 5. Instead of making a scan of the parameter space, we will just focus on a few representative benchmark points.
For m s = 125, v s = 205 GeV, we show in Fig.6 the isocontours reproducing the correct relic density for three reference values of λ hs (corresponding nucleation temperatures can be found in the Table 1). Firstly, C eff 1 for all three reference points. For λ hs = 0.424 the upper red curve corresponds to the case when DM production is dominated by the BE (bubble expansion) and the lower curve by FO. The steepness of the upper red curve (BE) comes from the fact that we are always in the region of parameter space where exp −M 2 φ /(2γ w v EW T nuc ) 1, leading to a very strong sensitivity on M φ mass. Physically this means that the model generically predicts large overproduction of DM in BE process unless the Boltzmann suppression exp −M 2 φ /(2γ w v EW T nuc ) is playing a role. For the other two reference points λ hs = 0.4242, 0.42424 we can see that there is an additional part of parameter space for the DM masses M φ ∼ 1 − 4 TeV, which corresponds to the region without the Boltzmann suppression exp −M 2 φ /(2γ w v EW T nuc ) ∼ 1. This is related to larger values of M MAX ∼ √ γ w v EW T nuc and smaller values of the nucleation temperature, reducing the excess of the DM abundance. On the right panel of Fig.6, we report similar plots for v s = 175, m s = 150 GeV. Finally, before closing this section, we comment about the possibility of considering the singlet s itself, in the limit of very precise Z 2 , as DM. After the phase transition T ∼ 40 GeV, the singlet is in thermal equilibrium and we can apply straightforwardly the freeze-out expression: From this estimate of the FO abundance for s and recalling that we considered λ hs ∼ 0.3 − 0.6 and M s (v EW , 0) ∼ 100 GeV, we conclude that the abundance of s produced in this fashion, today, is underproduced by one or two orders of magnitude to fit the observed amount of DM Ω today s, FO h 2 ≈ 0.1. Even in this underproduced case, there are severe bounds from the direct detection experiments except for the resonant region, where (57) is over-estimated. However, as we will discuss in the Appendix. D, we will have a Z 2 explicit breaking which makes s decay much before today.

Singlet portal DM
In this section we mention an alternate possibility of coupling the DM (φ) to the singlet field s via "singlet portal" (58) where the width of the singlet wall is similar to the length of the Higgs wall L w . Interestingly even though FOPT is from (0, v s ) → (v EW , 0), the singlet scattering of the wall can lead to the production of the φ field. Phenomenology of DM production is very similar to the Higgs portal case discussed in the previous section, with one main difference: in the false vacuum, the mass of the singlet is not small and the factor C eff introduced in the Eq.(54) plays an important role. The results are shown in Fig.7. For example, if we compare the curves for λ hs = 0.424, m s = 125, v s = 205 in Fig.7 and in Fig.6, we can see that, for the Higgs portal DM, the isocontour has the same shape, but with the larger values of DM M φ masses. As shown in Eq.(55), this is due to the proportionality between the DM relic abundance produced during the bubble expansion with ∝ C eff /M φ .
Singlet portal with additional field A slight modification of this scenario is to further introduce a light scalars. Then we can have, where again φ is the DM and we did not write down the mass terms for simplicity of notation. We assume thats is in the thermal bath before the PT. Then due to the field change of s in the bubble wall, the momentum conservation violating process s → φφ can occur (h → φφ may also occur if there is the h 2 φ 2 term.) In this model s →s + SM particles happen via the DM loop.

Fermion-mediated Dark Matter
In the previous section we noticed that the DM production during bubble expansion strongly depends on the mass of the incoming particle in the symmetric phase, due to the Boltzmann suppression factor. In case of an incoming scalar, generically this effect is relevant and crucially modifies the phenomenology, as we have seen in the section 6.2.2. In this section we construct a model where the incoming particle is a massless fermion in the symmetric phase, so that C eff ≡ 1 by definition. The model consists of a vector-like neutral fermion N which is a singlet under SM and a couple of Z 2 odd fields φ and χ: Here, L, H are SM lepton and Higgs doublets, respectively. The production mechanism works as follows: the heavy field N is produced during the phase transition L → N and it will subsequently decay into N → χφ, N → LH. The field N can be Majorana or Dirac (in the former case there will be a relation to neutrino masses and in the later it will be completely independent from neutrinos). In this model, heavy N are produced via L → N with a probability As a consequence, unstable heavy N accumulate behind the wall with initial density given by where v w = 1 − 1/γ 2 w , we expanded for large γ w and approximated the Fermi-Dirac distribution as a Boltzmann distribution. Compared to the original proposal in Ref. [2], the density of the heavy fields inside the bubble will be additionally enhanced by ∼ 16π 2 factor since 1 → 1 transitions are more effective than 1 → 2. Let us assume that M φ < M χ so that φ is the DM candidate, then DM production will happen via the following chain of processes: However, the heavy N has two channels of decay: toward the heavy dark sector φ, χ and back to the light L. The abundance of φ, χ after the transition is thus suppressed and given by and the final relic abundance redshifted to today thus reads For the freeze-out process in the symmetric phase, we have: φφ → LHLH by neglecting co-annihilation. The cross-section is highly phase space suppressed (closing a loop for a 2 to 2 annihilation gives a similar scaling): The abundance by taking account the supercooling is The total DM density today will be given by the sum of Eq.(65)- (66). Therefore, this scenario leads to the over-production of DM unless M φ , M χ 10 GeV. However, note that these equations are valid only for the heavy DM candidates which do not go back to equilibrium after the phase transition. Otherwise, the final density will be given by Eq.(66) only without T nuc /T reh 3 and we are going back to the normal freeze-out scenario.
Let us now investigate the regime M φ M χ , precisely |M φ − M χ | M φ /20, where the co-annihilation takes place. In this case we have the channel φχ → HL to decrease the abundance of φ as well as χ. The cross-section is Summing this estimate with the Ω today φ,BE h 2 in Eq.(65) we find that it becomes possible to reproduce the observed DM abundance. However we see that bubble expansion tends to overproduce the DM and the relic abundance in BE can be reproduced if only the factor exp[−M 2 N /(2v EW T nuc γ w )] starts playing a role in suppressing DM relic density (left boundary of Fig.8 is almost vertical).

Baryogenesis mechanism
Now we remind the scenario of Baryogenesis with relativistic bubble walls that was proposed in [1]. As a prototype, we worked with the following model, reminiscent of the toy model of Eq.(46), (omitting the kinetic terms) Thus, additionally to the SM and the singlet sector, the model contains a Majorana field χ and two vector-like B quarks with the masses M 1,2 ∼ m χ . η is a scalar field in the fundamental representation of QCD and with electric charge Q(η) = 1/3. We defined Q, u, d as the SM quark doublet and singlets respectively, and ignored the flavour indices. H is the SM Higgs. We first proved that the production mechanism, when computed up to one loop-level, indeed transforms a CP-violating phase into a chiral asymmetry in the abundances produced, in a fashion very similar to usual leptogenesis decay. Let us now sketch the mechanism in itself. First, b−quarks collide with the relativistic wall and produce B I , B c I via the mechanism explained above. Thus inside the bubble we have where n inc b is the number density of the bottom-type quark colliding with the wall from outside, n b and n B are the abundances inside of the bubble and is a loop suppressed coefficient which parametrizes the CP violating phase and the resonance between one-loop and tree-level diagrams. After the passage through the wall, the following asymmetric abundances are then generated From this expression we can read an apparent asymmetry in the bottom quark abundances. However, if the heavy B freshly produced were to decay back into b, the asymmetry would be washed out. This is however not the case if they decay in a dark sector containing χ, η, where the asymmetry is enhanced by the presence of a Majorana mass for χ. The final unsuppressed produced asymmetry is given by where s is the entropy at the moment of the production, g b is number of degrees of freedom of the bottom and g the number of relativistic degrees of freedom. The loop functions f IJ B have been computed in [1] and are controlled by the CP violating sector. Absence of strong wash-out conditions as well as experimental signatures (direct production in colliders, flavor violation, neutron oscillations) pushed the heavy particles to be M B,χ,η 2 × 10 3 GeV .
In the context of singlet extension with Z 2 that we studied, this opens up the range where the Baryogenesis mechanism proposed above is operative.

Impact of the heavy sector on the phase transition
The models we are considering by construction have new heavy fields coupled to the Higgs boson. These will lead to the finite corrections to the scalar parameters of the form (assuming a Yukawa type connection yBHb) where g N is the number of heavy degrees of freedom and M N is the typical mass of the heavy sector. One can wonder how these corrections can effect the tuning of the Higgs potential. However note that in our model the Higgs mass hierarchy problem is not addressed and generically we expect the size of m 2 h to be of the order of the cut off scale (M pl in SM). So the corrections in Eq. (75) do not make the tuning worse.
In case the Higgs hierarchy problem is solved at the scale of the heavy fields in Eq.(51),(68) the tuning in the Higgs potential will be roughly, We can combine this estimate with a tuning for low nucleation temperatures (see discussion in section 5.1) which are necessary for the heavy field production and the tuning estimate becomes: Using the estimates of the maximal values of γ w and the maximal mass of heavy particles which can be produced during the bubble-plasma collisions (see Eq. (31) and Eq. (61)) we get the following estimate for the maximal tuning in the model where we remind the reader that this estimate is valid only if the Higgs hierarchy problem is solved at the heavy fields scale.
The signal produced at the moment of the transition can be separated into different contributions: the bubble collision [98] contribution, the plasma sound waves [84] and finally the turbulence. Only the two first sources of GW are well understood. Another nice feature of those two sources is that they are expected to dominate in different physical situations; the bubble collision would dominate in case of runaway wall and the sound waves if the wall reaches a terminal velocity. We have already mentioned that the EWPT, if first order, will always happen in the regime of terminal velocity, because of the large number of strongly coupled vector bosons 6 . For GW produced by plasma sound wave, the peak frequency and amplitude are given by with z p ∼ 10, κ sw is the efficiency factor for the production of sound waves in the plasma [99], α and R have been defined in Eqs. (19) and (41) respectively, R ∼ O(10 −1 − 10 −4 )H −1 is the approximate size of the bubble at collision and all quantities (T, H, g ) have to be evaluated at reheating.
As we have seen in sections 6.3-6.2, for the baryogenesis and DM production we need relativistic walls with relatively low nucleation temperature 10 GeV. In this context, α 1 and κ sw → 1. The peak frequency and the signal amplitude are only function of the size of the bubbles at collision, which are reported in Table  2 and 3. We can observe that in this range β/H spans the value between [50, 10 4 ], with a preference for lower values. Going back to Eq.(79), emitted amplitude and frequencies will be of the order where we set z p = 10, g = 100. This range of frequencies and amplitude are largely in the expected sensitivity of the coming observer LISA [84,85], as expected for this class of models [34]. We thus conclude that strong GW signal in the LISA with spectrum controlled by the plasma sound waves is a generic prediction of Baryogenesis with relativistic bubble walls. This is in sharp opposition with the general expectation that usual EWBG demands slow walls, and thus suppressed signals. As a final comment, it should however be noticed that the current simulations do not directly provide a solutions for the regime of large α, and we only have an extrapolation of the numerical result. Thus, the conclusion above should be taken with a grain of salt.

Conclusion
In this study we have presented the first explicit realization of the baryogenesis and DM production during electroweak phase transition for ultra-relativistic bubble expansion. The work is based on the proposals in [1][2][3]29] where new heavy particles are produced in plasma−bubble wall collisions. We have shown that the model with SM extended by a real singlet with a Z 2 symmetry can indeed lead to ultrarelativistic bubbles, where the Lorentz factor γ w can reach the values ∼ 10 5−6 . Such fast bubbles can appear if the symmetry breaking occurs in two steps: first discrete Z 2 is spontaneously broken and in the second step electroweak symmetry breaking is accompanied by Z 2 restoration (0, 0) 0). We find that there exists a region of parameter space where the nucleation temperature can become as low as 1 − 2 GeV and the collision of the bubble wall with the plasma particles can lead to the non-thermal heavy particle production with the masses up to ∼ 10 TeV. Interestingly we find that the mechanism is most efficient for relatively low masses of the singlet field M s (v EW , 0) ∼ 70 − 100 GeV, close to the region excluded by the Higgs invisible decays. Subsequently, this region of parameter space will be probed by HL-LHC ( [70,100]) in the singlet production mediated by off-shell Higgs boson. By noting the slight Z 2 breaking, s, if produced, can decay into bb in collider experiments. Depending on the size of the breaking displaced vertices of bb may be probed. We find the typical bubble radius parameter of the order of R ∼ (10 −4 − 1)H −1 so that stochastic gravitational background signal becomes observable at GW experiments like LISA [84,85].
The model necessarily requires tuning ∝ (T nuc /m h ) 2 which numerically turns out to be of the order of 10 −4 − 10 −2 (using Giudice-Barbieri measure) for successful baryogenesis and DM production mechanism. In spite of this we believe it can provide a useful guidance for more appealing models where these hierarchies can appear naturally.

A The bounce in two dimensions
In this paper we studied numerically the phase transition from the minimum (0, v s ) (or in the vicinity of it) to (v EW , 0). The bounce computation can be done using existing codes for example FindBounce or CosmoTransition. However we have found that in the regime of long supercooling where the potential around the false vacuum is very flat and, the existing codes are often not stable and lead to numerical errors. Thus we have developed our own code (more stable for the flat potentials), following the procedure described in [101], while cross-checking the available values with FindBounce.

A.1 Computation of the bounce profile
In this Appendix, we briefly review the standard computation of the bounce action with only one field before going to describe the algorithm we used for the same computation but for the case of two fields PT. In order to compute the vacuum tunneling probability from the false vacuum to the true one in d dimensions, we need to minimize the Euclidean action given by It is known that the field configurations leading to the minimal action are the ones that exhibit an O(d) spherical symmetry, then the so-called bounce solution is the solution of the following Cauchy problem where we have chosen the false minimum to be at φ = 0. If we interpret the parameter r as a time and φ as a position, this problem becomes formally equivalent to the evolution of a mechanical ball in a potential −V [φ] undergoing a friction d−1 r dφ dr , released with vanishing velocity and stopping its evolution for r → ∞ at φ = 0. It is well known that this problem can be solved by applying numerically an overshoot/undershoot method on the position of the released point. Releasing the ball too close to the true vacuum would induce an overshoot configuration (the ball would continue after crossing φ = 0), we would thus shift the release point toward the false vacuum, while releasing it too close would end up in an undershoot configuration (the ball would never reach φ = 0 and starts oscillate around the minimum of −V [φ]) and we correct it by shifting the release point farther from the false vacuum. Iterating between those two situations, we are able to find the correct release point and obtain the bounce solution.
It is well known that the case of a PT triggered by temperature fluctuation at temperature T is formally equivalent to imposing a periodicity T −1 in the imaginary time t E , which imposes the following constraint on the field and the computation of thermally induced phase transition thus amounts to take d = 3 in the above equations.

A.2 Bounce action in two dimensions and path deformation
The problem complexifies when the transition involves many fields. Here there is no straightforward intuition for the path followed by the fields in field space during the tunneling. One can think that a straight line, connecting the two minima, could be a reasonable guess, but it turns out that it cannot be considered as a good approximation of the euclidean action 7 . Here we thus describe the algorithm [101] to find the right path in field space. In a multi-field case, the Eq.(83) becomes Since in this case an overshoot/undershoot procedure cannot be easily applied, the idea is to reduce the problem to one dimensional tunneling. In order to do so we start guessing the path, φ g (x), where x is now to be understood as the parameter that measure the distance along the path, i.e. the so-called curvilinear abscissa.
For the present case, if we parametrize the path in the field space as (h(t), s(t)) = (t, f (t)) ≡ (h, s(h)) it is defined as Here, we have been able to separate the dynamics along the parallel and perpendicular direction in such a way the first equation defines a new undershoot/overshoot problem, that we solve to obtain the value of the escape point, φ 0 (x ), and the Euclidean action corresponding to the potential along the path considered φ g , as in Fig.9. On the other hand, the second equation can be seen as a condition that the bounce solution has to satisfy and can be thought as a force field acting on the path, defined as following The right path will be the one where N is vanishing. The algorithm proceeds iteratively: first we guess a path, the straight line connecting the two minima, then we find the bounce solution along this path, we compute the normal force and deform the guessed path according to it. In practice, to define the path at the step n, φ n , we need to solve for the bounce profile for the path at φ n−1 , extract the escape point x ,n−1 , that is (h(x ,n−1 ), s(x ,n−1 )) in field space, we then discretize the path in the interval x ∈ [0, x ,n−1 ], creating a series ( φ n−1 ) j , for j = 1, ..., N and a series of values for the normal force ( N n−1 ) j . We then shift each point of the discretized path by ( φ n ) j = ( φ n−1 ) j + ρ( N n−1 ) j j = 1, . . . , N. In the end, we fit a path φ n along the shifted points from ( φ n−1 ) j . The procedure of deformation of the path will produce a series of paths φ i [x], over which we compute the Euclidean action according to Eq.(82) at each step of the deformation, like in Fig.10. The algorithm stops when the difference in the bounce action, S 3 , between two successive iterations is below some imposed precision. At a definite temperature T , we start by identifying the two minima, the false and the true ones and will keep the false minimum fixed during the whole procedure of deformation. Generally, especially when we have a sizable amount of supercooling, the escape point is just behind the barrier, so the escape point (v (T n ), v s, (T n )) will be different from the, zero-temperature, EWSB vacuum, but when the tunneling happens the system will classically roll down towards the global minimum, as we can see from Fig.10. We do not track the evolution of the fields profile after the tunneling.

B Supplemental numerical results
In this Appendix, we present all our supplemental numerical results. First, tough we focused mostly on a more weakly coupled part of the parameter space, we would like to compare our findings with the ones in the Ref. [47] and argue that we observed only small changes, due to the inclusion of loop-corrections and Daisy resummation. On Fig.11 we make a reproduction of the scan of the Fig.6 of [47] using our potential and emphasize the close similarities. The relations between the parameters κ, η and  Figure 11: Here is presented the same results found in Fig. 6 of Ref. [47]. It has to be noted that these results are obtained with only the thermal potential and without Daisy resummation, i.e. without thermal masses. The relation with our parameters is (λ s , λ hs ) = (η, 2κ) and M s (v EW , 0) = 300 GeV.  Table 3: Same as Tables 1 and 2 In the main text, we studied specifically the case where the parameter m s = 125 GeV, we observed that for this value, the region of deep supercooling displayed small masses of the singlet in the real vacuum, being on the verge of detection due to h → ss at M s 62 GeV. We also concluded in section 5 that this region was closing around M s ≈ 75 GeV. We could wonder if this conclusion would change if we modify the value of the parameter m s , and if so in which direction. On Fig.12 we show similar plots than in Fig.1 and 3 for the case of m s = 150 GeV. Thus, increase the value of m s pushes the deep supercooling region to M s ≈ 90 GeV, at the price of increasing the portal coupling λ hs . However, we can observe on the last plot of Fig.12 that the typical tuning remains roughly the same and that we can still trust our naive (T nuc /m h ) 2 for an order-of-magnitude estimate of the tuning.
On the other hand, we also observed that decreasing the parameter m s to ≈ 100 GeV was pushing all the deep supercooling region inside M s 62 GeV, which is thus strongly disfavored by colliders. We hope that this trend can be extrapolated to larger values of m s , until we hit perturbativity bounds for λ hs .
We could also wonder what happens at the upper boundary of the deep supercooling region, as we have observed on Fig.1 a sharp decrease in the supercooling allowed around v s 200 GeV (for m s = 125 GeV). This transition regime can be understood if we plot the explicit S 3 /T functions on Fig.13. Comparing the plot in Fig.13 with the one in Fig.5, we see that as we decrease v s , the full pattern of S 3 /T is shifted toward smaller values. At some critical point around v s ≈ 200 GeV, the nucleation becomes controlled by the first minimum in the function S 3 /T and not by the disappearance of the barrier. This largely suppresses the possibility for large supercooling.
Finally, in Table 3, we provide the value of the velocity, reheating and nucleation temperature for m s = 150 GeV and v s = 175 GeV that was used in Fig.6.

C The coefficient of NLO pressure
In this appendix we will review the calculation of the friction coefficient for the NLO pressure for EW phase transition. We will follow closely the discussion in [67] and report the quantity abc ν a g a β c C abc (92) where ν a = 1(3/4) for a a boson (fermion), β c ≡ Mc M Z and C abc stands for the couplings appearing in the vertex. Normalization of the C abc coefficient is the following: for a chiral fermion coupled to the vector field the amplitude for the process ψ → ψA sof t is equal to Where in the relation Eq. (93) is written only for one polarization of the vector field. Similarly for the scalar field

D Domain wall collapse
Our main discussion was focused on the two step phase transition (0, 0) → (0, v s ) → (v h , 0) where the first phase transition is Z 2 breaking. Obviously during such a phase transition domain walls will be formed which can drastically modify the cosmology of the system. We can avoid the stable domain walls if we assume some small Z 2 breaking, however in this case the question rises about the timescale for the stability of the domain walls. This is particularly important since recently it was shown [102] that for singlet extension of the SM the domain walls (if still present) will become seeds of the secondary phase transition (0, v s ) → (v h , 0) and will dominate the phase transition. We will follow closely the discussion in the section 4.2 using only the tree level potential and the thermal corrections to the masses. Then the Z 2 breaking phase transition will occur at the temperatures which is a temperature of the domain wall formation. The domain wall mediated transition will happen at the temperature T w which is found to be order one different from T Z 2 . The exact mechanism of the transition depends on the values of the couplings and can proceed either with the classical rolling or 2D bounces localized on the domain wall. The temperature when the classical rolling can start is reported in Ref. [102] and is equal to The nucleation temperature (T 2D w ) of 2D bounces should be found numerically (Ref. [102]) however it will be obviously smaller than T crit (of (v s , 0) → (v h , 0) phase transition). At this point we can safely ignore the seeded phase transition effects if all of the domain walls annihilate in the interval of temperatures where T * is the temperature when the seeded phase transition will be completed and it is obviously less than T crit of EW phase transition. Let us estimate how strong should be the bias ∆V between the potential energies of the two minima of Z 2 potential so that all of the walls can disappear. For these estimates it is sufficient to assume that there is order one difference between T * and T Z 2 , which is generically the case. The critical radius (above which) areas with true vacuum will start to expand is roughly where σ is surface energy density of the wall. So the domain walls will exist on the time scale of where u is velocity of the wall motion. The change of the temperature during the wall annihilation will be roughly ∆T ∼ T H∆t w .
So that if ∆t w H 1 ⇒ ∆T T the wall annihilation happens almost instantaneously. Assuming σ ∼ T 3 crit and H ∼ Balancing the pressure against the friction forces ∆V ∼ uT 4 crit we can estimate the velocity and then the condition for the quick wall annihilation becomes which is not restrictive at all.