Leak-in Dark Matter

We introduce leak-in dark matter, a novel out-of-equilibrium origin for the dark matter (DM) in the universe. We provide a comprehensive and unified discussion of a minimal, internallythermalized, hidden sector populated from an out-of-equilibrium, feeble connection to the hotter standard model (SM) sector. We emphasize that when this out-of-equilibrium interaction is renormalizable, the colder sector undergoes an extended phase of non-adiabatic evolution largely independent of initial conditions, which we dub “leak-in.” We discuss the leak-in phase in generality, and establish the general properties of dark matter that freezes out from a radiation bath undergoing such a leak-in phase. As a concrete example, we consider a model where the SM has an out-of-equilibrium B − L vector portal interaction with a minimal hidden sector. We discuss the interplay between leak-in and freezein processes in this theory in detail and demonstrate regions where leak-in yields the full relic abundance. We study observational prospects forB−L vector portal leak-in DM, and find that despite the requisite small coupling to the SM, a variety of experiments can serve as sensitive probes of leak-in dark matter. Additionally, regions allowed by all current constraints yield DM with self-interactions large enough to address small-scale structure anomalies. ARXIV EPRINT: nnnn.nnnnn a rX iv :1 90 9. 04 67 1v 1 [ he pph ] 1 0 Se p 20 19


Introduction
Despite overwhelming gravitational evidence for the existence of dark matter (DM), the particle properties of DM remain mysterious. Historically, one of the best-motivated candidates for particle DM has been a weakly-interacting massive particle (WIMP), or more generally, a particle that was in thermal equilibrium with the standard model (SM) plasma in the early universe, but froze out as number-changing interactions with the SM, e.g. annihilations DM DM → SM SM, departed from equilibrium. One major appealing feature of this class of models is that the DM relic abundance is directly tied to its couplings to the SM, giving rise to definite and accessible experimental targets. Owing to the spectacular success of experiments searching for DM-in direct, collider, and indirect probes-the surviving WIMP parameter space is rapidly shrinking. Other scenarios for the origin of particle dark matter, and their resulting experimental signatures, are thus of high interest.
One broad and generic scenario for the origin of DM is that its relic abundance can be determined by interactions within an internally thermalized hidden sector (HS), with minimal direct involvement of SM fields [1][2][3][4][5][6][7][8][9]. Such self-interacting hidden sectors open many avenues for addressing longstanding mysteries in both particle and astrophysics, and can predict qualitatively novel signatures. More broadly, internally thermalized hidden sectors are a simple and generic possible source for the DM of our universe, and it is worth addressing in some generality how the possible cosmological origin stories for such hidden sectors impact the dynamics and signatures of the DM they produce.
One minimal and predictive way to populate a thermal dark radiation bath in the early universe is by producing it directly from interactions with the SM radiation bath. In the simplest scenarios, these interactions are sufficiently strong to bring the hidden sector into thermal equilibrium with the SM. In this paper we will focus on the regime where the leading interaction between the two sectors never enters equilibrium.
In this scenario, feeble interactions allow energy to leak from the hot SM radiation bath into the colder hidden sector. When the leading interaction is non-renormalizable, the energy injection from the SM rapidly becomes negligible as the universe expands. The population of the hidden sector is thus dominated by a limited span of UV temperatures, after which the hidden sector evolves adiabatically [10]. By contrast, when the leading interaction is renormalizable, energy injection from the SM becomes more and more important as the universe cools. In this latter case, the hidden sector radiation bath undergoes an extended phase of non-adiabatic evolution that we dub "leak-in," which realizes a quasi-static equilibrium between the energy injection from the SM and the dilution from the expansion of the universe. The aim of this paper is to investigate the properties of DM that freezes out of a hidden radiation bath in this quasi-static leak-in phase, which we dub "leak-in dark matter" (LIDM). This scenario is distinct from freezein DM [11,12], where the DM itself is the hidden particle produced from the SM. The primary difference for leak-in is that the hidden sector is internally thermalized and acquires its own temperature, which fixes the abundances of particles, including dark matter, within the hidden sector.
Despite the feeble coupling to the SM, there are many potential experimental handles on leak-in dark matter. In particular, while LIDM annihilation cross-sections are typically suppressed relative to standard WIMP benchmarks, indirect detection signals are still within reach of a variety of cosmic ray experiments, such as Fermi, AMS-02, H.E.S.S., HAWC, CTA, and others [13][14][15][16][17]. Additionally, observations of the cosmic microwave background (CMB) can place stringent constraints on LIDM annihilations during recombination. Direct detection can also be a promising avenue for detecting LIDM, with complementary sensitivity to indirect detection, and XENON1T [18] is currently probing the edge of LIDM parameter space in the benchmark model we will consider later in this work. There can also be meaningful constraints on the mediator itself, from, for example, stellar cooling or fifth force experiments [19]. Additionally, regions of the LIDM parameter space realize sizable DM selfinteraction cross-sections. Very large DM self-interactions are constrained by dwarf structure and ellipticity, but somewhat smaller self-interactions may be favored by various small-scale structure anomalies [20].
We begin by discussing the general properties of a radiation bath populated by out-of-equilibrium renormalizable interactions in Sec. 2. We establish the general properties of DM that freezes out during the resulting leak-in phase in Sec. 2.3. In Sec. 3, we introduce a concrete model of a minimal hidden sector, consisting of a feebly-coupled B − L vector boson together with dark matter, and discuss the mechanisms governing the DM relic abundance in detail. Sec. 4 examines the observable signals of the model, with the viable regions of parameter space collected in Sec. 4.5. We conclude in Sec. 5. Appendices include criteria for attaining internal thermalization in the hidden sector in App. A, some B − L model-building considerations in Sec. B, and details of the energy transfer between SM and hidden sectors in App. C.

Leak-in: the out-of-equilibrium population of a hidden radiation bath
We begin by discussing the out-of-equilibrium population of a dark radiation bath from the SM in some generality. Throughout this work, we will denote hidden sector (SM) quantities with (without) a tilde. The Boltzmann equations describing the temperature evolution of two internally thermalized radiation baths areρ where C E is the collision term describing the energy transfer between sectors, ρ andρ are the SM and hidden sector energy densities, respectively, and we have made the simplifying approximation of neglecting the contribution of non-relativistic species to the energy density. Assuming that interactions within each sector keep the sectors in internal thermal equilibrium at separate temperatures T,T , these equations can be solved to obtain the dependence of T andT on the scale factor a.
The form of the collision term, and in particular its dependence on T andT , depend on the nature and structure of the leading interaction linking the two sectors. When the leading interaction is renormalizable, the collision term falls off more slowly with temperature than the Hubble term: in other words, scattering through renormalizable interactions becomes more important in the late universe relative to the early universe. This IR-dominance has the useful consequence of making the properties of hidden sectors populated through a renormalizable interaction with the SM relatively insensitive to the unknown reheating temperature of the universe.

Populating a cold sector through renormalizable interactions
There is a limited suite of possible renormalizable operators that allow SM particles to interact with a particle that is a total SM singlet. These operators include a dark fermion ψ coupling through the neutrino portal, O ν = ψHL, a dark scalar S coupling through the Higgs portal, O h = µS|H| 2 + S 2 |H| 2 , a dark vector boson Z D coupling through kinetic mixing with hypercharge, O Y = B µν Z µν D , and a dark vector boson Z D coupled to the SM through one of the anomaly free currents: either Each of these interactions together with those of the SM generate tree-level 2 → 2 scattering processes, which are, at zero temperature, independent of the dark particle mass in the E CM m limit. Dimensional analysis then suggests that the scattering rate in the early universe should scale like T , for T m, an expectation borne out in explicit kinetic theory calculations [21]. However, properly accounting for the contribution to the thermal self-energies in the dense medium of the radiation bath can in some cases parametrically alter this expectation [19,22,23]. In particular, when a dark species X couples to the SM entirely through mixing with another state in the plasma, there is a parametric suppression of the production rate of X from the SM plasma as m X /T → 0 [19]. If X can mix with a SM state A in medium, the propagating degrees of freedom can be found by diagonalizing the 2 × 2 propagation matrix, where Π IJ are (1PI) thermal self-energies, and for simplicity we have taken A to be massless (in many examples of interest it is the photon). The observation of Ref. [19] is that the off-diagonal entries of this matrix provide important corrections to the finite temperature propagator, and therefore to the net production rate of X from the SM plasma. Taking X to be coupled to the SM plasma through a parameter 1, we can write Π XX ≡ Π XX dk + Π XX SM , where Π XX SM is of order O( 2 ) and the O( 0 ) piece of the self-energy, Π XX dk , accounts for possible contributions from other dark species that may be in the plasma (with no direct coupling to the SM). Working to leading order in and absorbing Π XX dk into an effective mass for the dark state,m 2 X , the eigenmodes of Eq. 2.4 are (2.6) The production rate of Xs from the SM plasma is then given in terms of the imaginary part of this self energy, (2.7) Now, suppose that X inherits all its couplings to the SM from mixing with A. Then we can write (to lowest nontrivial order in ) With this relation, the term in parentheses in Eq. 2.7 can be expanded in the Π AA m 2 X limit to obtain which is directly proportional tom 2 X and vanishes in them/T → 0 limit, contrary to the naive expectation from kinetic theory, Γ ∝ T .
On the other hand, if the tight relationship of Eq. 2.8 doesn't hold, so that Π AA = C 1 , Π AX = C 2 , Π XX SM = 2 C 3 for generically ∼ O(1) differences between the various C i , then the cancellation of the leading terms in Eq. 2.9 does not occur, and the generic scaling Γ ∼ T does hold. Thus one expects the cosmological production rate of a kinetically-mixed dark photon in the early universe to be parametrically different from that of a B − L gauge boson, which has a distinct coupling structure from the photon. Another interesting case is a dark Higgs boson S [24], which can have unsuppressed thermal production in two ways. The interaction Lagrangian coupling S to the SM Higgs does directly give S unique couplings to the SM plasma through its interactions with the Higgs. Once the Higgs boson leaves the plasma shortly after electroweak symmetry breaking, the dark Higgs inherits all of its couplings to species remaining in the SM plasma from mixing with H; however, as the Higgs boson itself is gone, the SM production of S can still be unsuppressed.
Our focus in this paper will be on the case where the leading thermal scattering rates between HS and SM particles are unsuppressed in the m/T → 0 limit. This is partly for computational tractability, as it makes kinetic theory calculations a reliable guide to the temperature-dependence of the theory, and partly because these scenarios allow us to reveal some novel cosmological behavior. In these models, the production of dark states from the SM is dominated by 2 → 2 scattering, e.g. f g → Xf , where, for example, X may be a B − L dark vector boson or a Higgs-mixed dark scalar. Importantly, this particle X is not the dark matter, and is typically unstable on cosmological time-scales. These 2 → 2 processes have scattering rates that generically scale as Γ ∼ n σv ∝ T , and dominate the interactions between sectors when T m X [21]. The collision term C E describing the energy transferred between the two sectors through these scattering processes is given by the thermal average of the scattering amplitude weighted by the energy carried by the dark particle. If SM particles 1 and 2 scatter to SM particle 3 and a dark particle, labeled 4, the collision term can be expressed as Here in the last line we have introduced separate notation for the collision term governing forward scattering, C f E , which deposits energy into the dark sector, and the backward scattering term C b E , which transfers energy from the hidden sector back to the SM. In evaluating these collision terms, we will use classical (Maxwell-Boltzmann) statistics for simplicity. As these collision terms are important at energies where T m i for all particles involved, a priori the full dependence on quantum statistics should be retained. Fortunately, dark mediator production from the SM thermal bath is typically dominated by semi-fermionic processes such as f g → Xf , for which empirically we find that  Figure 1: Evolution of SM temperature (red) and hidden sector temperature (purple, green, blue) as a function of scale factor a when the two sectors are linked by the renormalizable 2 → 2 interaction of Eqs. 2.11 -2.12. Different HS temperature solutions follow from different initial conditions at a 0 . The grey line shows the attractor 'leak-in' solution,T ∝ a −3/4 , of Eq. 2.14. Solutions with initial conditions below the leak-in attractor rapidly converge to it, while the solution that starts at a higher temperature than the leak-in solution redshifts until it matches onto the attractor. The two sectors equilibrate near a = 0.1.
Maxwell-Boltzmann statistics provide a reasonable approximation to the full result, accurate to within a factor of 2 (see also [21,25,26]). To gain some quantitative intuition for the behavior of a leak-in sector, consider a toy model where the leading process transferring energy between the SM and HS radiation baths is described by a constant matrix element M = , neglecting all particle masses. In this case, the collision term describing forward energy transfer takes the simple form while the backward energy transfer from the reverse process is well-approximated as The resulting temperature evolution for both sectors is shown in Fig. 1, in the approximation that the SM dominates the Hubble expansion. Before the two sectors equilibrate, the hidden sector temperature exhibits a characteristic 'leak-in' phase, which realizes a quasi-static equilibrium between the energy injection from the SM and the energy dilution from the expansion of the universe. Hidden sectors that have a small initial reheat temperature rapidly rise up to reach the leak-in solution, as seen in the green and blue curves in Fig. 1. Meanwhile, if the hidden sector has a reheat temperature higher than the temperature of the leak-in phase, as for the purple curve in Fig. 1, it redshifts like a standard adiabatic radiation bath (T ∝ 1/a) until its temperature reaches the leak-in solution, at which point the energy injection from the SM is no longer negligible. The leak-in phase is thus an attractor solution, and in particular, at any given value of the SM temperature, the hidden sector temperature during leak-in is completely determined by the resulting energy transfer rate. This cosmology is thus IR-dominated, i.e., once the leak-in phase is attained, there is no remaining dependence on the initial conditions in the hidden sector. This ensures that the properties of DM freezing out during the leak-in phase are independent of T RH .

Essential properties of the leak-in phase
We can obtain several useful properties of the leak-in phase by analytically solving Eqs. 2.1-2.3 in the regime whereT T , and therefore • the energy of the universe is dominated by the SM radiation bath, H ∝ T 2 /M P ; • we can neglect the backward energy transfer rate into the SM; and • the SM entropy is approximately conserved, T ∝ 1/a.
With these assumptions, the evolution of the hidden sector energy density is given simply bẏ Let us now take C f E = c E T 5 , where c E is a numerical constant, as dimensional analysis requires when all masses are negligible. This will be a good approximation to the collision term in realistic models away from mass thresholds. In the toy model of Eqs. 2.11 and 2.12, c E = 2 /(64π 5 ). With ρ = (π 2 /30)g * (T )T 4 , Eq. 2.13 can be easily solved to obtaiñ (2.14) This expression for the HS temperature lets us observe two important things. Firstly,T ∝ T 3/4 ∝ a −3/4 -the HS radiation bath redshifts as if it were matter, and in particular dilutes less slowly than an adiabatic radiation bath. Secondly, the HS temperature is completely determined by the SM temperature and the strength of the leak-in interaction c E , so that it scales with the small portal coupling asT ∝ c 1/4 E ∝ 1/2 . It is worth noting that the scalingT ∝ 1/2 requires only that C E (T,T ) ≈ C f E (T ), i.e., it does not depend on the specific functional dependence on the SM temperature T , but is a direct consequence of taking the hidden sector cold compared to the SM. WhenT T , the Boltzmann equation describing the hidden sector evolution, Eq. 2.2, can be written as where both C f E and H are functions of the SM temperature only in this limit. But then, as C f E ∝ 2 , it is clear that all dependence can be scaled out by sendingT →T / √ .

Leak-in, freezeout
We would now like to understand what happens to DM that freezes out of a hidden sector radiation bath during a leak-in phase. As the leak-in phase is an attractor solution, freezeout during leak-in is a generic possibility, and does not require any fine-tuning of mass scales. We will begin with some analytic estimates to establish the main features of dark matter freezeout from a leak-in phaseor, for short, "leak-in dark matter" (LIDM) and highlight how it differs from a thermal relic in an adiabatic hidden sector. As a warm-up, we begin with a reminder of DM freezeout in a decoupled, but adiabatic, hidden sector [5], i.e., hidden sector freezeout whereT = ξT for a constant ξ. The sudden freezeout approximation, n(x f ) σv = H(x f ), lets us estimate the DM relic abundance as a function of its annihilation cross-section. Here we have definedx = m/T and x = m/T , with m the DM mass. Assuming that the SM energy dominates Hubble, the sudden freezeout approximation implies where H(m) is evaluated at x = 1 and d χ is the number of degrees of freedom for the dark matter. Iteratively solving this equation forx f yields the approximate solutionx f = ln A+ 1 2 ln ln A. To facilitate comparison with the canonical WIMP it is convenient to consider the yield Y ∞ ≈ n(x f )/s(x f ), where x f ≡ ξx f is the value of the SM temperature at DM freezeout, Necessarily this reduces to the standard result when ξ → 1. To leading order, obtaining the correct relic abundance for DM freezing out in a cold adiabatic HS requires the annihilation cross-section to be rescaled as σv → ξ σv , as the dependence ofx f on ξ is only logarithmic. Now let us repeat this exercise for DM freezing out of a leak-in radiation bath. In this case, we can read off from Eq. 2.14 thatx is related to x through which lets us express the (SM-dominated) Hubble rate in terms ofx, The sudden freezeout condition forx f is then which has the approximate solutionx f = lnÃ + 7 6 ln lnÃ. Comparing Eqs. 2.16 and 2.20 we can recognize that b is serving as a "coldness" parameter analogously to a fixed constant ξ, while the different fractional power ofx f reflects the different temperature evolution with redshift. However, b is not given by the temperature ratio between the two sectors at freezeout, which is rather ξ(x f ) = bx Using Eq. 2.20, we can derive the yield (recall that the SM entropy is approximately conserved during the leak-in phase). Thus we can again parametrically expect σv = bx where σv W denotes the annihilation cross-section for a standard thermal WIMP. In particular note the annihilation cross-section necessary to obtain the desired relic abundance scales with as b ∝ 2/3 . Since here DM freezeout depends on the hidden sector temperature, we typically expectx f ∼ 15, while x f 1 is possible.

Region of interest for the portal coupling
The LIDM mechanism can account for the observed DM abundance for a bounded range of portal couplings . At sufficiently large values of the dimensionless coupling , the dark radiation bath will thermalize with the SM, yielding a WIMP next door [21]. We refer to this transition as the "equilibration floor." For small enough , however, the dark sector never reaches a high enough co-moving dark matter number density to account for the observed DM relic abundance. This "absolute coupling floor" for leak-in dark matter can be straightforwardly estimated by requiring that the maximum value attained by the equilibrium leak-in DM yield should equal the observed relic abundance, where Ω DM ρ c,0 is the present-day energy density of DM and s 0 is the present-day entropy of the CMB. UsingT ∝ a −3/4 and T ∝ a −1 , maximizing Y eq (a) with respect to a tells us that the maximum yield is obtained atT (a max ) = 2m 5 . (2.25) Using Eq. 2.18, the maximum equilibrium yield obtained is then Requiring that this maximum yield is greater or equal to the observed relic abundance, Eq. 2.24, places a condition on the strength of the interaction with the SM, If c E is below this critical value, then even if the sector were to internally thermalize, there would never have been a large enough dark matter number density to correspond to the observed relic abundance today. Recalling c E ∝ 2 , we immediately observe that the absolute minimum value of consistent with the leak-in scenario is independent of the DM mass (although logarithmic dependence on the mass may enter through the collision term). This requirement defines an absolute coupling floor, below which leak-in cannot produce enough dark matter to reproduce the observed relic abundance. Of course, within any given model, the portal coupling will be subject to many terrestrial, astrophysical, and cosmological constraints that depend on the specific properties of the mediator X. Cosmological constraints on the mass and lifetime of X arise due to the relic hidden sector radiation bath in the early universe, which can lead to constraints through either its gravitational imprint on the early universe or the decays of X into the SM.
It is worth noting that these cosmological constraints have some model-dependence, even under the assumption that X is the lightest species in the dark sector. In the absence of other dark states that X can interact with, cosmological constraints on are dominated by "freezein" constraints on X, i.e., constraints on the out-of-equilibrium population of mediators produced in the early universe thanks to their couplings to the SM. This population is dominated by the production of Xs at SM temperatures T ∼ m X . However, when X is part of a larger dark sector that was once in internal thermal equilibrium, there is a separate population of Xs resulting from the relic radiation bath. After the HS bath leaves equilibrium, the freeze-in population will not be able to equilibrate with the relic bath population. The hidden radiation bath may leave equilibrium long before the late-time injection of freezein X, or-depending on the dark particle content-possibly not until freezein has effectively stopped. These two different scenarios lead to two very different phase space distributions for the final X population, and thus to different potential signatures. In minimal models, such as the one we will discuss below, the number density of X in a relic bath population can easily exceed the number density in the freeze-in population, and therefore may dominate any constraints arising from the decays of X. This is one example of a general theme: once we depart from thermal equilibrium, the details of which processes go out of equilibrium first can lead to rich behavior even within a simple model, e.g. [24,[27][28][29][30][31][32][33][34][35][36]. To go further and work out the observational consequences for leak-in DM, we will need to be more concrete and specify a model. In the next section, we build on this discussion, extending this toy model of leak-in DM to a more complete picture of dark matter production in a specific out-of-equilibrium hidden sector.

Dark matter relic abundance in an out-of-equilibrium hidden sector
In the previous section we developed a general analytic guide to the properties of DM that freezes out during the "leak-in" evolution of a hidden radiation bath, which we refer to as leak-in DM. The same interactions that separately govern leak-in and freezeout will typically also yield out-of-equilibrium production of DM directly from the SM, i.e., freezein [11,12]. Although direct production of DM from the SM will generally give a sub-leading contribution to the total energy density of the HS radiation bath, it has the potential to substantially affect the final DM number density, and thus the relic abundance. The leak-in mechanism is controlled by the HS temperature, and is governed by the properties of the HS radiation bath atT ∼ m/15. Freezein production, on the other hand, is dominated by SM temperatures near the dark matter mass, x f i ∼ 3 − 5. Depending on the coldness of the HS relative to the SM, direct production of DM from the SM may thus dominantly occur either prior to HS freezeout, i.e., x f i < x f , in which case its ultimate impact is negligible, or post-HS freezeout, x f i > x f , in which case it can sometimes, but not always, dominate the final DM abundance (see Fig. 2).
In this section, we introduce a specific model of a minimal hidden sector for concreteness, consisting of Dirac fermion dark matter χ coupled to a (massive) dark vector Z D that couples to the SM through B − L charges. We detail the resulting interplay of leak-in, freezein, and "reannihilation" [28] in determining the DM relic abundance when the coupling between Z D and the SM fields via the B − L interaction is too small to allow the dark sector to achieve equilibrium with the SM. We assume throughout that the dark sector is internally thermalized; criteria for attaining internal thermalization are discussed in Appendix A.

A minimal B − L vector portal leak-in hidden sector
We consider a minimal hidden sector consisting of a Dirac fermion DM candidate, χ, and a massive dark vector, Z D . This dark vector is the gauge boson for a U (1) symmetry, and interacts with the SM by coupling to the B − L current [37][38][39] where Q f is ±1 for leptons and ±1/3 for quarks. For simplicity and minimality, we consider a Stückelberg origin for the dark vector mass [40,41]. Since we are interested in dark sectors that never attain thermal equilibrium with the SM, the B − L portal coupling is assumed to be very small. This model thus assumes a large hierarchy between the couplings of Z D to DM and to the SM, g D , which, while technically natural, does invite model-building questions. This hierarchy of couplings could originate from (e.g.) dark matter with a very large B − L charge, or from a U (1) B−L × U (1) D symmetry broken at a higher scale. In principle, UV model-building can introduce some modeldependence through the introduction of new particles in the UV. To insulate the discussion from this UV sensitivity, we will simply take g D throughout the discussion of the next two sections, but in Appendix B, we will provide some simple UV completions to this hierarchical B − L model and discuss their consequences. 1  This minimal hidden sector can be described by four independent free parameters, which we will take to be α D , , m χ and m Z D . So long as m χ 10 m Z D , such that Z D is relativistic at the time of DM freezeout, the DM relic abundance will be largely insensitive to the dark vector mass: both the DM annihilation cross-section (discussed below) and the temperature evolution of the radiation bath prior to and during freezeout are largely independent of the dark vector mass when the dark vector is relativistic. Throughout this paper, we will thus consider m Z D ≤ m χ /10 in order for this specific minimal hidden sector to serve as a useful illustration of the dynamics of a general dark sector in a leak-in phase.
We compute the energy transfer collision term, C E , by considering processes that produce dark vectors from interactions with the SM plasma. In particular, we sum up the contributions from gf → [21]. With the collision term in hand, the hidden sector temperature can be determined numerically. A particularly useful function is the ratio of hidden sector and SM temperatures, ξ =T /T . If the transfer of energy out of the hidden sector is negligible,

Interplay of leak-in dark matter and freezein processes
In addition to leak-in, below the equilibration floor there are two related processes that can govern the relic abundance, freezein [11,12] and reannihilation [28]. "Freezein" refers to an out-of-equilibrium dark matter population injected predominantly near T ∼ m χ /(2 − 5) with little subsequent evolution, while "reannihilation" occurs when the freezein mechanism injects much more dark matter than is needed, but a large coupling between the DM and a dark mediator allows for the excess to annihilate down to the correct relic abundance, with this depletion typically completing near T ∼ m χ /10. Example evolution of the DM number abundance with temperature is shown for all three processes in Fig. 3.
In this minimal B −L model, dark matter freezeout is governed by the χχ → Z D Z D annihilation process with cross-section where r ≡ m Z D /m χ and α D ≡ g 2 D /4π. Freezein, however, is dominated by the direct production of DM from the SM through s-channel Z D exchange, ff → χχ: where η f = 1/3 (1) for quarks (leptons). A DM particle produced via freezein will, in the presence of a dark radiation bath atT , rapidly attain kinetic equilibrium in the parameter space of interest, though not necessarily chemical equilibrium. Thus, givenT as a function of T ,T (T ; ), we can obtain the relic abundance of DM by solving the single Boltzmann equation [24,28] dY is the equilibrium number density as dictated by the hidden sector temperature, relative to the SM entropy, which we approximate as conserved, 2 and Y χ,eq (T ) ≡ Y χ,eq (T, T ). LIDM is realized when the second term in Eq. 3.5 is unimportant for determining the final relic abundance, i.e., neglecting the effect of that source term will have only small effects on the final dark matter population. This can happen in two separate regimes. The first regime occurs when the hidden sector temperature is relatively close to the SM temperature, such that DM produced by freezein can reach thermal equilibrium with the dark radiation bath prior to freezeout (see Fig. 2): we call this "late" LIDM. The second regime occurs at very small values of and α D , where freezeout occurs before freezein stops, but freezein processes are sufficiently feeble to contribute only a tiny fraction to the final DM abundance. We call this more weakly coupled regime "early" LIDM. When the second term in Eq. 3.5 is not negligible, we find that generically the DM relic abundance is obtained through reannihilation. Freezein occurs when the first term is entirely negligible in comparison to the second term, and is realized in a very limited fraction of parameter space. Fig. 4 shows a schematic of the viable parameter space and of the mechanisms yielding the correct DM relic abundance in the minimal B − L model. At large portal couplings above the purple line, the HS and SM sectors are in thermal equilibrium, yielding a WIMP next door scenario [21]. At small portal couplings, the co-moving number density of dark matter is never high enough to produce the correct relic abundance (2.27). In practice, the high multiplicity of the SM sector and size of α s result in C E ∼ {few} × 2 T 5 , placing the absolute coupling floor near ∼ 10 −13 . For values of slightly above this floor, the hidden sector does not attain internal thermal equilibrium (for any m Z D ). While we will assume internal thermal equilibrium in this subsection, we will establish the validity The different regions in the parameter space. Leak-in dark matter is shown in light blue, with the narrow slice near the equilibration floor (purple line) corresponding to late LIDM, and the smaller epsilon region corresponding to early LIDM. Above the equilibration floor is the WIMP next door [21]. Reannihilation [28] is shown with shaded light red. The narrow slice of parameter space above m χ ∼ 3 TeV in shaded green indicates the three-solution region where three different choices of α D can produce a mostly LIDM, freezein, or reannihilation solution (see Fig. 3). At high masses, the solid red denotes where α D > 4π and non-perturbative couplings are clearly required to produce the right relic abundance. At very small , the model is below the absolute coupling floor, and not enough dark matter is produced to reach the relic abundance (dark green). The dark red shaded region above this indicates where the internal thermalization conditions are not satisfied. At small masses, constraints from the CMB, direct detection, and DM self-interactions together with constraints on the B − L boson forbid any valid DM solutions (shown in olive), as will be discussed in Sec. 4.
of this assumption in Appendix A. At high DM masses, the requisite α D becomes non-perturbative, while small DM masses are excluded by a combination of constraints on B − L vector bosons, CMB distortions, and DM self-interactions, as detailed below in Sec. 4. Late LIDM governs the DM abundance near the equilibration floor, while early LIDM governs the region of parameter space at small , and the two regimes transition smoothly into each other at smaller values of the DM mass. At larger masses and intermediate values of , re-annihilation governs the DM abundance. The boundary between reannihilation and late LIDM is set by requiring that the coupling as determined in a "leak-in-only" solution (i.e., σ ff →χχ v → 0 in Eq. 3.5) differs from the full solution by less than 10%. The boundary between reannihilation and early LIDM occurs  in practice when the value of α D , , and m χ cross the point where a "freezein-only" solution (i.e., σ χχ→Z D Z D v → 0 in Eq. 3.5) would produce the observed relic abundance. There is a narrow slice of parameter space at high DM mass and small , shown here in green, where a dominantly leak-in, a dominantly freezein, or a reannihilation solution can be achieved for different choices of α D . One such point is shown in Fig. 3. To distinguish between dominantly leak-in and dominantly freezein solutions, we consider the co-moving number density immediately after leak-in freezeout occurs, and ask whether it is greater (leak-in) or less (freezein) than 50% of the observed value. In all cases, near transitions both source terms in Eq. 3.5 are important for obtaining the final DM abundance. For a given m χ and , there is usually a unique value of α D that realizes the correct DM relic abundance, shown in the left panel of Fig. 5. Within the three-solution region, we display the α and σv values for the mostly leak-in solution. The corresponding annihilation cross-sections are displayed in the right panel of Fig. 5, where the wedge of the reannihilation region is clearly visible at high DM mass. In the absence of the freeze-in term in Eq. 3.5, the annihilation cross-section would display the simple scaling with expected from Eq. 2.22. However, as Fig. 5 shows, the presence of the freeze-in term instead leads to reannihilation and its larger annihilation cross-sections controlling the phenomenology. The net annihilation cross-sections are thus only slightly suppressed compared to expectations for a traditional WIMP over much of parameter space, with correspondingly better prospects for detectability; of course, DM in this model can also be much heavier than a traditional WIMP. At lower DM masses, where the leak-in solution dominates, the numerous mass thresholds of 10 -8 In Fig. 6, the possible solutions that provide the correct relic abundance in α D vs parameter space are shown for three different choices of m χ . At high mass, an appropriate choice of α D and could realize any one of the solutions. At moderate masses of m χ 3 TeV, there is no longer a dominantly freezein solution (unless there was never a formation of the dark vector plasma). For smaller masses m χ 10 GeV, there is no valid reannihilation solution because the dark matter abundance produced through freezein processes is either too small to account for the dark matter density or injected into a dark vector plasma where the dark matter is still in equilibrium. In this case, the early leak-in solution smoothly joins the late leak-in solution.
Generically, there is only one viable solution for α D for a given m χ and . However, at very high mass there is a region where different values of α D can provide a mostly late LIDM, mostly freezein, or reannihilation solution; these multiple solutions are manifest in Fig. 6 where the curve for m χ = 100 TeV becomes non-monotonic. In the sliver of parameter space where all three solutions are valid, α reann α F I α LI . This three-solution region is the only place in the parameter space that a mostly freezein solution can be found. However, in this construction we have implicitly assumed that the reheat temperature of the universe is large enough thatT RH , as dictated by the attractor solution, is larger than the dark matter freezeout temperature. If the reheat temperature was too low, or the hidden sector did not internally thermalize, freezein solutions can occur.

Signals of B − L vector portal LIDM
Despite the small size of the B − L portal coupling , there are many experimental handles on vector portal leak-in dark matter. In this section, we will discuss current limits on and potential future sensitivities to this parameter space.

Indirect detection
The same process that allows LIDM to freezeout can facilitate dark matter annihilation throughout the universe's history, including today. Indirectly detecting dark matter through these annihilation products is one of the most promising ways to probe LIDM models as the annihilation cross-section, which is s-wave in the B −L vector portal model, does not depend directly on the very small coupling to the SM particles .
Additionally, the exchange of light mediators can enhance the tree-level annihilation cross-section from Eq. 3.3 via the Sommerfeld effect [42][43][44][45]. The s-wave cross-section can be expressed as [6,21,46,47] where r = m Z D /m χ , v c is some characteristic dark matter velocity for the system of interest, and the Sommerfeld enhancement factor for a Hulthén potential (a good approximation to a Yukawa potential with nicer analytic properties [48]) is [21] S(α D , r, v) = 2πα D v sinh 6v πr cosh 6v πr − cosh 2) to place constraints on the parameter space. We do not consider the influence of the Sommerfeld effect on freezeout, as this would primarily affect only the large α D region which corresponds not to leak-in, but reannihilation, discussed in Sec. 3.2, that produces the bulk of the relic abundance. Some of the most stringent constraints on annihilating dark matter come from the detailed measurements of the CMB power spectrum [49,50]. Injection of energetic charged particles and photons into the plasma can distort the CMB anisotropies. Planck, SPT, ACT, and WMAP results restrict the power injected into the CMB from DM annihilation, per DM mass, to satisfy f ef f (m χ ) σv /m χ < 14 pb c / TeV [50], which allows for robust bounds to be placed on dark matter models. The effective energy deposition efficiency f ef f (m χ ) [50,51] depends on the branching fractions into specific annihilation channels, but it is 0.4 − 0.6 for electron-and photon-enriched annihilations, small for neutrinos, and typically ∼ 0.2 for everything else in the SM. Despite the smallness of the energy deposition efficiency for neutrinos, at high DM masses neutrino-induced energy deposition into the CMB can be large enough that even dark vectors that are only able to decay to neutrinos are excluded.
Despite the current excess above predictions [14], the observed positron flux at the AMS-02 experiment can be used to constrain dark matter annihilations that result in positrons. We follow [52] in choosing to bound σv × Br(Z D → e + e − ), which is most stringent for dark vectors in the range 2m e < m Z D < 2m µ .
The Fermi Large Area Telescope experiment (Fermi-LAT) [13] has observed gamma ray spectra for many dwarf galaxies [53], including many ultra-faint dwarf galaxies observed by the Dark Energy Survey (DES) [54]. As several dwarfs have low noise and large astrophysical J-factors, these observations can severely constrain dark matter annihilations [55]. We use log-likelihood-ratios (LLR) provided by the Fermi collaboration for 24 energy bins for the 41 dwarf galaxies within the nominal sample of Ref. [53]. We approximate the effect of correlated systematics by applying a 0.5σ downward shift to the J-factors for each dwarf galaxy, as this was determined to closely replicate the limits placed by Fermi [21]. The 41 dwarf LLRs are combined within each energy bin. We generate the gamma ray spectra from dark matter annihilation in Pythia 8 [56] at each point in m χ vs m Z D . All 24 bins are combined to form the χ 2 with one degree of freedom to place limits on the annihilation rate σv .
Additionally, the observations of the galactic center performed with the High Energy Stereoscopic System (H.E.S.S.) experiment [15] can place tighter constraints on heavier dark matter (m χ 1 TeV) than Fermi due to the very large J-factors expected at the galactic center. 3 As H.E.S.S. does not provide their data, we follow Refs. [66,67] and use the 112 hour data from a gamma ray study [68] to obtain an observed gamma ray spectra for the signal and background regions and simply scale these results up to 254 hours of data to project fairly conservatively what an updated study could achieve. We again use Pythia 8 [56] to generate the annihilation signal gamma ray spectra at each point in m χ vs m Z D . The effective area was collected from [69]. With the statistical procedure outlined in [70], we derive a χ 2 with one degree of freedom to approximate the limits on the annihilation rate σv that H.E.S.S. would be able to find. While this procedure allows us to place approximate limits on the model, firmer statements would be possible if H.E.S.S. were to provide the tools required to reliably recast their results, e.g. by providing the LLRs for a signal + background hypothesis for each energy bin as a function of injected signal strength.
The interplay of these constraints in the parameter space is illustrated in Fig. 7. For four fixed values of (at 10 −9 , 10 −10 , 10 −11 , and 10 −12 ), we show CMB constraints (blue), positrons from AMS-02 (pink), gamma-rays from dwarf galaxies at Fermi-LAT (red), and gamma rays from the galactic center at H.E.S.S. (green). In the lower left corner of each figure, the dark vector is below 2m e and decays to neutrinos. While most indirect detection constraints considered here are insensitive to neutrinos, the CMB power spectrum can be sufficiently distorted by very energetic neutrinos that arise from heavy dark matter annihilations [50]. The sharp transition near 100 TeV in the lower right panel occurs where the model moves from leak-in dark matter to reannihilation and the cross-section jumps due to the larger α D needed for the correct relic abundance.
Finally, Fermi-LAT observations of the smooth galactic halo may place more stringent constraints than dwarfs [71]. The Cherenkov Telescope Array (CTA) [17] would greatly enhance the sensitivity to heavy dark matter [72]. A full treatment of these (potential) limits is beyond the scope of this work.

Direct detection
Despite the smallness of the B −L portal coupling , direct detection experiments can be an important probe of leak-in dark matter. Dark vector exchange contributes to the non-relativistic, spin-averaged amplitude-squared for DM-nucleus scattering, which can be written as where E R is the nuclear recoil energy, m N and A are the mass and mass number, respectively, for the target nucleus (Xenon, in the case of interest), F 2 (E R ) is the nuclear form factor, for which we take the Helm form factor [73,74]. When m Z D v χ m χ m N /(m χ + m N ), the recoil energy dependence in the propagator is necessary to properly track the transition into a long-range interaction. When the amplitude is independent of the DM velocity, the event detection rate of the experiment per unit detector mass can be written as [75,76] where ρ χ = 0.3 GeV/cm 3 is the local DM density, (E R ) is the selection efficiency specific to the experiment, and the mean inverse speed η is defined by [76] η for which we use the expression in Ref. [74] to match the experiments (and not the more accurate expression found in Ref. [76]). If it is reasonable to approximateM N R (E R ) →M N R (0) in Eq. 4.4 (as is typical for contact interactions), then the particle physics inputs may be factorized from the astrophysical and experimental inputs. Most direct detection results are presented using a crosssection that has been both factorized in this manner and posed in terms of an effective cross-section per nucleon. Defining the reduced mass of the nucleon-DM system as µ χn = m χ m n /(m χ + m n ), the per-nucleon-DM cross-section in this model is (4.6) As can be seen from Eq. 4.3, the assumption of recoil energy independence breaks down when m 2 χN v 2 χ . We will determine the excluded cross-section via 10 -2 In red, we show where the model is required to be at or below the late leak-in dark matter scenario. In green, the model is required to be in reannihilation or below, while the narrow purple sliver is in early leak-in. See Fig. 4 for more details.
in order to correctly account for this important effect at low mediator masses. The latest XENON1T limits [18] place the tightest constraints in the parameter space. We show the current limits and regions where direct detection forces the model below the equilibration floor in Fig. 8. Recent limits from DarkSide-50, CRESST-III, and EDELWEISS [77][78][79] probe lighter masses, but are currently not sensitive enough to place meaningful constraints below the equilibration floor. The sensitivity scales as 1/m 4 Z D , but saturates when m 2 Interestingly for dark matter above 100 GeV, the limits from direct detection are nearly independent of the dark matter mass, and set the same constraint across the m Z D vs plane. This is because the dark matter flux drops as 1/m χ , while σ 0 χn ∝ α D , which for reannihilation also scales roughly as m χ . Several recent experiments, notably SENSEI [80], DAMIC [81], XENON10 [82], SuperCDMS [83], and DarkSide-50 [84], have constrained very light dark matter scattering off of electrons. The relevant cross-section for this is simply [85]  (4.9) As we will uncover in the next two sections, ∼ 10 −10 is the largest value possible and ∼100 keV is lightest that a B − L vector can be for dark matter in this mass range. With current constraints in thē σ χe ∼ 10 fb range, there are several orders of magnitude further to probe before these electron recoil experiments could have sensitivity to this model, sensitivity which may be achievable by the proposed DAMIC-M [85].

Constraints on the B − L vector boson
Independent of the nature of the dark matter, there is a wide variety of experimental searches for a U (1) B−L gauge boson. In the region of dark vector masses and coupling strengths of interest for LIDM, 10 −7 g B−L ≡ 10 −14 , the most important constraints come from fifth force experiments [19,[86][87][88][89], N ef f constraints on the number of relativistic species present during BBN due to the dark B − L vector maintaining thermal equilibrium with the neutrinos after decoupling from the electronphoton plasma, resulting in heating of the neutrino sector [90,91], the cooling of Supernova 1987A [91,92], the electron beam dump E137 [93,94], the neutrino experiment LSND, interpreted as a proton beam dump [95,96], and especially stellar cooling through emission of dark vectors in the sun (Sun), horizontal branch stars (HB), and red giants (RG) [19,91]. The current limits on weakly coupled B − L vector bosons are summarized in Fig. 9.
It is possible there are additional constraints both from SN1987A, where a B − L dark vector decays to positrons that at late times contribute to 511 keV gamma ray signal [97], and BBN, where B − L dark vectors that are produced in the early universe, survive until BBN, then decay causing photo-disintegration of nuclei [98]. Derivation of the specific constraints for a B − L dark vector is beyond the scope of this work.

Dark matter self-interactions
B−L vector portal LIDM can have sizable self-interactions, especially in the regime r = m Z D /m χ 1. The most stringent limits on DM self-interactions in this model come from the measured ellipticity of galaxy halos, or (less precisely) through the generation of very large cross-sections on dwarf scales, far in excess of those cross-sections that yield acceptable dwarf galaxy properties in simulations.
α D r, there is an analytic closed-form expression for the transfer cross-section valid for both χχ → χχ and χχ → χχ (calculated with Package-X [99]), In the classical Rutherford regime, the scattering is long-range, so that the momentum transfer is large compared to the mediator mass v r, and non-perturbative since α D r. In terms of β ≡ 2α D rv 2 , the transfer cross-section can be approximated as [100]   for an attractive interaction and  for a repulsive interaction. Since the symmetric DM in this model is composed of an equal number of particles and antiparticles, we will take σ T = 1 2 (σ + T + σ − T ). Between the Born regime and the classical regime is the resonant regime, characterized by α D r and v ∼ r, where the transfer cross-section has a complicated velocity dependence. The transfer cross-section here can be calculated by summing up contributions from a sufficiently large number of partial waves [47]. This procedure is computationally expensive, however, and in much of the parameter space we will be able to bypass computations in the resonant regime by employing a bounding method, described below.
Self-interaction cross-section on dwarf scales. If dark matter has too large of a transfer crosssection, then the galactic properties produced in N -body simulations do not match observations. Comparison of N -body simulations with observations suggest upper bounds on σ T /m χ in dwarf systems of order 10 cm 2 /g ≈ 20 barn/GeV or below, both in constant cross-section models [101,102] and long-range models [103].
Meanwhile simulations (albeit of constant self-interaction cross-sections) indicate that crosssections in excess of 50 cm 2 /g ≈ 100 barn/GeV begin to exhibit core collapse in dwarf galaxies [104,105].
To evaluate σ T , we construct two separate transfer cross-sections thermally averaged over a Maxwellian velocity distribution with v RM S = 30 km/s, one assuming the Born cross-section σ B T , and the other assuming the classical Rutherford cross-section σ R T . For a given value of m χ and (and thus α D from the relic abundance condition), we solve for the value of m Z D that realizes σ X T = 50 cm 2 /g × m χ for both cases (X = B, R). After solving for this minimum allowed m Z D in both regimes, we check that the solution is self-consistent (specifically, we require α D /r < 10 −1 in the Born regime and in the Rutherford regime both α D /r > 1 and v/r > 100). This procedure then produces a curve in the (m Z D , ) plane indicating (for a given m χ ) where the specified thermally averaged transfer cross-section is obtained. In Figs. 10 and 11 we show curves for both 10 cm 2 /g and 50 cm 2 /g, indicating where transfer cross-sections begin to exceed the values where N -body simulations accord with observations. Parameter points with larger self-interaction cross-sections on dwarf scales are disfavored.
However, when the thermal averaging is performed in the classical regime, this procedure is not completely accurate, as the classical expression for σ T is only valid when v r. For sufficiently small relative velocities, scattering occurs in the resonant regime instead. In order to overcome this issue, we consider the following bounding method which will also let us largely bypass the necessity of calculating the transfer cross-section in the resonant regime. Consider the thermally averaged cross-section in the non-perturbative regime: where we split the integral into classical and resonant contributions. Since both terms are nonnegative, we obtain the following lower bound: In order to construct an upper bound, we note from [100] that the classical cross-section is an overestimate in the resonant regime. Hence, σ R T ≥ σ resonant T , so that As long as the scattering process is indeed non-perturbative, this method gives us a bounding region for the constraint curve. The upper and lower bounds constructed in this manner often either closely coincide, or both lie deeply within excluded regions. Only for DM masses around m χ = 10-100 GeV do we need to explicitly evaluate the resonant contribution to the thermally-averaged transfer crosssection. In Figs. 10 and 11 we indicate with the blue hatched region the bound from resonant and/or classical (Rutherford) scattering, and with the red hatched region the bound from Born scattering.
Ellipticity. DM self-interactions will tend to increase isotropy within galaxy haloes. In particular, the measured ellipticity of the gravitational potential of the galaxy NGC720 [106] places a bound on DM self-interactions [107]. We here use a simple treatment of the ellipticity bound based on estimating the timescale τ e for isotropizing the velocity dispersion in a halo and requiring that it exceed the age of the universe [107], where the average DM energy E is given in terms of the velocity dispersion v 2 0 by (we take the velocity distribution to be locally given by a Maxwell-Boltzmann distribution). Meanwhile the average energy transferred in a DM-DM collision is given by For simplicity we evaluate Eq. 4.18 with ρ = 2.1 GeV/cm 3 and v 0 = 260 km/s, corresponding to the middle of the range of values reported in Ref. [107].
For this model, the integral in Eq. 4.20 is regulated by the finite dark vector mass. However, when the dark vector mass is sufficiently small compared to the momentum transfer, the integral will first be cut off by the net charge neutrality of the dark plasma, i.e., by requiring that the maximum impact parameter be smaller than the inter-particle spacing λ pp = (m χ /ρ χ ) 1/3 [108]. Thus τ e becomes independent of m Z D when m 2 where y = mχv 2 0 αχ λ pp . We use Eq. 4.18 as the constraint, which underestimates the time required to attain an isotropic distribution as it does not take into account the reduction in the energy transfer rate as initially anisotropic populations approach equilibration [108]. The resulting constraints are shown in brown in Figs. 10 and 11. This bound is conservative for the purposes of identifying clearly allowed regions, but (as argued in Ref. [108]) there are several ambiguities in translating the measured ellipticity of galaxy haloes into bounds on DM self-interactions, making it hard to conclude that the shaded regions to the left of this bound are definitively excluded.

Allowed parameter space
The regions of dark vector parameter space consistent with LIDM are shown for several different values of DM mass in Figs. 10 and 11. For a fixed DM mass m χ , there is a specific region in the (m Z D , ) plane consistent with the LIDM mechanism. For sufficiently large , the SM and the HS attain thermal equilibrium before DM freezeout, while for sufficiently small , DM will never obtain a a large enough co-moving number density to account for the relic abundance observed today. Internal thermalization (see Appendix A) provides a more stringent, but less robust, condition than underabundance; we will show lower bounds from both thermalization and absolute abundance on the plots below. Meanwhile, the upper bound on m Z D simply reflects the requirement that r = m Z D /m χ ≤ 0.1, so that (in the minimal model) the dark vector constitutes a relativistic radiation bath at DM freezeout.
From Figs. 10 and 11, we identify two distinct regions of parameter space consistent with dark vector constraints. First is an "invisible" region where the dark vector mass lies in the narrow window between stellar cooling bounds and CMB constraints on DM annihilations, 100 keV m Z D < 2m e . In this regime, the dark vector decays entirely to neutrinos, rendering DM annihilation (largely) invisible to cosmic ray searches. The second, "visible", region of parameter space occurs where m Z D > 2m e , and DM annihilation produces visible cosmic ray signals. Stringent constraints on very light dark B − L gauge bosons, combined with the excessively large DM self-interactions generated when m Z D ≪ m χ , disfavor values of m Z D below tens of keV.
For m χ 100 MeV (Fig. 10), the combination of dark vector bounds, the restriction m Z D ≤ m χ /10, the requirement of internal thermalization (light red), and CMB constraints on DM annihilation (green) leave only the small invisible region available. The narrow window of surviving  (Fig. 9). The horizontal purple dotted line shows the equilibration floor, the horizontal green dash-dotted line shows the absolute floor, and the vertical dashed black line indicates m Z D = 0.1m χ . The light red region indicates regions that do not attain full internal thermal equilibrium. Direct detection exclusions are shown in yellow (Fig. 8), and CMB constraints in green (Fig. 7). The red hatched region shows regions with σ T /m χ > 50 cm 2 /g (left curve) and 10 cm 2 /g (right curve), self-consistently computed in the Born regime (Sec space, where indirect detection searches have no reach. In the top left plot, for m χ = 10 GeV, we perform a full resonant calculation for both σ T /m χ and ellipticity. In the top right plot, for m χ = 100 GeV, the left (right) solid blue line and left (right) dashed blue line form the brackets for the classical regime calculation saturating σ T /m χ at 50 cm 2 /g (10 cm 2 /g). The light blue shaded region highlights the bracketing region for 50 cm 2 /g. Similarly, the solid and dashed brown lines bracket the ellipticity constraint, with the light brown shaded region highlighting the bracketed region. As the DM mass increases, Sommerfeld-enhanced indirect detection signals become increasingly effective at probing the parameter space. This is unsurprising, as the bulk of the high-mass parameter space is in the reannihilation regime, where the relatively large values of α D and σv accordingly yield interesting indirect detection signals. The remaining unexcluded territory is predominantly in the more weakly coupled leak-in regime, where indirect detection signals are much fainter.
For m χ 100 GeV, we obtain viable parameter space realizing LIDM with DM self-interactions σ T /m χ ∼ few cm 2 /g in dwarf systems, i.e., in the range of interest for addressing small-scale puzzles in galaxy formation, that are not obviously in tension with ellipticity constraints.

Summary and conclusions
In this paper we have examined in detail the properties of leak-in dark matter: dark matter that freezes out of a hidden sector evolving in a non-adiabatic leak-in phase. The quasi-static equilibrium leak-in phase, in which the energy density of the hidden sector redshifts like matter, is a generic behavior that emerges when a cold hidden sector is dominantly populated through a dimension-four interaction with the hotter SM. We provide analytic methods for consistently treating the out-of-equilibrium evolution of the hidden sector temperature in the presence of a known collision term.
We present a detailed study of DM freezing out of a leak-in radiation bath and the resulting observational consequences. The renormalizable nature of the interaction feeding the hidden sector radiation bath ensures that the cosmological evolution of the hidden sector is minimally sensitive to details of the unknown physics of reheating in our universe. This class of DM models are thus sharply predictive, and have a bounded parameter space. The strength of the interaction cannot be too large, in which case the interaction will reach equilibrium, or too small, in which case the dark sector will never reach a high enough internal temperature to produce the observed DM relic abundance. Meanwhile, the DM mass is bounded from above by the requirement of perturbativity, and from below by a (model-dependent) combination of terrestrial, astrophysical, and cosmological constraints. In an outof-equilibrium hidden sector, the DM relic abundance is determined by an interplay of freezeout and freezein processes, resulting in a rich solution space.
To establish some concrete constraints on and predictions from LIDM, we specialize to a particular model, where the dark sector consists of fermionic DM together with a dark vector boson that couples to the SM via the B − L current. Despite the smallness of the portal coupling , there are many experimental probes of this B − L LIDM model. While the DM annihilation cross-section is suppressed compared to standard WIMP scenarios thanks to the relative coldness of the hidden sector, indirect detection signals do not depend directly on the small portal coupling , and provide excellent sensitivity to large regions of the parameter space. In particular, this model can realize very large DM masses (m χ ∼ 10s -100s of TeV) with striking cosmic ray signals of DM annihilation, detectable due to sizable Sommerfeld enhancements in the late universe from the relatively large dark coupling constant. Additionally, the enhanced cross-sections obtained from light mediator exchange enable direct detection experiments to probe the cosmic history, and not just the particle content, of thermal dark sectors. In fact, XENON1T now provides the leading constraints on the very weakly coupled LIDM regime when m Z D < 2m e and indirect detection signals are suppressed.
Portions of the LIDM parameter space can realize very large DM self-interaction cross-sections. The combination of (i) stringent constraints on low-mass B −L gauge bosons, (ii) enormous DM self-interaction cross-sections, and (iii) the requirement of internal thermalization eliminates all parameter space where the B −L boson lies below the constraints from stellar cooling, ∼ 100 keV. Astrophysical tests of DM self-interactions could potentially provide a unique observational handle on the low-mass regions of LIDM parameter space, where neither direct nor indirect detection are sensitive. Viable parameter space at high masses, m χ ∼ 10 -100 GeV, can have DM self-interaction cross-sections that fall in the astrophysically interesting range σ T /m χ ∼ few cm 2 /g compatible with small-scale structure anomalies in dwarf systems.
Leak-in dark matter represents a simple, generic, and sharply predictive class of models for the origin of dark matter in our universe. For that reason, exploring the signature space of both this B − L model and other realizations of LIDM, coming from other choices of leading interactions between the SM and the dark sector, is an important aspect of broadening the search for DM.
Note added: While this work was nearing completion, the works Refs. [109][110][111] appeared, containing related but not identical material.

A Attaining internal thermalization
In order for the dynamics described here to be an accurate description of the hidden sector, the dark radiation bath must have sufficiently rapid self-interactions to attain internal thermal equilibrium. This criterion depends on the properties of the dark radiation bath itself, and is therefore necessarily somewhat model-dependent. In this subsection we will present an approximate criterion for internal thermalization of the minimal B − L vector portal hidden sector.
For the hidden sector to attain internal thermal equilibrium, processes that change the numbers of individual dark species must be efficient on cosmological timescales. At leading order, such a process is provided by the elastic scattering Z D Z D → χχ. Given a number density n Z D of "hard", pre-thermalized dark vectors, it is straightforward to estimate the rate Γ el for this process. The number density of dark photons in the absence of subsequent scattering within the hidden sector can be obtained by solving the Boltzmann equatioṅ under the simplifying assumptions that H ∝ T 2 depends only on the SM temperature, backward contributions to the collision term can be neglected, and the SM temperature simply redshifts as T ∝ 1/a. The collision term can be estimated as where σ qg→qZ D v ∼ 2 α s /(24T 2 ) is the spin-and color-averaged cross-section. Solving Eq. A.1 yields Comparing this result for n Z D to the analogous estimate for the energy density injected into the HS (see Sec. 2), we can see that (as expected) the typical energy carried by one of these hard dark vectors is ∼ T . The corresponding rate for initial production of DM particles from the primordial dark vector population is then It is worth observing that the essential parametrics of this elastic rate hold for any elastic 2 → 2 process occurring among the initial hard population of particles in the dark sector. Given Γ el , we can quickly estimate whether elastic scattering is sufficiently rapid to thermalize the hidden sector by requiring that Γ el > H at some temperature T > m χ . This estimate indicates that elastic scattering suffices to thermalize much but not all of the leak-in parameter space. However, inelastic scattering, χX → χX + Z D , is more effective than elastic scattering at thermalizing the hidden sector over much of the parameter space of interest. The importance of inelastic scattering in thermalizing a sector containing gauge interactions is well-known [112][113][114][115]. While the inelastic scattering process is higher-order in α D , it can be sufficiently enhanced by the region of low momentum transfer to more than compensate for the additional α D suppression. Our estimate of thermalization through this inelastic process will be parametric, and largely follows the related treatment in [115].
The inelastic scattering rate is approximately given by where µ is the effective IR scale that regulates the t-channel Z D propagator, n χ indicates the number density of hard DM particles produced directly from the SM, and we have temporarily neglected the possible complications that arise when the timescale for emitting a soft vector boson in the final state becomes longer than the timescale between hard scatterings, i.e., the Landau-Pomeranchuk-Migdal (LPM) effect [116,117]. The number density of hard χ particles, in the absence of subsequent scattering within the hidden sector, can be obtained analogously to the estimate for n Z D above. We can estimate σv ≈ ( f g f (Q f B−L ) 2 πα D 2 /T 2 for ff → χχ, giving the collision term Then, solving the Boltzmann equation for n χ yields 10 -4 10 -3 10 -2 10 -1 1 10 100 10 3 10 4 10 5 10 6 10 7 10 -13 10 -12 There are three possibilities for the effective IR scale µ that cuts off the momentum transfer in Eq. A.5. First is simply the (vacuum) dark vector mass itself, m Z D . Second is H, reflecting that the horizon is the largest range of physical interest for the dark interaction. Finally, in the medium, the dark vector propagator receives corrections from its interactions with the plasma. The screening scale in the non-equilibrium dark plasma can be estimated as [118] in terms of the hard DM population n χ . For the B−L dark vector, we should in principle also consider the contribution to its effective mass from interactions with the SM plasma, m SM,T ∼ T . Over our parameter range of interest, we find that both Hubble and the SM contribution to the dark vector's effective mass are always negligible in comparison with µ sc and m Z D . These possible screening scales have varying dependence on T , , and α D ; at any given temperature, the largest is the one that is physically relevant. Now, when the timescale for emitting a soft dark vector is larger than the typical timescale between 2 → 2 collisions, the inelastic 2 → 3 scattering can no longer be discussed in isolation. The result of multiple 2 → 2 scatterings occurring during the so-called "formation time" governing the 1 → 2 splitting is known as the LPM effect, and can be formally understood in an effective Boltzmann treatment by defining an effective splitting function that resums specific contributions to the amplitude from successive scatterings [118,119]. Destructive interference among these contributions results in a suppression of the brehmsstrahlung rate. Thus we need to correct the estimate of the inelastic rate for 2 → 3 scattering in Eq. A.5 with a factor f LP M ≤ 1 to account for this suppression, We use the estimate of [115] (see also [114]) for f LP M in the Abelian plasma: , the LPM suppression is not operative, so f LP M = 1. When the LPM effect is operable, i.e., f LP M < 1, the net inelastic rate is given by Γ inel = α 2 D n χ /T . If m Z D is small, then the LPM effect is active everywhere in the parameter space of interest. It is worth noting that this estimate for f LP M assumes an adiabatic evolution of n χ in estimating the evolution of the formation timescale. This is an underestimate of the non-adiabatic population of hard n χ , and therefore an underestimate of Γ inel . While this treatment could in principle be improved, it is a conservative choice, and further refinement is beyond the scope of this paper.
In Fig. 12, we show where the internal thermalization conditions are not satisfied, i.e. when Γ el + Γ inel < H at the freezeout temperature. In practice, 2 → 2 processes are more important at lighter DM masses, while for higher DM masses the 2 → 3 process are more important. In Figs B Alternative UV models for hierarchical B − L charges Throughout this paper, we introduced an extremely large B − L charge for DM to create a disparity in the B−L vector boson's coupling to SM fields compared to dark matter. This model has the advantage of having clear predictions and no UV sensitivity; however, the large DM charge invites modelbuilding questions. In this appendix, we present two simple models that provide an explanation for the hierarchical couplings of the dark U (1) gauge boson to DM and the SM B − L current, and briefly sketch the impact of the added states on the DM signatures. Both models involve a U (1) B−L ×U (1) D symmetry, with the first introducing kinetic mixing and the second introducing a Higgs state. In both models, if the reheat temperature is too high, there is a danger that the B − L vector could thermalize the SM and hidden sectors. where J µ B−L is the SM B − L current, and χ is the dark matter. In other words, we start with a model where, in the gauge basis, the dark gauge boson talks only to dark matter, and will inherit its couplings to the SM B − L current through kinetic mixing with a new B − L gauge boson. We assume that this B − L gauge boson gets a large mass through spontaneous symmetry breaking, m X,0 (the origin of this mass term, Higgs or Stückelberg, is unimportant).
Making the customary field redefinition and redefining g D =ĝ D / √ 1 −ˆ 2 , yields diagonal kinetic terms for the gauge bosons, and couplings to matter of the form Here we have defined Given masses m 2 Z D ,0 and m 2 X,0 forẐ D andX, the resulting mass-squared matrix forZ D andX is where δ 2 ≡ m 2 Z D ,0 /m 2 X,0 . We will be interested in δ 2 1. This mixing matrix is diagonalized by Expressing δ 2 in terms of the eigenmass m Z D , we have where in the last step we expanded to leading order inˆ (assuming m 2 D m 2 X ). Thus the mixing angle can be written as giving the two eigenstate couplings to matter (to leading order inˆ ) For δ,ˆ 1, the effective Z D coupling to the SM B − L current is then the product of the underlying g B−L and two independent small parameters, which is the small portal coupling used throughout this work. Importantly in this model, the heavy B − L vector couples dark matter to the SM particles at the same order as the lighter dark vector, which results in a cancellation of the leading amplitude for direct detection processes: which is suppressed by 2m N E R /m 2 Z D relative to Eq. 4.3.

B.2 Dark mixed Higgs
As before, this model has three U (1) factors in the UV: a U (1) D gauge boson, a separate U (1) B−L , and SM hypercharge. Additionally, we introduce a scalar field φ that has charges {Q D,φ , Q B−L,φ } under the U (1) D and U (1) B−L symmetries. The terms in our Lagrangian important for this discussion are where kinetic mixing is assumed to be absent. The mass for the vectorX µ could arise from a Stückelberg or Higgs mechanism, but this origin is unimportant. The gauge bosons couple to matter, notably φ, through covariant derivatives of the form In standard fashion, V (φ) results in a VEV for φ, φ = w, so that our low-energy mass matrix has the form M 2 V = m 2 X,0 1 + κ 2 δκ δκ δ 2 (B. 16) where κ = g B−L Q B−L,φ w/m X,0 and δ = g D Q D,φ w/m X,0 , and κ, δ 1. Diagonalizing this matrix gives masses that are simply m 2 X ≈ m 2 X,0 (1 + κ 2 ) and m 2 Z D = g 2 D Q 2 D w 2 + O w 4 /m 2 X,0 , and a mixing angle, sin θ ∼ δκ. After this the two eigenstates X µ and Z µ D couple to matter as Unlike the previous model, the heavy B − L vector contributions to dark matter -SM interactions are unimportant in the IR. In principle, the remaining scalar degree of freedom from φ could affect the model in a few ways. It could be in the plasma, which could affect bothg * and rates relevant for internal thermalization. One way to reduce phenomenological consequences from φ would be to introduce a fairly small Q D,φ , which can allow for a very large separation between m φ and m Z D . With m φ m χ , φ is effectively removed from the low-energy theory.

C Collision term
This Appendix collects details concerning the calculation of the energy transfer collision term C E governing the temperature evolution of the hidden sector.

C.1 Away from the equilibration floor
The hidden sector temperatureT can be numerically determined as a function of the SM temperature T by the following procedure. For an internally thermalized hidden sector, the energy density stored there defines its temperatureρ which holds provided that (i) the SM dominates the entropy in the universe and (ii) g * S is slowly varying, so that T 3 a 3 = const holds to a good approximation (near the QCD phase transition, this assumption will not be good). A particularly useful variable is ξ =T /T , the ratio of hidden sector to SM temperatures. Noting that dT /dT = T dξ/dT + ξ, we can express Eq. 2.2 as dξ dT = 30C E (T,T ) 4π 2 H(T )g * (ξT )ξ 3 T 5 . (C.4) Assuming again thatg * is constant in the region of interest and the hidden sector is sufficiently cold so that the transfer of energy out of the hidden sector is negligible, C E (T,T ) ∼ C f E (T ), we can solve Eq. C.4 to obtain Here we have used that ξ(T i ) ξ(T ) which is always true in this model of interest for a sufficiently high value of T i .

C.2 Near the equilibration floor
Near the equilibration floor, the collision term in Eq. C.4 can be expanded in terms of a parameter δ = 1 − ξ that goes to 0 when the two sectors are equilibrated, where we have used that C b E,0 (T ) = C f E (T ). Given the functions C b E,n (T ), the resulting equation, can be straightforwardly numerically integrated near the equilibration floor where the backward collision term becomes important. To derive the functions C b E,n (T ), we will (as throughout) use Maxwell-Boltzmann statistics. The backward scattering piece of the collision term (2.10) can be written as The collision term can be related to the cross-section σ(s) for the given process using is the dimensionless two-body kinematic factor. We can thus write (C.8) as where we have defined the dimensionless quantity R ≡ m 2 4 − m 2 3 /s. It is possible to integrate over E − in Eq. C.11 analytically to yield At this point, we define δ = 1 −T /T and expand the above expression in powers of δ. The terms in the resulting series can individually be integrated over E + analytically. As κ ≡ 1 − R appears frequently in these expressions, it is convenient to replace R with κ below. For any given scattering process, the backward collision term can thus be approximated using where the functions A