Dark Matter production from relativistic bubble walls

In this paper we present a novel mechanism for producing the observed Dark Matter(DM) relic abundance during the First Order Phase Transition (FOPT) in the early universe. We show that the bubble expansion with ultra-relativistic velocities can lead to the abundance of DM particles with masses much larger than the scale of the transition. We study this non-thermal production mechanism in the context of a generic phase transition and the electroweak phase transition. The application of the mechanism to the Higgs portal DM as well as the signal in the Stochastic Gravitational Background are discussed.


Introduction
Cosmological observations conspire to suggest the existence of a massive, undetected, dark component permeating the universe [1], this is the Dark Matter (DM) phenomenon. One of the earliest candidate for this DM, the celebrated WIMP component, demands that the Standard Model (SM) is coupled to the DM, whose stability is guaranteed by a symmetry. This interaction leads to quick thermalisation between the DM and the SM. In this mechanism, known as thermal Freeze-Out (FO), thermal relic density is naturally fixed via the decoupling of the SM-DM sectors, when the rate of interaction can not compete any more with the expansion of the universe [2,3,4]. The requirement that this relic density matches the observed abundance imposes a relation between the DM-SM coupling and the mass of the DM candidate. In this context, the surprising and exciting coincidence that weak coupling and TeV scale DM candidate are consistent with observed DM abundance is known as the WIMP miracle. Moreover, unitarity considerations on the coupling governing the scattering of DM provide an upper bound on the mass of the DM candidate [5], the Griest-Kamionkowski (GK) bound of O(100) TeV. However, today, many WIMP models have been excluded due to the bounds on the DM-nucleon scattering set by the direct detection experiments [6,7,8,9,10].
To diversify the range of possibilities inside the (coupling-mass) parameter space, many alternatives to FO have been proposed, as for example; freeze-in [11,12,13], forbidden freezein [14], super-heavy particles decay [15,16]. Several proposals also take advantage of the possibility of an early First-Order Phase Transitions (FOPT) occurring in the universe, with many different consequences on DM abundance [17]. Phase transitions offer a way to fix the final relic abundance via the VEV flip-flop mechanism [18,19,20], by modifying the stability of DM candidate [21,22], through the injection of entropy [23,24,25,26] or also via non-thermal production mechanism [27]. More recently, the mechanism of bubble filtering (BF) [28,29,30] was proposed as a way to go around the GK bound and produce ultra-heavy DM candidate with the observed abundance.
In this paper, we would like to present a new mechanism of DM production, occurring during strong FOPT's with ultra-relativistic walls and effective when DM is connected via portal coupling to the sector with FOPT. In [31], authors showed that an ultra-relativistic wall, with Lorentz factor γ w 1, sweeping through the plasma can excite degrees of freedom of mass up to M ∼ √ γ w T nuc , possibly producing out-of-equilibrium particles, mechanism that we call Bubble Expansion (BE) production. In this paper, we would like to show that those produced particles can be stable and thus constitute viable DM candidates. In addition to the possibility of evading the GK bound and thus possibly providing ultra-massive DM candidate, the relic density of these particles is set by the hierarchy between the mass of the DM and the scale of the transition and thus evades the exponential sensitivity typically showing up in the relic abundance controlled by FOPT's. In this context, a simple model for the DM sector perhaps is a real singlet scalar field stabilized by a Z 2 symmetry coupled via the portal coupling to the scalar field (Higgs) undergoing FOPT. In this minimal setting if the Higgs is SM field [32,33,34,35,36,37], FO mechanism is under strong constraints and the direct detection experiment have excluded most of the parameter region below the TeV range.
A similar production mechanism, the Bubble Collision mechanism, takes advantage of the large excursion of Higgs vacuum expectation value (VEV) during the collision of relativistic bubbles. It was first hinted in [38], predicting a production of particles as massive as M ∼ γT . This was shown to be too optimistic in [27] where only the vector and fermion DM candidate were considered as promising DM candidate, however for the scalar DM we find that the mechanism of production of DM via bubble collision of [27] is completely subdominant compared to BE, presented in this paper.
We will show that our production mechanism can proceed even with very massive DM candidate, thus possibly evading the direct detection experiment bounds, even if the coupling to SM is strong. However an irreducible prediction of the mechanism, which takes advantage of a strong FOPT, is the large imprint left in the Stochastic Gravitational Waves Background (SGWB). Such SGWB signal could be detected in forthcoming GW detectors such as LISA, advanced LIGO, BBO, DECIGO, etc, offering an alternative way to study DM production. This paper is organised as follows: in Section 2 we present the production mechanism and the amount of relics produced after the passage of the wall. In Section 3, we present first the maximal amount of DM abundance that can be produced via BE mechanism, and then discuss three ways of accommodating the parameter space to the observed DM abundance; 3.1, we discuss how annihilation can modify the early relics abundance, in Section 3.2 we discuss how some amount of supercooling modifies the relative FO and BE abundances and, finally, in Section 3.3, we discuss the case of very massive DM candidate in the absence of FO relics. In Section 4, we specialize to the Electroweak Phase transition (EWPT) and discuss the allowed range of parameter providing the observed relic abundance. In Section 5, we expose the unavoidable gravitational signature expected by such mechanism. Finally, in Section 6 we conclude.

DM production in the Bubble Expansion
Let us introduce the Lagrangian for the minimal model which suffices for the illustration of the advertised effect where h is a complex scalar field obtaining a non-vanishing VEV via the phase transition and V (h) is its potential. We will not specify the form of V (h) in this paper, but will assume that it leads to the first order phase transition in the early universe. This field h can be the physical Higgs, and thus the phase transition(PT) is electroweak (EWPT), or a new Dark Higgs, and then the transition happens only in the Dark Sector (DS). On the top of it, we introduce a DM candidate φ, that for simplicity we take to be only a single scalar field stabilized by a Z 2 symmetry, with Lagrangian of the form We have assumed that DM candidate is coupled to the symmetry breaking sector via the portal coupling which is the simplest and most natural non-gravitational connection between the symmetry breaking sector and the DM candidate (for review on portal DM, see [37]). We will also assume λ > 0 in order to make sure the potential is bounded from below. We will be mostly interested in masses of the DM candidate φ much larger than the Higgs scale, M φ m h . As a consequence, the abundance of φ in the plasma at the moment of the transition is Boltzmann-suppressed and can be ignored in the dynamics of the transition. We thus neglect the quartic part of φ potential in the discussion as well as the change of M φ due to the transition with v the VEV of the Higgs-like field in the zero-temperature true vacuum, v ≡ h . The hierarchy M φ m h , v introduces the usual tuning of the scalar mass into the model if λM 2 φ /(16π 2 ) m 2 h , v 2 (similar to the SM Higgs mass hierarchy problem), but in this paper we will not try to present a model where this hierarchy can be obtained naturally.

Dynamics of the bubble wall after nucleation
Let us now turn to the dynamics of the transition triggered by the Higgs-like field h. As already stated above, we will focus on the regime with ultra-relativistic bubble wall expansion with where v w is the wall velocity at the bubble center frame. This regime is favoured when the transition is strong enough to develop at least some amount of supercooling. Indeed the condition for the acceleration of the wall is fulfilled if the release of energy ≡ ∆V (using the zero-temperature potential) is larger than the pressure ∆P (computed using the zero-temperature minima) exerted on the wall by the plasma. In the relativistic limit, at the leading order (LO), the pressure is equal to [39,40] with g i the number of degrees of freedom (d.o.f.) contained in the plasma 1 , ∆m 2 = m 2 broken − m 2 symmetric and c i = 1(1/2) for bosons (fermions) and T nuc , the nucleation temperature, is the temperature when there is roughly one bubble per Hubble volume. Eq.(3) can be considered as an upper bound on the pressure [41] and the bubble becomes relativistic if [42] > ∆P LO (Relativistic wall condition).
As pressure scales like P LO ∝ v 2 T 2 and release of energy like ∝ v 4 , supercooled transitions, like in nearly conformal dynamics [43,44,45], drive the wall to ultra-relativistic regimes. Note that if no other contribution is present, the bubbles satisfying Eq.(4) become runaway (permanently accelerating). If some gauge field acquires a mass during the phase transition, it is known that the Next-To-Leading order (NLO) correction to the pressure [46], due to the emission of ultra-soft vector bosons, scales like γ w where g gauge is the gauge coupling and g i counts the number of degrees of freedom. This pressure will stop the acceleration of the wall and wield a terminal velocity with final boost factor γ w,MAX . Before proceeding further let us estimate the maximal velocity (or γ w,MAX factor) the bubble wall will reach before the bubble collisions. As we have seen from Eqs. (3) and (5), the discussion changes depending on the presence of phase-dependent vectors.
1. Runaway regime: When the PT does not involve phase-dependent vectors, there is no NLO pressure and the wall keeps accelerating until collision. The γ w at collision is [47,31] where R is an estimate for the bubble size at collision and R 0 is the bubble size at nucleation, β is the inverse duration parameter of the transition and M pl ≈ 2.4×10 18 GeV the reduced Planck mass.
2. Terminal velocity regime: When the PT gives a mass to vectors, the pressure becomes dominated by the emission of ultra-soft bosons and quickly wield a terminal velocity of the form where in the last step we have to take the minimal of the two values, since the bubble collision can happen before the terminal velocity regime is reached.
The last source of pressure is provided by the production of heavy particles [31] including DM itself Here M Heavy is the typical mass of the heavy particles. This additional contribution can stop as well the bubbles from being in the runaway regime (see for examples [31]). At last before we will proceed to the calculation of DM production, let us define a reheating temperature after the completion of the phase transition, which is approximately equal to where is the latent heat released during the FOPT. Generically we expect T reh ∼ T cr ∼ v, with T cr the critical temperature when the two minima are degenerate. Note that in the regime of large supercooling α 1 there will be a hierarchy between the nucleation and reheating temperatures T reh T nuc .

Production of DM via the bubble wall
After those preliminaries, we can now go to the production mechanism itself. In the wall frame, h particles hit the wall with typical energy and momentum E h ∼ p h z ∼ γ w T nuc . The VEV of the h h = v(z), varying along the wall, induces a new trilinear coupling of the form λv(z)hφ 2 that did not existed on the symmetric side of the wall. It was shown in [31] that, in such a situation, the transition from light to heavy states h → φφ has a probability of the form (see also Appendix A for the details of the computation) L w is the width of the wall which is approximately is a difference of momenta between final and initial state particles in the direction orthogonal to the wall. Immediately after the production, the typical energy of each φ in the bubble center frame isĒ if p h z M φ , see Appendix A. This is much larger than either m h or M φ . As a consequence, inside of the bubble, a non-vanishing density of non-thermal φ accumulates. Thus, this "Bubble Expansion (BE)" produced density of φ, in the wall rest frame, takes the following form v w = 1 − 1/γ 2 w is the velocity of the wall, and f h (p) is the equilibrium thermal distribution of h outside of the bubble. This writes , as the Higgs-like field should be at equilibrium with SM.
Using Boltzmann distribution as a simplifying assumption, we can perform the integral in Eq. (12), obtaining , the density in the plasma frame, in the limit of fast walls, becomes The (12). We can see that in the limit the exponential goes to one and the density becomes independent of the velocity of the wall v w . The step function is an approximation of the transition function which depends on the exact shape of the wall. We report it for different wall ansatzs in Appendix A. It is important to note that in the regime ∆p z L w 1 the step function presents a good approximation and the results are independent of the wall shape as expected from the Heisenberg uncertainty principle. However, if the inequality Eq. (15) is not satisfied and we are in the regime then the wall shape effects start to become important. We discuss this wall suppression for the tanh and gaussian walls in the Appendix A. We find that generically the deviations from the naive step function are exponentially suppressed, so that expression in Eq.(13) can be used as an estimate in the transition regime Eq. (16). At last, for the particle production gets additional suppression by the usual Boltzmann factor. From now we will keep working with expression (12), keeping in mind possible departure from pure exponential suppression behaviour.
The final number density of heavy non-thermal DM, in the unsuppressed region is of the form From the previous discussion, we see that an ultra-relativistic wall of FOPT sweeping through the plasma will produce heavy states, via portal coupling of Eq.(2). Assuming no subsequent reprocessing (thermalisation, annihilation, dilution by inflation ...) of the relic abundance, the nowadays abundance of Bubble Expansion (BE) produced DM is given by where T 0 is the temperature today, ρ c is the critical energy density and g S0 (g S (T reh )) is the entropy number of d.o.f. today (at the reheating temperature). As a consequence, plugging the expression Eq.(18), the final relic abundance today writes and we see that the produced relic abundance is controlled by the quantities So far we have shown that a bubble with Lorentz factor γ w sweeping through the plasma can produce massive states up to mass M 2 φ γ w T nuc /L w , where L w ∼ 1/v is the width of the wall. The maximal value of the γ w factor depends on the particle content of the theory (particularly the presence of the gauge fields) which influences the largest DM mass which can be produced. We can estimate this maximal mass by considering two generic cases of the bubble expansion.
1. Runaway regime: According to this maximal boost factor in Eq.(6), the maximal mass M MAX φ that can be produced, by the sweeping of the wall, scales like We will study Dark Sectors of this type in Section 3.

2.
Terminal velocity regime: Similar considerations from Eq. (7) give where we assumed, as in the remaining of this paper, that g i g 3 gauge ∼ O(1). Above this maximal mass, the production of DM becomes exponentially suppressed according to 2vTnucγw , as we have seen in Eq. (14). We will study a transition of this type in the context of EWPT in the Section 4. and Ω today BE,φ = Ω today obs define 4 regions. In I, BE abundance is dominant and FO is not enough to account for the observation. In II, FO is too large, but BE is still dominant. In III, both BE and FO are too large, but FO is dominant. Finally, in IV, FO is dominant, and BE is not enough to account for Ω today obs .
The final relic abundance produced during BE has to compete with the relic abundance coming from FO, which provides a final relic abundance roughly of the form Notice that this component exists if the reheating temperature of the Universe after inflation is higher than M φ and if φ couples to the thermal bath not too weakly so that φ is produced from the thermal scatterings. We assume this component in most parts of this paper. However, we will remove this assumption in section 3.3. The ratio of the nucleation temperature T nuc over the reheating temperature T reh in Eq. (20), originates from the fact that the heavy particles are actually produced at the nucleation temperature, but that the release of energy reheat the plasma at T reh , providing the new initial condition for the evolution of the universe. Strong FOPT's are often accompanied by long supercooling and thermal inflation [48,49], leading to the hierarchy between T nuc and T reh and strong suppression of the abundance. We will see that this new suppression factor can be useful in the range of parameters where the final relic abundance is overproduced, as illustrated on Fig 1: in the region II, where the BE abundance is dominant over FO, but both of them are too large to account for Ω today obs and I, where FO is not large enough. In this range, dilution related to thermal inflation can reduce the overproduced relic abundance to Ω today obs .

Dark Sector PT production of DM
In the previous section, we have presented a new mechanism of DM production. However it is important whether this mechanism can lead to the observed relic abundance. In order to consider the phenomenological relevance of our mechanism we will use the toy model presented in Eq.(1), which can perfectly constitute a viable model of DM. We consider the field h as some scalar field experiencing the phase transition at some scale v. Let us look at the nowadays relic abundance presented in Eq. (20). The results are presented on the Fig. 1 for v = 200 and 2×10 4 GeV. Generically we can define four regions as follows: in region I, the abundance is underproduced via FO, but largely overproduced via BE. The region IV is the symmetric situation, where the BE is small but FO is very large. In region II(III), both FO and BE are overproduced, but BE (FO) production dominates over FO (BE): Very naively these equations indicate that none of the regions leads to a viable phenomenology. However we have not yet taken into account few possibilities on the initial conditions as well as the evolution of ρ φ /s which can make some parts of those regions viable.
To be more precise, we will study three possibilities; in the regions where DM is overproduced annihilation processes can reduce the DM density back to the observed relic abundance, as this can be for example the case in region I. We discuss this possibility in the Section 3.1. As we already hinted above, another process which can reduce the DM density is a brief period of inflation during the FOPT, which happens if the nucleation temperature is significantly lower than the reheating temperature. This leads to the reduction of the overall DM density and as a result opens up some parameter space, typically inside of region I and II of Fig. 1. We discuss this effect in the section 3.2. At last in the case that the thermal history begins with a reheating scale below the FO temperature 2 , φ never reaches thermal equilibrium after the reheating and is (almost) not produced via FO. We discuss this possibility in the Section 3.3.

Late time annihilation
In the previous Section 2, we showed that if a relativistic bubble goes through the plasma, it can produce DM relics, possibly very over-abundant. On Fig. 1 we saw that, in region I, the FO contribution was not large enough to account for the observed DM abundance, but that on the contrary, BE production was extremely large. As a consequence we would like to track the evolution of the number of DM particles after the initial production. We will see that, as long as the DM density produced is very large, the final density does not depend on the initial density. Thus the physics of this part does not change even if φ is produced enormously from other dynamics e.g. inflaton/moduli decay. 3 Due to this reason, in the following, we make a general discussion which is not specific to the BE production unless otherwise stated. We assume for simplicity that the production happens instantaneously during the radiation domination epoch at T [t ini ] = T reh , and assume that the density just after the production is much larger than that for the observed DM abundance (which is the case of the region I of Fig. 1).
The annihilation cross-section for the process φφ → hh is well approximated as when φ is non-relativistic. Here v rel is the relativistic velocity, and is the average over the distribution functions of φ and h.
for h being the SM Higgs and a real singlet Dark Higgs, respectively. (In the real singlet Dark Higgs case we should take L scalar ⊃ −λφ 2 h 2 /4.) In calculating the average, we have assumed that just after the production, the DM velocity v φ soon slows down due to the scattering with the ambient plasma, and we further assume h soon decays into the SM plasma. When h is the SM Higgs, the assumptions are easily satisfied. The mean-free path in the thermal environment is set by the inverse of Γ MFP ∼ y 2 q λ 2 v 2 E 2 φ T reh where y q is the quark-Higgs Yukawa coupling (This expression is valid in the broken phase. In the case of symmetric phase, the scattering is with Higgs multiplet and the rate is larger.) Here T reh is comparable or larger than the mass of the quark q. Γ MFP is easily larger than the Hubble parameter unless E φ is extremely large. When the dominant annihilation product is a dark Higgs boson, we can still have a sub-dominant portal coupling between the DM and the SM Higgs, via which the kinetic equilibrium can be easily reached. The typical velocity of φ in the kinetic equilibrium is Thus a simple criterion to assess the stability of DM relics is the competition between the expansion rate of the universe, and the rate of annihilation Γ ann . A rough stability condition thus writes Γ ann ∼ σ φφ v rel n φ < H (Stability condition).
If this condition is violated the annihilation gradually takes place even if T reh is below the FO temperature ∼ M φ /20, as discussed in the Wino and Higgsino DM cases [52,53].
To evaluate the final abundance after the annihilation, we can solve the integrated Boltzmann equation (by assuming kinetic equilibrium as in the case of the WIMP): n eq (M φ T /(2π)) 3/2 exp (−M φ /T ) is the number density in the equilibrium, and the annihilation rate is given by being the thermal average and S eff is the Sommerfeld enhancement factor, i.e. the boost factor from the interacting long-range force. We assume the force potential between the φ pair distanced by r as where α med (m med ) is the messenger coupling (mass). For the Higgs-mediated force discussed in this section, we have The analytic approximation of the enhancement factor is given by [54,55] (See also Refs [56,57,58,59]) To solve numerically the Boltzmann equation, we set the initial condition of n φ [t ini ] 0.2 eV × s/M φ , i.e. much larger than the corresponding value of the observed DM number density. Here s is the entropy density. The Boltzmann equation can be solved to give Fig.2 where we plot the time evolution of the number density with n φ [t ini ] ∼ 40 eV × s/M φ . Indeed, we find that even when initially there is too large number density, with large enough coupling (and thus large annihilation rate), the number density decreases significantly within one Hubble time. We obtain suppressed abundance in the end (right top panel). On the bottom (left top) panel we can see that if the coupling is not very large this is not the case (if M φ < T reh /20, φ is thermalized soon and FO happens).
In Fig.3 with h being the real singlet Dark Higgs, we represent the numerical result giving Ω φ h 2 = Ω DM h 2 ≈ 0.1 [60] by taking v = T reh = 50, 100, 200, 400 GeV from top to bottom, with the initial condition set as n φ [t ini ] = 40 eV × s/M . We see that at lower mass range the predictions do not depend on T reh , which represents that the FO takes over since T reh > T FO ∼ M φ /20. The FO prediction is displayed on Fig.3 (and 4) by the dotted orange line. On the larger mass range, the late time annihilation becomes important and reduces the abundance relevant to T reh .
In fact, we can explain the final number density, n φ , in this region from the condition This condition is similar to the freeze-out condition for the ordinary WIMP: the annihilation should end when the rate becomes comparable to the Hubble expansion rate. We obtain λ = λ ann ∼ 0.53(g 4 rS eff ) −1/2 √ C g (T reh ) 103.5 from the condition that the φ abundance composes an r fraction of the observed dark matter abundance, Ω φ = rΩ DM (and we are now focusing on r = 1.) Notice again that to use Eq. (34) we have assumed T FO > T reh , otherwise the DM is thermalized and then usual FO takes place after a certain redshift. From the numerical fit by solving the Boltzmann equations, we obtain C = [0.1 − 1] depending on the initial condition. If the initial n φ [t ini ] is larger C becomes larger approaching to 1. In particular for our bubble wall scenario, we may have a very large n φ (t ini ) and, in this peculiar case, C is almost 1. So far we have been agnostic regarding the coupling of the DM to the SM sectors. We just have assumed that DM couples to the scalar field h to which it annihilates into. However, to be in kinetic equilibrium, the DM should somehow couple to the SM plasma. This leads to the possibility of detecting DM with direct and indirect detection experiments. In particular, when h is the SM Higgs boson, the coupling to nucleons is controlled by the coupling λ. The case where h is the SM Higgs multiplet is shown in Fig. 4, where the difference from Fig. 3 is that we fixed v = 174 GeV, g 4 = 1 and m h = 125 GeV. We adopt the bound XENON1T experiment [61] from [37] (The Purple region above the purple solid line), which is extrapolated by us to multi-TeV range. The green dashed and blue dotted lines represent the future reaches of the XENONnT [62] and DARWIN [63], respectively, which are also adapted and extrapolated from [37]. The Cerenkov Telescope Array (CTA) reach (by assuming the NFW distribution of DM) is adopted from [64] and also extrapolated by us. Consequently, the predicted parameter region can be fully covered in the future direct detection and indirect detection experiments such as XENONnT, DARWIN and CTA. Interestingly, since the predicted black lines are parallel to the direct detection reaches in the late time annihilation region, T reh corresponds to the DM-Nucleon interaction rate. If the DM is detected in the direct detection experiments, which implies the interaction rate is measured, we can tell the reheating temperature assuming late time annihilation.
Here we notice that the contribution of the Sommerfeld enhancement may be as large as S eff − 1 = O(10%) when the mass is large. Usually in the (SM) Higgs portal dark matter model, the Sommerfeld enhancement is negligible due to the small Higgs dark matter coupling, α med ∝ ( λv M φ ) 2 suppressed by the heavy dark matter mass. In the late annihilation scenario, since we need larger λ than conventional FO and smaller v φ , we have larger S eff . and DARWIN [63], respectively. The lines are adopted from [37]. The Cerenkov Telescope Array (CTA) reach (by assuming the NFW distribution of DM) is adopted from [64].
As a conclusion of this section, let us, finally, come back to the BE production. We have seen on Fig. 1 that in the region of parameter with large coupling and DM mass in the TeV range, the FO is subdominant and BE is largely over-produced, this was the region I of Fig. 1. This is exactly the setting we studied in this section and the result displayed on Fig. 3 can be used for the dark sector PT. Also Fig.4 can be straightforwardly extended to the EWPT, if we assume that some modification of the SM wield a strong enough EWPT. We will discuss this possibility further in Section 4.

Dilution by supercooling
In Section 3.1 we saw that even if the DM is over-produced by the wall, the relic abundance can be reduced by the reaction φφ → hh. For the case of v ≈ 174 GeV, this opened up the range of values M φ ∈ [1, 10] TeV and λ ∈ [0. 3,10], which is normally with too small abundance in usual FO. In this section, we would like to account for a second effect, which is the dilution induced by some amount of supercooling. Indeed, if some low-scale thermal inflation [48,49] occurs due to the supercooling, a possibly large hierarchy between the reheating temperature and the nucleation temperature can occur.
During the thermal inflation [43], the expansion factor scales like a ∝ e Ht and the temperature like T rad ∝ e −Ht , the FO abundance is a non-relativistic fluid scaling like Ω FO ∝ T 3 . As a consequence, the FO abundance receives a further Tnuc T reh 3 suppression factor with respect to usual cosmology evolution. Summing both FO and BE contributions the total relic abundance will be approximately given by (we are assuming M φ 20T reh ) When BE contribution and FO contribution are small, the thermal production may become dominant, especially with T reh 1/20M φ (see Eq. (29)). Assuming an instantaneous reheating after bubble collision and negligible non-thermal production of φ via bubble collision [27], the additional contribution from thermal production takes the form This formula agrees well with the numerical simulation by taking C ∼ 0.9 − 1. Since, around the T FO , this contribution changes exponentially with temperature via n eq , the range of C may be slightly wider, which depends on the detailed process of the bubble collision. Let us also mention that, insisting on dominant BE production (second term of Eq.(35) larger than first term and thermal production in Eq.(36)), perturbativity λ < 4π, maximal mass Eq. (22) and finally current bound on the relativistic species at BBN, impose the following constraints on the broken symmetry VEV of the (Dark-)Higgs: MeV v 10 8 GeV, (scale range).
The upper bound is due to the quadratic dependence of the BE production on the VEV v while the lower bound comes from stringent BBN bound on the number of relativistic species, which demands that our transition happens before T ∼ 1 MeV. On Fig. 5, we display the values of M φ and λ providing the observed amount of DM relics for the various values of the reheating (T reh ) and nucleation (T nuc ) temperatures for the fixed scale v = 2000 GeV. We have also assumed that the bubble wall could reach runaway regime due to suppressed plasma pressure (no phase dependent gauge fields), so that the upper bound for the DM mass in Eq.(22) becomes ∼ 10 8 GeV. These curves were obtained by numerical solution of the Boltzmann equations but qualitatively we can understand the shape of the isocontours as follows: • Let us start with the top left plot on the Figure 5. The orange dashed line corresponds to the usual DM freeze-out. As we can see, it is the case if the DM is lighter than roughly 20T reh and, in this case, the physics of the phase transition plays no role in the final DM relic abundance.
• For heavier masses the isocontours are given by the red dot-dashed triangles. The sides Ω BE,FO of the triangles are fixed by Eq. (35) and correspond to the cases when either Ω BE or Ω FO dominates the total relic abundance. Almost vertical side at M ∼ 20T reh is given by Eq. (36) and corresponds to the thermal production of DM during reheating after bubble collision. Inside the triangle the DM is under-produced and outside, it is over-produced.
• Let us move on to the other plots on the Fig. 5. Multiple triangles correspond to the different values of supercooling. Finally the origin of the black line (continuation of the dashed orange line) can be traced back to the discussion in Section 3.1. In this case the DM is produced by BE mechanism, however the large coupling leads to an efficient annihilation and the final relic abundance is set by Eq.(34) .

Super-Heavy Dark Matter candidate
Another possibility to suppress the freeze-out (FO) density is to assume that the usual inflation reheating temperature T R is too low and inflaton does not decay into the dark matter, so that φ is not produced by reheating and thermal scattering process. 4 At this point, we can GeV. The Red lines correspond to contributions from FO and BE providing the observed DM abundance and that do not undergo annihilation after the transition. The black line is the result of DM annihilation, as in Section 3.1. Roughly when M φ < 20T reh , the DM comes back to equilibrium after the transition and the final parameters compatible parameters are given by the orange dotted line. Let us also emphasize that we assumed runaway regime bubble, with the maximal DM mass given by Eq. (22) completely decouple FO contribution and we are left only with the BE production, so the region of parameter space with large masses M φ or small coupling λ opens up. This condition writes Going back to Eq. (35) and assuming T reh ≈ v, we see that, in this scenario, the final relic abundance is now simply given by the BE contribution with four controlling parameters: v, T nuc , M φ and λ. Assuming vanishing supercooling in order to compute the maximal mass that can be produced, DM with mass as high as GeV can provide the observed DM abundance, Ω BE = Ω obs . The second line was obtained by placing perturbativity bounds on the coupling, λ < 4π. Let us emphasize that this maximal mass has nothing to do with the previously computed maximal mass in Eqs. (22) and (23), where the production was suppressed by wall effects. In this case, the maximal mass originates from the fact that even in the unsuppressed region, the production scales as ∝ 1 Of course, those very large masses can only be activated by the transition if it does not contain gauge boson, according to (22). As a consequence, this possibility most probably can not be realised in the context of EWPT, as the wall quickly reaches a terminal velocity.
Fixing v = 2 × 10 2 GeV, and considering vanishing supercooling, the observed relic abundance is displayed, in the space (M φ , λ) on Fig.1 by the red line dubbed Ω BE = Ω obs .

BE production in EWPT
So far, with the exception of Fig.4, we have been general in our analysis and assumed that h is a generic field undergoing a very strong FOPT. Let us now specialize to the case of EWPT with v ≈ 200 GeV and assume that the transition is strong enough to induce a relativistic wall. During the SM-Higgs transition, gauge bosons W and Z receive a mass and thus contribute to the pressure at NLO order. Thus the wall will inevitably reach a terminal velocity, which puts an upper bound on the maximal DM mass M MAX φ , above which the DM production starts to become exponentially suppressed. In Eq. (23), we have seen that this maximal mass increases with the supercooling : As a consequence, we will study the φ relic abundance in the range We set the lower bound M MIN φ ∼ TeV, below which the usual FO takes over again after reheating if T reh ∼ 100 GeV, and the sub-TeV WIMP Miracle is mostly excluded as mentioned in the introduction. We will also assume that T nuc Λ QCD , otherwise QCD effects can become important and trigger themselves phase transition (see for example [65,66]), so that the longest supercooling will be roughly ∼ T reh Tnuc 10 3 . These assumptions confine the DM candidate mass to be in the range to TeV M MAX φ 10 3 TeV, thus leaving us with a generous range of exploration. However this setting renders the scenario of Section 3.3, with very massive DM, difficult, so we will not consider it. In this Section, we will only consider the two mechanisms of Section 3.1 and 3.2. The coupling λ in the Eq.(2) become the Higgs portal coupling and leads to the direct detection possibilities. Plotting the isocontours in the (λ, M φ ) space similarly to the Figure 5 we have checked the current bounds and future prospects for direct DM detection on the Fig.6. We can see that parts of the parameter space where the annihilation of DM (Black line of 6) plays a role is already probed by XENON1T experiment and parts of parameter space with BE production mechanism will be tested by the future DARWIN and XENONnT experiments, at least partially. The red-shaded region displays the regions of parameter space where the DM is under-produced, while outside of it, DM is over-produced and the observed DM abundance corresponds to the red line boundary. It is instructive to compare these results with the results of the Fig.5 where we have assumed that γ w → ∞ ⇒ M MAX φ → ∞. On left panel of Fig.6, for T reh Tnuc = 15 we can observe two islands of under-production: one at low mass and low coupling, which is exactly the same as on the Fig.5 and the one for the high masses and high couplings. In the later region the DM production from BE receives an additional suppression of the form e − M 2 φ 2vTnucγw , according with the Eq. (14). On the right panel we present a similar plot for T reh Tnuc = 30, however in this case two islands with and without exponentially suppressed DM production are joined.
Note that in our analysis we have included only the factor e −

Observable signatures
It is well known that an unavoidable signature of strong FOPT's, with very relativistic wall, is large a Stochastic Gravitational Waves Background (SGWB) signal, with peak frequency controlled by the scale of the transition f peak ∼ 10 −3 T reh GeV mHz. As an example, the EWPT signal is expected to peak in the mHz range, which is the optimal range of sensitivity of the forthcoming LISA detector. Then the constraint Eq.(37) turns into a constraint on the frequency of the signal 10 −6 mHz f peak 100 Hz (Frequency range) We can also more or less constrain the model parameters for a given reheating temperature or peak frequency. In Fig. 7, we show T reh (and thus f peak by assuming f peak = 10 −3 T reh GeV ) vs the mass range. The parameter region satisfies the constraints of correct DM abundance Eq. (35)≈ 0.1, the dominant BE production (second term of (35) dominant, suppressed thermal production T reh < 1/20M φ ), Eq. (22), perturbativity (λ < 4π), and consistency conditions T reh ≥ T nuc , v ≥ T reh . For the late time annihilation, we can read the relation for mass, λ, and T reh from Figs. 3 and 4. These imply that the observation of the SGWB provides a probe of the parameter range.
Theoretically, two different sources of GW are well understood; the bubble collision [67], dominating the signal in the case of runaway walls (theories with no gauge bosons), and the plasma sound wave [68], dominating in the case of terminal velocity walls, (theories with gauge bosons). Those two contributions have peak intensity and peak frequency of the form [67,68] 5 Ω peak collision h 2 ∼ 5 × 10 −8 100 g 1/3 κ wall α  Hz with z p ∼ 10, k sw is the efficiency factor for the production of sound waves in the plasma, α and β have been defined in Eqs. (9) and (6) respectively, R ∼ v w /β ∼ O(10 −2 − 10 −3 )H −1 is the approximate size of the bubble at collision, and all quantities (T, H, g ) have to be evaluated at reheating. The specific values of the parameters κ wall and κ sw depend on the regime of the bubble expansion: • Runaway wall A large fraction of the energy is stored in the wall of the bubble: [42] This regime produces GW via bubble collision and sound waves mechanism, with bubble collision dominating the signal.
• Terminal velocity In this case, most of the energy of the transition goes to the plasma motion and we have As a consequence GW are dominated by sound wave production. We can see that these two scenarios are quite exclusive: runaway behaviour is dominated by bubble component and terminal velocity -by sound waves. This difference in principle allows discrimination between the two bubble expansion scenarios.
Strong signals are obtained for: 1)large α, which is the consequence of long supercooling and large latent heat, 2)small β, which are obtained for slow transitions and thus large bubbles at collision, and 3)relativistic walls v w → 1. Thus, the same conditions necessary for the BE production of Dark Matter will induce the strongest GW signal. In Fig.8, we present the signal induced by four benchmark point, each representative of a specific regime: P1 (runaway α = 1, β = 100), P2 (runaway α = 0.1, β = 1000), P3 (terminal velocity α = 1, β = 100), P4 (terminal velocity α = 0.1, β = 1000) with T reh = v = 200 GeV. We also represent the GW signal with several T reh in the range corresponding to Fig. 7 by fixing α = 1 and β = 100. As we expect the scaling α∞ α ∝ Tnuc v 2 , we set a suppressed α ∞ = 0.001, due to quite large supercooling that we considered in most of our scenarios. We can see that generically BE mechanism for DM production leads to the stochastic gravitational wave signature in the frequency range Eq. (43), which is well in the reach of the current and future experimental studies

Conclusion
In this paper we have presented a novel mechanism of the DM production. We have shown that the ultra relativistic expansion of the bubbles during the first order phase transition in the early universe can produce a significant amount of the cold relics even if the mass of the DM candidate is much larger than the scale of the phase transition. This, as a consequence, "brings back to life" components that, due to Boltzmann suppression, did not belong to the plasma any more. We illustrate this mechanism on a simple renormalizable model where DM is a scalar . We also took α ∞ = 0.001. The signal-to-noise ratio and the sensitivity curves can be build following the recommendations of [69,70,71,72,73,74,75,76]. Right-The runaway GW signal with fixed α = 1, β = 100 are shown with T reh = 10 −2 , 10, 10 4 , 10 8 GeV corresponding to the parameter range given in Fig. 7.
coupled via portal coupling to the field experiencing the phase transition. When the bubble wall reaches velocities γ w > M 2 φ v 2 the exponential suppression of the heavy particle production disappears and BE mechanism can become very significant in large ranges of parameter space. Thus the produced DM density can be easily dominant. In the simple model presented in the paper both BE and FO contributions to the DM relic density were controlled by the same coupling, however this does not have to be the case for more complicated models, where additional interactions can suppress FO contribution further.
In the absence of FO produced relics, BE mechanism also provides the possibility of supermassive strongly coupled DM candidate, which is a scenario similar to the baby-zillas of [27].
We showed that there are parameter regions where the BE production dominates over the FO production and explains the observed amount of DM in the universe. This opened up the range of Multi-TeV DM with large coupling, thus being more detectable at direct detection (like forthcoming XENONnT and DARWIN) experiments and indirect detection (like the CTA) experiments than the usual FO mechanism.
Our mechanism is also characterized by an unavoidable and possibly observable imprint in the SGWB, with peak frequency controlled by the scale of the transition. The shape of the spectrum can then discriminate between runaway or terminal velocity bubble wall behaviour. Let us also emphasize that if the DM belong to a totally decoupled DS, SGWB signal is the only unavoidable imprint.
we define the kinematics as p h = (p 0 , 0, 0, p 2 0 − m 2 h (z)) Putting together the relevant pieces, the final matrix element squared is With those tools in hand, we can now compute the probability of 1 to 2 splitting. The expression for the probability of transition generically takes the form and putting together Eq. (55), (54), using the kinematics (48) and the large velocity approxi- where the Θ(γ w T nuc − 2M φ ) function appears from the trivial requirement that we need enough energy to produce the two heavy states and Θ(γ w T nuc − M 2 φ L w ) comes from the behaviour of the function sin α/α, suppressing the transition probability for large α.
One can also estimate the typical energy of the produced φ in the bubble center frame.
Here in the last approximation we have used that p h 0 ∼ γ w (1 + v w )T nuc and v w = 1 − γ −2 w .
A.1 Consequences of the shape of the wall So far, we have been assuming that the wall has a linear shape. This provided us with a suppression factor sin α α 2 . However given a wall shape we would have different type of suppression.
The wall shape depends on the Higgs potential and the specific interactions between the wall and the plasma. Here let us assume two types of the wall shape to explicitly calculate the suppression factor from Eq. (51). To this end, we use tanh wall shape and gaussian wall shape to perform the integral.
In general, we can perform the integral by using partial integration of where we have neglect the surface term. For the tanh (z/L w ) case, V tanh (z) = vλ/(2L w cosh 2 (z/L w )). By noting that z integral becomes the summation of residues at poles z = πiL w /2, 3πiL w /2, · · · for ∆p z > 0 or z = −πiL w /2, −3πiL w /2, · · · for ∆p z < 0, we obtain M tanh = sign[∆p z ]πiλvL w ∞ n=0 e −Lw|∆pz|(n+1/2)π = πiλvL w 2 sinh Lw∆pzπ 2 . (61) One finds that this has the exactly same behavior at ∆p z 1/L w but the suppression is rather exponential, ∝ e −Lw∆|pz| when L w |∆p z | 1. This implies that the linear approximation is good when L w |∆p z | 1 but may not be good enough when L w |∆p z | 1.
In the case of Eq. (59) similarly we obtain, where we have dropped the surface term. Again we have the same form as the linear approximation with ∆p z 1/L w but the suppression factor is gaussian.