We propose a novel scenario to obtain the correct relic abundance for thermally underproduced dark matter. This scenario utilizes a strongly firstorder phase transition at temperature Trm PT that gives rise to dark matter mass m. Freezeout in the broken phase can yield the desired abundance in the entire region currently allowed by observational bounds and theoretical constraints for 102 Trm PT lesssim m lesssim 104 Trm 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 multimessenger 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.
['High Energy Physics - Phenomenology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Theory']
Enhanced Dark Matter Abundance in FirstOrder Phase Transitions
{'Enhanced Dark Matter Abundance in First-Order Phase Transitions': 'Rouzbeh Allahverdi , 1, ∗ Cash Hauptmann , 2, † and Peisi Huang 2, ‡  \n1 Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87131, USA 2 Department of Physics and Astronomy, University of Nebraska, Lincoln, NE 68588, USA  \nWe 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.', 'I. INTRODUCTION': "There are various lines of evidence for the existence of dark matter (DM) in our universe [1]. 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 [2-5]. 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 [6] and newly discovered Milky Way satellites [7] have placed limits that rule out thermal DM with a mass below 20 GeV in a model-independent way [8] (barring special cases with p -wave annihilation or coannihilation).  \nThe 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) [9]. 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 [10-14] or production at the end of EMD [15, 16]. In the thermally under-produced case, where ⟨ σv ⟩ > 3 × 10 -26 cm 3 s -1 , the correct abundance can be obtained by a combina-  \ntion and its residual annihilation at the end of EMD [17, 18].  \nHere 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.  \nOur 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 [19-23], Peccei-Quinn axion models [24, 25], strongly coupled theories [26], and extra dimensional theories [27, 28]. 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 [29].", 'II. 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ˆıtreRobertson-Walker spacetime:  \nd n d t = -3 Hn -⟨ σv ⟩ [ n 2 -n 2 eq ] . (1)  \nHubble's parameter H ≡ (d a/ d 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 d n ) / ( ∫ d n ). In deriving Eq. 1 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 + | p | 2 is the energy of DM with momentum p , and T is its temperature. The equilibrium density n eq is given by an exact Maxwell-Boltzmann distribution:  \nn eq = ∫ e -E/T g (2 π ) 3 d 3 p (2)  \nwhere g is the number of internal degrees of freedom for the DM particle. This study uses g = 2 for concreteness.  \nGiven 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  \nΩ thermal = s 0 ρ c, 0 √ 45 π √ g ∗ g ∗ s m M Pl T f ⟨ σv ⟩ (3)  \n(see Appendix A). 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 [30] 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.  \nIn the standard freeze-out analysis, Eq. 3 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. 3 cannot be applied. This has been studied in the context of non-thermal DM production from the decay of inflaton fields [31] or moduli fields in supersymmetric models [3234]. 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  \nn c ( T ) = 3 H ⟨ σv ⟩ , (4)  \nwhich is defined for times after freeze-out ( T < T f ) when the DM has fallen out of thermal equilibrium  \nand therefore the equilibrium rate ⟨ σv ⟩ n 2 eq in Eq. 1 is negligible. This critical number density is found by equating the remaining expansion and annihilation rates: 3 Hn c = ⟨ σv ⟩ n 2 c .  \nIf 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 = mn ( T P ) /ρ c, 0 . The more interesting case occurs when DM is produced in excess of its critical value in Eq. 4.  \nWhen 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:  \nΩ DM = Ω thermal T f T P (5)  \nwhich is no longer a function of T f according to Eq. 3. 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. 5. Rapid annihilation drives the initial abundances toward a common attractor.", 'III. DARK MATTER PARAMETER SPACE': "For a given model of DM with mass m and thermal cross section ⟨ σv ⟩ , Eq. 5 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 Ω DM h 2 = 0 . 12 [5] are shown in Fig. 1. 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. [35], we use a unitarity bound of ⟨ σv ⟩ ≲ 4 π/ ( m 2 ⟨ v ⟩ ) where ⟨ v ⟩ ≈ √ 6 m/T f (see Appendix B). As for observational bounds, Ref. [8] combines limits on DM annihilation coming from Planck satellite measurements of the cosmic microwave background [40], Fermi-LAT measurements of dwarf spheroidal galaxies in the Milky Way [6, 7], and the Alpha Magnetic Spectrometer measurements of cosmic rays [41, 42]. 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  \nFIG. 1. DM parameter space with cross section upper bounds on s -wave annihilation to visible final states [8] and unitarity [35]. Additional bounds are provided by searches of τ -τ + DMannihilation channels from HAWC [36], LHAASO [37], and HESS [38]. The dashed curve provides the expected sensitivity of CTA [39] in the same channel. All τ -τ + searches assume cusped DM density profiles in the Milky Way (Einasto profiles for HAWC, HESS, and CTA; Navarro-Frenk-White profile for LHAASO).  \n<!-- image -->  \nlead to an increase in its density, this region remains excluded for such T P < T f models.  \nWith the remaining DM parameter space of Fig. 1, we include single-channel bounds from searches of DM annihilating to τ -τ + states from HAWC [36], LHAASO [37], and HESS [38] as this channel consistently provides the most stringent bounds across different experiments. Projected sensitivities of CTA [39] 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.", 'IV. 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 sub-  \ngroup. 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.  \nIn 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 ϕ :  \nL ⊃ yϕ ∗ χ R χ L +h . c . , (6)  \nwhere χ R and χ L denote the right and left chiral Weyl components of χ respectively. Our calculations use g = 2 and therefore apply to Majorana DM ( χ R = χ L ). 1 After the phase transition the scalar field develops an asymmetric vacuum expectation value (VEV) ⟨ ϕ ⟩ resulting in a mass term m = y ⟨ ϕ ⟩ for the DM.  \nThe 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  \ngap 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 gT 3 PT . 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  \nn n fv ⟨ σv ⟩ T PT ≥ √ 2 3 8 π 5 / 2 ζ (3) M Pl √ g ∗ g . (7)  \nThe situation can be easily extended to other massive DM production mechanisms in FOPTs, such as scalar decay outlined in Appendix C.  \nOne 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 [43, 44]). 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.", 'V. PENETRATING BUBBLE WALLS': "Eq. 7 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  \ntrue 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 p and mass m fv in the false vacuum, it could penetrate only if | p | 2 ≥ m 2 tv -m 2 fv 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.  \nIn 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:  \nf fv = 1 e [ γ w | p | (1 -v w cos θ ) -µ ] /T -( -1) 2 s , (8)  \nwhere γ w = (1 -v 2 w ) -1 / 2 is the Lorentz factor, µ is the chemical potential, and π -θ is the angle between the particle's three-momentum p and the wall's local velocity. The number density of fermions in the false vacuum is then  \nn fv = ∫ f fv g (2 π ) 3 d 3 p = g 2 π 2 γ w [ 3 2 ζ (3) + π 2 6 µ T +ln(2) µ 2 T 2 ] T 3 (9)  \nwhere 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  \nJ = ∫ | p | cos θ | p | f fv Θ( | p | cos θ -m ) g (2 π ) 3 d 3 p = g (2 π ) 2 [ 1 2 T 2 (1 -v w ) 2 γ 3 w { ( m 2 (1 -v w ) 2 γ 2 w -µ 2 ) ( -iπ -2 tanh -1 [ 1 + 2 e m (1 -vw ) γw -µ T ] +ln [ 1 + e -m (1 -vw ) γw + µ T ]) -2 mT (1 -v w ) γ w Li 2 ( -e -m (1 -vw ) γw + µ T ) -2 T 2 Li 3 ( -e -m (1 -vw ) γw + µ T )}] T 3 (10)  \nwhere Θ is the Heaviside step function and Li q 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:  \nn n fv = J v w n fv . (11)  \nNote that, as shown in Fig. 2, this fraction is a function only of wall speed v w and x ≡ m/T whenever µ = 0, which holds for Majorana fermions whose particle number freely changes 2 . We use this relationship to find the  \nminimum 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. 5 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. 7 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. Choosing the minimum penetrating fraction n/n fv as outlined above is not necessary to the DM production mechanism-this fraction can take any value satisfying Eq. 7 to overcome the critical DM number density. This is in contrast to mechanisms using 'filtered DM' [45, 46], which rely in part on inefficient DM penetration.  \nWe 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 [47-49] or, if the DM is fermionic, an intermediate stage of Fermi-ball formation [50-53] 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.", 'VI. 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  \nFIG. 2. The fraction of particles able to penetrate into true vacuum and acquire mass for different bubble wall speeds (line labels). A chemical potential of µ = 0 is assumed.  \n<!-- image -->  \nand 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 β :  \nα = 1 ρ r ( T ) { [ V eff ( ϕ false , T ) -V eff ( ϕ true , T )] -T [ d V eff ( ϕ false , T ) d T -d V eff ( ϕ true , T ) d T ]}∣ ∣ ∣ ∣ T = T PT β = -d( S 3 /T ) d t ∣ ∣ ∣ ∣ t = t PT . (12)  \nHere, ρ 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.  \nIn 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. [54], we approximate the total GW background to be a linear combination of these stochastic sources:  \nΩ GW = Ω sw +Ω col +Ω turb (13)  \nwith Ω = ρ/ρ c the corresponding fraction of the critical energy density, however GW signal spectra are usually presented in terms of Ω( f ) ≡ dΩ / dln( f ) [55, 56], 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 ).  \nThe sound wave contribution is numerically fitted as  \nh 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 × H ( T PT ) τ sw , 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 (14)  \nwhere κ 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. [57].) The effective number of degrees of freedom coming from relativistic SM particles changes with temperature [30]. Redshifting the peak frequencies to their expected values today is accomplished with the h ∗ factor which assumes adiabatic expansion of a radiation-dominated universe. For conservative estimates, a suppression factor of Hτ sw is used. This factor accounts for the finite lifetime of sound waves as sources of GWs, and is given by [58, 59]  \nHτ sw = min[1 , Hτ sh ] ,  \nHτ sh ≈ (8 π ) 1 / 3 v w ( β H ) -1 ( 3 4 κ sw α 1 + α ) -1 / 2 . (15)  \nA numerical fit for the contributions from magnetohydrodynamic turbulence is found to be [58]  \nh 2 Ω turb ( f ) = 3 . 35 × 10 -4 [ H ( T PT ) β ][ κ sw α 1 + α ] 3 / 2 × [ 100 g ∗ ( T PT ) ] 1 / 3 v w [1 -H ( T PT ) τ sw ] × [ ( f/f turb ) 3 [1 + ( f/f turb )] 11 / 3 [1 + 8 πf/h ∗ ] ] , f turb = 1 . 64 v w [ β H ( T PT ) ] h ∗ . (16)  \nFinally, for a strongly supercooled phase transition, we use lattice simulation results from Ref. [60] to provide  \na numerical fit for the contributions from collisions of bubble walls:  \nh 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 ∗ . (17)  \nHere, κ 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.  \nTo meet relic abundance requirements, Fig. 1 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 ≤ O (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 α 3 . In taking α ≫ 1 limits, we also calculate κ sw = κ sw ( α, v w ) using numerical fits found in Ref. [61]. 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. V. In the remaining GW spectra parameters κ col and κ turb , simulations [62] 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.  \nFIG. 3. Peak values of GW signal for different β/H ( T PT ). Each peak is colored to its corresponding DM mass. Disks are points compatible with HESS bounds while crosses are not. Additionally, disks outlined in black are within CTA's projected sensitivities. Colored lines show the projected sensitivities of the following proposed experiments: A+ [63], ALIA [64], aLIGO [65, 66], BBO [67, 68], CE [69], DECIGO [70, 71], EPTA [72], ET [73], Gaia [74], LISA [75-77], NANOGrav [78, 79], SKA [80], Taiji [81], THEIA [82], TianQin [83].  \n<!-- image -->", 'VII. RESULTS': 'To calculate GW spectra, we take points in the m -⟨ σv ⟩ plane, shown in Fig. 1, 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. 3 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. 1 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  \nby projected sensitivities of CTA.  \nFinally, in Fig. 4 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 four 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, hence our omission of points with β/H ( T PT ) = 10 6 from Fig. 4. As in Fig. 3, disks denote points compatible with HESS bounds while crosses do not; disks with black outlines can be further probed by CTA. The three distinct (although overlapping) sets of points correspond to different β/H ( T PT ) values in the FOPTs, labeled near the bottom of the plot.', 'VIII. 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  \nFIG. 4. Mixed space of DM particle properties and GW parameters. Point shapes give the same designations as Fig. 3.  \n<!-- image -->  \nscenaio 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.  \nIn 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.  \nAs shown in Fig. 3, 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 multimessenger test of our scenario. We see in Fig. 4 that almost the entire region of the ⟨ σv ⟩ -m plane that is allowed by the FERMI-LAT constraints can be probed by GWdetectors. Future indirect searches like CTA can access the bulk of this region that include masses between the weak scale and 10 4 GeV.  \nThe crucial ingredient of our scenario is a strongly FOPT at temperatures below a few TeV which gives  \nrise 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.', 'ACKNOWLEDGMENTS': '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.', 'A. 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 3 Hn and ⟨ σv ⟩ n 2 , respectively. Equilibrium interactions at high temperatures help maintain the number density at a rate ⟨ σv ⟩ n 2 eq . To factor out effects  \nFIG. 5. Typical solutions to Eq. 19 with the boundary condition Y (1) = Y eq (1). Departure from Y eq marks the onset of freeze-out, which occurs at later times for higher ⟨ σv ⟩ .  \n<!-- image -->  \nof 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. 1 to be rewritten as  \nd Y d t = -⟨ σv ⟩ s [ Y 2 -Y 2 eq ] , Y ≡ n/s (18)  \nwith 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 d x/ d t = Hx and Eq. 18 becomes  \nd Y d x = -⟨ σv ⟩ s Hx [ Y 2 -Y 2 eq ] , x ≡ m T . (19)  \nThe Friedmann equations also yield the following expressions for the Hubble parameter and entropy density:  \nH 2 = 8 π 3 M 2 Pl g ∗ ( T ) π 2 30 T 4 , s = 2 π 2 45 g ∗ s ( T ) T 3 (20)  \nwith 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 [30]. This study approximates these degrees of freedom to be equal: g ∗ = g ∗ s . The sum of Standard Model and DM relativistic degrees of freedom gives the total number of degrees of freedom in the plasma, that is g ∗ = g SM + g .  \nTypical solutions of the Boltzmann equation are displayed in Fig. 5 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.  \nEventually 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 [5]. Here, h is determined by measurement of today's value of the Hubble parameter H 0 = 100 h kms -1 Mpc -1 . With knowledge of the frozen-in DM abundance, we expect to find  \nΩ DM = mY ( x today ) s 0 ρ c, 0 (21)  \nfor nonrelativistic DM. Although Eq. 19 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. [84], for example.  \nThe 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  \n⟨ σv ⟩ ≡ σ 0 ( T/m ) k = σ 0 x -k { k = 0 ⇐⇒ s -wave annihilation k = 1 ⇐⇒ p -wave annihilation . (22)  \nOur study assumes s -wave annihilation is the dominant contribution and thus k = 0. This parameterization allows Eq. 19 to be rewritten as  \nd Y d x = -λx -k -2 [ Y 2 -Y 2 eq ] , λ ≡ xσ 0 s H = [ σ 0 s H ] x =1 = √ π 45 g ∗ s √ g ∗ M Pl mσ 0 (23)  \nwhere 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. 23 in terms of ∆ ≡ Y -Y eq to find  \nd Y d x = -λx -k -2 ∆(2 Y eq +∆) . (24)  \nAs Y eq decreases exponentially with x , we expect ∆ ≈ Y with cold relics for times after freeze-out, meaning Eq. 24 is approximately  \nd∆ d x ≈ -λx -k -2 ∆ 2 ( x ≫ x f ) . (25)  \nDespite 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  \n∫ ∆( x ) ∆( x f ) d∆ ∆ 2 ≈ -λ ∫ x x f d x ' x ' k +2 [ 1 ∆( x ) -1 ∆( x f ) ] ≈ -λ k +1 [ 1 x k +1 -1 x k +1 f ] ∆( x ) ≈ k +1 λ [ 1 x k +1 f -1 x k +1 ] -1 ∆( x ) ≈ k +1 λ x k +1 f . (26)  \nThe third line of Eqs. 26 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  \nY st ( x today ) ≈ k +1 λ x k +1 f . (27)  \nIf instead the DM is produced after its characteristic freeze-out time, we can simply integrate Eq. 25 from a lower bound given by the time of production x P . With x today ≫ x P > x f , a similar result is found:  \nY ( x today ) ≈ k +1 λ x k +1 P . (28)", 'B. Thermally Averaged Møller Speed': 'The Møller speed between two particles with fourmomenta p µ 1 = ( E 1 , p 1 ) and p µ 2 = ( E 2 , p 2 ) is  \nv ≡ √ ( p µ 1 p 2 µ ) 2 -( m 1 m 2 ) 2 E 1 E 2 (29)  \nwhere m 2 i ≡ 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  \nantiparallel), which we will call E and | p | . We then find  \n⟨ v ⟩ 2 = ( E 2 + | p | 2 ) 2 -m 4 E 4 = 1 + ( | p | E ) 4 +2 ( | p | E ) 2 -( m E ) 4 = 1 + ( 3 T γm ) 2 +2 ( 3 T γm ) -1 γ 4 = 6 x √ 1 -v rel + 9 x 2 (1 -v rel ) + v 2 rel (2 -v 2 rel ) . (30)  \nThe relativistic equipartition theorem of 3 T = | p | 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 .', 'C. 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  \nL ⊃ 1 2 m 2 ∆ ∆ 2 -1 2 λ ∆ 2 | ϕ | 2 (31)  \nThe 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  \n1 2 λ ∆ 2 | ϕ + ⟨ ϕ ⟩| 2 = 1 2 λ ∆ 2 ( | ϕ | 2 + ⟨ ϕ ⟩ 2 +2 ⟨ ϕ ⟩| ϕ | ) . (32)  \nHere, 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.  \nIn 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 |T ϕ → 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  \nY = Y ϕ Y ϕ fv n ϕ ( T PT ) s ( T PT ) N |T ϕ → N ∆ | 2 . (33)  \nPrior 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  \nY = Y ϕ Y ϕ fv 45 ζ (3) g ϕ 2 π 4 g ∗ ,s N |T ϕ → N ∆ | 2 . (34)  \nTo meet the condition of DM abundance to be en-  \n- [1] G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405 , 279 (2005), arXiv:hep-ph/0404175.\n- [2] E. Komatsu and C. L. Bennett (WMAP Science Team), PTEP 2014 , 06B102 (2014), arXiv:1404.5415 [astroph.CO].\n- [3] P. A. R. Ade et al. (Planck), Astron. Astrophys. 571 , A16 (2014), arXiv:1303.5076 [astro-ph.CO].\n- [4] N. Aghanim et al. (Planck), Astron. Astrophys. 641 , A1 (2020), arXiv:1807.06205 [astro-ph.CO].\n- [5] N. Aghanim et al. (Planck), Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].\n- [6] M. Ackermann et al. (Fermi-LAT), Phys. Rev. Lett. 115 , 231301 (2015), arXiv:1503.02641 [astro-ph.HE].\n- [7] A. Albert et al. (Fermi-LAT, DES), Astrophys. J. 834 , 110 (2017), arXiv:1611.03184 [astro-ph.HE].\n- [8] R. K. Leane, T. R. Slatyer, J. F. Beacom, and K. C. Y. Ng, Phys. Rev. D 98 , 023016 (2018), arXiv:1805.10305 [hep-ph].\n- [9] H. Baer, K.-Y. Choi, J. E. Kim, and L. Roszkowski, Phys. Rept. 555 , 1 (2015), arXiv:1407.0017 [hep-ph].\n- [10] G. F. Giudice, E. W. Kolb, and A. Riotto, Phys. Rev. D 64 , 023508 (2001), arXiv:hep-ph/0005123.\n- [11] A. L. Erickcek, Phys. Rev. D 92 , 103505 (2015), arXiv:1504.03335 [astro-ph.CO].\n- [12] M. Nemevsek, G. Senjanovic, and Y. Zhang, JCAP 07 , 006 (2012), arXiv:1205.0844 [hep-ph].\n- [13] M. Nemevˇsek and Y. Zhang, Phys. Rev. Lett. 130 , 121002 (2023), arXiv:2206.11293 [hep-ph].\n- [14] M. Nemevˇsek and Y. Zhang, Phys. Rev. D 109 , 056021 (2024), arXiv:2312.00129 [hep-ph].\n- [15] G. B. Gelmini and P. Gondolo, Phys. Rev. D 74 , 023510 (2006), arXiv:hep-ph/0602230.\n- [16] R. Allahverdi, B. Dutta, and K. Sinha, Phys. Rev. D 83 , 083502 (2011), arXiv:1011.1286 [hep-ph].\n- [17] M. Kawasaki, T. Moroi, and T. Yanagida, Phys. Lett. B 370 , 52 (1996), arXiv:hep-ph/9509399.\n- [18] T. Moroi and L. Randall, Nucl. Phys. B 570 , 455 (2000), arXiv:hep-ph/9906527.\n- [19] S. Iso, N. Okada, and Y. Orikasa, Phys. Rev. D 80 , 115007 (2009), arXiv:0909.0128 [hep-ph].\n- [20] S. Iso, N. Okada, and Y. Orikasa, Phys. Lett. B 676 , 81 (2009), arXiv:0902.4050 [hep-ph].\n- [21] R. Jinno and M. Takimoto, Phys. Rev. D 95 , 015020 (2017), arXiv:1604.05035 [hep-ph].\n- [22] C. Marzo, L. Marzola, and V. Vaskonen, Eur. Phys. J. C 79 , 601 (2019), arXiv:1811.11169 [hep-ph].\n- [23] P. Huang and K.-P. Xie, JHEP 09 , 052 (2022), arXiv:2206.04691 [hep-ph].\n- [24] L. Delle Rose, G. Panico, M. Redi, and A. Tesi, JHEP  \nhanced from its standard freeze-out value, we require Y ≥ Y c ( T PT ). Combining this condition with Eqs. 4 and 34 yields  \nY ϕ Y ϕ fv ⟨ σv ⟩ T PT ≥ 2 √ 6 π 5 / 2 ζ (3) M Pl √ g ∗ g ϕ 1 N |T ϕ → N ∆ | 2 . (35)  \n04 , 025 (2020), arXiv:1912.06139 [hep-ph].  \n- [25] B. Von Harling, A. Pomarol, O. Pujol'as, and F. Rompineve, JHEP 04 , 195 (2020), arXiv:1912.07587 [hep-ph].\n- [26] P. Baratella, A. Pomarol, and F. Rompineve, JHEP 03 , 100 (2019), arXiv:1812.06996 [hep-ph].\n- [27] K. Agashe, P. Du, M. Ekhterachian, S. Kumar, and R. Sundrum, JHEP 02 , 051 (2021), arXiv:2010.04083 [hep-th].\n- [28] L. Randall and G. Servant, JHEP 05 , 054 (2007), arXiv:hep-ph/0607158.\n- [29] T. Abe, K. S. Babu, and A. Kaladharan, (2023), arXiv:2312.17203 [hep-ph].\n- [30] S. Borsanyi et al. , Nature 539 , 69 (2016), arXiv:1606.07494 [hep-lat].\n- [31] K. Kaneta, Y. Mambrini, and K. A. Olive, Phys. Rev. D 99 , 063508 (2019).\n- [32] M. Fujii and K. Hamaguchi, Phys. Rev. D 66 , 083501 (2002), arXiv:hep-ph/0205044.\n- [33] B. S. Acharya, P. Kumar, K. Bobkov, G. Kane, J. Shao, and S. Watson, JHEP 06 , 064 (2008), arXiv:0804.0863 [hep-ph].\n- [34] K. Sinha, AIP Conference Proceedings 1604 , 105 (2014), https://pubs.aip.org/aip/acp/articlepdf/1604/1/105/11710526/105 1 online.pdf.\n- [35] K. Griest and M. Kamionkowski, Phys. Rev. Lett. 64 , 615 (1990).\n- [36] A. Albert et al. (HAWC), JCAP 12 , 038 (2023), arXiv:2305.09861 [astro-ph.HE].\n- [37] X.-J. Bi et al. , Chin. Phys. C 46 , 030005 (2022).\n- [38] H. Abdalla et al. (H.E.S.S.), Phys. Rev. Lett. 129 , 111101 (2022), arXiv:2207.10471 [astro-ph.HE].\n- [39] A. Acharyya et al. (CTA), JCAP 01 , 057 (2021), arXiv:2007.16129 [astro-ph.HE].\n- [40] P. A. R. Ade et al. (Planck), Astron. Astrophys. 594 , A13 (2016), arXiv:1502.01589 [astro-ph.CO].\n- [41] L. Accardo et al. (AMS Collaboration), Phys. Rev. Lett. 113 , 121101 (2014).\n- [42] M. Aguilar et al. (AMS Collaboration), Phys. Rev. Lett. 113 , 121102 (2014).\n- [43] D. Borah, S. Mahapatra, N. Sahu, and V. S. Thounaojam, Phys. Rev. D 108 , 083038 (2023), arXiv:2308.06172 [hep-ph].\n- [44] A. Adhikary, D. Borah, S. Mahapatra, I. Saha, N. Sahu, and V. Thounaojam, (2024), arXiv:2405.17564 [hep-ph].\n- [45] D. Chway, T. H. Jung, and C. S. Shin, Phys. Rev. D 101 , 095019 (2020), arXiv:1912.04238 [hep-ph].\n- [46] M. J. Baker, J. Kopp, and A. J. Long, Phys. Rev. Lett. 125 , 151102 (2020), arXiv:1912.02830 [hep-ph].\n- [47] C. Gross, G. Landini, A. Strumia, and D. Teresi, JHEP 09 , 033 (2021), arXiv:2105.02840 [hep-ph].\n- [48] M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht, (2021), arXiv:2105.07481 [astro-ph.CO].\n- [49] M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht, (2021), arXiv:2110.00005 [astro-ph.CO].\n- [50] K. Kawana and K.-P. Xie, Phys. Lett. B 824 , 136791 (2022), arXiv:2106.00111 [astro-ph.CO].\n- [51] J.-P. Hong, S. Jung, and K.-P. Xie, Phys. Rev. D 102 , 075028 (2020), arXiv:2008.04430 [hep-ph].\n- [52] P. Huang and K.-P. Xie, Phys. Rev. D 105 , 115033 (2022), arXiv:2201.07243 [hep-ph].\n- [53] K.-P. Xie, (2024), arXiv:2405.01227 [hep-ph].\n- [54] C. Caprini et al. , JCAP 04 , 001 (2016), arXiv:1512.06239 [astro-ph.CO].\n- [55] M. Maggiore, Phys. Rept. 331 , 283 (2000), arXiv:grqc/9909001.\n- [56] B. Allen and J. D. Romano, Phys. Rev. D 59 , 102001 (1999).\n- [57] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D 96 , 103520 (2017), [Erratum: Phys.Rev.D 101, 089902(E) (2020)], arXiv:1704.05871 [astro-ph.CO].\n- [58] J. Ellis, M. Lewicki, and J. M. No, JCAP 07 , 050 (2020), arXiv:2003.07360 [hep-ph].\n- [59] H.-K. Guo, K. Sinha, D. Vagie, and G. White, JCAP 01 , 001 (2021), arXiv:2007.08537 [hep-ph].\n- [60] M. Lewicki and V. Vaskonen, Eur. Phys. J. C 80 , 1003 (2020), arXiv:2007.04967 [astro-ph.CO].\n- [61] J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, JCAP 06 , 028 (2010), arXiv:1004.4187 [hep-ph].\n- [62] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D 92 , 123009 (2015), arXiv:1504.03291 [astro-ph.CO].\n- [63] T. L. S. Collaboration, (2014).\n- [64] X. Gong et al. , J. Phys. Conf. Ser. 610 , 012011 (2015), arXiv:1410.7296 [gr-qc].\n- [65] D. Shoemaker (LIGO Scientific), (2019), arXiv:1904.03187 [gr-qc].\n- [66] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32 , 074001 (2015), arXiv:1411.4547 [gr-qc].\n- [67] V. Corbin and N. J. Cornish, Class. Quant. Grav. 23 , 2435 (2006), arXiv:gr-qc/0512039.\n- [68] K. Yagi and N. Seto, Phys. Rev. D 83 , 044011 (2011), [Erratum: Phys.Rev.D 95, 109901(E) (2017)], arXiv:1101.3940 [astro-ph.CO].  \n- [69] B. P. Abbott et al. (LIGO Scientific), Class. Quant. Grav. 34 , 044001 (2017), arXiv:1607.08697 [astroph.IM].\n- [70] S. Kawamura et al. , PTEP 2021 , 05A105 (2021), arXiv:2006.13545 [gr-qc].\n- [71] S. Kawamura et al. , Class. Quant. Grav. 23 , S125 (2006).\n- [72] R. van Haasteren et al. (EPTA), Mon. Not. Roy. Astron. Soc. 414 , 3117 (2011), [Erratum: Mon.Not.Roy.Astron.Soc. 425, 1597 (2012)], arXiv:1103.0576 [astro-ph.CO].\n- [73] M. Punturo et al. , Class. Quant. Grav. 27 , 194002 (2010).\n- [74] A. Brown, A. Vallenari, T. Prusti, J. De Bruijne, C. Babusiaux, C. Bailer-Jones, M. Biermann, D. W. Evans, L. Eyer, F. Jansen, et al. , Astronomy & astrophysics 616 , A1 (2018).\n- [75] C. Caprini et al. , JCAP 03 , 024 (2020), arXiv:1910.13125 [astro-ph.CO].\n- [76] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, et al. , arXiv preprint arXiv:1702.00786 (2017).\n- [77] T. Robson, N. J. Cornish, and C. Liu, Class. Quant. Grav. 36 , 105011 (2019), arXiv:1803.01944 [astroph.HE].\n- [78] Z. Arzoumanian et al. (NANOGrav), Astrophys. J. Lett. 905 , L34 (2020), arXiv:2009.04496 [astro-ph.HE].\n- [79] T. J. W. Lazio, S. Bhaskaran, C. Cutler, W. M. Folkner, R. S. Park, J. A. Ellis, T. Ely, S. R. Taylor, and M. Vallisneri, IAU Symp. 337 , 150 (2017), arXiv:1801.02898 [astro-ph.IM].\n- [80] G. Janssen et al. , PoS AASKA14 , 037 (2015), arXiv:1501.00127 [astro-ph.IM].\n- [81] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Int. J. Mod. Phys. A 35 , 2050075 (2020), arXiv:1807.09495 [gr-qc].\n- [82] A. Vallenari, Frontiers in Astronomy and Space Sciences 5 , 11 (2018).\n- [83] J. Luo et al. (TianQin), Class. Quant. Grav. 33 , 035010 (2016), arXiv:1512.02076 [astro-ph.IM].\n- [84] E. W. Kolb and M. S. Turner, The Early Universe , Vol. 69 (1990)."}