Dark Matter Freeze-in from a $Z^\prime$ Reheaton

Abstract

We consider the Standard Model (SM) extended by a secluded $U(1)_D$ gauge sector encompassing a Dirac fermion ($χ$) dark matter (DM), an abelian gauge boson $Z^\prime$ and a SM-singlet complex-scalar field $Φ$, whose radial component drives cosmic inflation. When the Higgs portal coupling is small, the $Z^\prime$ then acts as a {\it ``reheaton''}, dominating the energy budget of the Universe before finally yielding the SM bath, with reheating temperature $< O(10)$ TeV, through the gauge portal interaction. We explore the possibility that DM freezes-in via non-thermal $Z^\prime$ decays before reheating ends, giving rise to substantial viable parameter space. We account for non-perturbative effects, relevant during the initial stages of reheating, using lattice simulations. We additionally show how the cosmological gravitational wave (GW) background produced by preheating and inflation allow for a direct probe of the reheating mechanism.

Figures (8)

• Figure 1: Left: we summarise the inflationary predictions for the generalised model of Higgs inflation in the $r$-$n_s$ plane. We use dashed black lines to illustrate the dependence on the e-folds of inflationary expansion ($N$) before the Hubble exit of the pivot scale ($k_*$) and the effective non-minimal coupling of the inflaton ($\xi$). We draw some contours in red to exhibit the dependence on the e-folds of primordial matter-domination, which is generally non-zero in our setup. The light and dark shaded regions correspond to 95% and 68% confidence intervals at the pivot scale $k_* = 0.05\ \text{MPc}^{-1}$; the model is in excellent agreement with the combined Planck and BK18 fit Planck:2018jriBICEP:2021xfz, as well as SPT SPT-3G:2025bzu and ACT ACT:2025fju, but there is a small tension with the inclusion of DESI DESI:2025zgx data by Ref. ACT:2025fju. Right: the normalisation of the scalar power spectrum is achieved by the parameter hierarchy $\lambda \ll \xi$, with the exact functional dependence plotted numerically above. As in Figure \ref{['fig:rnsplot']}, with contours that show a negligible dependence on $N_{m.d.}$.
• Figure 2: Power spectra at the end of inflation for the rescaled Goldstone modes in Coulomb gauge $\delta\tilde{\phi}_{I\mathbf{k}}$ (left) and the corresponding longitudinal polarisation modes of the $Z'$ in unitary gauge (right). The vectorial $\propto k^2$ contribution to the energy density per logarithmic interval is manifest in both cases, resulting in dramatic isocurvature suppression in the long wavelength limit Graham:2015rvaOzsoy:2023gnl.
• Figure 3: Left: the occupation number $n_k = \rho_k/\omega_k$ spectra for the real and imaginary components of the $\Phi$ field, with bluer contours for $t \sim t_{\text{end}}$ becoming redder at later times. See text for details. Right upper: the oscillating field zero mode compared with the root-mean-squared of the field fluctuations, which are amplified by preheating, restoring the potential minimum to the origin. Right lower: we compare estimates for the energy fraction in the background ($\rho_{bg}$) as well as fluctuations in the inflatons and Goldstones using $\delta\rho_X \sim(\dot{\delta\phi_X})^2$.
• Figure 4: Upper left: the kinetic ($\rho_K$), gradient ($\rho_G$) and potential ($\rho_V$) energy fractions (of the total $\rho$) are plotted over time (transparent), along with their time-averages (opaque), which approach constant values by the end of the simulation. (Note that the residual $\rho_V$ is likely subject to finite-size effects.) Lower left: the time-averaged equation of state is consistent with radiation ($w \simeq \frac{1}{3}$), even after fragmentation. Right: we plot the time evolution of the equation of state, from a simulation with $v_D = 3\times 10^{16}$ GeV ($k_{IR} = 0.08 \sqrt{\lambda_\Phi}\varphi_{\text{end}}$, $N= 128$, $dt = 10^{-3}\sqrt{\lambda_\Phi}\varphi_{\text{end}}$), during and after the non-thermal phase transition, together with a natural interpolation in yellow (leading to an extrapolation in blue) explained in the text. We rescale with a fiducial scale factor in order to use these fits for smaller $v_D$ consistent with (\ref{['eq:paramspace']}).
• Figure 5: Here we present the region of parameter space spanned by $\epsilon-m_{Z^\prime,0}$ where our proposed cosmology encompassing a $Z^\prime$ reheaton can be achieved (in white). In the green shaded region at the top, the $Z^\prime$ thermalise with the SM plasma, so that freeze-in is instead dominantly thermal. The brown shaded region is properly described by chiral perturbation theory (beyond the scope of this work). Several existing constraints from SN1987A Chang:2016ntp, E137 Marsicano:2018krp and BBN Berger:2016vxi are shown by the blue, red and light-blue shaded regions, respectively. The projected sensitivity reach of SHiP SHiP:2020vbd is also shown by the dark red dashed line. The light yellow shaded vertical region on the right indicates the parameter region where the near-global limit of $U(1)_D$ is invalid. Within the allowed region the purple contours represent the contours along which DM freeze-in relic density is satisfied for different values of the DM $U(1)_D$ charge $Q_\chi$, the gray dash-dotted lines denote the contours of different reheating temperatures, while the gray dotted lines represent the contours of different values of $N_{md}$. This plot is obtained for representative values of $v_D = 10^9\,{\rm GeV}$ and $m_\chi = 10\,{\rm keV}$, while the yellow star marks a suitably chosen benchmark point used in our discussions. See the text for details.