#  MS-TP-21-09 A multi-temperature universe can allow a sub-MeV dark photon
dark matter

Amin<EMAIL_ADDRESS>, Wan-Zhe<EMAIL_ADDRESS>, Pran<EMAIL_ADDRESS>and Zhu-Yao
<EMAIL_ADDRESS>

aInstitut für Theoretische Physik, Westfälische Wilhelms-Universität Münster,
Wilhelm-Klemm-Straße 9, 48149 Münster, Germany
bCenter for Joint Quantum Studies and Department of Physics,
School of Science, Tianjin University, Tianjin 300350, PR. China
cDepartment of Physics, Northeastern University, Boston, MA 02115-5000, USA

###### Abstract

An analysis of sub-MeV dark photon as dark matter is given which is achieved
with two hidden sectors, one of which interacts directly with the visible
sector while the second has only indirect coupling with the visible sector.
The formalism for the evolution of three bath temperatures for the visible
sector and the two hidden sectors is developed and utilized in solution of
Boltzmann equations coupling the three sectors. We present exclusion plots
where the sub-MeV dark photon can be dark matter. The analysis can be extended
to a multi-temperature universe with multiple hidden sectors and multiple heat
baths.

###### Contents

  1. 1 Introduction
  2. 2 A model for a multi-temperature universe
  3. 3 Boltzmann equations for yields with three bath temperatures
  4. 4 Temperature evolution in the dark sectors versus in the visible sector
  5. 5 Thermalization and dark freeze-out
  6. 6 Conclusions
  7. A Canonical normalization of extended $G_{\text{S}M}\times U(1)_{X_{1}}\times U(1)_{X_{2}}$ Lagrangian with kinetic and Stueckelberg mass mixings
  8. B Dark photon $\gamma^{\prime}$ and dark $Z^{\prime}$ couplings with Standard Model particles
  9. C Deduction of three temperature Boltzmann equations
  10. D Dark photon and dark fermion scattering cross sections and $Z^{\prime}$ decay width

##  1 Introduction

Supergravity and strings models typically contain hidden sectors with gauge
groups including $U(1)$ gauge group factors. These hidden sectors with $U(1)$
gauge groups can interact feebly with the visible sector and interact feebly
or with normal strength with each other. The fields in the visible and hidden
sectors in general will reside in different heat baths and the universe in
this case will be a multi-temperature universe. The multi-temperature nature
of the universe becomes a relevant issue if the observables in the visible
sector are functions of the visible and the hidden sector heat baths. Such is
the situation if dark matter (DM) resides in the hidden sector but interacts
feebly with the visible sector. In this case an accurate computation of the
relic density requires thermal averaging of cross sections and decay widths
which depend on temperatures of both the visible and the hidden sector heat
baths. In this work we develop a theoretical formalism which can correlate the
evolution of temperatures of the hidden sector and of the visible sectors (for
the specific case of two hidden sectors) in an accurate way.

The formalism noted above is used in the investigation of a dark photon and
dark fermions of hidden sectors as possible DM candidates. There exists a
considerable literature in the study of dark photons [1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16, 17, 18] (for review see [19, 20, 21]) to which
the interested reader is directed. While axions and dark photons in the light
to ultralight mass region (from keV to $10^{-22}$ eV) have been investigated
[11, 22, 23, 13, 24, 25, 26], the sub-MeV dark photon mass range appears
difficult to realize. The problem arises in part because with the visible
sector interacting with a hidden sector via kinetic mixing, the twin
constraints that the dark photon has a lifetime larger than the age of the
universe, and also produce a sufficient amount of DM to populate the universe
are difficult to satisfy. In addition to the relic density constraint there is
also the constraint on dark photon lifetime which needs to be larger than the
age of the universe as well as a constraint from BBN on the light degrees of
freedom. For the case of one hidden sector, these constraints are difficult to
satisfy. Specifically, $\Delta N_{\rm eff}$ at BBN time depends on the ratio
$(T_{\text{hid}}/T_{\gamma})^{4}$. The BBN temperature is typically $\sim 0.1$
MeV and, as will be seen later, for the case of one hidden sector the ratio
$(T_{\text{hid}}/T_{\gamma})^{4}\sim 1$ which gives a contribution to $\Delta
N_{\rm eff}$ at BBN time in excess of the current experimental constraint
which from the combined data from BBN, BAO and CMB [27] is $\Delta N_{\rm
eff}<0.214$. On the other hand, for the case of two hidden sectors it is
possible to satisfy all the current experimental constraints and for that
reason we will focus on the two hidden sector model which is the minimal
extension of one hidden sector model as discussed below.

In this study we show that a sub-MeV dark photon as DM can indeed be realized
in a simple extension of the Standard Model (SM) where the hidden sector is
constituted of two sectors $X_{1}$ and $X_{2}$ where the sector $X_{1}$ has
kinetic mixing with the visible sector while the sector $X_{2}$ has kinetic
and mass mixings only with sector $X_{1}$, c.f., Fig. 1. We assume that the
hidden sector $X_{1}$ has a dark fermion $D$ and its gauge boson $Z^{\prime}$
decays before the BBN while the hidden sector $X_{2}$ has only a dark photon
$\gamma^{\prime}$ which has a lifetime larger than the age of the Universe.
This set up is theoretically more complex because here one has three heat
baths and the computation of the relic density thus depends on three
temperatures, i.e., the temperature $T$ of the visible sector, the temperature
$T_{1}$ of the hidden sector $X_{1}$ and the temperature $T_{2}$ of the hidden
sector $X_{2}$. In the following we develop a formalism that allows one to
compute temperatures of all three heat baths in terms of one common
temperature, which can be chosen to be $T$, $T_{1}$ or $T_{2}$. In the
analysis below it is found convenient to choose the reference temperature to
be $T_{1}$. The outline of the rest of the paper is as follows: in section 2
we discuss the particle physics model used for our multi-temperature universe
and in sections 3 and 4 we write the coupled Boltzmann equations with three
temperatures and derive the temperature evolution of $T$ and $T_{2}$ relative
to $T_{1}$. Thermalization between the different sectors and dark freeze-out
are explained in section 5 with a numerical analysis and a discussion of the
astrophysical constraints. Conclusions are given in section 6. Further
analytical results are given in Appendices A through D.

Figure 1: The exhibition of the model we consider in the paper. The standard
model has a direct coupling with hidden sector $X_{1}$ with strength
proportional to $\delta_{1}$, whereas hidden sector $X_{2}$ only interacts
directly with $X_{1}$ with strength proportional to $\delta_{2}$, and thus
$X_{2}$ interacts with the standard model only indirectly.

##  2 A model for a multi-temperature universe

As mentioned in the introduction, supergravity and string models contain
hidden sectors. These hidden sectors in general would have both abelian and
non-abelian gauge groups and some of them could interact feebly with the
visible sector while others may interact with each other as shown in Fig. 2.
Thus, for example, in D-brane models one gets $U(N)$ gauge groups where
$U(N)\to SU(N)\times U(1)$. These extra $U(1)$ factors in general can acquire
kinetic mixing with the $U(1)_{Y}$ of the visible sector. Further, the gauge
bosons of the extra $U(1)$’s can acquire mass via the Stueckelberg mechanism
and also have Stueckelberg mass mixings with the hypercharge gauge boson.
Additionally if there are several hidden sector $U(1)$’s they can have also
gauge and Stueckelberg mass mixings among themselves. Interestingly, the
possible existence of these hidden sectors can have significant effect on
model building in the visible sector. As an example one phenomenon which is
deeply affected by the existence of hidden sectors is dark matter which we
discuss in further detail below.

In general the dark sectors with a gauge symmetry will contain gauge fields as
well as matter, but, as noted, typically they will have feeble interactions
with the SM particles and likely also with the inflaton. This means that these
particles would not be thermally produced in the reheat period after inflation
but would then acquire their relic density via annihilation and decay of the
SM particles. Thus in general the temperatures of the visible and the hidden
sectors will be different from the visible sector as well as from each other.
This means that their relic densities will be governed by a set of coupled
Boltzmann equations which depend on different temperatures, i.e., temperature
of the visible sector and those for the hidden sectors. One of the central
items in understanding of how to deal with such coupled systems with sectors
involving different temperatures is to understand fully how the temperatures
of the hidden sectors grow relative to the visible sector temperature. The
formalism of how to correlate the hidden and the visible sector temperatures
was worked out for the case of the visible sector interacting with one hidden
sector in [5]. However, a general framework does not exist. Here we discuss
the case where there are two hidden sectors $X_{1}$ and $X_{2}$ where the
hidden sector $X_{1}$ interacts with the visible sector, while the hidden
sector $X_{2}$ interacts only with the hidden sector $X_{1}$ as shown in Fig.
1. In this case two functions $\eta^{-1}=\xi=T_{1}/T$ and $\zeta=T_{2}/T_{1}$
enter in the coupled Boltzmann equations and we derive differential equations
for their evolution. The above setup has a direct application in achieving a
sub-MeV dark photon as dark matter as we show later in this work. However, a
consistent analysis of the coupled dynamics of the visible sector and two
hidden sectors is significantly more complex. This work develops the necessary
machinery to do so and which can be extended to multiple hidden sectors.

Figure 2: A schematic diagram exhibiting the coupling of the visible sector
with multiple dark sectors and of the dark sectors among themselves. The
visible sector may have direct couplings with some of the dark sectors, or
indirect couplings with others via interactions among the entire hidden
sector.

We assume that the two sectors $X_{1}$ and $X_{2}$ have $U(1)_{X_{1}}$ and
$U(1)_{X_{2}}$ gauge symmetries and that the field content of $X_{1}$ is
$(C_{\mu},D)$ where $C_{\mu}$ is the gauge field and $D$ is a dark fermion,
and the field content of $X_{2}$ is the gauge field $D_{\mu}$ and there is no
dark fermion in the sector $X_{2}$. We invoke a kinetic mixing [28, 29]
between the hypercharge field $B_{\mu}$ of SM and $C_{\mu}$ and a kinetic
mixing between $C_{\mu}$ and $D_{\mu}$ as well as a Stueckelberg mass growth
[30, 31, 32, 33, 34] for all the gauge fields as well as a Stueckelberg mass
mixing between the fields $C_{\mu}$ and $D_{\mu}$. The extended part of the
Lagrangian including both the kinetic and mass mixings is

$\mathcal{L}=\mathcal{L}_{\rm SM}+\mathcal{L}_{\rm kin}+\mathcal{L}_{\rm
mass}+g_{X}J_{X_{1}}^{\mu}C_{\mu},$ (2.1)

where $\mathcal{L}_{\rm SM}$ contains the SM terms and the kinetic part is
given by

$\mathcal{L}_{\text{kin}}=-\frac{1}{4}C^{\mu\nu}C_{\mu\nu}-\frac{1}{4}D^{\mu\nu}D_{\mu\nu}-\frac{\delta_{1}}{2}B^{\mu\nu}C_{\mu\nu}-\frac{\delta_{2}}{2}C^{\mu\nu}D_{\mu\nu},$
(2.2)

and in the unitary gauge the mass Lagrangian is given by

$\mathcal{L}_{\text{mass}}=-\frac{1}{2}(M_{1}C_{\mu}+M_{2}B_{\mu})^{2}-\frac{1}{2}(M_{3}C_{\mu}+M_{4}D_{\mu})^{2}-m_{D}\bar{D}D.$
(2.3)

The $D$ fermion is assumed charged under $U(1)_{X_{1}}$ with interaction
$g_{X}\bar{D}\gamma^{\mu}DC_{\mu}$. Canonical normalization of Eqs. (2.2) and
(2.3) is carried out in Appendix A which gives the mass eigenstates
$\gamma^{\prime},Z^{\prime},Z,\gamma$. The neutral current Lagrangian
contained in $\mathcal{L}_{\rm SM}$ for the mass eigenstates
$\gamma^{\prime},Z^{\prime},Z,\gamma$ describing the couplings between the
vector bosons $\gamma^{\prime},Z^{\prime},Z,\gamma$ with the visible sector
fermions is given by

$\displaystyle\mathcal{L}_{\rm NC}^{\rm v}$
$\displaystyle=\frac{g_{2}}{2\cos\theta_{w}}\bar{f}\gamma^{\mu}[(v_{f}-a_{f}\gamma_{5})Z_{\mu}+(v^{\prime}_{f}-a^{\prime}_{f}\gamma_{5})Z^{\prime}_{\mu}+(v^{\prime\prime}_{f}-a^{\prime\prime}_{f}\gamma_{5})A^{\gamma^{\prime}}_{\mu}]f$
$\displaystyle~{}~{}~{}+e\bar{f}\gamma^{\mu}Q_{f}A^{\gamma}_{\mu}f.$ (2.4)

Here $\theta_{w}$ is the weak angle and $e$ is defined as

$\frac{1}{e^{2}}=\frac{1}{g_{2}^{2}}+\frac{1+\epsilon_{1}^{2}-2\epsilon_{1}\delta_{1}+\epsilon_{1}^{2}\epsilon_{2}^{2}-2\epsilon_{1}^{2}\epsilon_{2}\delta_{2}}{g^{2}_{Y}},$
(2.5)

where $\epsilon_{1}=M_{2}/M_{1}$ and $\epsilon_{2}=M_{3}/M_{4}$. As seen from
Eqs. (2.1)$-$(2.3), the framework of the model allows for the inclusion of
both kinetic and mass mixing between the hidden and visible sectors. However,
in the analysis presented in this work we will only take kinetic mixing
between the visible and the hidden sectors, so that $M_{2}=0$. This is done
because in this case the $4\times 4$ neutral vector boson mass square matrix
factors into a block diagonal form consisting of two $2\times 2$ matrices as
shown in Eq. (A.4). In this case we can carry out a set of eight
$GL(2,\mathbb{R})$ transformations to put the kinetic energy of the visible
and the hidden sectors in a canonical form and at the same time to lowest
order the mass matrix is also in a canonical form. This allows us to use
perturbation theory around stable minima of the standard model mass matrix and
therefore to deduce the couplings between the hidden sectors and the SM
particles as given in Tables 2 and 3. As seen in Appendix A, the analysis is
rather non-trivial and significantly more involved than for the case of one
hidden sector. With the inclusion of $M_{2}$, the analysis becomes
analytically intractable and an exhibition of results corresponding to those
of Tables 2 and 3 is difficult. However, we note that even with $M_{2}=0$, we
still have both kinetic and mass mixing in the hidden sector. Thus
$\delta_{2}$ takes account of kinetic mixing and $M_{3}$ takes account of mass
mixing between the hidden sectors 1 and 2. So in summary in this model the
visible sector has only kinetic mixing with the hidden sector 1 while the
hidden sector 1 has both kinetic and mass mixing with hidden sector 2. It is
seen that the mass mixing from $M_{3}$ does have significant effect on the
model predictions, e.g., on the mass of the dark photon and on the relic
density as seen in model point (d) in Table 1.

With $M_{2}=0$, the neutral current Lagrangian for coupling to the hidden
sector fermion is given by

$\displaystyle\mathcal{L}^{\rm h}_{\rm NC}$
$\displaystyle=(c_{\gamma^{\prime}}A^{\gamma^{\prime}}_{\mu}+c_{Z}Z_{\mu}+c_{Z^{\prime}}Z^{\prime}_{\mu}+c_{\gamma}A^{\gamma}_{\mu})\bar{D}\gamma^{\mu}D,$
(2.6)

where

$\displaystyle c_{\gamma^{\prime}}\simeq$ $\displaystyle
g_{X}\dfrac{m^{2}_{\gamma^{\prime}}}{m^{2}_{Z^{\prime}}-m^{2}_{\gamma^{\prime}}}\delta_{2},~{}c_{Z}\simeq
g_{X}\delta_{1}\sin\theta_{w}(1+\epsilon_{z}^{2}),$ (2.7) $\displaystyle
c_{Z^{\prime}}\simeq$ $\displaystyle g_{X},~{}c_{\gamma}\simeq-
g_{X}\delta_{1}\delta_{2}\left(\dfrac{m_{\gamma^{\prime}}}{m_{Z^{\prime}}}\right)^{2}\sin\beta\cos\theta_{w},$
(2.8) $\displaystyle\tan 2\beta$
$\displaystyle=\frac{2M_{3}M_{4}}{M^{2}_{4}-M^{2}_{1}-M^{2}_{3}},~{}\epsilon_{z}=m_{Z^{\prime}}/m_{Z}.$
(2.9)

The vector and axial-vector couplings with SM fermions appearing in Eq. (2.4)
are given by Eqs. (B.1)$-$(B.6). Those couplings along with the ones in Eq.
(2.7) and Eq. (2.8) are calculated after a proper diagonalization and
normalization of the kinetic and mass square matrices. The complete analysis
is given in Appendix A.

##  3 Boltzmann equations for yields with three bath temperatures

The relic densities of the dark photon and of the dark fermion arise in part
from a freeze-in mechanism [35, 36, 37, 38, 39]. In general the visible sector
and the dark sectors will have different temperatures [40, 41, 42, 43, 44, 45]
(a similar setup has been considered in Ref. [45] but with a different
particle content, couplings and no explicit multi-temperature evolution. Also
the DM candidate was not the dark photon as in our case). As mentioned above,
we consider three different temperatures corresponding to the temperatures of
the visible sector $T$ and of the two hidden sectors, $T_{1}$ for $X_{1}$ and
$T_{2}$ for $X_{2}$. Defining the yield $Y=n/s$, where $n$ is the number
density and $s$ is the entropy density, and the bath functions $\eta$ and
$\zeta$ so that $T=\eta T_{1}$ and $T_{2}=\zeta T_{1}$, we write the Boltzmann
equations for the yields as

$\displaystyle\frac{dY_{D}}{dT_{1}}=$
$\displaystyle-\frac{d\rho/dT_{1}}{4H\rho}s\mathcal{J}_{D},$ (3.1)
$\displaystyle\frac{dY_{Z^{\prime}}}{dT_{1}}=$
$\displaystyle-\frac{d\rho/dT_{1}}{4H\rho}s\mathcal{J}_{Z^{\prime}},$ (3.2)
$\displaystyle\frac{dY_{\gamma^{\prime}}}{dT_{1}}=$
$\displaystyle-\frac{d\rho/dT_{1}}{4H\rho}s\mathcal{J}_{\gamma^{\prime}},$
(3.3)

where $H$ is the Hubble parameter given by

$H^{2}=\frac{8\pi G_{N}}{3}(\rho_{v}+\rho_{1}+\rho_{2}),$ (3.4)

with $\rho=\rho_{v}+\rho_{1}+\rho_{2}$ and $s$ being the energy and entropy
densities given by

$\displaystyle\rho$ $\displaystyle=\frac{\pi^{2}}{30}\left(g_{\rm
eff}^{v}T^{3}+g_{1\rm eff}T_{1}^{3}+g_{2\rm eff}T_{2}^{3}\right),$ (3.5)
$\displaystyle s$ $\displaystyle=\frac{2\pi^{2}}{45}\left(h_{\rm
eff}^{v}T^{3}+h_{1\rm eff}T_{1}^{3}+h_{2\rm eff}T_{2}^{3}\right),$ (3.6)

and the quantities $\mathcal{J}_{D}$, $\mathcal{J}_{Z^{\prime}}$ and
$\mathcal{J}_{\gamma^{\prime}}$ are defined in terms of the collision terms as

$\displaystyle\mathcal{J}_{D}=\frac{C_{D}}{s^{2}},~{}~{}\mathcal{J}_{Z^{\prime}}=\frac{C_{Z^{\prime}}}{s^{2}},~{}~{}\mathcal{J}_{\gamma^{\prime}}=\frac{C_{\gamma^{\prime}}}{s^{2}}.$
(3.7)

The collision terms $C_{D}$, $C_{Z^{\prime}}$ and $C_{\gamma^{\prime}}$ are
given by Eqs. (C.7)$-$(C.9) which allow us to write
$\mathcal{J}_{D},\mathcal{J}_{Z^{\prime}}$ and $\mathcal{J}_{\gamma^{\prime}}$
in terms of the yield as

$\displaystyle\mathcal{J}_{D}$ $\displaystyle=\langle\sigma
v\rangle_{i\bar{i}\to D\bar{D}}(\eta T_{1})Y_{i}^{2}(\eta
T_{1})-\frac{1}{2}\langle\sigma v\rangle_{D\bar{D}\to
i\bar{i}}(T_{1})Y_{D}^{2}+\langle\sigma v\rangle_{Z^{\prime}Z^{\prime}\to
D\bar{D}}(T_{1})Y^{2}_{Z^{\prime}}$ $\displaystyle-\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to
Z^{\prime}Z^{\prime}}(T_{1})Y_{D}^{2}-\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to\gamma^{\prime}\gamma^{\prime}}(T_{1})Y_{D}^{2}+\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to D\bar{D}}(\zeta
T_{1})Y_{\gamma^{\prime}}^{2}$ $\displaystyle-\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to Z^{\prime}\gamma^{\prime}}(T_{1})Y_{D}^{2}+\langle\sigma
v\rangle_{Z^{\prime}\gamma^{\prime}\to D\bar{D}}(T_{1},\zeta
T_{1})Y_{Z^{\prime}}Y_{\gamma^{\prime}},$ (3.8)
$\displaystyle\mathcal{J}_{Z^{\prime}}$ $\displaystyle=\langle\sigma
v\rangle_{i\bar{i}\to Z^{\prime}Z^{\prime}}(\eta T_{1})Y_{i}^{2}(\eta
T_{1})-\langle\sigma v\rangle_{Z^{\prime}Z^{\prime}\to
i\bar{i}}(T_{1})Y_{Z^{\prime}}^{2}-\langle\sigma
v\rangle_{Z^{\prime}Z^{\prime}\to D\bar{D}}(T_{1})Y^{2}_{Z^{\prime}}$
$\displaystyle+\langle\sigma v\rangle_{i\bar{i}\to Z^{\prime}}(\eta
T_{1})Y^{2}_{i}(\eta T_{1})+\frac{1}{2}\langle\sigma v\rangle_{D\bar{D}\to
Z^{\prime}Z^{\prime}}(T_{1})Y_{D}^{2}+\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to Z^{\prime}\gamma^{\prime}}(T_{1})Y_{D}^{2}$
$\displaystyle-\langle\sigma v\rangle_{Z^{\prime}\gamma^{\prime}\to
D\bar{D}}(T_{1},\zeta
T_{1})Y_{Z^{\prime}}Y_{\gamma^{\prime}}-\frac{1}{s}\langle\Gamma_{Z^{\prime}\to
i\bar{i}}\rangle(T_{1})Y_{Z^{\prime}},$ (3.9)
$\displaystyle\mathcal{J}_{\gamma^{\prime}}$
$\displaystyle=\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to\gamma^{\prime}\gamma^{\prime}}(T_{1})Y^{2}_{D}-\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to D\bar{D}}(\zeta
T_{1})Y^{2}_{\gamma^{\prime}}+\langle\sigma
v\rangle_{i\bar{i}\to\gamma^{\prime}}(\eta T_{1})Y_{i}^{2}(\eta T_{1})$
$\displaystyle+\langle\sigma
v\rangle_{i\bar{i}\to\gamma^{\prime}\gamma^{\prime}}(\eta T_{1})Y_{i}^{2}(\eta
T_{1})-\langle\sigma v\rangle_{\gamma^{\prime}\gamma^{\prime}\to
i\bar{i}}(\zeta
T_{1})Y_{\gamma^{\prime}}^{2}-\frac{1}{s}\langle\Gamma_{\gamma^{\prime}\to
i\bar{i}}\rangle(\zeta T_{1})Y_{\gamma^{\prime}}$
$\displaystyle+\frac{1}{2}\langle\sigma v\rangle_{D\bar{D}\to
Z^{\prime}\gamma^{\prime}}(T_{1})Y_{D}^{2}-\langle\sigma
v\rangle_{Z^{\prime}\gamma^{\prime}\to D\bar{D}}(T_{1},\zeta
T_{1})Y_{Z^{\prime}}Y_{\gamma^{\prime}}.$ (3.10)

The thermally averaged cross sections are calculated using Eq. (C.11) for a
single initial state temperature and with Eq. (C.12) in the case of two
different initial state temperatures, while the thermal averaging of the width
is given by Eq. (C.10). Note that the SM effective energy and entropy degrees
of freedom ($g^{v}_{\rm eff}$ and $h^{v}_{\rm eff}$) are read from tabulated
values [46, 47], while those pertaining to the hidden sectors $(g_{1\rm
eff},g_{2\rm eff},h_{1\rm eff},h_{2\rm eff})$ are calculated using Eqs.
(C.16), (C.17) and (C.18). Since the above Boltzmann equations are dependent
on the parameters $\eta$ and $\zeta$, then one must consider the evolution of
those parameters with temperature. We discuss the formalism in the next
section.

##  4 Temperature evolution in the dark sectors versus in the visible sector

In this section we derive the evolution equations for temperatures $T_{1}$ and
$T_{2}$ in the dark sectors and $T$ in the visible sector, i.e., $T/T_{1}$ and
$T_{2}/T_{1}$ as a function of $T$. However, for numerical integration
purposes it is found more convenient to use $T_{1}$ as the reference
temperature. Thus we are interested in deriving the evolution equations for
$d\eta/dT_{1}$ and $d\zeta/dT_{1}$ where we recall that $\eta$ and $\zeta$ are
defined so that

$\displaystyle T\equiv\eta T_{1},~{}T_{2}\equiv\zeta T_{1}.$ (4.1)

To this end we look at the equations for the energy densities in the visible
and the hidden sectors. In the analysis, we encounter the quantity
$d\rho/dT_{1}$. The main difficulty in computing the quantity $d\rho/dT_{1}$
is that $\rho$ is constituted of three parts,
$\rho=\rho_{v}+\rho_{1}+\rho_{2}$ which depends on three temperatures i.e.,
$\rho_{v}$ is controlled by $T$, $\rho_{1}$ is controlled by $T_{1}$ and
$\rho_{2}$ is controlled by $T_{2}$. Thus we need to express
$d\rho_{v}/dT_{1}$ in terms of $d\rho_{v}/dT$ and $d\rho_{2}/dT_{1}$ in terms
of $d\rho_{2}/dT_{2}$. Using the definitions of $\eta$ and $\zeta$ we can
write

$\displaystyle\frac{d\rho_{v}}{dT_{1}}=\left(\eta+T_{1}\frac{d\eta}{dT_{1}}\right)\frac{d\rho_{v}}{dT},~{}~{}~{}\text{and}~{}~{}~{}\frac{d\rho_{2}}{dT_{1}}=\left(\zeta+T_{1}\frac{d\zeta}{dT_{1}}\right)\frac{d\rho_{2}}{dT}.$
(4.2)

This means that a determination of $d\rho_{v}/dT_{1}$ and $d\rho_{2}/dT_{1}$
requires $d\eta/dT_{1}$ and $d\zeta/dT_{2}$. Next, we derive the evolution
equations for these quantities.

We note that $\rho_{v},\rho_{1},\rho_{2}$ satisfy the following evolution
equations

$\displaystyle\frac{d\rho_{v}}{dt}+4\rho_{v}H=j_{v},$ (4.3)
$\displaystyle\frac{d\rho_{1}}{dt}+4\rho_{1}H=j_{1},$
$\displaystyle\frac{d\rho_{2}}{dt}+4\rho_{2}H=j_{2},$

where $j_{v},j_{1},j_{2}$ are the corresponding sources. Instead of time we
will use temperature so we will need to convert derivatives with respect to
time to derivatives with respect to temperature. We note now that for any
given temperatures $T_{i}$ the time derivative of temperature is given by

$\displaystyle\frac{dT_{i}}{dt}=-\frac{4H\rho}{\frac{d\rho}{dT_{i}}}\,.$ (4.4)

As discussed above, we choose $T_{1}$ to be the reference temperature and the
evolution equation for $\rho_{v}$ in this case can be written as

$\displaystyle j_{v}-4\rho_{v}H$
$\displaystyle=-\frac{4H\rho}{\frac{d\rho}{dT_{1}}}\frac{d\rho_{v}}{dT_{1}}\,.$
(4.5)

From Eq. (4.5) we can deduce that

$\displaystyle\frac{d\rho_{v}}{dT_{1}}$
$\displaystyle=\frac{4\rho_{v}H-j_{v}}{4H(\rho_{1}+\rho_{2})+j_{v}}\left(\frac{d\rho_{1}}{dT_{1}}+\frac{d\rho_{2}}{dT_{1}}\right).$
(4.6)

In a similar fashion starting with the equation for $d\rho_{2}/dt$ we can
deduce

$\displaystyle\frac{d\rho_{2}}{dT_{1}}$
$\displaystyle=\frac{4\rho_{2}H-j_{2}}{4H(\rho_{1}+\rho_{v})+j_{2}}\left(\frac{d\rho_{1}}{dT_{1}}+\frac{d\rho_{v}}{dT_{1}}\right).$
(4.7)

Eqs. (4.6) and (4.7) are two coupled equations involving $d\rho_{v}/dT_{1}$
and $d\rho_{2}/dT_{1}$ which give the solution

$\displaystyle\frac{d\rho_{v}}{dT_{1}}$
$\displaystyle=\frac{(AB+A)}{(1-AB)}\frac{d\rho_{1}}{dT_{1}},~{}~{}~{}\text{and}~{}~{}~{}\frac{d\rho_{2}}{dT_{1}}=\frac{(AB+B)}{(1-AB)}\frac{d\rho_{1}}{dT_{1}},$
(4.8)

where

$A=\frac{4\rho_{v}H-j_{v}}{4H(\rho_{1}+\rho_{2})+j_{v}},~{}~{}~{}B=\frac{4\rho_{2}H-j_{2}}{4H(\rho_{v}+\rho_{1})+j_{2}}.$
(4.9)

Using Eqs. (4.2) and (4.8) one can then obtain the relations

$\displaystyle\frac{d\eta}{dT_{1}}$
$\displaystyle=-\frac{\eta}{T_{1}}+\frac{(AB+A)}{(1-AB)}\frac{d\rho_{1}/dT_{1}}{T_{1}\frac{d\rho_{v}}{dT}},~{}~{}~{}\text{and}~{}~{}~{}\frac{d\zeta}{dT_{1}}=-\frac{\zeta}{T_{1}}+\frac{(AB+B)}{(1-AB)}\frac{d\rho_{1}/dT_{1}}{T_{1}\frac{d\rho_{2}}{dT_{2}}},$
(4.10)

where

$\displaystyle\frac{d\rho}{dT_{1}}$
$\displaystyle=\left(\frac{4H\rho}{4H\rho_{1}-j_{1}}\right)\frac{d\rho_{1}}{dT_{1}},~{}~{}~{}\frac{d\rho_{v}}{dT}=\frac{\pi^{2}}{30}\left(\frac{dg_{\rm
eff}^{v}}{dT}\eta^{4}T_{1}^{4}+4g_{\rm eff}^{v}\eta^{3}T_{1}^{3}\right),$
(4.11) $\displaystyle\frac{d\rho_{1}}{dT_{1}}$
$\displaystyle=\frac{\pi^{2}}{30}\left(\frac{dg_{1\rm
eff}}{dT_{1}}T_{1}^{4}+4g_{1\rm
eff}T_{1}^{3}\right),~{}~{}~{}\frac{d\rho_{2}}{dT_{2}}=\frac{\pi^{2}}{30}\left(\frac{dg_{2\rm
eff}}{dT_{2}}\zeta^{4}T_{1}^{4}+4g_{2\rm eff}\zeta^{3}T_{1}^{3}\right).$
(4.12)

Using the fact that $j_{v}+j_{1}+j_{2}=0$, eliminating $j_{v}$ in favor of
$j_{1}$ and $j_{2}$ and inserting Eq. (4.9) in Eq. (4.10), one can further
simplify Eq. (4.10) to cast $d\eta/dT_{1}$ and $d\zeta/dT_{1}$ in their final
form as

$\displaystyle\frac{d\eta}{dT_{1}}=$
$\displaystyle-\frac{\eta}{T_{1}}+\left(\frac{4H\rho_{v}+j_{1}+j_{2}}{4H\rho_{1}-j_{1}}\right)\frac{d\rho_{1}/dT_{1}}{T_{1}\frac{d\rho_{v}}{dT}},$
(4.13) $\displaystyle\frac{d\zeta}{dT_{1}}=$
$\displaystyle-\frac{\zeta}{T_{1}}+\left(\frac{4H\rho_{2}-j_{2}}{4H\rho_{1}-j_{1}}\right)\frac{d\rho_{1}/dT_{1}}{T_{1}\frac{d\rho_{2}}{dT_{2}}}.$
(4.14)

The source terms $j_{1}$ and $j_{2}$ are given by

$\displaystyle j_{1}=\sum_{i}[$ $\displaystyle 2Y_{i}(T)^{2}J(i\bar{i}\to
D\bar{D})(T)+2Y_{i}(T)^{2}J(i\bar{i}\to
Z^{\prime}Z^{\prime})(T)+Y_{i}(T)^{2}J(i\bar{i}\to Z^{\prime})(T)$
$\displaystyle+2Y_{\gamma^{\prime}}^{2}J(\gamma^{\prime}\gamma^{\prime}\to
D\bar{D})(T_{2})-\frac{1}{2}Y_{D}^{2}J(D\bar{D}\to\gamma^{\prime}\gamma^{\prime})(T_{1})-\frac{1}{2}Y_{D}^{2}J(D\bar{D}\to
Z^{\prime}\gamma^{\prime})(T_{1})$
$\displaystyle+Y_{Z^{\prime}}Y_{\gamma^{\prime}}J(Z^{\prime}\gamma^{\prime}\to
D\bar{D})(T_{1},T_{2})]s^{2}-Y_{Z^{\prime}}J(Z^{\prime}\to i\bar{i})(T_{1})s,$
(4.15) $\displaystyle j_{2}=\sum_{i}[$ $\displaystyle
Y_{i}(T)^{2}J(i\bar{i}\to\gamma^{\prime})(T)+Y_{D}^{2}J(D\bar{D}\to\gamma^{\prime}\gamma^{\prime})(T_{1})-Y_{\gamma^{\prime}}^{2}J(\gamma^{\prime}\gamma^{\prime}\to
D\bar{D})(T_{2})$
$\displaystyle+2Y_{i}(T)^{2}J(i\bar{i}\to\gamma^{\prime}\gamma^{\prime})(T)-Y_{\gamma^{\prime}}^{2}J(\gamma^{\prime}\gamma^{\prime}\to
i\bar{i})(T_{2})+\frac{1}{2}Y_{D}^{2}J(D\bar{D}\to
Z^{\prime}\gamma^{\prime})(T_{1})$ $\displaystyle-
Y_{Z^{\prime}}Y_{\gamma^{\prime}}J(Z^{\prime}\gamma^{\prime}\to
D\bar{D})(T_{1},T_{2})]s^{2}-Y_{\gamma^{\prime}}J(\gamma^{\prime}\to\nu\bar{\nu})(T_{2})s,$
(4.16)

with

$\displaystyle n_{i}(T)^{2}J(i\bar{i}\to D\bar{D})(T)$
$\displaystyle=\frac{T}{32\pi^{4}}\int_{4m^{2}_{D}}^{\infty}ds~{}\sigma_{D\bar{D}\to
i\bar{i}}s(s-s_{i})K_{2}(\sqrt{s}/T),$ (4.17) $\displaystyle
n_{i}(T)^{2}J(i\bar{i}\to Z^{\prime})(T)$
$\displaystyle=\frac{T}{32\pi^{4}}\int_{4m^{2}_{i}}^{\infty}ds~{}\sigma_{i\bar{i}\to\gamma^{\prime}}s(s-s_{i})K_{2}(\sqrt{s}/T),$
(4.18) $\displaystyle
n_{D}(T_{1})^{2}J(D\bar{D}\to\gamma^{\prime}\gamma^{\prime})(T_{1})$
$\displaystyle=\frac{n_{D}(T_{1})^{2}}{8m^{4}_{D}T_{1}K_{2}^{2}(m_{D}/T_{1})}\int_{4m^{2}_{D}}^{\infty}ds~{}\sigma_{D\bar{D}\to\gamma^{\prime}\gamma^{\prime}}s(s-4m^{2}_{D})K_{2}(\sqrt{s}/T),$
(4.19) $\displaystyle n_{Z^{\prime}}J(Z^{\prime}\to i\bar{i})(T_{1})$
$\displaystyle=n_{Z^{\prime}}m_{Z^{\prime}}\Gamma_{Z^{\prime}\to i\bar{i}}\,.$
(4.20)

One should not confuse the variable $s$ in Eqs. (4.15) and (4.16) with the one
in Eq. (4.20). The former is the entropy density while the latter corresponds
to the Mandelstam variable.

##  5 Thermalization and dark freeze-out

In the analysis we make certain that the relic density of the dark relics is
consistent with the Planck data [27],

$\Omega h^{2}=0.1198\pm 0.0012,$ (5.1)

along with a $2\sigma$ corridor from theoretical calculations. Contribution to
the relic density arise from $\gamma^{\prime}$ and $D$ while $Z^{\prime}$
decays before BBN and is removed from the spectrum. In Table 1 we present four
benchmarks which satisfy all the experimental constraints.

Model | $m_{D}$ | $M_{1}$ | $M_{3}$ | $M_{4}$ | $\delta_{1}$ | $\delta_{2}$ | $m_{Z^{\prime}}$ | $m_{\gamma^{\prime}}$ | $\Omega h^{2}$
---|---|---|---|---|---|---|---|---|---
(a) | 1.00 | 4.50 | 0.0 | 0.43 | $4.0\times 10^{-10}$ | 0.40 | 4.90 | 0.43 | 0.124
(b) | 0.50 | 4.50 | 0.0 | 0.47 | $6.5\times 10^{-11}$ | 0.40 | 4.90 | 0.47 | 0.103
(c) | 0.05 | 4.50 | 0.0 | 0.45 | $5.6\times 10^{-12}$ | 0.40 | 4.91 | 0.45 | 0.102
(d) | 0.62 | 4.50 | -5.0 | 0.45 | $4.0\times 10^{-10}$ | 0.05 | 6.76 | 0.30 | 0.108
Table 1:  Benchmarks used in this analysis where $g_{X}=0.95$ and masses are
in MeV except $m_{D}$ which is in GeV. In the analysis of this table and in
the rest of the numerical analysis we choose $M_{2}=0$.

The relic density shown is that of $\gamma^{\prime}$ while that of $D$ is only
$\mathcal{O}(10^{-6})$ or less and thus negligible. We note here that in
addition to the particle physics interactions generating the dark photon relic
density, one could have in addition gravitational production [16, 17, 18].
However, such a production is highly dependent on the reheat temperature which
is model dependent. From Eq. (46) of Ref. [16], for the case when the dark
photon mass is $\sim 1$ MeV, one finds that the gravitational production of
dark photon would be suppressed when the reheat temperature $H_{I}<10^{11}$
GeV. So one may think of our model being valid for this restricted class of
inflationary models.

The dark photon is long-lived with a decay width to two neutrinos given by

$\displaystyle\Gamma_{\gamma^{\prime}\to\nu\bar{\nu}}$
$\displaystyle=\frac{g_{2}^{2}\delta_{1}^{2}(\delta_{2}-\sin\beta)^{2}\epsilon_{\gamma^{\prime}}^{4}}{8\pi}m_{\gamma^{\prime}}\tan^{2}\theta_{w}\,,$
(5.2)

where $\epsilon_{\gamma^{\prime}}=m_{\gamma^{\prime}}/m_{Z}$, and the partial
decay width of $\gamma^{\prime}$ to three photons reads [23, 48]

$\displaystyle\Gamma_{\gamma^{\prime}\to 3\gamma}$
$\displaystyle=\frac{17\alpha^{3}\alpha^{\prime}}{2^{7}3^{6}5^{3}\pi^{3}}\frac{m_{\gamma^{\prime}}^{9}}{m_{e}^{8}}\approx
4.70\times
10^{-8}\alpha^{3}\alpha^{\prime}\frac{m_{\gamma^{\prime}}^{9}}{m_{e}^{8}},$
(5.3)

where $\alpha=e^{2}/4\pi$, $\alpha^{\prime}=(ke)^{2}/4\pi$ and
$k=-\delta_{1}(\delta_{2}-\sin\beta)\cos\theta_{w}$. The dark photon’s
lifetime is larger than the age of the universe and this is illustrated by
model point (d) which gives $\tau_{\gamma^{\prime}\to\nu\bar{\nu}}\sim
8.4\times 10^{21}$ yrs and $\tau_{\gamma^{\prime}\to 3\gamma}\sim 5.3\times
10^{15}$ yrs.

Calculation of the relic density requires determining the yields by
numerically solving the five stiff coupled equations, Eqs. (3.1)$-$(3.3),
(4.13) and (4.14). The resulting yields for $D$, $Z^{\prime}$ and
$\gamma^{\prime}$ as a function of the hidden sector temperature $T_{1}$ are
shown in Fig. 3 for benchmarks (a) and (b). As the universe cools, the number
densities of $D$, $Z^{\prime}$ and $\gamma^{\prime}$ increase gradually. At
around $T_{1}=100$ GeV, the dark fermions $D$ start to freeze-out, and the
blue curve becomes flat around $T_{1}=0.1$ GeV. After the dark fermion
decouples, the dynamics of $Z^{\prime}$ and $\gamma^{\prime}$ is affected
mainly by SM particles freeze-in processes after $T_{1}=0.1$ GeV. However,
since $Z^{\prime}$ is unstable its density depletes to zero at $T_{1}\sim
10^{-4}$ GeV. The only particles that contribute to the relic density then are
$D$ (blue curve) and $\gamma^{\prime}$ (yellow curve). As noted the analysis
gives dark photon as the dominant component of DM. We note that the number
changing processes in hidden sector 1 are driven by
$D\bar{D}\longleftrightarrow Z^{\prime}Z^{\prime}$ owing to the sizable value
of the coupling $g_{X}$. Since $m_{Z^{\prime}}\ll m_{D}$, the reaction
$Z^{\prime}Z^{\prime}\to D\bar{D}$ shuts off early on as the temperature drops
while the reverse reaction remains active. This causes a significant drop in
$Y_{D}$ and as a consequence $Y_{Z^{\prime}}$ rises sharply as shown in Fig.
3. This is followed by a dramatic drop in $Y_{Z^{\prime}}$ due to the decay of
$Z^{\prime}$ to SM fermions.

Figure 3: The yields for the dark fermion $D$ and the dark bosons $Z^{\prime}$
and $\gamma^{\prime}$ as a function of the hidden sector temperature $T_{1}$
for benchmarks (a) (upper panel) and (b) (bottom panel). The horizontal dashed
line corresponds to the observed relic density which matches the freeze-out
yield of $\gamma^{\prime}$. Note that at dark freeze-out $Y_{D}\ll
Y_{\gamma^{\prime}}$.

Figure 4: Evolution of $\xi$ (upper panel) and $\kappa$ (bottom panel) as a
function of the visible sector temperature $T$ and that of $\zeta$ (middle
panel) as a function of $T_{1}$ for three benchmarks (a), (b) and (c) of Table
1.

The upper panel of Fig. 4 gives the evolution of $\xi=T_{1}/T$ as a function
of $T_{1}$ which shows that $\xi$ rises until it thermalizes with the visible
sector, i.e., $\xi\sim 1$. The middle panel of this figure gives the evolution
$\zeta=T_{2}/T_{1}$ as a function of $T_{1}$ while the bottom panel of Fig. 4
gives the evolution of $\kappa=T_{2}/T$ as function of $T$. We note that
$X_{2}$ does not thermalize with $X_{1}$. This happens because the energy
injection from $X_{1}$ into $X_{2}$ is not efficient enough. Consequently
$T_{2}/T\ll 1$ which also has implications for $\Delta N_{\rm eff}$ as we
explain later.

We note in passing that even though hidden sector 1 thermalizes with the
visible sector, there is a distinction between how that happens for the case
of the hidden sector versus the visible sector. In the presence of a coupling
induced either by kinetic mixing or by mass mixing, the visible sector and the
hidden sectors will eventually thermalize as long as the sectors do not
thermally decouple according to the second law of thermodynamics. This is what
happens in the top left panel of Fig. 4. However, we note that the speed with
which a dark photon thermalizes is much slower relative to a visible sector
particle such as a quark which has almost instantaneous thermalization with
the photon background. Further, the thermalization will cease once the
particles in the hidden sector fully decouple from the visible sector or from
each other as seen in the right top panel of Fig. 4. This is the case for
hidden sector 2.

In the left panel of Fig. 5 we show $n\langle\sigma v\rangle$ and the
thermally averaged $Z^{\prime}$ decay width as a function of $T_{1}$ for
benchmark (a). Also shown is the Hubble parameter $H(T_{1})$. As evident,
while $Z^{\prime}$ can enter into equilibrium with the visible sector for a
period of time, the dark photon barely does so. We indicate by arrows the
point at which the dark freeze-out of $D$ and $\gamma^{\prime}$ occurs. The
dark photons decouple earlier followed by the dark fermions which is also
evident in Fig. 3. We note that $\langle\Gamma_{Z^{\prime}}\rangle$ overtakes
$H(T_{1})$ at lower temperatures contributing to the depletion of $Z^{\prime}$
number density. It is of interest to ask how thermal equilibrium of dark
photons occurs once they are produced. Such an equilibrium can be achieved in
our model by considering a massless complex scalar field $\phi$ in the second
hidden sector with interactions with the dark photon of the type
$(\kappa\phi^{\dagger}\partial^{\mu}\phi
A_{\mu}^{\gamma^{\prime}}+\text{h.c.})$. Elastic scattering between
$\gamma^{\prime}$ and $\phi$, $\gamma^{\prime}\phi\to\gamma^{\prime}\phi$, can
keep the dark photon in local thermal equilibrium. We follow the method of
Ref. [49] to determine the temperature of kinetic decoupling of dark photons.
For $\kappa\sim 10^{-4}$, we find that kinetic decoupling occurs at 10 keV
which is much later than chemical decoupling which happens around 0.1 GeV. The
right panel of Fig. 5 shows the temperature of kinetic decoupling. It is to be
noted that the small value of $\kappa$ has a minimal effect on the relic
density of $\gamma^{\prime}$.

Figure 5: Left panel: A plot of $n\langle\sigma v\rangle$ for the dominant
processes in the hidden sector along with the Hubble parameter and the
thermally averaged decay width of $Z^{\prime}$. Right panel: the temperature
of kinetic decoupling of $\gamma^{\prime}$ for benchmarks (a) and (d). The
dark photon temperature traces that of the thermal bath before decoupling at
around 10 keV.

We note that the parameter space of $Z^{\prime}$ and $\gamma^{\prime}$ is
constrained by experiments such as BaBar, CHARM and other beam-dump
experiments as well as by astrophysical data from Supernova SN1987A and
stellar cooling. Those limits become even stronger when the dark photon is
assumed to be the dark matter particle. Thus, measurements of heating rates of
the Galactic center cold gas clouds [50], the temperature of the diffuse X-ray
background [51] as well as that of the intergalactic medium at the time of
He++ reionization [52, 53, 54] are affected by early $\gamma^{\prime}\to
3\gamma$ decays. Further constraints can be derived from energy injection
during the dark ages [55] and spectral distortion of the CMB [54]. The
presence of a long-lived sub-MeV particle species can contribute to the
relativistic number of degrees of freedom $\Delta N_{\rm eff}$ during BBN and
recombination [56]. All those constraints can exclude a sub-MeV dark photon
down to a kinetic mixing coefficient $\mathcal{O}(10^{-13})$.

In the model discussed here, the dark photon resides in a hidden sector
$X_{2}$ that does not interact directly with the visible sector. Instead, the
direct interaction is between the two hidden sectors $X_{1}$ and $X_{2}$ via
kinetic and mass mixings. Since $X_{1}$ mixes kinetically with the visible
sector, the interaction between $X_{2}$ and the visible sector becomes doubly
suppressed and all coupling will be proportional to
$\delta_{1}(\delta_{2}-\sin\beta)$. The quantity $\sin\beta$ is due to the
mass mixing between the hidden sector and such a term can impart millicharges
to the $D$ fermions. However, the coupling between the photon and the $D$
fermions is not only suppressed by $\delta_{1}\delta_{2}\sin\beta$ but also by
the mass ratio as evident from the expression of $c_{\gamma}$ in Eq. (2.8).
Therefore even for a modest value of $\sin\beta\sim 10^{-2}$, the millicharges
are very small and do not constitute a significant constraint on the model.
This doubly suppressed coupling between the dark photon and the SM can
alleviate the present constraints mainly from $\gamma^{\prime}\to 3\gamma$ as
seen in Fig. 6, which still removes a part of the parameter space of our
model. This constraint is derived from measurements of the intergalactic
diffuse photon background. Another decay channel for the dark photon is to two
neutrinos. This leads to a possible neutrino flux but experiments have not
reached the required sensitivity to probe masses in the sub-MeV region.
Experiments such as IceCube [57] have constrained only very heavy dark matter
decays. The region of the parameter space which would produce a dark photon
relic density within $2\sigma$ of the experimental value is shown in both
panels of Fig. 6 and labeled ‘Freeze-in’. For the case when a dominant
component of the relic density arises from gravitational production, the
parameter space of our model will be enlarged. The enlarged regions which are
represented by the hatched area in Fig. 6. This area accommodates for a dark
photon relic density down to $\sim 10^{-4}$.

It is argued in Ref. [51] that a dark photon with direct kinetic mixing with
the SM can only give a subdominant contribution to the relic density and that
such an observation can be dismissed if another production mechanism is in
effect. The model discussed here presents exactly this counter argument
required to produce a dominant dark photon dark matter. The main production
mechanism for the dark photon in the current analysis is not via the freeze-in
mechanism from the visible sector, $i\bar{i}\to\gamma^{\prime}$ and
$i\bar{i}\to\gamma^{\prime}\gamma^{\prime}$, because of the doubly suppressed
coupling (see left panel of Fig. 7) but rather from interactions between the
hidden sector particles. Thus, processes such as
$D\bar{D}\to\gamma^{\prime}\gamma^{\prime}$ have cross-sections proportional
to $g_{X}\delta_{2}$. As shown in Fig. 6, the sizes of $g_{X}$ and
$\delta_{2}$ are in the required ranges to produce a dark photon relic density
which dominates that of $D$. A smaller value of $g_{X}$ will reduce the dark
photon yield as expected (see right panel of Fig. 7). We note in passing that
the gauge coupling $g_{X}$ is constrained by the main annihilation channel
$D\bar{D}\to Z^{\prime}Z^{\prime}\to 4e$ from the Planck experiment [58, 59]
and for $m_{D}<1$ GeV, $g_{X}>0.1$ is excluded. However, this constraint does
not exist for our model since the relic abundance of our $D$ fermions is
negligible. The effect of the forward process
$D\bar{D}\to\gamma^{\prime}\gamma^{\prime}$ can be clearly seen in Fig. 3
where the drop in the $D$ fermions yield at a certain temperature is followed
by a rise in the yield of $\gamma^{\prime}$. It is worth mentioning that the
reverse process $\gamma^{\prime}\gamma^{\prime}\to D\bar{D}$ works on reducing
the number density of $\gamma^{\prime}$ on the expense of $D$, but this
process shuts off early on as shown in Fig. 7 allowing
$D\bar{D}\to\gamma^{\prime}\gamma^{\prime}$ to completely take over for lower
temperatures.

Figure 6: Exclusion limits from terrestrial and astrophysical experiments on
dark photon which has kinetic mixings with the SM sector. Excluded regions are
due to constraints from experiments which include electron and muon $g-2$
[60], BaBar [61], CHARM [62, 63], NA48 [64], E137 [65, 66], NA64 [67, 68],
E141 [69] and $\nu$-CAL [70, 71, 63]. The limits are obtained from darkcast
[72]. The strongest constraints on a dark photon (mass less than 1 MeV) come
from Supernova SN1987A (including a robustly excluded region and systematic
uncertainties) [73], stellar cooling [74] and from decay to $3\gamma$ on
cosmological timescales [75, 51]. The islands in pink are constraints from
BBN. The region where the freeze-in relic density is satisfied within
$2\sigma$ of the experimental constraint is shown in different shades of blue
corresponding to different choices of $m_{D}$. The hatched area represents an
enlargement to the $2\sigma$ region allowing a relic density as low as $\sim
10^{-4}$. In the upper panel, $m_{Z^{\prime}}=(5-20)m_{\gamma^{\prime}}$ and
$\delta_{2}=0.4$ while in the lower panel $m_{Z^{\prime}}=3$ MeV and
$\delta_{2}=0.01$.

Finally, we check the number of relativistic degrees of freedom generated by
the dark photon (and possibly a complex scalar) at BBN time. The SM gives
$N_{\rm eff}=3.046$. The dark photon contribution is given by

$\displaystyle\Delta N_{\rm
eff}\simeq\frac{12}{7}\left(\frac{11}{4}\right)^{\frac{4}{3}}\left(\frac{T_{2}}{T_{\gamma}}\right)^{4},$
(5.4)

where $T_{\gamma}=T$. Using the ratio of the temperatures $T_{2}/T<0.1$ from
Fig. 4, one finds that $\Delta N_{\rm eff}$ from dark photons is
$\mathcal{O}(10^{-4})$ which makes a negligible contribution to the SM $N_{\rm
eff}$. With the inclusion of a complex scalar, the model has now five new
bosonic degrees of freedom but this still gives a small contribution and does
not violate the bound on $\Delta N_{\rm eff}$. We note that the suppressed
value of $T_{2}/T$ which arises from the non-thermalization of sectors 1 and 2
is due to the low density of dark fermions in sector 1. This can be seen by
the yield for the dark fermion in Fig. 3.

Figure 7: Left panel: A plot of $n\langle\sigma v\rangle$ for the processes
contributing to the dark photon number density as a function of $T_{1}$. Right
panel: The yields of the hidden sector particles showing a diminishing
$Y_{\gamma^{\prime}}$ due to a smaller $g_{X}$ for benchmark (d).

##  6 Conclusions

In this work we discussed the possibility that DM in the universe is
constituted of sub-MeV dark photons which reside in the hidden sector. In this
case a proper analysis of the relic density requires a solution to coupled
Boltzmann equations which depend on multiple bath temperatures including the
bath temperature for the visible sector and those for all the hidden sectors
that are feebly coupled with the visible sector. In this work we discussed a
model where the visible sector couples with two hidden sectors $X_{1}$ and
$X_{2}$ and where the particles in the hidden sector consist of a dark
fermion, a dark $Z^{\prime}$ and a dark photon. The dark $Z^{\prime}$ decays
and disappears from the spectrum while the dark photon is a long lived relic.
We show that the relic density of the dark photon depends critically on
temperatures of both the visible and the hidden sectors. We present exclusion
plots where a sub-MeV dark photon can exist consistent with all the current
experimental constraints. We also show that the existence of a dark photon is
consistent with the constraints on $N_{\rm eff}$ from BBN. Thus a sub-MeV dark
photon is a viable candidate for DM within a constrained parameter space of
mass and kinetic couplings. The formalism developed here of correlated
evolution of bath temperatures in the visible and hidden sectors may find
application for a wider class of phenomena involving hidden sectors.

The research of AA was supported by the BMBF under contract 05H18PMCC1. The
research of WZF was supported in part by the National Natural Science
Foundation of China under Grant No. 11905158 and No. 11935009. The research of
PN and ZYW was supported in part by the NSF Grant PHY-1913328.

##  Appendix A Canonical normalization of extended $G_{\text{S}M}\times
U(1)_{X_{1}}\times U(1)_{X_{2}}$ Lagrangian with kinetic and Stueckelberg mass
mixings

Consider the Lagrangian

$\mathcal{L}=-\frac{1}{4}V_{\mu\nu}^{T}K_{E}V^{\mu\nu}-\frac{1}{2}V^{T}M^{2}V,$
(A.1)

where $V=(V_{1},V_{2},V_{3},V_{4})^{T}$ which we choose to be $(D,C,B,A_{3})$.
This Lagrangian can be put in the canonical form by appropriate
transformations on $K_{E}$ and $M^{2}$. For the case when there is one dark
sector it was done analytically in [33] and the basic reason which allows that
to happen is that one of the eigenvalues of $M^{2}$ is zero corresponding to
the photon which effectively reduces the analysis to two massive modes which
can be handled analytically. In the present case since we have two hidden
sectors, we have a $4\times 4$ matrix, and while one of the eigenvalues
corresponding to the photon is zero, one still has to deal with a cubic
equation which, although possible to solve analytically, quickly becomes
intractable in the presence of both kinetic and Stueckelberg mass mixing.
However, because the kinetic mixings are typically small, it is possible to
get accurate results by expanding couplings of the dark particles with the SM
particles in powers of the kinetic mixings. In this case the relevant
couplings can be recovered easily. However, such an expansion must occur
around stable minima. This means that we must first diagonalize the SM mass
squared matrix for the gauge bosons, and compute the kinetic mixing in this
basis. We can then put the kinetic term in the canonical form. This step
requires several $GL(2,\mathbb{R})$ transformations because of several mixings
of the hidden sector with the visible sector and the mixing of the two hidden
sectors. After the kinetic energy is put in the canonical form, we must write
the mass square matrix of the gauge bosons in the same basis which then
undiagonalizes the said matrix. However, because of the smallness of the
kinetic mixings, we can carry out a perturbation expansion of the mass square
matrix where the zeroth order mass square matrix is diagonal and the
perturbations are proportional to the kinetic mixings and are small. We make
this analysis concrete in the formalism below.

Let us consider an orthogonal transformation $V=RV^{(1)}$ such that
$R^{T}M^{2}R=M_{D}^{2}$ where $M^{2}_{D}$ is a diagonal matrix. In the
$V^{(1)}$ basis the kinetic energy has the form $K_{E}^{\prime}=R^{T}K_{E}R$.
Next, let us make a transformation $K$ such that $V^{(1)}=KV^{(n)}$, which
could be a product of several sub-transformations, such that the kinetic
energy is in the canonical form, i.e.,

$K^{T}K_{E}^{\prime}K=\mathbb{1}.$ (A.2)

In our case we will have $n=8$ (as discussed below). In the $V^{(8)}$ basis,
while the kinetic energy is in the canonical form, the mass matrix
$M^{2}=K^{T}M^{2}_{D}K$ is not. However, as explained above since the kinetic
mixings are small we can expand $M^{2}$ around $\delta_{1}=0=\delta_{2}$ so
that

$\displaystyle K^{T}M_{D}^{2}K$ $\displaystyle=M_{D}^{2}+\Delta M^{2}.$ (A.3)

Now since the kinetic mixing is supposed to be small, $K$ differs from a unit
matrix only by a small amount and thus $\Delta M^{2}$ is small relative to
$M_{D}^{2}$ and one may carry out perturbation expansion in $\Delta M^{2}$ to
arrive at the kinetic and mass mixing effects in the physical processes. To
compute $\Delta M^{2}$ we need $K$ defined by Eq. (A.2). The computation of
$K$ is significantly more complicated than for the case of one hidden sector.
Below we give its computation in some detail.

While the procedure outlined above is general we will discuss the specific
case where $M^{2}$ is block diagonal so that

$M^{2}=\begin{pmatrix}M_{4}^{2}&M_{3}M_{4}&0&0\\\
M_{3}M_{4}&M_{3}^{2}+M_{1}^{2}&0&0\\\
0&0&\frac{1}{4}g_{Y}^{2}v^{2}&-\frac{1}{4}g_{Y}g_{2}v^{2}\\\
0&0&-\frac{1}{4}g_{Y}g_{2}v^{2}&\frac{1}{4}g_{2}^{2}v^{2}\\\ \end{pmatrix},$
(A.4)

where the upper right $2\times 2$ matrix is for the hidden sector and the
lower left $2\times 2$ matrix is for the case of the standard model in the
basis $V^{T}=(D_{\mu},C_{\mu},B_{\mu},A_{3\mu})$. $M^{2}$ can be diagonalized
by $R$ where

$R=\begin{pmatrix}R_{\beta}&0\\\ 0&R_{w}\\\ \end{pmatrix},\quad
R_{\beta}=\begin{pmatrix}\cos\beta&-\sin\beta\\\ \sin\beta&\cos\beta\\\
\end{pmatrix},~{}~{}~{}~{}R_{w}=\begin{pmatrix}\cos\theta_{w}&-\sin\theta_{w}\\\
\sin\theta_{w}&\cos\theta_{w}\\\ \end{pmatrix},$ (A.5)

with $\theta_{w}$ being the weak mixing angle. Here $R_{\beta}$ diagonalizes
the hidden sector mass squared matrix while $R_{w}$ diagonalizes the standard
model mass squared matrix, and the diagonalization gives

$R^{T}M^{2}R=M_{D}^{2}\equiv\text{diag}(m_{\gamma^{\prime}}^{2},m_{Z^{\prime}}^{2},0,m_{Z}^{2}).$
(A.6)

Since the mass square matrix is now diagonal, it is a good starting point to
diagonalize and normalize the kinetic energy matrix. This is a bit non-trivial
and requires several steps which we outline below. In the basis
$(D,C,B,A_{3})$, the kinetic Lagrangian is given by

$\mathcal{L}_{\rm
KE}=-\frac{1}{4}(D^{2}+C^{2}+B^{2}+A_{3}^{2})-\frac{1}{4}(2\delta_{1}BC+2\delta_{2}CD),$
(A.7)

where we use an abbreviated notation so that $B^{2}=B_{\mu\nu}B^{\mu\nu}$,
$BC=B_{\mu\nu}C^{\mu\nu}$, etc. Next we write $\mathcal{L}_{\rm KE}$ in the
basis $(D^{(1)},C^{(1)},B^{(0)},A_{3}^{(0)})$ in which the mass square matrix
of the gauge bosons is diagonal. Thus, using

$\begin{pmatrix}D\\\ C\\\ \end{pmatrix}=R_{\beta}\begin{pmatrix}D^{(1)}\\\
C^{(1)}\\\ \end{pmatrix},\quad\begin{pmatrix}B\\\ A_{3}\\\
\end{pmatrix}=R_{w}\begin{pmatrix}B^{(0)}\\\ A_{3}^{(0)}\\\ \end{pmatrix},$
(A.8)

allows us to write $\mathcal{L}_{\rm KE}$ as

$\displaystyle\mathcal{L}_{\rm KE}=$
$\displaystyle-\frac{1}{4}(D^{(1)2}+C^{(1)2}+B^{(0)2}+A_{3}^{(0)2})$ (A.9)
$\displaystyle-\frac{1}{2}\delta_{1}(B^{(0)}\cos\theta_{w}-A_{3}^{(0)}\sin\theta_{w})(D^{(1)}\sin\beta+C^{(1)}\cos\beta)$
$\displaystyle-\frac{1}{2}\delta_{2}(D^{(1)}\sin\beta+C^{(1)}\cos\beta)(D^{(1)}\cos\beta-C^{(1)}\sin\beta).$

The diagonal kinetic terms for $D^{(1)}$ and $C^{(1)}$ have the form

$-\frac{1}{4}D^{(1)^{2}}(1+\delta_{2}\sin
2\beta)-\frac{1}{4}C^{(1)^{2}}(1-\delta_{2}\sin 2\beta).$ (A.10)

To normalize them to unity we make a transformation from the basis
$V^{(1)^{T}}=(D^{(1)},C^{(1)},B^{(0)},A_{3}^{(0)})$ to
$V^{(2)^{T}}=(D^{(2)},C^{(2)},B^{(0)},A_{3}^{(0)})$ so that

$V^{(1)}=K_{1}V^{(2)},~{}~{}K_{1}=\begin{pmatrix}\frac{1}{\sqrt{1+\delta_{2}\sin
2\beta}}&0&0&0\\\ 0&\frac{1}{\sqrt{1-\delta_{2}\sin 2\beta}}&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix},$ (A.11)

where

$\bar{\delta_{2}}=\frac{\delta_{2}\cos
2\beta}{\sqrt{1-\delta_{2}^{2}\sin^{2}2\beta}}.$ (A.12)

After the transformation, $\mathcal{L}_{\rm KE}$ in the $V^{(2)}$ basis has
the form

$\displaystyle\mathcal{L}_{\rm KE}=$
$\displaystyle-\frac{1}{4}(D^{(2)^{2}}+C^{(2)^{2}}+B^{(0)^{2}}+A_{3}^{(0)^{2}})-\frac{1}{2}\bar{\delta_{2}}C^{(2)}D^{(2)}$
(A.13)
$\displaystyle-\frac{1}{2}(\delta_{1}^{+}D^{(2)}+\delta_{1}^{-}C^{(2)})(B^{(0)}\cos\theta_{w}-A_{3}^{(0)}\sin\theta_{w}),$

where

$\delta_{1}^{+}=\frac{\delta_{1}\sin\beta}{\sqrt{1+\delta_{2}\sin
2\beta}},~{}~{}~{}\delta_{1}^{-}=\frac{\delta_{1}\cos\beta}{\sqrt{1-\delta_{2}\sin
2\beta}}.$ (A.14)

We now note that there is a $C^{(2)}D^{(2)}$ mixing term in Eq. (A.13) which
can be removed by a $GL(2,\mathbb{R})$ transformation. We do this by going
from the basis $V^{(2)}$ to $V^{(3)T}=(D^{(3)},C^{(3)},B^{(0)},A_{3}^{(0)})$
so that

$\displaystyle V^{(2)}=K_{2}V^{(3)},\quad
K_{2}=\begin{pmatrix}1&-s_{\bar{\delta_{2}}}&0&0\\\
0&c_{\bar{\delta_{2}}}&0&0\\\ 0&0&1&0\\\ 0&0&0&1\end{pmatrix},$ (A.15)

where

$s_{\bar{\delta_{2}}}=\frac{\bar{\delta_{2}}}{\sqrt{1-\bar{\delta_{2}}^{2}}},~{}~{}~{}c_{\bar{\delta_{2}}}=\frac{1}{\sqrt{1-\bar{\delta_{2}}^{2}}}.$
(A.16)

In the $V^{(3)}$ basis, $\mathcal{L}_{\rm KE}$ becomes

$\displaystyle\mathcal{L}_{\rm
KE}=-\frac{1}{4}(D^{(3)^{2}}+C^{(3)^{2}}+B^{(0)^{2}}+A_{3}^{(0)^{2}})$
$\displaystyle-\frac{1}{2}\delta_{1}^{-}C^{(3)}(B^{(0)}\cos\theta_{w}-A_{3}^{(0)}\sin\theta_{w})$
(A.17)
$\displaystyle-\frac{1}{2}\delta^{\prime}D^{(3)}(B^{(0)}\cos\theta_{w}-A_{3}^{(0)}\sin\theta_{w}),$

where

$\displaystyle\delta^{\prime}\equiv\delta_{1}^{+}c_{\bar{\delta_{2}}}-\delta_{1}^{-}s_{\bar{\delta_{2}}}.$
(A.18)

We note that while there are no kinetic mixing terms between $C^{(3)}$ and
$D^{(3)}$, there are kinetic mixing terms between them and the fields
$B^{(0)}$ and $A_{3}^{(0)}$. The mixing term between $D^{(3)}$ and $B^{(0)}$
can be removed by the transformation

$\displaystyle
V^{(3)}=K_{3}V^{(4)},~{}~{}K_{3}=\begin{pmatrix}1&0&-s_{\delta_{3}}&0\\\
0&1&0&0\\\ 0&0&c_{\delta_{3}}&0\\\ 0&0&0&1\end{pmatrix},$ (A.19)

where $V^{(4)^{T}}=(D^{(4)},C^{(3)},B^{(1)},A_{3}^{(0)})$, $s_{\delta_{3}}$
and $c_{\delta_{3}}$ are defined similar to Eq. (A.16) and $\delta_{3}$ is
defined by

$\delta_{3}\equiv\delta^{\prime}\cos\theta_{w}\,.$ (A.20)

In the $V^{(4)}$ basis the Lagrangian takes the form

$\displaystyle\mathcal{L}_{\rm KE}=$
$\displaystyle-\frac{1}{4}(D^{(4)^{2}}+C^{(3)^{2}}+B^{(1)^{2}}+A_{3}^{(0)^{2}})-\frac{1}{2}\delta_{1}^{-}C^{(3)}(B^{(1)}c_{\delta_{3}}\cos\theta_{w}-A_{3}^{(0)}\sin\theta_{w})$
(A.21)
$\displaystyle-\frac{1}{2}\delta^{\prime}D^{(4)}(B^{(1)}c_{\delta_{2}}\cos\theta_{w}-A_{3}^{(0)}\sin\theta_{w})+\frac{1}{2}\delta^{\prime}s_{\delta_{3}}B^{(1)}(-B^{(1)}s_{\delta_{3}}\cos\theta_{w}-A_{3}^{(0)}\sin\theta_{w}).$

A mixing term between $D^{(4)}$ and $A_{3}^{(0)}$ exists which can be removed
by the transformation

$\displaystyle
V^{(4)}=K_{4}V^{(5)},~{}~{}K_{4}=\begin{pmatrix}1&0&0&-s_{\delta_{4}}\\\
0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&c_{\delta_{4}}\end{pmatrix},$ (A.22)

where $V^{(5)}$ is given by
$V^{(5)^{T}}=(D^{(5)},C^{(3)},B^{(1)},A_{3}^{(1)})$ and where $s_{\delta_{4}}$
and $c_{\delta_{4}}$ are defined as in Eq. (A.16) and

$\delta_{4}\equiv-\delta^{\prime}\sin\theta_{w}\,.$ (A.23)

In the $V^{(5)}$ basis, $\mathcal{L}_{\rm KE}$ takes the form

$\displaystyle\mathcal{L}_{\rm
KE}=-\frac{1}{4}(D^{(5)^{2}}+C^{(3)^{2}}+B^{(1)^{2}}+A_{3}^{(1)^{2}})$
$\displaystyle-\frac{1}{2}(B^{(1)}c_{\delta_{3}}\cos\theta_{w}-A_{3}^{(1)}c_{\delta_{4}}\sin\theta_{w})\delta_{1}^{-}C^{(3)}$
$\displaystyle-\frac{1}{2}\sin\theta_{w}\delta^{\prime}\delta_{3}c_{\delta_{4}}A_{3}^{(1)}B^{(1)}.$
(A.24)

Next we look at the kinetic mixing of $(C^{(3)},B^{(1)})$. This mixing can be
removed by the transformation

$\displaystyle V^{(5)}=K_{5}V^{6},~{}~{}~{}K_{5}=\begin{pmatrix}1&0&0&0\\\
0&1&-s_{\delta_{5}}&0\\\ 0&0&c_{\delta_{5}}&0\\\ 0&0&0&1\end{pmatrix}.$ (A.25)

Here $V^{(6)^{T}}=(D^{(5)},C^{(4)},B^{(2)},A_{3}^{(1)})$ and

$\delta_{5}=\delta_{1}^{-}c_{\delta_{3}}\cos\theta_{w}\,,$ (A.26)

where $s_{\delta_{5}}$ and $c_{\delta_{5}}$ are defined as in Eq. (A.16).
After the transformation the kinetic energy Lagrangian in the $V^{(5)}$ basis
has the form

$\displaystyle\mathcal{L}_{\rm KE}=$
$\displaystyle-\frac{1}{4}(D^{(5)^{2}}+C^{(4)^{2}}+B^{(2)^{2}}+A_{3}^{(1)^{2}})+\frac{1}{2}\sin\theta_{w}\delta_{1}^{-}c_{\delta_{4}}A_{3}^{(1)}(C^{(4)}-s_{\delta_{5}}B^{(2)})$
(A.27)
$\displaystyle-\frac{1}{2}\sin\theta_{w}\delta^{\prime}\delta_{3}c_{\delta_{4}}c_{\delta_{5}}A_{3}^{(1)}B^{(2)}.$

Next we examine the kinetic mixing of the fields $(C^{(4)},A_{3}^{(1)})$. This
mixing term can be eliminated by the transformation

$\displaystyle V^{(6)}=K_{6}V^{(7)},~{}~{}~{}K_{6}=\begin{pmatrix}1&0&0&0\\\
0&1&0&-s_{\delta_{6}}\\\ 0&0&1&0\\\ 0&0&0&c_{\delta_{6}}\end{pmatrix},$ (A.28)

where $V^{(7)^{T}}=(D^{(5)},C^{(5)},B^{(2)},A_{3}^{(2)})$ and where
$\delta_{6}$ is defined by

$\delta_{6}=-\delta_{1}^{-}c_{\delta_{4}}\sin\theta_{w},$ (A.29)

and $s_{\delta_{6}}$ and $c_{\delta_{6}}$ are defined as usual.

After the transformation the kinetic energy Lagrangian in the $V^{(7)}$ basis
has the form

$\mathcal{L}_{\rm
KE}=-\frac{1}{4}(D^{(5)^{2}}+C^{(5)^{2}}+B^{(2)^{2}}+A_{3}^{(2)^{2}})+\frac{1}{2}\delta_{7}B^{(2)}A_{3}^{(2)},$
(A.30)

where

$\delta_{7}=\sin\theta_{w}\delta_{1}^{-}c_{\delta_{4}}c_{\delta_{6}}s_{\delta_{5}}+\sin\theta_{w}\delta^{\prime}c_{\delta_{4}}c_{\delta_{6}}c_{\delta_{5}}s_{\delta_{3}}.$
(A.31)

We are now left with the last kinetic mixing term involving the fields
$B^{(2)}$ and $A_{3}^{(2)}$. To eliminate this mixing we make the final
transformation

$\displaystyle V^{(7)}=K_{7}V^{(8)},\quad K_{7}=\begin{pmatrix}1&0&0&0\\\
0&1&0&0\\\ 0&0&1&-s_{\delta_{7}}\\\ 0&0&0&c_{\delta_{7}}\end{pmatrix},$ (A.32)

where $V^{(8)^{T}}=(D^{(5)},C^{(5)},B^{(3)},A_{3}^{(3)})$. In the basis
$V^{(8)}$ the kinetic energy for all the gauge fields is in the canonical form
so that

$\mathcal{L}_{\rm
KE}=-\frac{1}{4}(D^{(5)^{2}}+C^{(5)^{2}}+B^{(3)^{2}}+A_{3}^{(3)^{2}}).$ (A.33)

The free Lagrangian in the $V^{(8)}$ basis is then

$\displaystyle\mathcal{L}$
$\displaystyle=-\frac{1}{4}V^{(8)^{T}}V^{(8)}-\frac{1}{2}V^{(8)^{T}}K^{T}M_{D}^{2}KV^{(8)},$
(A.34)

with

$K\equiv K_{1}K_{2}K_{3}K_{4}K_{5}K_{6}K_{7}.$ (A.35)

As discussed in the beginning of this section we now make the expansion of Eq.
(A.3). Below we exhibit $K$ and $\Delta M^{2}$ in the limit
$\delta_{1},\delta_{2}\ll 1$. In this case
$s_{\delta_{1}}\sim\delta_{1},c_{\delta_{1}}\sim 1$, etc. and $K$ and $\Delta
M^{2}$ have the following form

$K\sim\begin{pmatrix}1&-\bar{\delta_{2}}&-\delta_{3}&-\delta_{4}\\\
0&1&-\delta_{5}&-\delta_{6}\\\ 0&0&1&-\delta_{7}\\\
0&0&0&1\end{pmatrix},~{}~{}\Delta
M^{2}=\begin{pmatrix}0&-m_{\gamma^{\prime}}^{2}\bar{\delta_{2}}&-m_{\gamma^{\prime}}^{2}\delta_{3}&-m_{\gamma^{\prime}}^{2}\delta_{4}\\\
-m_{\gamma^{\prime}}^{2}\bar{\delta_{2}}&0&-m_{Z^{\prime}}^{2}\delta_{5}&-m_{\gamma^{\prime}}^{2}\delta_{6}\\\
-m_{\gamma^{\prime}}^{2}\delta_{3}&-m_{Z^{\prime}}^{2}\delta_{5}&0&0\\\
-m_{\gamma^{\prime}}^{2}\delta_{4}&-m_{\gamma^{\prime}}^{2}\delta_{6}&0&0\end{pmatrix}.$
(A.36)

The interactions relevant for our computation arise from $\Delta M^{2}$ and
the relation

$\displaystyle V^{(1)}=KV^{(8)}.$ (A.37)

Eqs. (A.36) and (A.37) and non-degenerate perturbation theory is utilized in
the computation of couplings of the visible sector with the hidden sector.
This is discussed in the next section.

##  Appendix B Dark photon $\gamma^{\prime}$ and dark $Z^{\prime}$ couplings
with Standard Model particles

The couplings of the dark photon and $Z^{\prime}$ to the SM particles are
given by the Lagrangian Eq. (2.4). To compute the couplings proportional to
$\delta_{1}$ and $\delta_{2}$ that arise due to the kinetic and Stueckelberg
mass mixings, we use first order non-degenerate perturbation theory using
$\Delta M^{2}$ given in Eq. (A.36) as the perturbation. Thus, to first order
perturbation in $\Delta M^{2}$, the neutral currents of Eq. (2.4) which
involve the vector and axial-vector couplings of the dark photon, the dark
$Z^{\prime}$ and the SM gauge gauge bosons are given by

$\displaystyle v_{f}$ $\displaystyle=T_{3f}-2Q_{f}\sin^{2}\theta_{w},$ (B.1)
$\displaystyle a_{f}$ $\displaystyle=T_{3f},$ (B.2) $\displaystyle
v^{\prime}_{f}$ $\displaystyle=-Q_{f}\sin
2\theta_{w}\cos\theta_{w}(1+\epsilon_{z}^{2})\delta_{1}\left[1-\left(1-\frac{T_{3f}}{2Q}\right)\frac{m_{Z^{\prime}}^{2}}{m^{2}_{W}}\right],$
(B.3) $\displaystyle a^{\prime}_{f}$
$\displaystyle=-\delta_{1}T_{3f}\sin\theta_{w}\epsilon^{2}_{z}(1+\epsilon^{2}_{z}),$
(B.4) $\displaystyle v^{\prime\prime}_{f}$ $\displaystyle=Q_{f}\sin
2\theta_{w}\cos\theta_{w}(1+\epsilon_{\gamma^{\prime}}^{2})\left[1-\left(1-\frac{T_{3f}}{2Q}\right)\frac{m_{\gamma^{\prime}}^{2}}{m^{2}_{W}}\right]\delta_{1}(\delta_{2}-\sin\beta),$
(B.5) $\displaystyle a^{\prime\prime}_{f}$
$\displaystyle=T_{3f}\sin\theta_{w}\epsilon_{\gamma^{\prime}}^{2}(1+\epsilon_{\gamma^{\prime}}^{2})\delta_{1}(\delta_{2}-\sin\beta),$
(B.6)

where $\epsilon_{z}=m_{Z^{\prime}}/m_{Z}$ and
$\epsilon_{\gamma^{\prime}}=m_{\gamma^{\prime}}/m_{Z}$. The relevant couplings
with the visible sector are summarized in Tables 2 and 3. Thus Table 2 gives
the couplings of $Z$ and $Z^{\prime}$ to the visible sector fermions
$f\bar{f}$ and Table 3. gives the coupling of $\gamma^{\prime}$ to $f\bar{f}$.

$f$ | $Q_{f}$ | $v_{f}$ | $a_{f}$ | $v^{\prime}_{f}$ | $a^{\prime}_{f}$
---|---|---|---|---|---
$\nu_{e},\nu_{\mu},\nu_{\tau}$ | 0 | $\frac{1}{2}$ | $\frac{1}{2}$ | $-\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{z}(1+\epsilon^{2}_{z})\delta_{1}$ | $-\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{z}(1+\epsilon^{2}_{z})\delta_{1}$
$e,\mu,\tau$ | $-1$ | $-\frac{1}{2}+2\sin^{2}\theta_{w}$ | $-\frac{1}{2}$ | $\sin 2\theta_{w}\cos\theta_{w}\left(1-\frac{3m^{2}_{Z^{\prime}}}{4m^{2}_{W}}\right)(1+\epsilon^{2}_{z})\delta_{1}$ | $\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{z}(1+\epsilon^{2}_{z})\delta_{1}$
$u,c,t$ | $\frac{2}{3}$ | $\frac{1}{2}-\frac{4}{3}\sin^{2}\theta_{w}$ | $\frac{1}{2}$ | $-\frac{2}{3}\sin 2\theta_{w}\cos\theta_{w}\left(1-\frac{5m^{2}_{Z^{\prime}}}{8m^{2}_{W}}\right)(1+\epsilon^{2}_{z})\delta_{1}$ | $-\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{z}(1+\epsilon^{2}_{z})\delta_{1}$
$d,s,b$ | $-\frac{1}{3}$ | $-\frac{1}{2}+\frac{2}{3}\sin^{2}\theta_{w}$ | $-\frac{1}{2}$ | $\frac{1}{3}\sin 2\theta_{w}\cos\theta_{w}\left(1-\frac{m^{2}_{Z^{\prime}}}{4m^{2}_{W}}\right)(1+\epsilon^{2}_{z})\delta_{1}$ | $\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{z}(1+\epsilon^{2}_{z})\delta_{1}$

Table 2:  The $Z\to f\bar{f}$ and $Z^{\prime}\to f\bar{f}$ vertices. In the
above, $\epsilon_{z}=m_{Z^{\prime}}/m_{Z}$ and $m_{W}$ is the $W$ boson mass.

$f$ | $Q_{f}$ | $v^{\prime\prime}_{f}$ | $a^{\prime\prime}_{f}$
---|---|---|---
$\nu_{e},\nu_{\mu},\nu_{\tau}$ | 0 | $\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{\gamma^{\prime}}(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(\delta_{2}-s_{\beta})$ | $\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{\gamma^{\prime}}(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(\delta_{2}-s_{\beta})$
$e,\mu,\tau$ | $-1$ | $\sin 2\theta_{w}\cos\theta_{w}\left(1-\frac{3m^{2}_{\gamma^{\prime}}}{4m^{2}_{W}}\right)(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(s_{\beta}-\delta_{2})$ | $-\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{\gamma^{\prime}}(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(\delta_{2}-s_{\beta})$
$u,c,t$ | $\frac{2}{3}$ | $\frac{2}{3}\sin 2\theta_{w}\cos\theta_{w}\left(1-\frac{5m^{2}_{\gamma^{\prime}}}{8m^{2}_{W}}\right)(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(\delta_{2}-s_{\beta})$ | $\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{\gamma^{\prime}}(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(\delta_{2}-s_{\beta})$
$d,s,b$ | $-\frac{1}{3}$ | $\frac{1}{3}\sin 2\theta_{w}\cos\theta_{w}\left(1-\frac{m^{2}_{\gamma^{\prime}}}{4m^{2}_{W}}\right)(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(s_{\beta}-\delta_{2})$ | $-\frac{1}{2}\sin\theta_{w}\epsilon^{2}_{\gamma^{\prime}}(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(\delta_{2}-s_{\beta})$

Table 3:  The $\gamma^{\prime}\to f\bar{f}$ vertices. In the above,
$\epsilon_{\gamma^{\prime}}=m_{\gamma^{\prime}}/m_{Z}$.

The triple gauge boson couplings of $\gamma^{\prime},Z^{\prime},Z$ are given
by are

$\displaystyle
WW\gamma^{\prime}:~{}-ig_{2}\cos\theta_{w}\sin\theta_{w}(1+\epsilon^{2}_{\gamma^{\prime}})\delta_{1}(\delta_{2}-\sin\beta).$
(B.7) $\displaystyle
WWZ^{\prime}:~{}ig_{2}\cos\theta_{w}\sin\theta_{w}(1+\epsilon^{2}_{z})\delta_{1},$
(B.8) $\displaystyle WWZ:~{}-ig_{2}\cos\theta_{w}.$ (B.9)

Couplings in the limit of large $\delta_{2}$:  Some of the processes such as
the lifetime of the dark photon require only that the product
$\delta_{1}\delta_{2}$ be small which could be achieved by $\delta_{1}$ being
small while $\delta_{2}$ is $\mathcal{O}(1)$ size. Thus we list below the
vector and axial-vector couplings in the limit of small $\delta_{1}$ and
$\beta$ while $\delta_{2}$ is not necessarily small

$\displaystyle v_{f}$
$\displaystyle=T_{3f}-2Q_{f}\sin^{2}\theta_{w}+\frac{Q_{f}\delta_{1}\delta_{2}\sin\theta_{w}\sin
2\theta_{w}(\delta_{2}-\sin\beta+\delta_{2}^{2}\sin\beta)}{\sqrt{1-\delta_{2}^{2}}},$
(B.10) $\displaystyle a_{f}$ $\displaystyle=T_{3f},$ (B.11) $\displaystyle
v^{\prime}_{f}$ $\displaystyle=-Q_{f}\sin
2\theta_{w}\cos\theta_{w}(1+\epsilon_{z}^{2}+\epsilon^{2}_{\gamma^{\prime}}\delta_{2}^{2})\delta_{1}\left[\frac{m^{2}_{Z^{\prime}}}{m^{2}_{Z^{\prime}}+m^{2}_{\gamma^{\prime}}\delta_{2}^{2}}-\left(1-\frac{T_{3f}}{2Q_{f}}\right)\frac{m_{Z^{\prime}}^{2}}{m^{2}_{W}}\right]$
$\displaystyle\hskip
113.81102pt\times\left[(1+\delta_{2}\sin\beta)+\frac{m^{2}_{\gamma^{\prime}}}{m^{2}_{Z^{\prime}}}\frac{\delta_{2}^{2}(1+\delta_{2}\sin\beta)-\delta_{2}\sin\beta}{\sqrt{1-\delta_{2}^{2}}}\right],$
(B.12) $\displaystyle a^{\prime}_{f}$
$\displaystyle=-\delta_{1}T_{3f}\sin\theta_{w}(1+\epsilon^{2}_{z}+\delta_{2}^{2}\epsilon^{2}_{\gamma^{\prime}})\left[\epsilon^{2}_{z}(1+\delta_{2}\sin\beta)+\epsilon^{2}_{\gamma^{\prime}}\frac{\delta_{2}^{2}(1+\delta_{2}\sin\beta)-\delta_{2}\sin\beta}{\sqrt{1-\delta_{2}^{2}}}\right],$
(B.13) $\displaystyle v^{\prime\prime}_{f}$ $\displaystyle=Q_{f}\sin
2\theta_{w}\cos\theta_{w}(1+\epsilon_{\gamma^{\prime}}^{2})\left[1-\left(1-\frac{T_{3f}}{2Q_{f}}\right)\frac{m_{\gamma^{\prime}}^{2}}{m^{2}_{W}}\right]\frac{\delta_{1}(\delta_{2}-\sin\beta+\delta_{2}^{2}\sin\beta)}{\sqrt{1-\delta_{2}^{2}}},$
(B.14) $\displaystyle a^{\prime\prime}_{f}$
$\displaystyle=T_{3f}\sin\theta_{w}\epsilon_{\gamma^{\prime}}^{2}(1+\epsilon_{\gamma^{\prime}}^{2})\frac{\delta_{1}(\delta_{2}-\sin\beta+\delta_{2}^{2}\sin\beta)}{\sqrt{1-\delta_{2}^{2}}}.$
(B.15)

The couplings with the $D$ fermions become

$\displaystyle\gamma^{\prime}D\bar{D}:~{}g_{X}\dfrac{m^{2}_{\gamma^{\prime}}\delta_{2}}{m^{2}_{Z^{\prime}}-m^{2}_{\gamma^{\prime}}(1-\delta_{2}^{2})},$
(B.16) $\displaystyle Z^{\prime}D\bar{D}:~{}g_{X},$ (B.17)
$\displaystyle\gamma
D\bar{D}:~{}-g_{X}\dfrac{m^{2}_{\gamma^{\prime}}\delta_{1}\delta_{2}\cos\theta_{w}}{m^{2}_{Z^{\prime}}+m^{2}_{\gamma^{\prime}}\delta_{2}^{2}}\left[\delta_{2}(1+\beta\delta_{2})+\dfrac{\beta-\delta_{2}-\beta\delta_{2}^{2}}{\sqrt{1-\delta_{2}^{2}}}\right],$
(B.18) $\displaystyle
ZD\bar{D}:~{}g_{X}\delta_{1}\sin\theta_{w}(1+\epsilon_{z}^{2}+\delta_{2}^{2}\epsilon^{2}_{\gamma^{\prime}})\left[(1+\beta\delta_{2})(1+\delta_{2}^{2}\epsilon^{2}_{\gamma^{\prime}})+\dfrac{\delta_{2}\epsilon^{2}_{\gamma^{\prime}}}{\sqrt{1-\delta_{2}^{2}}}(\beta-\delta_{2}-\beta\delta_{2}^{2})\right],$
(B.19)

and the triple gauge boson couplings take the form

$\displaystyle
WW\gamma^{\prime}:~{}-ig_{2}\cos\theta_{w}\sin\theta_{w}(1+\epsilon^{2}_{\gamma^{\prime}})\frac{\delta_{1}(\delta_{2}-s_{\beta}+\beta\delta_{2}^{2})}{\sqrt{1-\delta_{2}^{2}}},$
(B.20) $\displaystyle
WWZ^{\prime}:~{}ig_{2}\cos\theta_{w}\sin\theta_{w}\frac{\delta_{1}(1+\epsilon^{2}_{z}+\epsilon^{2}_{\gamma^{\prime}}\delta_{2}^{2})}{m^{2}_{Z^{\prime}}+m^{2}_{\gamma^{\prime}}\delta_{2}^{2}}\left[m^{2}_{Z^{\prime}}(1+\delta_{2}\beta)+\frac{m^{2}_{\gamma^{\prime}}\delta_{2}(\beta\delta_{2}^{2}+\delta_{2}-\beta)}{\sqrt{1-\delta_{2}^{2}}}\right],$
(B.21) $\displaystyle
WWZ:~{}-ig_{2}\cos\theta_{w}\left[1+\frac{\sin\theta_{w}\tan\theta_{w}\delta_{1}\delta_{2}(\beta\delta_{2}^{2}+\delta_{2}-\beta)}{\sqrt{1-\delta_{2}^{2}}}\right].$
(B.22)

##  Appendix C Deduction of three temperature Boltzmann equations

Next we give a deduction of the Boltzmann equations for the case of three heat
baths. We will use $T_{1}$ as the reference temperature. Let us consider a
generic number density $n_{i}$. In this case $n_{i}R^{3}$ is conserved during
the expansion if there is no injection and we have
$\frac{d(n_{i}R^{3})}{dt}=0$ where $R$ is the scale factor, while in the
presence of injection one has

$\displaystyle\frac{dn_{i}}{dt}+3Hn_{i}=C_{i}\,,$ (C.1)

where $C_{i}$ represent the integrated collision terms. Next, if $S=sR^{3}$ is
the total entropy, it is conserved which implies that

$\displaystyle\frac{ds}{dt}+3Hs=0.$ (C.2)

Using Eqs. (C.1) and (C.2), and the fact that $n_{i}=sY_{i}$, one finds

$\displaystyle\frac{dY_{i}}{dt}=\frac{1}{s}C_{i}\,.$ (C.3)

We can convert this equation to one that uses temperature $T_{1}$ rather than
time which gives

$\displaystyle\frac{dY_{i}}{dT_{1}}$
$\displaystyle=-\frac{d\rho/dT_{1}}{4H\rho}\frac{1}{s}C_{i}\,,$ (C.4)

where $d\rho/dT_{1}$ is given by

$\displaystyle\frac{d\rho}{dT_{1}}$
$\displaystyle=(\eta+T_{1}\eta^{\prime})\frac{d\rho_{v}}{dT}+\frac{d\rho_{1}}{dT_{1}}+(\zeta+T_{1}\zeta^{\prime})\frac{d\rho_{2}}{dT_{2}}.$
(C.5)

The Boltzmann equations for $Y_{D},Y_{Z^{\prime}},Y_{\gamma^{\prime}}$ may now
be written as

$\displaystyle\frac{dY_{D}}{dT_{1}}=-\frac{d\rho/dT_{1}}{4H\rho}\frac{1}{s}C_{D},$
(C.6)
$\displaystyle\frac{dY_{Z^{\prime}}}{dT_{1}}=-\frac{d\rho/dT_{1}}{4H\rho}\frac{1}{s}C_{Z^{\prime}},$
$\displaystyle\frac{dY_{\gamma^{\prime}}}{dT_{1}}=-\frac{d\rho/dT_{1}}{4H\rho}\frac{1}{s}C_{\gamma^{\prime}},$

where $C_{D}$, $C_{Z^{\prime}}$ and $C_{\gamma^{\prime}}$ are given by

$\displaystyle C_{D}=$ $\displaystyle n_{i}^{2}(T)\langle\sigma
v\rangle_{i\bar{i}\to D\bar{D}}(T)+n_{Z^{\prime}}^{2}\langle\sigma
v\rangle_{Z^{\prime}Z^{\prime}\to
D\bar{D}}(T_{1})-\frac{1}{2}n_{D}^{2}\langle\sigma v\rangle_{D\bar{D}\to
i\bar{i}}(T_{1})$ $\displaystyle-\frac{1}{2}n_{D}^{2}\langle\sigma
v\rangle_{D\bar{D}\to
Z^{\prime}Z^{\prime}}(T_{1})-\frac{1}{2}n_{D}^{2}\langle\sigma
v\rangle_{D\bar{D}\to\gamma^{\prime}\gamma^{\prime}}(T_{1})+n_{\gamma^{\prime}}^{2}\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to D\bar{D}}(T_{2})$
$\displaystyle-\frac{1}{2}n_{D}^{2}\langle\sigma v\rangle_{D\bar{D}\to
Z^{\prime}\gamma^{\prime}}(T_{1})+n_{Z^{\prime}}n_{\gamma^{\prime}}\langle\sigma
v\rangle_{Z^{\prime}\gamma^{\prime}\to D\bar{D}}(T_{1},T_{2}),$ (C.7)
$\displaystyle C_{Z^{\prime}}=$ $\displaystyle n_{i}^{2}(T)\langle\sigma
v\rangle_{i\bar{i}\to Z^{\prime}}(T)+n_{i}^{2}(T)\langle\sigma
v\rangle_{i\bar{i}\to
Z^{\prime}Z^{\prime}}(T)+\frac{1}{2}n_{D}^{2}\langle\sigma
v\rangle_{D\bar{D}\to Z^{\prime}Z^{\prime}}(T_{1})$ $\displaystyle-
n_{Z^{\prime}}^{2}\langle\sigma v\rangle_{Z^{\prime}Z^{\prime}\to
i\bar{i}}(T_{1})-n_{Z^{\prime}}^{2}\langle\sigma
v\rangle_{Z^{\prime}Z^{\prime}\to
D\bar{D}}(T_{1})-\frac{1}{2}n_{D}^{2}\langle\sigma v\rangle_{D\bar{D}\to
Z^{\prime}\gamma^{\prime}}(T_{1})$ $\displaystyle-
n_{Z^{\prime}}n_{\gamma^{\prime}}\langle\sigma
v\rangle_{Z^{\prime}\gamma^{\prime}\to
D\bar{D}}(T_{1},T_{2})-n_{Z^{\prime}}\langle\Gamma_{Z^{\prime}\to
i\bar{i}}\rangle(T_{1}),$ (C.8) $\displaystyle C_{\gamma^{\prime}}=$
$\displaystyle n_{i}^{2}(T)\langle\sigma
v\rangle_{i\bar{i}\to\gamma^{\prime}}(T)+\frac{1}{2}n_{D}^{2}\langle\sigma
v\rangle_{D\bar{D}\to\gamma^{\prime}\gamma^{\prime}}(T_{1})-n_{\gamma^{\prime}}\langle\Gamma_{\gamma^{\prime}\to
i\bar{i}}\rangle(T_{2})-n_{\gamma^{\prime}}^{2}\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to D\bar{D}}(T_{2})$
$\displaystyle+n_{i}^{2}(T)\langle\sigma
v\rangle_{i\bar{i}\to\gamma^{\prime}\gamma^{\prime}}(T)-n_{\gamma^{\prime}}^{2}\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to
i\bar{i}}(T_{2})+\frac{1}{2}n_{D}^{2}\langle\sigma v\rangle_{D\bar{D}\to
Z^{\prime}\gamma^{\prime}}(T_{1})$ $\displaystyle-
n_{Z^{\prime}}n_{\gamma^{\prime}}\langle\sigma
v\rangle_{Z^{\prime}\gamma^{\prime}\to D\bar{D}}(T_{1},T_{2}).$ (C.9)

In Eqs. (C.7), (C.8) and (C.9) one encounters thermally averaged decay width
and thermally averaged cross sections. The thermally averaged decay width is
given by

$\langle\Gamma_{a\to bc}\rangle=\Gamma_{a\to
bc}\frac{K_{1}(m_{a}/T)}{K_{2}(m_{a}/T)},$ (C.10)

and the thermally averaged cross-section is given by

$\langle\sigma v\rangle^{a\bar{a}\to
bc}(T)=\frac{1}{8m^{4}_{a}TK^{2}_{2}(m_{a}/T)}\int_{4m_{a}^{2}}^{\infty}ds~{}\sigma(s)\sqrt{s}\,(s-4m_{a}^{2})K_{1}(\sqrt{s}/T).$
(C.11)

$K_{1}$ and $K_{2}$ are the modified Bessel functions of the second kind and
degrees one and two, respectively. For the case when the annihilating
particles have different masses $m_{1}$ and $m_{2}$ and are at different
temperatures $T_{1}$ and $T_{2}$, the thermally averaged cross-section becomes

$\displaystyle\langle\sigma v\rangle_{12\to 34}(T_{1},T_{2})=$
$\displaystyle\frac{1}{4m_{1}^{2}m_{2}^{2}K_{2}(m_{1}/T_{1})K_{2}(m_{2}/T_{2})}\int_{(m_{1}+m_{2})^{2}}^{\infty}ds\,\sigma(s)\sqrt{s(s-(m_{1}+m_{2})^{2})}I(s),$
(C.12)

where

$I(s)=\frac{1}{T_{1}-T_{2}}\int_{\sqrt{s}}^{\infty}dx\,e^{-a_{+}x/2}\sinh\left(\frac{a_{-}}{2}\sqrt{1-\frac{(m_{1}+m_{2})^{2}}{s}}\sqrt{x^{2}-s}\right),$
(C.13)

and where

$a_{+}=\frac{T_{1}+T_{2}}{T_{1}T_{2}},~{}~{}~{}a_{-}=\frac{T_{1}-T_{2}}{T_{1}T_{2}}.$
(C.14)

Note that in the limit $T_{1}\to T_{2}$,
$I(s)\to\dfrac{\sqrt{s-(m_{1}+m_{2})^{2}}}{2T}K_{1}(\sqrt{s}/T)$ which for
$m_{1}=m_{2}$ allows us to recover Eq. (C.11) using Eq. (C.12). The
equilibrium yield of species $i$ is given by

$Y_{i}=\frac{n_{i}^{\rm
eq}}{s}=\frac{g_{i}}{2\pi^{2}s}m_{i}^{2}TK_{2}(m_{i}/T).$ (C.15)

The hidden sectors degrees of freedom is given by

$\displaystyle g_{1\rm eff}$ $\displaystyle=g^{Z^{\prime}}_{\rm
eff}+\frac{7}{8}g^{D}_{\rm eff},~{}~{}~{}\text{and}~{}~{}~{}h_{1\rm
eff}=h^{Z^{\prime}}_{\rm eff}+\frac{7}{8}h^{D}_{\rm eff},$ (C.16)
$\displaystyle g_{2\rm eff}$ $\displaystyle=g^{\gamma^{\prime}}_{\rm
eff},~{}~{}~{}\text{and}~{}~{}~{}h_{2\rm eff}=h^{\gamma^{\prime}}_{\rm eff},$
(C.17)

where

$\displaystyle g^{V}_{\rm eff}$
$\displaystyle=\frac{45}{\pi^{4}}\int_{x_{V}}^{\infty}\frac{\sqrt{x^{2}-x_{V}^{2}}}{e^{x}-1}x^{2}dx,~{}~{}~{}\text{and}~{}~{}~{}h^{V}_{\rm
eff}=\frac{45}{4\pi^{4}}\int_{x_{V}}^{\infty}\frac{\sqrt{x^{2}-x_{V}^{2}}}{e^{x}-1}(4x^{2}-x_{V}^{2})dx,$
(C.18) $\displaystyle g^{D}_{\rm eff}$
$\displaystyle=\frac{60}{\pi^{4}}\int_{x_{D}}^{\infty}\frac{\sqrt{x^{2}-x_{D}^{2}}}{e^{x}+1}x^{2}dx,~{}~{}~{}\text{and}~{}~{}~{}h^{D}_{\rm
eff}=\frac{15}{\pi^{4}}\int_{x_{D}}^{\infty}\frac{\sqrt{x^{2}-x_{D}^{2}}}{e^{x}+1}(4x^{2}-x_{D}^{2})dx.$

In Eq. (C.18), $V=Z^{\prime},\gamma^{\prime}$ and we take
$g_{\gamma^{\prime}}=g_{Z^{\prime}}=3$ and $g_{D}=4$.

##  Appendix D Dark photon and dark fermion scattering cross sections and
$Z^{\prime}$ decay width

The calculation of the relic densities of the dark photon and dark fermion
require solving the coupled Boltzmann equations which contain a variety of
cross sections involving the standard model and dark sector particles. We list
these below.

  1. 1.

Processes: $D\bar{D}\to Z,Z^{\prime},\gamma^{\prime}\to f\bar{f}$

$\displaystyle\sigma^{D\bar{D}\to f\bar{f}}(s)$
$\displaystyle=\frac{g_{X}^{2}g_{2}^{2}N_{c}}{12\pi\cos^{2}\theta_{w}}\left(1+\frac{2m_{D}^{2}}{s}\right)\sqrt{\frac{s-4m_{f}^{2}}{s-4m_{D}^{2}}}\Bigg{\\{}\frac{[a_{f}^{\prime\prime
2}(s-4m^{2}_{f})+v_{f}^{\prime\prime
2}(s+2m^{2}_{f})]\delta_{2}^{2}}{\kappa^{2}[(s-m^{2}_{\gamma^{\prime}})^{2}+m^{2}_{\gamma^{\prime}}\Gamma^{2}_{\gamma^{\prime}}]}$
$\displaystyle+\frac{a_{f}^{\prime 2}(s-4m^{2}_{f})+v_{f}^{\prime
2}(s+2m^{2}_{f})}{(s-m^{2}_{Z^{\prime}})^{2}+m^{2}_{Z^{\prime}}\Gamma^{2}_{Z^{\prime}}}+\frac{[a_{f}^{2}(s-4m^{2}_{f})+v_{f}^{2}(s+2m^{2}_{f})]\delta_{1}^{2}\sin^{2}\theta_{w}}{(s-m^{2}_{Z})^{2}+m^{2}_{Z}\Gamma^{2}_{Z}}$
$\displaystyle+\frac{2\delta_{2}[a_{f}^{\prime}a_{f}^{\prime\prime}(s-4m^{2}_{f})+v_{f}^{\prime}v_{f}^{\prime\prime}(s+2m^{2}_{f})]}{\kappa[(s-m_{Z^{\prime}}^{2})^{2}+m^{2}_{Z^{\prime}}\Gamma^{2}_{Z^{\prime}}][(s-m_{\gamma^{\prime}}^{2})^{2}+m^{2}_{\gamma^{\prime}}\Gamma^{2}_{\gamma^{\prime}}]}G(s,m_{Z^{\prime}},m_{\gamma^{\prime}})$
$\displaystyle+\frac{2\delta_{1}\delta_{2}[a_{f}a_{f}^{\prime\prime}(s-4m^{2}_{f})+v_{f}v_{f}^{\prime\prime}(s+2m^{2}_{f})]\sin\theta_{w}}{\kappa[(s-m_{Z}^{2})^{2}+m^{2}_{Z}\Gamma^{2}_{Z}][(s-m_{\gamma^{\prime}}^{2})^{2}+m^{2}_{\gamma^{\prime}}\Gamma^{2}_{\gamma^{\prime}}]}G(s,m_{Z},m_{\gamma^{\prime}})$
$\displaystyle+\frac{2\delta_{1}[a_{f}a_{f}^{\prime}(s-4m^{2}_{f})+v_{f}v_{f}^{\prime}(s+2m^{2}_{f})]\sin\theta_{w}}{[(s-m_{Z^{\prime}}^{2})^{2}+m^{2}_{Z^{\prime}}\Gamma^{2}_{Z^{\prime}}][(s-m_{Z}^{2})^{2}+m^{2}_{Z}\Gamma^{2}_{Z}]}G(s,m_{Z},m_{Z^{\prime}})\Bigg{\\}},$
(D.1)

Here $s$ is the Mandelstam variable which gives the square of the total energy
in the CM system. Further, the notation used above is as follows:
$f=e,\mu,\tau$: $T_{3f}=-1/2$ and $Q_{f}=-1$ and for $f=u,c,t$:
$T_{f}^{3}=1/2$ and $Q_{f}=2/3$ and for $f=d,s,b$: $T_{f}^{3}=-1/2$ and
$Q_{f}=-1/3$, $N_{c}$ is the color number and

$\displaystyle\kappa$
$\displaystyle=-(1-m^{2}_{Z^{\prime}}/m^{2}_{\gamma^{\prime}}),$ (D.2)
$\displaystyle G(s,m_{1},m_{2})$
$\displaystyle=(s-m^{2}_{1})(s-m^{2}_{2})+\Gamma_{1}\Gamma_{2}m_{1}m_{2}.$
(D.3)

  2. 2.

Processes: $D\bar{D}\to Z,Z^{\prime},\gamma^{\prime}\to\nu\bar{\nu}$

$\displaystyle\sigma^{D\bar{D}\to\nu\bar{\nu}}(s)=$
$\displaystyle\frac{g_{X}^{2}g_{2}^{2}\delta_{1}^{2}}{8\pi}\frac{(s+4m_{D}^{2})\tan^{2}\theta_{w}}{(1-4m^{2}_{D}/s)^{1/2}}\left[\frac{A}{(s-m^{2}_{Z^{\prime}})(s-m^{2}_{\gamma^{\prime}})}-\frac{1}{(s-m^{2}_{Z})}\right]^{2},$
(D.4)

where

$A=\epsilon^{2}_{z}(s-m^{2}_{\gamma^{\prime}})+\frac{\delta_{2}\epsilon^{2}_{\gamma^{\prime}}(s-m^{2}_{Z^{\prime}})(\delta_{2}-\sin\beta)}{1-m^{2}_{Z^{\prime}}/m^{2}_{\gamma^{\prime}}}.$
(D.5)

  3. 3.

Processes: $f\bar{f}\to Z,Z^{\prime},\gamma^{\prime}\to D\bar{D}$

$(s-4m^{2}_{D})\sigma^{D\bar{D}\to
f\bar{f}}(s)=N_{c}^{2}(s-4m^{2}_{f})\sigma^{f\bar{f}\to D\bar{D}}(s).$ (D.6)

  4. 4.

Processes: $\nu\bar{\nu}\to Z,Z^{\prime},\gamma^{\prime}\to D\bar{D}$

$\sigma^{\nu\bar{\nu}\to
D\bar{D}}(s)=4\left(1-\frac{4m^{2}_{D}}{s}\right)\sigma^{D\bar{D}\to\nu\bar{\nu}}(s).$
(D.7)

  5. 5.

Process: $D\bar{D}\to Z^{\prime}Z^{\prime}$

$\displaystyle\sigma^{D\bar{D}\to Z^{\prime}Z^{\prime}}(s)$
$\displaystyle=\frac{g_{X}^{4}}{8\pi
s}\Bigg{\\{}-\sqrt{\frac{s-4m^{2}_{Z^{\prime}}}{s-4m^{2}_{D}}}\left[\frac{m^{2}_{D}s+2m^{4}_{Z^{\prime}}+4m^{4}_{D}}{(s-4m^{2}_{Z^{\prime}})m^{2}_{D}+m^{4}_{Z^{\prime}}}\right]$
$\displaystyle+\frac{s^{2}+4m^{2}_{D}(s-2m^{2}_{Z^{\prime}})+4m^{4}_{Z^{\prime}}-8m^{4}_{D}}{(s-2m^{2}_{Z^{\prime}})(s-4m^{2}_{D})}\log
B\Bigg{\\}},$ (D.8)

where

$B=\frac{s-2m^{2}_{Z^{\prime}}+\sqrt{(s-4m^{2}_{Z^{\prime}})(s-4m^{2}_{D})}}{s-2m^{2}_{Z^{\prime}}-\sqrt{(s-4m^{2}_{Z^{\prime}})(s-4m^{2}_{D})}}.$
(D.9)

  6. 6.

Process: $D\bar{D}\to\gamma^{\prime}\gamma^{\prime}$

$\sigma^{D\bar{D}\to\gamma^{\prime}\gamma^{\prime}}(s)=\left(\frac{\delta_{2}M^{2}_{\gamma^{\prime}}}{M^{2}_{Z^{\prime}}-M^{2}_{\gamma^{\prime}}}\right)^{4}\sigma^{D\bar{D}\to
Z^{\prime}Z^{\prime}}(s)\Big{|}_{m_{Z^{\prime}}\longleftrightarrow
m_{\gamma^{\prime}}}.$ (D.10)

  7. 7.

Processes: $D\bar{D}\to VV$ with $V=Z^{\prime},\gamma^{\prime}$

$8(s-4m^{2}_{D})\sigma^{D\bar{D}\to VV}(s)=9(s-4m^{2}_{V})\sigma^{VV\to
D\bar{D}}(s).$

  8. 8.

Processes $D\bar{D}\rightarrow Z^{\prime}\gamma^{\prime}$

$\displaystyle\sigma^{D\bar{D}\rightarrow Z^{\prime}\gamma^{\prime}}(s)$
$\displaystyle=\frac{\delta_{2}^{2}g_{X}^{4}m_{\gamma^{\prime}}^{4}}{4\pi(m_{Z^{\prime}}^{2}-m_{\gamma^{\prime}}^{2})s(s-4m_{D}^{2})}$
$\displaystyle\times\bigg{\\{}\frac{4m_{D}^{4}s-2m_{Z^{\prime}}m_{\gamma^{\prime}}s+m_{D}^{2}\left[(m_{Z^{\prime}}^{2}-m_{\gamma^{\prime}}^{2})^{2}+s^{2}\right]}{m_{Z^{\prime}}^{2}m_{\gamma^{\prime}}^{2}s+m_{D}^{2}\left[m_{\gamma^{\prime}}^{4}+(s-m_{Z^{\prime}}^{2})^{2}-2m_{\gamma^{\prime}}^{2}(s+m_{Z^{\prime}}^{2})\right]}E$
$\displaystyle+\left[-8m_{D}^{4}+(m_{Z^{\prime}}^{2}+m_{\gamma^{\prime}}^{2})^{2}-4m_{D}^{2}(m_{Z^{\prime}}^{2}-m_{\gamma^{\prime}}^{2}-s)+s^{2}\right]\log
F\bigg{\\}},$ (D.11)

where

$\displaystyle E$
$\displaystyle=\sqrt{\left(1-\frac{4m_{D}^{2}}{s}\right)\left[m_{\gamma^{\prime}}^{4}+(s-m_{Z^{\prime}}^{2})^{2}-2m_{\gamma^{\prime}}^{2}(s+m_{Z^{\prime}}^{2})\right]},$
(D.12) $\displaystyle F$
$\displaystyle=\frac{m_{\gamma^{\prime}}^{2}+m_{Z^{\prime}}^{2}-s+E}{m_{\gamma^{\prime}}^{2}+m_{Z^{\prime}}^{2}-s-E}.$
(D.13)

  9. 9.

Processes $Z^{\prime}\gamma^{\prime}\rightarrow D\bar{D}$

$\sigma^{Z^{\prime}\gamma^{\prime}\rightarrow
D\bar{D}}(s)=\frac{4s(s-4m_{D}^{2})}{9[(m_{\gamma^{\prime}}^{2}+m_{Z^{\prime}}^{2}-s)^{2}]-4m_{\gamma^{\prime}}^{2}m_{Z^{\prime}}^{2}}\sigma^{D\bar{D}\rightarrow
Z^{\prime}\gamma^{\prime}}(s)$ (D.14)

  10. 10.

Processes: $VV\to f\bar{f}$ with $V=Z^{\prime},\gamma^{\prime}$

$\displaystyle\sigma^{VV\to f\bar{f}}(s)$
$\displaystyle=\frac{g_{2}^{4}N_{c}}{9\pi
m^{4}_{V}s(s-4m^{2}_{V})\cos^{4}\theta_{w}}\Bigg{\\{}\frac{\sqrt{(s-4m^{2}_{f})(s-4m^{2}_{V})}}{m^{4}_{V}+m^{2}_{f}(s-4m^{2}_{V})}$
$\displaystyle\times\Big{(}c_{A}^{4}[-2m_{V}^{8}+m_{f}^{2}m_{V}^{4}(s+4m_{V}^{2})+2m_{f}^{4}(8m_{V}^{4}-8m_{V}^{2}s+s^{2})]$
$\displaystyle~{}~{}~{}~{}+2c_{A}^{2}c_{V}^{2}m_{V}^{4}[8m_{f}^{4}-6m_{V}^{4}+m_{f}^{2}(22m_{V}^{2}-7s)]$
$\displaystyle~{}~{}~{}~{}-m_{V}^{4}(4m_{f}^{4}+2m_{V}^{4}+m_{f}^{2}s)c_{V}^{4}\Big{)}$
$\displaystyle+\frac{\log
C}{(s-2m^{2}_{V})}\Big{(}c_{A}^{4}[4m_{f}^{4}(4m_{V}^{2}-s)s+m_{V}^{4}(4m_{V}^{4}+s^{2})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+4m_{f}^{2}m_{V}^{2}(-4m_{V}^{4}-3m_{V}^{2}s+s^{2})]+2c_{A}^{2}m_{V}^{2}[16m_{f}^{4}m_{V}^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+3m_{V}^{2}(s^{2}+4m_{V}^{4})+2m_{f}^{2}(s^{2}-10m_{V}^{2}s-10m_{V}^{4})]c_{V}^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+m_{V}^{4}[s^{2}+4m_{f}^{2}(s-2m_{V}^{2})-8m_{f}^{4}+4m_{V}^{4}]c_{V}^{4}\Big{)}\Bigg{\\}},$
(D.15)

where

$C=\frac{s-2m^{2}_{V}+\sqrt{(s-4m^{2}_{V})(s-4m^{2}_{f})}}{s-2m^{2}_{V}-\sqrt{(s-4m^{2}_{V})(s-4m^{2}_{f})}},$
(D.16)

and $c_{A}=a_{f}^{\prime}$, $c_{V}=v_{f}^{\prime}$ for $V=Z^{\prime}$ and
$c_{A}=a_{f}^{\prime\prime}$, $c_{V}=v_{f}^{\prime\prime}$ for
$V=\gamma^{\prime}$.

  11. 11.

Process: $VV\to\nu\bar{\nu}$ with $V=Z^{\prime},\gamma^{\prime}$

$\displaystyle\sigma^{VV\to\nu\bar{\nu}}(s)$
$\displaystyle=\frac{g_{2}^{4}\delta_{1}^{4}}{6\pi
s}\tan^{4}\theta_{w}\Bigg{[}\frac{s^{2}+4m^{4}_{V}}{(s-4m^{2}_{V})(s-2m^{2}_{V})}\log
D-2\left(1-\frac{4m^{2}_{V}}{s}\right)^{-1/2}\Bigg{]}$ $\displaystyle\hskip
85.35826pt\times\begin{cases}\epsilon_{z}^{8},&\text{for }V=Z^{\prime}\\\
(\delta_{2}-\sin\beta)^{4}\epsilon_{\gamma^{\prime}}^{8},&\text{for
}V=\gamma^{\prime}\end{cases},$ (D.17)

where

$D=\frac{s-2m^{2}_{V}+\sqrt{s(s-4m^{2}_{V})}}{s-2m^{2}_{V}-\sqrt{s(s-4m^{2}_{V})}}.$
(D.18)

  12. 12.

Processes $f\bar{f},\nu,\bar{\nu}\to VV$ with $V=Z^{\prime},\gamma^{\prime}$

$\displaystyle 9(s-4m^{2}_{V})\sigma^{VV\to f\bar{f}}$
$\displaystyle=8(s-4m^{2}_{f})\sigma^{f\bar{f}\to VV},$ (D.19) $\displaystyle
9(s-4m^{2}_{V})\sigma^{VV\to\nu\bar{\nu}}$
$\displaystyle=2s\,\sigma^{\nu\bar{\nu}\to VV}.$ (D.20)

  13. 13.

Process: $f\bar{f}\to Z^{\prime}$

$\displaystyle\sigma^{f\bar{f}\to Z^{\prime}}(s)$ $\displaystyle=\frac{\pi
g_{2}^{2}\delta_{1}^{2}m_{Z^{\prime}}^{2}}{2s\sqrt{s-4m^{2}_{f}}N_{c}}\Bigg{[}\left(1+\frac{2m^{2}_{f}}{m_{Z^{\prime}}^{2}}\right)\left(Q_{f}+\frac{m_{Z^{\prime}}^{2}}{2m_{W}^{2}}(T_{3f}-2Q_{f})\right)^{2}\sin^{2}2\theta_{w}$
$\displaystyle\hskip
99.58464pt+\left(1-\frac{4m^{2}_{f}}{m^{2}_{Z^{\prime}}}\right)\epsilon_{z}^{4}T^{2}_{3f}\tan^{2}\theta_{w}\Bigg{]}\delta(\sqrt{s}-m_{Z^{\prime}}).$
(D.21)

  14. 14.

Process: $\nu\bar{\nu}\to Z^{\prime}$

$\displaystyle\sigma^{\nu\bar{\nu}\to Z^{\prime}}(s)$
$\displaystyle=\frac{3\pi
g_{2}^{2}\delta_{1}^{2}\epsilon_{z}^{4}m^{2}_{Z^{\prime}}}{s^{3/2}}\tan^{2}\theta_{w}\,\delta(\sqrt{s}-m_{Z^{\prime}}).$
(D.22)

  15. 15.

Process: $f\bar{f},\nu\bar{\nu}\to\gamma^{\prime}$

$\displaystyle\sigma^{f\bar{f}\to\gamma^{\prime}}(s)=(\delta_{2}-\sin\beta)^{2}\,\sigma^{f\bar{f}\to
Z^{\prime}}(s)\Big{|}_{m_{Z^{\prime}}\longleftrightarrow
m_{\gamma^{\prime}}}.$ (D.23)

Same applies to $\sigma^{\nu\bar{\nu}\to\gamma^{\prime}}(s)$.

  16. 16.

Process: $Z^{\prime}\to f\bar{f},\nu\bar{\nu}$

The decay width of $Z^{\prime}$ to SM fermions is given by

$\displaystyle\Gamma_{Z^{\prime}\to f\bar{f}}$
$\displaystyle=\frac{g_{2}^{2}\delta_{1}^{2}N_{c}}{12\pi}m_{Z^{\prime}}\sqrt{1-\left(\frac{2m_{f}}{m_{Z^{\prime}}}\right)^{2}}\Bigg{[}\left(1-\frac{4m^{2}_{f}}{m^{2}_{Z^{\prime}}}\right)T^{2}_{3f}\epsilon^{4}_{z}\tan^{2}\theta_{w}$
$\displaystyle+\left(1+\frac{2m^{2}_{f}}{m^{2}_{Z^{\prime}}}\right)\left(Q_{f}+\frac{m_{Z^{\prime}}^{2}}{2m_{W}^{2}}(T_{3f}-2Q_{f})\right)^{2}\sin^{2}2\theta_{w}\Bigg{]},$
(D.24)

and its invisible decay is

$\displaystyle\Gamma_{Z^{\prime}\to\nu\bar{\nu}}$
$\displaystyle=\frac{g_{2}^{2}\delta_{1}^{2}\epsilon_{z}^{4}}{8\pi}m_{Z^{\prime}}\tan^{2}\theta_{w}.$
(D.25)

## References

  * [1] M. R. Buckley and P. J. Fox, Phys. Rev. D 81, 083522 (2010) doi:10.1103/PhysRevD.81.083522 [arXiv:0911.3898 [hep-ph]].
  * [2] A. Loeb and N. Weiner, Phys. Rev. Lett. 106, 171302 (2011) doi:10.1103/PhysRevLett.106.171302 [arXiv:1011.6374 [astro-ph.CO]].
  * [3] M. Kaplinghat, S. Tulin and H. B. Yu, Phys. Rev. Lett. 116, no.4, 041302 (2016) doi:10.1103/PhysRevLett.116.041302 [arXiv:1508.03339 [astro-ph.CO]].
  * [4] L. Sagunski, S. Gad-Nasr, B. Colquhoun, A. Robertson and S. Tulin, JCAP 01, 024 (2021) doi:10.1088/1475-7516/2021/01/024 [arXiv:2006.12515 [astro-ph.CO]].
  * [5] A. Aboubrahim, W. Z. Feng, P. Nath and Z. Y. Wang, Phys. Rev. D 103, no.7, 075014 (2021) doi:10.1103/PhysRevD.103.075014 [arXiv:2008.00529 [hep-ph]].
  * [6] K. Kaneta, H. S. Lee and S. Yun, Phys. Rev. Lett. 118, no.10, 101802 (2017) doi:10.1103/PhysRevLett.118.101802 [arXiv:1611.01466 [hep-ph]].
  * [7] K. Kaneta, H. S. Lee and S. Yun, Phys. Rev. D 95, no.11, 115032 (2017) doi:10.1103/PhysRevD.95.115032 [arXiv:1704.07542 [hep-ph]].
  * [8] R. T. Co, A. Pierce, Z. Zhang and Y. Zhao, Phys. Rev. D 99, no.7, 075002 (2019) doi:10.1103/PhysRevD.99.075002 [arXiv:1810.07196 [hep-ph]].
  * [9] J. A. Dror, K. Harigaya and V. Narayan, Phys. Rev. D 99, no.3, 035036 (2019) doi:10.1103/PhysRevD.99.035036 [arXiv:1810.07195 [hep-ph]].
  * [10] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi and F. Takahashi, Phys. Lett. B 801, 135136 (2020) doi:10.1016/j.physletb.2019.135136 [arXiv:1810.07188 [hep-ph]].
  * [11] A. J. Long and L. T. Wang, Phys. Rev. D 99, no.6, 063529 (2019) doi:10.1103/PhysRevD.99.063529 [arXiv:1901.03312 [hep-ph]].
  * [12] G. Alonso-Álvarez, T. Hugle and J. Jaeckel, JCAP 02, 014 (2020) doi:10.1088/1475-7516/2020/02/014 [arXiv:1905.09836 [hep-ph]].
  * [13] Y. Nakai, R. Namba and Z. Wang, JHEP 12, 170 (2020) doi:10.1007/JHEP12(2020)170 [arXiv:2004.10743 [hep-ph]].
  * [14] G. Choi, T. T. Yanagida and N. Yokozaki, JHEP 01, 057 (2021) doi:10.1007/JHEP01(2021)057 [arXiv:2008.12180 [hep-ph]].
  * [15] C. Delaunay, T. Ma and Y. Soreq, JHEP 02, 010 (2021) doi:10.1007/JHEP02(2021)010 [arXiv:2009.03060 [hep-ph]].
  * [16] P. W. Graham, J. Mardon and S. Rajendran, Phys. Rev. D 93, no.10, 103520 (2016) doi:10.1103/PhysRevD.93.103520 [arXiv:1504.02102 [hep-ph]].
  * [17] Y. Ema, K. Nakayama and Y. Tang, JHEP 07, 060 (2019) doi:10.1007/JHEP07(2019)060 [arXiv:1903.10973 [hep-ph]].
  * [18] A. Ahmed, B. Grzadkowski and A. Socha, JHEP 08, 059 (2020) doi:10.1007/JHEP08(2020)059 [arXiv:2005.01766 [hep-ph]].
  * [19] S. Tulin and H. B. Yu, Phys. Rept. 730, 1-57 (2018) doi:10.1016/j.physrep.2017.11.004 [arXiv:1705.02358 [hep-ph]].
  * [20] J. Alexander, M. Battaglieri, B. Echenard, R. Essig, M. Graham, E. Izaguirre, J. Jaros, G. Krnjaic, J. Mardon and D. Morrissey, et al. [arXiv:1608.08632 [hep-ph]].
  * [21] M. Fabbrichesi, E. Gabrielli and G. Lanfranchi, doi:10.1007/978-3-030-62519-1 [arXiv:2005.01515 [hep-ph]].
  * [22] I. M. Bloch, R. Essig, K. Tobioka, T. Volansky and T. T. Yu, JHEP 06, 087 (2017) doi:10.1007/JHEP06(2017)087 [arXiv:1608.02123 [hep-ph]].
  * [23] M. Pospelov, A. Ritz and M. B. Voloshin, Phys. Rev. D 78, 115012 (2008) doi:10.1103/PhysRevD.78.115012 [arXiv:0807.3279 [hep-ph]].
  * [24] J. E. Kim and D. J. E. Marsh, Phys. Rev. D 93, no.2, 025027 (2016) doi:10.1103/PhysRevD.93.025027 [arXiv:1510.01701 [hep-ph]].
  * [25] L. Hui, J. P. Ostriker, S. Tremaine and E. Witten, Phys. Rev. D 95, no.4, 043541 (2017) doi:10.1103/PhysRevD.95.043541 [arXiv:1610.08297 [astro-ph.CO]].
  * [26] J. Halverson, C. Long and P. Nath, Phys. Rev. D 96, no.5, 056025 (2017) doi:10.1103/PhysRevD.96.056025 [arXiv:1703.07779 [hep-ph]].
  * [27] N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].
  * [28] B. Holdom, Phys. Lett. B 166, 196-198 (1986) doi:10.1016/0370-2693(86)91377-8
  * [29] M. Dutra, M. Lindner, S. Profumo, F. S. Queiroz, W. Rodejohann and C. Siqueira, JCAP 03, 037 (2018) doi:10.1088/1475-7516/2018/03/037 [arXiv:1801.05447 [hep-ph]].
  * [30] B. Kors and P. Nath, Phys. Lett. B 586, 366-372 (2004) doi:10.1016/j.physletb.2004.02.051 [arXiv:hep-ph/0402047 [hep-ph]].
  * [31] K. Cheung and T. C. Yuan, JHEP 03, 120 (2007) doi:10.1088/1126-6708/2007/03/120 [arXiv:hep-ph/0701107 [hep-ph]].
  * [32] D. Feldman, Z. Liu and P. Nath, JHEP 11, 007 (2006) doi:10.1088/1126-6708/2006/11/007 [arXiv:hep-ph/0606294 [hep-ph]].
  * [33] D. Feldman, Z. Liu and P. Nath, Phys. Rev. D 75, 115001 (2007) doi:10.1103/PhysRevD.75.115001 [arXiv:hep-ph/0702123 [hep-ph]].
  * [34] A. Aboubrahim and P. Nath, Phys. Rev. D 99, no.5, 055037 (2019) doi:10.1103/PhysRevD.99.055037 [arXiv:1902.05538 [hep-ph]].
  * [35] L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 03, 080 (2010) doi:10.1007/JHEP03(2010)080 [arXiv:0911.1120 [hep-ph]].
  * [36] A. Aboubrahim, W. Z. Feng and P. Nath, JHEP 02, 118 (2020) doi:10.1007/JHEP02(2020)118 [arXiv:1910.14092 [hep-ph]].
  * [37] A. Aboubrahim, W. Z. Feng and P. Nath, JHEP 04, 144 (2020) doi:10.1007/JHEP04(2020)144 [arXiv:2003.02267 [hep-ph]].
  * [38] S. Koren and R. McGehee, Phys. Rev. D 101, no.5, 055024 (2020) doi:10.1103/PhysRevD.101.055024 [arXiv:1908.03559 [hep-ph]].
  * [39] Y. Du, F. Huang, H. L. Li and J. H. Yu, JHEP 12, 207 (2020) doi:10.1007/JHEP12(2020)207 [arXiv:2005.01717 [hep-ph]].
  * [40] J. L. Feng, H. Tu and H. B. Yu, JCAP 10, 043 (2008) doi:10.1088/1475-7516/2008/10/043 [arXiv:0808.2318 [hep-ph]].
  * [41] X. Chu, T. Hambye and M. H. G. Tytgat, JCAP 05, 034 (2012) doi:10.1088/1475-7516/2012/05/034 [arXiv:1112.0493 [hep-ph]].
  * [42] L. Ackerman, M. R. Buckley, S. M. Carroll and M. Kamionkowski, doi:10.1103/PhysRevD.79.023519 [arXiv:0810.5126 [hep-ph]].
  * [43] R. Foot and S. Vagnozzi, Phys. Rev. D 91, 023512 (2015) doi:10.1103/PhysRevD.91.023512 [arXiv:1409.7174 [hep-ph]].
  * [44] R. Foot and S. Vagnozzi, JCAP 07, 013 (2016) doi:10.1088/1475-7516/2016/07/013 [arXiv:1602.02467 [astro-ph.CO]].
  * [45] T. Hambye, M. H. G. Tytgat, J. Vandecasteele and L. Vanderheyden, Phys. Rev. D 100, no.9, 095018 (2019) doi:10.1103/PhysRevD.100.095018 [arXiv:1908.09864 [hep-ph]].
  * [46] P. Gondolo and G. Gelmini, Nucl. Phys. B 360, 145-179 (1991) doi:10.1016/0550-3213(91)90438-4
  * [47] G. B. Gelmini, P. Gondolo and E. Roulet, Nucl. Phys. B 351, 623-644 (1991) doi:10.1016/S0550-3213(05)80036-7
  * [48] S. D. McDermott, H. H. Patel and H. Ramani, Phys. Rev. D 97, no.7, 073005 (2018) doi:10.1103/PhysRevD.97.073005 [arXiv:1705.00619 [hep-ph]].
  * [49] T. Bringmann, New J. Phys. 11, 105027 (2009) doi:10.1088/1367-2630/11/10/105027 [arXiv:0903.0189 [astro-ph.CO]].
  * [50] A. Bhoonah, J. Bramante, F. Elahi and S. Schon, Phys. Rev. D 100, no.2, 023001 (2019) doi:10.1103/PhysRevD.100.023001 [arXiv:1812.10919 [hep-ph]].
  * [51] J. Redondo and M. Postma, JCAP 02, 005 (2009) doi:10.1088/1475-7516/2009/02/005 [arXiv:0811.0326 [hep-ph]].
  * [52] A. Caputo, H. Liu, S. Mishra-Sharma and J. T. Ruderman, Phys. Rev. Lett. 125, no.22, 221303 (2020) doi:10.1103/PhysRevLett.125.221303 [arXiv:2002.05165 [astro-ph.CO]].
  * [53] A. A. Garcia, K. Bondarenko, S. Ploeckinger, J. Pradler and A. Sokolenko, JCAP 10, 011 (2020) doi:10.1088/1475-7516/2020/10/011 [arXiv:2003.10465 [astro-ph.CO]].
  * [54] S. J. Witte, S. Rosauro-Alcaraz, S. D. McDermott and V. Poulin, JHEP 06, 132 (2020) doi:10.1007/JHEP06(2020)132 [arXiv:2003.13698 [astro-ph.CO]].
  * [55] S. D. McDermott and S. J. Witte, Phys. Rev. D 101, no.6, 063030 (2020) doi:10.1103/PhysRevD.101.063030 [arXiv:1911.05086 [hep-ph]].
  * [56] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, JCAP 06, 013 (2012) doi:10.1088/1475-7516/2012/06/013 [arXiv:1201.5902 [hep-ph]].
  * [57] M. G. Aartsen et al. [IceCube], Eur. Phys. J. C 78, no.10, 831 (2018) doi:10.1140/epjc/s10052-018-6273-3 [arXiv:1804.03848 [astro-ph.HE]].
  * [58] P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016) doi:10.1051/0004-6361/201525830 [arXiv:1502.01589 [astro-ph.CO]].
  * [59] T. R. Slatyer, Phys. Rev. D 93, no.2, 023527 (2016) doi:10.1103/PhysRevD.93.023527 [arXiv:1506.03811 [hep-ph]].
  * [60] M. Endo, K. Hamaguchi and G. Mishima, Phys. Rev. D 86, 095029 (2012) doi:10.1103/PhysRevD.86.095029 [arXiv:1209.2558 [hep-ph]].
  * [61] J. P. Lees et al. [BaBar], Phys. Rev. Lett. 113, no.20, 201801 (2014) doi:10.1103/PhysRevLett.113.201801 [arXiv:1406.2980 [hep-ex]].
  * [62] F. Bergsma et al. [CHARM], Phys. Lett. B 157, 458-462 (1985) doi:10.1016/0370-2693(85)90400-9
  * [63] Y. D. Tsai, P. deNiverville and M. X. Liu, Phys. Rev. Lett. 126, no.18, 181801 (2021) doi:10.1103/PhysRevLett.126.181801 [arXiv:1908.07525 [hep-ph]].
  * [64] J. R. Batley et al. [NA48/2], Phys. Lett. B 746, 178-185 (2015) doi:10.1016/j.physletb.2015.04.068 [arXiv:1504.00607 [hep-ex]].
  * [65] S. Andreas, C. Niebuhr and A. Ringwald, Phys. Rev. D 86, 095019 (2012) doi:10.1103/PhysRevD.86.095019 [arXiv:1209.6083 [hep-ph]].
  * [66] J. D. Bjorken, R. Essig, P. Schuster and N. Toro, Phys. Rev. D 80, 075018 (2009) doi:10.1103/PhysRevD.80.075018 [arXiv:0906.0580 [hep-ph]].
  * [67] D. Banerjee et al. [NA64], Phys. Rev. Lett. 120, no.23, 231802 (2018) doi:10.1103/PhysRevLett.120.231802 [arXiv:1803.07748 [hep-ex]].
  * [68] D. Banerjee et al. [NA64], Phys. Rev. D 101, no.7, 071101 (2020) doi:10.1103/PhysRevD.101.071101 [arXiv:1912.11389 [hep-ex]].
  * [69] E. M. Riordan, M. W. Krasny, K. Lang, P. De Barbaro, A. Bodek, S. Dasu, N. Varelas, X. Wang, R. G. Arnold and D. Benton, et al. Phys. Rev. Lett. 59, 755 (1987) doi:10.1103/PhysRevLett.59.755
  * [70] J. Blumlein, J. Brunner, H. J. Grabosch, P. Lanius, S. Nowak, C. Rethfeldt, H. E. Ryseck, M. Walter, D. Kiss and Z. Jaki, et al. Z. Phys. C 51, 341-350 (1991) doi:10.1007/BF01548556
  * [71] J. Blumlein, J. Brunner, H. J. Grabosch, P. Lanius, S. Nowak, C. Rethfeldt, H. E. Ryseck, M. Walter, D. Kiss and Z. Jaki, et al. Int. J. Mod. Phys. A 7, 3835-3850 (1992) doi:10.1142/S0217751X9200171X
  * [72] P. Ilten, Y. Soreq, M. Williams and W. Xue, JHEP 06, 004 (2018) doi:10.1007/JHEP06(2018)004 [arXiv:1801.04847 [hep-ph]].
  * [73] J. H. Chang, R. Essig and S. D. McDermott, JHEP 01, 107 (2017) doi:10.1007/JHEP01(2017)107 [arXiv:1611.03864 [hep-ph]].
  * [74] H. An, M. Pospelov and J. Pradler, Phys. Lett. B 725, 190-195 (2013) doi:10.1016/j.physletb.2013.07.008 [arXiv:1302.3884 [hep-ph]].
  * [75] R. Essig, E. Kuflik, S. D. McDermott, T. Volansky and K. M. Zurek, JHEP 11, 193 (2013) doi:10.1007/JHEP11(2013)193 [arXiv:1309.4091 [hep-ph]].