Regurgitated Dark Matter

We present a new paradigm for the production of the dark matter (DM) relic abundance based on the evaporation of early Universe primordial black holes (PBHs) themselves formed from DM particles. As a concrete realization, we consider a minimal model of the dark sector in which a first-order phase transition results in the formation of Fermiball remnants that collapse to PBHs, which then emit DM particles. We show that the regurgitated DM scenario allows for DM in the mass range $\sim1$ GeV $- \,10^{16}$ GeV, thereby unlocking parameter space considered excluded.

Introduction.-Darkmatter (DM) constitutes ∼ 85% of all matter in the Universe, as determined by a multitude of astronomical observations (for reviews, see e.g.[1,2]).Despite decades of efforts to detect its non-gravitational interactions, the nature of DM remains mysterious.
A significant focus has been on weakly interacting massive particles (WIMPs) as the DM paradigm [3,4], with typical masses in the GeV to multi-TeV range that often appear in theories that address the hierarchy problem of the standard model (SM).The scenario of decoupling from the thermal bath in the early Universe, which successfully explains cosmological observations such as the light element abundances [5], suggests that WIMPs with typical electroweak-scale masses and annihilation cross sections can readily account for the observed DM relic abundance through thermal freeze-out (the so-called "WIMP miracle").Sensitive experimental searches (e.g.[6,7]) significantly constrain the parameter space of minimal WIMP scenarios.DM scenarios based on additional or number-changing DM particle interactions, such as the strongly interacting massive particle miracle [8], provide alternative approaches to achieving the DM relic abundance.
In this work we present a novel paradigm, which we call regurgitated dark matter (RDM), for producing the DM relic abundance.It is based on the evaporation of early Universe primordial black holes (PBHs), themselves constituted by DM particles.PBHs scramble and re-emit DM particles with altered energy-momentum and abundance distributions, distinct from that of the original DM particles in the thermal bath that formed the PBHs, i.e., these properties are not determined by direct DM interactions as in conventional DM production mechanisms.While particle DM emission from evaporating PBHs has been studied (e.g.[9][10][11]), in the RDM scenario, PBHs originate from the same DM particles that later constitute the DM relic abundance and not from some distinct mechanism; for PBH production mechanisms see e.g.[12][13][14][15][16][17][18][19][20][21][22][23][24][25].As we demonstrate, RDM can open a new window in the WIMP parameter space with either fermion or scalar DM particles.
Model.-We illustrate an elegant realization of RDM in the context of a first-order phase transition (FOPT) in an asymmetric dark sector that produces Fermiball remnants composed of dark sector particles that subsequently collapse to PBHs.The dark sector particles that form the PBHs are emitted by these PBHs through Hawking evaporation.Depending on the particle mass and PBH mass, RDM particles may be relativistic at the epoch of Big Bang nucleosynthesis (BBN) or contribute to warm DM.We note, however, that our RDM mechanism is general and may be realized in the context of other scenarios, such as the collapse of solitonic macroscopic objects to PBHs or in models with additional dark sector forces and SM portals beyond the Higgs portal.We leave the exploration of such possibilities for future work.
We consider the model of Refs.[26][27][28][29][30] given by where L SM is the standard model (SM) Lagrangian, dark sector fermions χ, χ and scalar ϕ interact via an attractive Yukawa force with coupling y χ , and ϕ interacts with the SM sector through the Higgs H doublet portal coupling κ.On receiving thermal corrections, the potential V (ϕ) becomes V (ϕ, T ) which triggers a FOPT below the critical temperature, forming Fermiball remnants that collapse to PBHs.We remain agnostic about the details of the potential, allowing a general discussion of RDM.We assume the existence of an asymmetry between the number density of χ and χ with η χ = (n χ − n χ)/s(T ⋆ ), where s(T ⋆ ) = 2π 2 g(T ⋆ )T 3 ⋆ /45 is the entropy number density with relativistic degrees of freedom (d.o.f.) g(T ) at temperature T ⋆ of the FOPT; we neglect the small contribution of the dark sector d.o.f. to g(T ).Such an asymmetry can be realized in a variety of asymmetric DM mechanisms [31][32][33], with Fermi-degenerate remnants dominated by χ.
Formation of Fermiballs.-We consider the following thermal history of cosmology for production of RDM, illustrated schematically in Fig. 1.
At first, the dark sector and SM particles are in thermal equilibrium after inflationary reheating, which can occur either due to the Higgs portal coupling κ or inflaton sector couplings.As the Universe expands and temperature decreases below the electroweak phase transition at T ∼ 160 GeV, the dark sector decouples from the visible sector at T ∼ m h = 125 GeV due to the diminishing interactions between ϕ and the Higgs.Hereafter, the SM and dark sector temperatures T SM and T , respectively, evolve separately.The effective number of relativistic degrees of freedom (d.o.f.) of the SM at dark sector decoupling is g(T dec SM ).At a temperature T ⋆ , a FOPT is triggered by the dark scalar potential V (ϕ, T ⋆ ).The phase transition changes the expectation value of ϕ, from ⟨ϕ⟩ = 0 in the false vacuum to ⟨ϕ⟩ = v ⋆ in the true vacuum.The FOPT proceeds through bubble nucleation (with expanding bubble wall speed v w ) and can be characterized by the following parameters [34]: β is the inverse duration of the FOPT, and α D ≃ 30∆V /π 2 g D T 4 ⋆ quantifies its strength.Here, g D = 4.5 is the d.o.f. of the dark sector, and ∆V is the potential energy difference between the false and true vacua.
The FOPT can readily induce a significant mass gap ∆m ϕ , ∆m χ ≫ T ⋆ between the vacua, with massless fermions acquiring mass m χ = y χ v ⋆ via the Yukawa coupling.For the particles to be trapped, the dark sector particle masses in the true vacuum must exceed the FOPT temperature: This condition can be fulfilled in supercooled scenarios with large v ⋆ /T ⋆ [35][36][37][38][39], or with y χ ≫ 1 [40,41].Then the mass of dark sector particles is significantly larger than their thermal kinetic energy and cannot penetrate into the true vacuum bubbles.As true vacuum bubbles expand, dark sector particles are efficiently trapped in contracting regions of false vacuum.In much of the parameter space, explicit calculations confirm that the trapping fraction of dark sector particles in the false vacuum remnants is ∼ 1 if Eq. ( 2) is satisfied.The remnants get compressed to form non-topological solitonic Fermiball remnants [26-30, 42, 43].We assume that the number of d.o.f.does not significantly vary between decoupling and the FOPT, which will not affect our conclusions.
Fermiball cooling.-Thedark sector temperature of the Fermiballs will be T 1 = (90∆V /π 2 g D ) 1/4 , during the slow remnant cooling process [29].The Fermiballs cool via SM particle production through the Higgs portal.A detailed analysis of Fermiball cooling for our regimes of interest can be found in Supplemental Material [44].The asymmetry η χ ensures that a population of χ particles survives after annihilation to ϕ's.As the Fermiball cools, these particles dominate until the Fermiball collapses into a black hole.
The dominant cooling channel for Fermiballs depends on the dark sector temperature T 1 .For T 1 below the electroweak scale, the cooling rate is suppressed by the Higgs mass.For small values of κ, we have verified that volumetric cooling occurs with a rate Ċ = n 2 ⟨2E⟩σv rel (see Supplemental Material [44] for details).The scalar ϕ follows a thermal Bose-Einstein distribution with number density n = (ζ(3)/π 2 )T 3 and average energy ⟨E⟩ ≃ 2.7T 1 .
In this regime, Fermiball remnants will predominantly cool through ϕϕ − → f f emission, where f is the heaviest available fermion with mass m f .The remnant is initially supported by thermal pressure, but as it cools, it becomes supported by the Fermi degeneracy pressure of the asymmetric χ population.In the volumetric cooling regime, this transition happens at a temperature, For T 1 above the electroweak scale, direct Higgs production, ϕϕ → HH, can occur.For small κ, we have volu-metric cooling with the transition temperature, For larger κ in both temperature ranges, the mean free path of the particles can become shorter than the Fermiball radius, and blackbody surface cooling dominates.In this case, the transition time is negligible, i.e., T tr SM ∼ T ⋆ (see Supplemental Material [44] for details).Hence, the resulting DM abundance is determined primarily by the black hole evaporation timescale.
Black hole formation.-Blackholes are formed from the cooling Fermiballs when the length scale associated with attractive Yukawa force ∼ 1/m ϕ is comparable with the mean separation of χ (χ) inside.In practice, this collapse occurs shortly after the transition to Fermi degeneracy pressure, at temperature T tr SM [29,30].(Note that gravitational collapse is also possible [45].)This instability is ensured for α D > 0.01 [30], which is readily satisfied for α D > 1/3 so that the initial remnant shrinks efficiently.Then, the average mass of PBHs formed from the collapse of Fermiballs is [29] M PBH ≃ (6.61 × 10 14 g) α The number density of these PBHs is with g s (T SM ) the number of entropic d.o.f. of the visible sector.
PBHs will eventually evaporate through Hawking radiation [46] at a time t PBH after formation, as reviewed in Supplemental Material [44].If PBHs come to dominate the matter density, then evaporation will reheat the Universe to a temperature, where g H,SM is the number of Hawking d.o.f. for the SM sector, which we elaborate on below.If there is no PBH-dominated era, then the prefactor in Eq. ( 7) becomes T evap SM ∼ 43 MeV.The time at which evaporation occurs is the sum of the three timescales -the PBH lifetime, the formation (cooling) time, and t ⋆ .Because of the hierarchy of scales, the evaporation occurs at ∼ min(T tr SM , T RH SM , T ⋆ ).
There are several key timescales in our setup.We are primarily interested in PBH evaporating into sufficient quantities of RDM before BBN to maintain the successful predictions of the light element abundances.This sequence includes the production of Fermiballs, their subsequent cooling and collapse into PBHs, and PBH evaporation.The relevant timescales for Fermiball formation and PBH collapse are related to the temperatures T ⋆ and T tr SM respectively, and the evaporation temperature is T RH SM (T evap SM ).To evade BBN constraints, we require that T ⋆ , T tr SM , T The relative Hawking emission rates of dark sector particles to SM particles is given by the ratio of their effective Hawking d.o.f.g H,D /g H,SM , with g H,D = 5.82 for the dark sector and g H,SM = 108 for the SM sector [49,50] (see Supplemental Material [44]).Most emitted particles have masses smaller than the Hawking temperature T PBH of the PBH.So, T PBH ≳ m χ , m ϕ ≫ T 1 , T ⋆ holds for dark sector particle emission.While both ϕ and χ are emitted, the scalar ϕ develops a mixing with the SM Higgs after the electroweak and dark sector phase transitions, allowing for the decay of ϕ into SM particles.In the relevant regions of allowed parameter space, this decay timescale is shorter than the lifetime of the Universe.Hence, in this particular model realization of the RDM paradigm, ϕ does not constitute a significant DM component.In the following, we focus on the abundance of the stable fermion χ.Although ϕ is unstable and does not contribute to the current mass density, it's parameters can still be constrained by considerations of the cooling rate and the Higgs invisible decay.In these instances, we assume m ϕ = m χ .
When the initial Hawking temperature is smaller than the particle mass, the dark sector particles are emitted once the Hawking temperature reaches their mass in the true vacuum, ϵ em T PBH ≃ m χ where ϵ em = 2.66, 4.53, 6.04 for spin s = 0, 1/2, 1 particles, respectively [51].This higher temperature corresponds to a lighter PBH mass below which these heavy particles can be emitted: The bulk of emitted particles will be nonrelativistic, with the initial density at emission reduced by a factor of (M em PBH /M PBH ), such that ρ χ /ρ SM = Ω DM = 0.264 and T1 ≃ T⋆ ≃ 0.05mχ.Also displayed are the most restrictive bounds from several experiments (labeled "Exp." in the left panel), including the 95% CL bound from the invisible Higgs decay branching fraction [47], and the 90% CL bounds from XENONnT [6] and LUX-ZEPLIN (LZ) [48].The solid black curve corresponds to conventional Higgs portal WIMP freezeout, and the solid and dashed gray curves correspond to thermal production and gravitational overproduction of superheavy WIMPZillas, respectively.In the lower shaded region labeled "Cool.",PBHs form after BBN.
From the emission spectrum of PBHs, the average energy for heavy particles is E = 2m χ .Thus, the present day DM mass density of such primarily nonrelativistic emitted DM particles is where T 0 SM = 2.73 K and g s (T 0 SM ) = 3.9 are the present day values.An additional factor of ρ PBH /ρ SM appears if the PBHs are subdominant (see Supplemental Material [44]).
If the particles are emitted relativistically, the present day mass density contribution is reduced by the redshifting of these particles until they become nonrelativistic.The bulk of the Hawking radiation emitted particles will have a Lorentz boost factor γ ≃ ϵ em T PBH /m χ .The resulting density of these initially relativistic dark sector particles is The behavior is opposite to the nonrelativistic case, with lighter particles being less dense.Dark matter detection.-Asdark matter, χ can be observed in direct detection experiments by interactions through the Higgs portal coupling in Eq. ( 1).The resulting elastic scattering cross section of χ on nucleons N is given by [47,52] (see Supplemental Material [44] for details) where m N ≃ 1 GeV is the nucleon mass and f N ∼ 0.3 is Higgs-nucleon interaction parameter.Employing Eq. ( 11), we recast existing bounds on the spin-independent scattering cross section into constraints on the coupling κ.In Fig. 2, we display constraints on κ (upper shaded region) as well as predictions (solid colored lines) for χ to constitute all of the DM abundance for different values of β/H and η χ .In the lower shaded region labeled "Cool.",PBHs form after BBN via Fermiball cooling.Clearly, the regurgitated χ can saturate the DM relic abundance for a wide range of masses.The final abundance has a κ dependence if the cooling time is longer than the lifetime of the PBH and the transition time t ⋆ .In the WIMP mass range, the cooling rate depends on fermion channels that open sequentially with increasing m χ = 20T 1 .Here, g(T ⋆ ) also changes similarly and affects the PBH number density through Eq. (6).
The Yukawa interaction between χ and ϕ keeps χ in thermal equilibrium after it decouples from the SM.Therefore, we calculate χ's freeze-out abundance by as-suming the two particles freeze-out together when ϕ decouples from the SM plasma.We interpret the results for thermal freeze-out production of Ref. [53] for ϕϕ → f f in the nonrelativistic limit, (see Supplemental Material [44]), and plot the result as the solid black curve in Fig. 2. While WIMP masses are strongly constrained, RDM can be efficiently produced in the unconstrained parameter space below m χ ∼ TeV.
For m χ < m h /2, stringent bounds on κ arise from the invisible Higgs decay branching fraction to DM particles Br(H → inv) constrained by Large Hadron Collider (LHC) data at 95% confidence level [47].We show this bound in Fig. 2 with m ϕ = m χ .Direct detection constraints on κ weaken for heavier m χ since In the mass range 10 9 GeV ≲ m χ ≲ 10 14 GeV, gravitational overproduction of WIMPZillas [54][55][56][57] can restrict χ from being a viable DM candidate depending on the Hubble rate H e at the end of inflation [58]; see the gray dashed lines in Fig. 2.However, if inflation occurs at a lower energy scale the abundance of WIMPZillas can be significantly suppressed.
If the WIMPZilla has additional interactions, a thermal relic can be realized.For the Higgs portal scenario [58], we show the case of H e = 10 13 GeV with the reheating temperature T RH SM = 10 12 GeV.For even heavier DM masses, m χ ≳ 3 × 10 16 GeV, stringent bounds originate from DEAP-3600 and mica searches [59], but they fall in a region with no reliable scaling relation.Furthermore, a robust theoretical bound for pointlike DM [60] cannot be meaningfully applied to κ, since this would imply that the nuclear scattering cross section of χ plateaus to a maximum value even for an arbitrarily large κ.
Gravitational waves (GWs) from the FOPT can provide a correlated signature with DM detection in scattering experiments.
Conclusions.-We proposed a novel paradigm of regurgitated DM, stemming from the emission of evaporating PBHs formed from the DM particles themselves.This is distinct from conventional particle DM production mechanisms, since the resulting DM relic abundance is not determined by particle interactions.Intriguingly, as we demonstrate with a concrete realization, this paradigm can produce the inferred abundance of DM in a very broad mass range ∼ 1 GeV − 10 16 GeV, and opens up parameter space previously thought to be excluded.A stochastic background of GWs is a possible correlated signature of the scenario.

SUPPLEMENTAL MATERIAL Regurgitated Dark Matter
TaeHun Kim, Philip Lu, Danny Marfatia, Volodymyr Takhistov We provide additional details of RDM production.Specifically, we discuss Hawking evaporation, Fermiball cooling rates, the conditions for PBH domination, and the interactions after the FOPT.

PBH EVAPORATION
The lifetime of a PBH depends sensitively on M PBH as t PBH ∝ M 3 PBH , with a PBH of initial mass ∼ 4 × 10 8 g evaporating around the time of BBN.The Hawking temperature of a PBH is [46] T PBH = 1.06 × 10 5 GeV M PBH 10 8 g −1 . (S1) The black hole emission further depends on the spin of the particles, with the relative rates of emission per d.o.f. with respect to spin-1/2 fermions given by [49,50] where w i is the number of spin states of each particle.The total emission rate is The SM d.o.f.contribute a total of g H ≃ 108.Using Eq. (S2), the dark sector fermions (four spin-1/2 d.o.f.) and scalar (one spin-0 d.o.f.) contribute ≃ 5.82, resulting in dark sector emission contributing approximately ∼ 5.1% of the total emission once the PBH temperature becomes sufficiently high for efficient emission of the dark sector particles.From Eq. (S1), the shape of the integrated emission spectrum can be obtained.Differentiating, dM PBH /dT PBH ∝ T −2 PBH so that EdN/dE ∝ T −2 PBH or dN/dE ∝ T −3 PBH for the high energy tail.The scaling relations for the emitted particle number density are given by [51] with M i being the initial mass of the PBH.We can construct an approximate integrated spectrum by connecting the two scaling relations at the (instantaneous) peak energy E = ϵ em T PBH .The spectrum can be normalized to the PBH mass, as each primary DM particle (before decay) should represent g H,D /(108 + g H,D ) of the total emission.For dark fermions, we use ϵ em = 4.53 and g H,D = 5.82.Then, We make the common assumption that the emission takes place nearly instantaneously on cosmological time scales.

FERMIBALL COOLING
Here we provide additional details for Fermiball cooling through the Higgs portal.For energies and masses below the electroweak scale, the cross section for ϕϕ − → f f through a Higgs propagator is where √ s is the center of mass energy.The cross section is dominated by the heaviest fermion kinematically allowed.With the hierarchy m f ≪ m ϕ ≪ m h , we find In our freeze-out calculations, we take the limit of heavy nonrelativistic ϕ.
(S9) Here, ln(R 1 /R tr ) −1/2 is an O(1) factor that depends on the ratio of the initial Fermiball radius R 1 to the transition time radius R tr approximately as R 1 /R tr ∼ η −1/3 χ .For the case when T 1 is above the electroweak scale, we have direct Higgs production through ϕϕ − → HH with cross section The corresponding cooling rate can be calculated as which yields the corresponding transition temperature, However, for large κ the mean free path of the Higgs becomes shorter than the size of the Fermiball.Surface cooling happens when nσR 1 ∼ O(1).Above the electroweak scale, this occurs for The cooling approaches blackbody radiation so that the transition temperature becomes [30] T tr SM ≃ (7.29 TeV) ) This can also occur below the electroweak scale for sufficiently large κ.Under these approximations, the transition temperature can seem to be higher than the phase transition temperature.In practice, this means the cooling is very rapid and the transition happens within a Hubble time.The PBH is quickly formed, and the transition temperature is ∼ T ⋆ .
A subdominant population of PBHs will emit a correspondingly smaller fraction of DM particles.The density in the case of a subdominant PBH population is similar to the PBH dominant case, but with an extra factor of Eq. (S21).For the interesting mass range from GeV to TeV, the masses of the dark sector particles are always below the Hawking temperature of the evaporating PBHs.The resulting DM density is obtained by combining of Eqs. ( 10) and (S21):

INTERACTIONS AFTER THE DARK SECTOR AND ELECTROWEAK PHASE TRANSITIONS
After the FOPT and electroweak symmetry breaking, the fields can be expanded around the new minima as ϕ → ϕ + v ⋆ and H = (0, (v h + h)) T / √ 2, giving the interaction Lagrangian, The masses include corrections arising from the vacuum expectation values, the Higgs mass includes cancellations from other contributions, and m χ = y χ v ⋆ .The fields ϕ, h, and χ are in the flavor basis.The bilinear term in ϕ and h implies a mixing of the form, where φ and h are the mass eigenstates.The Lagrangian is diagonalized for the mixing angle Hereafter, we assume that naturalness dictates m ϕ,0 ∼ v ⋆ and m h,0 ∼ v h .Then, θ in Eq. (S25) is suppressed by the hierarchy between v ⋆ and v h .Since κ < √ 4π by unitarity, θ ≪ 1 except for m ϕ,0 ≃ m h,0 , a case that we simply discard.Then, the masses of the mass eigenstates are m ϕ ≃ m ϕ,0 ∼ v ⋆ and m h ≃ m h,0 ∼ v h , and the leading order Lagrangian with all self-interaction terms neglected is where the mixing angle is We discuss several relevant interactions between the Higgs, ϕ and χ.In doing so, we neglect any process involving two or more particles in the initial state, assuming that they are too sparse for such an interaction to take place.

Decay of ϕ
We first calculate the decay rate of ϕ and check whether it can be a stable DM candidate.For m ϕ > 2m h , the leading process is the direct decay into a Higgs pair via the interaction term, giving the decay rate Assuming this decay channel is dominant, the lifetime of the ϕ particle is (S30) This imposes strong constraints on the Higgs portal coupling.Even the minimum requirement of τ ϕ ≳ t U ≈ 14 Gyr gives κ ≲ 10 −22 for m ϕ = 10 3 GeV and κ ≲ 10 −25 for m ϕ = 10 9 GeV.Such small values of κ do not give sufficient cooling rates for ϕ to constitute a significant component of RDM.
For m ϕ ≤ 2m h , decay into Higgs pairs is not possible and the dominant channel for ϕ decay is into SM fermions through their mixing with the flavor eigenstate h.For example, the decay rate into quarks is [71]

FIG. 1 .
FIG. 1. Cosmological thermal history of RDM production.The dark sector particles in the Fermiball are re-emitted at a higher temperature after black hole formation.

9 )FIG. 2 .
FIG.2.Constraints on the Higgs portal coupling κ as a function of mχ (assuming m ϕ = mχ) over a wide mass range[left]   and for the WIMP mass range[right].On the solid iso-ΩDM = 0.264 contours, χ constitutes the entire DM abundance for specific choices of the inverse duration of the FOPT β/H and dark sector fermion asymmetry ηχ.The PBHs dominate the Universe's energy density on the blue contours and are subdominant on the red contours.We fix αD = 1, vw = 0.67, gD, * = 4 and T1 ≃ T⋆ ≃ 0.05mχ.Also displayed are the most restrictive bounds from several experiments (labeled "Exp." in the left panel), including the 95% CL bound from the invisible Higgs decay branching fraction[47], and the 90% CL bounds from XENONnT[6] and LUX-ZEPLIN (LZ)[48].The solid black curve corresponds to conventional Higgs portal WIMP freezeout, and the solid and dashed gray curves correspond to thermal production and gravitational overproduction of superheavy WIMPZillas, respectively.In the lower shaded region labeled "Cool.",PBHs form after BBN.
RH SM ≥ 5 MeV.The condition for PBH domination before they evaporate is M PBH n PBH (T RH SM ) > ρ SM .Regurgitated dark matter.-Thedark sector particles are regurgitated during PBH evaporation and can constitute the main component of DM either in the WIMP mass range, 1 GeV ≲ m χ ≲ 1 TeV, or the superheavy mass range with m χ ≫ 1 TeV.The energy density of RDM depends on the PBH mass fraction at evaporation.If the PBHs evaporate during the radiation dominated era, the density of regurgitated particles is suppressed by ρ PBH /(ρ PBH + ρ SM ) ≃ ρ PBH /ρ SM .The ratio is unity if PBHs dominate the matter density.