$Z^\prime$-mediated dark matter freeze-in at stronger coupling
Giorgio Arcadi David Cabo-Almeida Oleg Lebedev


September 9, 2024

 



a,b]Giorgio Arcadi,
a,b,c,d]David Cabo-Almeida,
e]Oleg Lebedev 

[a]Dipartimento di Scienze Matematiche e Informatiche, Scienze Fisiche e Scienze della Terra, 
 Universita degli Studi di Messina, Viale Ferdinando Stagno d'Alcontres 31, I-98166 Messina, Italy
[b]INFN Sezione di Catania, Via Santa Sofia 64, I-95123 Catania, Italy
[c]Departament de Física Quàntica i Astrofísica, Universitat de Barcelona,
Martí i Franquès 1,
 E08028 Barcelona, Spain
[d]Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona,
Martí i Franquès 1,
 E08028 Barcelona, Spain
[e]Department of Physics and Helsinki Institute of Physics,
Gustaf Hällströmin katu 2a, 

FI-00014 Helsinki, Finland

giorgio.arcadi@unime.it
david.cabo@ct.infn.it
oleg.lebedev@helsinki.fi
 


 
We study freeze-in production of fermionic dark matter mediated by a Z^' gauge boson. In particular, we explore the regime of Boltzmann-suppressed production, when the Standard Model (SM) thermal bath
temperature never exceeds the dark matter mass. The corresponding  gauge coupling is then required to be significant, up to order one. As a result, this class of freeze-in models can be probed by the current and future direct  
dark matter detection experiments. 









Z^'-mediated dark matter 
 freeze-in  
           at stronger coupling
    [
    September 9, 2024
======================================================================




 



  

§ INTRODUCTION

   
   Cosmological dark matter remains one of the most mysterious entities in modern physics. The most common approach to this problem is to postulate the existence of a stable particle which has no quantum numbers with respect to the Standard Model
   gauge symmetries and whose gravitational effects we observe. Its further properties depend on the coupling to the observed particles as well as itself. If such couplings are very small, the dark matter field never reaches thermal equilibrium and its abundance accumulates over time due to 
   its slow production by the SM thermal bath <cit.>. This  mechanism is known as “freeze-in” <cit.>.
   
   The freeze-in mechanism is predictive only if other production channels are negligible. Generally, there are further non-thermal DM production mechanisms which are operative in the Early Universe <cit.>. Most importantly, all particles are produced by gravity itself during and after inflation. For example, particle production immediately after inflation, during the inflaton oscillation epoch, can be very efficient due to the presence of the Planck-suppressed operators that couple the inflaton  ϕ to DM such as <cit.>
   
    1 M_ Pl^2  ϕ^4 s^2  , 1 M_ Pl ϕ^2 χ̅χ ,

where s is scalar DM and χ is fermion DM.
   These operators are generated by both classical and quantum gravitational effects. The resulting abundance of DM exceeds the observed value by many orders of magnitude unless the corresponding Wilson coefficients happen to be very small or DM is extremely light
    <cit.>.
   
   One way to address this problem is to allow for a long inflaton-dominated expansion period before reheating. The  DM quanta produced immediately after inflation  are relativistic and their energy density gets diluted by expansion in a non-relativistic background. At the reheating stage, the 
   resulting DM abundance  is then small  <cit.>.
   This implies  that the
   reheating temperature T_R is  relatively low, being only limited from below, T_R > 4MeV,  by observations <cit.>. It opens the possibility that the SM bath temperature has never been higher than the DM mass and the freeze-in production is Boltzmann suppressed <cit.>.
   The SM-DM coupling required for producing the correct DM relic abundance can then be as large as order one, which facilitates direct searches for dark matter   via DM-nucleon scattering and  in collider experiments. 
   
   
   The purpose of our current work is to explore this idea in the context of a Z^' extension of the Standard Model <cit.>. The role of dark matter is played by a Dirac fermion which couples to a Z^' and whose production by the SM thermal bath is mediated by the
   Z^'. We delineate parameter space leading to the correct DM relic abundance, while satisfying   the current direct detection and collider constraints. We find that the Z^'-mediated freeze-in can be probed by direct DM detection experiments as long as
   the Z^' has significant “vector” couplings to the SM fermions. 
   
   
   

 §.§ The model

   
   We study an extension of the Standard Model with a massive Z^' boson and a Dirac fermion χ, which has no SM quantum numbers and plays the role of dark matter. The relevant couplings are
   
    ℒ⊃-m_χχ̅χ-1/2 M^2_Z^' Z_μ^' Z^'μ+χ̅γ^μ(V_χ-A_χγ_5) χ Z_μ^'+∑_f f̅γ^μ(V_f-A_f γ_5) f  Z_μ^' ,

Here f represents the SM fermions;  V_f,χ and A_f,χ are the vectorial and axial couplings. These are related to the Z^' gauge coupling  g_Z^'  and the charges of the left- and right-handed fermions   X_f_L,X_f_R      as 
 V_f=g_Z^'(X_f_L+X_f_R)/2 and A_f=g_Z^'(X_f_L-X_f_R)/2.  

   
   We assume that the Z^' is the only portal between the SM and dark matter, focussing on the TeV scale M_Z^'. 
      We also  take its mixing with the SM gauge bosons to be negligible for our purposes. The U(1) extensions of the SM are subject to the anomaly cancellation constraint, which generally necessitates the presence of additional fermions. Such fermions, however, can be vector-like with respect to the SM gauge charges and hence  have large masses making their impact on the DM phenomenology insignificant. A similar approach was taken in related work <cit.> (see also <cit.>).
      
   
   In what follows, we study freeze-in production of dark matter χ mediated by the Z^'. Our main assumption is that the SM thermal bath temperature T is always below m_χ. This means, in particular, that the reheating temperature T_R is sufficiently low,
   allowing for dilution of the gravitationally produced relics, as explained above.
   To compute the DM production, we resort to the instant reheating approximation, e.g. take  T to increase abruptly from zero to T_R and then decrease as usual, in accordance with the SM entropy conservation. 
   This gives a good estimate of the DM abundance in cosmological models with a flat SM temperature profile before reheating <cit.>. 
   Using appropriate rescaling and replacing T_R with the maximal SM sector temperature, we can extend our results to a broader  class of models.             
   Earlier work on the Z^'-mediated freeze-in, which assumes a high reheating temperature, can be found in <cit.>.
  
   



§ DARK MATTER  RELIC ABUNDANCE

   
   In this section, we study freeze-in dark matter production in the Boltzmann-suppressed regime T≪ m_χ. 
   Dark matter  is produced via annihilation of the SM thermal bath particles, e.g.  f̅ f → Z^'→χ̅χ, and its density n= n_χ + n_χ̅ is controlled by the Boltzmann equation,
   
    ṅ + 3Hn =2 Γ (  SM→χ̅χ) - 2 Γ ( χ̅χ→ SM )  ,

   where the factor of 2 accounts for production of 2 DM states, χ and χ̅, and Γ is the reaction rate per unit volume.
   
    Since dark matter is non-relativistic at T ≪ m_χ, the a→ b reaction rate is given by 

    Γ_a→ b = ∫( ∏_i∈ ad^3  p_i  (2 π)^3 2E_i f(p_i)) ( ∏_j∈ bd^3  p_j  (2 π)^3 2E_j)
    | M_a→ b|^2   (2π)^4 δ^4(p_a-p_b)  ,

where p_i and p_j are the initial and final state momenta, respectively, and f(p) is the momentum distribution function. 
M_a→ b is the  QFT  a → b  transition amplitude, 
in which 
  both
 the  initial and final state phase space symmetry factors are absorbed.
For the 2→ 2 reactions,  energy conservation combined with the Boltzmann statistics f(p)=e^-E/T implies

    f (p_1) f(p_2) = f(p_3) f (p_4) .

   Therefore, the DM production rate via thermal SM states can be written as the  thermal DM annihilation rate into SM quanta,
   
    Γ (  SM→χ̅χ) = Γ ( χ̅χ|_ therm→ SM)  .

 The latter can be computed following Gelmini and Gondolo
   <cit.>,
 
    Γ (χ̅χ|_ therm→ SM )   =   ⟨σ(χ̅χ|_ therm→ SM  )    v_r ⟩  n_χ n_χ̅ = 2^2 (2π)^6 ∫σ v_r   e^-E_1/T e^-E_2/T  d^3 p_1 d^3 p_2  
       =   2^3π^2 T (2π)^6 ∫_4m_χ^2^∞  d s σ ( s- 4 m_χ^2) √( s)  K_1 (√( s)/T)  ,

 where σ is the   χ̅χ→ SM    cross section; v_r is the relative velocity of the colliding quanta with energies E_1,E_2 and momenta p_1,p_2; n_χ is the χ  number density;
 ⟨ ... ⟩ denotes a thermal average;
 s is the Mandelstam variable, and
 K_1(x) is the modified Bessel function of the first kind.   We have explicitly factored out the spin d.o.f. factor 2^2 specific to fermion annihilation  (n_χ= 2  e^-E/T).
  
   
   Let us  focus for now on a heavy Z^' regime, M_Z'  > 2m_χ. It is instructive to consider in detail the limit M_Z' ≫ m_χ, m_f, where f represents the SM fermions. This allows us to obtain simple analytical approximations in the pure freeze-in limit. In our numerical studies however, we do not resort to this approximation and use the exact (tree-level) results together with the backreaction term in (<ref>).
In what follows, we  consider the vector and axial Z^' cases separately.
   
   
  

 §.§ Vector coupling

  
  In most of the parameter space, the main production/annihilation channel is 
  
    χ̅χ↔f̅ f  ,

  where f is an SM fermion and contributions of all the fermions lighter than χ should be summed.
  The DM annihilation via a vector Z^' is allowed already at the s-wave level and the corresponding cross section for m_χ≫ m_f is 
  
    σ (χ̅χ→f̅ f) ≃V_f^2 V_χ^2  m_χ^2  2π M_Z'^4 √(s  s-4m_χ^2) .

  The reaction rate is found from Eq. <ref>.
  In the regime m_χ≫ T, we may use the asymptotic form K_1 (√( s)/T) ≃√(π 2)T^1/2 s^1/4 e^-√(s) /T. 
   When computing the integral over s, one can use the following approximation: since the integrand peaks sharply close to s= 4m_χ^2, one may replace
  s→  4m_χ^2 in “slow” functions of s, while keeping the factors √(  s-4m_χ^2) and e^-√(s) /T as they are. 
   The resulting integral reduces to the Gamma function such that
   
    ∫_4m_χ^2^∞  d s √(  s-4m_χ^2)  e^-√(s) /T≃√(π) 2 (4 m_χ T)^3/2 e^-2m_χ /T .

   The resulting DM production rate via SM fermion annihilation  according to (<ref>) is
   
    Γ (f̅ f →χ̅χ) = V_f^2 V_χ^2    m_χ^5 T^3  2π^4 M_Z'^4  e^-2m_χ /T .

   
   In the pure freeze-in regime,  the Boltzmann equation (<ref>) reduces to ṅ + 3nH =2  Γ (f̅ f →χ̅χ). It can be solved analytically for simple reheating scenarios. For definiteness, 
   we resort to the instant reheating approximation, that is, we assume that the SM sector temperature increases abruptly from zero to T_R. After that it is determined by entropy conservation, as usual.
   Integrating the Boltzmann equation from T=T_R to T=0, one finds that n∝ T^3 at late times and, due to the exponential suppression of the reaction rate,  the production is dominated by the initial moments when T ≃ T_R.
   Using
   
    s_ SM= 2π^2  45 g_* T^3 ,  H = √(g_* π^2  90 ) T^2  M_ Pl ,

  where g_* is the effective number of the  SM d.o.f. and M_ Pl= 2.4 × 10^18 GeV, 
   we find
   
    Y ≡n s_ SM≃∑_f 45 √(90) 4V_f^2 V_χ^2 π^7 g_*^3/2 m_χ^4 M_ Pl  M_Z'^4 T_R  e^-2m_χ /T_R ,

   where the sum runs over all SM Dirac fermions f  with the color multiplicity properly included.[The neutrinos contribute 1/2 of the Dirac fermion contribution. For V_f = V_χ =λ and m_χ > m_t,
   the sum over all the SM fermions amounts to a factor of 22.5.]
   Here we have assumed that   g_* stays approximately constant in the regime of interest.
   For the universal couplings V_f = V_χ =λ, the observed Y=4.4 × 10^-10  GeV/m_χ requires
   
    λ≃  10^-6 M_Z' T_R^1/4 m_χ^5/4  e^m_χ/(2T_R)

   for g_* ≃ 107.
   Since m_χ≫ T_R, the exponential factor can be  large making λ as large as O (1).  Nevertheless, the produced dark matter density is low and it does not thermalize, thereby justifying the assumption of 
   the  freeze-in regime.  We note that the size of the  coupling is determined primarily by   m_χ /T_R and rather insensitive to a specific charge assignment, thus making the  universal charge limit justified.
   
   
   
   The above  result is obtained in the instant reheating approximation. As shown in <cit.>,  this gives a good estimate of the required coupling within  a larger class of models. If the SM temperature is constant before reheating, 
     Eq. <ref> stands as it is, up to a small (percent level) shift in T_R. More generally, DM production peaks at the maximal temperature T_ max such that one replaces T_R → T_ max in the above formulas
     and, to account for entropy production between T=T_ max and T=T_R, rescales the coupling with the factor (T_R/T_ max)^κ, where κ is model dependent. 
   
   The common feature of freeze-in models at stronger coupling is that the DM abundance exhibits the factor e^-2m_ DM /T_R times some power of the dark matter coupling <cit.>. 
   As in the other models, we find that an order one coupling requires m_ DM∼ 20  T_R for typical parameter values.
   
   
   
      
   
   
   
    

 §.§ Axial coupling

   
   Consider now DM production via  f̅ f →χ̅χ mediated by an axial Z^' <cit.>.
   To compute the reaction rate, we consider non-relativistic DM annihilation. In contrast to the previous case, in the m_f → 0 limit, the process is p-wave and hence less efficient.
   Neglecting the fermion masses, one finds
    
    σ (χ̅χ→f̅ f) ≃A_f^2 A_χ^2   12 π M_Z'^4 √(s(  s-4m_χ^2)) .

   This approximation applies for heavy DM, m_χ≫ m_t, or  in the regime where the top-quark channel is not available, m_χ < m_t.
   The reaction rate calculation proceeds as before, except the different velocity dependence leads to the integral 
    
    ∫_4m_χ^2^∞  d s   (s-4m_χ^2)^3/2  e^-√(s) /T≃3√(π) 4 (4 m_χ T)^5/2 e^-2m_χ /T

   instead of (<ref>).
   The result is
    
    Γ (f̅ f →χ̅χ) = A_f^2 A_χ^2    m_χ^4 T^4  2π^4 M_Z'^4  e^-2m_χ /T .

   Note that the rate is now proportional to T^4, unlike that in the vector case.
   The DM abundance is given by
   
    Y ≃∑_f 45 √(90) 4A_f^2 A_χ^2 π^7 g_*^3/2 m_χ^3 M_ Pl  M_Z'^4   e^-2m_χ /T_R .

For A_f = A_χ =λ and g_* ≃ 107, the required coupling is  

    λ≃   10^-6 M_Z' m_χ  e^m_χ/(2T_R) .


      
     
     
     
             

 §.§ Additional channels

        
      For heavier DM, m_χ∼ M_Z', the decay and t-channel processes 
      
    Z^'→χ̅χ ,  Z^' Z^'→χ̅χ

        become important.
  The Z^' number density receives  Boltzmann suppression similar  to e^-2m_χ /T, hence processes with the Z^' on-shell cannot be neglected.
      The Z^' gauge boson  maintains  thermal equilibrium  with the SM bath due to its coupling to light fermions. The reaction
      
    f̅ f ↔ Z^'

      is very efficient for the coupling range of interest and faster than the rate of expansion of the Universe (see e.g. analogous calculations in <cit.>), thereby implying thermal equilibrium. 
      
      The decay channel Z^'→χ̅χ dominates  around M_Z'∼ 2m_χ, while the efficiency of the t-channel mode Z^' Z^'→χ̅χ depends on the nature of
      the Z^' couplings. In case of the axial couplings, the fermion annihilation mode is suppressed by m_f or DM velocity. On the other hand, the processes involving the longitudinal component of a Z^'
      at high energy are enhanced by E/M_Z'. Hence, one expects  Z^' Z^'→χ̅χ to dominate for large m_χ. 
      Fig. <ref> (left) shows  the relative reaction rate contributions of the different production channels for the axial coupling case, produced with   micrOMEGAs  <cit.>.
      As expected, we observe that the t-channel annihilation becomes dominant for heavy dark matter.  In the vector coupling case, on the other hand, 
       the SM fermion annihilation mode is efficient and we find that it remains the main  production channel at large m_χ.
      
      
      
      
       
 


   


    
   
   
   

 §.§ Dark matter thermalization

   
   
   
    The above  considerations help us understand the qualitative behaviour of  λ (m_χ, T_R)  at low DM densities, when the term Γ ( χ̅χ→ SM ) in Eq. <ref> can be neglected.
    However, in our numerical analysis we use the micrOMEGAs tool  <cit.>
      which takes into account all the channels 
      as well as the DM annihilation effects. 
      The latter lead to a qualitative change in the relic abundance calculations at larger coupling.
        As one increases λ, the DM density grows thus enhancing DM annihilation and eventually leading to its thermalization.
          This manifests itself in the freeze-in lines  λ (m_χ, T_R)   merging with the freeze-out curve
            λ (m_χ), as we show in the next section. Such thermalization   has been studied in detail in <cit.> and entails a smooth FIMP-WIMP transition <cit.>.
        

Fig. <ref> (right)  illustrates the effect of the coupling increase on the eventual DM abundance. While at small λ, Y grows with the coupling,  this ceases to be the case above a certain critical coupling, which is around λ∼ 0.02 in this example.
We observe that the abundance decreases over time due to DM annihilation.
At yet larger couplings, Y(T) follows its equilibrium value Y_ eq(T) for an extended period of time, before freeze-out. This signals DM thermalization.


We note that, in order to compute  Γ ( χ̅χ→ SM ), one needs to know the momentum distribution of dark matter. This reaction becomes important at large enough coupling and, thus, it is safe to assume that the SM-DM system is in kinetic equilibrium at that stage. 
   Indeed, the  elastic scattering 
reaction        SM+χ→ SM+χ        is more efficient than the annihilation one due to the higher density of the SM states. Hence, kinetic equilibrium sets in at smaller couplings.  This is also the assumption adopted in micrOMEGAs .
 




§ CONSTRAINTS AND PARAMETER SPACE ANALYSIS


In this section, we delineate parameter space of the model  producing the correct DM relic abundance and discuss various constraints on the “universal”  coupling λ. These include  bounds on the Z^' interactions  from collider experiments as well as constraints  on dark matter from direct and indirect searches.











 §.§ LEP constraints


  A heavy Z^' can be integrated out leading to a set of 4-fermion interactions. 
 The strictest constraints on the contact interactions are   imposed by the 4-lepton operators ℓℓℓℓ, while the lepton-quark operators    ℓℓ  q q  <cit.> are less important.
 The most significant observables are the muon decay constant, which is affected by the muon-electron operators, as well as the LEP lepton measurements <cit.>. 
 The relevant constraints for the flavor-universal case can be read off from Fig. 3 (right) of <cit.>. 
 The vector coupling     (κ_L=κ_R)    is constrained by
 
    λ M_Z' < 0.12  TeV
 
 at 95 % CL. The constraint on the axial coupling (κ_L=-κ_R) is somewhat weaker,
  
    λ M_Z' < 0.19  TeV
 
 at 95 % CL.
The corresponding excluded regions  are shown in Fig. <ref>  in light orange, marked “LEP”.
   
   

 §.§ LHC constraints

   
   The LHC experiments are conducting direct searches for a new gauge boson Z^' in the channel pp →ℓ^+ ℓ^- X, 
    where X represents
the beam fragment jets. The constraints are obtained on the quantity 
    σ (pp → Z^')   BR(Z^'→ℓ^+ ℓ^- ), employing
      the narrow width approximation. We use the current data from ATLAS <cit.> based on LHC Run 2 at √(s)=13TeV  
   with a total integrated luminosity of L = 139  fb^-1.
 The theory prediction for 
   σ (pp → Z^') is  computed using 
  MadGraph <cit.>. As the first step, we reproduce the ATLAS limits on   the Z_ψ^' model <cit.>. Then, we derive the bounds on the universal coupling
   λ in our model. The results are presented in  Fig. <ref>: the light blue area marked Z^'→ℓℓ is excluded by the LHC.
   
   The typical bound on the Z^' mass for the SM-like couplings is around 5 TeV <cit.>. The constraint fades away very fast as M_Z' increases and above 6-7 TeV the LHC 
    has no sensitivity to the extra gauge boson, at least at the perturbative level. This is unlike the other bounds (LEP, direct DM detection) which scale as a function of  λ/M_Z'. 
    
    
   
    

 §.§ Direct DM detection

    
    
    
    Direct DM detection experiments traditionally set the strongest constraints on dark matter models (see, e.g. a recent discussion in <cit.>). These are based on non-observation of any significant interaction of dark matter with nucleons. 
    In our case, such an interaction is mediated by the Z^'.
    The DM-nucleon scattering occurs at low energy, hence the Z^' can be integrated out resulting in the effective Lagrangian 
  
    ℒ_eff    ⊃   1/M_Z^'^2[V^χ(2 V^u+V^d) χ̅γ^μχ p̅γ_μ p. 
          .+  A^χ(Δ_u^p  A^u+(Δ_d^p+Δ_s^p) A^d) χ̅γ^μγ^5 χ p̅γ_μγ^5 p],

   where p represents the proton, and a similar expression applies to the neutron n. Δ_i^p,n stands for the spin content of  quark i in the proton/neutron.
   According to <cit.>,

    Δ_u^p=0.842,   Δ_d^p=-0.427,   Δ_s^p=-0.085  ,

and the corresponding neutron quantities are obtained by the isospin symmetry. 
 The resulting cross section for spin-independent scattering of DM on a nucleus with charge Z and atomic weight A is
 
    σ_χ N^SI=μ_χ N^2 (   V^χ)^2 /π M_Z^'^4[V^u(1+Z/A)+V^d(2-Z/A)]^2  ,
  
  where μ_χ N  is the reduced mass of the DM-nucleon system. 
  This cross section must lie below the LZ 2022  and XENON1T bounds <cit.>, which imposes a significant constraint on the model parameters.
 For the case of the universal coupling λ, the direct detection constraint scales as 
  λ^4/M_Z'^4, while σ_χ N^SI remains independent of the DM mass as long as m_χ≫ 1GeV.
  
  The axial couplings lead to the spin-dependent DM-nucleon scattering, which is only weakly constrained
  and the resulting bounds are unimportant for our analysis.  The LHC searches impose the leading constraints in this case, at least for the M_Z' range of interest.

  
  In Fig. <ref>, we show the  direct detection constraints on the universal coupling λ with the help of  the micrOMEGAs tool.  
  The grey regions marked by “LZ+XENON1T” are excluded.[ The  LZ collaboration has updated its 2022 direct detection bound. According to the preliminary (unpublished) estimate, the bound has improved by a factor of 5 - 10, depending on the DM mass. This translates into a stronger constraint on λ by a factor of 1.5 - 1.8.]
  We observe that for a heavy Z^', these constraints supersede the other bounds, while for a lighter  Z^', the LHC
searches  set the strictest bounds. 
The figures also display the “neutrino fog”  which represents the neutrino background for the direct DM detection experiments <cit.>. The standard techniques would not be able to distinguish the DM scattering from that of neutrinos,
if λ falls below the neutrino fog curve.  While the corresponding parameter space could be explored with innovative techniques, DM detection is challenging for such low couplings. 
Nevertheless, we see that large regions of the freeze-in parameter space, especially for a heavy Z^', lie above the neutrino background and thus can be probed by the current and future direct DM detection experiments
such as XENONnT <cit.> and DARWIN <cit.>.[This property is shared by a class of  freeze-in models  with a light mediator <cit.>
as well as the Higgs-portal-type  WIMP models with a first order phase transition <cit.>.  ] 

 

 

 §.§ Indirect DM detection

 
 Indirect dark matter detection is based on possible signatures of its annihilation in regions with significant DM density, e.g. the Galactic center. 
 The strictest bounds come from the 
  Fermi-LAT   <cit.> observation of 30 dSphs for 14.3 years.
  Interpreting these constraints in our model using  micrOMEGAs, we find that they are superseded by the bounds discussed above and thus are insignificant.
   
  

 §.§ Discussion and summary

  
 The correct DM relic abundance is reproduced 
  along the colored lines in Fig. <ref>. Each of them has a different reheating temperature T_R, which we also identify with the maximal SM bath temperature.  
  We observe that, at sufficiently large coupling, all  the freeze-in    lines   merge with the thermal  black curve. This signals DM thermalization such that the relic abundance is controlled by the 
  usual freeze-out. The corresponding parameter space is, however, ruled out experimentally, leaving only a very narrow resonance region m_χ≃ M_Z'/2.
  
  
  In addition to the  constraints discussed above, we impose the perturbativity bound which excludes the red hatched area.
    In this   region, the Z^' resonance becomes too broad such that Γ_Z' > 0.5   M_Z' and it ceases to be a “particle” in the conventional sense.
  While for a light Z^' this bound is superseded by other constraints, in the case of a heavy Z^', it becomes important, especially for the axial-type couplings. 
  
   
   The axial coupling case is constrained mostly by collider observables, whereas  the direct detection bounds are very weak. This also implies that an axial Z^'-induced freeze-in can hardly be probed experimentally.
   On the other hand, the vector case is more interesting and significant parts of the freeze-in parameter space can be probed by the direct detection experiments. 
   For M_Z'=1TeV, the LHC constraints dominate, yet they leave a modestly sized region above the ν-fog line, which can therefore be probed by direct DM detection. The size of this area grows as M_Z' increases
   and at M_Z'=7TeV, the LHC bounds fade away. 
     
    We observe that, as long as the Z^' couplings have a significant vector component, the direct DM detection prospects are good in the entire range of the DM mass considered, 10-10^4 GeV.  
     The required gauge coupling is in the range      10^-3-10^-1 for a Z^' mass between 1 and 10 TeV.
  
   
     
   


§ CONCLUSIONS

   
   We have analyzed Z^'-mediated fermion dark matter production in the freeze-in regime at stronger coupling, when the SM bath temperature is below the DM mass.
   Models with a low reheating temperature are motivated by the problem of gravitational particle production which mars  the usual freeze-in mechanism. In particular, Planck–suppressed
   operators coupling dark matter to the inflaton are very efficient in particle production immediately after inflation, which leads to a large initial abundance of dark matter. This problem can be solved by 
   allowing for an extended  period of inflaton-dominated expansion thereby diluting the initial DM abundance.  As a result, the reheating temperature in this framework  is relatively low, depending on further details.
 
   Using both analytic estimates and more sophisticated numerical tools, we find that the correct DM relic abundance can be produced for a broad range of the DM mass m_χ, assuming a TeV-scale Z^'. The main factor is the reheating temperature
   T_R and, for m_χ / T_R ∼ 20, the Z^' couplings can be as large as O(1).  We distinguish the vector and axial coupling cases, which exhibit different phenomenology. While the DM production rates are similar in both cases, 
   the   constraints on the parameter space differ substantially.   In particular, the direct DM detection results  set strict bounds on  the vector coupling, whereas the axial coupling remains essentially unconstrained. 
   In our freeze-in framework, this implies that the vector case can be probed by the current and future direct detection experiments, which will reach the sensitivity at the level of the “neutrino fog”.  The axial case, on the other hand, is very difficult to test due to the suppression of the DM-nucleus cross section. Further bounds on the parameter space are imposed by the direct Z^' searches at the LHC, LEP measurements of the lepton production and
   constraints from the muon decay. 
   
   As the Z^' coupling increases, the DM production grows more efficient and the inverse reaction becomes significant.  At some critical value of the coupling (depending on T_R), both reactions equilibrate and dark matter thermalizes. 
   We observe this explicitly as the freeze-in relic abundance curves merge with the standard freeze-out line. In the vector coupling case, the corresponding region of the parameter space is ruled out by the direct DM detection.
   However, at lower couplings, the Boltzmann-suppressed freeze-in is still operative and produces the correct relic abundance while evading 
    such constraints. For a range of the vector couplings, typically of order 10^-3-10^-1, the DM-nucleon cross section lies above the neutrino fog and thus can be tested in the near future.  Needless to say, this does not require the Z^'
to be purely “vectorial”, it only sets a lower bound on the vector component of the Z^' coupling. 
    
    Our main conclusion is that, as long as the vector coupling of the Z^' is substantial,  the  freeze-in dark matter can be probed further (and possibly discovered) by direct DM detection experiments such as XENONnT and DARWIN.  
   
  

   Acknowledgments
   The authors thank Francesco Costa for fruitful discussions.
   D.C.A. acknowledges funding from the Spanish MCIN/AEI/10.13039/501100011033 through grant PID2022-136224NB-C21.


   

§ GENERAL CROSS SECTIONS AND DECAY WIDTHS 
   
   
   In this appendix, we provide some analytical formulas for the general case of the vector and axial couplings being present simultaneously. 
   The DM annihilation cross section 
    into a Standard Model Dirac fermion f is given by   
   

    σ_χ̅χ→f̅ f   =   1/12 π s[(s-M_Z^'^2)^2+M_Z^'^2 Γ_Z'^2]√(1-4 m_f^2 / s/1-4 m_χ^2 / s)[A_ f ^  2 V_χ^2(s-4 m_f^2)(2 m_χ^2+s). 
          +A_ f^2A_χ^2(4 m_χ^2[m_f^2(7-6 s/M_Z^'^2+3 s^2/M_Z^'^4)-s]
    +s(s-4 m_f^2))
          .+V_f^2(2 m_f^2+s)(A_χ^2(s-4 m_χ^2)+V_χ^2(2 m_χ^2+s))] ,

 where s is the Mandelstam variable. This result agrees with the corresponding cross section presented in <cit.>.
  The Z^' decay width is
  

    Γ_Z'=Γ_Z^'→χ̅χ+∑_fΓ_Z^'→f̅ f  ,


with


    1.2Γ_Z^'→χ̅χ=√(M_Z'^2-4 m_χ^2) (A_χ^2 (M_Z'^2-4 m_χ^2)+V_χ^2 (2 m_χ^2+M_Z'^2))/12 π  M_Z'^2 ,


    1.2Γ_Z^'→f̅ f=√(M_Z'^2-4 m_f^2) (A_f^2 (M_Z'^2-4 m_f^2)+V_f^2 (2
       m_f^2+M_Z'^2))/12 π  M_Z'^2 .

 In the universal coupling case with a heavy Z^', the sum over the SM fermions amounts to a factor of 22.5, which includes 3×6 quark contributions, 3 charged lepton contributions and 3×1/2 neutrino terms.  
   The branching ratio for the Z^'→ e^+ e^- decay approaches 1/23.5 in this limit.

The t-channel Z^' Z^'→χ̅χ analytical results can be found in <cit.>, which we have also reproduced.   




 
99



Dodelson:1993je 
  S. Dodelson and L. M. Widrow,
  
  Phys. Rev. Lett.  72, 17 (1994).
   
  


 
Hall:2009bx
L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West,

JHEP 03, 080 (2010).




 
Lebedev:2022cic
O. Lebedev,

JCAP 02, 032 (2023).



 
Lebedev:2022ljz
O. Lebedev and J. H. Yoon,

JCAP 07, no.07, 001 (2022).


 
Koutroulis:2023fgp
F. Koutroulis, O. Lebedev and S. Pokorski,

JHEP 04, 027 (2024).




Hannestad:2004px
S. Hannestad,

Phys. Rev. D 70 (2004), 043506.




 
Cosme:2023xpa
C. Cosme, F. Costa and O. Lebedev,

Phys. Rev. D 109, no.7, 075038 (2024).




Langacker:2008yv
P. Langacker,

Rev. Mod. Phys. 81 (2009), 1199-1228










Arcadi:2013qia
G. Arcadi, Y. Mambrini, M. H. G. Tytgat and B. Zaldivar,

JHEP 03 (2014), 134.


Lebedev:2014bba
O. Lebedev and Y. Mambrini,

Phys. Lett. B 734 (2014), 350-353.


Arcadi:2014lta
G. Arcadi, Y. Mambrini and F. Richard,

JCAP 03 (2015), 018








Cosme:2024ndc
C. Cosme, F. Costa and O. Lebedev,

JCAP 06 (2024), 031.






Cosme:2021baj
C. Cosme, M. Dutra, S. Godfrey and T. R. Gray,

JHEP 09 (2021), 056






















 

Gondolo:1990dk
P. Gondolo and G. Gelmini,

Nucl. Phys. B 360 (1991), 145-179






Koivunen:2024vhr
N. Koivunen, O. Lebedev and M. Raidal,

[arXiv:2403.15533 [hep-ph]].


Arcadi:2024wwg
G. Arcadi, F. Costa, A. Goudelis and O. Lebedev,

JHEP 07 (2024), 044.




 
Belanger:2018ccd
G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar,

Comput. Phys. Commun. 231, 173-186 (2018).


Alguero:2023zol
G. Alguero, G. Belanger, F. Boudjema, S. Chakraborti, A. Goudelis, S. Kraml, A. Mjallal and A. Pukhov,

Comput. Phys. Commun. 299, 109133 (2024).






Silva-Malpartida:2024emu
J. Silva-Malpartida, N. Bernal, J. Jones-Pérez and R. A. Lineros,

[arXiv:2408.08950 [hep-ph]].



















Cheung:2001wx
K. m. Cheung,

Phys. Lett. B 517 (2001), 167-176




 
 
 
ALEPH:2006bhb
J. Alcaraz et al. [ALEPH, DELPHI, L3, OPAL and LEP Electroweak Working Group],

[arXiv:hep-ex/0612034 [hep-ex]].



ALEPH:2013dgf
S. Schael et al. [ALEPH, DELPHI, L3, OPAL and LEP Electroweak],

Phys. Rept. 532 (2013), 119-244






 
Falkowski:2015krw
A. Falkowski and K. Mimouni,

JHEP 02 (2016), 086









ATLAS:2019erb
G. Aad et al. [ATLAS],

Phys. Lett. B 796 (2019), 68-87






Alwall:2014hca
J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro,

JHEP 07 (2014), 079.



Bandyopadhyay:2018cwu
T. Bandyopadhyay, G. Bhattacharyya, D. Das and A. Raychaudhuri,

Phys. Rev. D 98 (2018) no.3, 035027









Arcadi:2024ukq
G. Arcadi, D. Cabo-Almeida, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, J. P. Neto, M. Pierre, S. Profumo and F. S. Queiroz,

[arXiv:2403.15860 [hep-ph]].





HERMES:2006jyl
A. Airapetian et al. [HERMES],

Phys. Rev. D 75 (2007), 012007




 
 
 
 

LZ:2022lsv
J. Aalbers et al. [LZ],

Phys. Rev. Lett. 131 (2023) no.4, 041002






XENON:2018voc
E. Aprile et al. [XENON],

Phys. Rev. Lett. 121, no.11, 111302 (2018).




Billard:2021uyg
J. Billard, M. Boulay, S. Cebrián, L. Covi, G. Fiorillo, A. Green, J. Kopp, B. Majorovits, K. Palladino and F. Petricca, et al.

Rept. Prog. Phys. 85, no.5, 056201 (2022).



 
XENON:2020kmp
E. Aprile et al. [XENON],

JCAP 11, 031 (2020).
 
   
 
DARWIN:2016hyl
J. Aalbers et al. [DARWIN],

JCAP 11, 017 (2016).


 Hambye:2018dpi
T. Hambye, M. H. G. Tytgat, J. Vandecasteele and L. Vanderheyden,

Phys. Rev. D 98, no.7, 075017 (2018).



Boddy:2024vgt
K. K. Boddy, K. Freese, G. Montefalcone and B. Shams Es Haghi,

[arXiv:2405.06226 [hep-ph]].



Wong:2023qon
X. R. Wong and K. P. Xie,

Phys. Rev. D 108 (2023) no.5, 055035










 


McDaniel:2023bju
A. McDaniel, M. Ajello, C. M. Karwin, M. Di Mauro, A. Drlica-Wagner and M. A. Sánchez-Conde,

Phys. Rev. D 109 (2024) no.6, 063024







 
 
 


Berlin:2014tja
A. Berlin, D. Hooper and S. D. McDermott,

Phys. Rev. D 89 (2014) no.11, 115022




 

Klasen:2016qux
M. Klasen, F. Lyonnet and F. S. Queiroz,

Eur. Phys. J. C 77 (2017) no.5, 348




 




