Enhanced Dark Matter Abundance in First-Order Phase Transitions
Rouzbeh Allahverdi Cash Hauptmann Peisi Huang



rouzbeh@unm.edu
Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87131,
USA


chauptmann2@huskers.unl.edu
Department of Physics and Astronomy, University of Nebraska, Lincoln, NE 68588, USA


peisi.huang@unl.edu
Department of Physics and Astronomy, University of Nebraska, Lincoln, NE 68588, USA



§ ABSTRACT 
    We propose a novel scenario to obtain the correct relic abundance for thermally under-produced dark matter. This scenario utilizes a strongly first-order phase transition at temperature T_ PT that gives rise to dark matter mass m. Freeze-out in the broken phase can yield the desired abundance in the entire region currently allowed by observational bounds and theoretical constraints for 10^2 T_ PT≲ m ≲ 10^4 T_ PT. We show that the accompanying gravitational waves are strong enough to be detected by many upcoming and proposed experiments. This, in tandem with dark matter indirect searches, provides a multi-messenger probe of such models. Positive signals in the future can help reconstruct the potential governing the phase transition and shed light on an underlying particle physics realization.   
    
    


Enhanced Dark Matter Abundance in First-Order Phase Transitions
    Peisi Huang 0000-0003-3360-2641
    Received ; accepted 
===============================================================





§ INTRODUCTION


There are various lines of evidence for the existence of dark matter (DM) in our universe <cit.>. However, the nature of DM remains a major open problem at the interface of cosmology and particle physics, and an important question concerns the relic abundance of DM that is precisely inferred from cosmological measurements <cit.>. An interesting and widely studied mechanism to address this question is thermal freeze-out in a radiation-dominated universe. The correct relic abundance can be obtained in this framework if the thermally averaged DM annihilation rate ⟨σ v ⟩ takes the nominal value of 

3 × 10^-26 cm^3 s^-1 at the time of freeze-out (coined “WIMP miracle"). Indirect detection experiments have put ⟨σ v ⟩ under increasing scrutiny. For example, Fermi-LAT’s
results from observations of dwarf spheroidal galaxies <cit.>
and newly discovered Milky Way satellites <cit.> have placed limits that rule out thermal DM with a mass below 20 GeV in a model-independent way <cit.> (barring special cases with p-wave annihilation or coannihilation).

The correct relic abundance may be obtained for larger and smaller values of ⟨σ v ⟩ if DM does not drop out of chemical equilibrium in a radiation-dominated universe. An important such scenario is DM production in nonstandard thermal histories with an epoch of early matter domination (EMD) prior to big bang nucleosynthesis (BBN) <cit.>. In the case with ⟨σ v ⟩ < 3 × 10^-26 cm^3 s^-1, which leads to thermal overproduction, entropy generation at the end of EMD can yield an acceptable relic abundance from freeze-out during EMD <cit.> or production at the end of EMD <cit.>. In the thermally under-produced case, where ⟨σ v ⟩ > 3 × 10^-26 cm^3 s^-1, the correct abundance can be obtained by a combination of DM production and its residual annihilation at the end of EMD <cit.>.     

Here we propose an alternative scenario that can yield the desired relic abundance for large values of ⟨σ v ⟩. Instead of relying on a period of EMD, this scenario employs a first-order phase transition (FOPT) within the standard thermal history. The DM particles remain light at temperatures above that of the FOPT, denoted by T_ PT, and hence stay in chemical equilibrium with the thermal bath until then. The FOPT bumps the DM mass to values m ≫ T_ PT, as a result of which DM particles undergo annihilation. For suitable values of T_ PT

one can obtain the correct DM relic abundance. Additionally, the FOPT is expected to form gravitational waves (GWs) that could be probed by upcoming detectors. This allows us to test this scenario via its GW signal in tandem with the current experimental bounds and future projections from indirect DM searches.             







Our main goal is to examine viability of the proposed scenario and its observable signatures. Therefore we remain agnostic to the specifics of underlying particle physics realizations of this scenario. That said, our results for the GW signal fit well within models accommodating supercooled FOPTs such as classically conformal B-L models <cit.>, Peccei-Quinn axion models <cit.>, strongly coupled theories <cit.>, and extra dimensional theories <cit.>. Moreover, an increase in DM mass after a FOPT may happen for any model with a fermionic DM candidate that is charged under a spontaneously broken chiral gauge symmetry <cit.>.   





















§ FREEZE-OUT AND ABUNDANCE ENHANCEMENT


The following Boltzmann equation provides the dynamics of the particle number (proper/physical) density n of DM as it evolves in a Friedmann-Lemaître-Robertson-Walker spacetime:

    nt = -3Hn - ⟨σ v ⟩[ n^2 - n_eq^2 ]   .

Hubble's parameter H ≡ ( a /  t) / a is given by the cosmic scale factor a and its rate of change with respect to cosmic time t, σ is the DM + DM annihilation cross section, and v is the Møller speed between two interacting DM particles. Angular brackets denote a thermal average: ⟨σ v ⟩≡ (∫σ v   n)/(∫ n). In deriving Eq. <ref> the DM is assumed to have a Maxwell-Boltzmann-like distribution function f ∝ e^-E/T in a Lorentz frame where E^2 = m^2 + |𝐩|^2 is the energy of DM with momentum 𝐩, and T
is its temperature. The equilibrium density n_eq is given by an exact Maxwell-Boltzmann distribution:

    n_eq = ∫ e^-E/Tg/(2π)^3 ^3 p

where g is the number of internal degrees of freedom for the DM particle. This study uses g=2 for concreteness.




Given enough time, the cosmic expansion rate dominates the Boltzmann equation and DM densities are too sparse for number-changing processes to occur. This is called freeze-out and happens at some temperature T_f. The final amount of DM can be characterized by its relic abundance Ω≡ρ_ DM / ρ_ c,0; the ratio of DM energy density to today's critical energy density. Asymptotic solutions to the Boltzmann equation find

    Ω_ thermal = s_0/ρ_c,0√(45/π)√(g_*)/g_*sm/M_Pl T_f ⟨σ v ⟩






(see Appendix <ref>). Today's entropy density is denoted with s_0, while g_* and g_*s respectively give the total relativistic degrees of freedom in the SM + DM plasma and entropy density <cit.> in a radiation-dominated universe (we take g_* = g_*s = 106.75 + 2). The Planck mass is M_ Pl = 1.22 × 10^19 GeV. This study assumes s-wave annihilation is dominant.

In the standard freeze-out analysis, Eq. <ref> is derived under the assumption that the DM number density closely tracks its equilibrium values for early times before freeze-out (called thermal DM). Given sufficient interactions between DM and the SM, this assumption is justified for DM produced in these early times, as the particles will quickly reach thermal equilibrium with the dense universe.

However, if no appreciable number of massive DM particles exist 

for T ≥ T_f then Eq. <ref> cannot be applied. This has been studied in the context of non-thermal DM production from the decay of inflaton fields <cit.> or moduli fields in supersymmetric models <cit.>. For DM produced after what would have been considered its characteristic freeze-out (had it been in equilibrium in the early universe), two cases must be considered, demarcated by some critical DM number density n_c

given by

    n_c(T) = 3H/⟨σ v ⟩  ,

which is defined for times after freeze-out (T < T_f) when the DM has fallen out of thermal equilibrium and therefore the equilibrium rate ⟨σ v ⟩ n_eq^2 in Eq. <ref> is negligible. This critical number density is found by equating the remaining expansion and annihilation rates: 3Hn_c = ⟨σ v ⟩ n_c^2.

If there is too little DM produced at T_P < T_f and n(T_P) < n_c(T_P), the cosmic expansion dominates annihilation. In this case, the relic abundance is already frozen in at T_P and is determined by the DM abundance at its production: Ω_DM = m   n(T_P) / ρ_c,0.




The more interesting case occurs when DM is produced in excess of its critical value in Eq. <ref>.











When DM is produced at some temperature T_P < T_f with density n(T_P) > n_c(T_P), the annihilation rate dominates the cosmic expansion rate and rapid annihilation will occur until the two rates are equal and the DM density begins to stabilize. The same techniques of standard freeze-out analysis can be applied to find the subsequent Boltzmann evolution of this delayed, non-thermal production of DM. The resulting relic abundance finds an enhancement from the standard thermal DM scenario:

    Ω_ DM = Ω_ thermalT_f/T_P








which is no longer a function of T_f according to Eq. <ref>.


















As long as the criterion n(T_P) > n_c(T_P) is met, the exact value of the initial DM density at its production plays no role in determining the asymptotic values, as seen in Eq. <ref>. Rapid annihilation drives the initial abundances toward a common attractor.





§ DARK MATTER PARAMETER SPACE




For a given model of DM with mass m and thermal cross section ⟨σ v ⟩, Eq. <ref> gives the necessary temperature T_P at which the DM is produced to yield some relic abundance Ω_DM. The production temperatures which yield today's relic abundance of Ω_DMh^2 = 0.12 <cit.> are shown in Fig. <ref>. Much of the DM parameter space has been ruled out by both theoretical and observational bounds. On the theoretical side, non-composite DM must have upper bounds on ⟨σ v ⟩ to maintain unitarity in 2 → 2 interactions. Following arguments first made in Ref. <cit.>, we use a unitarity bound of ⟨σ v ⟩≲ 4π / (m^2 ⟨ v ⟩) where ⟨ v ⟩≈√(6m/T_f) (see Appendix <ref>). As for observational bounds, Ref. <cit.> combines limits on DM annihilation coming from Planck satellite measurements of the cosmic microwave background <cit.>, Fermi-LAT measurements of dwarf spheroidal galaxies in the Milky Way <cit.>, and the Alpha Magnetic Spectrometer measurements of cosmic rays <cit.>. These bounds only consider “visible” final states consisting of photons, leptons, or hadrons. Lastly, we carve out the parameter space which yields an overabundance of DM; models within this region find freeze-out times T_f giving Ω_DM > 0.12. As the delayed production of DM above its critical abundance can only lead to an increase in its density, this region remains excluded for such T_P < T_f models.

With the remaining DM parameter space of Fig. <ref>, we include single-channel bounds from searches of DM annihilating to τ^- τ^+ states from HAWC <cit.>, LHAASO <cit.>, and HESS <cit.> as this channel consistently provides the most stringent bounds across different experiments. Projected sensitivities of CTA <cit.> are shown by the dashed curve. Note that a significant portion of the DM parameter space allowed by HESS is expected 

to be probed by CTA searches, yet masses near 20 GeV or 100 TeV 

will still remain out of reach.






§ MASSIVE DARK MATTER PRODUCTION IN PHASE TRANSITIONS


We now outline a simple mechanism which can produce massive DM in a cosmological FOPT occurring after its characteristic freeze-out time. Cosmological phase transitions are present in any particle physics model exhibiting symmetry breaking from an 

underlying (gauge) group. As this happens, the particles of the universe transition from a phase of interactions described by the original symmetry of the full (gauge) group to a phase of interactions given by some 

spontaneously broken subgroup. 

FOPTs are characterized by the broken phase nucleating in pockets that are stochastically distributed throughout the symmetric phase. As the FOPT proceeds, 

bubbles of broken phase expand 

until the universe is completely enveloped by the broken phase.

In this mechanism of DM “production", the DM exists prior to the phase transition, however a change in its mass occurs. The usual scenario begins with massless 

fermionic species χ (eventually the DM candidate) 

with a Yukawa 

coupling to some scalar field ϕ:

    ℒ⊃ - y ϕ^* χ_R χ_L +  h.c. ,

where χ_R and χ_L denote the right and left chiral components of χ respectively. After the phase transition

the scalar field develops an asymmetric vacuum expectation value (VEV) ⟨ϕ⟩ 







resulting in a mass term m = y⟨ϕ⟩ for the DM.



The massive DM abundance in the true vacuum is determined by the amount of massless 

χ in the false vacuum that have enough momentum to overcome the mass gap across the bubble wall. Assuming thermal equilibrium at the time of the FOPT, the number density of massless 

χ is given by n_fv(T_PT) = (3 / 4π^2) ζ(3) γ_w g T_PT^3. Therefore, if some fraction n / n_fv of 

χ particles penetrates into the true vacuum, the condition for abundance enhancement reads (n / n_fv) n_fv≥ n_c(T_PT), or

    n/n_fv⟨σ v ⟩ T_PT≥√(2/3)8 π^5/2/ζ(3) M_Pl√(g_*)/g  .


The situation can be easily extended to other massive DM production mechanisms in FOPTs, such as scalar decay outlined in Appendix <ref>.

One comment is in order at this point. 

In general, mass of the particle that mediates interactions of DM may also change during a FOPT leading to an increase in ⟨σ v ⟩ (for example, see <cit.>). However, this will not affect the final DM abundance in our scenario. This is because the number density of χ particles can acquire and track the equilibrium value from scatterings in the thermal bath even if ⟨σ v ⟩ is much (but not extremely) smaller than 3 × 10^-26 cm^3 s^-1 at T ≫ T_ PT. All that matters is that ⟨σ v ⟩ has reached its final value when the phase transition ends.   




§ PENETRATING BUBBLE WALLS
 

Eq. <ref> provides a lower bound on the abundance of DM particles in the true vacuum required for abundance enhancement. However, not all particles can penetrate into the true vacuum through the expanding bubble walls. When particles acquire a larger mass upon entering the true vacuum due to a nonzero VEV, only false vacuum particles with energies (relative to a bubble wall) above the mass difference can penetrate. More accurately, for a particle with three-momentum 𝐩 and mass m_fv in the false vacuum, it could penetrate only if |𝐩|^2 ≥ m_tv^2 - m_fv^2 where m_tv is its mass in the true vacuum. We will take m_fv = 0 in the following. To calculate the fraction of particles that successfully transition to true vacuum, we assume the false vacuum populations are relativistic and in thermal equilibrium with the matter content of the universe, and we find the flux of particles with enough momentum to penetrate the walls.

In the frame of reference of a bubble wall with relative speed v_w, a relativistic particle of spin s is described by the following one-particle phase space distribution function when in equilibrium:

    f_fv = 1/e^[γ_w |𝐩| (1 - v_w cosθ) - μ]/T - (-1)^2s  ,

where γ_w = (1-v_w^2)^-1/2 is the Lorentz factor, μ is the chemical potential, and π - θ is the angle between the particle's three-momentum 𝐩 and the wall's local velocity. The number density of fermions in the false vacuum is then

    n_fv   = ∫ f_fvg/(2π)^3^3p 
       = g/2 π^2γ_w [ 3/2ζ(3) + π^2/6μ/T + ln(2) μ^2/T^2] T^3

where g is the particle's number of internal degrees of freedom and ζ is the Riemann zeta function. The flux density of these particles which can penetrate a bubble wall is given by


    J    = ∫|𝐩| cosθ/|𝐩| f_fvΘ(|𝐩| cosθ - m) g/(2π)^3^3p 
       = g/(2π)^2[ 1/2T^2 (1-v_w)^2 γ_w^3{( m^2 (1-v_w)^2 γ_w^2 - μ^2 ) ( -iπ - 2tanh^-1[ 1+2e^m(1-v_w)γ_w-μ/T] . . . 
           . . . + ln[ 1+e^-m(1-v_w)γ_w+μ/T] ) - 2mT(1-v_w)γ_w Li_2( -e^-m(1-v_w)γ_w+μ/T) - 2T^2 Li_3( -e^-m(1-v_w)γ_w+μ/T)}] T^3


where Θ is the Heaviside step function and Li_2 is the polylogarithm function of order q. As J is the number of particles per area per time that can penetrate bubble walls, the number density of particles that pass into the true vacuum is n = J / v_w.




One can now find the fraction of particles which penetrate into the true vacuum and acquire mass:

    n/n_fv = J/v_w n_fv  .

Note that, as shown in Fig. <ref>, this fraction is a function only of wall speed v_w and x ≡ m / T whenever μ = 0 (which is assumed here). We use this relationship to find the minimum wall speed for a given FOPT to produce today's relic abundance of DM in the following way. For DM produced after its characteristic freeze-out, requiring Ω_DM h^2 = 0.12 allows Eq. <ref> to find the necessary time x_P = m / T_P at which the DM must be produced. Yet, the abundance enhancement only applies if Eq. <ref> is satisfied, requiring a minimum fraction of penetrating particles. Thus, a given x_P and minimum n / n_fv can identify a minimum wall speed.

We end this section with a discussion on the possibility of forming compact objects from DM trapped in the false vacuum. As a FOPT nears completion, bubbles of the true vacuum dominate the universe and domains or remnants of the false vacuum shrink as bubbles continue to expand. These remnants may contain those DM particles which lack the sufficient three-momenta to overcome the mass gap across bubble walls. If regions within these remnants become sufficiently dense, compact DM objects may form under their own gravitational interactions. For example, these compact objects could be black holes <cit.> or, if the DM is fermionic, an intermediate stage of Fermi-ball formation <cit.> may take place before collapsing into black holes. Whatever the mechanism, DM annihilation can severely preclude the formation of compact objects. To keep DM from completely annihilating within false vacuum remnants, particle models must have either an asymmetric distribution of DM to anti-DM populations, or some suppression of the annihilation mechanism via small Yukawa couplings or kinematic constraints. Since our analysis has no requirement for small Yukawa couplings between the DM and the FOPT-inducing scalar field, we do not consider scenarios with suppressed DM annihilation, and therefore assume no appreciable formation of compact objects before DM annihilates within (or escapes from) domains of the false vacuum.








§ GRAVITATIONAL WAVES



Particularly strong FOPTs could displace enough energy to produce a stochastic GW background detectable by future experiments which aim to measure the amplitudes and frequencies of GWs. FOPTs that release more latent heat produce GWs of greater amplitudes and are thus easier to detect. Given the stochastic nature of FOPTs, their GW signal is expected to yield a wide range of frequencies peaked around some frequency heavily influenced by the characteristic timescale of the phase transition. These two properties of a FOPT, its latent heat and timescale, are characterized by α and β:

    α   = 1/ρ_r(T){ [V_eff(ϕ_false, T) - V_eff(ϕ_true, T)] . 
       . - T [ V_eff(ϕ_false, T)T - V_eff(ϕ_true, T)T] }|_T=T_PT
    β   = - (S_3 / T)t|_t=t_PT
              .

Here, ρ_r(T) = π^2 g_*(T) T^4 / 30 is the energy density of the radiation bath at the time of the FOPT with g_* the relativistic degrees of freedom in the plasma. Subscripts “true” and “false” denote quantities evaluated in the true and false vacuua. More useful in GW calculations is β / H(T_PT) where H is the Hubble parameter. T_PT is a characteristic temperature at which the FOPT occurs. In what follows, we also take this to be the temperature at which massive DM is produced (T_PT = T_P). It should be noted that T_PT is a post-supercooling temperature and all reheating effects take place before this time.

In a FOPT, the three dominant contributions to a GW background come from (1) sound waves in the universal plasma, (2) collisions of true vacuum bubbles, and (3) magnetohydrodynamic turbulence in the plasma. Following Ref. <cit.>, we approximate the total GW background to be a linear combination of these stochastic sources:

    Ω_GW = Ω_sw + Ω_col + Ω_turb

with Ω = ρ / ρ_c the corresponding fraction of the critical energy density, however GW signal spectra are usually presented in terms of Ω(f) ≡dΩ / dln(f) <cit.>, which is the notation used here. To determine each of these contributions, one needs knowledge of the fractional amount of available latent heat that transfers into (1) the bulk motion of the fluid (κ_sw), (2) the bubble wall kinetic energy (κ_col), and (3) turbulence (κ_turb). The sound wave contribution is numerically fitted as

    h^2 Ω_sw(f)    = 3 × 2.65 × 10^-6[ H(T_PT)/β] [ κ_swα/1+α]^2 
       ×[ 100/g_*(T_PT)]^1/3 v_w [ f/f_sw]^3 [ 7/4 + 3(f/f_sw)^2]^7/2 , 
    
            f_sw   = 1.15/v_w[ β/H(T_PT)] h_*    , 
    
            h_*    = 1.65 × 10^-5 Hz[ T_PT/100  GeV] [ g_*(T_PT)/100]^1/6

where κ_sw is the fraction of latent heat transferred into the bulk motion of the plasma, and g_*(T_PT) is the number of relativistic degrees of freedom for all fields when the FOPT occurs. (The prefactor of 3 is found in the erratum of Ref. <cit.>.) The effective number of degrees of freedom coming from relativistic SM particles changes with temperature <cit.>. Redshifting the peak frequencies to their expected values today is accomplished with the h_* factor which assumes adiabatic expansion of a radiation-dominated universe. A numerical fit for the contributions from magnetohydrodynamic turbulence is found to be

    h^2 Ω_turb(f)    = 3.35 × 10^-4[ H(T_PT)/β] [ κ_turbα/1+α]^3/2
       ×[ 100/g_*(T_PT)]^1/3 v_w 
       ×[ (f/f_turb)^3/[1+(f/f_turb)]^11/3 [1+8π f/h_*]] , 
    
            f_turb   = 1.64/v_w[ β/H(T_PT)] h_*

where κ_turb gives the fraction of latent heat transferred into plasma turbulence. Finally, for a strongly supercooled phase transition, we use lattice simulation results from Ref. <cit.> to provide a numerical fit for the contributions from collisions of bubble walls: 

    h^2 Ω_col(f)    = 1.67 × 10^-5[ H(T_PT)/β]^2 [ κ_colα/1+α]^2 
       ×[ 100/g_*(T_PT)]^1/3[ 1 + (2π f/0.13 β)^d-a/1 + (2π f_col/0.13 β)^d-a] 
       ×3.63 × 10^-2× (a+b)^c/[b(f/f_col)^-a/c + a(f/f_col)^b/c]^c, 
    
            f_col   = 0.81/2π[ β/H(T_PT)] h_*   .

Here, κ_col giving the fraction of latent heat converted into kinetic energy of the bubble walls. The numerical parameters are a=2.54, b=2.24, c=2.30, and d=0.93.

To meet relic abundance requirements, Fig. <ref> shows that the DM must be produced at times x_P ≡ m/T_P ∼ [10^2, 10^4]. If DM mass is proportional to the scalar field's VEV, high x_P indicates the FOPTs should occur at temperatures well below the VEV, given that 

y ≤𝒪(1) in the perturbative regime. 

This suggests the FOPT is supercooled, meaning it completes at some temperature well below its critical temperature (the temperature, typically on the same energy scale of the VEV, at which the effective potential acquires degenerate minima). With supercooling present, at the end of the FOPT a large difference in the effective potential obtains between the true and false vacuum phases, releasing much more energy than a non-supercooled transition. For this analysis, we thus assume α≫ 1 meaning the contributions to Ω_GW(f) are practically independent of α[For α≈ 1, the GW signal strength is only reduced by a factor of about four. Even with this reduction, Fig. <ref> shows GW signals are still well within projected detector sensitivities. On the other hand, α should not be allowed arbitrarily large values to avoid an extended period of vacuum energy domination, requiring adjustments to the GW signal calculations.]. In taking α≫ 1 limits, we also calculate κ_sw = κ_sw(α, v_w) using numerical fits found in Ref. <cit.>. For values of v_w required in these fits and in Ω_GW(f), we use the minimum wall speed required for abundance enhancement outlined in Sec. <ref>. In the remaining GW spectra parameters κ_col and κ_turb, simulations <cit.> have shown power in transverse modes of plasma fluid velocity to be at most 5-10% of that in longitudinal modes. We thus conservatively take κ_turb = 0.05 κ_sw. Finally, to find κ_col we enforce κ_sw + κ_turb + κ_col = 1 which assumes all latent heat is dispersed among sound waves, turbulence, and bubble collisions.





§ RESULTS


To calculate GW spectra, we take points in the m-⟨σ v ⟩ plane, shown in Fig. <ref>, which yield the correct DM relic abundance. Furthermore, we do not consider the points excluded by visibles and unitarity bounds. Within this study, the only free parameter left in calculating GW spectra is β / H(T_PT), which we take to be greater than one to ensure bubble nucleation rates exceed universal expansion rates and the FOPT can complete. Peak values of GW signals from FOPTs are plotted in Fig. <ref> to compare against expected sensitivities of proposed GW detection experiments. It is clear that supercooled FOPTs with α≫ 1 help provide GW signals strong enough to reach expected sensitivities for a wide range of β / H(T_PT). Although the correlation is not exact, we can also see that in models endowing larger masses to DM fields, the resulting GW signals peak at higher frequencies; suggesting these models require more rapid FOPTs. As mentioned above, each signal peak is determined by a point in Fig. <ref> which is not excluded by bounds on visible annihilation channels, unitarity, and overabundance. The peaks are further distinguished by their shape: disks (∘) are points permitted by HESS bounds while crosses (×) are not. In addition to evading HESS bounds, disks with black outlines are probeable by projected sensitivities of CTA.



Finally, in Fig. <ref> we offer a fuller look at possible DM candidates in a parameter space mixing particle characteristics and GW data from the accompanying FOPT. In this plot, the DM particle mass is given by the vertical axis while thermally averaged cross sections color each point. The horizontal axis provides the frequencies of peak GW signals, and we include three frequency bands demarcating the projected sensitivity ranges of THEIA, LISA, ALIA, and ET. Note that although a given point may fall within one of these frequency bands, its GW signal strength will determine whether the point might be probed by some detector. As in Fig. <ref>, disks denote points compatible with HESS bounds while crosses do not; disks with black outlines 

can be further probed by CTA. The four distinct (although overlapping) sets of points correspond to different β / H(T_PT) values in the FOPTs, labeled near the bottom of the plot.






§ CONCLUSION


In this paper, we presented a new scenario for obtaining the correct relic abundance in the case of thermally under-produced DM, i.e., ⟨σ v ⟩ > 3 × 10^-26 cm^3 s^-1. This scenaio utilizes a FOPT in a radiation-dominated universe as a result of which DM acquires its mass m. DM particles maintain a constant comoving density above FOPT temperature T_ PT and their freeze-out in the broken phase yields the correct abundance if T_ PT is sufficiently lower than the freeze-out temperature in the standard scenario T_ f∼ m/20. 

In particular, we found that 10^2 T_ PT≲ m ≲ 10^4 T_ PT within the mass range allowed by current experimental limits from indirect searches and theoretical bounds from unitarity. The hierarchy between m and T_ P suggests that the FOPT must be supercooled thereby leading to potentially significant GW signals. We found that although the GWs sourced by such a FOPT would likely not be detectable by current experiments, they are within projected sensitivities of many next-generation GW detectors. 


As shown in Fig. <ref>, for reasonable values of PT parameters, the GW signal is within the reach of experiments like ET, LISA, ALIA, and THEIA. This, in tandem with indirect detection experiments, provides a multi-messenger test of our scenario. We see in Fig. <ref> that almost the entire region of the ⟨σ v ⟩-m plane that is allowed by the FERMI-LAT constraints can be probed by GW detectors. Future indirect searches like CTA can access the bulk of this region that include masses between the weak scale and 10^4 GeV.         

The crucial ingredient of our scenario is a strongly FOPT at temperatures below a few TeV which gives rise to DM mass. We largely focus on general conditions for its success and prospects of its multi-messenger probes. A natural direction for future investigations involves model-specific studies. Detection of GWs combined with a signal from indirect DM searches can help us reconstruct the potential governing the dynamics of the FOPT. This can then be used to build explicit particle physics models employing a FOPT to give mass to the DM. 






§ ACKNOWLEDGEMENTS


The work of R.A. is supported in part by NSF Grant No. PHY-2210367. The work of C.H. and P.H. is supported by the National Science Foundation under grant number PHY-2112680. 





§ APPENDIX




 §.§ Freeze-Out of Cold Relics




The Boltzmann equation reveals three processes responsible for the changing DM number density. Cosmic expansion and self-annihilation deplete the number density at rates of 3Hn and ⟨σ v ⟩ n^2, respectively. Equilibrium interactions at high temperatures help maintain the number density at a rate ⟨σ v ⟩ n_eq^2. To factor out effects of cosmic expansion, one can define the yield of DM as Y ≡ n / s where s is the entropy density of the universe. During cosmic expansion n and s scale equivalently with a^-3, allowing Eq. <ref> to be rewritten as

    Yt = - ⟨σ v ⟩ s [ Y^2 - Y_eq^2 ]   ,
            
        Y ≡ n/s

with Y_eq≡ n_eq / s. One can further adimensionalize this equation with the parameterization of x ≡ m / T. Assuming the DM to be in thermal equilibrium with a radiation-dominated universe, the Friedmann equations give x /   t = Hx and Eq. <ref> becomes

    Yx = - ⟨σ v ⟩ s/Hx[ Y^2 - Y_eq^2 ]   ,
            
        x ≡m/T  .

The Friedmann equations also yield the following expressions for the Hubble parameter and entropy density:

    H^2 = 8 π/3 M_Pl^2 g_*(T) π^2/30 T^4   ,
            
        s = 2π^2/45 g_*s(T) T^3

with Planck mass M_Pl = 1.22 × 10^19 GeV. In these expressions g_* and g_*s give the total degrees of freedom in the universal plasma and entropy density, respectively <cit.>. This study approximates these degrees of freedom to be equal: g_* = g_*s. The sum of SM and DM relativistic degrees of freedom gives the total number of degrees of freedom in the plasma, that is g_* = g_SM + g.



Typical solutions of the Boltzmann equation are displayed in Fig. <ref> where Y = Y_eq is taken as a boundary condition for early times (small x). The universe is dense and hot throughout these times, keeping the DM yield near its thermal equilibrium value Y_eq, while self-annihilation decreases the number of DM particles. Eventually the expansion of the universe “freezes out" DM annihilation processes in that the number of DM particles within a comoving volume becomes too low for annihilation rates to significantly affect Y(x). The approximate time at which this occurs is denoted by x_f. With a boundary condition specified for some unique Boltzmann equation, a characteristic freeze-out for the system is determined. For x > x_f, the DM comoving abundance is said to be “frozen in" as the number of DM particles per comoving volume remains relatively constant after freeze-out. As long as freeze-out occurs well into the past, the frozen-in DM abundance determines its relic abundance as measured today: Ω_DM h^2 = 0.12 <cit.>. Here, h is determined by measurement of today's value of the Hubble parameter H_0 = 100   h   km s^-1 Mpc^-1. With knowledge of the frozen-in DM abundance, we expect to find

    Ω_DM = m Y(x_today) s_0/ρ_c, 0

for nonrelativistic DM. Although Eq. <ref> has no known closed-form solutions, one can find approximations for today's yield Y(x_today) for two types of DM relics. Hot relics are models of DM wherein freeze-out occurs while the species is still relativistic, that is x_f ≪ 3. Models with nonrelativistic DM at freeze-out (x_f ≫ 3) are known as cold relics. The next section follows a standard analysis that can be found in Ref. <cit.>, for example.


The first simplifying assumption is made on the thermally averaged annihilation cross section ⟨σ v ⟩. It is expected that σ v ∝ v^k/2 where k=0 corresponds to s-wave annihilation and k=1 corresponds to p-wave annihilation. As ⟨ v ⟩∝ T^1/2, we parameterize the cross section as

    ⟨σ v ⟩   ≡σ_0 (T/m)^k 
       = σ_0 x^-k  
            k = 0    ⟺ s-wave    annihilation
    
            k = 1    ⟺ p-wave    annihilation  .

Our study assumes s-wave annihilation is the dominant contribution and thus k=0. This parameterization allows Eq. <ref> to be rewritten as

    Yx   = - λ x^-k-2[ Y^2 - Y_eq^2 ]   , 
    λ≡x σ_0 s/H   = [σ_0 s/H]_x=1 
        = √(π/45)g_*s/√(g_*) M_Pl m σ_0

where a new parameter λ has been defined which is approximately independent of x.






To track the yield's deviation from its equilibrium value, we express the right-hand side of Eq. <ref> in terms of Δ≡ Y - Y_eq to find

    Yx = -λ x^-k-2Δ (2 Y_eq + Δ)   .

As Y_eq decreases exponentially with x, we expect Δ≈ Y with cold relics for times after freeze-out, meaning Eq. <ref> is approximately

    Δx≈ -λ x^-k-2Δ^2
            
        (x ≫ x_f)   .

Despite the validity of this approximation holding for x ≫ x_f, the standard method is to integrate the above equation from x_f out to some x ≫ x_f to find

    ∫_Δ(x_f)^Δ(x)Δ/Δ^2   ≈ -λ∫_x_f^x x^'/x^'^k+2
    [ 1/Δ(x) - 1/Δ(x_f)]    ≈-λ/k+1[ 1/x^k+1 - 1/x_f^k+1] 
    Δ(x)    ≈k+1/λ[ 1/x_f^k+1 - 1/x^k+1]^-1
    Δ(x)    ≈k+1/λ x_f^k+1  .

The third line of Eqs. <ref> ignores 1/Δ(x_f) since the yield at freeze-out is typically larger than its value at some later time. The last line utilizes x ≫ x_f and rids the dependence on x. With Δ(x) ≈ Y(x), the standard approximation for today's DM yield is finally given by

    Y_st(x_today) ≈k+1/λ x_f^k+1  .

If instead the DM is produced after its characteristic freeze-out time, we can simply integrate Eq. <ref> from a lower bound given by the time of production x_P. With x_ today≫ x_P > x_f, a similar result is found:

    Y(x_today) ≈k+1/λ x_P^k+1  .




 §.§ Thermally Averaged Møller Speed


The Møller speed between two particles with four-momenta p_1^μ = (E_1, 𝐩_1) and p_2^μ = (E_2, 𝐩_2) is

    v ≡√((p_1^μ p_2 μ)^2 - (m_1 m_2)^2)/E_1 E_2

where m_i^2 ≡ p_i^μ p_i μ. Consider a colinear collision between particles of the same mass: m_1 = m_2 ≡ m. Furthermore, in taking a thermal average, we should expect each particle to have the same energy and momentum (although antiparallel), which we will call E and |𝐩|. We then find

    ⟨ v ⟩^2    = (E^2 + |𝐩|^2)^2 - m^4/E^4
       = 1 + ( |𝐩|/E)^4 + 2( |𝐩|/E)^2 - ( m/E)^4 
       = 1 + ( 3T/γ m)^2 + 2( 3T/γ m) - 1/γ^4
       = 6/x√(1-v_ rel) + 9/x^2 (1-v_ rel) + v_ rel^2 (2 - v_ rel^2)   .

The relativistic equipartition theorem of 3T = |𝐩|^2 / E (along with E = γ m) was used in the third line, and γ≡ 1 / √(1-v_ rel) is the Lorentz factor between the two rest frames of the particles. We thus see that for cold relics with x_f ≫ 3 and nonrelativistic relative speeds v_ rel≪ 1, the Møller speed at freeze-out is ⟨ v ⟩≈√(6/x_f).



 §.§ Scalar Decay



Consider scalar DM particles Δ produced inside the bubble walls of the developing broken phase via decay of a symmetry-breaking scalar field ϕ. For scalar DM, the relevant terms in the Lagrangian density are 

    ℒ⊃ - 1/2 m_Δ^2 Δ^2 - 1/2λΔ^2 |ϕ|^2

The interaction term develops a three-point vertex when ϕ acquires a nonzero VEV, inducing a phase transition which we take to be of first order. Letting ϕ→ϕ + ⟨ϕ⟩, the interaction term becomes

    1/2λΔ^2 |ϕ + ⟨ϕ⟩|^2 = 
        1/2λΔ^2 ( |ϕ|^2 + ⟨ϕ⟩^2 + 2 ⟨ϕ⟩ |ϕ| )   .

Here, the first term describes a quartic interaction between the DM and scalar as in the original Lagrangian, the second term provides the DM with additional mass, and the third term is a trilinear coupling between DM and scalar. This last term is responsible for the decay of ϕ into DM within bubble walls where ⟨ϕ⟩≠ 0.

In this mechanism, the abundance of DM in the true vacuum phase is determined by the ϕ→Δ decay process, as well as how many scalars are able to penetrate the bubble walls. Suppose the scalar field has some yield Y^ϕ_fv = n^ϕ_fv / s and decays into N DM particles with probability |𝒯_ϕ→ N Δ|^2 at any given time. If some fraction Y^ϕ / Y^ϕ_fv of scalar particles penetrate the bubble walls into the true vacuum, the yield of DM in the true vacuum is given by

    Y = Y^ϕ/Y^ϕ_fvn^ϕ(T_PT)/s(T_PT) N |𝒯_ϕ→ N Δ|^2   .

Prior to the FOPT, we assume the scalar field ϕ to be in thermal equilibrium with the field content of the universe so that n^ϕ(T) = ζ(3) γ_w g_ϕ T^3 / π^2 which gives

    Y = Y^ϕ/Y^ϕ_fv45 ζ(3) g_ϕ/2π^4 g_*,s N |𝒯_ϕ→ N Δ|^2   .

To meet the condition of DM abundance to be enhanced from its standard freeze-out value, we require Y ≥ Y_c(T_PT). Combining this condition with Eqs. <ref> and <ref> yields

    Y^ϕ/Y^ϕ_fv⟨σ v ⟩ T_PT≥2√(6)π^5/2/ζ(3) M_Pl√(g_*)/g_ϕ1/N |𝒯_ϕ→ N Δ|^2  .





