aainstitutetext: School of Physical Sciences, Indian Association for the
Cultivation of Science,
2A & 2B Raja S.C. Mullick Road, Jadavpur, Kolkata 700 032, India

# Secluded Dark Sector and Muon $(g-2)$ in the Light of Fast Expanding
Universe

Sougata Ganguly a Sourov Roy a Ananya Tapadar<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The lack of information before Big Bang Neucleosynthesis (BBN) allow us to
assume the presence of a new species $\phi$ whose energy density redshifts as
$a^{-(4+n)}$ where $n>0$ and $a$ is the scale factor. This non-standard
cosmological setup facilitates a larger portal coupling $(\epsilon)$ between
the dark and the visible sectors even when the two sectors are not in thermal
equilibrium. Here, we have considered $U(1)_{L_{\mu}-L_{\tau}}\otimes
U(1)_{X}$ gauge extension of the Standard Model (SM) and studied different
phases of the cosmological evolution of a thermally decoupled dark sector such
as leak-in, freeze-in, reannihilation, and late-time annihilation in the
presence of fast expansion. Due to the tree level kinetic mixing between
$U(1)_{X}$ and $U(1)_{L_{\mu}-L_{\tau}}$ gauge bosons, the dark sector couples
with the $\mu$ and $\tau$ flavored leptons of the SM. We show that in our
scenario it is possible to reconcile the dark matter relic density and muon
$(g-2)$ anomaly. In particular, we show that for $2\times
10^{-4}\lesssim\epsilon\lesssim 10^{-3}$, $5.5{\rm MeV}\lesssim
m_{Z^{\prime}}\lesssim 200{\rm MeV}$, $n=4$, and $1{\rm TeV}\lesssim
m_{\chi}\lesssim 10{\rm TeV}$ relic density constraint of dark matter,
constraint from muon $(g-2)$ anomaly, and other cosmological, astrophysical
constraints are satisfied.

##  1 Introduction

The Standard Model (SM) of particle physics is an extremely successful theory
to explain observable phenomena at the microscopic level. Despite huge
successful predictions of SM, there are some unsolved issues which cannot be
explained by SM. The existence of dark matter (DM), revealed from various
cosmological and astrophysical observations (such as measurement of cosmic
microwave background(CMB) anisotropy by WMAP WMAP:2012nax , Planck 2018
Planck:2018vyg , galaxy rotation curve Sofue:2000jx , bullet cluster
observation Clowe:2003tk ; Clowe:2006eq  etc.), is one such unsolved issue
which SM cannot explain. According to the latest data from Planck 2018
observations Planck:2018vyg , the abundance of DM at the present epoch is
$\Omega_{\rm{DM}}h^{2}=0.1200\pm 0.0012$ Planck:2018vyg . Even after decades
of first observational evidences, the origin and nature of DM are still
elusive. None of the SM particles is suitable as a DM candidate and thus the
physics beyond SM is indispensable.

Weakly interacting massive particle (WIMP) has been a well motivated DM
candidate for a long time Srednicki:1988ce ; Gondolo:1990dk ; Jungman:1995df ;
Feng:2010gw . Inspite of having many attractive theoretical features, WIMP
paradigm is not suitable to explain the negative results of direct, indirect
as well as laboratory searches Roszkowski:2017nbc ; Lin:2019uvt . As a result,
the allowed parameter space of the WIMP scenario has been stringently
constrained.

On the other hand, recent measurement of the anomalous magnetic moment of muon
by Fermilab Muong-2:2021ojo  reveals that the tension between the theoretical
prediction Czarnecki:2002nt ; Melnikov:2003xd ; Aoyama:2012wk ;
Gnendiger:2013pva ; Kurz:2014wya ; Colangelo:2014qya ; Masjuan:2017tvw ;
Colangelo:2017fiz ; Davier:2017zfy ; Keshavarzi:2018mgv ; Hoferichter:2018kwz
; Colangelo:2018mtw ; Aoyama:2019ryr ; Gerardin:2019vio ; Hoferichter:2019mqg
; Bijnens:2019ghy ; Davier:2019can ; Colangelo:2019uex ; Keshavarzi:2019abf ;
Blum:2019ugy ; Aoyama:2020ynm  and experimental observations of muon $(g-2)$
is at $4.2\sigma$. Combining the result of Fermilab with the old measurement
by BNL E821 experiment Muong-2:2006rrc , the deviation of
$a_{\mu}(\equiv(g-2)/2)$ from the SM prediction is given by Muong-2:2021ojo

$\displaystyle\Delta a_{\mu}$ $\displaystyle=$ $\displaystyle\left(251\pm
59\right)\times 10^{-11}\,\,.$

Anomaly free $U(1)_{L_{\mu}-L_{\tau}}$ gauge extension of the SM He:1990pn ;
He:1991qd  is a very well motivated model to explain the mismatch between
theoretical prediction and experimental observation of muon $(g-2)$
Baek:2001kca ; Harigaya:2013twa ; Altmannshofer:2016brv .
$U(1)_{L_{\mu}-L_{\tau}}$ gauge extension has also been extensively studied in
the Supersymmetric framework to explain the anomalous magnetic moment of muon
Ma:2001md ; Banerjee:2018eaf ; Banerjee:2020zvi . In the context of neutrino
masses and mixings $U(1)_{L_{\mu}-L_{\tau}}$ model has been studied in
Ma:2001md ; Baek:2015mna ; Chen:2017gvf .

In the light of the above discussion, we have considered
$U(1)_{L_{\mu}-L_{\tau}}\otimes U(1)_{X}$ gauge extension of SM to study the
cosmological evolution of a thermally decoupled dark sector and muon $(g-2)$
anomaly. The dark sector contains the gauge boson corresponding to the
$U(1)_{X}$ gauge symmetry and the DM $\chi$ which is a Dirac fermion, charged
only under $U(1)_{X}$ gauge symmetry. In order to connect dark and SM sector,
we have considered tree level kinetic mixing between $U(1)_{L_{\mu}-L_{\tau}}$
and $U(1)_{X}$ gauge bosons and the dark sector particles couple with the
$\mu$ and $\tau$ flavored SM leptons through the portal coupling $\epsilon$.
Let us note that, other possible tree level kinetic mixings between $U(1)_{Y}$
with $U(1)_{L_{\mu}-L_{\tau}}$ and $U(1)_{X}$ gauge bosons are not considered
here. However, these kinetic mixings can be generated at one loop level.

For a thermally decoupled dark sector, there exists an upper limit on the
portal coupling $\epsilon$ and the temperature evolution of both the sectors
are very much different hambye:2019 ; Evans:2019vxr ; Tapadar:2021kgw . Since
both the sectors are thermally decoupled, the value of $\epsilon$ is very
small. Due to the small value of the portal coupling, although the dark vector
boson $Z^{\prime}$, having a coupling with muons, cannot explain muon $(g-2)$
anomaly. However, the situation may improve if we modify the cosmological
history before Big Bang Neucleosynthesis (BBN). The modification of the
cosmology in the pre-BBN era helps ameliorate the theoretical prediction of
muon $(g-2)$ in the following manner. If we assume the existence of a new
species $\phi$ whose energy density redshifts as $a^{-(4+n)}$ where $n>0$
before BBN, then the Universe expands faster in comparison to the standard
radiation dominated Universe DEramo:2017gpl ; DEramo:2017ecx . The key feature
of the fast expanding scenario is the requirement of large coupling constant
to satisfy the relic density constraint. Therefore, to study muon $(g-2)$
anomaly and phases of a thermally decoupled dark sector in a unified
framework, we have utilised the idea of fast expanding Universe to uplift the
upper limit of the portal coupling $\epsilon$. In addition to that, we have
also investigated the constraints from Borexino Bellini:2011rx , BABAR
experiments BaBar:2016sci , BBN, white dwarf cooling, self interaction of
DM,and neutrino trident production at CCFR CCFR:1991lpl .

The rest of the paper is structured as follows. In section 2 we have discussed
our model briefly. The fast expanding scenario, thermalization and temperature
evolution of the dark sector in the fast expanding Universe has been discussed
in section 3. Analytical calculations for the DM relic density in fast
expanding Universe has been performed in section 4. We provide a detail
analysis of DM number density evolution in section 5. Consequence of
$Z^{\prime}$ on muon anomalous magnetic moment, relevant constraints on portal
coupling and dark vector boson mass are discussed in section 6. Finally, we
conclude in section 7. In appendix A, we show the detailed analytical
calculations of freeze-in scenario in case of fast expanding Universe.

##  2 The model

We consider $U(1)_{L_{\mu}-L_{\tau}}\otimes U(1)_{X}$ gauge extension of the
SM to study the cosmological evolution of the dark sector and muon $(g-2)$
anomaly simultaneously. We introduce a Dirac fermion $\chi$ which is charged
only under $U(1)_{X}$ gauge symmetry. In our framework, $\chi$ acts as a DM
candidate and it only interacts with the gauge boson of $U(1)_{X}$ gauge
symmetry $\hat{Z^{\prime}}$. All the SM fields are singlet under $U(1)_{X}$
gauge symmetry. We have considered the tree level kinetic mixing between
$\hat{Z^{\prime}}$ and $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson
$\hat{Z}_{\mu\tau}$ and because of this, the dark sector couples with the
$\mu$ and $\tau$ flavored SM leptons. One of the interesting features of this
model is that the dark gauge boson i.e. the gauge boson corresponding to the
$U(1)_{X}$ gauge symmetry can contribute to the anomalous magnetic moment of
muon if it is sufficiently light. In our model we have assumed that the masses
of the gauge bosons are generated via Stückelberg mechanism
Stueckelberg:1938hvi ; Ruegg:2003ps . As mentioned in the introduction, tree
level kinetic mixing between $U(1)_{Y}$ gauge boson with
$U(1)_{L_{\mu}-L_{\tau}}$ and $U(1)_{X}$ gauge bosons are neglected here.
However, we would like to mention that these mixings can be generated at one
loop level and we have considered the effect of loop induced kinetic mixing in
our discussion on the constraints of dark gauge boson.

The Lagrangian of our model is given by Tapadar:2021kgw

$\displaystyle\mathcal{L}$ $\displaystyle=$ $\displaystyle\mathcal{L}_{\rm
SM}+\bar{\chi}\left(i\not{\partial}-m_{\chi}\right)\chi-\dfrac{1}{4}\hat{X}^{\rho\sigma}\hat{X}_{\rho\sigma}-\dfrac{1}{4}\hat{F}_{\mu\tau}^{\rho\sigma}\hat{F}_{\mu\tau_{\rho\sigma}}+\dfrac{\sin\delta}{2}\hat{F}_{\mu\tau}^{\rho\sigma}\hat{X}_{\rho\sigma}$
(1) $\displaystyle-
g_{X}\bar{\chi}\gamma^{\rho}\chi{\hat{Z^{\prime}}}_{\rho}-g_{\mu\tau}\left(\bar{\mu}\gamma_{\rho}\mu+\bar{\nu}_{\mu}\gamma_{\rho}P_{L}\nu_{\mu}-\bar{\tau}\gamma_{\rho}\tau-\bar{\nu}_{\tau}\gamma_{\rho}P_{L}\nu_{\tau}\right)\hat{Z}_{\mu\tau}^{\rho}$
$\displaystyle+\dfrac{1}{2}\hat{m}_{\mu\tau}^{2}\hat{Z}_{\mu\tau}^{\rho}\hat{Z}_{\mu\tau_{\rho}}+\dfrac{1}{2}\hat{m}^{\prime
2}\hat{Z^{\prime}}^{\rho}\hat{Z^{\prime}}_{\rho}\,\,.$

Here,
$\hat{X}^{\rho\sigma}=\partial^{\rho}\hat{Z^{\prime}}^{\sigma}-\partial^{\sigma}\hat{Z}^{\prime}{}^{\rho}$
and
$\hat{F}_{\mu\tau}^{\rho\sigma}=\partial^{\rho}\hat{Z}_{\mu\tau}^{\sigma}-\partial^{\sigma}\hat{Z}_{\mu\tau}^{\rho}$
are the field strength tensors of $U(1)_{X}$ and $U(1)_{L_{\mu}-L_{\tau}}$
gauge symmetry respectively. $m_{\chi}$ is the mass of the DM $\chi$,
$\hat{m}^{\prime}$ and $\hat{m}_{\mu\tau}$ are the gauge bosons’ mass
parameters of the theory. The gauge couplings corresponding to $U(1)_{X}$ and
$U(1)_{L_{\mu}-L_{\tau}}$ are denoted by $g_{X}$ and $g_{\mu\tau}$
respectively. The last term in the first line of Eq. 1 denotes the tree level
kinetic mixing between $\hat{Z}^{\prime}$ and $\hat{Z}_{\mu\tau}$.

To put the Lagrangian in canonical form, we have performed a non-orthogonal
transformation from “hatted” to “tilde” basis and the transformation is

$\displaystyle\begin{pmatrix}\hat{Z}^{\prime}_{\rho}\\\
\hat{Z}_{\mu\tau_{\rho}}\end{pmatrix}=\begin{pmatrix}\sec\delta&0\\\
\tan\delta&1\end{pmatrix}\begin{pmatrix}\tilde{Z^{\prime}}_{\rho}\\\
\tilde{Z}_{\mu\tau_{\rho}}\end{pmatrix}\,\,.$ (2)

Removal of kinetic mixing term generates mass mixing between
$\tilde{Z^{\prime}}$ and $\tilde{Z}_{\mu\tau}$ bosons. This can be removed by
transforming “tilde” basis into mass basis by the following orthogonal
transformation

$\displaystyle\begin{pmatrix}\tilde{Z}_{\mu\tau_{\rho}}\\\
\tilde{Z}^{\prime}_{\rho}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\cos\beta&-\sin\beta\\\
\sin\beta&\cos\beta\end{pmatrix}\begin{pmatrix}Z_{\mu\tau_{\rho}}\\\
Z^{\prime}_{\rho}\end{pmatrix}\,\,\,,$ (3)

where $\beta$ is the $Z^{\prime}-Z_{\mu\tau}$ mixing angle and in the limit
$\delta\ll 1$, it is given by

$\displaystyle\tan\beta$ $\displaystyle\simeq$
$\displaystyle\dfrac{\tan\delta}{1-\hat{r}^{2}}\,\,,\ $ (4)

where $\hat{r}=m^{\prime 2}/m_{\mu\tau}^{2}$, $m^{\prime}$ and $m_{\mu\tau}$
are the masses of $Z^{\prime}$ and $Z_{\mu\tau}$ respectively. In the limit
$\delta,\,\hat{r}\ll 1$, $\hat{Z}$ and $\hat{Z}_{\mu\tau}$ can be expressed in
the following manner

$\displaystyle\hat{Z}_{\mu\tau_{\rho}}$ $\displaystyle\simeq$ $\displaystyle
Z_{\mu\tau_{\rho}}-\hat{\epsilon}Z^{\prime}_{\rho}\,\,,$
$\displaystyle\hat{Z}^{\prime}_{\rho}$ $\displaystyle\simeq$ $\displaystyle
Z^{\prime}_{\rho}+\tan\delta Z_{\mu\tau_{\rho}}\,\,,$ (5)

where $\hat{\epsilon}=\hat{r}^{2}\rm{tan}\,\delta$. Using Eq. 5 we can rewrite
the interaction term of dark vector boson $Z^{\prime}$ as follows

$\displaystyle\mathcal{L}\supset-g_{X}\bar{\chi}\gamma^{\rho}\chi
Z^{\prime}_{\rho}+\epsilon\left(\bar{\mu}\gamma^{\rho}\mu+\bar{\nu}_{\mu}\gamma^{\rho}P_{L}\nu_{\mu}-\bar{\tau}\gamma^{\rho}\tau-\bar{\nu}_{\tau}\gamma^{\rho}P_{L}\nu_{\tau}\right)Z^{\prime}_{\rho}\,\,\,.$
(6)

Here, $\epsilon=g_{\mu\tau}\hat{\epsilon}$ is the portal coupling between dark
and the visible sector. Note that the condition $\hat{r}\ll 1$ tells us that
$Z_{\mu\tau}$ will have negligible contribution to muon $(g-2)$ anomaly and we
have only considered the contribution of $Z^{\prime}$ to the muon $(g-2)$
anomaly in the rest of our discussion.

##  3 Fast expanding Universe

###  3.1 Cosmology of fast expanding Universe

Prediction of light element abundance at the time of BBN tells us that our
Universe was radiation dominated at the time of BBN. However, it is possible
that the Universe was not radiation dominated before BBN. This implies that a
new species $\phi$ might have existed at the pre-BBN era whose energy density
redshifts faster than radiation DEramo:2017gpl ; DEramo:2017ecx . In that
case, total energy density of the Universe can be expressed as

$\displaystyle\rho(T)$ $\displaystyle=$
$\displaystyle\rho_{\phi}(T)+\rho_{r}(T)+\rho_{r}^{\prime}(T^{\prime})\,\,,$
(7)

where $\rho_{r}(T)=\dfrac{\pi^{2}}{30}g_{\rho}(T)T^{4}$ is the radiation
energy density of the SM bath at SM temperature $T$,
$\rho_{r}^{\prime}(T^{\prime})=\dfrac{\pi^{2}}{30}g_{\rho^{\prime}}(T^{\prime}){T^{\prime}}^{4}$
is the radiation energy density of dark sector (DS) and $T^{\prime}$ is the
temperature of the DS. $\rho_{\phi}(T)$ is the energy density of the new
species $\phi$ and it redshifts as $a^{-(4+n)}$ where $a$ is the scale factor
and $n>0$. Here, $g_{\rho}(T)$ and $g_{\rho^{\prime}}(T^{\prime})$ are
effective number of relativistic degrees of freedom of SM sector and DS
respectively, contributing to the energy density. Here, we have assumed
$T^{\prime}\ll T$ and we may neglect DS sector energy density compared to SM
bath. Therefore the total energy density can be expressed as

$\displaystyle\rho(T)$ $\displaystyle=$
$\displaystyle\rho_{\phi}(T)+\rho_{r}(T).$ (8)

Motivated by the constraints from BBN, one can define a temperature $T_{r}$ at
which111The abundances of the light elements at the time of BBN remain
unaltered if we choose $T_{r}\gtrsim(15.4)^{1/n}\rm MeV$ where $n>0$
DEramo:2017gpl . $\rho_{\phi}(T_{r})=\rho_{r}(T_{r})$. Using this definition
of $T_{r}$ and entropy conservation of the Universe, we can write
$\rho_{\phi}(T)$ as

$\displaystyle\rho_{\phi}(T)$ $\displaystyle=$
$\displaystyle\rho_{r}(T)\left(\dfrac{g_{\rho}(T_{r})}{g_{\rho}(T)}\right)\left(\dfrac{g_{*s}(T)}{g_{*s}(T_{r})}\right)^{\frac{4+n}{3}}\left(\dfrac{T}{T_{r}}\right)^{n}\,\,.\,$
(9)

Therefore, total energy density can be written as

$\displaystyle\rho(T)$ $\displaystyle=$
$\displaystyle\rho_{r}(T)\left[1+\left(\dfrac{g_{\rho}(T_{r})}{g_{\rho}(T)}\right)\left(\dfrac{g_{*s}(T)}{g_{*s}(T_{r})}\right)^{\frac{4+n}{3}}\left(\dfrac{T}{T_{r}}\right)^{n}\right]\,\,.\,$
(10)

It is clear from Eq. 10 that at early times, i.e. $T\gg T_{r}$, energy density
of $\phi$ dominates the energy budget of the Universe and the Universe expands
faster than radiation dominated era. Similarly, for $T\ll T_{r}$, one can see
$\rho(T)\simeq\rho_{r}(T)$. We would like to mention that throughout our
analysis we have chosen $T_{r}=20\rm MeV$ and the choice of $T_{r}$ is
consistent with the BBN constraint as discussed earlier.

Now, in presence of a new species $\phi$, Hubble parameter is given by

$\displaystyle H(T)$ $\displaystyle=$ $\displaystyle\sqrt{\dfrac{8\pi}{3M_{\rm
Pl}^{2}}\rho(T)}\,\,,$ (11)

where $M_{\rm{Pl}}=1.22\times 10^{19}\rm GeV$ is the Planck mass. At high
temperature limit i.e. $T\gg T_{r}$, $H(T)$ can be approximated as

$\displaystyle H(T)$ $\displaystyle\simeq$ $\displaystyle\dfrac{\pi}{3M_{\rm
Pl}}\sqrt{\dfrac{4\pi}{5}}\sqrt{g_{\rho}(T)}T^{2}\left(\dfrac{T}{T_{r}}\right)^{n/2}\,\,.\
$ (12)

###  3.2 Thermalisation of visible sector and dark sector

Figure 1:  Allowed regions in $m_{\chi}-\epsilon$ plane from the
thermalisation criterion for $n=1$ (red), $n=2$ (violet), $n=3$ (green), $n=4$
(yellow) are shown. The constraint considering radiation dominated Universe is
shown with cyan color.

Thermalisation of DS and SM bath not only depends on the interaction strength
of the portal coupling, but also on the expansion rate of the Universe. In
standard radiation dominated Universe for large portal coupling, there would
be sufficient amount of energy exchange between the two sectors and they will
be in thermal equilibrium with each other. In this framework DM freezes-out as
the universe expands and DM relic abundance is determined by its annihilation
into dark vector boson, $\chi\bar{\chi}\to Z^{\prime}Z^{\prime}$ if
$m_{\chi}>m_{Z^{\prime}}$. This is known as secluded sector freeze out
Pospelov:2007mp ; Feng:2008mu ; Cheung:2010gj ; Chu:2011be ; Evans:2017kti .

In our study, the central idea of DM number density evolution is focused on
the fact that two sectors are not in thermal equilibrium with each other.
Therefore, first we identify the parameter space in $m_{\chi}-\epsilon$ plane
in which the dark sector is not in thermal equilibrium with the visible
sector. To identify this region, initially we assume that the two sectors are
in thermal equilibrium. We calculate reaction rate $\left(\Gamma\right)$ of
all processes producing $Z^{\prime}$ from SM bath (the processes are mentioned
in section 3.3). If $\Gamma(T)_{{\rm SM}\rightarrow Z^{\prime}}\gtrsim H(T)$
at the time of DM freeze-out, then two sectors are in thermal contact with
each other, at the time of DM freeze-out. This region is shown in Fig.1 for
different values of expansion parameter $n$ which controls the expansion rate
of the Universe.

In this scenario, at early Universe in the absence of any mass thresholds,
reaction rates go as $\Gamma_{{\rm SM}\rightarrow
Z^{\prime}}\sim\epsilon^{2}T$. Now, in case of radiation dominated Universe,
the ratio $\Gamma_{{\rm SM}\rightarrow Z^{\prime}}/H(T)$, at the time of DM
freeze-out, goes as $\sim\epsilon^{2}/m_{\chi}$. Thus we need large value of
$\epsilon$ as we increase the value of $m_{\chi}$ for the dark sector to be
thermalised with the SM bath and we can see this from Fig. 1. However in case
of fast expanding Universe, $H(T)$ is much larger in comparison to the
radiation dominated Universe as one can see from Eq. 12. As a result, the
required value of $\epsilon$ to establish thermal equilibrium between dark and
visible sectors will be larger in comparison to the radiation dominated
Universe. Thus for a fixed value of $m_{\chi}$, the upper limit of $\epsilon$
for thermally decoupled dark sector increases as we increase the value of $n$.
In Fig. 1, we show the parameter space in $m_{\chi}-\epsilon$ plane for
thermally decoupled dark sector and one can clearly see from the figure that
the required value of $\epsilon$ for the thermalisation of dark sector
increases with the increase in the value of $n$.

Modification of the cosmological history of the early Universe opens up the
possibility of explaining muon $(g-2)$ anomaly while the dark sector is still
thermally decoupled. In particular, for $\epsilon\sim 10^{-3}$, the dark
sector is still thermally decoupled if $n\geq 3$ and $m_{\chi}\gtrsim 30\rm
TeV$. However, in this region of parameter space, the dark gauge boson
$Z^{\prime}$ can elevate the theoretical prediction of muon $(g-2)$ if
$m_{Z^{\prime}}\sim\mathcal{O}(100\rm MeV)$. In the subsequent sections, we
will discuss the DM production mechanisms for the thermally decoupled dark
sector in fast expanding Universe.

###  3.3 Temperature evolution of dark sector

To study the temperature evolution of the DS, our next step is to calculate DS
temperature $T^{\prime}$. In our scenario, dark sector is populated from $\mu$
and $\tau$ flavored SM leptons involving annihilations processes in the early
Universe and the DS temperature depends on interaction strength of the portal
coupling $(\epsilon)$ between DS and VS, $T_{r}$, and $n$. To calculate
$T^{\prime}$, we have considered all possible $Z^{\prime}$ production
processes from SM bath such as
$l\left(\bar{l}\right)+W^{+}\left(W^{-}\right)\rightarrow\nu_{l}\left(\bar{\nu_{l}}\right)+Z^{\prime}$,
$\nu_{l}\left(\bar{\nu_{l}}\right)+W^{-}\left(W^{+}\right)\rightarrow
l\left(\bar{l}\right)+Z^{\prime}$,
$l\left(\bar{l}\right)+\nu_{l}\left(\bar{\nu_{l}}\right)\rightarrow
W^{-}\left(W^{+}\right)+Z^{\prime}$,
$l+\bar{l}\rightarrow\gamma\left(Z,h\right)+Z^{\prime}$,
$\nu_{l}+\bar{\nu_{l}}\rightarrow Z\left(h\right)+Z^{\prime}$,
$l+\gamma\left(Z,h\right)\rightarrow l+Z^{\prime}$,
$\bar{l}+\gamma\left(Z,h\right)\rightarrow\bar{l}+Z^{\prime}$,
$\nu_{l}+Z\left(h\right)\rightarrow\nu_{l}+Z^{\prime}$,
$\bar{\nu_{l}}+Z\left(h\right)\rightarrow\bar{\nu_{l}}+Z^{\prime}$,
$l+\bar{l}\rightarrow Z^{\prime}$, $\nu_{l}+\bar{\nu_{l}}\rightarrow
Z^{\prime}$ where, $l=\\{\mu,\tau\\}$. Now the DS temperature evolution can be
obtained from the Boltzmann equation for the DS energy density
$\rho^{\prime}$, and it is given by,

$\displaystyle\dfrac{d\rho^{\prime}}{dt}+4H\rho^{\prime}$
$\displaystyle\simeq$ $\displaystyle\mathcal{C}_{{\rm SM}\to
Z^{\prime}}\left(T\right)\,\,,\ $ (13)

where $\mathcal{C}_{{\rm SM}\to Z^{\prime}}$ is the collision term for above-
mentioned processes. In our calculation of the collision term, we have only
considered the processes which are proportional to $\epsilon$ at the amplitude
level. Since $T^{\prime}\ll T$, energy exchange from DS to VS is very small,
and therefore we have not considered the collision term due to dark sector
annihilation into SM. Now, we have defined $\xi$ as
$\xi\equiv\dfrac{T^{\prime}}{T}$ and using Eq. 13, the solution of $\xi(T)$ as
a function of $T$ is as follows

$\displaystyle\xi(T)$ $\displaystyle=$
$\displaystyle\left[\int_{T}^{T_{0}}\dfrac{30\,\mathcal{C}_{{\rm SM}\to
Z^{\prime}}(\tilde{T})}{g_{\rho^{\prime}}\pi^{2}H(\tilde{T})\tilde{T}^{5}}d\tilde{T}\right]^{1/4}\,\,.\,$
(14)

Here, $T_{0}$ is the initial temperature of the Universe and we take
$T^{\prime}\simeq 0$ at $T=T_{0}$, $g_{\rho^{\prime}}$ is the relativistic
degrees of freedom contributing to the radiation bath of the dark sector,
$H(T)$ is the expansion rate of the Universe as given in Eq. 11.

Figure 2:  Evolution of $\xi(T)$ (solid lines) as a function of $T$ for
$n=1$(black) for $n=2$ (red), $n=3$ (green), $n=4$ (blue). The dashed lines
with same color codes denote the behavior of semi-analytic expression of
$\xi(T)$ as a function of $T$ (given in Eq. 15). Here, we have considered
$m_{Z^{\prime}}=100\rm MeV$ and $\epsilon=10^{-10}$.

In Fig. 2, we show the variation of $\xi(T)$ as a function of $T$ for
$m_{Z^{\prime}}=100\rm MeV$ and $\epsilon=10^{-10}$ for four different values
of $n$. One can clearly see that at $T\sim m_{Z^{\prime}}$ energy injection
stops and after that $\xi(T)$ remains constant. Another important feature is
that the value of $\xi(T)$ at a fixed value of $T$ decreases as we increase
$n$. This is because with the increase in $n$, the energy injection processes
decouple earlier and as a result amount of energy injection will be less.

In our scenario, the collision terms are proportional to the portal coupling
$\epsilon^{2}$ and the portal interaction is renormalisable. Since the portal
interaction is renormalisable, using Eq. 14 one can write a semi-analytic
expression for $\xi$ as follows Evans:2019vxr ; Tapadar:2021kgw

$\displaystyle\xi$ $\displaystyle\simeq$
$\displaystyle\dfrac{\zeta(m_{Z^{\prime}})\sqrt{\epsilon}}{(1+\frac{n}{2})^{\frac{1}{4}}}\,T^{-\frac{1}{4}}\left(\dfrac{T_{r}}{T}\right)^{\frac{n}{8}}\,\,{\rm
for}\,\,T>T_{r}$ (15) $\displaystyle\simeq$
$\displaystyle\zeta(m_{Z^{\prime}})\sqrt{\epsilon}T^{-\frac{1}{4}}\,\,\,{\rm
for}\,\,T\leq T_{r},$

where $\zeta(m_{Z^{\prime}})$ will be extracted from our numerical analysis.
Since the effect of non-standard cosmology will enter into the evolution of
$\xi$ only through Eq. 11, therefore $\zeta(m_{Z^{\prime}})$ does not depend
on $n$. In Fig. 3, we show the variation of $\zeta(m_{Z^{\prime}})$ for
different values of $m_{Z^{\prime}}$ and one can see that for $T\gtrsim
10m_{Z^{\prime}}$, $\zeta(m_{Z^{\prime}})\simeq 10^{4}$. Thus during the phase
of energy injection, one can use Eq. 15 for the evolution of $T^{\prime}$. The
dashed lines of Fig. 2 denote the behaviour of Eq. 15 as a function of $T$ for
different values of $n$.

Figure 3:  Variation of $\zeta(m_{Z^{\prime}})$ as a function $T$ for
$m_{Z^{\prime}}=100\rm MeV$ (blue dotted line), $1\rm GeV$ (black dashed
line), and $100\rm GeV$ (red solid line). Here we have considered
$\epsilon=10^{-10}$ and $n=1$.

##  4 Properties of non-adiabatic evolution in fast expanding Universe

Dark matter freeze-out through non-adiabatic evolution in fast expanding
Universe is quite different from either of non-adiabatic evolution scenario as
discussed earlier in Evans:2017kti ; Tapadar:2021kgw  or dark matter freeze-
out in non-standard cosmological scenario DEramo:2017gpl . We shall start with
some analytical calculation of DM freeze-out in this scenario and also put
some bounds on portal coupling.

Dark matter freeze-out within DS through non-adiabatic evolution is referred
to as “leak-in” dark matter Evans:2019vxr . Difference between non-adiabatic
evolution of DM number density with dark sector freeze-out is that in the
former scenario, $\xi$ is a dynamical quantity which depends on the portal
coupling as well as the SM temperature whereas $\xi$ is constant in the latter
scenario. Now, the relevant Boltzmann equation for dark sector freeze-out is
given by

$\displaystyle\dfrac{dY_{\chi}}{dx}$ $\displaystyle=$
$\displaystyle\dfrac{\langle\sigma
v\rangle^{T^{\prime}}s(x)}{xH(x)}\left({Y_{\chi}^{\rm
eq}}(x^{\prime})^{2}-Y_{\chi}^{2}\right)\,\,,\ $ (16)

where, $x=m_{\chi}/T$, $Y_{\chi}^{\rm eq}(x^{\prime})=n_{\chi}^{\rm
eq}(x^{\prime})/s(x)$, where $n_{\chi}^{\rm eq}(x^{\prime})$ is the
equilibrium DM number density, $s(x)$ is the entropy density of the Universe,
and $x^{\prime}=m_{\chi}/T^{\prime}$. $\langle\sigma v\rangle^{T^{\prime}}$ is
the thermally averaged cross-section of DM annihilation, calculated at
$T^{\prime}$ and $H(x)$ is the Hubble parameter which is given in Eq. 11.

Initially, DM comoving number density follows $Y_{\chi}^{\rm eq}(x^{\prime})$
and at later times it decouples from dark sector thermal bath and freezes-out
at $x^{\prime}_{f}=\dfrac{m_{\chi}}{T^{\prime}_{f}}$ where, $T^{\prime}_{f}$
is the dark sector temperature at the time of DM freeze-out. In case of fast
expanding Universe, due to the higher expansion rate, number changing
processes for DM decouple much earlier in comparison to the radiation
dominated Universe. Therefore, in this scenario $x_{f}^{\prime}$ appeared to
be smaller than the radiation dominated Universe. Exact value of
$x_{f}^{\prime}$ can be obtained by solving the following equation
DEramo:2017gpl

$\displaystyle\dfrac{e^{x^{\prime}_{f}}{\sqrt{x_{f}^{\prime}}}}{\xi\langle\sigma
v\rangle}\left(\dfrac{x_{r}}{\xi x^{\prime}_{f}}\right)^{\frac{n}{2}}$
$\displaystyle=$
$\displaystyle\dfrac{c\left(c+2\right)}{\left(c+1\right)}\dfrac{3\sqrt{5}}{4\sqrt{2}\pi^{3}}\dfrac{g_{\chi}}{\sqrt{g_{*s}}}m_{\chi}M_{\rm{Pl}}\,\,,\
$ (17)

where $g_{\chi}$ is the degrees of freedom of DM, which is $2$ in our
scenario, $x_{r}=m_{\chi}/T_{r}$, $g_{*s}$ is effective degrees of freedom
contributing to entropy density, and $c$ is a constant which is numerically
found to be equal to $c(c+2)=\left(\frac{n}{2}+1\right)$ Scherrer:1985zt . For
$s$-wave DM annihilation cross-section Eq. 17 can be approximated as,

$\displaystyle x^{\prime}_{f}\simeq\ln A-\dfrac{n}{2}\ln x_{r}\,\,,\ $ (18)

where
$A=\dfrac{c\left(c+2\right)}{\left(c+1\right)}\dfrac{3\sqrt{5}}{4\sqrt{2}\pi^{3}}\dfrac{g_{\chi}}{\sqrt{g_{*s}}}\langle\sigma
v\rangle m_{\chi}M_{\rm Pl}$.

To get an analytical expression of final relic abundance $Y_{\infty}$, we have
considered $s$-wave dominated annihilation cross section. As DM freeze-out
temperature from DS thermal bath is $T^{\prime}_{f}$ then SM bath temperature
at the time of DM freeze-out is $T^{\prime}_{f}/\xi$. At late time after
freeze-out $Y_{\chi}^{\rm eq}(T^{\prime})\simeq 0$ compared to $Y_{\infty}$
which is the final value of DM comoving number density. Therefore, final
comoving number density of DM is given by,

$\displaystyle\dfrac{1}{Y_{\infty}}$ $\displaystyle=$
$\displaystyle\dfrac{1}{Y_{f}\left(\xi x^{\prime}_{f}\right)}+\langle\sigma
v\rangle_{0}\int_{\xi
x^{\prime}_{f}}^{x_{r}}\dfrac{s(x)}{xH(x)}dx+\int_{x_{r}}^{\infty}\dfrac{s(x)}{xH_{\rm{rad}}(x)}dx\,\,,\
$ (19)

where $\langle\sigma v\rangle_{0}$ is the s-wave thermally averaged DM
annihilation cross section, $H(x)$ is the Hubble parameter in fast expanding
Universe as approximated in Eq. 12. For $T<T_{r}$, Universe is radiation
dominated and the Hubble parameter in this regime is defined by
$H_{\rm{rad}}(x)=\dfrac{1}{M_{\rm Pl}}\sqrt{\dfrac{8\pi}{3}\rho_{r}(T)}$.

Therefore, by integrating Eq. 19, approximated expression for the comoving
number density is given by,

$\displaystyle Y_{\infty}$ $\displaystyle=$
$\displaystyle\sqrt{\dfrac{45}{\pi}}\dfrac{\sqrt{g_{\rho}(x)}}{g_{*s}(x)}\dfrac{2x_{r}}{\langle\sigma
v\rangle m_{\chi}M_{\rm Pl}}\left[\dfrac{2}{\xi
x^{\prime}_{f}}\left(\dfrac{\xi
x^{\prime}_{f}}{x_{r}}\right)^{\frac{n}{2}-1}+\dfrac{1}{n-2}\left(1-\left(\dfrac{\xi
x^{\prime}_{f}}{x_{r}}\right)^{\frac{n}{2}-1}\right)+1\right]^{-1}\,,\ n\neq
2$ (21) $\displaystyle=$
$\displaystyle\sqrt{\dfrac{45}{\pi}}\dfrac{\sqrt{g_{\rho}(x)}}{g_{*s}(x)}\dfrac{x_{r}}{\langle\sigma
v\rangle m_{\chi}{M_{\rm Pl}}}\left[\dfrac{2}{\xi x^{\prime}_{f}}+{\rm
ln}\left(\frac{x_{r}}{\xi x^{\prime}_{f}}\right)+1\right]^{-1}\,,\ n=2\,\,.\ $

Eq. 21 reduces to standard WIMP scenario as soon as we put $\xi=1$ and $n=0$.
Thus, it is clear that apart from the interactions within DS, final relic
abundance for internally thermalized DS in non-standard cosmology also depends
on $n$, portal interaction between DS and VS and $T_{r}$.

DM production through this mechanism can take place within a certain range of
portal coupling. Leak-in mechanism is only possible when dark and visible
sectors are not in thermal equilibrium, but there is small exchange of energy
between the two sectors. If portal coupling is sufficiently large then both
the sectors equilibrate and hence no leak-in mechanism is observed. On the
other hand, if portal coupling is very small then there would never be
sufficient DM produced at early Universe that could reach observed relic
abundance of DM after freeze-out.

Analytically, we can obtain a lower limit on portal coupling by comparing
observed relic density with the maximum relic abundance of DM in leak-in
scenario. To do this, we first need to find the maximum value attained by
equilibrium yield,

$\displaystyle Y_{\chi}^{\rm
eq}=\dfrac{g_{\chi}}{g_{*s}T^{3}}\dfrac{45}{2\pi^{2}}\left(\dfrac{m_{\chi}T^{\prime}}{2\pi}\right)^{3/2}e^{-m_{\chi}/T^{\prime}}\,\,.\
$ (22)

We maximize $Y_{\chi}^{\rm eq}$ with respect to $T^{\prime}$ to obtain the
dark sector temperature $T^{\prime}_{\rm max}$ at which the equilibrium yield
is maximum. The value of $T^{\prime}_{\rm max}$ is given by

$\displaystyle T^{\prime}_{\rm
max}=\dfrac{2m_{\chi}\left(6-n\right)}{3n+30}\,\,.\ $ (23)

One can notice that Eq. 23 reduces to standard result of $T^{\prime}_{\rm
max}=\frac{2m_{\chi}}{5}$ Evans:2019vxr  as we choose $n=0$. Another point to
note is that as $n$ increases, $T^{\prime}_{\rm max}$ decreases. Using Eq. 23,
the equilibrium yield at $T^{\prime}_{\rm max}$ is given by,

$\displaystyle{Y_{\chi}^{\rm eq}}_{\rm max}$ $\displaystyle=$
$\displaystyle\dfrac{g_{\chi}}{g_{*s}}\dfrac{45}{2\pi^{2}}\left(\dfrac{m_{\chi}}{2\pi}\right)^{3/2}e^{-m_{\chi}/T^{\prime}_{\rm
max}}\dfrac{{T^{\prime}_{\rm max}}^{-(30+3n)/2(6-n)}}{K^{-24/(6-n)}}\,\,,\ $
(24)

where
$K=\frac{10^{4}\sqrt{\epsilon}}{\left(1+\frac{n}{2}\right)^{1/4}}T_{r}^{n/8}$.

Now, in case of fast expanding Universe, the DM can be produced via leak-in
mechanism if ${Y_{\chi}^{\rm eq}}_{\rm max}\geq Y_{\rm obs}$ where $Y_{\rm{\rm
obs}}=4.355\times 10^{-10}\left(1{\rm GeV}/m_{\chi}\right)$. From this
condition, one can derive an upper limit on $\epsilon$ and it is given by

$\displaystyle\epsilon\gtrsim 2.616\times 10^{-12}\left(3.561\times
10^{8}\right)^{n/12}e^{n/8}\sqrt{\left(2+n\right)}\,\
\left(\dfrac{g_{\chi}}{g_{*s}}\right)^{\dfrac{n-6}{12}}\left(\dfrac{2\left(6-n\right)}{3\left(10+n\right)}\right)^{\dfrac{10+n}{8}}T_{r}^{-n/4}m_{\chi}^{n/3}\,\,.\
$ (25)

Thus, Eq. 25 tells us that unlike DM freeze-out in non-adiabatic scenario in
standard cosmology, lower bound on $\epsilon$ is proportional to $m_{\chi}$.
Thus, for leak-in mechanism, we need larger value of $\epsilon$ as we increase
the mass of the DM. The upper limit given in Eq. 25 is known as “Absolute
coupling floor” below which it is not possible to get correct relic density
via leak-in mechanism.

##  5 Numerical result for relic density

In this section, we discuss the evolution of DM comoving number density in
fast expanding Universe where visible sector and dark sector are not in
thermal equilibrium. The temperature evolution of DS is different from the
evolution of $T$ and we use the semi-analytic expression of $\xi$, given in
Eq. 15 to calculate $T^{\prime}$.

###  5.1 Boltzmann equation

Evolution of the DM number density depends on the production of DM from SM
fields and annihilation of DM into $Z^{\prime}$. Relevant Boltzmann equation
for the evolution of total number density of DM is given by,

$\displaystyle\dfrac{dY_{\chi_{\rm tot}}}{dx}$ $\displaystyle=$
$\displaystyle\dfrac{s(x)h_{\rm
eff}(x)}{H(x)x}\left[\dfrac{1}{2}\left[\langle\sigma
v\rangle^{T^{\prime}}_{\chi\bar{\chi}\to
Z^{\prime}Z^{\prime}}\left(Y_{\chi_{\rm tot}}^{\rm
eq}(x,x^{\prime})^{2}-Y_{\chi_{\rm
tot}}^{2}\right)\right]+2\sum_{f}\langle\sigma
v\rangle^{T}_{{f}\bar{f}\to{\chi}\bar{\chi}}Y_{f}^{\rm eq}(x)^{2}\right]\,,$

where $Y_{\chi_{\rm tot}}^{\rm eq}(x,x^{\prime})=\dfrac{n_{\chi_{\rm
tot}}^{\rm eq}(x^{\prime})}{s(x)}$ is the total equilibrium comoving number
density of DM where, $n_{\chi_{\rm tot}}^{\rm eq}(x^{\prime})=n_{\chi}^{\rm
eq}(x^{\prime})+n_{\bar{\chi}}^{\rm eq}(x^{\prime})$ is the total equilibrium
DM number density at $T^{\prime}$. $Y_{f}^{\rm eq}(x)=\dfrac{n_{f}^{\rm
eq}(x)}{s(x)}$ where, $n_{f}^{\rm eq}(x)$ is the equilibrium number density of
SM fermion $f$ at SM temperature $T$, $h_{\rm
eff}(x)=\left(1-\dfrac{1}{3}\dfrac{d\ln g_{*s}(x)}{d\ln x}\right)$, and $H(x)$
is the Hubble parameter, defined in Eq. 11.

The first term inside the square bracket on the right hand side (RHS) of Eq.
5.1 denotes the contribution to the change in $Y_{\chi_{\rm tot}}$ due to the
annihilation of DM into $Z^{\prime}$ and we call this term as “annihilation
term” in the rest of our discussion. $\langle\sigma
v\rangle^{T^{\prime}}_{\chi\bar{\chi}\to Z^{\prime}Z^{\prime}}$ is the
thermally averaged cross section of $\chi\bar{\chi}\to Z^{\prime}Z^{\prime}$,
calculated at $T^{\prime}$ and the $s$-wave approximated form of
$\langle\sigma v\rangle^{T^{\prime}}_{\chi\bar{\chi}\to Z^{\prime}Z^{\prime}}$
is given by

$\displaystyle\langle\sigma v\rangle_{{\chi}\bar{\chi}\to
Z^{\prime}Z^{\prime}}$ $\displaystyle\simeq$
$\displaystyle\dfrac{4\pi\alpha_{X}^{2}}{m_{\chi}^{2}}\dfrac{m_{\chi}\left(m_{\chi}^{2}-m_{Z^{\prime}}^{2}\right)^{3/2}}{\left(2m_{\chi}^{2}-m_{Z^{\prime}}^{2}\right)^{2}}\,\,.$
(27)

The second term inside the square bracket on the RHS of Eq. 5.1 denotes the
production of DM from the SM bath and we call this term as “freeze-in” term in
the rest of our discussion. Here, $\langle\sigma
v\rangle^{T}_{{f}\bar{f}\to{\chi}\bar{\chi}}$ is the thermally averaged cross
section of $f\bar{f}\to\chi\bar{\chi}$, calculated at temperature $T$. The
cross sections of the relevant processes for the freeze-in are given by

$\displaystyle\sigma_{{f}\bar{f}\to{\chi}\bar{\chi}}$ $\displaystyle=$
$\displaystyle\dfrac{\kappa\alpha_{X}\epsilon^{2}}{3\hat{s}}\sqrt{\dfrac{\hat{s}-4m_{\chi}^{2}}{\hat{s}-4m_{f}^{2}}}\dfrac{(\hat{s}+2m_{f}^{2})(\hat{s}+2m_{\chi}^{2})}{(\hat{s}-m_{Z^{\prime}}^{2})^{2}}\,\,,$
(28)

where $m_{f}$ is the mass of the SM fermion $f$, $\kappa=1$ for $f=\mu,\,\tau$
and $\kappa=1/2$ for $f=\nu_{\mu},\,\nu_{\tau}$.

In the following section, we discuss the effect of annihilation and freeze-in
terms on the DM relic density in the fast expanding Universe. The numerical
results related to the DM relic density has also been discussed in the next
section.

###  5.2 Outcome of numerical analysis

In this section, we present our numerical analysis by solving Boltzmann
equation given in Eq. 5.1 in the presence of fast expanding Universe. In our
analysis, we use the semi-analytic expression of $\xi(T)$, given in Eq. 15, to
calculate $T^{\prime}$. Here, we have considered $T_{r}=20\rm MeV$ as
discussed earlier and $m_{Z^{\prime}}=100\rm MeV$. Let us note that the choice
of $m_{Z^{\prime}}$ is motivated by the explanation of muon $(g-2)$ anomaly.
However, as far as the relic density analysis is concerned, it is independent
on the choice of $m_{Z^{\prime}}$ for $m_{\chi}/m_{Z^{\prime}}\gtrsim 10$.

Figure 4: Evolution of relative abundance of DM as a function of $x$ for
$n=1$, $\epsilon=3\times 10^{-9}$, and $m_{Z^{\prime}}=100\rm MeV$. Left
panel: Here we consider $m_{\chi}=10^{3}\rm GeV$ and for
$\alpha_{X}=2.35\times 10^{-5}$ (red and blue solid line) and
$\alpha_{X}=0.144$ (black solid line), DM is produced via freeze-in and
reannihilation mechanisms respectively. For $\alpha_{X}=5\times 10^{-3}$
(brown dot-dashed line), final abundance of DM is determined by the leak-in
mechanism but it overcloses the Universe. Right panel: Here, we consider
$m_{\chi}=10^{5}\rm GeV$. For $\alpha_{X}=2.24\times 10^{-2}$ DM is produced
by freeze-in mechanism and it is shown by blue solid line. For $\alpha_{X}=1$,
(black solid line) and $4\pi$ (red solid line), the amount of produced DM
overcloses the Universe. In both panels, black dotted and black dashed lines
denote the evolution of $\Omega h^{2}_{\rm eq}/\Omega h^{2}_{\rm obs}$ and
$\Omega h^{2}_{\rm QSE}/\Omega h^{2}_{\rm obs}$ respectively.

In Fig. 4, we have shown the variation of the relative abundance of DM
(defined as $\Omega h^{2}/\Omega h^{2}_{\rm obs}$, where $\Omega
h^{2}=2.755\times 10^{8}(m_{\chi}/1\rm GeV)Y_{\chi_{\rm tot}}$ and $\Omega
h^{2}_{\rm obs}=0.12$) as a function of $x$ for $n=1$, $\epsilon=3\times
10^{-9}$ and two different values of $m_{\chi}$. In the left panel of Fig. 4,
we have considered $m_{\chi}=10^{3}\rm GeV$. Here, one can see that for three
different values of $\alpha_{X}$, DM production mechanisms are different. When
$\alpha_{X}$ is very small, DS is not internally thermalized and DM can be
produced by both $Z^{\prime}Z^{\prime}\rightarrow\chi\bar{\chi}$ and
$f\bar{f}\rightarrow\chi\bar{\chi}$ processes. At large temperature i.e. when
$(T^{\prime}\gtrsim m_{\chi})$, DM production is dominated by
$Z^{\prime}Z^{\prime}\rightarrow\chi\bar{\chi}$ process and it is shown by the
red solid line. Now, due to the presence of the freeze-in term in the
Boltzmann equation, DM can also be produced from the annihilation of SM
fermions and this production will stop at $T\sim m_{\chi}$. The solid blue
line indicates the production of DM from SM annihilation. Therefore, for
$\alpha_{X}=2.35\times 10^{-5}$, DM is produced via freeze-in mechanism. Now,
if we increase the value of $\alpha_{X}$, one of the possible production
mechanisms leak-in is observed, which is shown by brown dot-dashed line. One
can notice that DM produced through leak-in mechanism is over abundant. Now,
for $\alpha_{X}=0.144$, variation of the relative abundance of DM is shown by
black solid line. As we can see, initially the relative abundance follows its
equilibrium value i.e. the evolution of $\Omega h^{2}_{\rm eq}/\Omega
h^{2}_{\rm obs}$ where $\Omega h^{2}_{\rm eq}=2.755\times 10^{8}(m_{\chi}/1\rm
GeV)Y^{\rm eq}_{\chi_{\rm tot}}$ and it is depicted by black dotted line. When
$T^{\prime}<m_{\chi}$, number density of DM started to fall off exponentially
but due to the presence of freeze-in term $Y_{\chi_{\rm tot}}$ starts to
increase and at this time it follows a quasi static equilibrium. Quasi static
equilibrium is established when DM annihilation rate due to annihilation term
and DM production rate due to freeze-in term are equal Cheung:2010gj . DM
comoving number density in this quasi static equilibrium is denoted by
$Y_{\rm{QSE}}$ and it is independent of the expansion parameter $n$. The
analytic form of $Y_{\rm QSE}$ is given by

$\displaystyle Y_{\rm{QSE}}$ $\displaystyle=$
$\displaystyle\dfrac{\langle\sigma
v\rangle^{T}_{f\bar{f}\to\chi\bar{\chi}}n_{f}^{\rm eq}(T)^{2}}{\langle\sigma
v\rangle^{T^{\prime}}_{\chi\bar{\chi}\to Z^{\prime}Z^{\prime}}s(T)^{2}}\,\,.\
$ (29)

In this figure we show the variation of the quantity $\Omega h^{2}_{\rm
QSE}/\Omega h^{2}_{\rm obs}$ by the dashed black line and $\Omega h^{2}_{\rm
QSE}=2.755\times 10^{8}(m_{\chi}/1{\rm GeV})Y_{\rm QSE}$. Here, one can see
that nearly at $x\simeq 10$, DM departs from quasi static equilibrium and
freezes out. This phenomena is known as reannihilation. In the right panel of
Fig. 4, we have shown that for $m_{\chi}=10^{5}\rm GeV$, only possible way to
produce correct relic abundance of DM is freeze-in and it is shown by solid
blue line. For $\alpha_{X}=1$ (black solid line) and $\alpha_{X}=4\pi$ (red
solid line), excess amount of DM produced from the SM bath cannot annihilate
sufficiently into $Z^{\prime}$ to give correct relic density. Thus for $n=1$,
$m_{\chi}=10^{3}\rm GeV$ and $10^{5}\rm GeV$ reannihilation and freeze-in are
the only possible way to produce observed relic density. Let us note that, to
satisfy the relic density constraint, we need larger values of $\alpha_{X}$ in
comparison to the standard freeze-in scenario Hall:2009bx .

For $n=2$, we show the evolution of relative abundance of DM for
$\epsilon=10^{-7}$ in Fig. 5. In the left panel of Fig. 5, we have considered
$m_{\chi}=10^{3}\rm GeV$. Here, we can see that for $\alpha_{X}=6.8\times
10^{-5}$, DM is produced via freeze-in and it is shown by solid blue line. The
red line over blue line indicates the DM production from $Z^{\prime}$
annihilation in the early Universe. For $\alpha_{X}=0.8$, DS can internally
thermalize and one of the possible production mechanisms is reannihilation,
shown by black solid line. Evolution of $\Omega h^{2}_{\rm QSE}/\Omega
h^{2}_{\rm obs}$ and $\Omega h^{2}_{\rm eq}/\Omega h^{2}_{\rm obs}$ are shown
by the black dashed and black dotted lines respectively. DM annihilation into
$Z^{\prime}$ increases if we increase $\alpha_{X}$ and the magenta line
represents the evolution of relative abundance of DM for $\alpha_{X}=4\pi$.
Thus, possible ways of DM production in this case are freeze-in and
reannihilation. In the right panel of Fig. 5, we have considered
$m_{\chi}=10^{5}\rm GeV$. One can see from this figure that initially the DM
does not follow the equilibrium comoving number density. This is because the
annihilation cross section $\chi\bar{\chi}\to Z^{\prime}Z^{\prime}$ decreases
as we increase $m_{\chi}$ and as a result the DM annihilation cross section is
not sufficient to thermalise the DM. Here, we have shown that for
$\alpha_{X}=6.9\times 10^{-3}$ and $1$, relic abundance of DM can be obtained
via freeze-in mechanism and they are denoted by solid blue and dashed black
lines respectively. For $\alpha_{X}=4\pi$, excess amount of DM is produced
from the SM bath and it is represented by red dot-dashed line. Now, due to
large value of $\alpha_{X}$, excess DM produced from the SM bath, starts to
annihilate into $Z^{\prime}$. We call this phenomenon as “late-time
annihilation”. However, within the perturbative limit of $\alpha_{X}$,
annihilation cross section of $\chi\bar{\chi}\to Z^{\prime}Z^{\prime}$ is not
sufficiently large to give correct relic abundance of DM. Here, one can see
that for $\alpha_{X}=1$ and $4\pi$, DM is overabundant.

Figure 5:  Left panel: Evolution of relative abundance of DM for
$m_{\chi}=10^{3}\rm GeV$. For $\alpha_{X}=6.8\times 10^{-5}$ and $0.8$ DM
production mechanisms are freeze-in (blue solid line) and reannhilation (black
solid line) respectively. Red line over blue line indicates DM production from
$Z^{\prime}$. For $\alpha_{X}=4\pi$, DM is underabundant and it is shown with
solid magenta line. The evolution of $\Omega h^{2}_{\rm eq}/\Omega h^{2}_{\rm
obs}$ and $\Omega h^{2}_{\rm QSE}/\Omega h^{2}_{\rm obs}$ are denoted by black
dotted and black dashed lines respectively Right panel: Evolution of DM
relative abundance for $m_{\chi}=10^{5}\rm GeV$. For $\alpha_{X}=6.9\times
10^{-3},1,\,\text{and }4\pi$, DM production mechanisms are freeze-in (blue
solid line), late-time annihilation (black dashed line), and late-time
annihilation (red dot-dashed line) respectively. Here we choose $n=2$,
$\epsilon=10^{-7}$, and $m_{Z^{\prime}}=100\rm MeV$ in both the figures.

Figure 6: Relative abundance of DM as a function of $x$ for
$m_{\chi}=10^{3}\rm GeV$, $\epsilon=10^{-6}$, and $m_{Z^{\prime}}=100\rm MeV$.
Left panel: We consider $n=3$. For $\alpha_{X}=1.9\times 10^{-3}$ (solid black
line) and $\alpha_{X}=1.58$ (red dashed line) DM is produced via freeze-in and
late-time annihilation mechanisms respectively. Right panel: DM is produced
via freeze-in and late-time annihilation mechanisms for $\alpha_{X}=4.83\times
10^{-2}$ (black solid line) and $1.93$ (red dashed line) respectively. Here we
have considered the expansion parameter $n=4$. In both panels
$x_{r}=m_{\chi}/T_{r}$.

As $n$ increases, parameter space for internally thermalized DS reduces. As a
result, the possible ways of producing DM are freeze-in and late-time
annihilation. As we can see in the left panel of Fig. 6, for
$m_{\chi}=10^{3}\rm{GeV}$, one of the possible mechanisms here is freeze-in
when $\alpha_{X}=1.9\times 10^{-3}$, and it is shown in black solid line. The
red dashed line depicts a new possible way of producing DM for $n>2$. Here we
see that for $\alpha_{X}=1.58$, the dark sector is not internally thermalized
and DM is produced from SM bath via freeze-in mechanism. Due to the large
value of $\alpha_{X}$, the excess amount of DM produced via freeze-in can
annihilate into $Z^{\prime}$ and finally DM freezes out and the final
abundance of DM is set by the late-time annihilation. Furthermore, due to the
presence of the fast expanding component $\phi$ in the Hubble parameter, DM
keeps annihilating and as a result the freeze-out process is prolonged. This
is known as “relentless DM” DEramo:2017gpl . In right panel of Fig. 6 similar
features are shown for $n=4$. Here, we can see a plateau region for large
$\alpha_{X}$, this is because freeze-in mechanism stops when $x\simeq 1$ and
at this time annihilation rate of DM is much less than $H$. However, it is
expected that the annihilation of DM will start again due to the presence of
the fast expanding component in the Hubble parameter. We have numerically
checked that $\Gamma_{\chi\bar{\chi}\to Z^{\prime}Z^{\prime}}\simeq H$ near
$x\simeq 10^{3}$ and DM starts to annihilate into $Z^{\prime}$ and final relic
density of DM is determined by the late time annihilation process.

Figure 7:  Evolution of relative abundance of DM as a function of $x$ for
$\epsilon=7\times 10^{-4}$ and $m_{Z^{\prime}}=100\rm MeV$, and $n=4$ for two
different set of values of $(m_{\chi},\,\alpha_{X})$. The evolution for
$(10^{4}{\rm GeV},\,9.8\times 10^{-6})$ and $(10^{6}{\rm GeV},\,9.8\times
10^{-2})$ are denoted by red dashed line and solid black line respectively.
Figure 8: Constraints in the $m_{\chi}-\epsilon$ plane. In the yellow region,
dark and visible sectors are in thermal equilibrium. The relic density
satisfied region is depicted by the light red region. Cyan region denotes the
absolute coupling floor for the leak-in scenario. DM self interaction
constraint is shown by the blue region. The horizontal solid lines denote the
required value of $\epsilon$ to satisfy the muon $(g-2)$ anomaly within
$2\sigma$ limit for $m_{Z^{\prime}}=100\rm MeV$ (violet line) and
$m_{Z^{\prime}}=50\rm MeV$ (black line). In this figure, we have considered
the value of the expansion parameter $n=4$.

In section 3.2, we mentioned a certain region of parameter space for a
thermally decoupled dark sector in which it is possible to explain muon
$(g-2)$ anomaly. Motivated by this result, in Fig. 7, we show the evolution of
relative DM abundance for different values of $m_{\chi}$ and $\alpha_{X}$ and
the choices of $m_{Z^{\prime}}$ and $\epsilon$ are such that it can enhance
the theoretical prediction of muon $(g-2)$. One can see that in this figure,
the DM is produced via freeze-in mechanism. This is because for large values
of $n$, internal thermalization of dark sector is not possible within the
perturbative limit of $\alpha_{X}$. An important point to note is that in
standard cosmology, DM relic density produced via freeze-in mechanism, is
independent of DM mass and it only depends on coupling strength Hall:2009bx .
However, for fast expansion, the relic density depends on $m_{\chi}$,
$\epsilon$, and $\alpha_{X}$. Therefore, for fixed $\epsilon$ we need larger
value of $\alpha_{X}$ as we increase $m_{\chi}$. The calculations related to
freeze-in in fast expanding Universe is discussed in appendix A.

In Fig. 8 we have shown DM relic density satisfied region and different
constraints in case of $n=4$ in $m_{\chi}-\epsilon$ plane. As we have
discussed earlier, DM production mechanism for $n=4$ is freeze-in and late
time annihilation. Relic satisfied region is shown by the light red colored
region. Here, one can see that the required value of $\epsilon$ increases as
we increase DM mass and this can be easily understood from Eq. 41. In the
yellow shaded region of the figure, dark and visible sectors are in thermal
equilibrium. The absolute coupling floor for leak-in scenario is shown by the
cyan region. The horizontal violet and black solid lines denote the required
values of $\epsilon$ for $m_{Z^{\prime}}=100\rm MeV$ and $50\rm MeV$ to
satisfy the muon $(g-2)$ anomaly respectively. We have also studied the
constraint from self interaction of DM which is discussed later and the
constrained region in $m_{\chi}-\epsilon$ plane is shown by the blue shaded
region of Fig. 8.

##  6 Relevant constraints on model

###  6.1 Self interacting dark matter (SIDM)

DM self interaction cross-section per unit DM mass can be constrained from
bullet cluster observation and the upper limit of $\sigma^{T}/m_{\chi}$ is
$1.25\,\rm cm^{2}g^{-1}$ Randall:2008ppe  where $\sigma^{T}$ is the momentum
transfer cross section of DM elastic scattering processes. In our model, DM
elastic scattering processes are mediated by dark vector boson $Z^{\prime}$.
We have considered $m_{Z^{\prime}}=100\,\rm MeV$ and choice of
$m_{Z^{\prime}}$ is motivated by the explanation of muon $(g-2)$ anomaly. We
calculate total $\sigma^{T}$ and it is defined as

$\displaystyle\sigma^{T}=\dfrac{1}{4}\left[\sigma^{T}_{\chi\bar{\chi}\to\chi\bar{\chi}}+\sigma^{T}_{\chi\chi\to\chi\chi}+\sigma^{T}_{\bar{\chi}\bar{\chi}\to\bar{\chi}\bar{\chi}}\right]\,.$
(30)

The blue region of Fig. 8 depicts the excluded region of parameter space from
DM self interaction. In deriving the SIDM bound, we have expressed
$\alpha_{X}$ in terms of $m_{\chi}$ and $\epsilon$ using Eq. 41 and the relic
density constraint $\Omega h^{2}_{\rm obs}=0.12$.

###  6.2 Muon anomalous magnetic moment

In our model, the dark vector boson $Z^{\prime}$ can contribute to the
anomalous magnetic moment of muon if its coupling strength with muon is $\sim
10^{-3}$ and $m_{Z^{\prime}}\sim 100\rm MeV$. As we discuss in the earlier
sections, it is possible to satisfy the relic density constraint for
$\epsilon\sim 10^{-3}$ if we consider fast expanding Universe. Thus in our
framework, one can have simultaneous explanations of DM relic density and muon
$(g-2)$.

The contribution of $Z^{\prime}$ to $\Delta a_{\mu}$ is given by Baek:2001kca
; Ma:2001md ; Escudero:2019gzq ; Banerjee:2020zvi

$\displaystyle\Delta a_{\mu}$ $\displaystyle=$
$\displaystyle\dfrac{\epsilon^{2}}{4\pi^{2}}\int_{0}^{1}\dfrac{m_{\mu}^{2}z\left(1-z\right)^{2}}{m_{\mu}^{2}\left(1-z\right)^{2}+m_{Z^{\prime}}^{2}z}dz\,\,.\
$ (31)

The present value of $\Delta a_{\mu}=(251\pm 59)\times 10^{-11}$
Muong-2:2021ojo  and the yellow region of Fig. 9 depicts the muon $(g-2)$
satisfied region within $2\sigma$ limit in $m_{Z^{\prime}}-\epsilon$ plane.

Figure 9: Constraints in $m_{Z^{\prime}}-\epsilon$ plane. The cyan region
denotes the constraint from the neutrino trident production at CCFR
CCFR:1991lpl ; Altmannshofer:2014pba . BBN Constraint is shown by the gray
colored region. Brown region depicts the constraint from white dwarf cooling.
The observed deviation from the SM prediction of muon $(g-2)$ can be explained
in the yellow region within $2\sigma$ limit. In the black hatched region,
depending on the expansion parameter $n$ one can satisfy the relic density
constraint. The left side of the black dashed line represents the region where
$Z^{\prime}$ is in thermal equilibrium neutrino bath. The constraints from
Borexino Bellini:2011rx  and BABAR BaBar:2016sci  experiments discussed in
Gninenko:2020xys , are shown by violet and green regions respectively
Gninenko:2020xys .

###  6.3 White dwarf cooling

Several astrophysical observations constrain $L_{\mu}-L_{\tau}$ portal. For
example, SN1987A constrains very low portal coupling Fradette:2014sza ;
Chang:2016ntp ; Escudero:2019gzq  whereas presence of $Z^{\prime}$ inside
stellar core affects the rate of stellar cooling An:2013yfc ; Hardy:2016kme
and the parameter space of very low mass $Z^{\prime}$ is severely constrained
from the stellar cooling constraint. For the mass range of $Z^{\prime}$ in
which we are interested in, most relevant constraint comes from white dwarf
(WD) cooling.

The WD cooling constraint in our scenario can be derived from the following
effective electron-neutrino interaction Lagrangian

$\displaystyle\mathcal{L}^{\rm NP}_{\nu e}\supset-C_{\rm
WD}\left(\bar{\nu}_{l}\gamma_{\mu}P_{L}\nu_{l}\right)(\bar{e}\gamma^{\mu}e)\,\,,$
(32)

where

$\displaystyle C_{\rm
WD}=\frac{\epsilon}{m_{Z^{\prime}}^{2}}\sqrt{4\pi\alpha_{\rm{em}}}\,\Pi(0)\,\,.$
(33)

The $\gamma-Z^{\prime}$ kinetic mixing at one loop level $\Pi(q^{2})$ is given
by Araki:2017wyg ; Banerjee:2018mnw

$\displaystyle
i\Pi(q^{2})=-\frac{ie\epsilon}{2\pi^{2}}\int_{0}^{1}x(1-x)\ln\left(\dfrac{m_{\mu}^{2}-q^{2}x(1-x)}{m_{\tau}^{2}-q^{2}x(1-x)}\right)dx\,\,\,.$
(34)

In the low momentum transfer limit i.e. $q^{2}\to 0$, Eq. 34 takes the
following form

$\displaystyle
i\Pi(0)=-\frac{ie\epsilon}{12\pi^{2}}\ln\left(\dfrac{m_{\mu}^{2}}{m_{\tau}^{2}}\right)\,\,.$
(35)

From the rate of cooling of WD, the Wilson coefficient $C_{\rm WD}$ of
neutrino-electron interaction is constrained as Dreiner:2013tja ;
Arguelles:2022xxa

$\displaystyle\dfrac{1.12\times
10^{-5}}{\rm{GeV}^{2}}<C_{\rm{WD}}<\dfrac{4.5\times
10^{-3}}{\rm{GeV}^{2}}\,\,.$ (36)

We have shown the constraint coming from WD cooling by brown colored region in
Fig. 9.

###  6.4 BBN constraint

Phenomenology of early Universe near $1\rm{MeV}$ is highly constrained from
BBN predictions Mukhanov:2003xs ; Mangano:2011ar . Due to presence of light
$Z^{\prime}$, it is possible that $Z^{\prime}$ can thermalise with the SM
neutrinos via $Z^{\prime}\leftrightarrow\bar{\nu}_{i}\nu_{i}$ ($i=\mu,\,\tau$)
and contribute to the effective number of neutrino species ($N_{\nu}$) at the
time of BBN.

To check the thermalisation of $Z^{\prime}$ at the time of BBN, we have
compared the reaction rate of $\bar{\nu}_{i}\nu_{i}\leftrightarrow Z^{\prime}$
with the Hubble parameter at the time of BBN (i.e. at $T=1\rm MeV$) and found
that the $Z^{\prime}$ is in thermal equilibrium with the neutrino bath for
$10^{-4}\leq\epsilon\leq 10^{-2}$ and ${\rm 1MeV}\leq m_{Z^{\prime}}\lesssim
40\rm MeV$. This region is outlined by the black dashed line in Fig. 9. In
this region of parameter space, depending on the mass of $Z^{\prime}$, it may
give additional contribution to the effective number of neutrino species
($\Delta N_{\nu}\equiv N_{\nu}-3$) whose upper limit is 0.168 at 95% C.L
Fields:2019pfx . We have calculated the constraint on the $Z^{\prime}$ mass
from $\Delta N_{\nu}$ and found that $Z^{\prime}$ of mass greater than $5.5\rm
MeV$ is allowed from the $\Delta N_{\nu}$ constraint even if it is in thermal
equilibrium with the neutrino bath. The grey region of Fig. 9 depicts the
constraint on $m_{Z^{\prime}}$ from $\Delta N_{\nu}$. Let us note in passing
that in our parameter space of interest, non-thermal production $Z^{\prime}$
can happen for $m_{Z^{\prime}}\gtrsim 40\rm MeV$ and it does not contribute to
the $\Delta N_{\nu}$. This is because i) the energy density of the non-
thermally produced $Z^{\prime}$s are smaller compared to the thermally
produced $Z^{\prime}$ and they are non-relativistic at the time of BBN, and
ii) the decay lifetime of $Z^{\prime}$ in our parameter space of interest is
much smaller than $1s$. Thus the produced $Z^{\prime}$ decays much earlier
than SM neutrino decoupling and hence it cannot heat up the neutrino bath.

##  7 Summary and conclusions

In this work we have considered two sectors, dark and visible sectors in non-
standard cosmological background. Here, by non-standard cosmology, we mean the
presence of a new species $\phi$ which redshifts faster than the radiation and
total energy density in the early Universe is dominated by $\phi$. We have
analyzed DM evolution in the fast expanding Universe when dark and visible
sectors are not in thermal equilibrium. To provide a detailed discussion, we
have considered $U(1)_{L_{\mu}-L_{\tau}}\otimes U(1)_{X}$ gauge extension of
SM. A Dirac fermion $\chi$ is the DM candidate which is charged only under
$U(1)_{X}$ gauge symmetry. We have assumed tree level kinetic mixing between
$U(1)_{X}$ and $U(1)_{L_{\mu}-L_{\tau}}$ which generates an interaction at
tree level between $\mu$, $\tau$ flavored SM leptons with the dark gauge boson
$Z^{\prime}$. In standard cosmological scenario i.e. in the radiation
dominated Universe, there exists an upper limit of $\epsilon$ for a fixed
value of $m_{\chi}$ above which both the sectors will be in thermal
equilibrium. In the presence of $\phi$, the upper limit of the portal coupling
is larger compared to the radiation dominated Universe. We have shown that for
$n=1,\,2$, correct relic abundance of DM can be produced either via
reannihilation or by freeze-in mechanisms whereas for $n=3,\,4$, DM production
mechanisms are freeze-in and late time annihilation. In case of late time
annihilation, excess DM is produced from freeze-in mechanism and annihilates
into $Z^{\prime}$.

Furthermore by considering the fast expanding Universe, it is possible to
reconcile DM relic density and muon $(g-2)$ anomaly in our framework. We have
also investigated the constraints from CCFR, Borexino, BABAR, WD cooling and
BBN and found that $m_{Z^{\prime}}-\epsilon$ plane is stringently constrained
from these observations. In particular, we have shown that for $n=4$,
reconciliation of DM relic density for a thermally decoupled dark sector and
muon anomalous magnetic moment is possible while satisfying the other
cosmological and astrophysical constraints if $2\times
10^{-4}\lesssim\epsilon\lesssim 10^{-3}$, $5.5{\rm MeV}\lesssim
m_{Z^{\prime}}\lesssim 200{\rm MeV}$, and $1{\rm TeV}\lesssim m_{\chi}\lesssim
10{\rm TeV}$.

##  8 Acknowledgement

SG acknowledges University Grants Commission (UGC), Government of India for
providing senior research fellowship. AT would like to acknowledge the
financial support provided by the Indian Association for the Cultivation of
Science (IACS), Kolkata.

##  Appendix A Freeze-in production in fast expanding Universe

In this appendix, we derive approximate DM yield for pure freeze-in mechanism
from scattering when Universe redshifts faster than radiation. To derive this,
we need to calculate the collision term for
$f(P_{1})\bar{f}(P_{2})\rightarrow\chi(P_{3})\bar{\chi}(P_{4})$ process, and
it is given by,

$\displaystyle\mathcal{C_{E}}$ $\displaystyle=$
$\displaystyle\int\prod_{i=1}^{4}d\Pi_{i}\left(2\pi\right)^{4}\delta^{4}\left(P_{1}+P_{2}-P_{3}-P_{4}\right)\overline{|\mathcal{M}|^{2}}\,f_{1}^{\rm
eq}(E_{1},T)f_{2}^{\rm eq}(E_{2},T)\,\,.\ $ (37)

Here, $P_{i}$, $\vec{p}_{i}$, and $E_{i}$ are four momentum, magnitude of
three momentum, and energy of the $i^{th}$ species. The Lorentz invariant
phase space measure is
$d\Pi_{i}=\dfrac{g_{i}d^{3}\vec{p}_{i}}{(2\pi)^{3}2E_{i}}$ where $g_{i}$ is
the internal degrees of freedom of the $i^{th}$ species. $f_{i}(E_{i},T)$ is
the energy distribution function of the $i^{th}$ species at temperature $T$
and $f_{i}(E_{i},T)=\exp(-E_{i}/T)$. $\overline{|\mathcal{M}|^{2}}$ is the
matrix amplitude square of $f\bar{f}\to\chi\bar{\chi}$ and it is averaged over
the degrees of freedom of initial and final state particles. Now, one can
write $\mathcal{C_{E}}$ as follows Gondolo:1990dk

$\displaystyle\mathcal{C_{E}}$ $\displaystyle=$
$\displaystyle\dfrac{T}{128\pi^{4}}\int_{s_{min}}^{\infty}ds\,\sqrt{s}\,\,\tilde{\sigma}\left(s\right)K_{1}\left(\dfrac{\sqrt{s}}{T}\right)\,\,.\
$ (38)

Here, $s_{min}={\rm Min}\left[4m_{f}^{2},4m_{\chi}^{2}\right]$, $m_{f}$ is the
mass of fermion, $K_{1}(\sqrt{s}/T)$ is the modified Bessel function of second
kind of order one and $\tilde{\sigma}$ is given by

$\tilde{\sigma}=4g_{f}^{2}\left(s-4m_{f}^{2}\right)\sigma\left(f\bar{f}\to\chi\bar{\chi}\right)$

To estimate the relic density of $\chi$, we calculate the annihilation cross
section in the limit $s\gg m_{\chi}^{2},m_{f}^{2}$ and in this limit
$\sigma\left(\chi\bar{\chi}\rightarrow
f\bar{f}\right)\simeq\dfrac{\alpha_{X}\epsilon^{2}}{16\,s}$ and one can write
$\mathcal{C_{E}}$ as follows

$\displaystyle\mathcal{C_{E}}$ $\displaystyle\simeq$
$\displaystyle\dfrac{\alpha_{X}\epsilon^{2}}{16\pi^{4}}T^{3}m_{\chi}K_{1}\left(\dfrac{2m_{\chi}}{T}\right)\,\,.\
$ (39)

To obtain final relic density of DM, we need to solve the following Boltzmann
equation,

$\displaystyle\dfrac{dY}{dx}$ $\displaystyle=$
$\displaystyle\dfrac{1}{H(x)s(x)x}\,\mathcal{C_{E}}(x)\,\,.\ $ (40)

To get an analytical result, we assume that the production of DM from SM bath
stops at $T_{\rm fi}$ and $T_{\rm fi}\gg T_{r}$. Thus, we can use the
functional form of $H$ as given in Eq. 12 and the final comoving abundance of
DM is given by,

$\displaystyle Y_{\infty}$ $\displaystyle\simeq$
$\displaystyle\dfrac{45}{128\,\pi^{7}}\sqrt{\dfrac{5}{4\pi}}\dfrac{\alpha_{X}\epsilon^{2}}{\sqrt{g_{\rho}\left(x_{r}\right)}g_{*s}}\,\dfrac{3\,M_{\rm
Pl}}{m_{\chi}\,x_{r}^{n/2}}\,\Gamma\left(\dfrac{2+n}{4}\right)\Gamma\left(\dfrac{6+n}{4}\right)\,\,,\
$ (41)

where we have neglected variation of degrees of freedom of entropy $g_{*s}$
with SM bath temperature $T$. Now, the relic density of DM can be calculated
from $\Omega h^{2}=2.755\times 10^{8}(m_{\chi}/1\rm GeV)Y_{\infty}$.

## References

  * (1) WMAP collaboration, _Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Parameter Results_ , _Astrophys. J. Suppl._ 208 (2013) 19 [1212.5226].
  * (2) Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * (3) Y. Sofue and V. Rubin, _Rotation curves of spiral galaxies_ , _Ann. Rev. Astron. Astrophys._ 39 (2001) 137 [astro-ph/0010594].
  * (4) D. Clowe, A. Gonzalez and M. Markevitch, _Weak lensing mass reconstruction of the interacting cluster 1E0657-558: Direct evidence for the existence of dark matter_ , _Astrophys. J._ 604 (2004) 596 [astro-ph/0312273].
  * (5) D. Clowe, M. Bradac, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones et al., _A direct empirical proof of the existence of dark matter_ , _Astrophys. J. Lett._ 648 (2006) L109 [astro-ph/0608407].
  * (6) M. Srednicki, R. Watkins and K. A. Olive, _Calculations of Relic Densities in the Early Universe_ , _Nucl. Phys. B_ 310 (1988) 693.
  * (7) P. Gondolo and G. Gelmini, _Cosmic abundances of stable particles: Improved analysis_ , _Nucl. Phys. B_ 360 (1991) 145.
  * (8) G. Jungman, M. Kamionkowski and K. Griest, _Supersymmetric dark matter_ , _Phys. Rept._ 267 (1996) 195 [hep-ph/9506380].
  * (9) J. L. Feng, _Dark Matter Candidates from Particle Physics and Methods of Detection_ , _Ann. Rev. Astron. Astrophys._ 48 (2010) 495 [1003.0904].
  * (10) L. Roszkowski, E. M. Sessolo and S. Trojanowski, _WIMP dark matter candidates and searches—current status and future prospects_ , _Rept. Prog. Phys._ 81 (2018) 066201 [1707.06277].
  * (11) T. Lin, _Dark matter models and direct detection_ , _PoS_ 333 (2019) 009 [1904.07915].
  * (12) Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * (13) A. Czarnecki, W. J. Marciano and A. Vainshtein, _Refinements in electroweak contributions to the muon anomalous magnetic moment_ , _Phys. Rev. D_ 67 (2003) 073006 [hep-ph/0212229].
  * (14) K. Melnikov and A. Vainshtein, _Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment revisited_ , _Phys. Rev. D_ 70 (2004) 113006 [hep-ph/0312226].
  * (15) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, _Complete Tenth-Order QED Contribution to the Muon g-2_ , _Phys. Rev. Lett._ 109 (2012) 111808 [1205.5370].
  * (16) C. Gnendiger, D. Stöckinger and H. Stöckinger-Kim, _The electroweak contributions to $(g-2)_{\mu}$ after the Higgs boson mass measurement_, _Phys. Rev. D_ 88 (2013) 053005 [1306.5546].
  * (17) A. Kurz, T. Liu, P. Marquard and M. Steinhauser, _Hadronic contribution to the muon anomalous magnetic moment to next-to-next-to-leading order_ , _Phys. Lett. B_ 734 (2014) 144 [1403.6400].
  * (18) G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, _Remarks on higher-order hadronic corrections to the muon g $-$2_, _Phys. Lett. B_ 735 (2014) 90 [1403.7512].
  * (19) P. Masjuan and P. Sanchez-Puertas, _Pseudoscalar-pole contribution to the $(g_{\mu}-2)$: a rational approach_, _Phys. Rev. D_ 95 (2017) 054026 [1701.05829].
  * (20) G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, _Dispersion relation for hadronic light-by-light scattering: two-pion contributions_ , _JHEP_ 04 (2017) 161 [1702.07347].
  * (21) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, _Reevaluation of the hadronic vacuum polarisation contributions to the Standard Model predictions of the muon $g-2$ and ${\alpha(m_{Z}^{2})}$ using newest hadronic cross-section data_, _Eur. Phys. J. C_ 77 (2017) 827 [1706.09436].
  * (22) A. Keshavarzi, D. Nomura and T. Teubner, _Muon $g-2$ and $\alpha(M_{Z}^{2})$: a new data-based analysis_, _Phys. Rev. D_ 97 (2018) 114025 [1802.02995].
  * (23) M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold and S. P. Schneider, _Dispersion relation for hadronic light-by-light scattering: pion pole_ , _JHEP_ 10 (2018) 141 [1808.04823].
  * (24) G. Colangelo, M. Hoferichter and P. Stoffer, _Two-pion contribution to hadronic vacuum polarization_ , _JHEP_ 02 (2019) 006 [1810.00007].
  * (25) T. Aoyama, T. Kinoshita and M. Nio, _Theory of the Anomalous Magnetic Moment of the Electron_ , _Atoms_ 7 (2019) 28.
  * (26) A. Gérardin, H. B. Meyer and A. Nyffeler, _Lattice calculation of the pion transition form factor with $N_{f}=2+1$ Wilson quarks_, _Phys. Rev. D_ 100 (2019) 034520 [1903.09471].
  * (27) M. Hoferichter, B.-L. Hoid and B. Kubis, _Three-pion contribution to hadronic vacuum polarization_ , _JHEP_ 08 (2019) 137 [1907.01556].
  * (28) J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, _Short-distance constraints for the HLbL contribution to the muon anomalous magnetic moment_ , _Phys. Lett. B_ 798 (2019) 134994 [1908.03331].
  * (29) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, _A new evaluation of the hadronic vacuum polarisation contributions to the muon anomalous magnetic moment and to $\mathbf{\bm{\alpha}(m_{Z}^{2})}$_, _Eur. Phys. J. C_ 80 (2020) 241 [1908.00921].
  * (30) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, _Longitudinal short-distance constraints for the hadronic light-by-light contribution to $(g-2)_{\mu}$ with large-$N_{c}$ Regge models_, _JHEP_ 03 (2020) 101 [1910.13432].
  * (31) A. Keshavarzi, D. Nomura and T. Teubner, _$g-2$ of charged leptons, $\alpha(M^{2}_{Z})$ , and the hyperfine splitting of muonium_, _Phys. Rev. D_ 101 (2020) 014029 [1911.00367].
  * (32) T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung et al., _Hadronic Light-by-Light Scattering Contribution to the Muon Anomalous Magnetic Moment from Lattice QCD_ , _Phys. Rev. Lett._ 124 (2020) 132002 [1911.08123].
  * (33) T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , _Phys. Rept._ 887 (2020) 1 [2006.04822].
  * (34) Muon g-2 collaboration, _Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL_ , _Phys. Rev. D_ 73 (2006) 072003 [hep-ex/0602035].
  * (35) X. G. He, G. C. Joshi, H. Lew and R. R. Volkas, _NEW Z-prime PHENOMENOLOGY_ , _Phys. Rev. D_ 43 (1991) 22.
  * (36) X.-G. He, G. C. Joshi, H. Lew and R. R. Volkas, _Simplest Z-prime model_ , _Phys. Rev. D_ 44 (1991) 2118.
  * (37) S. Baek, N. G. Deshpande, X. G. He and P. Ko, _Muon anomalous g-2 and gauged L(muon) - L(tau) models_ , _Phys. Rev. D_ 64 (2001) 055006 [hep-ph/0104141].
  * (38) K. Harigaya, T. Igari, M. M. Nojiri, M. Takeuchi and K. Tobe, _Muon g-2 and LHC phenomenology in the $L_{\mu}-L_{\tau}$ gauge symmetric model_, _JHEP_ 03 (2014) 105 [1311.0870].
  * (39) W. Altmannshofer, C.-Y. Chen, P. S. Bhupal Dev and A. Soni, _Lepton flavor violating Z’ explanation of the muon anomalous magnetic moment_ , _Phys. Lett. B_ 762 (2016) 389 [1607.06832].
  * (40) E. Ma, D. P. Roy and S. Roy, _Gauged L(mu) - L(tau) with large muon anomalous magnetic moment and the bimaximal mixing of neutrinos_ , _Phys. Lett. B_ 525 (2002) 101 [hep-ph/0110146].
  * (41) H. Banerjee, P. Byakti and S. Roy, _Supersymmetric gauged U(1) ${}_{L_{\mu}-L_{\tau}}$ model for neutrinos and the muon $(g-2)$ anomaly_, _Phys. Rev. D_ 98 (2018) 075022 [1805.04415].
  * (42) H. Banerjee, B. Dutta and S. Roy, _Supersymmetric gauged $U(1)_{L_{\mu}-L_{\tau}}$ model for electron and muon $g-2$ anomaly_, _JHEP_ 03 (2021) 211 [2011.05083].
  * (43) S. Baek, H. Okada and K. Yagyu, _Flavour Dependent Gauged Radiative Neutrino Mass Model_ , _JHEP_ 04 (2015) 049 [1501.01530].
  * (44) C.-H. Chen and T. Nomura, _Neutrino mass in a gauged $L_{\mu}-L_{\tau}$ model_, _Nucl. Phys. B_ 940 (2019) 292 [1705.10620].
  * (45) T. Hambye, M. H. Tytgat, J. Vandecasteele and L. Vanderheyden, _Dark matter from dark photons: A taxonomy of dark matter production_ , _Physical Review D_ 100 (2019) .
  * (46) J. A. Evans, C. Gaidau and J. Shelton, _Leak-in Dark Matter_ , _JHEP_ 01 (2020) 032 [1909.04671].
  * (47) A. Tapadar, S. Ganguly and S. Roy, _Non-adiabatic evolution of dark sector in the presence of $U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry_, _JCAP_ 05 (2022) 019 [2109.13609].
  * (48) F. D’Eramo, N. Fernandez and S. Profumo, _When the Universe Expands Too Fast: Relentless Dark Matter_ , _JCAP_ 05 (2017) 012 [1703.04793].
  * (49) F. D’Eramo, N. Fernandez and S. Profumo, _Dark Matter Freeze-in Production in Fast-Expanding Universes_ , _JCAP_ 02 (2018) 046 [1712.07453].
  * (50) G. Bellini et al., _Precision measurement of the 7Be solar neutrino interaction rate in Borexino_ , _Phys. Rev. Lett._ 107 (2011) 141302 [1104.1816].
  * (51) BaBar collaboration, _Search for a muonic dark force at BABAR_ , _Phys. Rev. D_ 94 (2016) 011102 [1606.03501].
  * (52) CCFR collaboration, _Neutrino tridents and W Z interference_ , _Phys. Rev. Lett._ 66 (1991) 3117.
  * (53) E. C. G. Stueckelberg, _Interaction energy in electrodynamics and in the field theory of nuclear forces_ , _Helv. Phys. Acta_ 11 (1938) 225.
  * (54) H. Ruegg and M. Ruiz-Altaba, _The Stueckelberg field_ , _Int. J. Mod. Phys. A_ 19 (2004) 3265 [hep-th/0304245].
  * (55) M. Pospelov, A. Ritz and M. B. Voloshin, _Secluded WIMP Dark Matter_ , _Phys. Lett. B_ 662 (2008) 53 [0711.4866].
  * (56) J. L. Feng, H. Tu and H.-B. Yu, _Thermal Relics in Hidden Sectors_ , _JCAP_ 10 (2008) 043 [0808.2318].
  * (57) C. Cheung, G. Elor, L. J. Hall and P. Kumar, _Origins of Hidden Sector Dark Matter I: Cosmology_ , _JHEP_ 03 (2011) 042 [1010.0022].
  * (58) X. Chu, T. Hambye and M. H. G. Tytgat, _The Four Basic Ways of Creating Dark Matter Through a Portal_ , _JCAP_ 05 (2012) 034 [1112.0493].
  * (59) J. A. Evans, S. Gori and J. Shelton, _Looking for the WIMP Next Door_ , _JHEP_ 02 (2018) 100 [1712.03974].
  * (60) R. J. Scherrer and M. S. Turner, _On the Relic, Cosmic Abundance of Stable Weakly Interacting Massive Particles_ , _Phys. Rev. D_ 33 (1986) 1585.
  * (61) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080 [0911.1120].
  * (62) S. W. Randall, M. Markevitch, D. Clowe, A. H. Gonzalez and M. Bradac, _Constraints on the Self-Interaction Cross-Section of Dark Matter from Numerical Simulations of the Merging Galaxy Cluster 1E 0657-56_ , _Astrophys. J._ 679 (2008) 1173 [0704.0261].
  * (63) M. Escudero, D. Hooper, G. Krnjaic and M. Pierre, _Cosmology with A Very Light L μ $-$ Lτ Gauge Boson_, _JHEP_ 03 (2019) 071 [1901.02010].
  * (64) W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, _Neutrino Trident Production: A Powerful Probe of New Physics with Neutrino Beams_ , _Phys. Rev. Lett._ 113 (2014) 091801 [1406.2332].
  * (65) S. Gninenko and D. Gorbunov, _Refining constraints from Borexino measurements on a light Z’-boson coupled to L $\mu$-L$\tau$ current_, _Phys. Lett. B_ 823 (2021) 136739 [2007.16098].
  * (66) A. Fradette, M. Pospelov, J. Pradler and A. Ritz, _Cosmological Constraints on Very Dark Photons_ , _Phys. Rev. D_ 90 (2014) 035022 [1407.0993].
  * (67) J. H. Chang, R. Essig and S. D. McDermott, _Revisiting Supernova 1987A Constraints on Dark Photons_ , _JHEP_ 01 (2017) 107 [1611.03864].
  * (68) H. An, M. Pospelov and J. Pradler, _New stellar constraints on dark photons_ , _Phys. Lett. B_ 725 (2013) 190 [1302.3884].
  * (69) E. Hardy and R. Lasenby, _Stellar cooling bounds on new light particles: plasma mixing effects_ , _JHEP_ 02 (2017) 033 [1611.05852].
  * (70) T. Araki, S. Hoshino, T. Ota, J. Sato and T. Shimomura, _Detecting the $L_{\mu}-L_{\tau}$ gauge boson at Belle II_, _Phys. Rev. D_ 95 (2017) 055006 [1702.01497].
  * (71) H. Banerjee and S. Roy, _Signatures of supersymmetry and $L_{\mu}-L_{\tau}$ gauge bosons at Belle-II_, _Phys. Rev. D_ 99 (2019) 035035 [1811.00407].
  * (72) H. K. Dreiner, J.-F. Fortin, J. Isern and L. Ubaldi, _White Dwarfs constrain Dark Forces_ , _Phys. Rev. D_ 88 (2013) 043517 [1303.7232].
  * (73) C. A. Argüelles et al., _Snowmass White Paper: Beyond the Standard Model effects on Neutrino Flavor_ , in _2022 Snowmass Summer Study_ , 3, 2022, 2203.10811.
  * (74) V. F. Mukhanov, _Nucleosynthesis without a computer_ , _Int. J. Theor. Phys._ 43 (2004) 669 [astro-ph/0303073].
  * (75) G. Mangano and P. D. Serpico, _A robust upper limit on $N_{\rm eff}$ from BBN, circa 2011_, _Phys. Lett. B_ 701 (2011) 296 [1103.1261].
  * (76) B. D. Fields, K. A. Olive, T.-H. Yeh and C. Young, _Big-Bang Nucleosynthesis after Planck_ , _JCAP_ 03 (2020) 010 [1912.01132].