Revisiting sterile neutrino dark matter in gauged $U(1)_{B-L}$ model

Shintaro Eijima [email protected] ICRR, The University of Tokyo, Kashiwa, Chiba 277-8582, Japan Osamu Seto [email protected] Department of Physics, Hokkaido University, Sapporo 060-0810, Japan Takashi Shimomura [email protected] Faculty of Education, Miyazaki University, Miyazaki 889-2192, Japan Department of Physics, Kyushu University, Fukuoka 819-0395, Japan

Abstract

We reexamine sterile neutrino dark matter in gauged $U(1)_{B-L}$ model. Improvements have been made by tracing and careful evaluation of the evolution of the number densities of sterile neutrinos $N$ and extra neutral gauge bosons $Z^{\prime}$ . As a result, the cosmologically-interesting gauge coupling of $U(1)_{B-L}$ for freeze-in sterile neutrinos turns out to be smaller than the values reported in the literature. This avoids the overproduction of $Z^{\prime}$ so that it is consistent with the big bang nucleosynthesis and the cosmic microwave background constraints on the effective number of neutrino species. Similarly, the free-streaming length constraints exclude a large parameter space derived in previous studies. In addition to known freeze-in pair production of $N$ from the standard model fermion pairs, we find the case that $N$ is dominantly produced from a pair of $Z^{\prime}$ at the temperature characterized by the $B-L$ breaking scalar mass. Thus, the naive truncation of the $U(1)_{B-L}$ scalar contribution made in the literature is not valid.

^†^†preprint: UME-PP-021^†^†preprint: EPHOU-22-013^†^†preprint: KYUSHU-HET-244

I Introduction

Not only the existence of nonbaryonic dark matter (DM) but also nonvanishing neutrino masses are major open questions in particle physics and cosmology. Those are also hints and evidence for new physics beyond the Standard Model (SM). One of the simple extensions of the SM to explain those two problems is to introduce right-handed (RH) neutrinos, which are singlets under the SM gauge group. With Majorana masses of RH neutrinos, the tiny neutrino masses can be naturally generated through the seesaw mechanism Minkowski (1977); Yanagida (1979); Gell-Mann et al. (1979); Mohapatra and Senjanovic (1980). The resultant heavier-mass eigenstates compared with the active neutrinos responsible for neutrino oscillations are called sterile neutrinos. They slightly mix with the left-handed (LH) components via the active-sterile mixings. The lightest of them is a good candidate for (warm) DM Dolgov and Hansen (2002) if the lifetime is long enough Drewes et al. (2017); Boyarsky et al. (2019). The $\nu$ MSM Asaka et al. (2005); Asaka and Shaposhnikov (2005) is known as an extension with three RH neutrinos including such a sterile neutrino dark matter.

Three generations of RH neutrinos can be theoretically verified, once the SM is extended by introducing an extra $U(1)$ gauge symmetry under which RH neutrinos are charged. Since RH neutrinos are chiral, their charges and generations are determined by gauge anomaly cancellation. A well-studied example is gauged $U(1)_{B-L}$ symmetry Davidson (1979); Mohapatra and Marshak (1980); Marshak and Mohapatra (1980).

If RH neutrinos interact through an extra $U(1)$ gauge interaction, the sterile neutrino DM can be generated effectively. The Dodelson-Widrow (DW) mechanism Dodelson and Widrow (1994), where sterile neutrino DM is produced through the active-sterile mixings, conflicts with the observations in x-ray Boyarsky et al. (2006a, b); Boyarsky et al. (2007, 2008); Yuksel et al. (2008) and for Lyman- $\alpha$ forest (See, e.g., Boyarsky et al. (2019) and references therein). Thus, the production from scatterings through new mediators Khalil and Seto (2008); Kaneta et al. (2017); Biswas and Gupta (2016); Seto and Shimomura (2020); De Romeri et al. (2020); Lucente (2021); Bélanger et al. (2021) is an alternative promising scenario. This class of nonthermal production is sometime called “freeze-in” Hall et al. (2010). For a review, see, e.g., Refs. Baer et al. (2015); Shakya (2016). Nonthermal productions of an extra $U(1)$ gauge interacting sterile neutrino DM have been investigated for a large gauge coupling Khalil and Seto (2008), and for a small gauge coupling in Refs. Kaneta et al. (2017); Biswas and Gupta (2016); Fileviez Pérez et al. (2019); Heeba and Kahlhoefer (2020); Okada et al. (2020, 2021); Iwamoto et al. (2022).

The purpose of this paper is to reexamine the freeze-in production of $U(1)_{B-L}$ gauge-interacting sterile neutrino DM. After we examine the production processes of the sterile neutrino $N$ and the extra neutral gauge boson $Z^{\prime}$ in detail and trace the evolution of the number densities of $N$ and $Z^{\prime}$ , we evaluate cosmological constraints on the model from the viewpoint of light degrees of freedom and the structure formation. We find that the inverse decay production of $Z^{\prime}$ from SM particles in equilibrium is dominant in cases with a very small $U(1)_{B-L}$ gauge coupling. Due to the contribution, the resultant abundances of the sterile neutrino DM and light $Z^{\prime}$ as dark radiation (DR) are larger than previously estimated. Therefore, the impacts of this increase on the big bang nucleosynthesis (BBN) and the cosmic microwave background (CMB) are estimated, in which constraints on the gauge coupling are obtained. We also find that the contribution from the on shell $Z^{\prime}$ decay is dominant in the scattering production process through $s$ -channel $Z^{\prime}$ exchange. Since the produced $N$ from the $Z^{\prime}$ decay is relativistic, the free-streaming length constraints limit the viable parameter space. Moreover, we find cases where the sterile neutrino DM is dominantly produced by the pair production from a pair of $Z^{\prime}$ . The cross section of this mode depends on the mass of $U(1)_{B-L}$ breaking scalar. Thus, the naive truncation of the scalar contribution made in literature cannot be validated.

This paper is organized as follows. We describe Langrangian of the model, and read vertices of interactions in Sec. II. After we summarize Boltzmann equations which need to be solved in Sec. III, we consider three production scenarios of the sterile neutrino DM and estimate the DM abundance by taking into account cosmological constraints in Sec. IV. We shortly discuss the implication to dark matter detection in Sec. V. We discuss the conclusion in Sec. VI.

II Model

We consider the extension of the SM by gauging $U(1)_{B-L}$ symmetry, where $B$ and $L$ are the baryon and lepton number, respectively. Under the $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times U(1)_{B-L}$ gauge groups, three generations of RH neutrinos ( $\nu_{R}^{i}$ with $i$ running over $1,2,3$ ) have to be introduced for the anomaly cancellation. The scalar sector of the SM is also extended by introducing one complex scalar $\Phi_{B-L}$ , which is charged under the $U(1)_{B-L}$ , to break the extra gauge symmetry spontaneously. This is the minimal extension regarding the gauged $U(1)_{B-L}$ symmetry. In Table 1, $Q$ $(u_{R},~{}d_{R})$ and $L$ $(e_{R},~{}\nu_{R})$ denote the LH (RH) quarks and leptons, respectively.

	$SU(3)_{C}$	$SU(2)_{L}$	$U(1)_{Y}$	$U(1)_{B-L}$
$Q^{i}$	$\mathbf{3}$	$\mathbf{2}$	$\frac{1}{6}$	$\frac{1}{3}$
$u_{R}^{i}$	$\mathbf{3}$	$\mathbf{1}$	$\frac{2}{3}$	$\frac{1}{3}$
$d_{R}^{i}$	$\mathbf{3}$	$\mathbf{1}$	$-\frac{1}{3}$	$\frac{1}{3}$
$L^{i}$	$\mathbf{1}$	$\mathbf{2}$	$-\frac{1}{2}$	$-1$
$e_{R}^{i}$	$\mathbf{1}$	$\mathbf{1}$	$-1$	$-1$
$\nu_{R}^{i}$	$\mathbf{1}$	$\mathbf{1}$	$0$	$-1$
$\Phi_{H}$	$\mathbf{1}$	$\mathbf{2}$	$\frac{1}{2}$	$0$
$\Phi_{B-L}$	$\mathbf{1}$	$\mathbf{1}$	$0$	$2$

Table 1: In addition to the SM particle content, three RH neutrinos

\nu_{R}^{i}

(

i=1,2,3

) and one

U(1)_{B-L}

Higgs field

\Phi_{B-L}

are introduced.

II.1 Lagrangian

The Lagrangian is given by

	$\displaystyle\mathcal{L}$	$\displaystyle=\mathcal{L}_{\text{SM}}+\mathcal{L}_{\nu_{R}}+V(\Phi_{H},\Phi_{B-L}),$		(1)
	$\displaystyle\mathcal{L}_{\nu_{R}}$	$\displaystyle=\overline{\nu_{R}^{i}}(iD_{\mu}\gamma^{\mu})\nu_{R}^{i}-y_{\nu}{}_{ij}\overline{L^{i}}\tilde{\Phi}_{H}\nu_{R}^{j}-\frac{1}{2}y_{\nu_{R}^{i}}\Phi_{B-L}\overline{\nu_{R}^{i~{}C}}\nu_{R}^{i}+\mathrm{H.c.},$		(2)

where $\mathcal{L}_{\text{SM}}$ , $\mathcal{L}_{\nu_{R}}$ , and $V(\Phi_{H},\Phi_{B-L})$ are Lagrangian of the SM, RH neutrinos, and the scalar potential of this model, respectively. In Eq. (2), the superscript $C$ denotes the charge conjugation of $\nu^{i}_{R}$ , and $\tilde{\Phi}_{H}\equiv\epsilon\Phi_{H}^{\dagger}$ is the conjugation of the Higgs field $\Phi$ with $\epsilon$ being the antisymmetric tensor. Yukawa couplings with LH lepton doublets and among RH neutrinos are denoted by $y_{\nu}$ and $y_{\nu_{R}^{i}}$ respectively, in which $i$ and $j$ are indices of flavor or generation. We work on the diagonal basis of $y_{\nu_{R}^{i}}$ without loss of generality.

Due to the $U(1)_{B-L}$ symmetry, the covariant derivative is modified as

\displaystyle D_{\mu}=\partial_{\mu}-ig_{2}W_{\mu}-ig_{1}YB_{\mu}-ig_{B-L}q_{B-L}X_{\mu},

(3)

where $W,~{}B$ , and $X$ represent the gauge fields, and $g_{2},~{}g_{1}$ , and $g_{B-L}$ are the gauge coupling constants of $SU(2)_{L},~{}U(1)_{Y}$ , and $U(1)_{B-L}$ , respectively. The $U(1)_{Y}$ and $U(1)_{B-L}$ charges listed in Table 1 are denoted as $Y$ and $q_{B-L}$ . We omit any symbol about the $SU(3)_{C}$ color interaction.

There may exist a gauge kinetic mixing term

\displaystyle\mathcal{L}_{\varepsilon}=\frac{\varepsilon}{2}B_{\mu\nu}X^{\mu\nu},

(4)

where $B_{\mu\nu}$ and $X_{\mu\nu}$ are the gauge field strength of $U(1)_{Y}$ and $U(1)_{B-L}$ gauge field respectively, and $\varepsilon$ is a mixing parameter. The importance of this term has been studied intensively in the context of the so-called dark photon. Since we are interested in the $U(1)_{B-L}$ gauge interaction on RH neutrinos, we concentrate on the cases where the effect of the gauge kinetic mixing is negligible and set $\varepsilon$ vanishing in this paper.¹¹1If the gauge kinetic mixing effects are more dominant than the direct gauge interaction, the model would reduce to a dark photon model. For a recent review on the dark photon, see, e.g., Refs. Caputo et al. (2021).

II.2 Scalar potential

II.2.1 Mass eigenstates and Higgs mixing

We derive the masses of introduced particles with the following scalar potential,

	$\displaystyle V(\Phi_{H},\Phi_{B-L})=$	$\displaystyle\frac{1}{2}\lambda_{1}\left(\|\Phi_{H}\|^{2}-\frac{v^{2}}{2}\right)^{2}+\frac{1}{2}\lambda_{2}\left(\|\Phi_{B-L}\|^{2}-\frac{v_{B-L}^{2}}{2}\right)^{2}$
		$\displaystyle+\lambda_{3}\left(\|\Phi_{H}\|^{2}-\frac{v^{2}}{2}\right)\left(\|\Phi_{B-L}\|^{2}-\frac{v_{B-L}^{2}}{2}\right),$		(5)

where all parameters, $\lambda_{1},\lambda_{2},\lambda_{3},v_{B-L}$ , and $v\simeq 246$ GeV, are real and positive. At the electroweak and $B-L$ broken vacuum, $\Phi_{H}$ and $\Phi_{B-L}$ fields can be expanded around those vacuum expectation values (VEVs) $v$ and $v_{B-L}$ , respectively.

In this vacuum, the $U(1)_{B-L}$ gauge boson $X$ absorbs the corresponding Nambu-Goldstone mode and becomes the massive extra neutral gauge boson $Z^{\prime}$ with the mass

\displaystyle m_{Z^{\prime}}^{2}=4g_{B-L}^{2}v_{B-L}^{2},

(6)

and three RH neutrinos also obtain their Majorana masses

\displaystyle m_{\nu_{R}^{i}}=\frac{y_{\nu_{R}^{i}}}{\sqrt{2}}v_{B-L}.

(7)

Then, the nonvanishing neutrino masses can be generated through the type-I seesaw mechanism as

\displaystyle m_{\nu}{}_{ij}

\displaystyle=-m_{D}{}_{ik}\frac{1}{m_{\nu_{R}^{k}}}m_{D}^{T}{}_{kj},

(8)

with

\displaystyle m_{D}{}_{ik}

\displaystyle=\frac{y_{\nu ik}v}{\sqrt{2}}.

(9)

The mass eigenstates of the scalars, $h$ and $\phi$ , are obtained from those physical fluctuations $\phi_{H}$ and $\phi_{B-L}$ in $\Phi_{H}$ and $\Phi_{B-L}$ , respectively, as

\displaystyle\left(\begin{array}[]{c}\phi_{H}\\ \phi_{B-L}\\ \end{array}\right)=\left(\begin{array}[]{cc}\cos\alpha&\sin\alpha\\ -\sin\alpha&\cos\alpha\\ \end{array}\right)\left(\begin{array}[]{c}h\\ \phi\\ \end{array}\right),

(16)

with the mixing angle $\alpha$ . Their masses are given by


	$\displaystyle m_{h}^{2}=\frac{1}{2}\left(\lambda_{1}v^{2}+\lambda_{2}v_{B-L}^{2}+\frac{\lambda_{1}v^{2}-\lambda_{2}v_{B-L}^{2}}{\cos(2\alpha)}\right),$		(17a)
	$\displaystyle m_{\phi}^{2}=\frac{1}{2}\left(\lambda_{1}v^{2}+\lambda_{2}v_{B-L}^{2}-\frac{\lambda_{1}v^{2}-\lambda_{2}v_{B-L}^{2}}{\cos(2\alpha)}\right).$		(17b)

At the $\alpha\rightarrow 0$ limit, $h$ is reduced to the SM Higgs boson $\phi_{H}$ . For a small $\alpha\ll 1$ , $h$ and $\phi$ are identified with the SM-like Higgs boson and the singlet-like scalar, respectively. Thus, we take $m_{h}\simeq 125$ GeV. The mixing angle $\alpha$ can be expressed in terms of $\lambda_{3}$ as

\displaystyle\sin(2\alpha)

\displaystyle\simeq\frac{2vv_{B-L}}{m_{\phi}^{2}-m_{h}^{2}}\lambda_{3}=\frac{vm_{Z^{\prime}}}{m_{\phi}^{2}-m_{h}^{2}}\frac{\lambda_{3}}{g_{B-L}},

(18)

where we have used Eq. (6) in the last equality.

In terms of $h$ and $\phi$ , the scalar potential (5) is rewritten as

\displaystyle V=

\displaystyle\frac{1}{2}m_{h}^{2}h^{2}+\frac{1}{2}m_{\phi}^{2}\phi^{2}+\frac{1}{6}C_{hhh}h^{3}+\frac{1}{2}C_{hh\phi}h^{2}\phi+\frac{1}{2}C_{h\phi\phi}h\phi^{2}+\frac{1}{6}C_{\phi\phi\phi}\phi^{3}+\cdots,

(19)

where the ellipsis denotes quartic terms that are irrelevant for our following analysis. The scalar trilinear couplings are given by


$\displaystyle C_{hhh}$	$\displaystyle=3\frac{m_{h}^{2}\left(v_{B-L}\cos^{3}\alpha-v\sin^{3}\alpha\right)}{vv_{B-L}},$	(20a)
$\displaystyle C_{hh\phi}$	$\displaystyle=\frac{\sin(2\alpha)(2m_{h}^{2}+m_{\phi}^{2})(v\sin\alpha+v_{B-L}\cos\alpha)}{2vv_{B-L}},$	(20b)
$\displaystyle C_{h\phi\phi}$	$\displaystyle=\frac{\sin(2\alpha)(m_{h}^{2}+2m_{\phi}^{2})(v_{B-L}\sin\alpha-v\cos\alpha)}{2vv_{B-L}},$	(20c)
$\displaystyle C_{\phi\phi\phi}$	$\displaystyle=3\frac{m_{\phi}^{2}\left(v\cos^{3}\alpha+v_{B-L}\sin^{3}\alpha\right)}{vv_{B-L}},$	(20d)

where $C_{hh\phi}$ and $C_{h\phi\phi}$ are suppressed for the small mixing angle $\alpha$ .

Similarly, each Yukawa coupling of the SM fermions $f$ and RH neutrinos $\nu_{R}$ with $h$ and $\phi$ is suppressed due to the Higgs mixing with the following factor


	$\displaystyle C_{hf}=\cos\alpha,$		(21a)
	$\displaystyle C_{\phi f}=\sin\alpha,$		(21b)
	$\displaystyle C_{h\nu_{R}}=-\sin\alpha,$		(21c)
	$\displaystyle C_{\phi\nu_{R}}=\cos\alpha,$		(21d)

and the mixing suppression factors to the gauge couplings of the SM gauge bosons $V(W,Z,A)$ and $Z^{\prime}$ to $h$ and $\phi$ are given by


	$\displaystyle C_{hV}=\cos\alpha,$		(22a)
	$\displaystyle C_{\phi V}=\sin\alpha,$		(22b)
	$\displaystyle C_{hZ^{\prime}}=-\sin\alpha,$		(22c)
	$\displaystyle C_{\phi Z^{\prime}}=\cos\alpha.$		(22d)

From the above equations, we can read the interactions of $Z^{\prime}$ and $\nu_{R}^{i}$ are affected by the Higgs mixing. We will see later that the mixing is essential for the production of the sterile neutrino DM in some cases.

II.2.2 Decays of scalars and range of Higgs mixing

Due to the Higgs mixing, decay rates of the SM-like Higgs boson into the SM particles are multiplied by $\cos^{2}\alpha$ in our model. In addition, if the mass of the singletlike scalar $\phi$ is smaller than one half of the SM-like Higgs boson mass, another decay mode $h\rightarrow\phi\phi$ is possible. Thus, the partial decay rates of $h$ are given by


$\displaystyle\Gamma_{h}(h\rightarrow\mathrm{SM})$	$\displaystyle=\cos^{2}\alpha\Gamma_{h_{\mathrm{SM}}}(h_{\mathrm{SM}}\rightarrow\mathrm{SM}),$	(23a)
$\displaystyle\Gamma_{h}(h\rightarrow\phi\phi)$	$\displaystyle=\frac{\sqrt{m_{h}^{2}-4m_{\phi}^{2}}}{16\pi m_{h}^{2}}\|C_{h\phi\phi}\|^{2},$	(23b)
$\displaystyle\Gamma_{h}(h\rightarrow Z^{\prime}Z^{\prime})$	$\displaystyle=\frac{\sqrt{m_{h}^{2}-4m_{Z^{\prime}}^{2}}}{32\pi m_{h}^{2}}\sin^{2}\alpha\frac{m_{h}^{4}-4m_{h}^{2}m_{Z^{\prime}}^{2}+12m_{Z^{\prime}}^{4}}{v_{B-L}^{2}},$	(23c)
$\displaystyle\Gamma_{h}(h\rightarrow NN)$	$\displaystyle=\sum_{i}\frac{1}{16\pi m_{h}}\left(1-\frac{4m_{N_{i}}^{2}}{m_{h}^{2}}\right)^{3/2}\sin^{2}\alpha\frac{m_{h}^{2}m_{N_{i}}^{2}}{v_{B-L}^{2}},$	(23d)

where $h_{\mathrm{SM}}\rightarrow\mathrm{SM}$ stands for the decay processes of SM Higgs boson into all final states in the SM model and its decay rate is $\Gamma_{h_{\mathrm{SM}}}(h_{\mathrm{SM}}\rightarrow\mathrm{SM})\simeq 4$ MeV. Here, $m_{N_{i}}$ is the mass of sterile neutrinos and $m_{N_{i}}\simeq m_{\nu_{R}^{i}}$ for small active-sterile mixings.

The total decay rate $\Gamma_{h}$ is given by the sum of Eqs. (23). For $m_{\phi}<m_{h}/2$ , the current constraints on exotic decay of the SM Higgs boson as $\text{Br}(h\rightarrow\text{invisible})\lesssim 19~{}\%$ Zyla et al. (2020) can be recast as $\sin\alpha\lesssim 0.2$ . For $m_{\phi}>m_{h}/2$ , the obtained range $\alpha<\mathcal{O}(0.1)$ is also consistent with measurements in the LHC Aad et al. (2020).

In the mass spectra of our interest, the singletlike scalar $\phi$ decays dominantly into pairs of the SM fermions through the Higgs mixing with the rate

\Gamma_{\phi}(\phi\rightarrow\mathrm{SM})=\sin^{2}\alpha\Gamma_{h_{\mathrm{SM}}}(\phi\rightarrow\mathrm{SM}),

(24)

where $\Gamma_{h_{\mathrm{SM}}}(\phi\rightarrow\mathrm{SM})$ expresses the decay rate of $\phi$ with the same SM interactions of $h$ , and the scale of running parameters is taken at $m_{\phi}$ in the calculation. The following decay modes with partial decay rates


$\displaystyle\Gamma_{\phi}(\phi\rightarrow hh)$	$\displaystyle=\frac{\sqrt{m_{\phi}^{2}-4m_{h}^{2}}}{16\pi m_{\phi}^{2}}\|C_{hh\phi}\|^{2},$	(25a)
$\displaystyle\Gamma_{\phi}(\phi\rightarrow Z^{\prime}Z^{\prime})$	$\displaystyle=\frac{\sqrt{m_{\phi}^{2}-4m_{Z^{\prime}}^{2}}}{16\pi m_{\phi}^{2}}\cos^{2}\alpha\frac{m_{\phi}^{4}-4m_{\phi}^{2}m_{Z^{\prime}}^{2}+12m_{Z^{\prime}}^{4}}{v_{B-L}^{2}},$	(25b)
$\displaystyle\Gamma_{\phi}(\phi\rightarrow NN)$	$\displaystyle=\sum_{i}\frac{1}{16\pi m_{\phi}}\left(1-\frac{4m_{N_{i}}^{2}}{m_{\phi}^{2}}\right)^{3/2}\cos^{2}\alpha\frac{m_{\phi}^{2}m_{N_{i}}^{2}}{v_{B-L}^{2}},$	(25c)

can also open depending on the mass spectrum. These, however, are negligible compared with Eq. (24); $\Gamma_{\phi}(\phi\rightarrow hh)$ is due to $v_{B-L}\gg v$ , and the others are also with the suppression by $v_{B-L}$ unless we take $\alpha\rightarrow 0$ . Therefore, we find typically $\Gamma_{\phi}\sim\sin^{2}\alpha$ MeV. Bounds on the Higgs mixing between a light scalar and the SM-like Higgs boson have been derived from the LEP experiments similarly Barate et al. (2003), and the range $\alpha<\mathcal{O}(0.1)$ coincides with the bounds as well.

If the singletlike scalar is lighter than about a few GeV, the constraints from meson decays by the LHCb Aaij et al. (2015, 2017) and CHARM Bergsma et al. (1985) are more stringent than ATLAS, CMS, and the LEP. In such a mass range, the bound $10^{-5}\lesssim\sin\alpha\lesssim 10^{-4}$ has been obtained Winkler (2019) where the lower bound is set by demanding that the lifetime of $\phi$ must be shorter than $\mathcal{O}(0.1)$ seconds so as not to affect the BBN Fradette and Pospelov (2017).

III The Boltzmann equation

In this section we describe our Boltzmann equations to calculate the evolution and abundance of the sterile neutrino DM and $Z^{\prime}$ via freeze-in production. Here and hereafter the DM is denoted with $N$ , which is a suitable one among the three sterile neutrinos. The Boltzmann equations for the number density of $N$ and $Z^{\prime}$ are given by


$\displaystyle\frac{dn_{N}}{dt}+3Hn_{N}=$	$\displaystyle\sum_{i,j}\langle\sigma v(ij\rightarrow NN)\rangle(n_{i}n_{j}-n_{N}^{2})+\sum_{i}\langle\Gamma(i\rightarrow NN)\rangle n_{i},$	(26a)
$\displaystyle\frac{dn_{Z^{\prime}}}{dt}+3Hn_{Z^{\prime}}=$	$\displaystyle\sum_{i,j}\langle\sigma v(ij\rightarrow Z^{\prime}Z^{\prime})\rangle(n_{i}n_{j}-n_{Z^{\prime}}^{2})+\sum_{i}\langle\Gamma(i\rightarrow Z^{\prime}Z^{\prime})\rangle n_{i}$
	$\displaystyle+\sum_{i,j,k}\langle\sigma v(Z^{\prime}i\rightarrow jk)n_{i}\rangle\left(n_{Z^{\prime}}-n^{\mathrm{eq}}_{Z^{\prime}}\right)-\sum_{i,j}\langle\Gamma(Z^{\prime}\rightarrow ij)\rangle(n_{Z^{\prime}}-n^{\mathrm{eq}}_{Z^{\prime}}),$	(26b)

where $i,j$ , and $k$ are possible initial and final states in reactions, and $n_{i}$ is the number density of $i$ -particle. Note that $N$ and $Z^{\prime}$ productions from the decays of intermediate particles produced on-pole in the first term are treated properly to avoid double-counting in Eqs. (26). The cosmic expansion rate $H$ in the radiation-dominated (RD) universe is given by

	$\displaystyle H^{2}$	$\displaystyle\equiv\left(\frac{\dot{a}}{a}\right)^{2}=\frac{1}{3M_{P}^{2}}\rho_{r},$		(27)
	$\displaystyle\rho_{r}$	$\displaystyle=\frac{\pi^{2}g_{*}}{30}T^{4},$		(28)

where $a$ is the scale factor of the universe and dot denotes derivative with respect to the cosmic time $t$ . The reduced Planck mass is $M_{P}\simeq 2.4\times 10^{18}$ GeV, $\rho_{r}$ is the energy density of radiation with the temperature $T$ and $g_{*}$ denotes the number of relativistic degrees of freedom.

A thermally averaged product of the scattering cross section and the relative velocity in Eqs. (26a) and (26b) are given by Gondolo and Gelmini (1991)

\displaystyle\langle\sigma v\rangle n_{i}n_{j}

\displaystyle=\frac{T}{32\pi^{4}}\sum_{i,j}\int^{\infty}_{(m_{i}+m_{j})^{2}}dsg_{i}g_{j}p_{ij}4E_{i}E_{j}\sigma vK_{1}\left(\frac{\sqrt{s}}{T}\right),

(29)

and

$\displaystyle 4E_{i}E_{j}\sigma v$	$\displaystyle\equiv\prod_{f}\int\frac{d^{3}p_{f}}{(2\pi)^{3}}\frac{1}{2E_{f}}\overline{\|\mathcal{M}\|^{2}}(2\pi)^{4}\delta^{(4)}(p_{i}+p_{j}-\sum p_{f})$
	$\displaystyle=\frac{1}{16\pi}\frac{2\|q_{f}\|}{\sqrt{s}}\int\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$	(30)
$\displaystyle 2\|q_{f}\|$	$\displaystyle=\sqrt{s-4m_{N}^{2}},$	(31)
$\displaystyle p_{ij}$	$\displaystyle\equiv\frac{\sqrt{s-(m_{i}+m_{j})^{2}}\sqrt{s-(m_{i}-m_{j})^{2}}}{2\sqrt{s}},$	(32)

where $i$ is an initial state with mass $m_{i}$ , energy $E_{i}$ , and internal degrees of freedom $g_{i}$ . The center-of-mass energy squared is given by $s=(E_{i}+E_{j})^{2}$ and three-momentum of a final-state particle is denoted by $q_{f}$ . $K_{i}(z)$ is the modified Bessel function of the $i$ th kind.

We consider the vanishing limit of the active-sterile mixings unless otherwise stated in the following analyses. This is because our purpose is to investigate the production of $N$ dominantly through the $g_{B-L}$ gauge interaction. For the very small $g_{B-L}$ we are interested in, the first term in right-handed side of Eq. (26b) representing the pair productions $ij\rightarrow Z^{\prime}Z^{\prime}$ is actually negligible compared with the other terms. This is because the cross sections for such pair productions are suppressed by $g_{B-L}^{4}$ while processes described in other terms are suppressed by only $g_{B-L}^{2}$ .

The second terms in the right-handed side of Eqs. (26a) and (26b) represent the production by the decay of an $i$ -particle, and

\displaystyle\langle\Gamma(i\rightarrow NN\,\mathrm{or}\,Z^{\prime}Z^{\prime})\rangle=\frac{K_{1}\left(\frac{m_{i}}{T}\right)}{K_{2}\left(\frac{m_{i}}{T}\right)}\Gamma(i\rightarrow NN\,\mathrm{or}\,Z^{\prime}Z^{\prime}),

(33)

is the thermal averaged partial decay rate of the $i$ -particle which is suppressed for a high temperature, $T\gg m_{i}$ , by the time dilation.

The third term in right-handed side of Eq. (26b) principally denotes the processes of $f\bar{f}\leftrightarrow Z^{\prime}\gamma$ , $f\gamma\leftrightarrow fZ^{\prime}$ and $\bar{f}\gamma\leftrightarrow\bar{f}Z^{\prime}$ . The thermal averaging of $\sigma vn$ is defined as

\displaystyle\langle\sigma vn_{i}\rangle n_{Z^{\prime}}^{\mathrm{eq}}

\displaystyle\equiv\frac{T}{32\pi^{4}}\int_{(m_{i}+m_{Z^{\prime}})^{2}}^{\infty}dsg_{i}g_{Z^{\prime}}p_{i{Z^{\prime}}}(4E_{i}E_{Z^{\prime}}\sigma v)K_{1}\left(\frac{\sqrt{s}}{T}\right),

(34)

with

\displaystyle n_{Z^{\prime}}^{\mathrm{eq}}

\displaystyle=\frac{T}{2\pi^{2}}g_{Z^{\prime}}m_{Z^{\prime}}^{2}K_{2}\left(\frac{m_{Z^{\prime}}}{T}\right),

(35)

where the superscript “ $\mathrm{eq}$ ” stands for the equilibrium value. Those turn out to be actually negligible compared with the fourth term. This can be understood from the fact that the cross sections of those $\gamma-Z^{\prime}$ scatterings are suppressed by $g_{B-L}^{2}\alpha_{\mathrm{em}}$ with $\alpha_{\mathrm{em}}=e^{2}/(4\pi)$ , while the following fourth term is suppressed by only $g_{B-L}^{2}$ .

The fourth term in Eqs. (26b) denotes the decay and inverse decay of $Z^{\prime}$ . The extra neutral gauge boson $Z^{\prime}$ decays into all fermions charged under the $U(1)_{B-L}$ . The partial decay widths of $Z^{\prime}$ are given by


$\displaystyle\Gamma_{Z^{\prime}}(Z^{\prime}\rightarrow f\bar{f})$	$\displaystyle=\sum_{f}\frac{g_{B-L}^{2}q_{B-L}^{2}N_{c}}{12\pi m_{Z^{\prime}}^{2}}\left(2m_{f}^{2}+m_{Z^{\prime}}^{2}\right)\left(m_{Z^{\prime}}^{2}-4m_{f}^{2}\right)^{1/2},$	(36a)
$\displaystyle\Gamma_{Z^{\prime}}(Z^{\prime}\rightarrow NN)$	$\displaystyle=\sum_{i}\frac{g_{B-L}^{2}}{24\pi m_{Z^{\prime}}^{2}}\left(m_{Z^{\prime}}^{2}-4m_{N_{i}}^{2}\right)^{3/2},$	(36b)

where the number of color $N_{c}=3$ is for quark final states. In this paper, we consider the situation where $h$ and $\phi$ are enough heavier than $Z^{\prime}$ . Then, the total decay rate is given by the sum of those. The inverse decay is the most efficient process to thermalize $Z^{\prime}$ disregarded in previous studies. This is our new observation in this work. Because of this efficient thermalization of $Z^{\prime}$ , we find the new freeze-in scenario where the $Z^{\prime}Z^{\prime}\rightarrow NN$ is the dominant production mode for $m_{Z^{\prime}}\gtrsim 1$ MeV. On the other hand, for $m_{Z^{\prime}}\lesssim 1$ MeV, the magnitude of the $B-L$ gauge coupling $g_{B-L}$ turns out to be smaller than $\mathcal{O}(10^{-12})$ to avoid the overproduction of light $Z^{\prime}$ as dark radiation.

IV Abundance of sterile neutrino DM

In this section we present the parameter region where cosmological constraints are satisfied and the observed DM abundance is explained with $N$ . The following three mass spectra are considered. Dominant processes of the DM production depend on the spectra, and hence the obtained regions are different in the three spectra. Before we show the results in detail, we briefly summarize those as

A

Heavy gauge boson region $m_{Z^{\prime}}>2m_{N}$ : The DM $N$ and $Z^{\prime}$ are not thermalized in this case. The $Z^{\prime}$ are dominantly produced on shell by dominantly its inverse decay. The DM is produced by the subsequent nonthermal decay $Z^{\prime}\to NN$ . By taking the free-streaming length constraints into account for $m_{Z^{\prime}}\lesssim 100$ GeV, the allowed mass range of $N$ turns out to be $m_{N}\gtrsim 1$ MeV and $g_{B-L}>10^{-12}$ .
B

Light gauge boson region $2m_{N}>m_{Z^{\prime}}>1$ MeV : $Z^{\prime}$ is thermalized by its inverse decay. The DM is dominantly produced via pair annihilation $Z^{\prime}Z^{\prime}\to NN$ . In this case, the DM abundance depends on the mass of $\phi$ due to the $s$ -channel exchange of $\phi$ . Taking $m_{N}=1$ GeV and $1$ GeV $<m_{\phi}<100$ GeV, the allowed region is found in $10^{-10}<g_{B-L}<10^{-6}$ and $10^{-3}$ GeV $<m_{Z^{\prime}}<2$ GeV.
C

Very light gauge boson region $2m_{N}>1$ MeV $>m_{Z^{\prime}}$ : The $B-L$ gauge coupling must be $g_{B-L}\lesssim 10^{-12}$ to avoid the BBN and the CMB constraints. The DM must be dominantly produced from the scatterings of the SM particles and $\phi$ via $\phi/h$ $s$ -channel exchange.

IV.1 Heavy gauge boson region $m_{Z^{\prime}}>2m_{N}$

The freeze-in DM production by the mediation of the extra gauge boson can be effective for $m_{Z^{\prime}}>2m_{N}$ . Under this mass spectrum, the production of $N$ can be dominated by the nonthermal decay of ${Z^{\prime}}$ . Since not only $N$ but also ${Z^{\prime}}$ cannot be thermalized, we need to simultaneously solve the Boltzmann Eqs. (26a) and (26b), which are rewritten as


$\displaystyle\left(\frac{dx}{dt}\right)\frac{dY_{N}}{dx}=$	$\displaystyle\langle\Gamma(Z^{\prime}\rightarrow NN)\rangle Y_{Z^{\prime}},$	(37a)
$\displaystyle\left(\frac{dx}{dt}\right)\frac{dY_{Z^{\prime}}}{dx}=$	$\displaystyle\langle\Gamma(\phi\rightarrow Z^{\prime}Z^{\prime})\rangle Y_{\phi}-\sum\langle\sigma v(Z^{\prime}i\leftrightarrow jk)n_{i}\rangle(Y_{Z^{\prime}}-Y^{\mathrm{eq}}_{Z^{\prime}})$
	$\displaystyle-\sum\langle\Gamma(Z^{\prime}\leftrightarrow ij)\rangle(Y_{Z^{\prime}}-Y^{\mathrm{eq}}_{Z^{\prime}}),$	(37b)

where $x\equiv M/T$ with $M$ being a mass scale for the normalization is a dimensionless variable. The yield abundance $Y_{i}\equiv n_{i}/s$ is defined as the ratio of the number density to the entropy density

\displaystyle s=\frac{2\pi^{2}g_{*S}}{45}T^{3},

(38)

with $g_{*S}$ being the total relativistic degrees of freedom for the entropy. Here, the pair production of $N$ by scattering $ij\rightarrow NN$ is dominated by the resonant processes of $s$ -channel $Z^{\prime}$ mediation from $f\bar{f}$ initial states. Since we have included the inverse decay of $Z^{\prime}$ , $f\bar{f}\rightarrow Z^{\prime}$ , and the decay of $Z^{\prime}$ into $NN$ , we have discarded the term for $f\bar{f}\rightarrow NN$ to avoid the double counting.

Refer to caption — Figure 1: Evolution of the yields of ${Z^{\prime}}$ and $N$ for $m_{Z^{\prime}}=100$ GeV, $g_{B-L}=5\times 10^{-10}$ , $m_{N}=1$ MeV, $\alpha=0$ .

We show, in Fig. 1, the typical evolution of $Y_{Z^{\prime}}$ and $Y_{N}$ for $m_{Z^{\prime}}=100$ GeV, $g_{B-L}=5\times 10^{-10}$ , $m_{N}=1$ MeV, and $\alpha=0$ .²²2This condition $\alpha=0$ is taken to suppress scalar mediated processes and is not necessarily satisfied exactly. However, when $\alpha$ is sufficiently large, the scalar mediated processes easily dominate over the $Z^{\prime}$ mediated processes and it is reduced to the Higgs portal freeze-in Majorana DM model. The orange curve is the thermal equilibrium yield value of the $Z^{\prime}$ boson. The blue and green solid curves represent values of $Y_{Z^{\prime}}$ and $Y_{N}$ , respectively. Here, we have confirmed that the $\gamma-Z^{\prime}$ scatterings are negligible compared to the inverse decay of $Z^{\prime}$ , as mentioned above. Once $Y_{N}(x\rightarrow\infty)$ is obtained, the present relic density is evaluated as

\displaystyle\Omega h^{2}=\frac{m_{N}}{\rho_{\mathrm{crit}}/s_{0}}Y_{N}

(39)

where $(\rho_{\mathrm{crit}}/s_{0})^{-1}=2.8\times 10^{8}/\mathrm{GeV}$ is given by the present entropy density $s_{0}$ and the critical density $\rho_{\mathrm{crit}}=3M_{P}^{2}H_{0}^{2}$ with $H_{0}$ being the present Hubble parameter.

Since the DM $N$ are produced by the decay of out-of-equilibrium $Z^{\prime}$ , under the mass spectrum of $m_{Z^{\prime}}\gg m_{N}$ , $N$ could be generated with a large momentum of about a half of the $Z^{\prime}$ mass. Such an energetic $N$ can have a large free-streaming length and erase small scale structures. The resultant comoving free-streaming scale can be calculated as Kolb and Turner (1990)

	$\displaystyle R_{f}$	$\displaystyle=\int^{t_{\mathrm{mre}}}_{t_{\mathrm{dec}}}\frac{v(t^{\prime})}{a(t^{\prime})}dt^{\prime}$
		$\displaystyle\simeq\frac{1}{2a_{\mathrm{mre}}H_{\mathrm{mre}}}\frac{a_{\mathrm{nr}}}{a_{\mathrm{mre}}}\left(\log\left(1+\frac{1}{\sqrt{1+\left(\frac{a_{\mathrm{nr}}}{a_{\mathrm{mre}}}\right)^{2}}}\right)-\log\left(1-\frac{1}{\sqrt{1+\left(\frac{a_{\mathrm{nr}}}{a_{\mathrm{mre}}}\right)^{2}}}\right)\right),$		(40)

with the velocity $v$ . The three-momentum of produced DM normalized by the mass

\displaystyle\frac{p_{N}}{m_{N}}=u=\frac{v}{\sqrt{1-v^{2}}},

(41)

whose initial value

\displaystyle u(t_{\mathrm{dec}})=\frac{\sqrt{m_{Z^{\prime}}^{2}-4m_{N}^{2}}}{2m_{N}},

(42)

is given at the time of the $Z^{\prime}$ decay , $t_{\mathrm{dec}}=1/\Gamma_{Z^{\prime}}$ , red-shifts inversely proportional to the scale factor as $u\propto a(t_{\mathrm{dec}})/a(t)=a_{\mathrm{dec}}/a(t)$ . The scale factor $a_{\mathrm{nr}}$ is one at the time $t_{\mathrm{nr}}$ when $N$ becomes nonrelativistic, i.e., $u=1$ , which is evaluated as

\displaystyle u(t_{\mathrm{dec}})\frac{a_{\mathrm{dec}}}{a_{\mathrm{nr}}}=1.

(43)

At the time of the matter-radiation equality, $t_{\mathrm{mre}}$ ,

	$\displaystyle\frac{a_{0}}{a_{\mathrm{mre}}}$	$\displaystyle=1+z_{\mathrm{mre}}=\frac{\Omega_{m}h^{2}}{\Omega_{r}h^{2}},$		(44)
	$\displaystyle H_{\mathrm{mre}}$	$\displaystyle=\sqrt{2\Omega_{m}}H_{0}(1+z_{\mathrm{mre}})^{3/2},$		(45)

where $\Omega_{r(m)}$ is the density parameter of the radiation(matter), and $a_{0}$ and $H_{0}$ are the present scale factor and Hubble parameter, respectively. We find

\displaystyle\frac{a_{\mathrm{nr}}}{a_{\mathrm{mre}}}=\frac{a_{\mathrm{dec}}}{a_{\mathrm{mre}}}u(t_{\mathrm{dec}})=\sqrt{\frac{t_{\mathrm{mre}}}{t_{\mathrm{dec}}}}u(t_{\mathrm{dec}}),

(46)

from Eq. (43) and $a\propto\sqrt{t}$ in the RD universe. By substituting Eqs. (42), (44), (45), and (46) into Eq. (40), we can evaluate the free-streaming length as $\lambda_{\mathrm{fs}}=a_{0}R_{f}$ . While the non-negligible free-streaming length would be interesting from the viewpoint of small-scale problems in cold dark matter model, it should not be larger than sub-Mpc. The recent constraints on warm DM or the free-streaming scale have been reported in Ref. Iršič et al. (2017); Nadler et al. (2021).

The contours of the observed DM density $\Omega h^{2}=0.12$ in Fig. 2 appear as blueish curves for $m_{N}=1$ MeV, $100$ MeV, and $10$ GeV, from top to bottom, respectively. The dark and light blue curves correspond to $\lambda_{\mathrm{fs}}<0.01$ Mpc, and $0.01$ Mpc $<\lambda_{\mathrm{fs}}<0.1$ Mpc, respectively. The left endpoint of the light blue curve for $M_{N}=1~{}\mathrm{MeV}$ corresponds to $\lambda_{\mathrm{fs}}=0.1$ Mpc. From the line for $M_{N}=100~{}\mathrm{MeV}$ , we can read that the free-streaming length is shortened not only for $m_{Z^{\prime}}\gtrsim 2m_{N}$ but also for $m_{Z^{\prime}}\gg m_{N}$ . In the latter case, the $Z^{\prime}$ decay happens relatively early and thus there is enough time to be redshifted for the momentum of $N$ . As long as we consider $m_{Z^{\prime}}\lesssim 100$ GeV, the mass of the DM produced by this mechanism must be larger than about MeV to have a small enough free-streaming length. The region where the lifetime of the $Z^{\prime}$ boson is longer than $0.1$ second is shaded in gray. It is ruled out because the decay of such long-lived abundant $Z^{\prime}$ bosons produces energetic particles and destroys the light elements synthesized by big bang nucleosynthesis. We employed the SN1987A constraint from Ref. Croon et al. (2021) as a reference and the excluded region is colored magenta. For the constraint see also recent other discussions Shin and Yun (2022); Caputo et al. (2022). The region colored brown is excluded by electron and proton beam dump experiments (Feng et al. (2022); Asai et al. (2022) and references therein).³³3The excluded regions by beam dump experiments and SN1987A have been derived under the assumption that $Z^{\prime}$ does not decay into RH neutrinos. Therefore, although Fig. 2 is presented for the $m_{Z^{\prime}}>2m_{N}$ case, the shaded area would not be exact. We, however, show them for reference purposes, because the excluded region is very far away from the parameter region of our interest and the viable parameter space of DM is unaffected. Since sterile neutrino $N$ is also a very weakly-interacting light particle, $N$ in addition to $Z^{\prime}$ would contribute the energy loss of supernovae. However, while the production cross section of $Z^{\prime}$ is as $\sigma\propto g_{B-L}^{2}$ , that of $N$ pair is $\sigma\propto g_{B-L}^{8}$ because the vertex of $f\rightarrow fNN$ is induced by the one loop diagram running $N$ and $Z^{\prime}$ with $f$ being SM particles. Thus, the latter is negligible compared with the former, for $g_{B-L}\ll 1$ .

IV.2 Light gauge boson region : $2m_{N}>m_{Z^{\prime}}>1$ MeV

Next, we consider the mass spectrum of $2m_{N}>m_{Z^{\prime}}>1$ MeV, in which the lower limit corresponds to the typical temperature of BBN. If we consider a smaller mass of $Z^{\prime}$ or larger coupling $g_{B-L}$ than the parameter sets studied in Sec. IV.1, the $Z^{\prime}$ gauge boson is fully thermalized by its decay and inverse decay described by the fourth terms in the right-handed side of Eq. (26b).

As a result, $N$ is dominantly produced by the pair-production processes $Z^{\prime}Z^{\prime}\rightarrow NN$ via $t$ ( $u$ )-channel $N$ exchange processes and $s$ -channel $\phi$ and $h$ exchange processes. The remarkable feature is the $m_{\phi}$ dependence of the DM abundance. For $m_{Z^{\prime}}<m_{N}$ , the scattering cross section of this process grows with respect to $s$ at lower energy than the mediator mass scale due to the longitudinal mode of $Z^{\prime}$ . In fact, the leading part of invariant amplitude squared, whose full expression is noted in Appendix, at a large $s>m_{N}^{2}\gg m_{Z^{\prime}}^{2}$ is

\displaystyle\int\sum|\mathcal{M}|^{2}d\cos\theta\sim\frac{32g_{B-L}^{4}}{m_{Z^{\prime}}^{4}}\left(\frac{2m_{N}^{4}s^{2}}{m_{N}^{2}\left(s-4m_{Z^{\prime}}^{2}\right)+m_{Z^{\prime}}^{4}}-\frac{2m_{N}^{2}s^{3}}{\Gamma_{\phi}^{2}m_{\phi}^{2}+\left(m_{\phi}^{2}-s\right)^{2}}\right),

(47)

where $\theta$ is the scattering angle. The first term from the $t(u)$ -channel $N$ exchange processes grows linearly with respect to the center-of-mass energy $s$ , as $s$ becomes large. The second term from $s$ -channel $\phi$ exchange processes becomes comparable for $s\gtrsim m_{\phi}^{2}$ and cancels with the first term at $s\gg m_{\phi}^{2}$ . The $m_{\phi}$ dependence in the thermally averaged cross section for $m_{N}=1$ GeV is shown in the left panel of Fig. 3. Each curve is for $m_{\phi}=1$ GeV (gray dotted), $10$ GeV (black dashed), and $100$ GeV (black solid). The black dashed and solid curves for $m_{\phi}>m_{N}$ have a sharp peak at $T\sim m_{\phi}/5$ due to the resonance pole of $\phi$ and decrease at the high temperature region $T>\mathcal{O}(m_{\phi})$ . It should be emphasized here that the DM production takes place most effectively at not $T\sim m_{N}$ but $T\sim m_{\phi}/5$ , unlike most freeze-in scenarios with renormalizable couplings. Hence, the result is sensitive to the mass of mediator $m_{\phi}$ .

This characteristic dependence of $m_{\phi}$ can be seen concretely in an example of the evolution of $Y_{Z^{\prime}}$ and $Y_{N}$ for $m_{N}=1$ GeV, $m_{\phi}=100$ GeV, $m_{Z^{\prime}}=0.1$ GeV, $g_{B-L}=2\times 10^{-8}$ , $\alpha=0$ presented in right panel of Fig. 3. The blue and green curves represent values of $Y_{N}$ and $Y_{Z^{\prime}}$ , respectively. The $Z^{\prime}$ boson is thermalized and its yield follows the thermal value. The slight increase of $Y_{Z^{\prime}}$ around $m_{N}/T\sim 10$ is due to the change of $g_{*}$ at the quark hadron transition. We can see that $N$ is gradually generated until the temperature becomes as low as $\sim m_{\phi}/5$ not the DM mass $m_{N}$ . The energy density of $Z^{\prime}$ decreases due to the Boltzmann suppression. Thus, the energy density of $Z^{\prime}$ gets negligible by the onset of the BBN.

The contours of $\Omega h^{2}=0.12$ in this mass spectrum are shown in Fig. 4. The colored, excluded regions are the same as Fig. 2. The curves shift to the right or left when we increase or decrease $m_{N}$ . We also note that parameters to reproduce the observed DM abundance interestingly lie in the reach of future experiments for long-lived particle search or beam dump experiments.

IV.3 Very light gauge boson region : $2m_{N}>1$ MeV $>m_{Z^{\prime}}$

Finally, we consider the mass spectrum of $2m_{N}>1$ MeV $>m_{Z^{\prime}}$ . If such a light $Z^{\prime}$ boson is thermalized by its inverse decay, the contribution of $Z^{\prime}$ to the energy density of relativistic degrees of freedom at the BBN conflicts with the observation. To avoid this, the gauge coupling constant must be much smaller than about $10^{-10}$ . The smallness of the coupling makes all processes for sterile neutrino production via $g_{B-L}$ gauge interaction negligible. However, $f\bar{f}\rightarrow NN$ , $W^{+}W^{-}(ZZ)\rightarrow NN$ and $hh(h\phi,\phi\phi)\rightarrow NN$ of all the s-channel scalar ( $h$ or $\phi$ ) exchanges are relevant processes to produce the DM $N$ . Those processes get effective at the electroweak and $B-L$ broken vacuum, and thus only $b\bar{b}$ and $\tau\bar{\tau}$ initial states are dominant, $\phi\phi$ initial state could be non-negligible for some parameter sets, and the other initial states are negligible.

Since the evolution of abundance of $N$ and $Z^{\prime}$ are independent in this case, two Boltzmann equations (26a) and (26b) are decoupled. We can solve both individually. First, let us see the evolution of the abundance of the gauge boson $Z^{\prime}$ ,

\displaystyle\frac{dn_{Z^{\prime}}}{dt}+3Hn_{Z^{\prime}}=

\displaystyle-\sum\langle\Gamma(Z^{\prime}\leftrightarrow ij)\rangle(n_{Z^{\prime}}-n^{\mathrm{eq}}_{Z^{\prime}}),

(48)

where we omit negligible inverse processes and negligible terms $\sigma v(ij\leftrightarrow Z^{\prime}Z^{\prime})$ , $\Gamma(i\leftrightarrow Z^{\prime}Z^{\prime})$ , and the $\gamma-Z^{\prime}$ scattering. The term in the right-handed side is the decay and the inverse decay of $Z^{\prime}$ . The final abundance is determined by $g_{B-L}$ , because the magnitude of the rate is proportional to $g_{B-L}^{2}$ . In Fig. 5, the evolution of the energy density of $Z^{\prime}$ for $g_{B-L}=5\times 10^{-12}$ is shown and compared with that of one generation of neutrino. Because of the very long lifetime of $Z^{\prime}$ , the constraint from CMB is much more stringent than that from the BBN. As seen in Fig. 5, even if the abundance of $Z^{\prime}$ is sufficiently small at the time of the BBN $T_{\nu}\sim 1$ MeV, the energy fraction starts to increase for $T\lesssim m_{Z^{\prime}}$ because the energy density decreases as $a^{-3}$ after $Z^{\prime}$ becomes nonrelativistic. Thus, the extra radiation generated by $Z^{\prime}$ decay could be significant at the recombination epoch and affects the temperature anisotropy of CMB unless the energy density of $Z^{\prime}$ is small sufficiently at the BBN era.

While the evolution of the energy density of $Z^{\prime}$ can be followed as above, the abundance of $Z^{\prime}$ can be easily and directly estimated by the following single integration expression

	$\displaystyle Y_{Z^{\prime}}=$	$\displaystyle\int_{T_{0}}^{T_{R}}\frac{\sum\langle\Gamma(Z^{\prime}\leftrightarrow ij)\rangle n_{i}n_{j}}{sTH}dT$
	$\displaystyle=$	$\displaystyle\frac{135\sqrt{10}M_{P}}{2\pi^{3}}\int_{T_{0}}^{T_{R}}\frac{dT}{g_{S}\sqrt{g_{}}T^{6}}\sum_{i,j}\langle\Gamma(Z^{\prime}\leftrightarrow ij)\rangle n_{Z^{\prime}}^{\mathrm{eq}},$		(49)

where $T_{0}$ is a low temperature before the decay of $Z^{\prime}$ . The $Z^{\prime}$ decays into LH neutrino pairs at $t=1/\Gamma_{Z^{\prime}}$ . By comparing the energy density of neutrino of one species

\rho_{\nu}=g_{\nu}\frac{7}{8}\frac{\pi^{2}}{30}T_{\nu}^{4},

(50)

where $g_{\nu}=2$ is the internal degrees of freedom of neutrinos and $T_{\nu}$ is the temperature of neutrinos, the energy density of decay products from $Z^{\prime}$ can be parametrized as

\displaystyle\Delta N_{\mathrm{eff}}\equiv\left.\frac{\rho_{Z^{\prime}}}{\rho_{\nu}}\right|_{t=1/\Gamma_{Z^{\prime}}}.

(51)

Similarly, by integrating Eq. (26a) from a low temperature $T_{0}$ to the reheating temperature after inflation $T_{R}$ as

	$\displaystyle Y_{N}$	$\displaystyle=\int_{T_{0}}^{T_{R}}\frac{\langle\sigma v(ij\rightarrow NN)\rangle n_{i}n_{j}}{sTH}dT$
		$\displaystyle=\frac{135\sqrt{10}M_{P}}{64\pi^{7}}\int_{T_{0}}^{T_{R}}\frac{dT}{g_{S}\sqrt{g_{}}T^{5}}\sum_{i,j}\int^{\infty}_{(m_{i}+m_{j})^{2}}dsg_{i}g_{j}p_{ij}4E_{i}E_{j}\sigma vK_{1}\left(\frac{\sqrt{s}}{T}\right),$		(52)

in the second equality, we have used Eqs. (27) and (28). In this framework, there are five free parameters associated with the production of DM $N$ ; $g_{B-L},m_{Z^{\prime}},m_{N},m_{\phi}$ , and $\alpha$ . Although $\Omega h^{2}$ seems to be dependent on all the five parameters, it practically depends on only three of $g_{B-L}\sin(2\alpha)/m_{Z^{\prime}},m_{\phi}$ , and $m_{N}$ . We can find a simple scaling of the resultant abundance as

\displaystyle\Omega h^{2}

\displaystyle\propto\left(\frac{g_{B-L}\sin(2\alpha)}{m_{Z^{\prime}}}\right)^{2},

(53)

because the coupling vertex appears only in this combination for the Higgs portal main processes. As we discussed in Sec. II, the Higgs mixing is constrained as $\alpha\lesssim 0.1$ for $m_{\phi}\gtrsim 10$ GeV and $\alpha\lesssim 10^{-4}$ for $m_{\phi}\lesssim$ several GeV.

In Fig. 6 the contours of $\Omega h^{2}=0.12$ are shown with bluish curves for some $m_{N}$ and for two sets of $\alpha$ and $m_{\phi}$ . Two black solid and one black dashed curves are contours of $\Delta N_{\mathrm{eff}}=0.5,0.2$ , and $0.06$ from top to bottom, respectively. We note that extra relativistic degrees of freedom with $0.2\lesssim\Delta N_{\mathrm{eff}}\lesssim 0.5$ would be favored to relax the so-called Hubble tension Aghanim et al. (2020); Seto and Toda (2021a, b) and that $\Delta N_{\mathrm{eff}}=0.06$ is the expected reach of the future experiment CMB-S4 Abazajian et al. (2019). The excluded region by the horizontal branch star and the red giant star constraints are shaded with orange and brown, respectively Redondo and Raffelt (2013); Heeck (2014). In the left panel, three lines of $\Omega_{N}h^{2}\simeq 0.12$ are shown for $m_{N}=0.1,1$ , and $10$ GeV. The uneven intervals are due to the change of available modes. Namely, for a smaller $m_{N}$ , the pair-production mode with heavy SM fermion initial states are suppressed at $T\sim m_{N}$ . The right plot is an example with smaller $m_{\phi}$ and $\alpha$ . In this case, the larger $m_{N}>\mathcal{O}(10)$ GeV is required to reproduce the desired DM abundance, because the production cross section is strongly suppressed by the tiny mixing $\alpha$ .

V Implication of sterile neutrino DM detection

We have considered sterile neutrino DM whose mass is larger than MeV. The visible decay modes include $N\rightarrow\nu\gamma$ , three body decay $N\rightarrow\nu f\bar{f}$ through off-shell $W^{\pm}$ and $Z$ , and even hadronic mode.⁴⁴4For formulas of those decay rates, see for instance Refs. Atre et al. (2009); Ballett et al. (2017). While the usual keV scale sterile neutrino DM is searched for with x-ray lines induced by its radiative decay Pal and Wolfenstein (1982), the DM argued in this paper can be also probed by seeking other modes, such as the decay into $e^{-}e^{+}$ and the continuous spectrum of gamma rays from hadrons.

All the decays are induced by the SM processes through the active-sterile mixing for $m_{N}<m_{Z^{\prime}}$ . On the other hand, for $m_{N}>m_{Z^{\prime}}$ , a new decay mode $N\rightarrow\nu Z^{\prime}$ is additionally possible, and thus signals to search for $N$ are decided by the decay modes of $Z^{\prime}$ depending on $m_{Z^{\prime}}$ . In the very light gauge boson region, $Z^{\prime}$ can decay into only neutrinos, while decaying modes into other fermions opens as the mass increases.

Therefore, some of the decay rates depend on whether $Z^{\prime}$ is heavier than $N$ or not, for example, the rate of $N\rightarrow e^{-}e^{+}\nu$ . Since $N\rightarrow\nu Z^{\prime}$ is a tree-level process, the constraint on the active-sterile mixing may be more stringent than that from x-ray observations, especially for $m_{N}>m_{Z^{\prime}}$ cases. The detailed study is beyond the scope of this paper, and we will evaluate this issue elsewhere.

VI Conclusion

We have evaluated the freeze-in production of $U(1)_{B-L}$ gauge interacting sterile neutrino DM by taking into account processes overlooked in the literature; principally the inverse decay of $Z^{\prime}$ and the longitudinal mode effect in $Z^{\prime}Z^{\prime}\rightarrow NN$ scattering. We have found that the inverse decay of $Z^{\prime}$ indeed gives a non-negligible contribution to the production of $Z^{\prime}$ .

For $m_{Z^{\prime}}>2m_{N}$ cases, the final value of $\Omega_{N}h^{2}$ agrees with the previous estimation. Our finding for this mass spectrum case is that the constraint from the free-streaming length is more stringent than that has been thought. As the result, the mass of sterile neutrino DM under this mass spectrum must be larger than about one MeV as long as we assume $m_{Z^{\prime}}\lesssim 100$ GeV.

For the other spectrum $m_{Z^{\prime}}<2m_{N}$ , the gauge coupling $g_{B-L}=\mathcal{O}(10^{-6})$ independent from the $m_{Z^{\prime}}$ has been regarded as the viable parameter region to reproduce the desired DM abundance. We, however, have found that this is not correct. Since $Z^{\prime}$ is lighter than $N$ , when $N$ is produced by $Z^{\prime}Z^{\prime}\rightarrow NN$ , this cross section is enhanced by small $m_{Z^{\prime}}^{2}$ and increases with respect to the energy until $s\sim m_{\phi}^{2}$ , due to the longitudinal mode of $Z^{\prime}$ . Thus, the resultant DM abundance depends on $m_{Z^{\prime}}$ and $m_{\phi}$ . In addition, for such a large coupling of $g_{B-L}\sim 10^{-6}$ , $Z^{\prime}$ can be thermalized and gives $\Delta N_{\mathrm{eff}}\sim 1$ at the BBN epoch for $m_{Z^{\prime}}\lesssim 1$ MeV. Moreover, CMB gives a more stringent constraint on the gauge coupling as $g_{B-L}\lesssim 10^{-12}$ than BBN, because the energy density of $Z^{\prime}$ decreases slower than the background radiation after it becomes nonrelativistic. Thus, for $2m_{N}>1$ MeV $>m_{Z^{\prime}}$ , the freeze-in production by the $Z^{\prime}$ mediation is not available and the DM in this parameter region can be produced only by scalar portal scatterings.

Having the parameter space of consistent sterile neutrino DM mentioned above, the viable parameters for $m_{Z^{\prime}}>2m_{N}$ lie far beyond the reach of the near future experiments of long-lived particle searches as shown in Fig. 2, while a large interesting parameter space in the spectrum of $1$ MeV $<m_{Z^{\prime}}<2m_{N}$ is already constrained partially by the current experimental limits and will be probed more by the experiments in future. The case of the spectrum with $m_{Z^{\prime}}<1$ MeV $<2m_{N}$ will be examined by future measurements of $N_{\mathrm{eff}}$ .

Acknowledgments

We thank Masahiro Ibe for valuable comments on the IR divergence of the $Z^{\prime}-\gamma$ scattering. This work is supported, in part, by JSPS KAKENHI Grants No. JP19K03860 and No. JP19K03865 and MEXT KAKENHI Grant No. 21H00060 (O. S.), and JSPS KAKENHI Grants No. 18K03651, No. 18H01210, No. 22K03622 and MEXT KAKENHI Grant No. 18H05543 (T. S.).

Appendix A Amplitude of the sterile neutrino pair-production processes

We give explicit formulas of the invariant amplitude squared.

A.1 $f(p_{1})\bar{f}(p_{2})\rightarrow N(q_{1})N(q_{2})$

\displaystyle\overline{|\mathcal{M}|^{2}}d\cos\theta=\frac{m_{f}^{2}m_{N}^{2}N_{c}\sin^{2}(2\alpha)\left(s-4m_{f}^{2}\right)\left(s-4m_{N}^{2}\right)}{2v^{2}v_{B-L}^{2}}\left|\frac{-1}{s-m_{h}^{2}-i\Gamma_{h}m_{h}}+\frac{1}{s-m_{\phi}^{2}-im_{\phi}\Gamma_{\phi}}\right|^{2},

(54)

	$\displaystyle 4E_{i}E_{j}\sigma v(ij\rightarrow)=\frac{1}{16\pi}\sqrt{\frac{s-4m_{N}^{2}}{s}}\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$		(55)
	$\displaystyle p_{ij}=\frac{\sqrt{s-4m_{f}^{2}}}{2}.$		(56)

A.2 $W(p_{1})W(p_{2})\rightarrow N(q_{1})N(q_{2})$

	$\displaystyle\overline{\|\mathcal{M}\|^{2}}d\cos\theta=$	$\displaystyle\frac{g_{2}^{2}m_{N}^{2}\sin^{2}(2\alpha)\left(s-4m_{N}^{2}\right)\left(-4m_{N}^{2}\left(2m_{W}^{2}+s\right)+4m_{N}^{4}+16m_{W}^{4}+s^{2}\right)}{36m_{W}^{2}v_{B-L}^{2}}$
		$\displaystyle\times\left\|\frac{-1}{s-m_{h}^{2}-i\Gamma_{h}m_{h}}+\frac{1}{s-m_{\phi}^{2}-im_{\phi}\Gamma_{\phi}}\right\|^{2},$		(57)

	$\displaystyle 4E_{i}E_{j}\sigma v(ij\rightarrow)=\frac{1}{16\pi}\sqrt{\frac{s-4m_{N}^{2}}{s}}\int\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$		(58)
	$\displaystyle p_{ij}=\frac{\sqrt{s-4m_{W}^{2}}}{2}.$		(59)

A.3 $Z(p_{1})Z(p_{2})\rightarrow N(q_{1})N(q_{2})$

	$\displaystyle\overline{\|\mathcal{M}\|^{2}}d\cos\theta=$	$\displaystyle\frac{g_{2}^{2}m_{N}^{2}\sin^{2}(2\alpha)\left(s-4m_{N}^{2}\right)\left(-4m_{N}^{2}\left(2m_{Z}^{2}+s\right)+4m_{N}^{4}+16m_{Z}^{4}+s^{2}\right)}{36m_{Z}^{2}c_{W}^{2}v_{B-L}^{2}}$
		$\displaystyle\times\left\|\frac{-1}{s-m_{h}^{2}-i\Gamma_{h}m_{h}}+\frac{1}{s-m_{\phi}^{2}-im_{\phi}\Gamma_{\phi}}\right\|^{2},$		(60)

	$\displaystyle 4E_{i}E_{j}\sigma v(ij\rightarrow)=\frac{1}{16\pi}\sqrt{\frac{s-4m_{N}^{2}}{s}}\int\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$		(61)
	$\displaystyle p_{ij}=\frac{\sqrt{s-4m_{Z}^{2}}}{2}.$		(62)

A.4 $(h,\phi)(p_{1})(h,\phi)(p_{2})\rightarrow N(p_{1})N(p_{2})$

\displaystyle i\mathcal{M}=

\displaystyle-\frac{im_{N}\bar{u}(p_{1},m_{N})v(p_{2},m_{N})}{v_{B-L}}\left(\frac{C_{h}\sin\alpha}{(q_{1}+q_{2})^{2}-m_{h}^{2}+i\Gamma_{h}m_{h}}+\frac{C_{\phi}\cos\alpha}{(q_{1}+q_{2})^{2}-m_{\phi}^{2}+im_{\phi}\Gamma_{\phi}}\right),

(63)

where we have omitted $t(u)$ -channel $N$ exchange contributions because those are suppressed by $m_{N}^{2}/v_{B-L}^{2}$ and are thus negligible.

\displaystyle\overline{|\mathcal{M}_{s}|^{2}}d\cos\theta=\frac{4m_{N}^{2}(s-4m_{N}^{2})}{v_{B-L}^{2}}\left|\frac{C_{h}\sin\alpha}{s-m_{h}^{2}+i\Gamma_{h}m_{h}}+\frac{C_{\phi}\cos\alpha}{s-m_{\phi}^{2}+im_{\phi}\Gamma_{\phi}}\right|^{2},

(64)

	$\displaystyle(C_{h},C_{\phi})=$	$\displaystyle(C_{h\phi\phi},C_{\phi\phi\phi})\quad\mathrm{for}\quad(\phi,\phi),$
		$\displaystyle(C_{hh\phi},C_{h\phi\phi})\quad\mathrm{for}\quad(h,\phi),$
		$\displaystyle(C_{hhh},C_{hh\phi})\quad\mathrm{for}\quad(h,h),$

	$\displaystyle 4E_{i}E_{j}\sigma v(ij\rightarrow)=\frac{1}{16\pi}\sqrt{\frac{s-4m_{N}^{2}}{s}}\int\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$		(65)
	$\displaystyle p_{ij}=\frac{1}{2}\sqrt{\frac{\left(s-(m_{\phi_{i}}-m_{\phi_{j}})^{2}\right)\left(s-(m_{\phi_{i}}+m_{\phi_{j}})^{2}\right)}{s}}.$		(66)

A.5 $Z^{\prime}(q_{1})Z^{\prime}(q_{2})\rightarrow N(p_{1})N(p_{2})$

	$\displaystyle\int\sum\|\mathcal{M}\|^{2}d\cos\theta$
$\displaystyle=$	$\displaystyle\frac{32g_{B-L}^{4}}{m_{Z^{\prime}}^{4}}\left(\frac{4m_{N}^{2}\left(m_{\phi}^{2}-s\right)\left(-2m_{Z^{\prime}}^{2}s+4m_{Z^{\prime}}^{4}+s^{2}\right)}{\Gamma_{\phi}^{2}m_{\phi}^{2}+\left(m_{\phi}^{2}-s\right)^{2}}-\frac{2m_{N}^{2}\left(4m_{N}^{2}-s\right)\left(-4m_{Z^{\prime}}^{2}s+12m_{Z^{\prime}}^{4}+s^{2}\right)}{\Gamma_{\phi}^{2}m_{\phi}^{2}+\left(m_{\phi}^{2}-s\right)^{2}}\right.$
	$\displaystyle\left.+\frac{2m_{N}^{4}\left(-8m_{Z^{\prime}}^{2}s+8m_{Z^{\prime}}^{4}+s^{2}\right)+m_{N}^{2}m_{Z^{\prime}}^{4}\left(4m_{Z^{\prime}}^{2}+s\right)-2m_{Z^{\prime}}^{8}}{m_{N}^{2}\left(s-4m_{Z^{\prime}}^{2}\right)+m_{Z^{\prime}}^{4}}\right)$
	$\displaystyle+\frac{32g_{B-L}^{4}}{m_{Z^{\prime}}^{4}\sqrt{s-4m_{N}^{2}}\sqrt{s-4m_{Z^{\prime}}^{2}}}\left(\frac{8m_{N}^{2}\left(s-m_{\phi}^{2}\right)\left(m_{N}^{2}\left(-4m_{Z^{\prime}}^{2}s+8m_{Z^{\prime}}^{4}+s^{2}\right)-2m_{Z^{\prime}}^{6}\right)}{\Gamma_{\phi}^{2}m_{\phi}^{2}+\left(m_{\phi}^{2}-s\right)^{2}}\right.$
	$\displaystyle\left.-\frac{4m_{N}^{4}s\left(s-4m_{Z^{\prime}}^{2}\right)+4m_{N}^{2}m_{Z^{\prime}}^{2}\left(4m_{Z^{\prime}}^{2}-s\right)\left(m_{Z^{\prime}}^{2}+s\right)-m_{Z^{\prime}}^{4}\left(4m_{Z^{\prime}}^{4}+s^{2}\right)}{s-2m_{Z^{\prime}}^{2}}\right)$
	$\displaystyle\times\log\left(\frac{s-2m_{Z^{\prime}}^{2}+\sqrt{s-4m_{N}^{2}}\sqrt{s-4m_{Z^{\prime}}^{2}}}{s-2m_{Z^{\prime}}^{2}-\sqrt{s-4m_{N}^{2}}\sqrt{s-4m_{Z^{\prime}}^{2}}}\right),$	(67)

	$\displaystyle 4E_{i}E_{j}\sigma v(ij\rightarrow)=\frac{1}{16\pi}\sqrt{\frac{s-4m_{N}^{2}}{s}}\int\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$		(68)
	$\displaystyle p_{ij}=\frac{\sqrt{s-4m_{Z^{\prime}}^{2}}}{2}.$		(69)

Appendix B Amplitude of the $\gamma-Z^{\prime}$ scattering processes

We give explicit formulas of the invariant amplitude squared.

B.1 $Z^{\prime}(q_{1})\gamma(q_{2})\rightarrow f(p_{1})\bar{f}(p_{2})$

$\displaystyle\int\sum\|\mathcal{M}\|^{2}d\cos\theta=$	$\displaystyle(g_{B-L}q_{Xf}eq_{f})^{2}\frac{32}{(s-m_{Z^{\prime}}^{2})^{2}}$
	$\displaystyle\times\left(\frac{\sqrt{s}\left(4m_{f}^{2}\left(s-m_{Z^{\prime}}^{2}\right)-8m_{f}^{4}+m_{Z^{\prime}}^{4}+s^{2}\right)}{\sqrt{s-4m_{f}^{2}}}\log\left(\frac{\sqrt{s}+\sqrt{s-4m_{f}^{2}}}{\sqrt{s}-\sqrt{s-4m_{f}^{2}}}\right)\right.$
	$\displaystyle\left.-s\left(4m_{f}^{2}+s\right)-m_{Z^{\prime}}^{4}\right),$	(70)

where the IR divergence at $s=m_{Z^{\prime}}^{2}$ is due to the on-shell $t(u)$ -channel mediator and can be regulated by introducing thermal photon mass Redondo and Postma (2009).

	$\displaystyle 4E_{i}E_{j}\sigma v(ij\rightarrow)=\frac{1}{16\pi}\sqrt{\frac{s-4m_{f}^{2}}{s}}\int\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$		(71)
	$\displaystyle p_{ij}=\frac{1}{2}\frac{s-m_{Z^{\prime}}^{2}}{\sqrt{s}}.$		(72)

B.2 $Z^{\prime}(q_{1})f(q_{2})\rightarrow\gamma(p_{1})f(p_{2})$

$\displaystyle\int\sum\|\mathcal{M}\|^{2}d\cos\theta=$	$\displaystyle\frac{8(g_{B-L}q_{Xf}eq_{f})^{2}}{s\left(m_{f}^{2}-s\right)^{2}}\left(-m_{f}^{4}\left(m_{Z^{\prime}}^{2}+s\right)+m_{f}^{2}s\left(2m_{Z^{\prime}}^{2}+15s\right)+m_{f}^{6}+s^{2}\left(7m_{Z^{\prime}}^{2}+s\right)\right.$
	$\displaystyle\left.+\frac{2s^{2}\left(2m_{f}^{2}\left(m_{Z^{\prime}}^{2}-3s\right)-3m_{f}^{4}-2m_{Z^{\prime}}^{2}s+2m_{Z^{\prime}}^{4}+s^{2}\right)}{\sqrt{\left(s-(m_{f}-m_{Z^{\prime}})^{2}\right)\left(s-(m_{f}+m_{Z^{\prime}})^{2}\right)}}\right.$
	$\displaystyle\left.\times\log\left(\frac{m_{f}^{2}-m_{Z^{\prime}}^{2}+s+\sqrt{\left(s-(m_{f}-m_{Z^{\prime}})^{2}\right)\left(s-(m_{f}+m_{Z^{\prime}})^{2}\right)}}{m_{f}^{2}-m_{Z^{\prime}}^{2}+s-\sqrt{\left(s-(m_{f}-m_{Z^{\prime}})^{2}\right)\left(s-(m_{f}+m_{Z^{\prime}})^{2}\right)}}\right)\right),$	(73)

	$\displaystyle 4E_{i}E_{j}\sigma v(ij\rightarrow)=\frac{1}{16\pi}\sqrt{\frac{s-m_{f}^{2}}{s}}\int\overline{\|\mathcal{M}\|^{2}}d\cos\theta,$		(74)
	$\displaystyle p_{ij}=\frac{1}{2}\sqrt{\frac{\left(s-(m_{f}-m_{Z^{\prime}})^{2}\right)\left(s-(m_{f}+m_{Z^{\prime}})^{2}\right)}{s}}.$		(75)

B.3 $Z^{\prime}(q_{1})\bar{f}(q_{2})\rightarrow\gamma(p_{1})\bar{f}(p_{2})$

It is same as for $Z^{\prime}(q_{1})f(q_{2})\rightarrow\gamma(p_{1})f(p_{2})$ .

References

Minkowski (1977) P. Minkowski, Phys. Lett. B 67, 421 (1977).
Yanagida (1979) T. Yanagida, Conf. Proc. C 7902131, 95 (1979).
Gell-Mann et al. (1979) M. Gell-Mann, P. Ramond, and R. Slansky, Conf. Proc. C 790927, 315 (1979), eprint 1306.4669.
Mohapatra and Senjanovic (1980) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
Dolgov and Hansen (2002) A. Dolgov and S. Hansen, Astropart. Phys. 16, 339 (2002), eprint hep-ph/0009083.
Drewes et al. (2017) M. Drewes et al., JCAP 01, 025 (2017), eprint 1602.04816.
Boyarsky et al. (2019) A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens, and O. Ruchayskiy, Prog. Part. Nucl. Phys. 104, 1 (2019), eprint 1807.07938.
Asaka et al. (2005) T. Asaka, S. Blanchet, and M. Shaposhnikov, Phys. Lett. B 631, 151 (2005), eprint hep-ph/0503065.
Asaka and Shaposhnikov (2005) T. Asaka and M. Shaposhnikov, Phys. Lett. B 620, 17 (2005), eprint hep-ph/0505013.
Davidson (1979) A. Davidson, Phys. Rev. D 20, 776 (1979).
Mohapatra and Marshak (1980) R. N. Mohapatra and R. Marshak, Phys. Rev. Lett. 44, 1316 (1980), [Erratum: Phys.Rev.Lett. 44, 1643 (1980)].
Marshak and Mohapatra (1980) R. Marshak and R. N. Mohapatra, Phys. Lett. B 91, 222 (1980).
Dodelson and Widrow (1994) S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72, 17 (1994), eprint hep-ph/9303287.
Boyarsky et al. (2006a) A. Boyarsky, A. Neronov, O. Ruchayskiy, and M. Shaposhnikov, Mon. Not. Roy. Astron. Soc. 370, 213 (2006a), eprint astro-ph/0512509.
Boyarsky et al. (2006b) A. Boyarsky, A. Neronov, O. Ruchayskiy, M. Shaposhnikov, and I. Tkachev, Phys. Rev. Lett. 97, 261302 (2006b), eprint astro-ph/0603660.
Boyarsky et al. (2007) A. Boyarsky, J. Nevalainen, and O. Ruchayskiy, Astron. Astrophys. 471, 51 (2007), eprint astro-ph/0610961.
Boyarsky et al. (2008) A. Boyarsky, D. Iakubovskyi, O. Ruchayskiy, and V. Savchenko, Mon. Not. Roy. Astron. Soc. 387, 1361 (2008), eprint 0709.2301.
Yuksel et al. (2008) H. Yuksel, J. F. Beacom, and C. R. Watson, Phys. Rev. Lett. 101, 121301 (2008), eprint 0706.4084.
Khalil and Seto (2008) S. Khalil and O. Seto, JCAP 10, 024 (2008), eprint 0804.0336.
Kaneta et al. (2017) K. Kaneta, Z. Kang, and H.-S. Lee, JHEP 02, 031 (2017), eprint 1606.09317.
Biswas and Gupta (2016) A. Biswas and A. Gupta, JCAP 09, 044 (2016), [Addendum: JCAP 05, A01 (2017)], eprint 1607.01469.
Seto and Shimomura (2020) O. Seto and T. Shimomura, Phys. Lett. B 811, 135880 (2020), eprint 2007.14605.
De Romeri et al. (2020) V. De Romeri, D. Karamitros, O. Lebedev, and T. Toma, JHEP 10, 137 (2020), eprint 2003.12606.
Lucente (2021) M. Lucente (2021), eprint 2103.03253.
Bélanger et al. (2021) G. Bélanger, S. Khan, R. Padhan, M. Mitra, and S. Shil, Phys. Rev. D 104, 055047 (2021), eprint 2104.04373.
Hall et al. (2010) L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, JHEP 03, 080 (2010), eprint 0911.1120.
Baer et al. (2015) H. Baer, K.-Y. Choi, J. E. Kim, and L. Roszkowski, Phys. Rept. 555, 1 (2015), eprint 1407.0017.
Shakya (2016) B. Shakya, Mod. Phys. Lett. A 31, 1630005 (2016), eprint 1512.02751.
Fileviez Pérez et al. (2019) P. Fileviez Pérez, C. Murgui, and A. D. Plascencia, Phys. Rev. D 100, 035041 (2019), eprint 1905.06344.
Heeba and Kahlhoefer (2020) S. Heeba and F. Kahlhoefer, Phys. Rev. D 101, 035043 (2020), eprint 1908.09834.
Okada et al. (2020) N. Okada, S. Okada, and Q. Shafi, Phys. Lett. B 810, 135845 (2020), eprint 2003.02667.
Okada et al. (2021) H. Okada, Y. Orikasa, and Y. Shoji, JCAP 07, 006 (2021), eprint 2102.10944.
Iwamoto et al. (2022) S. Iwamoto, K. Seller, and Z. Trócsányi, JCAP 01, 035 (2022), eprint 2104.11248.
Caputo et al. (2021) A. Caputo, A. J. Millar, C. A. J. O’Hare, and E. Vitagliano, Phys. Rev. D 104, 095029 (2021), eprint 2105.04565.
Zyla et al. (2020) P. A. Zyla et al. (Particle Data Group), PTEP 2020, 083C01 (2020).
Aad et al. (2020) G. Aad et al. (ATLAS), Phys. Rev. D 101, 012002 (2020), eprint 1909.02845.
Barate et al. (2003) R. Barate et al. (LEP Working Group for Higgs boson searches, ALEPH, DELPHI, L3, OPAL), Phys. Lett. B 565, 61 (2003), eprint hep-ex/0306033.
Aaij et al. (2015) R. Aaij et al. (LHCb), Phys. Rev. Lett. 115, 161802 (2015), eprint 1508.04094.
Aaij et al. (2017) R. Aaij et al. (LHCb), Phys. Rev. D 95, 071101 (2017), eprint 1612.07818.
Bergsma et al. (1985) F. Bergsma et al. (CHARM), Phys. Lett. B 157, 458 (1985).
Winkler (2019) M. W. Winkler, Phys. Rev. D 99, 015018 (2019), eprint 1809.01876.
Fradette and Pospelov (2017) A. Fradette and M. Pospelov, Phys. Rev. D 96, 075033 (2017), eprint 1706.01920.
Gondolo and Gelmini (1991) P. Gondolo and G. Gelmini, Nucl. Phys. B 360, 145 (1991).
Kolb and Turner (1990) E. W. Kolb and M. S. Turner, The Early Universe, vol. 69 of Frontiers in physics (Addison-Wesley, 1990), ISBN 978-0-201-116038.
Iršič et al. (2017) V. Iršič et al., Phys. Rev. D 96, 023522 (2017), eprint 1702.01764.
Nadler et al. (2021) E. O. Nadler et al. (DES), Phys. Rev. Lett. 126, 091101 (2021), eprint 2008.00022.
Croon et al. (2021) D. Croon, G. Elor, R. K. Leane, and S. D. McDermott, JHEP 01, 107 (2021), eprint 2006.13942.
Shin and Yun (2022) C. S. Shin and S. Yun, JHEP 02, 133 (2022), eprint 2110.03362.
Caputo et al. (2022) A. Caputo, G. Raffelt, and E. Vitagliano, JCAP 08, 045 (2022), eprint 2204.11862.
Feng et al. (2022) J. L. Feng et al. (2022), eprint 2203.05090.
Asai et al. (2022) K. Asai, A. Das, J. Li, T. Nomura, and O. Seto (2022), eprint 2206.12676.
Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
Seto and Toda (2021a) O. Seto and Y. Toda, Phys. Rev. D 103, 123501 (2021a), eprint 2101.03740.
Seto and Toda (2021b) O. Seto and Y. Toda, Phys. Rev. D 104, 063019 (2021b), eprint 2104.04381.
Abazajian et al. (2019) K. Abazajian et al. (2019), eprint 1907.04473.
Redondo and Raffelt (2013) J. Redondo and G. Raffelt, JCAP 08, 034 (2013), eprint 1305.2920.
Heeck (2014) J. Heeck, Phys. Lett. B 739, 256 (2014), eprint 1408.6845.
Atre et al. (2009) A. Atre, T. Han, S. Pascoli, and B. Zhang, JHEP 05, 030 (2009), eprint 0901.3589.
Ballett et al. (2017) P. Ballett, S. Pascoli, and M. Ross-Lonergan, JHEP 04, 102 (2017), eprint 1610.08512.
Pal and Wolfenstein (1982) P. B. Pal and L. Wolfenstein, Phys. Rev. D 25, 766 (1982).
Redondo and Postma (2009) J. Redondo and M. Postma, JCAP 02, 005 (2009), eprint 0811.0326.

Revisiting sterile neutrino dark matter in gauged U​(1)B−LU(1)_{B-L} model

Abstract

I Introduction

II Model

II.1 Lagrangian

II.2 Scalar potential

II.2.1 Mass eigenstates and Higgs mixing

II.2.2 Decays of scalars and range of Higgs mixing

III The Boltzmann equation

IV Abundance of sterile neutrino DM

IV.1 Heavy gauge boson region mZ′>2​mNm_{Z^{\prime}}>2m_{N}

IV.2 Light gauge boson region : 2​mN>mZ′>12m_{N}>m_{Z^{\prime}}>1 MeV

IV.3 Very light gauge boson region : 2​mN>12m_{N}>1 MeV >mZ′>m_{Z^{\prime}}

V Implication of sterile neutrino DM detection

VI Conclusion

Acknowledgments

Appendix A Amplitude of the sterile neutrino pair-production processes

A.1 f​(p1)​f¯​(p2)→N​(q1)​N​(q2)f(p_{1})\bar{f}(p_{2})\rightarrow N(q_{1})N(q_{2})

A.2 W​(p1)​W​(p2)→N​(q1)​N​(q2)W(p_{1})W(p_{2})\rightarrow N(q_{1})N(q_{2})

A.3 Z​(p1)​Z​(p2)→N​(q1)​N​(q2)Z(p_{1})Z(p_{2})\rightarrow N(q_{1})N(q_{2})

A.4 (h,ϕ)​(p1)​(h,ϕ)​(p2)→N​(p1)​N​(p2)(h,\phi)(p_{1})(h,\phi)(p_{2})\rightarrow N(p_{1})N(p_{2})

A.5 Z′​(q1)​Z′​(q2)→N​(p1)​N​(p2)Z^{\prime}(q_{1})Z^{\prime}(q_{2})\rightarrow N(p_{1})N(p_{2})

Appendix B Amplitude of the γ−Z′\gamma-Z^{\prime} scattering processes

B.1 Z′​(q1)​γ​(q2)→f​(p1)​f¯​(p2)Z^{\prime}(q_{1})\gamma(q_{2})\rightarrow f(p_{1})\bar{f}(p_{2})

B.2 Z′​(q1)​f​(q2)→γ​(p1)​f​(p2)Z^{\prime}(q_{1})f(q_{2})\rightarrow\gamma(p_{1})f(p_{2})

B.3 Z′​(q1)​f¯​(q2)→γ​(p1)​f¯​(p2)Z^{\prime}(q_{1})\bar{f}(q_{2})\rightarrow\gamma(p_{1})\bar{f}(p_{2})