Gravitational wave signals of dark matter freeze-out

We study the stochastic background of gravitational waves which accompany the sudden freeze-out of dark matter triggered by a cosmological first order phase transition that endows dark matter with mass. We consider models that produce the measured dark matter relic abundance via (1) bubble filtering, and (2) inflation and reheating, and show that gravitational waves from these mechanisms are detectable at future interferometers.


Introduction
The identity of dark matter (DM) and its production mechanism are among the most important open questions in physics. In the weakly interacting massive particle (WIMP) paradigm, with its thermal freeze-out mechanism, the measured DM relic density requires a WIMP mass of O(10-10 3 ) GeV, and an electroweak-scale DM annihilation cross section. The decoupling temperature T dec is related to the DM mass m χ by T dec m χ /24. The vanilla version of this scenario had been challenged by the non-observation of WIMPs in DM direct detection searches. Alternative scenarios for DM production in the early universe often assume a DM sector that is out of thermal equilibrium with the standard model (SM) sector. For example, DM may be produced in the decays of a heavy particle [1,2]. DM may also be produced by freeze-in through the feeble annihilation of particles which are thermalized with the SM bath [3][4][5]. However, in all of the above scenarios, the DM mass is constant during DM production.
The discovery of the 125 GeV SM-like Higgs boson h at the Large Hadron Collider (LHC) [6,7] consolidates spontaneous symmetry breaking as the mechanism that gives the SM particles their mass. The Higgs mechanism gives the simple relation m f = y f · v EW between the fermion mass m f and its Yukawa coupling to the Higgs boson y f , where h ≡ v EW 246 GeV is the vacuum expectation value (VEV) of the SM Higgs. A picture of the universe going through an electroweak phase transition because finite temperature effects modify its scalar potential as the universe cools down, emerges. Before the phase transition, when all the SM particles are massless, the global minimum of the scalar potential is located at h = 0. After the phase transition, the global minima of the potential shift to nontrivial values h = 0, which gives mass to the SM particles. The SM Higgs quartic potential predicts a second order phase transition. However, since we do not fully understand the entire structure of the scalar potential of the 125 GeV Higgs boson, and since the existence of additional scalars is a possibility, it is unknown whether the early universe experienced a first or second order phase transition The DM mass may be generated by a similar mechanism [8][9][10][11][12]. The mass originates from its couplings to a scalar, which obtains a non-trivial VEV in the early universe, so that massless DM becomes massive during the phase transition. The scalar may or may not be the 125 GeV Higgs boson. We consider a first order phase transition (FOPT) in the early universe, with vacuum bubbles nucleated at temperature T , which ends with the expanding bubbles populating the entire universe. The symmetric and broken phases are located outside and inside the bubbles, respectively. The massless DM particles outside the bubbles become massive when they enter the bubbles. Only massless DM particles that carry kinematic energy larger than m χ can penetrate the bubble walls and become massive. DM inside the bubbles abruptly decouples from the thermal bath if T < T dec . The result is that the bubbles filter out a certain amount of DM and determine the DM relic abundance [11,12]. The massless DM outside the bubbles remains thermalized with SM radiation. It is also possible that all the massless DM particles enter the bubbles after being diluted by a period of inflation, which determines the relic abundance [9]. DM particles with insufficient kinetic energy to enter the bubbles, bounce back to the symmetric phase and slow down the bubble expansion by applying pressure on the bubble walls.
The value of m χ /T needed to produce the correct DM relic abundance depends on the velocity of the bubble walls v w . For instance, T m χ /30 for m χ = 1 TeV and v w = 0.01, which satisfies T < T dec . Note that DM freeze-out induced by a FOPT can easily accommodate DM masses above a PeV, which is beyond the current sensitivities of DM direct detection and LHC searches.
In this paper, we focus on gravitational wave (GW) signals of sudden DM freeze-out caused by a FOPT during which DM mass is generated. Because the power and frequency spectrum of the GW signal is model dependent, we choose two example models, i) Scalar Quartic Model [13][14][15][16] and ii) SU (2) X Model [8,9], to demonstrate that in parameter space regions that yield the observed DM relic abundance, a detection is possible at future GW interferometers. In the Scalar Quartic Model, the DM abundance is determined by bubble filtering, while in the SU (2) X Model, the DM abundance is set by inflation and reheating.
The paper is organized as follows. Bubble filtering is described in section 2, and computations of the bubble wall velocity are detailed in section 3. In section 4, we list the contributions to GW spectra from various processes. We calculate the GW signals for the two example models in section 5, and summarize in section 6.

Bubble filtering
During the FOPT and bubble expansion, massless (massive) DM particles are located outside (inside) the bubble, and momentum conservation must be satisfied at the bubble wall. An incident DM particle enters the bubble if it carries kinetic energy larger than its mass inside the bubble. Otherwise, the massless DM particle is reflected and stays outside the bubble. If a thermal flux of χ is incident on the wall, the number density of DM particles that enter the bubble is [12] n in χ = n in where γ w is the Lorentz boost factor of the wall in the rest frame of the plasma, g DM is the number of spin states of the DM particle, and the DM distribution has been approximated to be Boltzmann. In the non-relativistic limit, v w → 0, filtering strongly suppresses the DM number density inside the bubble as e −mχ/T . In the relativistic limit, m χ /(γ w T ) → 0, the number density ∼ e −mχ/(2γwT ) , so there is very little filtering and n in χ approaches the equilibrium number density outside the bubble, n eq χ | T =T = g DM T 3 /π 2 . If T is lower than the thermal decoupling temperature T dec , the DM inside the bubble is already decoupled from the thermal bath and makes up the DM relic abundance, On the other hand, if T > T dec , the DM filtered by the bubble wall remains in thermal equilibrium and the relic abundance is determined by standard thermal freeze-out with m χ /T dec 24. 1 The DM abundance today can be calculated by dividing n in χ +n in χ (at T ) by the entropy density s = (2π 2 /45)g S T 3 , where g S is the effective number of relativistic degrees of freedom associated with entropy, and normalizing to the critical density, ρ c = 3H 2 0 M 2 pl [11]: (2.2) Using Eq. (2.1), this can be simplified to For example, for v w → 1, taking m χ ≈ 1 TeV and g DM = 2, requires

Bubble wall velocity
We consider a fermioinic or bosonic DM particle χ that couples to a scalar η (that could be the SM Higgs or a new particle) with coupling g χ (not to be confused with g DM , the number of spin states). The scalar undergoes a FOPT at temperature T , during which the VEV jumps from η = 0 to η = v η . Nucleation starts at T , and the bubbles expand and merge until the entire universe is populated with the v η phase. During the bubble expansion two phases coexist. Inside the bubbles η = v η , and DM gets a mass m χ g χ v η . Outside the bubbles, χ is massless because η = 0. Bubble filtering occurs as described in the previous section.
DM particles that are reflected by the bubble wall exert pressure P on it, and slows down the bubble wall velocity, which is give by the equilibrium condition ∆V = P , where ∆V is the potential energy difference between the false and true vacua. The strength of the phase transition is defined by where the radiation energy density, ρ rad (T ) = π 2 g T 4 /30, with g the number of effectively massless degrees of freedom at temperature T . For the SM, far above the electroweak scale, g 106.75. Then the pressure on the bubble wall can be obtained from the difference in the number of light degrees of freedom inside and outside the bubble [12,17,18]: where the ratio of the number of light degrees of freedom is with g b i and g f i , the number of degrees of freedom of the bosons and fermions, respectively. If particle i of the thermal plasma gains mass M i inside the bubble and 0.2M i γ w T , then most of the i particles fail to penetrate the wall and instead exert pressure on it [12]. If i is fermionic DM with g DM = 2, then d n 0.032 including particle and antiparticle contributions. Therefore, once α is known from the scalar potential, v w can be obtained by solving the equation, ∆V = P : In the limit v w → 1, with d n = 0.032, we find α 0.085γ 2 ω from Eq. (3.2). Eliminating γ w from Eq. (2.5) yields the condition, to produce the measured relic abundance for m χ ≈ 1 TeV. If we assume g χ O(1), a large v η /T O(10) and small α O(0.1) is required.

Gravitational wave production
A FOPT generates GWs from three processes [19]: i) Bubble collisions. ii) Sound waves in the plasma following bubble collisions and before the kinetic energy is dissipated by bubble expansion. iii) Magnetohydrodynamic (MHD) turbulence in the plasma after the bubble collisions. The GW signal is given by the linear combination of the three processes, The parameters that control the signal are v w , T , the phase transition strength α, the inverse of the duration of the phase transition β/H in units of the Hubble parameter at T , all of which are model and scalar potential specific. Our calculations of the GW spectra follow the semi-analytic treatment in Refs. [17,19,20]. Here, we simply point out some aspects of the three contributions without regurgitating the equations used. Increasing the values of T and β/H increases the peak GW frequency, but the latter also suppresses the power of the GW signal. The power also decreases as v w is decreased. These properties are shared by all three GW contributions.
The GW contribution from bubble collisions can be calculated directly from the scalar field η in the envelope approximation. In this approximation, an important quantity is the fraction of latent heat transformed into scalar field gradient energy, κ η .
GWs are produced by the sound waves created during percolation. For values of v w not too close to the sound speed or speed of light, parametric fits to the numerically obtained GW spectrum can be found in Ref. [19]. These fits include an efficiency parameter κ v for the fraction of latent heat transformed into bulk motion of the fluid, that depends on the expansion mode of the bubble. The peak frequency of the contribution from sound waves is inversely propositional to v w .
The contribution from MHD turbulence arises when percolation transfers a κ turb fraction of the latent heat into turbulence in the plasma. This parameter is related to κ v via κ turb κ v , where represents the fraction of bulk motion that is turbulent. The value of is still under investigation, and we conservatively take = 0.05 [19], which makes the contribution from MHD turbulence small.

Models
We now investigate two example models to demonstrate that abrupt DM freeze-out produces a detectable stochastic GW background.
We consider the Scalar Quartic Model and SU (2) X Model. Both models have a quartic term as the highest order term in their scalar potentials. However, in the former model, the effective scalar potential is composed of only one scalar field η, and may be viewed as approximating a multi-field potential. There may be thermal or non-thermal contributions to the cubic term from new particles that are not heavy enough to be integrated out [14]. In this model, the DM candidate is unspecified. On the other hand, the SU (2) X Model has the SM gauge group with an extra SU (2) X , and the scalar potential at the Planck scale only permits quartic terms built from the SM scalar doublet H and a scalar doublet S under SU (2) X . The absence of quadratic terms renders the model dimensionless at tree level. The quadratic terms and electroweak scale are dynamically generated [21,22], and the SU (2) X vector bosons are automatically stable and are the DM candidates [23,24].

Scalar Quartic Model
The effective scalar potential at finite temperature is [13][14][15][16] where we have neglected non-thermal contributions to the cubic term. We set the zerotemperature VEV to the SM value, v η = v SM = 246 GeV. Then the critical temperature is [14] T where is the temperature when the potential barrier vanishes. The two minima are There are three independent parameters ξ, D, and λ in the above effective potential. For simplicity, we fix λ = 0.1 in following analysis. Many particle physics models such as the inert singlet, inert doublet, and minimal supersymmetry models, can be parametrized by the above finite-temperature effective potential. If η is the SM Higgs, the value of ξ should be small because the Higgs portal coupling is strongly constrained by current collider data. However, if η is not the SM Higgs, ξ can be much larger as in some hidden phase transition models [14].
The nucleation temperature T n is determined by requiring the bounce action S 3 (T n )/T n 142, when the vacuum tunneling rate equals the Hubble expansion rate [13]. We adopt the following analytic approximation from Ref. [16]: where δ ≡ λ(µ 2 + DT 2 )/(ξT ) 2 , and β 1 = 8.2983, β 2 = −5.5330 and β 3 = 0.8180 are the results of a numerical fit. The expression is valid for 0 < δ < 2 which corresponds to T 0 < T < T c . We choose T = T n to compute α and β/H = d(S 3 /T )/d(ln T )| T =Tn [13]. Figure 2 shows that in only narrow parameters region (for example around (D, ξ) (18, 0.85)), is v η /T large enough ( 25) and α small enough (< 0.1), as dictated by Eq.  Table 1. of D, ξ and g χ for which Ω DM h 2 = 0.11 are displayed in the right panel of Fig. 3. The three benchmark points marked with stars in Fig. 3 are listed in Table 1.
The GW spectra for the benchmark points are shown in Fig. 4 Table 1 and Fig. 3. and BBO have sensitivity to P2. Only BBO has sensitivity to P3 because the small α = 0.037 gives a small bubble wall velocity, v w = 0.24.

SU (2) X Model
In this dimensionless model, the SM gauge group is extended by an SU (2) X with gauge coupling g X , and a scalar S, which transforms as a doublet under SU (2) X and is a singlet under the SM gauge group [8,9]. The scalar potential at tree level is SU (2) X is spontaneously broken after η acquires a VEV η = v η . We treat the three vector bosons of SU (2) X cumulatively as a single DM candidate with g DM = 9 and mass m χ = g X v η /2. In this model, as the universe cools down, the universe remains trapped in the false vacuum (i.e., η = h = 0) during thermal inflation due to the thermal effects. Around this vacuum, all particles are massless. When the energy of the false vacuum exceeds the radiation energy (i.e., α > 1), thermal inflation begins at temperature T infl with Hubble constant H * , which are given by During this phase, all particles undergo supercooling, because the scale factor grows exponentially and the temperature falls inversely with the scale factor. Supercooling ends at temperature T end with a phase transition to the true vacuum at η = v η , h = v SM . Supercooling ends when the temperature falls to the nucleation temperature T n , or earlier at the QCD phase transition temperature T QCD if T QCD > T n : where h QCD 100 MeV. To compute the GW spectra we take T = T end . From Eq. (5.6), the energy in the false vacuum is ∆V 9m 4 χ /(128π 2 ), which implies that supercooling starts at [9] T infl m χ 8.5 and To calculate T n , we use the bounce action, , for g X ≥ 1.18 (5.10) which exactly reproduces the numerical result in Fig. 1 of Ref. [9]. Nucleation occurs when S 3 (T n )/T n 4 ln(M pl /m χ ) 142. After inflation ends, the universe is reheated by the transfer of vacuum energy ∆V from the scalars to the other particles. How quickly this occurs determines the reheating temperature T rh . If the scalars decay rapidly, T rh ∼ T infl , and if they oscillate and transfer energy at a rate Γ much slower than the Hubble rate before decaying, T rh is lower, i.e., We assume that the energy transfer rate is dominated by Higgs decay, so Γ Γ h sin 2 (v SM /v η ), where the Higgs decay rate is Γ h ≈ 4 MeV.  Table 2.

Dark matter abundance
Having calculated T infl , T end and T rh , we now consider the DM relic abundance in two regimes: T rh > T dec and T rh < T dec , where T dec is the decoupling temperature in the conventional freeze-out scenario. For T rh < T dec , the DM abundance is dictated by supercooling and by sub-thermal production via scattering. Although we account for bubble filtering, its effect is negligible. On the other hand, for T rh > T dec , the supercooled population is washed out, and the sub-thermal population reattains thermal equilibrium and produces the relic abundance as in the standard freeze-out scenario. The Ω DM h 2 = 0.11 contour in the upper-left corner of Fig. 5 corresponds to this case.
The DM abundance resulting from inflationary supercooling is where f in ≡ (n in χ | T end )/(n eq χ | T end ) quantifies the filtering effect with T = T end in Eq. (2.1). However, f in = 1 for most of the parameter space of this model because the bubble wall velocity is close to the speed of light and γ w T end m χ . The dilution from supercooling is significant for T infl /T end 1, and can lead to DM being under-produced; this corresponds to the white region in the lower-right corner of Fig. 5. The DM density today can be calculated by rescaling from T rh to the temperature today, 0.235 meV, and using Eq. (2.2).
We now consider sub-thermal DM production after supercooling. The decoupling temperature of this population is T dec m χ / ln λ, where λ ≡ M pl m χ σ ann v πg /45, and σ ann v is the thermal averaged DM annihilation cross section of the DMDM → ηη process [9]. The abundance of the sub-thermal population is obtained by solving the Boltzmann equation. For T rh < T dec , both supercooling and sub-thermal production contribute to the DM relic abundance, For T rh > T dec , the plasma thermalizes again, and the usual freeze-out mechanism yields the relic abundance, 14) The DM relic abundance is shown in Fig. 5. We mark seven benchmark points along the dashed curves (which indicate Ω DM h 2 0.11), and their values are listed in Table 2. For benchmark point BP7, T rh > T dec , so the DM abundance is produced by the usual thermal freeze-out. For BP2 and BP3 sub-thermal processes dominate. Dilution by supercooling fixes the DM abundances for BP1, BP4, BP5, and BP6. The end of supercooling occurs at the nucleation temperature for BP4, BP5 and BP6, and at the QCD phase transition temperature for BP1.

Gravitational wave signals
To calculate the GW spectra, we need the phase transition strength α, inverse phase transition duration β/H , v w , and T . We evaluate α and v w by following the procedure of Section 3 and taking T = T end . The bounce action is used to find β/H ≡ d( The values of these parameters are provided in Table 2 for the seven benchmark points. The points with extremely large values of α and γ w are representative of ultra supercooling, for which the pressure P cannot counter the vacuum energy ∆V , so that bubble expansion keeps accelerating, until the bubbles collide. In the α 1 case, GWs result mainly from bubble collisions, i.e., κ η 1, κ v 0, and Ω GW h 2 Ω co h 2 .
In Fig. 6, we display the GW spectra for the seven benchmark points and the sensitivities of the LIGO O2 and O5 observing runs [25], LISA [19,26], ET [27], BBO [28], and DECIGO [29], are provided for comparison. The peak frequencies of BP1 and BP6 fall in the frequency range of LIGO and ET, because the large β/H for BP1 and large  Table 2 and Fig. 5.

BP1
, the power is suppressed by β/H to an unobservable level. Signals of BP4, BP5 and BP6 can be easily detected by BBO and DECIGO. BP4, BP5, and BP7 produce strong signals at LISA, and BP2 is marginally detectable by LISA. Finally, the frequency of BP3 is too low to be observable at these interferometers.

Summary
We studied the sudden freeze-out of DM as an alternative to the continuous thermal freezeout mechanism. A necessary ingredient for sudden freeze-out is that a FOPT generates DM mass. DM mass is generated via the coupling to a scalar particle, whose potential is responsible for a FOPT. When the scalar field acquires a non-zero VEV, DM becomes massive. The DM relic abundance may be determined by bubble filtering or by inflation and reheating. Because a FOPT triggers sudden DM freeze-out, GWs offer a signature for sudden freeze-out not available for thermal freeze-out.
To assess the viability of GWs as a signal of sudden freeze-out, we considered two example models that produce a DM abundance either by bubble-filtering (Scalar Quartic Model) or by inflation and reheating (SU (2) X Model). We showed that the observed DM relic abundance can be realized in these models with detectable GW signals in future interferometers.
In the Scalar Quartic Model, the perturbativity condition, g χ √ 4π, forces the preferred parameter space to have a large v η /T 20 and small phase transition strength, α 0.1. To produce the DM relic abundance, the expanding bubbles must filter out most of the thermal DM in the symmetric phase via a large m χ /T and non-relativistic bubble wall velocity. In these parameter regions the GW spectra have peak frequencies O(10 −2 ) Hz, and powers large enough to be probed by LISA, DECIGO, and BBO.
In the SU (2) X Model, bubble filtering has a negligible effect on the DM number density, and the DM relic abundance is governed either by supercooling during thermal inflation or sub-thermal DM production. The parameter regions that give the DM relic abundance favor α 1, which corresponds to ultra supercooling. Therefore, GWs mainly originate from bubble collisions. The GW spectra have peak frequencies that span a wide range from 10 −7 -10 3 Hz, and enough power to be easily probed by LISA, DECIGO, BBO and ET. For m χ ≥ 10 5 GeV, the GW power can be as large as Ω GW h 2 10 −7 , which may even be detected by LIGO; see BP6 in Fig. 6.