Dark Matter Abundance from Sequential Freeze-in Mechanism

We present a thorough analysis of the sequential freeze-in mechanism for dark matter production in the early universe. In this mechanism the dark matter relic density results from pair annihilation of mediator particles which are themselves produced by thermal collisions of standard model particles. Below some critical value of the mediator coupling to standard model fields, this sequential channel dominates over the usual freeze-in where dark matter is directly produced from thermal collisions, even when the mediator is not in thermal equilibrium. The latter case requires computing the full non-thermal distribution of the mediators, for which finite temperature corrections are particularly important.


I. INTRODUCTION
The nature of the dark matter (DM) is perhaps the most acute open question in particle physics. Despite the strong observational evidence for an important DM component in the Universe, most of its properties remains unknown. Requiring that the DM be a thermal relic implies non gravitational interactions with ordinary matter. This nourishes hope to probe the DM in our local environment, either by detecting DM particles, directly in underground experiments or indirectly through the products of DM annihilation within our galactic neighborhood, or by producing them at colliders.
In particular, the hypothesis of a new DM particle around the electroweak scale has been tested extensively and the lack of evidence for DM from these searches triggered a growing interest in exploring a wider class of DM models. One possibility is that the DM and/or the mediator connecting it to the SM is below the GeV scale, thus leaving no traces in nuclear-scattering based direct detection experiments and colliders. This regime can be tested nonetheless with different experimental probes [1,2]. Besides, a mediator in the sub-GeV mass range also helps in resolving small-scale problems related to DM [3][4][5].
Another possibility is that the DM particle still lies above the GeV scale but signals in standard searches are suppressed because it interacts only very weakly with the SM. In this scenario, the DM is never in equilibrium with the SM in the early Universe and is instead produced by freeze-in through pair annihilation or decay of particles in the thermal bath [6,7]. The DM is generally assumed to be singlet under the SM group and part of a hidden sector that couples to the SM through renormalizable portal operators or with a mediator particle. Despite its tiny couplings with the SM fields the DM relic density often remains detectable in existing experiments [8], like direct detection when the mediator is light [9,10] or indirect detection for decaying DM [11,12]. Frozen-in DM could also be tested in cosmology [13][14][15] and at colliders through signatures involving displaced vertices and/or long-lived particles [16][17][18][19][20][21][22] including with detectors located some distance from the interaction point [23][24][25][26].
In this work we consider scenarios where dark sector particles are feebly coupled to the SM and where the DM relic is produced non-thermally through the freezein mechanism. We assume the mediator mass is small, at the 10 MeV scale, and is a scalar, for simplicity. 1 Lighter mediators also have interesting phenomenology but suffer severe constraints from Big Bang nucleosynthesis. The light mediator can potentially provide an explanation for anomalies observed at the cluster scale [29,30] and enhances the DM-nucleus scattering cross section, which offers the possibility to probe this scenario in direct detection. We will consider a simplified model with an hadrophilic scalar mediator which couples only to light quarks, thus alleviating several constraints that affect Higgs portal models where couplings to heavy quarks dominate [31][32][33]. This framework will be sufficient to illustrate the main phenomenological properties that are expected to be relevant for a larger set of models.
For freeze-in to take place, the product of mediator couplings to the SM and to DM must be very small, O(10 −12 − 10 −10 ), while their relative size remains a free parameter [7]. In models where DM is much heavier than the mediator, DM can generically be produced via pair annihilation of SM particles or mediators, assuming the latter are in equilibrium with the SM thermal bath. Here, we point out that even when the mediator coupling to the SM is too small for the mediator to ever reach equilibrium, a finite density of mediators can be produced through SM induced processes. The contribution of such non-thermal mediators to DM production parametrically dominates over that of pair annihilation of SM particles. This new phase of DM production, called sequential freeze-in [34], requires first solving for the momentum distribution of the mediator yield before using it for DM production. For this we solved the unintegrated Boltzmann equation rather than assuming φ to be in kinetic equilibrium with the thermal bath. Because the mediator is much lighter than DM, the tail of the mediator distribution is most relevant for DM production. Moreover we show that thermal effects, which we approximate by taking into account thermal masses for quarks, gluons and photons 2 play an important role in mediator production. Finally, our calculation takes into account the Bose-Einstein and Fermi-Dirac distributions for bosons and fermions, respectively, rather than simply using Maxwell-Boltzmann (MB) distributions. We show that these effects lead to large corrections in the calculation of the relic density [36].
The paper is organised as follows. We first present the simplified model and discuss at length the possible DM production mechanisms, giving an extensive description of the specific case where the mediator is not in thermal equilibrium (Section 3). With the complete calculation of DM production in hand (Section 4), we then determine the potential of current and future DD detectors to probe this model (Section 5) and examine numerous astrophysics and terrestrial constraints (Section 6). Our conclusions are presented in Section 7. The appendices contain details on the reaction rates for mediator production as well as an approximate analytical solution for the mediator distribution.

II. SIMPLIFIED DARK MATTER MODEL
Consider a simplified model for the dark sector which comprises a Dirac fermion χ (the DM candidate) interacting with SM quarks through a real scalar mediator φ with parity-preserving renormalizable couplings where the Lagrangian above is written below the scale of EW symmetry breaking. Both dark states are neutral under the SM gauge group and χ is assumed odd under a Z 2 symmetry and is thus stable. In the following, we will only consider a nonzero coupling for the first-generation up quark, q = u. 3 The question of the origin of the interactions in Eq. (1) might be raised. If φ is a SM singlet, we envisage two simple possibilities. For instance, φ and the SM could connect directly through the (renormalizable) Higgs portal, giving y f = m f /v × sin θ where v ≈ 246 GeV is the SM Higgs VEV breaking EW symmetry and θ is a mixing angle. In this case though, the mediator would couple to all charged fermions and its interactions with the third family would dominate. 4 In the absence of Higgs portal, φ could couple to SM fermions through interactions that involve additionnal states beyond the SM. For example, SM quarks could mix with new vector-like fermions that undergo Yukawa interactions with φ. In the limit that these additional states are heavy, characterized by a mass scale Λ, their dynamics will be captured by non-renormalizable operators like (x f /Λ)QHf R φ + h.c., with Q and f = u, d denoting respectively SU(2) L doublets and singlets. In this case y f ∼ x f v/Λ and coupling to the first family could dominate depending on the UV dynamics setting the flavor structure of the x f couplings. Moreover, taking Λ v would provide a simple rationale for the feeble couplings characterizing the freeze-in mechanism which we ought to consider in this article.
At energies below Λ QCD ≈ 200 MeV, quarks are no longer free and the φ interactions with the SM is better described in terms of hadronic resonances. We will limit ourselves here to protons, neutrons and pions. The lowenergy interaction Lagrangian is where the hadronic couplings are related to the fundamental quark coupling by matching. For first generation quarks, the coupling to nucleons is where the f N q 's are nuclear form factors whose values are extracted from matching nuclear data with lattice simulations [40,41]. For energies well below 4πf π , where f π ≈ 93 MeV is the pion decay constant, the coupling to pions can be derived from chiral perturbation theory, which gives at leading order for first generation quarks [42,43].
In order to retain the possibility of velocity-dependent cross sections for DM self-interactions required by clusters anomaly, we will consider the light mediator limit m φ m χ with m φ > 1 MeV. For mediators below the MeV scale, very strong constraints apply on scalar couplings to nucleons which basically excludes freeze-in production of DM [28].
Here the mediator is Z 2 -even and thus unstable. For m φ < 2m π the leading decay channel is in two photons through loops of charged hadrons.
which corresponds to a lifetime of τ φ ≈ 8.4×10 −11 sec/y 2 q for m φ = 1 MeV and a q = u coupling evaluated at the matching scale µ = 2 GeV. The decay of such a light, long-lived mediator would typically alter big bang nucleosynthesis (BBN), unless the mediator decouples and decays before it starts. At m φ = 1 MeV, the model is in tension with BBN for y q O(10 −5 ), greatly restricting the region of parameter space favored by the freeze-in mechanism. This strong constraint can be easily evaded by, for instance, introducing in L int an additional interaction of the mediator to neutrinos, −y ν φνν. The decay width into neutrinos Γ φ→νν = y 2 ν m φ /8π can be sufficiently large to avoid BBN constraints with a relatively small coupling y ν ∼ O(10 −10 ) which has no significant impact on the DM phenomenology. Note that the mediator has to decay into neutrinos sufficiently early so that most of them thermalize before neutrino decoupling at T ≈ few MeV.
In order to avoid strong constraints from BBN while maximizing the effect of the light mediator on DM phenomenology we conservatively set m φ = 10 MeV in the remainder of this article. Moreover, we focus on DM in the 1 − 100 GeV mass range where significant DM-nucleus scattering signals in next-generation direct detection experiments are expected.

III. MEDIATOR FREEZE-IN PRODUCTION
The mediator contribution to DM production γ φφ→χχ requires knownledge of the phase-space distribution of φ particles, f φ (p, T ). The latter is obtained from solving the (unintegrated) Boltzmann equation where p and E = (m 2 φ + p 2 ) 1/2 are respectively the 3momentum and energy of φ in the frame of the thermal bath, H is the Hubble rate and C[f φ ] is the collision term. Solving Eq. (6) is numerically challenging partly because of the ∂ p term which account for the momentum change due Hubble expansion. It is however possible to factor out this effect by introducing the dimensionless variable where h eff (T ) is the number of degrees of freedom contributing to the entropy density s(T ) = 2π 2 /45h eff T 3 and T 0 ≈ 2.3 × 10 −4 eV is the photon temperature today. The entropy ratio in Eq. (7) further accounts for the slowdown of the Hubble rate due to the decoupling of species across the cosmic history.
In terms of this reduced momentum variable Eq. (6) is brought to a one-derivative differential equation here expressed in terms of . This equation can be solved for fixed q.
The other complication lies in the form of the collision term C[f φ ]. Several interactions contribute to bring the mediator in thermal contact with the SM bath. The dominant contribution arises from QCD processes with a single φ in the final state. Those are gq → qφ and qq → gφ as well as g → qqφ once thermal corrections are included (see section III B). We also include subdominant electromagnetic processes qγ → qφ and qq → γφ which contributes at the O(10%) level. Pair production of φ is suppressed by a factor of O(4πy 2 q /α s ), where α s is the QCD coupling, and is therefore negligible in the limit y q 1 that is required for freeze-in. Hence, we have is the collision term associated with the process i. A sum over V = g, γ is implicit in Eq. (9) and the factor 2 multiplying the second term on the right-hand side accounts for the charge-conjugated processqV ↔qφ. Note that the last term is only sourced by thermal plasma effects.

A. Unintegrated collision rates
The collision terms for 2 → 2 and 1 → 3 processes are generally expressed as (upper/lower sign applies to fermions/bosons) where M denotes the scattering amplitude with initial and final state spins and colors summed over, the in (out) label denotes all the incoming (outgoing) particles other than φ with total momentum P in(out) , the index i runs over bath particles and processes. The first term inside the bracket represents the creation contribution from the process in → out + φ, while the second one accounts for depletion from its reverse counterpart out + φ → in. These two contributions, respectively denoted bŷ γ in→out+φ andγ out+φ→in , are related thanks to equilibrium conditions. Indeed, since particles 1, 2 and 3 are in where E is the energy of φ in the rest frame of the plasma. Moreover, in the absence of CP violation (as in the simplified model of interest) M in→out+φ = M out+φ→in , hence we havê Then, the generic expression for the collision terms in Eq. (8) simplifies tô where is the rate of the reaction in ↔ out + φ. Note that rates can be evaluated considering either the creation or depletion process of φ, thanks to equilibrium of SM particles. Finally, in the limit where this rate is much faster than the Hubble rate, R/H 1, the bracket in Eq. (13) goes to zero, meaning that φ reaches thermal equilibrium with the SM, f φ → (e E/T − 1) −1 .
For 2 → 2 scattering, it is more convenient to consider depletion processes (with initial state φ) to calculate the associated rates. Indeed, in this case, neglecting Pauli blocking and stimulated emission effects, ie. taking ( 1,γ 3φ→12 admits a simple expression in terms of the scattering cross-section σ 3φ→12 and, where v Møl is the Møller velocity and g i is the number of degrees of freedom of particle i. On the other hand, considering creation processes (with final state φ) is more convenient for decay. Within the same approximation, the collision term for the 1 → 3 process is expressed in terms of the differential (partial) decay width dΓ 1→23φ /d 3 p in the frame of the thermal bath, yielding We refer the interested reader to Appendix A for a fully detailed evaluation of the integrals in Eqs. (15) and (16). Note that integratingγ over the φ phase space yields d 3 p/(2π) 3γ 3φ→12 = σ 3φ→12 v n 3 n φ and d 3 p/(2π) 3γ 1→23φ = Γ 1→23φ n 1 , where · · · denotes thermal averaging and n the number density The reaction rates depend on temperature and rescaled momentum q, and they are typically faster for low momenta. For instance, as we show in Appendix A, the rate of 2 → 2 reactions approximately scales like 1/q at large q. As a result, energetic φ particles, whose momentum is larger than temperature, take more time to reach equilibrium relative to less energetic ones. For sake of illustration the rates of all relevant processes, including important plasma effects that we discuss in the next section, are shown in Fig. 1 for T = 5 GeV. Rates for different temperatures show similar behavior (see Appendix A).

B. Finite temperature corrections
Thermal corrections to the collision term C[f φ ] are implemented as follows. The perturbative formulation of gauge theory in vacuum (in powers of gauge coupling) breaks down in the presence of a hot medium due to the emergence of an external scale, the temperature T of the plasma. Gauge theory at finite temperature can still be formulated consistently only with a reorganized perturbative series where a certain class of diagrams needs to be resummed [44,45] (see also Ref. [46] for a recent review). We limit ourselves here to the so-called hardthermal-loop (HTL) approximation [47] which only resums the higher-order loop diagrams associated with soft momenta ∼ gT T where g denote gauge couplings. In this approximation, fermions and gauge bosons are quasiparticles with temperature dependent masses.
Gauge bosons at finite temperature have polarizarization-dependent dispersion relations [48]. However, the propagators of transverse and longitudinal polarizations have the same pole at zero momentum, which is given by the plasma frequency related to the Debye screening of the (chromo)electric field in a medium, and only develop small differences for non zero momentum. We neglect these differences here and in the calculation of scattering amplitudes we only replace the gauge propagator by a massive one with a pole mass given by the thermal Debye mass. To leading order in gauge coupling those are [49] for gluons and for photons, where n f (T ) and n ch (T ) are the number of active (namely satisfying m T ) quark flavors and charged particles in the plasma, respectively.
Quark dispersion relations are also modified at finite temperature with the emergence of hole excitations [50]. Nevertheless, particle and hole states are together well described by a quark propagator with a momentum independent thermal mass [51]. We therefore neglect these differences and simply add to the quark propagator the thermal mass [52], where Q q is the quark electric charge, in the calculation of scattering amplitudes. Finally, interaction vertices also receive finite temperature corrections. Those are captured to a very good approximation by renormalizing all coupling constants at the scale of the first Matsubara mode, µ R = 2πT , using renormalization group equations in vacuum [51].
There are several important implications of the plasma effects described above on mediator production. First of all, the quark thermal mass of O(g s T ) regulates the forward-enhancement of t-channel diagrams and thus strongly suppresses the production cross section compared to the zero-temperature limit. This is particularly visible at large q 1. In the opposite limit of small q 0.1 − 1, thermal masses affect the production rates in different and much more dramatic ways. For the φq → gq process, the large thermal mass of the gluon in the final state requires a highly energetic initial quark which is Boltzmann suppressed, causing the exponential drop below q ∼ 0.1. The φg → qq process however shows a strong enhancement relative to the zero-temperature limit below q ∼ 0.1. This is more conveniently understood considering the direct process qq → φg. If the light φ particle emitted from the intial quark states is sufficiently soft, it becomes possible, since m g > 2m q within the plasma, that the gluon produced from the annihilating qq pair go on-shell, which strongly enhances the 2 → 2 scattering amplitude. These effects of thermal masses are illustrated in Fig. 2 which shows the ratio of the reaction rates in Eq. 15, calculated with and without thermal masses. 5 Second, since the gluon thermal mass is always larger than twice that of the quark, opening a new production channel from the decay g → qqφ which is absent at zerotemperature. Note that the photon thermal mass, which emerges from QED interactions with charged particles in the plasma, is always too small to allow for the decay γ → qqφ.
The impact of plasma effects on the φ distribution resulting from Eq. (6) is illustrated on Fig. 3 for T = 5 GeV, where the solid black (orange) line represents the distribution with (without) including thermal corrections. For 2 → 2 processes the inclusion of thermal masses suppresses by a factor of O(10) the production of φ particles with momentum larger than temperature. At small momenta, however, thermal masses strongly enhances φ production due to significant emission of soft φ particles together with gluons. This enhancement in the production rate allows φ to reach equilibrium at small momenta much faster than in the absence of thermal corrections since, in this case, the production peaks at momenta of ∼ T . Consequently, the small momentum region is less populated, as shown with the orange line in Fig. ??. The gluon decay contribution is typically much smaller than the scattering contributions and quickly becomes inefficient in producing energetic φ particles due to phase space limitation. The latter contribution drops from about 30% at p ∼ T to roughly 10% at high momenta.

C. Simplifying assumptions
The full resolution of the Boltzmann equations is rather cumbersome. However a relatively accurate result for the relic density can be obtained upon making several simplifying assumptions.
First of all, we evaluate the impact of the choice of statistical distributions to describe particles in the plasma. The convenient assumption that particles follow a Maxwell-Boltzmann (MB) distribution is typically not justified for the freeze-in mechanism where most DM particles are produced from collisions of very relativistic particles. As shown in Fig. 3, making the approximation that all particles have MB distributions would overestimate the production of φ particles by more than a factor 2 for p ≈ T and by about 10% for much larger momenta. Note that using the MB distribution does not have a strong impact on the gluon decay contribution [36].
Second, a simple approximation would be to assume that φ is in kinetic equilibrium with the thermal bath [34]. In this case f φ /f eq is independent of momentum and simply given by the ratio n φ /n eq of the φ number density, thus avoiding having to solve the unintegrated Boltzmann equation in Eq. (8). This approximation is not justifed a priori, unless n φ ≈ n eq , because there is no elastic scattering rate between φ and SM particles that is faster than the Hubble rate. Moreover, since φ production rates are faster at low momentum, the kinetic equilibrium approximation largely overestimates (underestimates) f φ at high (low) momenta by several orders of magnitude, as shown in Fig. 3 with the horizontal red dot-dashed line. As argued in Appendix B the peak of DM production through fusion of out-of-equilibrium φ pairs occurs for one φ particle with a large momentum of O(m χ ) colliding with another one nearly at rest. Therefore, within the kinetic equilibrium assumption, there is a large compensation between the φ distribution at small and large q. As a result of this partial cancellation the kinetic approximation allows to estimate the DM relic density from out-of-equilibrium mediator fusion within an O(1) factor (see below).

IV. DM PRODUCTION
We assume negligible initial abundance for the dark sector at the end of inflation, n χ = n φ = 0 at T = T R with T R denoting the reheating temperature. In constrast with thermal production, the DM relic is produced by the so-called freeze-in mechanism [6,7] through feeble interactions with the thermal bath (during the radiationdominated era). There are two possible channels for DM production: qq → χχ (with φ in the s−channel) and φφ → χχ (with φ in the t, u− channels) where the φ density is produced from thermal collisions of SM fields (see below).
The DM yield Y χ ≡ n χ /s where n χ is the DM number density and s is the entropy density associated with the SM degrees of freedom, is governed by the following Boltzmann equation, being the Hubble parameter, and γ A↔B ≡ γ A→B −γ B→A . The γ's are the so-called (integrated) collision terms associated with the production processes described above (A = qq or φφ) and their depletion counterparts. The total DM energy density today is obtained from integrating Eq. (20) between T R and T 0 ∼ 2.7 K (the photon temperature today) where s 0 ≈ 2.89 × 10 9 /m −3 and ρ c ≈ 10.54h 2 GeV/m −3 are today's entropy and critical energy densities of the universe, respectively. h ≈ 0.674(5) [53] is related to the value of the Hubble parameter today as H 0 = 100h km/sec/Mpc. The value of T R is somewhat arbitrary. The simplified model under consideration being only valid below the EW scale we set T R = 100 GeV for consistency. Higher values of T R would require to embed the interaction Lagrangian in Eq. (1) into a specific UV complete theory respecting the SU(2) L ×U(1) Y invariance of the SM. Note that for m χ T R , DM is dominantly produced at much lower temperatures T ∼ m χ where the relevant dynamics is well described by Eq. (1) and the precise value of T R irrelevant. However, production of heavier DM particles would be strongly suppressed. Nevertheless, the freeze-in mechanism for m χ 50 GeV is well covered by direct detection [9] and most probably excluded by Xenon1T [54], see also section V.
We solve Eq. (20) neglecting reverse processes where DM annihilates back into qq and φφ. This is certainly a justified approximation for the hadronic channel since γ χχ→qq /γ qq→χχ ∼ O(n 2 χ /n 2 eq ) and n χ n eq at all times. The situation is less clear for mediator channel though, in particular because, as we show below, DM could be efficiently produced also in the case that φ is not in equilibrium with the thermal bath and n φ n eq . We verified numerically that the number densityn φ of φ particles with energy above m χ is larger than n χ by a factor of O(10 3 ) or more in regions of parameter space where φφ → χχ is dominant. Hence γ χχ→φφ /γ φφ→χχ ∼ O(n 2 χ /n 2 φ ) 1 whenever relevant and the reverse process is also negligible in this case.
For illustration, we show in Fig. 4 the "phase diagram" in the y q − y χ coupling plane resulting from the calculation of the relic density of 5 GeV mass DM in the model described in the previous section. One distinguishes three different regimes for DM production depending on the value of the quark-mediator coupling.
For relatively large values of y q the dominant DM production mechanism is directly from collisions of thermal SM particles, through qq → χχ whose cross section (as well as Ω χ ) scales like (y q y χ ) 2 .
As the quark coupling is decreased, SM collisions are less and less frequently producing DM particles and below a critical value of y crit q (the precise value slightly varies with m χ ) collisions of mediator particles become the dominant production mechanism. Since φ is assumed with negligible initial density, the efficiency of this process is determined by how much φ particles are produced from SM collisions. For values not too far below y crit q , the quark coupling is typically still sufficiently large so CHARM K→πϕ n-scat.

BBN SN1987a
10 -9 10 -7 10 -5 10 -3 10 -9 that the mediator reaches equilibrium with the thermal bath before DM production effectively starts at T ∼ m χ . In that case, the density of φ no longer depends on y q and Ω χ scales like y 4 χ . Hence the plateau in the phase diagram of Fig. 4.
For even smaller values of y q y eq q , where y eq q is the minimal coupling value needed to keep φ in thermal equilibrium with SM, the φ production rate is too slow, such that the mediator is out-of-equilibrium during DM production. In this case, the φφ → χχ rate is suppressed by the (square of the) small density of non-thermal φ. Nevertheless, this mechanism still dominates over direct production from SM collisions. This is understood as follows. Far from equilibrium, ie. for y q y eq q , the momentum distribution of φ is proportional to the production rate and scales as f φ ∼ (y q /y eq q ) 2 f eq since the rate is dominated by single production processes. As a result, the φφ → χχ contribution to the relic density scales as (y q y χ /y eq q ) 4 in this regime, and the ratio of collision terms in Eq. (20) is (roughly) γ φφ→χχ /γ qq→χχ ∼ (y q y χ ) 2 /(y eq q ) 4 . The qq-dominated freeze-in mechanism requires y q y χ ∼ 10 −11 in order to reproduce the observed DM relic density, while typically y eq q ≈ 10 −7 . Hence, the freeze-in production of DM from collisions of non-thermal φ particles dominates over the direct contribution from SM collisions by a factor of ∼ 10 5 .
To conclude this section, we stress the importance of using the full solution of the Boltzmann equation for φ in computing the DM relic density. Detailed comparison reveals that assuming a kinetic equilibrium distribution for the mediators overestimates Ω χ by a factor ∼ 2 for y q < 10 −9 while the discrepancy with the full calculation rapidly disappears for coupling values large enough so that φ approaches thermal equilibirum. As noted in Section III, using Maxwell-Boltzmann distributions leads to an overproduction of φ particles. For m χ = 5 GeV, this overestimates the relic density by about 50% in a regime where φ is out-of-equilibirum. The neglect of finite temperature corrections and plasma effects also leads to an overproduction of φ particles that yields an O(1) increase in the relic density. On the other hand, if φ particles are in thermal equilibrium, the relic density is roughly 40% higher when using the Bose-Einstein statistical distribution. Finally, when DM is mainly produced from u-quarks, its relic density increases more mildly, around 25%, when using a MB distribution.

V. PREDICTIONS FOR DIRECT DETECTION EXPERIMENTS
A relic of DM particles of mass above the GeV-scale can be directly detected by observing scattering events on heavy nuclei [55]. For a scalar mediator the expected signal is spin-independent (SI) with differential rate dR/dE R as a function of the nuclear recoil energy E R given by [56], where q 2 ≡ √ 2m N E R is the momentum transfered, ρ 0 = 0.3GeV/cm 3 is the DM energy density today, m N is the mass of the target nucleus, µ χN ≡ m χ m N /(m χ + m N ) is the reduced mass of the DM-nucleus system and N A is the Avogadro constant. F (q) is a nuclear form factor that describes the loss of coherence among nucleons at finite momentum transfer, while η(q 2 ) captures the dependence on the DM velocity distribution. Their explicit functions are given in Appendix C.
σ SI is the SI DM-nucleus scattering cross section evaluated at zero momentum transfer. In the limit of isospin symmetry it is related to the cross section on a single nucleon, say proton, asσ SI /µ 2 χN = A 2σ p SI /µ 2 χp where A is the total number of nucleons in the target and µ χp = m χ m p /(m χ + m p ) is the reduced mass of DM and the proton. In our simplified model, assuming m φ m χ , we haveσ where y p is defined in Eq. 3. Finally, the last term on the right-hand side of Eq. (22) parameterizes the t-channel propagator of the mediator. For q 2 m 2 φ , the DM-nucleus scattering is well described by a contact interaction, which is the implicit assumption behind the limit onσ p SI (as a function of m χ ) presented by all DD experiments. However for m φ q max ∼ O(GeV), the q 2 -dependence of the cross section is not negligible, and limits assuming contact interactions no longer apply. Nonetheless, the DD sensitivity for light mediators can be estimated by recasting existing limits based on the event rate expected from DM scattering in a given experiment. We have followed the recasting procedure of mi-crOMEGAs [57] for Xenon1T [54] and DarkSide50 [58]. In the low mass region, the latter is superseded by two analyses from Xenon1T using the S2 signal only [59] and taking advantage of the Migdal effect [60]. To estimate the projected sensitivity of SuperCDMS [61], the expected event rate is computed using where (E R ) denotes the detection efficiency. We assume that the efficiency vanishes below the nuclear energy threshold of 0.04keV and increases linearly to reach 85% at 2keV, for higher energies we take a constant efficiency [61]. Note that the exact shape of the efficiency curve at low nuclear recoil energies strongly affects the event rate, since the energy distribution for a light mediator peaks at low energies. The freeze-in prediction for the SI cross-section strongly depends on whether the relic abundance is dominated by qq-initiated or φφ-initiated collisions in the early universe. In the first case, both Ω χ andσ p SI depend on the same combination of couplings, (y q y χ ) 2 , such that the relic density uniquely determines the direct detection signal for a fixed DM mass. This prediction is represented by the upper black line on Fig. 5. In the second case, when mediator collisions dominate the freeze-in production of DM, the product of couplings y q y χ could be much smaller, and its value depends on whether the mediator is in thermal equilibrium with the SM or not, see Fig. 4. When freeze-in is dominated by non-thermal φ, Ω χ ∝ (y q y χ ) 4 and the product y q y χ is also fixed. Thus, the relic density also makes a unique prediction forσ p SI in this case, which is represented by the lower black line on Fig. 5. Conversely, when the mediator is in thermal equilibrium during DM production, Ω ∝ y 4 χ while the quark coupling is in the range y eq q y q y crit q . Hence, the direct detection cross section predicted by the relic density is not unique, but rather lies within the entire interval between the predicted value of qq-dominated regime (above) and that of the φφ-dominated one with non-thermal φ (below). This results in the gray-shaded band shown in Fig. 5 which approximately spans three to five orders of magnitude, depending on the DM mass. Note that theσ p SI range predicted by freeze-in is narrower for larger values of m χ , which is simply due to the fact that the mediator requires a larger coupling to the SM in order to reach equilibrium before DM production starts.
The predicted direct detection signals are already well covered by existing experiments. In particular, the current limits from Xenon1T rule out the parameter space consistent with freeze-in for DM masses above 30 GeV and partly covers the φφ-dominated freeze-in down to its threshold sensitivity, corresponding to m χ ≈ 6 GeV. For lower DM masses, the new analyses by Xenon1T based on the Migdal effect or exploiting the S2 signal only exclude the qq-dominated freeze-in as well as part of the parameter space of the φφ-dominated regime. This region is also excluded partly by DarkSide50. Moreover, future experiments will significantly improve the coverage for light DM. For instance, the projected reach for SuperCDMS will allow to probe a significant fraction of the freeze-in prediction below 5 GeV.
Finally, several terrestrial, astrophysics and cosmology constraints, which we summarise in the next section for completeness, can be imposed on a MeV-scale hadrophilic scalar. Imposing all these constraints at face value for m φ = 10 MeV severely restricts the range for the quark-mediator coupling, such that most of the mediatordominated regime would be excluded. The narrower band of the SI cross-sections that reproduce the relic abundance and satisfy these constraints corresponds to the upper part of the light grey area above the BBN line in Fig. 5. This region is significantly enlarged to the whole light grey area when the BBN bound is evaded in the presence of an additional decay channel for the mediator into neutrinos.

VI. OTHER CONSTRAINTS
Other constraints exist on individual couplings of a light mediator to SM quarks, both from terrestrial experiments and astrophysics and cosmological observations. See Refs. [28,37] for detailed reviews. For completeness we quickly decribe below the constraints shown in Figs. 4,5 that are relevant to our scenario.
Several laboratory experiments are sensitive to light particles coupled to quarks. Neutron-scattering experiments at low energies are sensitive to the coupling of φ to neutrons. The strongest constraint for a 10 MeV mass [62] arises from analyzing the momentum distribution of keV-scale neutrons scattered off lead nuclei, giving y n 1.5 × 10 −3 [63]. This bound translates to y q (2 GeV) 2.9 × 10 −4 for a scalar coupled to u-quark only [40,41].
Additional constraints come from rare meson decays. Light mediators coupled to quarks can be produced onshell in B → Kφ and K → πφ decays [64]. For the coupling values we envisage here φ is stable on collider scales and would appear as missing energy in the decays. Due to the heavy bottom mass, the B → K transition is induced at one-loop by an electroweak penguin diagram. This penguin is suppressed by the light u-quark mass and small CKM matrix elements V us V ub , which makes it negligible given the current experimental bound on such decay [65]. However, the K → π transition receives a less suppressed tree-level contribution from the chiral Lagrangian [37]. The limit on the BR(K + → π + νν) [66] provides a strong constraint from meson decays, yielding y q (2 GeV) 4.2 × 10 −6 for u-quark. The forthcoming NA62 experiment is expected to improve the sensitivity in the K + → π + νν channel by a factor of ∼ 3 [67]. Light mediators coupled to quarks are also constrained from proton beam dump experiments. In particular, the axion-like particle search at the 400 GeV SPS by the CHARM collaboration [68] can be used to constrain the process η → πφ where φ decay into two photons [37]. This search yields the strongest upper limit from meson decay, y q (2 GeV) 2.8 × 10 −6 for u-quark. Note however that this constraint can be relaxed if φ is allowed to decay into an invisible channel, like neutrinos.
Light bosons coupled to nucleons can be emitted in stars. Below a critical coupling value, the emitted bosons interact so weakly with the stellar medium that they escape the star without being reabsorbed, thus contributing to its cooling. Lack of evidence of such additional energy loss mechanisms in several stellar systems thus constrains the coupling of light bosons [69]. For large enough couplings, the new bosons are efficiently reabsorbed and trapped in the stellar medium, thus no longer contributing to energy losses. Horizontal branch and red giant stars are too cold to emit 10 MeV-scale bosons. Those can however be constrained from supernova 1987A (SN1987A) whose temperature reached T ∼ 30 MeV, excluding u-quark coupling values in the range 4.2 × 10 −11 y q (2 GeV) 1.4 × 10 −8 .
Light bosons with tiny coupling to SM fields typically live long enough to leave traces in well-understood late cosmological phenomena, such as BBN or the cosmic microwave background (CMB). If the mediator survives un-til BBN, its φ → γγ decay could inject electromagnetic energy in the thermal bath, hence increasing its entropy density and (if the decay products are sufficiently energetic) dissociating the freshly formed light elements. Given its very small scattering cross-section with the SM, the mediator will decouple relativistically from the thermal bath at T ∼ m φ = 10 MeV. Then, in order to avoid strong alteration of the standard BBN predictions for the abundances of light elements, its relic must decay away before the onset of the first nuclear reactions at t 1 sec. From Eq. (5), this implies y q (2 GeV) 2.9 × 10 −7 . Note that this constraint can be evaded by shortening the lifetime of the mediator through either increasing its mass or opening an additionnal decay channel into neutrinos. In the latter case, we checked, using the alterBBN [70] code, that the neutrinos produced from the decay of a 10 MeV mediator thermalize before neutrino decoupling and do not spoil BBN predictions.

VII. DISCUSSION
In the above analysis, we concentrated on the specific case of a 10 MeV scalar mediator that couples to DM and u-quarks. We also implicitly assumed a coupling to SM neutrinos whenever necessary to avoid cosmological constraints. However, the mechanism of DM freeze-in production from fusion of out-of-equilibrium mediators is more generic. First of all, our results for the DM relic density would equally apply in cases where φ couples to any of the light quark flavors. Moreover, DM phenomenology remain valid as long as the mediator mass is below ∼ 100 MeV. Indeed such a light mediator has little impact on the relic density, since m χ m φ , as well as on predictions and limits from direct detection, since the scaling of the cross sectionσ p SI ∝ m −4 φ holds in this mass range. The mass of the mediator and its possibility to decay into a neutrino channel do, however, affect other constraints, most notably those from BBN. In the absence of the neutrino channel, a mediator's lifetime τ φ 1 sec can be achieved with a large enough coupling y q , more precisely for y q (m φ /10 MeV) 3/2 > 2.9 × 10 −7 . While BBN constraints restrict the out-of-equilibrium regime for lighter mediators, m φ > 30 MeV allows for sufficiently smaller values of y q such that the out-ofequilibrium regime opens up significantly. This also implies a wider region in direct detection of light dark matter, with m χ 6 GeV, that is free of constraints. Note also that for m φ > 30 MeV, the supernova constraint disappears thus further relaxing constraints on the whole out-of-equilibrium window.
In summary, we performed a detailed calculation of DM production via the production of out-of-equilibrium mediators by solving the unintegrated Boltzmann equation for the latter, instead of making the kinetic equilibrium approximation, and including the effect of thermal masses and of quantum statistics. Each of these effect has a large impact on the prediction of the dark matter relic density. We also showed using a simplified model that while this mechanism faces cosmological constraint it can be probed by direct detection experiments. Increasing the sensitivity of direct detection experiments at low masses is however crucial to completely cover the pure freeze-in via out of equilibrium region in the future.
For fixed s and E, E 3 reaches its extremal values when particle 3 and φ are collinear. In this case, E 3 relate to s and E through a Lorentz boost transformation, where E 3 * = (s + m 2 3 − m 2 φ )/2 √ s and p 3 * = −p * are the energy and momentum of particle 3 in the center-of-mass frame and y is the rapidity of the center-of-mass, which satisfies E = E * cosh(y)+p * sinh(y) with E * = √ s−E 3 * . There are two independent solutions, where y + (y − ) is reached when the momenta of particle 3 and φ are parallel (antiparallel), and the energy of particle 3 in the plasma frame is E max 3 (E min 3 ).

Decay
The rate of the 3-body decay reaction 1 → 23φ is see Eq. (16), where Γ 1→23φ is the decay width of particle 1 in the plasma frame, p is the φ momentum in that frame and In order to determine dΓ 1→23φ /dE, we will first calculate the differential width in the rest frame of the decaying particle 1, dΓ * 1→23φ /dE * , and boost it to plasma frame. First of all, the width is boosted by a Lorentz factor dΓ 1→23φ = m 1 /E 1 × dΓ * 1→23φ . Then, one needs to find how many φ particle with energy E * in the rest frame of particle 1 wind up with energy E in the plasma frame. The energies E and E * are related by a boost transformation E = E * cosh(y) − p * cos θ * sinh(y) , with rapidity y = 1/2 log[(E 1 + p 1 )/(E 1 − p 1 )]. θ * is the angle between the φ momentum in the particle 1 rest frame and the boost direction, given by particle 1 momentum. Assuming the differential width is flat in cos θ * , which is the case for the decay of a scalar or an unpolarized particle with non-zero spin, a δ-distribution dΓ/dE * = Aδ(E * −E 0 * ) in the particle 1 rest frame yields a rectangular function in terms of E in the plasma frame, dΓ/dE = A/N [Θ(E−E min )−Θ(E−E max )] whose boundaries are defined by Eq. (A7) upon setting cos θ * = ±1 and normalization is rescaled by the size of the rectangle N = E max − E min = 2(E 2 0 * − m 2 φ ) 1/2 sinh(y). Since any generic spectrum in E * can be decomposed as a (infinite) set of δ-function peaking at different values of E * , the boosted spectrum is simply given by adding up the rectangles, yielding p 1 ∼ p 2 ∼ p + /2, or one particle carries all the total required momentum, p 1 ∼ p + and p 2 ∼ 0. In the first regime, η/η eq ∼ [1 − exp(−2cM Pl /p + )] 2 ∼ (2cM Pl /p + ) 2 , where in the last expression we used p + cM Pl , typically valid out-of-equilibrium. Conversely, in the second regime the suppression is less severe, η/η eq ∼ 1 − exp(−cM Pl /p + ) ∼ cM Pl /p + , because slow particles are close to equilibrium. Hence, production of heavy DM from much lighter, out-of-equilibrium φ fusion is dominated by fast particles colliding slow ones.