Non-thermal Dark Matter Production from the Electroweak Phase Transition: Multi-TeV WIMPs and"Baby-Zillas"

Particle production at the end of a first-order electroweak phase transition may be rather generic in theories beyond the standard model. Dark matter may then be abundantly produced by this mechanism if it has a sizable coupling to the Higgs field. For an electroweak phase transition occuring at a temperature T_EW ~ 50-100 GeV, non-thermally generated dark matter with mass M_X>TeV will survive thermalization after the phase transition, and could then potentially account for the observed dark matter relic density in scenarios where a thermal dark matter component is either too small or absent. Dark matter in these scenarios could then either be multi-TeV WIMPs whose relic abundace is mostly generated at the electroweak phase transition, or"Baby-Zillas"with mass M_GUT>>M_X>>v_EW that never reach thermal equilibrium in the early universe.


Introduction
The most popular paradigm for the origin of dark matter (DM) in the Universe is the thermal freeze-out. In that scenario, the dark matter particle with mass M X annihilates into matter with a cross section σ v thermal ∼ 3 × 10 −26 cm 3 /s. This ensures dark matter is in thermal equilibrium with the rest of the plasma in the early universe while T M but decouples when T ∼ M X /20, leaving the relic abundance in agreement with the value Ω X = 0.228±0.027 measured by WMAP [1]. Incidentally, σ v thermal is a generic cross section for a weak scale mass particle interacting with order one couplings, this fact being referred to as the WIMP miracle. In spite of these attractive features, non-thermal mechanisms of dark matter production have also received considerable attention. Examples include right-handed neutrinos produced by oscillations [2], axions produced by vacuum misalignment [3], winos produced from moduli decays [4], and super-massive dark matter (WIMP-zillas) produced during reheating after inflation [5]. These studies allow one to recognize a wider range of possible collider and astrophysical signals of dark matter than what would result from the thermal WIMP scenario.
In this paper we study the possibility of non-thermal dark matter production during a first-order electroweak (EW) phase transition. Bubble collisions at the end of the EW phase transition may give rise to abundant non-thermal particle production when a sizable amount of the energy budget of the transition is stored in the bubble walls, possibly leading to new and appealing scenarios. Many models of dark matter contain a direct coupling between the Higgs and the dark matter candidate fields (MSSM and its extensions, Little Higgs theories with T-parity and DM extensions of the standard model (SM) via the Higgs portal, to name a few). It is thus reasonable to expect that dark matter may be abundantly produced non-thermally at the end of a first-order EW phase transition. Note that, much like in the thermal WIMP case, dark matter would then be a particle with M X ∼ 10 GeV−10 TeV with significant coupling to the SM, thus being within reach of colliders and DM direct detection experiments.
There is however one generic problem with this scenario. Since the temperature of the Universe right after the EW phase transition is T EW ∼ 50 − 100 GeV (for strong transitions T EW may be somewhat lower than 100 GeV), thermalization will typically lead to a washout of the non-thermal abundance, thus rendering the particle production at the EW phase transition irrelevant for the subsequent evolution of the Universe. The wash-out process can nevertheless be avoided in a number of ways, each resulting in a scenario where nonthermal dark matter production is (fully or partially) responsible for the observed dark matter relic density. One possibility, recently outlined in [6], is to allow for a few e-foldings of inflation prior to the beginning of the transition (which can happen for very strong EW phase transitions), diluting the plasma and leaving the Universe partially empty. If the reheating temperature after the phase transition is low, T RH ≪ 100 GeV, it may be possible for a dark matter candidate with weak couplings to the Higgs field and mass M X ∼ 100 GeV to remain out of thermal equilibrium after the EW phase transition. In this paper we investigate other scenarios allowing for a survival of the non-thermal abundance.
One possibility corresponds to the case of relatively heavy (multi-TeV) dark matter: for M X 1 TeV, dark matter will be very non-relativistic when the EW phase transition takes place, and the decoupling/freeze-out temperature T FO will satisfy T FO ∼ M X /20 T EW .
Then, heavy dark matter produced non-thermally through bubble collisions may remain out of thermal equilibrium after the EW phase transition (or at least wash-out will be partially avoided). Another possibility is that bubble collisions produce super-heavy dark matter, M X ∼ 10 6 -10 8 GeV, a scenario we call "baby-zillas". We argue this may be possible for a very strong EW phase transition and dark matter having a large coupling the Higgs. In order for baby-zillas with M X ≫ v EW to be a viable dark matter candidate, they must have never reached thermal equilibrium in the early universe after inflation, since otherwise they would have over-closed the universe. This sets a relatively low upper bound on the reheating temperature after inflation in that scenario. Finally, asymmetric dark matter production might allow to avoid complete wash-out of the non-thermal abundance through thermalization after the EW phase transition. The paper is organized as follows: in Section 2 we review the formalism that describes particle production at the end of the EW phase transition for the case of very elastic bubble collisions [7,8] and extend it to the case of very inelastic ones, highlighting the differences between both scenarios [9]. Then, in Section 3 we explicitly compute the particle production efficiency of scalar, fermion, and vector boson particles coupled to the Higgs (either directly or via an effective Higgs portal). In Sections 4 and 5 we focus on dark matter production at the end of the EW phase transition. First we discuss in Section 4 the conditions for non-thermally produced dark matter to avoid subsequent wash-out and constitute the bulk of the present dark matter density, selecting heavy (multi-TeV) vector boson dark matter as a viable example. We go on to analyze in detail non-thermal dark matter production in that scenario and the subsequent evolution of the non-thermally generated abundance after the EW phase transition, including finally a discussion on the current XENON100 bounds and direct detection prospects. Then, in Section 5 we study the non-thermal production of very heavy (M X ≫ v EW ) vector boson dark matter, and outline the conditions under which these baby-zillas constitute a viable dark matter candidate. In the case of asymmetric non-thermal dark matter production, we find it difficult to avoid subsequent wash-out, and the discussion is left for an appendix. We summarize our results and conclude in Section 6.
2 Particle Production at the EW Phase Transition

Bubble Collisions in the EW Phase Transition
If the early Universe was hotter than T EW ∼ 100 GeV it must have undergone an EW phase transition at some point in its history. Within the SM, the EW phase transition is a smooth cross-over, however it is conceivable that new degrees of freedom beyond the SM modify the Higgs potential so as to make the transition first order. This is what we assume throughout this paper, without specifying the full theory that makes the first order transition possible. In that case, the EW phase transition proceeded through nucleation and expansion of bubbles of true Higgs vacuum, which eventually collided among each other completing the transition. As this was happening during the radiation dominated era, the bubble expansion process would then have taken place in a thermal environment (except for the case when a period of inflation would have preceded the phase transition).
For a first order phase transition occuring in a thermal environment, the study of the bubble expansion process reveals that the thermal plasma exerts some amount of friction on the expanding bubble wall, and this friction tends to balance the pressure difference on the bubble wall driving the expansion. In the usual picture, nucleated bubbles reach a stationary state after a very short period of acceleration, with a constant wall velocity depending specifically on the interactions of the bubble wall with the degrees of freedom in the plasma [10,11] and on the resulting fluid dynamics [12][13][14] (see [15] for a review). In this case, the amount of energy stored in the bubble walls at the time of the bubble collisions is negligible compared to the available energy of the transition, since most of this available energy gets converted into plasma bulk motion and thermal energy [16].
However, this picture was recently challenged in [17], where it was shown that the friction exerted by the plasma may saturate to a finite value for ultrarelativistic bubble walls. As a consequence the stationary state assumption will no longer be true when the pressure difference on the bubble wall is larger than the friction saturation value, which may happen for strongly first order phase transitions. In this scenario, if there are no hydrodynamic obstacles that prohibit the bubble walls to become highly relativistic in the first place (see however [18]), bubbles will expand in an accelerated way ('the so-called runaway bubbles), with almost all the energy of the transition being used to accelerate the bubble walls 1 [15]. By the end of the phase transition (when bubbles start colliding), these runaway bubbles may reach very large values of γ w : with v T the value of the Higgs VEV in the broken phase and β −1 ∼ (10 −3 − 10 −2 ) H −1 being the duration of the phase transition [19]. The estimate (2.1) follows from balancing the surface energy on the bubble wall (2.5) and the energy available inside the bubble. Once bubbles start colliding, the energy stored on the bubble walls will be liberated into the plasma. As argued above, for "runaway" bubbles this will correspond to a very large portion of the energy budget of the phase transition, and therefore this process can be very important. Under certain circumstances, this may also hold true for highly relativistic bubble walls (γ w ≫ 1) that reach a stationary state long before bubble collisions start (meaning that γ w ≪ γ max w ), in which case the amount of energy stored in the bubble walls will be very small compared to the available energy of the transition, but still important when released into the plasma at the end of the transition.
The process of bubble collisions in cosmological first order phase transitions is by itself a very complicated one. Consider a configuration of two planar bubble walls 2 initially far away from each other, that approach and collide [7,8,20]. Depending of the shape of the potential for the scalar field φ driving the transition (in our case, the Higgs field h), the bubble collision will be approximately elastic or partially inelastic [7,8] (see also [9]). In the first case, the walls reflect off one another after the collision, which reestablishes a region of symmetric phase between the bubble walls. For a perfectly elastic collision the field profile of the colliding walls in the limit of infinitely thin bubble walls (taken as step-functions) can be written as [8] where v w is the bubble wall velocity, the bubble walls move in the z-direction and the collision is assumed to happen at t = 0. Since we are ultimately interested in scenarios where γ w ≫ 1, we will take the ultrarelativistic limit v w → 1 in the rest of the section. The field profile (2.2) neglects the thickness of the bubble walls l w (generically, l w ∼ (10 − 30)/T EW , with T EW ∼ 50 − 100 GeV). To capture the wall thickness effects one can consider another ansatz for the profile of the colliding bubble walls: ( Figure 1 -Left), where it will perform small-amplitude oscillations. In this case the collision is approximately elastic as described above, with the bubble walls being effectively reflected as a region of symmetric phase is re-established between them. The walls move then away from each other until vacuum pressure makes them approach and collide again, repeating the process several times. In each collision some fraction of the energy stored in the walls is radiated into scalar waves and quanta of the fields coupled to h, until all of the energy in the walls is radiated away. In contrast to this scenario, for a potential V (h) with very non-degenerate minima ( Figure 1 -Right), the field oscillation after the "kick" in the region close to the collision point does not effectively drive the field over the potential barrier. As a consequence, the field stays in the basin of attraction of the broken minimum v T and performs relatively large-amplitude oscillations around it, giving rise to a large amount of energy radiated into scalar waves (as opposed to the previous scenario). In this case the collision is very inelastic.
Following [9], we compute the numerical solution for the field profile h(z, t) corresponding to the collision of two bubble walls, obtained from solving (2.4) with a toy potential V (h) of the form both in the case of nearly degenerate minima (Figure 1 -Left) and very non-degenerate minima (Figure 1 -Right). The results are shown in Figure 2 (similar plots appeared earlier in [21]). can obtain an analytic solution h(z, t) = h TI for the case of a "totally inelastic collision" (as opposed to the "perfectly elastic collision" described earlier), in which all the energy is radiated in the form of scalar waves after the bubble collision. For t < 0 (before the collision) we have which matches h lw (z, t < 0). In order to obtain h TI (z, t) for t > 0, we note that the field will not leave the basin of attraction of the broken minimum v T after the collision. We can then approximate the potential V (h) the field will feel for t > 0 as This allows to solve the equation of motion (2.4) explicitly for δh TI (z, t) ≡ h TI (z, t) − v T : where the boundary conditions follow from imposing continuity of δh TI (z, t) and ∂ t δh TI (z, t) at t = 0. From (2.9), we finally obtain Sinh π lw pz 2 γw Sin Notice that in the limit m h → 0, (2.9) becomes (∂ 2 t − ∂ 2 z ) δh TI (z, t) = 0 and (2.7) is also a solution for t > 0, case in which the two bubble walls would pass through each other without actually colliding.
The analysis for the dynamics of bubble collisions presented here may be extended to phase transitions involving multiple fields (see for example [20]), although in this case the analysis of the field evolution after the bubble collision becomes much more complicated (since the scalar potential is multidimensional and the field "excursion" at the moment of the bubble collision will involve several fields), and we will not attempt it here.

Particle Production Through Bubble Collisions
The bubble collision processes analyzed in the previous section allow to liberate into the plasma the energy contained in the bubble walls. This can happen either via direct particle production in the collisions or via radiation of classical scalar waves which will subsequently decay into particles. For bubble collisions taking place in a thermal environment, the number densities n α for the different particle species created during the collisions should very quickly approach the ones in thermal equilibrium n EQ α after the phase transition, thus rendering the particle production process irrelevant for the subsequent evolution of the Universe. However, as it has been briefly discussed in the introduction, under certain conditions fast thermalization of certain species after the phase transition may be avoided, which can make the particle production process very important in that case.
In order to study the particle production through bubble collisions, we will treat the scalar field configuration h(z, t) as a classical external field and the states coupled to it as quantum fields in the presence of this source. In doing so, we will neglect the back-reaction of particle production on the evolution of the bubble walls themselves throughout the collision, which should be a good approximation when the energy of the produced particles (for each species) is much less than the energy contained in the field configuration h(z, t). The probability of particle production is given by [8] is the generating functional of 1PI Green functions, and to the quadratic order in h being the 2-point 1PI Green function. In terms of its Fourier transformΓ (2) (p 2 ), and using (2.11) and (2.12) we get The last integral in (2.13) is just h (p) 2 , withh(p) being the Fourier transform of the Higgs field configuration h(x)h For a background field configuration h(z, t), its Fourier transform is given byh(p) = (2π) 2 δ(p x ) δ(p y )h(p z , ω). Then, using (2.13), we obtain the mean number of particles produced per unit area [8]: The physical interpretation of (2.15) is rather simple [8]: the scalar field configuration h(z, t), corresponding to the two bubble walls that approach and collide, can be decomposed into modes of definite four-momentum p 2 = ω 2 − p 2 z via the Fourier transform. Modes with p 2 > 0 represent propagating field quanta with mass squared m 2 = p 2 . Then, (2.15) integrates over the amount of field quanta of mass p 2 contained in the field configuration multiplied by the probability of those quanta to decay.
The Fourier transform of the background field configuration h(z, t) can be performed explicitly both for the case of a perfectly elastic collision and of a totally inelastic one analyzed in the previous section. For a perfectly elastic collision, in the limit of infinitely thin walls (h(z, t) = h ∞ ), we obtaiñ However, since the highest values of p z and ω available in the field configuration are naively expected to be of order γ w /l w (modes with p z , ω ≫ γ w /l w will be exponentially damped), the integration in (2.15) should in this case be cut-off for p z > γ w /l w and ω > γ w /l w . From (2.15) and (2.16) we then obtain Alternatively, when the thickness of the bubble walls is accounted for (h(z, t) = h lw ), the Fourier transform of (2.3) gives which automatically incorporates the exponential damping for ω, p z ≫ γ w /l w . The mean number of particles per unit area now reads For the opposite case of a totally inelastic collision (h(z, t) = h TI ), the Fourier transform is given bỹ The relative "−" sign between the two contributions in (2.20) can be easily understood noticing that in the limit m h → 0 the Fourier transform of h TI (z, t) should giveh(p z , ω) ∼ δ(ω ± p z ). From (2.20), the mean number of particles produced per unit area in the case of a totally inelastic collision is given by The expressions (2.17), (2.19) and (2.21) can be rewritten in a more compact form by making the change of variables χ = ω 2 − p 2 z , Ψ = ω 2 + p 2 z . After performing the integral in Ψ, the mean number of particles produced per unit area finally reads The function f (χ) encodes the details of the bubble collision process and quantifies the efficiency of particle production. For a perfectly elastic collision, in the limit of infinitely thin bubble walls, we have For a perfectly elastic collision, and for bubble walls with finite thickness, we have Finally, for a totally inelastic collision, we have In Figure 3 we compare the efficiency f (χ) for the various cases (2.23), (2.24) and (2.25). Notice that f TI (χ) diverges as χ → m 2 h . This divergence is artificial, due to considering h(z, t) over infinite time and space, and should be cut-off since our solution is not valid over distances larger than the bubble radius R B . Implementing this cut-off can be well approximated by replacing in (2.24) Defining χ min as the minimum value of χ for which particle production is possible (corresponding to the squared sum of the masses M α of the particles being produced), we immediately see from Figure 3 that for a totally inelastic collision, production of light particles (χ min < m 2 h ) may be very efficient, while production of heavy particles (χ min ≫ m 2 h ) will be extremely suppressed. For a perfectly elastic collision, however, the production of heavy particles may be relatively efficient (we will comment further on this point at the end of section 3). For the study of the efficiency of particle production in varios different scenarios in the next sections, we will use (2.23) for the case of an elastic collision, while for the case of a very inelastic one it is possible to show that (2.25) (together with (2.26)) can be approximated as Let us now turn to the evaluation of the imaginary part of the 2-point 1PI Green function's Fourier transformΓ (2) . Through the optical theorem, we can write: where M(h → α) 2 is the spin-averaged squared amplitude for the decay of h into a set of particles α with masses M α , χ min ≡ ( M α ) 2 is the minimum value of χ for which this decay is possible and dΠ α is the relativistically invariant n-body phase space element Then, the number of particles of a certain type α produced per unit area during the bubble collision will simply read from (2.22) and (2.28) The amount of energy produced per unit area in the form of particles α is obtained by weighting (2.30) by the energy of each decaying Fourier mode. This yields From (2.30) and (2.31), the non-thermally produced energy density ρ α (assuming that the produced particles quickly diffuse into the bubble interior) reads with A ∼ 4 π R 2 B being the total collision area and V the volume of the two colliding bubbles. From (2.32), and bearing in mind that R B ≃ β −1 , the non-thermally generated comoving energy density is with s(T EW ) the entropy density after the EW phase transition.

Particle Production via the Higgs Portal
The efficiency of particle production may strongly depend on the nature of the particles being produced. In this section we will analyze the particle production efficiency for scalars S, fermions f and vector bosons V µ coupled to the Higgs field. Apart from estimating the production of SM fermions and gauge bosons through this process, we will consider a simple Higgs-portal extension of the SM in order to study the production of other possible scalar, fermion or vector boson particles. Furthermore, we will restrict ourselves to Z 2 symmetric Higgs-portal scenarios, since we will ultimately be interested in dark matter analyses. We also comment on how to interpret the results in the case when the calculated particle production exceeds the energy available in the bubble wall.

Scalars
For the complex scalar S interacting with the SM via the Higgs portal, the relevant part of the lagrangian is given by (3.34) In this case, M(h → SS) 2 = λ 2 s v 2 T , and one immediately obtains  From Figure 4 it can be clearly seen that scalar particle production is quite suppressed for elastic collisions. For very inelastic collisions, heavy-scalar particle production is extremely suppressed, while production of light scalars turns out to be very efficient in this case. In fact, Figure 4 shows that for large values of λ s (λ s 1) the naively calculated energy of the produced particles E exceeds the amount of energy on the bubble walls E w . That inconsistency indicates that in these cases backreaction cannot be neglected. We will comment and expand on this issue in section 3.4.

Fermions
Turning now to fermionic particle production, in the presence of a tree-level Yukawa coupling between the Higgs and the fermions λ f H f f , the squared decay amplitude reads which, in the case of SM fermions, leads directly to The production of (SM) fermions will then be enhanced with respect to the one of Higgsportal S−scalars (specially in the limit of very elastic collisions, see Figure 5) due to the extra factor χ − 4 m 2 f in (3.37). Scenarios where the fermionic particle production might be important include (apart from the SM itself) the MSSM and its various extensions, due to the tree-level coupling between Higgses, Higgsinos and Gauginos 4 .
In the absence of a direct coupling, the interaction between the Higgs and the fermions will occur via an effective operator. This is the case for the so-called fermionic Higgs-portal: However, since bubble collisions may excite very massive Higgs field modes (p 2 ≫ T 2 EW ), particle production in this case may be sensitive to the UV completion of the Higgs-portal effective theory, making it unreliable to compute the particle production in the fermionic Higgs-portal via (3.38). Here we consider a simple UV completion for the fermionic Higgsportal, and compute the particle production in this case. We add a singlet scalar field S as a mediator between the Higgs field and the fermion f , the relevant part of the lagrangian being For simplicity, we will avoid a vev for S (it can be done through the addition of a linear term for S in (3.39)). For µ s = 0 the effective fermionic Higgs-portal operator |H| 2 f f will be generated at tree-level. The squared decay amplitude for h → f f will then be leading finally to When µ s = 0 the effective fermionic Higgs-portal operator is not generated at tree-level, but rather the decay h → f f occurs via a finite 1-loop diagram, yielding where F m 2 f , M 2 s , χ is a form factor that scales as Fermionic Higgs-portal particle production both in the µ s = 0 and µ s = 0 is shown in Figure 5, where it can be clearly seen that the production in the absence of a direct coupling between the Higgs and the fermions f differs from what would have been naively obtained using (3.38). As for the case of scalar particle production, under certain circumstances the estimate of fermionic particle production neglecting backreaction exceeds the amount of energy stored in the bubble walls (E > E w ), and in order to obtain a physically meaningful result backreaction should be included (We will expand on this issue in section 3.4).

Vector Bosons
Finally, we study the production of vector boson particles. In the presence of a tree-level coupling between the Higgs and the vector bosons λ V M V h V µ V µ , the squared decay amplitude reads Comparing (3.35), (3.37) and (3.46) we immediately observe the relative efficiency of particle production for scalars, fermions and vector bosons. While Im [Γ (2) (χ)] scales as χ 0 for scalars, and as χ for fermions, in the case of vector bosons it scales as χ 2 , thus greatly enhancing production of vector bosons with respect to scalars or fermions for very elastic collisions (see Figure 6). It is then expected that most of the available energy from the EW phase transition will go into W µ and Z µ gauge boson production and (possibly) other vector bosons coupled at tree-level to the Higgs in extensions of the SM 5 .
In the absence of a direct coupling, the interaction between the Higgs and the vector bosons may occur via an effective operator, as in the so-called vector Higgs-portal [23]: However (like for the fermionic Higgs-portal) an analysis of vector boson particle production in the context of the effective theory (3.47) will be unreliable due to very massive Higgs field modes (p 2 ≫ T 2 EW ) being excited during the bubble collisions. Vector boson particle production will then be sensitive to the way in which the effective operator |H| 2 V µ V µ is generated. One possible way of generating the effective operator at tree-level, being V µ a hidden U(1) gauge field, is by integrating out a U(1)-charged complex scalar S which has a Higgs portal coupling |H| 2 S * S, the relevant part of the lagrangian then being In this scenario, the vector boson V µ acquires a mass via the spontaneous breaking of the hidden U(1), through a vev v S for the S-scalar 6 . The squared decay amplitude for h → V µ V µ will then be with Γ s being the decay width of S. This leads to Vector boson effective Higgs-portal particle production is shown in Figure 6, resulting in a very suppressed particle production with respect to the case in which the vector bosons and the Higgs couple directly at tree-level, specially for very elastic collisions. From Figure 6 it is also clear that backreaction is most important for direct vector boson particle production (for which the production estimate yields E ≫ E w ).

Backreaction and Relative Efficiency
Clearly, for the present analysis of particle production to be physically meaningful it must be assumed that the total energy of the produced particles is less than the energy contained in the background field configuration h(z, t). Moreover, when the energy of the produced particles starts being comparable to the energy of the background field we expect backreaction on h(z, t) due to the particle production to be important. Then, in order for the previous analysis to be reliable, it is needed As it has been shown in the previous section, for fermion or vector boson particle production the previous condition (3.51) is not satisfied, and in some cases even E ≫ E w is obtained (Figure 6 LEFT), signaling the extreme importance of backreaction in those scenarios.
Since incorporating backreaction into the present analysis of particle production is extremely difficult and lies beyond the scope of this paper, we simply note that the relative efficiency in particle production for the different species in the present analysis should be roughly correct even when backreaction is important. Then, an estimate of the particle production in cases where some of the species are very efficiently produced may be obtained just by normalizing the production to the total energy in the bubble walls. For very elastic bubble collisions, it has been shown in section 3.3 that production of W µ and Z µ gauge bosons is extremely efficient, which will then leave very little energy left in the bubble walls for producing other particle species. The relative efficiencies (defined as ratios of energy in produced particles) of the different species for a perfectly elastic collision, normalized to the energy contained in the bubble walls (assuming that most of the available energy goes into producing W µ and Z µ ) is shown in Figure 7. A good estimate of the non-thermally generated comoving energy density (per particle species α) in this case may then be given by The fact that this is a reliable estimate of the particle production efficiency for the case of very elastic collisions is due to the high-p 2 modes of the bubble wall carrying almost all the energy of the bubble wall. The energy carried by the high-p 2 modes will then mostly go into vector boson production (their production efficiency at high p 2 is much larger than fermionic or scalar ones), result that holds even without incorporating backreaction into the analysis.
On the other hand, for very inelastic collisions the results from the previous section show that particle production is only effective for light particles (M X m h /2). Therefore, production of W µ and Z µ will be very suppressed in this case, along with any other heavy particle, and most of the available energy will go into production of SM fermions (mainly bottom quarks) and (possibly) new light scalars or fermions with sizable couplings to the Higgs.

Non-thermal Multi-TeV WIMP Dark Matter
In this section we focus on the case of relatively heavy dark matter, M X TeV, and explore the conditions under which the amount of non-thermally produced heavy dark matter can end-up accounting for a sizable part of the observed dark matter relic density (dark matter may nevertheless still have a thermal component coming from the usual freeze-out process). The first condition is clearly that bubble collisions have to be fairly elastic: it has been shown in sections 3 that for very inelastic bubble collisions only light (M X m h /2) particles are efficiently produced, while heavy particle production is extremely suppressed. Since fast thermalization of light species after the EW phase transition seems unavoidable 1 , for very inelastic bubble collisions either dark matter is not efficiently produced or it thermalizes immediately after the end of the EW phase transition, not having any influence on the subsequent evolution of the Universe.
For very elastic bubble collisions, the analysis from sections 3 and 3.4 shows that electroweak gauge bosons W µ and Z µ are most efficiently produced, and the relative production efficiency of heavy fermions and scalars is too low (for them to be able to account for a sizable part of the observed dark matter relic abundance, see Figure 7). This leaves heavy vector bosons with a direct coupling to the Higgs field as the only viable candidate for non-thermally produced dark matter during the EW phase transition.
In the following we perform an analysis of heavy vector boson dark matter coupled to the Higgs, including an overview of thermal freeze-out and direct detection constrains from XENON100 [22] (see [23,24] for more details), and a comparison between the amount of non-thermally produced dark matter and the amount of dark matter produced through thermal freeze-out. We also study the evolution of the non-thermally produced dark matter component after the EW phase transition.

Higgs-Vector Dark Matter Interplay
Consider a vector boson dark matter candidate with mass M V and a tree-level coupling to the Higgs [23,25], This coupling mediates the dark matter annihilation into Standard Model particles, as well as the elastic scattering on nucleons relevant for dark matter direct detection. Concerning the former process, the Higgs boson can mediate annihilation of dark matter into electroweak gauge bosons (for heavy dark matter they are the most important annihilation channel) through the couplings The spin-averaged amplitude squared for the annihilation process V µ V µ → W + µ W − µ in the limit s ≫ m 2 h is given by Given (4.55), the thermally averaged annihilation cross section is given by where z = M V /T , and K 1 (z), K 2 (z) are Bessel functions. For z ≫ 1, (4.56) reduces to The thermal cross section giving rise to the observed value of the relic density σv WMAP ≈ 2.6 · 10 −9 GeV −2 corresponds, for heavy dark matter M V ≫ m h and using (4.57), to Turning now to dark matter direct detection, the spin-averaged amplitude squared for Higgs-mediated dark matter elastic scattering on nucleons reads Here, m N ≈ 0.939 GeV is the proton/neutron mass and f N is the effective Yukawa coupling of the Higgs to nucleons which, following [24], we take f N = 0.326 based on the lattice estimate in [26]. In the last step we have taken the limit t ≪ m 2 N , m 2 h , M 2 V . The elastic scattering cross section then reads On the other hand, the XENON100 bound on the dark matter elastic scattering cross section on nucleons, for M V TeV is approximately Therefore, (4.58) and (4.61) leave a sizable window in the parameter space (M V , λ V ) for which the dark matter abundance obtained via thermal freeze-out is significantly smaller than the observed dark matter relic density, and still the value of λ V (as a function of M V ) is below the XENON100 bound (as shown in Figure 8).

Fate of Non-Thermally Produced Vector Dark Matter
Given the results from the previous section (summarized in Figure 8), it is fair to ask if, in the region of (M V , λ V ) parameter space in which the thermal component is not enough to account for the observed dark matter relic density, dark matter produced non-thermally at the EW phase transition could account for the extra needed amount. Using the results from production efficiency of heavy vector boson dark matter obtained in sections 3.3 and 3.4, we show in Figure 9 the value of λ V (as a function of M V ) for which the amount of non-thermal vector boson production equals the observed dark matter relic density (dashed-blue line). Then, for values of λ V above the thermal cross section giving the observed relic density σ v WMAP , non-thermal production of heavy vector bosons is so efficient as to generate amounts of dark matter much larger than the observed dark matter relic density.
Assuming that at the time of the EW phase transition vector boson dark matter is already frozen-out (T fo ≃ M V /20 > T EW ), we can study the evolution of the non-thermally generated dark matter abundance via a simple Boltzmann equation in which the comoving dark matter number density Y fulfills Y (z) ≫ Y EQ (z) (with Y EQ (z) being the equilibrium comoving number density), yielding with α = (4π 2 √ ξ g * )/45 ≃ 2.642 (g * ∼ 100 being the number of relativistic degrees of freedom in the thermal plasma and ξ ≡ 90/(32 π 3 )), M Pl = 1.2 × 10 19 GeV and y(z) = α σ v M Pl M V Y (z). Integration of (4.62) for z > z EW yields Then, given the fact that non-thermal vector boson dark matter production is much larger than the observed relic density in the (M V , λ V ) region of interest, we can take the limit y(∞) ≪ y(z EW ), obtaining From (4.64), we immediately obtain that the value of the annihilation cross section that will yield the observed dark matter relic density once the non-thermally generated dark matter evolves after the EW phase transition is simply given by The red lines in Figure 9 show the values of λ V yielding the correct "non-thermal" annihilation cross section (4.65) for several values of T EW .
This analysis shows that non-thermal production of multi-TeV vector boson dark matter at the EW phase transition (in (M V , λ V ) parameter space in which the amount of dark matter yielded by thermal freeze-out is not enough to account for the observed dark matter relic density) is efficient as to generate a dark matter amount much larger than the observed relic density. This results in a reactivation of thermalization processes that lead to partial wash-out of the non-thermally generated dark matter (wash-out is not complete due to the reactivation happening for T < T EW < T fo ), meaning that multi-TeV dark matter may have a thermal spectrum despite a large fraction of it having been produced non-thermally at the EW phase transition. As shown in Figure 9, in the presence of these non-thermally produced WIMPs, the relation between mass and coupling giving rise to the observed dark matter relic density gets modified with respect to the usual thermal freeze-out scenario, leading to better detection prospects in the multi-TeV region for future dark matter direct detection experiments.

Baby-zillas: Super-Heavy Dark Matter from the EW Phase Transition
In this section we study the production of super-heavy dark matter with a mass M X satisfying M GUT ≫ M X ≫ v EW in the bubble collisions at the end of a very strong EW phase transition. We call these dark matter particles baby-zillas because of many similarities (but smaller mass) to the WIMP-zilla scenario. From Figure 7, it can be inferred that for γ w ∼ 10 14 − 10 15 non-thermal heavy vector boson production in elastic bubble collisions can be so efficient as to generate the observed dark matter relic density even for very large dark matter masses M V ∼ 10 6 − 10 8 GeV and perturbative values of the coupling λ V . Using (3.52), we plot in Figure 10 the region in parameter space (M V , λ V ) for which non-thermal V µ production directly yields the observed dark matter relic density.

Bounds on the Reheating Temperature After Inflation
A stable particle with mass M V ∼ 10 5 − 10 8 GeV would yield a much larger relic abundance than the observed DM relic density. were it in thermal equilibrium at some stage after inflation. For such a massive species, the annihilation cross is always smaller than the one needed to yield the observed DM relic density through thermal freeze-out. It is then needed that this particle species never reached thermal equilibrium after the end of inflation. This sets an upper bound on the reheating temperature after inflation, specifically T RH < T fo (with T fo being the temperature below which the particle is decoupled from the thermal plasma). For a heavy vector boson V µ annihilating into SU(2) gauge bosons (the most important annihilation channel in this case) through the Higgs, T fo satifies where the thermally averaged annihilation cross section σv is given by (4.56). In Figure  11 we plot the minimum value of z (corresponding to the maximum allowed value of the reheating temperature T RH ) as a function of the mass M V for the range of λ V values giving rise to the observed dark matter relic abundance for γ w = 10 14 −10 15 (see Figure 10). We see that the upper bound on T RH is relatively insensitive to the precise value of γ w , and roughly scales as T max RH ∼ M V /10.  Figure 11: Bounds on the Reheating temperature after inflation for the requirement that dark matter never reaches thermal equilibrium after inflation, namely T RH ≤ T fo , as a function of the dark matter mass M V , and assuming λ V (M V ) for which non-thermal production yields the observed relic abundance (as shown in Figure 10).

Conclusions
Dark matter may have been efficiently produced at the end of a first order EW phase transition if it has a large coupling to the Higgs field. In this paper we investigated the conditions for this non-thermal production mechanism to account for most of dark matter in the Universe. We considered scalar, fermion and vector dark matter coupled to the SM through the Higgs (either via a direct, tree-level interaction or an effective Higgs-portal coupling), and found that production of vector bosons directly coupled to the Higgs is most efficient, while for scalar and fermions most of the energy stored in the bubble walls is bound to be released into production of SM particles. This analysis singles out vector dark matter in the present context. For very inelastic bubble collisions only dark with M X 100 GeV can be efficiently produced, while production of heavier dark matter is extremely suppressed. Unfortunately, for a dark matter mass in this range, we did not find a way to avoid subsequent thermalization and the wash-out of the non-thermal component, and therefore in this case dark matter production at the EW phase transition is irrelevant. The situation is quite different for highly elastic bubble collisions. In that case, dark matter with M X ≫ 100 GeV can be efficiently produced for the so-called runaway bubbles, that expand with a very large γfactor.
We have identified two scenarios where wash-out of dark matter produced at the EW phase transition can be naturally avoided. One has dark matter in the multi-TeV range, which makes it possible for non-thermally produced dark matter to remain out of thermal equilibrium after the EW phase transition. We determined the region in the parameter space of dark matter mass and coupling to the Higgs where the correct relic abundance is reproduced. For a given mass, the coupling has to be larger than in the usual thermal freezeout scenario for Higgs portal dark matter, which can be especially relevant for direct detection searches, as it opens the possibility of detecting a signal from multi-TeV non-thermal dark matter in the near future by XENON100 and LUX experiments. The other scenario is baby-zilla dark matter with M X ∼ 10 6 -10 8 GeV. Surprisingly enough, such super-heavy dark matter can be produced in important quantities at the end of a strongly first-order EW phase transition, provided the dark matter coupling to the Higgs is large, and the γ factor of the bubble walls is near its maximal value of γ w ∼ 10 15 . In order for the baby-zillas to be a viable dark matter candidate, they must have never reached thermal equilibrium, which then constrains the reheating temperature after inflation in this scenario.
For multi-component dark matter (X = X α ), an asymmetry in the number densities of X α and X α may be generated during the particle production. We will analyze in detail below the generation of this asymmetry for scalars (X α = S α ). Then, in section A.2 we study the evolution of the generated asymmetries after the EW phase transition.

A.1 Decay Asymmetries: Producing a Dark Matter Asymmetry
Let us consider a set of N i real scalars h i (that includes the field(s) involved in the EW phase transition) coupled to a set of N α complex scalars S α via a trilinear interaction. The relevant part of the lagrangian is where by hermiticity C iαβ = C * iβα (it follows that C iαα are real, but C iαβ with α = β can be complex), and the mass matrix for the scalars S α is taken to be diagonal without loss of generality. We also consider a possible term µ ij h i h 2 j appearing in V (h i ). The lagrangian (A.1) incorporates a Z 2 symmetry that makes the lightest of the scalars S α stable, which may then be a suitable dark matter candidate. In order for an asymmetry in the production of S α and S * α to be generated, we need a nonzero value for At tree level 3) and there is no asymmetry generated. At 1-loop we include the 1PI diagrams shown in Figure 12. Their contribution to the 1-loop decay amplitude is where the integrals I T andĨ T correspond to and can be computed in terms of the usual Passarino-Veltman 3-point scalar loop integral C 0 . The leading order difference between |M(i → α * β)| 2 and |M(i → α β * )| 2 is due to the interference between the tree level and 1-loop decay amplitudes. We obtain |M(h i → S * α S β )| 2 − M(h i → S * β S α ) 2 = 1 4 π 2 j,γ,δ µ ij Im C * iαβ C jαδ C jδβ Im Ĩ T + Im C * iαβ C iγδ C jαγ C jδβ Im [I T ] (A.6) An interplay between "weak" and "strong" phases (complex couplings and imaginary part of a 1-loop integral due to particles in the loop going on-shell) is then needed to get CP violation in the decay. For N i = 1 (i = j = 1) both terms on the right-hand side of (A.6) will vanish for N α < 3 but may be zonzero for N α ≥ 3. For the case of two fields h i (N i = 2), already with two scalars (N α ≥ 2) it is possible to obtain an asymmetry. Figure 12: Tree-level and 1PI 1-loop contributions to the decay h i → S * α S β To obtain the total combined production of S α and S * α particles, we will just consider the tree level contribution to h i → S * α S α , h i → S * α S β and h i → S * β S α . We then get Im Γ (2) (χ) α = |C iαα | 2 dΠ αα + |C iαβ | 2 dΠ αβ (A.7) with For the asymmetry in the production of S α and S * α particles we obtain Im Γ (2) (χ) with |M(i → α * β)| 2 − |M(i → α β * )| 2 given by (A.6).
leading to Ξ α → ∆ α . However, the various processes entering (A.13) will tend to erase the asymmetries ∆ α . In particular, the process X α + SM → X β + SM (responsible for kinetic equilibrium among the different species X α ) and the decay process X α → X β + SM, will, if active, wash-out the asymmetries very rapidly, being also quite insensitive to the temperature of the Universe after the EW phase transition (the annihilation process X α + X * β → SM also washes-out the asymmetries, but is suppressed for low T EW ). If these processes are in equilibrium after the EW phase transition, they will very rapidly drive which, together with (A.11), leads to ∆ α → 0. While it is possible for the decay process X α → X β + SM to be suppressed if the different species are quite degenerate, this automatically results in unsuppressed kinetic equilibrium. This seems to rule-out asymmetric production as a viable mechanism to avoid complete wash-out of the amount of light dark matter nonthermally produced during the EW phase transition, the reason being that it is not possible (at least in this simple scenario) to keep all the processes that tend to erase the asymmetries ∆ α out of equilibrium after the EW phase transition, while having an efficient particle production at the transition itself.