# Sub-GeV millicharge dark matter from the $U(1)_{X}$ hidden sector

Wan-Zhe Feng111Email<EMAIL_ADDRESS>Zi-Hui Zhang222Email:
<EMAIL_ADDRESS>Kai-Yu Zhang333Email<EMAIL_ADDRESS>
Center for Joint Quantum Studies and Department of Physics,
School of Science, Tianjin University, Tianjin 300350, PR. China

###### Abstract

We conduct a comprehensive study on the sub-GeV millicharge dark matter
produced through the freeze-in mechanism. We discuss in general the mixing
mechanism, encompassing both kinetic mixing and mass mixing, between the
$U(1)_{X}$ hidden sector and the standard model, which can generate
millicharge carried by the dark fermions from the hidden sector. We discuss in
depth how such millicharge is generated, and clarify several misunderstandings
regarding this subject in the literature. Without employing an effective field
theory approach, where the photon field directly mixed with the additional
$U(1)$, we analyze a general renormalizable model and investigate the complete
evolution of the hidden sector particles precisely. Due to the substantial
self-interactions among hidden sector particles, the evolution of the hidden
sector temperature plays a crucial role, which is addressed concurrently with
the number densities of hidden sector particles by solving a set of coupled
Boltzmann equations. We thoroughly examine eight benchmark models covering six
distinct cases. Some of our key findings from the analysis of these benchmark
models may be generalizable and applicable to broader freeze-in scenarios. We
also explore the possibility that the $\mathcal{O}$(keV) $U(1)_{X}$ dark
photon is a viable dark matter candidate, even though it can contribute at
most $\sim 5\%$ to the total observed dark matter relic density.

###### Contents

  1. 1 Introduction
  2. 2 The $U(1)$ extension of the SM and the generation of millicharge particle
    1. 2.1 Kinetic mixing between two massless $U(1)$’s
    2. 2.2 Kinetic mixing between a massless $U(1)$ and a massive $U(1)$
    3. 2.3 Kinetic mixing and mass mixing between two massive $U(1)$’s
    4. 2.4 Kinetic mixing and mass mixing in the full electroweak theory
      1. 2.4.1 $U(1)_{X}$ gauge boson achieves a mass via the Stueckelberg mechanism
      2. 2.4.2 $U(1)_{X}$ is broken via a dark Higgs
      3. 2.4.3 Kinetic mixing between a massless $U(1)$ and the hypercharge gauge field
  3. 3 Coupled Boltzmann equations with multiple temperatures
    1. 3.1 Derivation of the evolution of hidden sector particles
    2. 3.2 Coupled Boltzmann equations from two sectors
  4. 4 Evolution of matter in visible and hidden sectors
    1. 4.1 $U(1)_{X}$ sub-GeV dark matter via freeze-in
    2. 4.2 Constraints, benchmark models and phenomenology
    3. 4.3 Evolution of hidden sector particles
  5. 5 Conclusion
  6. A Kinetic mixing and mass mixing in the Stueckelberg extension
  7. B Plasmon Effect
  8. C Relevant scattering cross-sections

##  1 Introduction

The Dirac quantization condition explains the observed quantization of the
electric charge, although the magnetic monopole is not yet discovered. From an
experiment point of view, however, the pursuit of discovering new phenomena
and making new discoveries continues to be the foremost goal. The search for
new particle carrying tiny electric charge, normally referred to as the
millicharge,444From here and the following, we will use the term millicharge
referring to “milli-electric charge”. The tiny amount of other quantum charges
will be otherwise specified. is one of those interesting directions for the
new physics exploration. If millicharge particles exist, the lightest
millicharge particle must be stable and can serve as the dark matter
candidate. Experiments have searched for years and established stringent
bounds for the millicharge dark matter particle, for reviews see [1, 2]. In
this paper, we concentrate on millicharge dark matter in the sub-GeV mass
range, which has been extensively investigated by experimental searches [3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25]. From the theoretical perspective, the exploration of how a consistent
theory can generate such a millicharge, particularly at high energies,
represents an intriguing avenue of research.

The connection between the kinetic mixing and millicharge is traced back to
[26, 27]. The kinetic mixing between a massless $U(1)_{a}$ and a massless
$U(1)_{b}$ can generate a milli-$U(1)_{a}$ charge for a $U(1)_{b}$-charged
fermion, proportional to the kinetic mixing parameter. In a massive amount of
literature, the generation of millicharge particle boils down to a simplified
effective model that an extra $U(1)$ mixed with the photon field. After the
mixing fermions charged under the extra $U(1)$ will carry a millicharge.
However, as we will review in Section 2.2, when one of the $U(1)$’s gets
massive, say $U(1)_{b}$, the kinetic mixing with the massless $U(1)_{a}$ can
no longer generate a milli-$U(1)_{a}$ charge for a $U(1)_{b}$-charged fermion,
no matter how tiny the massive $U(1)_{b}$ is.

In the full electroweak theory, the photon field is a linear combination of
the hypercharge ($B$) and the neutral $SU(2)$ gauge fields $A_{3}$. Thus the
kinetic mixing between the extra $U(1)$ and the hypercharge gauge fields
should be treated simultaneously with the $\\{B,A_{3}\\}$ mass mixing matrix.
It has been shown that pure kinetic mixing between an extra $U(1)$ and the
hypercharge gauge field cannot generate a millicharge [28]. More precisely,
for a general discussion of the Standard Model (SM) extended by a $U(1)_{X}$
hidden sector, the millicharge carried by the $U(1)_{X}$ dark fermions can be
obtained from: (1) The dark particle carries a tiny amount of hypercharge as a
prior. (2) A kinetic mixing between a massless $U(1)$ and the hypercharge
gauge field, and the generated millicharge is proportional to the kinetic
mixing parameter. (3) The mass mixing between a massive $U(1)$ with the
hypercharge gauge field (massive in the initial gauge eigenbasis) [29, 28]. In
the last case the kinetic mixing does not play any role in generating the
millicharge, and the millicharge generated is proportional to the mass mixing
parameter. In this paper we will review all these possibilities and will focus
on the third case, exploring the details in a realistic model involving viable
freeze-in millicharge dark matter from a $U(1)_{X}$ hidden sector. For a
realistic model, it is essential to ensure that in the final mass eigenbasis,
the photon is precisely massless.

The nature of dark matter remains to be a puzzle in modern physics. In
particle physics, dark matter is typically assumed to reside in the hidden
sector, comprising new hypothetical particles and/or new interactions. Hidden
sectors are common in various grand unified theory models and in string
theory. In general, the Universe may encompass a large amount of hidden
sectors that can have direct or indirect, strong or weak connections with the
SM. With direct and indirect detections of dark matter imposing increasingly
stringent constraints on freeze-out scenarios, researchers have turned their
attention to the freeze-in mechanism [30, 31, 32, 33, 34, 35, 36, 37, 38, 39].
While the freeze-in mechanism demands a feeble coupling between the SM and the
dark matter particle, posing a puzzle that why such tiny coupling might emerge
in the theory. The kinetic mixing precisely provide such smallness for the
coupling constant between the SM and the hidden sector [40, 41, 42, 43, 44,
45, 46, 47, 48, 49, 50].

If a hidden sector is connected with the SM with only feeble couplings, the
hidden sector will evolve almost independently and will have its own
temperature. The temperature of the hidden sector may differ significantly
from the temperature of the observed Universe, and this can be analyzed using
the formalism established in [37]. Consequently, our Universe could contain
numerous hidden sectors with temperatures identical or different from the
temperature of the Universe.

When the hidden sector possesses strong interactions among its particles,
calculating the evolution of hidden sector particles becomes challenging,
making it difficult to determine the final relic density of dark matter. In
the literature, numerous assumptions and approximations are employed, and such
calculations cannot be carried out convincingly. The full evolution of a
$U(1)_{X}$ hidden sector with feeble connection to the SM, is first analyzed
in [37], where the number densities of dark particles as well as the hidden
sector temperature can be accurately computed, by solving a set of coupled
Boltzmann equations. In this paper, we concentrate on the extension of the SM
with an extra $U(1)_{X}$ hidden sector, interacting feebly through tiny
kinetic and mass mixings, which may produce sub-GeV millicharge dark matter
(some models also involve dark photon dark matter). We follow the approach
represented in [37], calculating the complete evolution of the $U(1)_{X}$
hidden sector and determining the final relic abundance of the dark matter
candidate(s). We discover several interesting features in the computations,
which could provide valuable insights for future research on hidden sectors
feebly interacting with the SM.

In Section 2, we discuss in general the kinetic and mass mixing effect,
including the simplified model where the photon field mixed directly with an
extra $U(1)$, as well as the general case that the extra $U(1)_{X}$ mixed with
the hypercharge gauge field. In particular, we discuss the kinetic mixing and
mass mixing effects in the full electroweak theory in Section 2.4, with a
detailed derivation shown in Appendix A for the case of Stueckelberg mass
mixing. In Section 3 we review the general formalism of computing the complete
evolution of the hidden sector, and write down the corresponding coupled
Boltzmann equations. A comprehensive investigation on the sub-GeV dark matter
from a $U(1)_{X}$ hidden sector is discussed in Section 4. To this end, we
explore six distinct cases with different evolution details, depending on the
masses of the dark fermion and the dark photon. We discussed various
experimental constraints and carefully chose eight illustrative benchmark
models to show the evolution of hidden sector particles and the hidden sector
temperature. In some of the models, dark photon can also be a valid dark
matter candidate, although it can only occupy at most $\sim 5\%$ of the total
observed dark matter relic density. For the freeze-in production of
millicharge dark particles, the decay from effective massive photon from hot
thermal bath (plasmon) is one of the important production channels and is
reviewed in Appendix B. All relevant cross-sections are summarized in Appendix
C. Finally we conclude in Section 5.

##  2 The $U(1)$ extension of the SM and the generation of millicharge
particle

The SM can be extended by a $U(1)_{X}$ hidden sector mixed with the SM via
gauge kinetic terms and mass terms. In this section, we discuss in general the
effect of both the kinetic mixing and the mass mixing. The simplest $U(1)_{X}$
extension consists of a dark $U(1)_{X}$ gauge boson $C_{\mu}$ and the dark
fermion $\chi$ carrying $U(1)_{X}$ charge $Q_{\chi}=+1$ but no SM quantum
number. In the final mass eigenbasis the dark $U(1)_{X}$ gauge boson is
referred to as the $Z^{\prime}$ gauge boson and is usually called a “dark
photon” if its mass is less than $\mathcal{O}$(GeV). The total Lagrangian of
the $U(1)_{X}$ extended SM is given by

$\mathcal{L}=\mathcal{L}_{{\rm SM}}+\mathcal{L}_{{\rm hid}}+\mathcal{L}_{{\rm
mix}}^{{\rm kin}}+\mathcal{L}_{{\rm mix}}^{{\rm mass}}\,,$ (1)

where $\mathcal{L}_{{\rm hid}}$ is given by

$\mathcal{L}_{{\rm
hid}}=-\frac{1}{4}F_{X\mu\nu}F_{X}^{\mu\nu}+g_{X}\bar{\chi}\gamma^{\mu}\chi
C_{\mu}+m_{\chi}\bar{\chi}\chi\,,$ (2)

where $g_{X}$ is the $U(1)_{X}$ gauge coupling, $C_{\mu}$ is the $U(1)_{X}$
gauge boson, $F_{X\mu\nu}$ is the field strength of $U(1)_{X}$, $\chi$ is the
dark fermion with mass $m_{\chi}$. $\mathcal{L}_{{\rm mix}}^{{\rm kin}}$ is
the kinetic mixing term, written as [26]

$\mathcal{L}^{{\rm kin}}_{{\rm
mix}}=-\frac{\delta}{2}F_{\mu\nu}F_{X}^{\mu\nu}\,,$ (3)

where $F_{\mu\nu}$ is the hypercharge field strength, $\delta$ is the kinetic
mixing parameter. If the $U(1)_{X}$ gauge boson acquires a mass via
Stueckelberg mechanism, $\mathcal{L}_{{\rm mix}}^{{\rm mass}}$ is written as
[28]

$\mathcal{L}_{{\rm mix}}^{{\rm
st}}=-\frac{1}{2}(M_{1}C_{\mu}+M_{2}B_{\mu}+\partial_{\mu}\sigma)^{2}\,,$ (4)

and we define the mass mixing parameter throughout the paper as
$\epsilon=M_{2}/M_{1}$, denoting the mass mixing between the $U(1)_{X}$ and
the hypercharge gauge field. In the Stueckelberg case the gauge field
$C_{\mu}$ is massive without breaking the $U(1)_{X}$ gauge symmetry.

If the $U(1)_{X}$ gauge boson achieves a mass via a dark Higgs $\phi$, the
covariant derivative reads

$D_{\mu}\phi=(\partial_{\mu}-ig_{X}C_{\mu})\,\phi\,,$ (5)

and the hidden $U(1)_{X}$ is broken at some new scale $u$. However, the mixing
between the dark Higgs $\phi$ and the SM Higgs field $H$ via the term
$(\phi^{*}\phi)(H^{\dagger}H)$ will not generate a mass mixing in the gauge
sector, i.e., the mixing of $U(1)_{X}$ with either the hypercharge gauge field
or the neutral $SU(2)$ gauge field is absent. The mass mixing in the gauge
sector can only be generated only if the hidden $U(1)_{X}$ gauge boson
$C_{\mu}$ also receives a tiny SM Higgs mass. One way to realize such scenario
is to gauge a linear combination of the hypercharge and any quantum number
$X$, i.e., $U(1)_{y+X}$ where $y$ is a very small fraction of the hypercharge
$Y$ [51]. Thus the dark particle, in this case, can carry a $U(1)$ charge $y$.
The covariant derivative acting on the SM Higgs doublet $H$ is then written as

$D_{\mu}H=\Big{(}\partial_{\mu}-ig_{2}T^{a}A_{\mu}^{a}-\frac{i}{2}g_{Y}YB_{\mu}-\frac{i}{2}g_{Y}yC_{\mu}\Big{)}\,H\,,$
(6)

where the hidden $U(1)_{X}$ also carries a tiny hypercharge $y$. After the
symmetry breaking and the mixing, the dark fermion originally carried a tiny
hypercharge proportional to $y$ now carries a millicharge. However, this
millicharge is not achieved by either the kinetic mixing or the mass mixing,
c.f., Section 2.4.2 for details. Though there seems to be no evident
explanation for the inclusion of milli-hypercharge in the hidden $U(1)_{X}$ as
a prior.

We outline our primary conclusion at the beginning that the millicharge
particle can be generated by mixing only in two ways: (1) kinetic mixing
between two massless $U(1)$’s [26], and (2) mass mixing between two massive
$U(1)$’s (though in the mass eigenbasis there can still be a massless $U(1)$)
[28]. We emphasize here that except for the two massless $U(1)$ case, the
kinetic mixing does not play any role in generating the millicharge carried by
the dark particle. For the case of two massive $U(1)$’s, if kinetic mixing and
mass mixing are both evoked, the millicharge carried by the fermion will only
depend on the mass mixing.

Since the photon is a linear combination of the hypercharge gauge field and
the neutral $SU(2)$ gauge field in the electroweak theory, the dark $U(1)_{X}$
should mix with the hypercharge at high energies rather than with the photon
field. For the analysis in the early universe, it is important to explore how
the effective kinetic mixing term of photon and dark photon is generated [52,
53], and further, how the millicharge carried by the dark particle is
generated at low energies. It is also intriguing to explore the stringy origin
of the kinetic mixing [54, 55, 56] as well as the production of the
millicharge [57, 58, 59]. After the discussion of this section, one will
observe that in the case of a massive $U(1)_{X}$ mixed with the hypercharge
gauge field in the full electroweak theory, the kinetic mixing effect cannot
generate a millicharge in any case. The millicharge can only be generated by
(1) the Stueckelberg mass mixing or (2) the dark particle originally carrying
a tiny hypercharge.

In the following we would like to summarize some of the known results
regarding the generation of millicharge particles, and we wish to clarify
several important issues as well as some misunderstandings regarding the
generation of the millicharge with kinetic and/or mass mixing effect. We will
first review the mixing effect for two $U(1)$’s and then discuss the mixing
effect of an extra $U(1)_{X}$ with the full electroweak sector.

###  2.1 Kinetic mixing between two massless $U(1)$’s

The simplest case of kinetic mixing arises from two massless $U(1)$’s [26,
27]. From an effective theory approach, at low energies one can consider the
case that an extra $U(1)_{X}$ mixed with the photon field. In the gauge
eigenbasis, the kinetic part and interaction part of the Lagrangian read

$\displaystyle\mathcal{L}_{{\rm kin}}$
$\displaystyle=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{4}F_{X\,\mu\nu}F_{X}^{\mu\nu}-\frac{\bm{\delta}}{2}F_{\mu\nu}F_{X}^{\mu\nu}\,,$
(7) $\displaystyle\mathcal{L}_{{\rm int}}$ $\displaystyle=eA_{\mu}J_{{\rm
em}}^{\mu}+g_{X}C_{\mu}J_{{\rm d}}^{\mu}\,,$ (8)

where $\bm{\delta}$ is the kinetic mixing parameter,555 In this paper we use
the boldface $\bm{\delta}$ to denote the kinetic mixing between the the extra
$U(1)_{X}$ and $U(1)_{\rm em}$ in the simplified model, and use the regular
$\delta$ to denote kinetic mixing between the extra $U(1)_{X}$ and the
hypercharge gauge field in the full electroweak theory. $e$ is the
electromagnetic coupling constant, i.e., the electric charge and $g_{X}$ is
the dark $U(1)_{X}$ gauge coupling constant, $J_{{\rm em}}^{\mu}$ and $J_{{\rm
d}}^{\mu}$ are the electromagnetic current and dark current respectively. To
keep the gauge kinetic terms in the canonical form, one performs a non-unitary
transformation

$K=\left(\begin{array}[]{cc}c_{\bm{\delta}}&0\\\
-s_{\bm{\delta}}&1\end{array}\right)\,,$ (9)

where $c_{\bm{\delta}}=1/\sqrt{1-\bm{\delta}^{2}}$ and
$s_{\bm{\delta}}=\bm{\delta}/\sqrt{1-\bm{\delta}^{2}}$. After the
diagonalization, the gauge kinetic terms are in the canonical form with the
mixing term eliminated, and the fields transform from the gauge eigenbasis
$\\{A,C\\}$ to the physical eigenbasis $\\{A_{\gamma},A_{X}\\}$ including the
physical photon and the dark photon as

$\left(\begin{array}[]{c}A\\\
C\end{array}\right)=\left(\begin{array}[]{cc}c_{\bm{\delta}}&0\\\
-s_{\bm{\delta}}&1\end{array}\right)\left(\begin{array}[]{c}A_{\gamma}\\\
A_{X}\end{array}\right)\,.$ (10)

In the physical eigenbasis, the interactions can be rewritten as

$\mathcal{L}_{{\rm int}}=ec_{\bm{\delta}}A_{\gamma}J_{{\rm
em}}+g_{X}(-s_{\bm{\delta}}A_{\gamma}+A_{X})J_{{\rm d}}\,.$ (11)

One can see now the dark current can couple to photon with strength
$s_{\bm{\delta}}g_{X}$, corresponding to a millicharge
$s_{\bm{\delta}}g_{X}/e$ carried by $\chi$ (we assume $\chi$ carries the
$U(1)_{X}$ charge $+1$). However, the transformation Eq. (9) is arbitrary, one
can reverse $A,C$ in Eq. (10) as

$\left(\begin{array}[]{c}C\\\
A\end{array}\right)=\left(\begin{array}[]{cc}c_{\bm{\delta}}&0\\\
-s_{\bm{\delta}}&1\end{array}\right)\left(\begin{array}[]{c}A_{X}\\\
A_{\gamma}\end{array}\right)\,.$ (12)

The interactions in this case are then written as

$\mathcal{L}_{{\rm int}}=g_{X}c_{\bm{\delta}}A_{X}J_{{\rm
d}}+e(-s_{\bm{\delta}}A_{X}+A_{\gamma})J_{{\rm em}}\,,$ (13)

and one can see now the SM fermions carry milli-$U(1)_{X}$ charges
$s_{\bm{\delta}}eQ_{i}/g_{X}$, where $Q_{i}$ denotes the electric charges
carried by SM fermions in the unit of $e$, though this transformation is
phenomenologically not interesting.

More generally, there is a rotation degree of freedom for the transformation
Eq. (9). One can always multiply an orthogonal matrix

$R=\left(\begin{array}[]{cc}\cos\alpha&-\sin\alpha\\\
\sin\alpha&\cos\alpha\end{array}\right)\equiv\left(\begin{array}[]{cc}c_{\alpha}&-s_{\alpha}\\\
s_{\alpha}&c_{\alpha}\end{array}\right)$ (14)

to the non-unitary transformation $K$, which gives rise to the interaction
terms as

$\mathcal{L}_{{\rm
int}}=e(c_{\bm{\delta}}c_{\alpha}A_{\gamma}-c_{\bm{\delta}}s_{\alpha}A_{X})J_{{\rm
em}}+g_{X}\big{[}(-c_{\alpha}s_{\bm{\delta}}+s_{\alpha})A_{\gamma}+(c_{\alpha}+s_{\bm{\delta}}s_{\alpha})A_{X}\big{]}J_{{\rm
d}}\,.$ (15)

The case $\alpha=0$ which corresponding to no rotation goes back to Eq. (11),
the case $\alpha=\pi/2$ corresponds to the reverse case Eq. (13) with the
redefinition of the physical fields $A_{X}\leftrightarrow A_{\gamma}$. In
practice, this rotation degree of freedom is redundant because the physical
photon is measured to be the gauge boson that couples to the EM current
$J_{{\rm em}}$. Namely one should redefine a real (measured) photon field
$A_{\gamma}^{{\rm real}}$ and $A_{X}^{{\rm real}}$ as

$\displaystyle A_{\gamma}^{{\rm real}}$ $\displaystyle\equiv
c_{\alpha}A_{\gamma}-s_{\alpha}A_{X}\,,$ (16) $\displaystyle A_{X}^{{\rm
real}}$ $\displaystyle\equiv s_{\alpha}A_{\gamma}+c_{\alpha}A_{X}\,.$ (17)

Under such field redefinition, one finds the interactions in the physical
eigenbasis as

$\mathcal{L}_{{\rm int}}=ec_{\bm{\delta}}A_{\gamma}^{{\rm real}}J_{{\rm
em}}+g_{X}(-s_{\bm{\delta}}A_{\gamma}^{{\rm real}}+A_{X}^{{\rm real}})J_{{\rm
d}}\,,$ (18)

which is exactly in the same form as Eq. (11).

In summary, the kinetic mixing of two massless $U(1)$ gauge fields can
generate the dark particle carrying a millicharge proportional to the kinetic
mixing parameter $\bm{\delta}$.

###  2.2 Kinetic mixing between a massless $U(1)$ and a massive $U(1)$

A more realistic case is kinetic mixing between the photon field with a
massive $U(1)$. In this case the Lagrangian reads

$\displaystyle\mathcal{L}_{{\rm eff}}$
$\displaystyle=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{4}F_{X\,\mu\nu}F_{X}^{\mu\nu}-\frac{\bm{\delta}}{2}F_{\mu\nu}F_{X}^{\mu\nu}-\frac{1}{2}M^{2}C^{2}\,,$
(19) $\displaystyle\mathcal{L}_{{\rm int}}$ $\displaystyle=eA_{\mu}J_{{\rm
em}}^{\mu}+g_{X}C_{\mu}J_{{\rm d}}^{\mu}\,.$ (20)

Now the only possible way of eliminating the kinetic mixing term is given by
the following transformation

$\left(\begin{array}[]{c}C\\\
A\end{array}\right)=\left(\begin{array}[]{cc}c_{\bm{\delta}}&0\\\
-s_{\bm{\delta}}&1\end{array}\right)\left(\begin{array}[]{c}A_{X}\\\
A_{\gamma}\end{array}\right)\,,$ (21)

similar to Eq. (12), to maintain the physical photon field exactly massless.
In other words, the nonzero mass of the massive $U(1)$ field fixes the angle
$\alpha$ to be $\pi/2$ and thus fixes the rotation matrix $R$ in Eq. (14). The
above transformation gives rise to the interaction in the physical eigenbasis
as

$\mathcal{L}_{{\rm int}}=g_{X}c_{\bm{\delta}}A_{X}J_{{\rm
d}}+e(-s_{\bm{\delta}}A_{X}+A_{\gamma})J_{{\rm em}}\,,$ (22)

and one can see in this case, the dark particle carries exactly zero electric
charge.

Hence, it is not appropriate to consider a kinetic mixing between the photon
and a massive dark photon, and still claiming the dark particles from the dark
$U(1)$ sector can carry a millicharge proportional to the kinetic mixing
parameter $\bm{\delta}$. More precisely, the case of kinetic mixing between
two massless $U(1)$’s, is distinct from the situation where a massless $U(1)$
kinetically mixed with a massive $U(1)$, even when the latter has an extremely
small but nonzero mass $M$. One cannot take the $M\to 0$ limit in this case
and use directly the result from the case of kinetic mixing between two
massless $U(1)$’s. In summary, for a dark photon with a non-negligible mass
(for example, a $\mathcal{O}$(keV–MeV) dark photon as commonly studied), even
though this mass is much smaller than the experimental scale one considered,
the kinetic mixing between a massless $U(1)$ and a massive $U(1)$ does not
generate a millicharge carried by the dark particle.

As long as one of the $U(1)$’s is massive, the only possible way of generating
the millicharge through mixing is from the mass mixing effect, and we will
review this mechanism in the next subsection.

###  2.3 Kinetic mixing and mass mixing between two massive $U(1)$’s

The Stueckelberg mass mixing effect of two massive $U(1)$ is discussed in
[28]. The full Lagrangian reads

$\displaystyle\mathcal{L}_{{\rm kin}}$
$\displaystyle=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{4}F_{X\,\mu\nu}F_{X}^{\mu\nu}-\frac{\bm{\delta}}{2}F_{\mu\nu}F_{X}^{\mu\nu}\,,$
(23) $\displaystyle\mathcal{L}_{{\rm mass}}$
$\displaystyle=-\frac{1}{2}M_{2}^{2}A^{2}-\frac{1}{2}M_{1}^{2}C^{2}-M_{1}M_{2}AC\,,$
(24) $\displaystyle\mathcal{L}_{{\rm int}}$ $\displaystyle=eA_{\mu}J_{{\rm
em}}^{\mu}+g_{X}C_{\mu}J_{{\rm d}}^{\mu}\,.$ (25)

To obtain the physical eigenbasis, one needs to diagonalize the kinetic mixing
matrix and mass mixing matrix simultaneously for the gauge eigenbasis
$V^{T}=(C,A)$,

$\mathcal{K}=\left(\begin{array}[]{cc}1&\bm{\delta}\\\
\bm{\delta}&1\end{array}\right)\,,\qquad M_{{\rm
St}}^{2}=\left(\begin{array}[]{cc}M_{1}^{2}&M_{1}M_{2}\\\
M_{1}M_{2}&M_{2}^{2}\end{array}\right)=M_{1}^{2}\left(\begin{array}[]{cc}1&\epsilon\\\
\epsilon&\epsilon^{2}\end{array}\right)\,,$ (26)

where we define the mass mixing parameter $\epsilon=M_{2}/M_{1}$.

Firstly we know that the transformation $K$ of Eq. (9) takes $\mathcal{K}$ to
be diagonal,

$K^{T}\mathcal{K}K=1\,,$ (27)

and sets a new basis as

$KV^{\prime}=V=(C,A)^{T}\,,\qquad{\rm or}\qquad V^{\prime}=K^{-1}V\,.$ (28)

The gauge kinetic terms then become canonical with the new basis $V^{\prime}$.
At the same time the new basis $KV^{\prime}$ also changes the mass terms to be

$V^{T}M_{{\rm St}}^{2}V\to(KV^{\prime})^{T}M_{{\rm
St}}^{2}(KV^{\prime})=V^{\prime T}K^{T}M_{{\rm St}}^{2}KV^{\prime}\,.$ (29)

Performing one further orthogonal transformation $O$ to diagonalize the matrix
$K^{T}M_{{\rm St}}^{2}K$ such that

$O^{T}K^{T}M_{{\rm St}}^{2}KO=D^{2}\,,$ (30)

where the matrix $D$ is the final diagonalized mass matrix. The mass terms
reduce to

$V^{\prime T}K^{T}M_{{\rm St}}^{2}KV^{\prime}=V^{\prime T}O(O^{T}K^{T}M_{{\rm
St}}^{2}KO)O^{T}V^{\prime}=(V^{\prime T}O)D^{2}(O^{T}V^{\prime})\equiv
E^{T}D^{2}E\,,$ (31)

where $E$ is the final mass eigenbasis. The kinetic terms will still stay in
the canonical form since $V^{\prime T}V^{\prime}=V^{\prime
T}OO^{T}V^{\prime}=E^{T}E$. Thus the transformation $R=KO$ diagonalizes the
kinetic and mass mixing matrices simultaneously, and it transforms the
original gauge eigenbasis $V^{T}=(C,A)$ to the mass eigenbasis
$E^{T}=(A_{X},A_{\gamma})$, which is

$V=RE=(KO)E\,.$ (32)

One should keep in mind that the physical photon must be exactly massless,
i.e., the mass term is zero for $A_{\gamma}$. One finds

$\tan\theta=\frac{\epsilon\sqrt{1-\bm{\delta}^{2}}}{1-\epsilon\bm{\delta}}\,,$
(33)

which leads to the following interaction terms

$\displaystyle\mathcal{L}_{{\rm int}}$ $\displaystyle=(g_{X}J_{{\rm
d}}^{\mu},eJ_{{\rm em}}^{\mu})V=(g_{X}J_{{\rm d}}^{\mu},eJ_{{\rm
em}}^{\mu})RE$
$\displaystyle=\frac{1}{\sqrt{1-2\epsilon\bm{\delta}+\epsilon^{2}}}\frac{1}{\sqrt{1-\bm{\delta}^{2}}}A_{X}\big{[}(1-\epsilon\bm{\delta})g_{X}J_{{\rm
d}}+eJ_{{\rm em}}(\epsilon-\bm{\delta})\big{]}$
$\displaystyle+\frac{1}{\sqrt{1-2\epsilon\bm{\delta}+\epsilon^{2}}}A_{\gamma}(eJ_{{\rm
em}}-\epsilon g_{X}J_{{\rm d}})\,.$ (34)

In this case, one can see clearly that the millicharge can be induced only
when the mass mixing, denoted by $\epsilon$, is non-vanishing, and is not
related to the kinetic mixing parameter $\bm{\delta}$ at all. One can also
see, for the case that $M_{2}=0$ or $\epsilon=0$, the above result reduces to
the result in the previous subsection, i.e., a massless $U(1)$ mixed with a
massive $U(1)$, where again no millicharge is produced.

###  2.4 Kinetic mixing and mass mixing in the full electroweak theory

In the full electroweak theory, the Higgs mechanism generates the mass mixing
between the hypercharge gauge field $B$ and the neutral $SU(2)$ gauge field
$A_{3}$, and the photon is a linear combination of them. Hence, the extra
$U(1)_{X}$ should mix with the hypercharge rather than the photon field at
high energies. For a realistic model beyond the SM, the extra $U(1)_{X}$ can
mix with the hypercharge via kinetic and mass terms. To obtain the
$\\{A^{\prime},A^{\gamma},Z\\}$ mass eigenstates, where $A^{\prime}$ is the
additional gauge boson in the theory, either called a $Z^{\prime}$ or a dark
photon, the $A_{3}$ gauge boson is always entangled in the mass matrix after
the electroweak phase transition when Higgs picks a vev.

In this subsection we discuss realist models that can generate a millicharge
carried by the dark particle. To this end, we discuss the kinetic mixing and
mass mixing effects of an extra $U(1)_{X}$ with the full electroweak sector.
This encompasses three different situations, and we will see that the kinetic
mixing between $U(1)_{Y}$ and a massive $U(1)_{X}$ alone cannot produce a
millicharge for the dark fermion.

####  2.4.1 $U(1)_{X}$ gauge boson achieves a mass via the Stueckelberg
mechanism

In this subsection we review the kinetic mixing and mass mixing in the
Stueckelberg extension of the SM, c.f., Eqs. (3) and (4), first discussed in
[28]. The Stueckelberg mass mixing effect is discussed in details in the
context of string theory in [58, 59]. Within the context of string theory, the
quantization of the electric charge, i.e., any charge should be an integral
(more accurately fractional number) multiple of the electric charge, can be
explained that the millicharge carried by the dark particles corresponds to
the ratio of integer or half-integer wrapping numbers of D-branes which wrap
on cycles in the internal six dimensional manifold in string theory.

In the gauge eigenbasis of the $U(1)_{X}$, hypercharge and the neutral $SU(2)$
gauge field $V^{T}=(C,B,A^{3})$, the mixing matrices can be written as

$\mathcal{K}=\left(\begin{array}[]{ccc}1&\delta&0\\\ \delta&1&0\\\
0&0&1\end{array}\right)\,,\qquad M_{{\rm
St}}^{2}=\left(\begin{array}[]{ccc}M_{1}^{2}&M_{1}M_{2}&0\\\
M_{1}M_{2}&M_{2}^{2}+\frac{1}{4}v^{2}g_{Y}^{2}&-\frac{1}{4}v^{2}g_{2}g_{Y}\\\
0&-\frac{1}{4}v^{2}g_{2}g_{Y}&\frac{1}{4}v^{2}g_{2}^{2}\end{array}\right)\,.$
(35)

Similar to the discussion for the two $U(1)$ case in previous subsection, a
simultaneous diagonalization is performed for the above two mixing matrices.
This tackles down to the problem of diagonalizing the matrix $K^{T}M_{{\rm
St}}^{2}K$ using an orthogonal transformation $O$ such that $O^{T}K^{T}M_{{\rm
St}}^{2}KO$ is diagonal, where

$K=\left(\begin{array}[]{ccc}c_{\delta}&0&0\\\ -s_{\delta}&1&0\\\
0&0&1\end{array}\right)\,,$ (36)

with $c_{\delta}=1/\sqrt{1-\delta^{2}}$ and
$s_{\delta}=\delta/\sqrt{1-\delta^{2}}$. The original gauge eigenbasis
$V^{T}=(C,B,A^{3})$ then transforms to the mass eigenbasis
$E^{T}=(A^{\prime},A^{\gamma},Z)$ following the relation

$V=\mathcal{R}E=(KO)E\,,$ (37)

the detailed calculation and the transformation matrix $\mathcal{R}$ can be
found in Appendix A. It will be useful to define the mass mixing parameter
$\epsilon=M_{2}/M_{1}$. The interaction of gauge bosons with fermions can be
obtained from

$\mathcal{L}_{{\rm
int}}=\left(g_{X}J_{X},g_{Y}J_{Y},g_{2}J_{3}\right)V=\left(g_{X}J_{X},g_{Y}J_{Y},g_{2}J_{3}\right)\mathcal{R}E\,.$
(38)

Here we summarize the neutral current couplings. The interactions of
$A^{\prime},Z,A^{\gamma}$ to SM fermions can be written as

$-\mathcal{L}_{A^{\prime},Z,A^{\gamma}}=\frac{1}{2}\bar{f}_{i}\gamma^{\mu}\left[\left(v_{i}^{\prime}-a_{i}^{\prime}\gamma^{5}\right)f_{i}A_{\mu}^{\prime}+\left(v_{i}-a_{i}\gamma^{5}\right)f_{i}Z_{\mu}\right]+eQ_{i}\bar{f}_{i}\gamma^{\mu}f_{i}A_{\mu}^{\gamma}\,,$
(39)

where

$\displaystyle v_{i}$
$\displaystyle=\left(g_{2}\mathcal{R}_{33}-g_{Y}\mathcal{R}_{23}\right)T_{i}^{3}+2g_{Y}\mathcal{R}_{23}Q_{i}\,,$
(40) $\displaystyle a_{f}$
$\displaystyle=\left(g_{2}\mathcal{R}_{33}-g_{Y}\mathcal{R}_{23}\right)T_{i}^{3}\,,$
(41) $\displaystyle v_{i}^{\prime}$
$\displaystyle=\left(g_{2}\mathcal{R}_{31}-g_{Y}\mathcal{R}_{21}\right)T_{i}^{3}+2g_{Y}\mathcal{R}_{21}Q_{i}\,,$
(42) $\displaystyle a_{f}^{\prime}$
$\displaystyle=\left(g_{2}\mathcal{R}_{31}-g_{Y}\mathcal{R}_{21}\right)T_{i}^{3}\,.$
(43)

In the above equations, $Q_{i}e$ is the electric charges of the SM fermions,
$g_{Y}$ is the original hypercharge gauge coupling in the electroweak theory.

The electric charge is altered by the mixing effects and is written as

$e=\frac{g_{2}g_{Y}^{\prime}}{\sqrt{g_{2}^{2}+g_{Y}^{\prime 2}}}\,,$ (44)

where the hypercharge gauge coupling is also modified to be
$g_{Y}^{\prime}=g_{Y}/\sqrt{1+\epsilon^{2}-2\delta\epsilon}$.

When the mass of the dark photon is less than twice of dark fermion mass and
is also less than twice of electron mass, the dark photon can only decay to
pair of neutrino and anti-neutrino, and to three photons, as shown in Fig. 1.
The coupling of the dark photon to neutrinos is calculated to be

$\displaystyle-\mathcal{L}_{A^{\prime}\bar{\nu}\nu}$
$\displaystyle=\frac{1}{2}a_{\nu}^{\prime}\bar{\nu}_{L}\gamma^{\mu}\nu_{L}A_{\mu}^{\prime}\approx-\frac{1}{2}g_{Y}^{\prime}\left(\delta-\epsilon\right)\epsilon_{0}^{2}\bar{\nu}_{L}\gamma^{\mu}\nu_{L}A_{\mu}^{\prime}\,,$
(45)

where $M_{0}=v\sqrt{g_{2}^{2}+g_{Y}^{\prime 2}}/2\approx M_{Z}$, almost
identical to the original $Z$ boson mass in the electroweak theory for small
$\delta,\epsilon$, and we define
$\epsilon_{0}=M_{\gamma^{\prime}}/M_{0}\approx M_{1}/M_{Z}$.

Figure 1:  A display of dark photon decay channels for
$M_{\gamma^{\prime}}<2m_{e}$, including the decay to pair of neutrino and
anti-neutrino, and to three photons. Considering various constraints, the
decay of the dark photon to neutrinos due to the mixing effect is always
suppressed compared to the three-photon decay channel.

The decay width of the dark photon to neutrinos is given by

$\Gamma_{\gamma^{\prime}\nu\bar{\nu}}\sim\frac{g_{Y}^{\prime
2}(\delta-\epsilon)^{2}\epsilon_{0}^{4}}{48\pi}M_{\gamma^{\prime}}\,.$ (46)

For a 100 keV dark photon where we have $\epsilon_{0}\sim 10^{-6}$, take the
largest value of $|\delta-\epsilon|\sim 10^{-12}$, one estimates that
$\Gamma_{\gamma^{\prime}\nu\bar{\nu}}\sim 10^{-55}$ GeV, which is much smaller
than the dark photon to three photon decay width, given by

$\Gamma_{\gamma^{\prime}\to
3\gamma}\sim\frac{(\delta-\epsilon)^{2}\alpha^{4}M_{\gamma^{\prime}}^{9}}{m_{e}^{8}}\,.$
(47)

The above decay width gives rise to the lifetime of the dark photon to be
roughly $10^{18}$ second, greater than the age of the Universe. While we are
going to see in Section 4 that even the dark photon is stable during the age
of the Universe, its couplings to the SM particles will be further constrained
by the isotropic diffuse photon background (IDPB) data.

The interactions of $A^{\prime},A^{\gamma},Z$ to the dark fermion $\chi$ can
be written as

$-\mathcal{L}_{\chi}=g_{X}Q_{\chi}\left(\mathcal{R}_{11}A_{\mu}^{\prime}+\mathcal{R}_{12}A_{\mu}^{\gamma}+\mathcal{R}_{13}Z_{\mu}\right)\bar{\chi}\gamma^{\mu}\chi\,.$
(48)

After the mixing, the dark fermion $\chi$ will couple to the photon with

$-\mathcal{L}_{\chi}^{{\rm em}}=Q_{\epsilon}\bar{\chi}\gamma^{\mu}\chi
A_{\mu}^{\gamma}=\epsilon g_{X}Q_{\chi}\cos\theta\,\bar{\chi}\gamma^{\mu}\chi
A_{\mu}^{\gamma}\,,$ (49)

where $Q_{\epsilon}\equiv\epsilon g_{X}Q_{\chi}\cos\theta$ is the millicharge
carried by $\chi$.

One can see clearly from the Eq. (49) that the millicharge carried by the dark
fermion $\chi$ is generated by the mass mixing effect, and thus only depends
on the mass mixing parameter $\epsilon$. Once the mass mixing is turned off
($\epsilon=0$), no millicharge will be generated. Thus the generation of
millicharge has nothing to do with the kinetic mixing between $U(1)_{X}$ and
the hypercharge gauge field.

Finally we summarize the order of the couplings induced by the kinetic mixing
and the mass mixing:

$\displaystyle A^{\gamma}i\bar{i}\sim eQ_{i}\,,\qquad
A^{\gamma}\chi\bar{\chi}\sim\epsilon g_{X}\,,$ (50) $\displaystyle
Zi\bar{i}\sim g_{2}\,,\qquad Z\chi\bar{\chi}\sim\delta g_{X}\,,$ (51)
$\displaystyle A^{\prime}i\bar{i}\sim|\delta-\epsilon|g_{Y}\,,\qquad
A^{\prime}\chi\bar{\chi}\sim g_{X}\,.$ (52)

####  2.4.2 $U(1)_{X}$ is broken via a dark Higgs

If the $U(1)_{X}$ is broken at some scale and achieves a mass via the dark
$U(1)_{X}$ Higgs $\phi$, the dark $U(1)_{X}$ breaking scale is important.
Prior to the $U(1)_{X}$ breaking, the massless $U(1)_{X}$ mixed with the
massless hypercharge gauge field via gauge kinetic terms. After the mixing the
dark fermion $\chi$ carrying the $U(1)_{X}$ charge, carries a milli-
hypercharge proportional to the kinetic mixing parameter $\delta$, as
described in Section 2.1. As the temperature of the universe drops down, and
considering the case that the $U(1)_{X}$ breaking scale is higher than the
electroweak phase transition, in the stage of $U(1)_{X}$ is broken while
before the electroweak phase transition, the massive $U(1)_{X}$ mixed with the
massless hypercharge via kinetic terms. As discussed in Section 2.2, the dark
particles carry zero hypercharge. When both of the $U(1)_{X}$ and electroweak
symmetry are broken, the $U(1)_{X}$ gauge boson will mix with the electroweak
sector and the mixing matrices become $3\times 3$, c.f., Eqs. (5) and (6)

$\mathcal{K}=\left(\begin{array}[]{ccc}1&\delta&0\\\ \delta&1&0\\\
0&0&1\end{array}\right)\,,\qquad
M^{2}=\left(\begin{array}[]{ccc}g_{X}^{2}u^{2}+\frac{1}{4}y^{2}g_{Y}^{2}v^{2}&\frac{1}{4}yg_{Y}^{2}v^{2}&-\frac{1}{4}yg_{2}g_{Y}v^{2}\\\
\frac{1}{4}yg_{Y}^{2}v^{2}&\frac{1}{4}g_{Y}^{2}v^{2}&-\frac{1}{4}g_{2}g_{Y}v^{2}\\\
-\frac{1}{4}yg_{2}g_{Y}v^{2}&-\frac{1}{4}g_{2}g_{Y}v^{2}&\frac{1}{4}g_{2}^{2}v^{2}\end{array}\right)\,.$
(53)

For the case of a general gauged $U(1)_{y+X}$ theory, where $y$ is a very
small fraction of the hypercharge, studied in for example [51], the dark
particles can carry a tiny amount of hypercharge proportional to $y$. Finally,
in the mass eigenbasis the dark particles then carry a tiny electric charge
proportional to $y$. In this case the kinetic mixing and mass mixing does not
play any role of generating such millicharge.

####  2.4.3 Kinetic mixing between a massless $U(1)$ and the hypercharge gauge
field

In this subsection we discuss the case that the kinetic mixing between a
massless $U(1)_{X}$ and the hypercharge gauge field, which can generate a
millicharge proportional to the kinetic mixing parameter. Now the kinetic
mixing matrix together with the mass matrix corresponding to the gauge
eigenbasis $V^{T}=\left(B,C,A_{3}\right)$ are given by

$\mathcal{K}=\left(\begin{array}[]{ccc}1&\delta&0\\\ \delta&1&0\\\
0&0&1\end{array}\right)\,,\qquad
M^{2}=\left(\begin{array}[]{ccc}\frac{1}{4}g_{Y}^{2}v^{2}&0&-\frac{1}{4}g_{2}g_{Y}v^{2}\\\
0&0&0\\\
-\frac{1}{4}g_{2}g_{Y}v^{2}&0&\frac{1}{4}g_{2}^{2}v^{2}\end{array}\right)\,.$
(54)

Using the formalism we have discussed in the previous subsection,
diagonalizing both of the above two matrices simultaneously, is equivalent to
diagonalizing a single matrix $K^{T}M^{2}K$, with matrix $K$ given in Eq.
(36). In this case, this is done by using an orthogonal transformation $O$
given by

$O=\left(\begin{array}[]{ccc}\frac{-\sqrt{1-\delta^{2}}g_{2}}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}&0&\frac{g_{Y}}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}\\\
0&1&0\\\
\frac{-g_{Y}}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}&0&\frac{-\sqrt{1-\delta^{2}}g_{2}}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}\end{array}\right)\,,$
(55)

such that $\left(KO\right)^{T}M^{2}KO$ is diagonal

$\left(KO\right)^{T}M^{2}KO=\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&0\\\
0&0&\frac{v^{2}}{4}\left(g_{2}^{2}+\frac{1}{1-\delta^{2}}g_{Y}^{2}\right)\end{array}\right)\,,$
(56)

which gives rise to only one nonzero eigenvalue for the $Z$ boson mass-square.

The full transformation matrix is then computed as

$\mathcal{R}^{0}=KO=\left(\begin{array}[]{ccc}-\frac{g_{2}}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}&0&\frac{g_{Y}}{\sqrt{1-\delta^{2}}\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}\\\
\frac{g_{2}\,\delta}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}&1&-\frac{g_{Y}\,\delta}{\sqrt{1-\delta^{2}}\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}\\\
-\frac{g_{Y}}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}&0&-\frac{g_{2}\sqrt{1-\delta^{2}}}{\sqrt{\left(1-\delta^{2}\right)g_{2}^{2}+g_{Y}^{2}}}\end{array}\right)\,.$
(57)

The interaction of three neutral gauge bosons to the dark fermion $\chi$ is
thus given by

$-\mathcal{L}_{\chi}=g_{X}Q_{\chi}\left(\mathcal{R}^{0}_{21}A_{\mu}^{\gamma}+\mathcal{R}^{0}_{22}A_{\mu}^{\prime}+\mathcal{R}^{0}_{23}Z_{\mu}\right)\bar{\chi}\gamma^{\mu}\chi\,.$
(58)

Specifically, the coupling of the dark fermion $\chi$ to the photon reads

$-\mathcal{L}_{\chi}^{\mathrm{em}}=g_{X}Q_{\chi}\mathcal{R}^{0}_{21}A_{\mu}^{\gamma}\bar{\chi}\gamma^{\mu}\chi\approx\frac{\delta
g_{2}g_{X}Q_{\chi}}{\sqrt{g_{2}^{2}+g_{Y}^{2}}}A_{\mu}^{\gamma}\bar{\chi}\gamma^{\mu}\chi\,,$
(59)

generating a millicharge proportional to the kinetic mixing parameter, carried
by the dark fermion.

##  3 Coupled Boltzmann equations with multiple temperatures

In Fig. 2, we graphically illustrate the model we discuss. Initially there is
negligible amount of matter in the hidden sector. The dark fermion $\chi$ and
the dark photon $\gamma^{\prime}$ are produced by the freeze-in production of
the SM particles. Once the matter inside the hidden sector accumulated, hidden
sector interactions start to take place, such as
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ which will
change the number densities of $\chi$ and $\gamma^{\prime}$. The hidden sector
interactions occur simultaneously with the freeze-in processes. The hidden
sector possesses a different temperature $T_{h}$, which is related to the
visible sector temperature $T_{v}$ (the temperature of the observed Universe,
and we will omit the subscript v from here), with the function $\eta(T_{h})$.
When the dark photon mass is greater than twice of electron mass, the dark
photon will eventually decay to (mostly) electron-positron pairs. When the
dark photon mass is less than twice of electron mass, it can only decay into
two neutrinos and to three photons and must be made stable through the age of
the Universe, see discussions in Section 2.4.1.

Figure 2: A graphic illustration of the model we discuss. The $U(1)_{X}$
hidden sector connects to the SM with tiny kinetic and mass mixings, and thus
it has only feeble interactions with SM particles. The hidden sector possesses
a different temperature $T_{h}$, related to the visible sector temperature
$T_{v}$ with the function $\eta(T_{h})$. The dark fermion $\chi$ and the dark
photon $\gamma^{\prime}$ are produced through the freeze-in of the SM
particles. Inside the hidden sector, $\chi$ and $\gamma^{\prime}$ can have
either strong or weak (maybe ultraweak) interactions such as
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$. The hidden
sector interactions and the freeze-in productions occur simultaneously and are
related to $T_{h}$ and $T_{v}$ respectively, and must be addressed jointly.

As shown in [37], for the cases that SM extended by one or multiple hidden
sectors interacting feebly with the visible sector, the hidden sectors will
evolve by themselves, especially they can possess temperatures distinct from
the visible sector temperature. In this section we will review the
temperature-dependent coupled Boltzmann equations of the millicharge dark
particle $\chi$ and the dark photon $\gamma^{\prime}$ for the model described
in Section 2.4.1, as well as the evolution equation of $\eta(T_{h})=T/T_{h}$,
with $T\,(T_{h})$ being the visible (hidden) sector temperature. By solving
these coupled differential equations simultaneously, one obtains the evolution
of hidden sector particle number densities and the hidden sector temperature.
In the analysis we notice that the entropies in the hidden and the visible
sectors are not separately preserved, but their sum is conserved, which is
imposed in the analysis. Since we are mostly interested in the evolution of
the hidden sector, we will use the hidden sector temperature $T_{h}$ as the
clock and the temperature in the visible sector is related to the hidden
sector via the function $\eta$.

###  3.1 Derivation of the evolution of hidden sector particles

In this subsection we review the derivation of the evolution of hidden sector
temperature [37]. We start from the two Friedmann equations for a flat
universe

$\displaystyle H^{2}$ $\displaystyle=\frac{8\pi G_{N}}{3}\rho\,,$ (60)
$\displaystyle\frac{\ddot{a}}{a}$ $\displaystyle=-\frac{4\pi
G_{N}}{3}(\rho+3p)\,,$ (61)

where $G_{N}$ is Newton’s gravitational constant, $\rho$ and $p$ are
respectively the energy density and pressure, give rise to the continuity
equation

$\frac{{\rm d}\rho}{{\rm d}t}+3H(\rho+p)=0\,.$ (62)

By using $T_{h}$ as the clock, we can then obtain the following relation from
Eq. (62)

$\frac{{\rm d}T_{h}}{{\rm d}t}=-\frac{4\zeta\rho}{\frac{{\rm d}\rho}{{\rm
d}T_{h}}}H\,,$ (63)

where $\zeta=\frac{3}{4}(1+p/\rho)$. Here $\zeta=1$ is for the radiation
dominated era and $\zeta=3/4$ for the matter dominated universe. We assume
that the hidden sector interacts with the visible sector feebly and the
continuity equations for the two sectors are modified to be

$\displaystyle\frac{{\rm d}\rho_{h}}{{\rm d}t}+3H\left(\rho_{h}+p_{h}\right)$
$\displaystyle=j_{h}\,,$ (64) $\displaystyle\frac{{\rm d}\rho_{v}}{{\rm
d}t}+3H\left(\rho_{v}+p_{v}\right)$ $\displaystyle=-j_{h}\,,$ (65)

where the subscripts $v$ and $h$ correspond to the visible and hidden sectors,
respectively, and $j_{h}$ is the source term in the hidden sector and arises
from freeze-in. Using the chain rule and together with Eqs. (63) and (64), one
obtains

$\rho\frac{{\rm d}\rho_{h}}{{\rm
d}T_{h}}=\left(\frac{\zeta_{h}}{\zeta}\rho_{h}-\frac{j_{h}}{4H\zeta}\right)\frac{{\rm
d}\rho}{{\rm d}T_{h}}\,,$ (66)

where $\zeta_{h}=\frac{3}{4}(1+p_{h}/\rho_{h})$ and interpolates between
$\zeta_{h}=1$ for radiation dominance and $\zeta_{h}=3/4$ for matter dominance
in the hidden sector. Since $\rho=\rho_{v}+\rho_{h}$, we have ${\rm
d}\rho/{\rm d}T_{h}={\rm d}\rho_{v}/{\rm d}T_{h}+{\rm d}\rho_{h}/{\rm d}T_{h}$
and together with Eq. (66), we can solve for ${\rm d}\rho_{v}/{\rm d}T_{h}$ in
terms of ${\rm d}\rho_{h}/{\rm d}T_{h}$ and get

$\frac{{\rm d}\rho_{v}}{{\rm
d}T_{h}}=\frac{\zeta\rho_{v}+\rho_{h}(\zeta-\zeta_{h})+j_{h}/(4H)}{\zeta_{h}\rho_{h}-j_{h}/(4H)}\frac{{\rm
d}\rho_{h}}{{\rm d}T_{h}}.$ (67)

In the analysis, we use the constraint that the total entropy $S=sR^{3}$ is
conserved which gives ${\rm d}s/{\rm d}t+3Hs=0$ with

$\displaystyle s$ $\displaystyle=\frac{2\pi^{2}}{45}\big{(}h_{{\rm
eff}}^{v}T^{3}+h_{{\rm eff}}^{h}T_{h}^{3}\big{)},$ (68)

where $h_{{\rm eff}}^{v}\,(h_{{\rm eff}}^{h})$ is the visible (hidden)
effective entropy degrees of freedom. The Hubble parameter also depends on
both $T$ and $T_{h}$ as can be seen from the Friedmann equation

$H^{2}=\frac{8\pi G_{N}}{3}\big{[}\rho_{v}(T)+\rho_{h}(T_{h})\big{]},$ (69)

where $\rho_{v}(T)$ and $\rho_{h}(T_{h})$ are the energy density in the
visible sector and in the hidden sector respectively and are given by

$\displaystyle\rho_{v}$ $\displaystyle=\frac{\pi^{2}}{30}g_{{\rm
eff}}^{v}T^{4},\leavevmode\nobreak\ \leavevmode\nobreak\
\rho_{h}=\frac{\pi^{2}}{30}g_{{\rm eff}}^{h}T_{h}^{4}.$ (70)

$g_{{\rm eff}}^{v},h_{{\rm eff}}^{v}$ are functions of $T$ and we use the fits
given in [60, 61, 62] to parameterize them while $g_{{\rm eff}}^{h},h_{{\rm
eff}}^{h}$ are functions of $T_{h}$ and we use temperature dependent integrals
given in [63] to parameterize them.

Combining equations discussed above, the evolution of the ratio of
temperatures for two sectors $\eta=T/T_{h}$ can be derived as

$\frac{{\rm d}\eta}{{\rm
d}T_{h}}=-\frac{A_{v}}{B_{v}}+\frac{\zeta\rho_{v}+\rho_{h}(\zeta-\zeta_{h})+j_{h}/(4H)}{B_{v}[\zeta_{h}\rho_{h}-j_{h}/(4H)]}\frac{{\rm
d}\rho_{h}}{{\rm d}T_{h}}\,,$ (71)

with

$\displaystyle A_{v}$ $\displaystyle=\frac{\pi^{2}}{30}\Big{(}\frac{{\rm
d}g_{{\rm eff}}^{v}}{{\rm d}T}\eta^{5}T_{h}^{4}+4g_{{\rm
eff}}^{v}\eta^{4}T_{h}^{3}\Big{)},$ (72) $\displaystyle B_{v}$
$\displaystyle=\frac{\pi^{2}}{30}\Big{(}\frac{{\rm d}g_{{\rm eff}}^{v}}{{\rm
d}T}\eta^{4}T_{h}^{5}+4g_{{\rm eff}}^{v}\eta^{3}T_{h}^{4}\Big{)}.$ (73)

Now $\rho_{h}$ and $p_{h}$, which enter in the definition of $\zeta_{h}$, are
determined in terms of
$\rho_{\gamma^{\prime}},p_{\gamma^{\prime}},\rho_{\chi},p_{\chi}$ so that
$\rho_{h}=\rho_{\gamma^{\prime}}+\rho_{\chi}$ and
$p_{h}=p_{\gamma^{\prime}}+p_{\chi}$, where $\rho_{\gamma^{\prime}}$ and
$p_{\gamma^{\prime}}$ are given by

$\displaystyle\rho_{\gamma^{\prime}}$
$\displaystyle=\frac{g_{\gamma^{\prime}}T_{h}^{4}}{2\pi^{2}}\int_{x_{\gamma^{\prime}}}^{\infty}\frac{x^{3}{\rm
d}x}{e^{x}-1}\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
p_{\gamma^{\prime}}=\frac{g_{\gamma^{\prime}}T_{h}^{4}}{6\pi^{2}}\int_{x_{\gamma^{\prime}}}^{\infty}\frac{\left(x^{2}-x_{\gamma^{\prime}}^{2}\right){\rm
d}x}{e^{x}-1}\,,$ (74) $\displaystyle\rho_{\chi}$
$\displaystyle=\frac{g_{\chi}T_{h}^{4}}{2\pi^{2}}\int_{x_{\chi}}^{\infty}\frac{x^{3}{\rm
d}x}{e^{x}+1}\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
p_{\chi}=\frac{g_{\chi}T_{h}^{4}}{6\pi^{2}}\int_{x_{\chi}}^{\infty}\frac{\left(x^{2}-x_{\chi}^{2}\right){\rm
d}x}{e^{x}+1}\,,$ (75)

where $g_{\gamma^{\prime}}=3,g_{\chi}=4$ are degrees of freedom of dark photon
and dark fermion respectively, and
$x_{\gamma^{\prime}}=m_{\gamma^{\prime}}/T_{h},x_{\chi}=m_{\chi}/T_{h}$. The
$H$ and $s$ are given by

$\displaystyle s$
$\displaystyle=\frac{2\pi^{2}}{45}\left(h_{\mathrm{eff}}^{v}\eta^{3}+h_{\mathrm{eff}}^{h}\right)T_{h}^{3}\,,$
(76) $\displaystyle H^{2}$
$\displaystyle=\frac{\pi^{2}}{90M_{\mathrm{Pl}}^{2}}\left(g_{\mathrm{eff}}^{v}\eta^{4}+g_{\mathrm{eff}}^{h}\right)T_{h}^{4}\,,$
(77)

and the ratio $s/H$ can be written as

$\displaystyle\frac{s}{H}$
$\displaystyle=\frac{2\sqrt{2}\pi}{\sqrt{45}}\frac{h_{{\rm
eff}}}{\sqrt{g_{{\rm eff}}}}M_{{\rm Pl}}T_{h}\,,$ (78)

where $M_{\mathrm{Pl}}\equiv\sqrt{\frac{1}{8\pi G_{N}}}=2.4\times
10^{18}\leavevmode\nobreak\ \mathrm{GeV}$.

In the dark sector, the effective degrees of freedom include those for the
dark photon and for the dark fermion

$\displaystyle g_{{\rm eff}}^{h}$ $\displaystyle=g_{{\rm
eff}}^{\gamma^{\prime}}+\frac{7}{8}g_{{\rm eff}}^{\chi},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ h_{{\rm eff}}^{h}=h_{{\rm
eff}}^{\gamma^{\prime}}+\frac{7}{8}h_{{\rm eff}}^{\chi}.$ (79)

At temperature $T_{h}$, $g_{{\rm eff}}$ and $h_{{\rm eff}}$ for the particles
$\gamma^{\prime}$ and $\chi$ are given by

$\displaystyle g_{{\rm eff}}^{\gamma^{\prime}}$
$\displaystyle=\frac{45}{\pi^{4}}\int_{x_{\gamma^{\prime}}}^{\infty}\frac{\sqrt{x^{2}-x_{\gamma^{\prime}}^{2}}}{e^{x}-1}x^{2}{\rm
d}x\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
h_{{\rm
eff}}^{\gamma^{\prime}}=\frac{45}{4\pi^{4}}\int_{x_{\gamma^{\prime}}}^{\infty}\frac{\sqrt{x^{2}-x_{\gamma^{\prime}}^{2}}}{e^{x}-1}(4x^{2}-x_{\gamma^{\prime}}^{2}){\rm
d}x\,,$ (80) $\displaystyle g_{{\rm eff}}^{\chi}$
$\displaystyle=\frac{60}{\pi^{4}}\int_{x_{\chi}}^{\infty}\frac{\sqrt{x^{2}-x_{\chi}^{2}}}{e^{x}+1}x^{2}{\rm
d}x\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\,\text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
h_{{\rm
eff}}^{\chi}=\frac{15}{\pi^{4}}\int_{x_{\chi}}^{\infty}\frac{\sqrt{x^{2}-x_{\chi}^{2}}}{e^{x}+1}(4x^{2}-x_{\chi}^{2}){\rm
d}x\,.$

When temperatures go above the mass of the particles, the particles become
relativistic and their effective degrees of freedom take constant values:
$g_{{\rm eff}}^{\gamma^{\prime}}=h_{{\rm eff}}^{\gamma^{\prime}}=3$ for
$T_{h}\gg m_{\gamma^{\prime}}$ and $g_{{\rm eff}}^{\chi}=h_{{\rm
eff}}^{\chi}=4$ for $T_{h}\gg m_{\chi}$.

###  3.2 Coupled Boltzmann equations from two sectors

With the above preparation, we are now ready to write down the Boltzmann
equations for the number densities of the dark fermions $\chi$ and of the dark
photons $\gamma^{\prime}$. For the model we discuss, we have the following
Boltzmann equations for the $\chi$ fermion and dark photon $\gamma^{\prime}$
number densities

$\displaystyle\frac{{\rm d}n_{\chi}}{{\rm d}t}+3Hn_{\chi}=\sum_{i\in{\rm SM}}$
$\displaystyle\Big{\\{}{\color[rgb]{0,1,0}(n_{\chi}^{{\rm
eq}})^{2}\langle\sigma v\rangle_{\chi\bar{\chi}\to
i\bar{i}}^{T}+n_{\gamma^{*}}\langle\Gamma\rangle_{\gamma^{*}\rightarrow\chi\bar{\chi}}^{T}}$
$\displaystyle+\theta(M_{\gamma^{\prime}}-2m_{\chi})\big{[}-n_{\chi}^{2}\langle\sigma
v\rangle_{\bar{\chi}\chi\rightarrow\gamma^{\prime}}^{T_{h}}+n_{\gamma^{\prime}}\langle\Gamma\rangle_{\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}\big{]}$
$\displaystyle-n_{\chi}^{2}\langle\sigma
v\rangle_{\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}}^{T_{h}}+n_{\gamma^{\prime}}^{2}\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}\Big{\\}}\,,$
(81) $\displaystyle\frac{{\rm d}n_{\gamma^{\prime}}}{{\rm
d}t}+3Hn_{\gamma^{\prime}}=\sum_{i\in{\rm SM}}$
$\displaystyle\Big{\\{}n_{\chi}^{2}\langle\sigma
v\rangle_{\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}}^{T_{h}}-n_{\gamma^{\prime}}^{2}\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}$
$\displaystyle+\theta(M_{\gamma^{\prime}}-2m_{\chi})\big{[}n_{\chi}^{2}\langle\sigma
v\rangle_{\bar{\chi}\chi\rightarrow\gamma^{\prime}}^{T_{h}}-n_{\gamma^{\prime}}\langle\Gamma\rangle_{\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}\big{]}$
$\displaystyle+{\color[rgb]{0,1,0}\theta(M_{\gamma^{\prime}}-2m_{i})\big{[}n_{i}^{2}\langle\sigma
v\rangle_{i\bar{i}\rightarrow\gamma^{\prime}}^{T}}-n_{\gamma^{\prime}}\langle\Gamma\rangle_{\gamma^{\prime}\to
i\bar{i}}^{T_{h}}\big{]}$
$\displaystyle+{\color[rgb]{0,1,0}n_{i}^{2}\langle\sigma
v\rangle_{i\bar{i}\rightarrow\gamma^{\prime}\gamma}^{T}+2n_{i}n_{\gamma^{\prime}}^{{\rm
eq}}\langle\sigma v\rangle_{i\gamma^{\prime}\rightarrow
i\gamma}^{T}}\Big{\\}}\,,$ (82)

where the upper subscript is the visible sector temperature $T$ or hidden
sector temperature $T_{h}$ that the corresponding process is depending on. The
freeze-in processes are marked by  green color to be distinguishable to other
hidden sector processes. Later we will use $T_{h}$ as the reference
temperature and replace $T$ by $T_{h}\eta$. In the above equations, $i$ and
$\bar{i}$ represent SM fermions and anti-fermions, $\gamma^{*}$ denotes the
plasmon which expresses photons with effective in-medium mass in the early
universe, c.f., Appendix B for a brief review. All the terms in the above
equations with subscript eq indicate that we are calculating the corresponding
process reversely for computation simplicity. In such calculations, we use the
equilibrium values of comoving number densities (depending on the visible
sector temperature $T$)

$Y_{a}^{{\rm
eq}}=\frac{n_{a}^{\mathrm{eq}}}{s}=\frac{g_{a}}{2\pi^{2}s}m_{a}^{2}TK_{2}\left(\frac{m_{a}}{T}\right)\,,\qquad{\rm
for\,particle\,species}\,a\,,$ (83)

which are usually not in equilibrium in the thermal bath. For example, the
freeze-in process $i\bar{i}\to\chi\bar{\chi}$ is computed via the inverse
process $\chi\bar{\chi}\to i\bar{i}$ using the comoving number density of
$\chi$ as if it is in equilibrium in the thermal bath, which is much easier to
calculate compared to the original process. The theta functions indicate
whether the three-point decay or combination can occur. If the mass of dark
photon is greater than twice of the dark fermion mass, i.e.,
$M_{\gamma^{\prime}}>2m_{\chi}$, three-point processes
$\gamma^{\prime}\leftrightarrow\chi\bar{\chi}$ are much stronger than the
four-point processes
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ which can be
safely ignored. However, we can conclude from our analysis in Section 4.3: for
the freeze-in production of the dark photon, the four-point freeze-in channels
must be considered even the three-point freeze-in channels are present. The
additional factor of 2 corresponding to the process $i\gamma\rightarrow
i\gamma^{\prime}$ also counts the anti-fermion contribution
$\bar{i}\gamma\rightarrow\bar{i}\gamma^{\prime}$.

All relevant freeze-in processes for the dark fermion $\chi$ are shown in Fig.
3. The freeze-in processes for the dark photon $\gamma^{\prime}$ are shown in
Fig. 4, if $M_{\gamma^{\prime}}$ is greater than twice of electron so that the
three-point processes via SM fermion combination $i\bar{i}\to\gamma^{\prime}$
can occur. If $M_{\gamma^{\prime}}<2m_{e}$, since the partial decay widths of
dark photon to neutrinos are highly suppressed, c.f., Eq. (45), the feasible
production channels can only be the four-point interactions, shown in Fig. 5.
We notice here that the production channel of
$i\bar{i}\to\gamma^{\prime}\gamma^{\prime}$ are suppressed, compared to
$i\bar{i}\to\gamma\gamma^{\prime}$. Four-point freeze-in production via
interactions $i\bar{i}\to Z\gamma^{\prime}$, $iZ\rightarrow i\gamma^{\prime}$
and $\bar{i}Z\rightarrow\bar{i}\gamma^{\prime}$ are also suppressed compared
with the same interactions with $Z$ replaced by the photon. We also collect
all the dark sector self-interactions in Fig. 6, depending on the masses of
the dark fermion and the dark photon.

Figure 3:  A display of dark fermion production channels, including the SM
fermion annihilation through $s$ and $t$ channels, as well as the plasmon
decay.  Figure 4:  A display of dark photon three-point production channel,
when the dark photon mass is heavier than at least twice of electrons. As the
mass of the dark photon increases, more combination channels involving SM
fermions will become accessible. Figure 5: A display of dark photon four-point
production channels, including $i\bar{i}\to\gamma^{\prime}+\gamma,Z$ and
$i(\bar{i})+\gamma,Z\to i(\bar{i})\gamma^{\prime}$, although interactions
involving the $Z$ boson are suppressed. Crucially, we noticed in our analysis
that even the three-point productions shown in Fig. 4 are available, the
contribution from four-point production channels are still quite significant
and should not be neglected. Figure 6:  A display of various hidden sector
self-interactions with particle number changing, including:
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}$,
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$. The hidden
sector will encounter dark matter freeze-out depending on the masses of the
dark particles, discussed explicitly in Section 4 for six different cases.

We then analyze the evolution of $n_{\chi}$, $n_{\gamma^{\prime}}$ and $\eta$
as a function of $T_{h}$. For the computation of the relic density, it is more
convenient to deal directly with the comoving number densities defined by the
number density of a species $a$ divided by the entropy $s$, i.e.,
$Y_{a}=n_{a}/s$. We assume there is no initial abundance for the dark
particles $\chi,\gamma^{\prime}$ and they are produced through freeze-in
processes from SM particles. The hidden sector particles can have rather
strong interactions such as
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ within the
hidden sector. Thus for the case that $U(1)_{X}$ hidden sector interacts with
the visible sector with feeble couplings and evolves almost independently with
a temperature $T_{h}=T/\eta$, the evolution equations of the dark fermion
$\chi$ and the dark photon $\gamma^{\prime}$ are given by the coupled
Boltzmann equations for the comoving number densities $Y_{\chi}$ and
$Y_{\gamma^{\prime}}$ and $\eta$:

$\displaystyle\frac{{\rm d}Y_{\chi}}{{\rm d}T_{h}}=$
$\displaystyle-\frac{s}{H}\frac{{\rm d}\rho_{h}/{\rm
d}T_{h}}{4\rho_{h}-j_{h}/H}\sum_{i\in{\rm
SM}}\Big{\\{}{\color[rgb]{0,1,0}{(Y_{\chi}^{{\rm eq}})^{2}\langle\sigma
v\rangle_{\chi\bar{\chi}\to
i\bar{i}}^{T_{h}\eta}+\frac{1}{s}Y_{\gamma^{*}}\langle\Gamma\rangle_{\gamma^{*}\rightarrow\chi\bar{\chi}}^{T_{h}\eta}}}$
$\displaystyle+\theta(M_{\gamma^{\prime}}-2m_{\chi})\big{[}-Y_{\chi}^{2}\langle\sigma
v\rangle_{\bar{\chi}\chi\rightarrow\gamma^{\prime}}^{T_{h}}+\frac{1}{s}Y_{\gamma^{\prime}}\langle\Gamma\rangle_{\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}\big{]}$
$\displaystyle-Y_{\chi}^{2}\langle\sigma
v\rangle_{\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}}^{T_{h}}+Y_{\gamma^{\prime}}^{2}\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}\Big{\\}}\,,$
(84) $\displaystyle\frac{{\rm d}Y_{\gamma^{\prime}}}{{\rm d}T_{h}}=$
$\displaystyle-\frac{s}{H}\frac{{\rm d}\rho_{h}/{\rm
d}T_{h}}{4\rho_{h}-j_{h}/H}\sum_{i\in{\rm
SM}}\Big{\\{}Y_{\chi}^{2}\langle\sigma
v\rangle_{\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}}^{T_{h}}-Y_{\gamma^{\prime}}^{2}\langle\sigma
v\rangle_{\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}$
$\displaystyle+\theta(M_{\gamma^{\prime}}-2m_{\chi})\big{[}Y_{\chi}^{2}\langle\sigma
v\rangle_{\bar{\chi}\chi\rightarrow\gamma^{\prime}}^{T_{h}}-\frac{1}{s}Y_{\gamma^{\prime}}\langle\Gamma\rangle_{\gamma^{\prime}\to\chi\bar{\chi}}^{T_{h}}\big{]}$
$\displaystyle+{\color[rgb]{0,1,0}\theta(M_{\gamma^{\prime}}-2m_{i})\big{[}Y_{i}^{2}\langle\sigma
v\rangle_{i\bar{i}\rightarrow\gamma^{\prime}}^{T_{h}\eta}}-\frac{1}{s}Y_{\gamma^{\prime}}\langle\Gamma\rangle_{\gamma^{\prime}\to
i\bar{i}}^{T_{h}}\big{]}$
$\displaystyle+{\color[rgb]{0,1,0}Y_{i}^{2}\langle\sigma
v\rangle_{i\bar{i}\rightarrow\gamma^{\prime}\gamma}^{T_{h}\eta}+2Y_{i}Y_{\gamma^{\prime}}^{{\rm
eq}}\langle\sigma v\rangle_{i\gamma^{\prime}\rightarrow
i\gamma}^{T_{h}\eta}}\Big{\\}}\,,$ (85) $\displaystyle\frac{{\rm d}\eta}{{\rm
d}T_{h}}=$
$\displaystyle-\frac{A_{v}}{B_{v}}+\frac{\rho_{v}+j_{h}/(4H)}{\rho_{h}-j_{h}/(4H)}\frac{\frac{{\rm
d}\rho_{h}}{{\rm d}T_{h}}}{B_{v}}\,,$ (86)

where the energy transfer density (the source term of the hidden sector) is
given by

$\displaystyle j_{h}=$ $\displaystyle\sum_{i\in{\rm
SM}}\Big{[}2Y_{i}^{2}s^{2}J_{i\bar{i}\rightarrow\chi\bar{\chi}}^{T}+2Y_{\gamma^{*}}sJ_{\gamma^{*}\rightarrow\chi\bar{\chi}}^{T}+Y_{i}^{2}s^{2}J_{i\bar{i}\rightarrow\gamma^{\prime}\gamma}^{T}+2Y_{i}Y_{\gamma^{\prime}}^{{\rm
eq}}s^{2}J_{i\gamma^{\prime}\rightarrow i\gamma}^{T}$
$\displaystyle+\theta(M_{\gamma^{\prime}}-2m_{i})\big{(}Y_{i}^{2}s^{2}J_{i\bar{i}\rightarrow\gamma^{\prime}}^{T}-Y_{\gamma^{\prime}}sJ_{\gamma^{\prime}\rightarrow
i\bar{i}}^{T_{h}}\big{)}\Big{]}\,,$ (87)

and

$\displaystyle
Y_{a}^{\mathrm{eq}}(T)Y_{b}^{\mathrm{eq}}(T)s^{2}J(ab\rightarrow cd)(T)$
$\displaystyle=\frac{T}{32\pi^{4}}\int_{s_{0}}^{\infty}ds\sigma_{cd\rightarrow
ab}s\left(s-s_{0}\right)K_{2}(\sqrt{s}/T)\,,$ (88) $\displaystyle
Y_{a}J_{a\rightarrow bc}$ $\displaystyle=Y_{a}m_{a}\Gamma_{a\rightarrow
bc}\,.$ (89)

The cross-section is given by

$\displaystyle\sigma_{ab\to cd}$
$\displaystyle=\frac{1}{g_{a}g_{b}}\frac{1}{2E_{a}2E_{b}\left|v_{a}-v_{b}\right|}\int\mathrm{d}\Pi_{c}\mathrm{d}\Pi_{d}\left|M\right|_{ab\to
cd}^{2}\left(2\pi\right)^{4}\delta^{4}\left(p_{a}+p_{b}-p_{c}-p_{d}\right)\,,$
(90)

where $\left|M\right|_{ab\to cd}^{2}$ is summed over initial and final spins,
the momentum phase space element is defined as

$\mathrm{d}\Pi_{i}\equiv\frac{\mathrm{d}^{3}\mathbf{p}_{i}}{\left(2\pi\right)^{3}2E_{i}}\,,$
(91)

and the thermally averaged cross-section reads

$\displaystyle\left\langle\sigma v\right\rangle_{ab\to cd}^{T}\equiv$
$\displaystyle\frac{\int\mathrm{d}\Pi_{a}\mathrm{d}\Pi_{b}\,\sigma
ve^{-E_{a}/T}e^{-E_{b}/T}}{\int\mathrm{d}\Pi_{a}\mathrm{d}\Pi_{b}e^{-E_{a}/T}e^{-E_{b}/T}}\,.$
(92)

We now inspect some specific details concerning coefficients that require
attention when solving the Boltzmann equations.

  * •

For processes involving SM fermions, one needs to consider all types of
fermions along with their corresponding anti-fermions. Take the process
$i\bar{i}\to\gamma^{\prime}$ as an example, in the calculation we need the
value of $n_{q_{i}}n_{\bar{q}_{i}}$ which is derived as

$\displaystyle
n_{q_{i}}n_{\bar{q}_{i}}=3g_{q}g_{\bar{q}}\left[\frac{m_{i}^{2}}{2\pi^{2}}TK_{2}\left(\frac{m_{i}}{T}\right)\right]^{2}=\left(3\times
2\times
2\right)\times\left[\frac{m_{i}^{2}}{2\pi^{2}}TK_{2}\left(\frac{m_{i}}{T}\right)\right]^{2}\,,$
(93)

where the factor 3 indicates three colors of quarks annihilate with anti-
quarks carrying the same color, for all six flavors of quarks. This
calculation is analogous for leptons without accounting for the color factor.

  * •

When calculating the energy density and pressure, contributions from both
$\chi$ and $\bar{\chi}$ should be considered. Hence in Eq. (75) we have

$g_{\text{$\chi$}}^{\mathrm{total}}=g_{\chi}+g_{\bar{\chi}}=4\,.$ (94)

  * •

Since we consider a model with symmetric component of dark matter, we have
$Y_{\chi}=Y_{\bar{\chi}}$ with $g_{\chi}=g_{\bar{\chi}}=2$. The evolution of
$\bar{\chi}$ is identical to the evolution of $\chi$ in the hidden sector, and
thus the total dark fermion number density in a comoving volume is
$2Y_{\chi}$.

All types of thermally averaged cross-sections for various processes are
summarized below:

  1. 1.

Fermion annihilation $a\bar{a}\rightarrow bc$. SM fermions or hidden sector
fermions can annihilate to a pair of fermion and anti-fermion, or to two gauge
bosons, which can be universally expressed as

$\langle\sigma v\rangle_{a\bar{a}\rightarrow
bc}^{T}=\frac{1}{8m_{a}^{4}TK_{2}^{2}\left(m_{a}/T\right)}\int_{4m_{a}^{2}}^{\infty}ds\sigma(s)\sqrt{s}\left(s-4m_{a}^{2}\right)K_{1}(\sqrt{s}/T)\,.$
(95)

  2. 2.

Fermion elastic scattering $ab\rightarrow cd$. A fermion scatters off a gauge
boson, producing a fermion with a gauge boson

$\displaystyle\langle\sigma v\rangle_{ab\rightarrow cd}^{T}$
$\displaystyle=\frac{1}{8m_{a}^{2}m_{b}^{2}TK_{2}\left(m_{a}/T\right)K_{2}\left(m_{b}/T\right)}\times$
(96)
$\displaystyle\quad\quad\quad\int_{\left(m_{a}+m_{b}\right)^{2}}^{\infty}ds\sigma(s)\sqrt{s}\left[s-\left(m_{a}+m_{b}\right)^{2}\right]K_{1}\left(\sqrt{s}/T\right)\,.$

  3. 3.

Neutral gauge boson decay $a\rightarrow b\bar{b}$. A neutral gauge boson
($Z,\gamma^{\prime}$ and plasmon $\gamma^{*}$) can decay into a pair of
fermion and anti-fermion

$\left\langle\Gamma\right\rangle_{a\rightarrow
b\bar{b}}^{T}=\Gamma_{a\rightarrow
b\bar{b}}\frac{K_{1}\left(m_{a}/T\right)}{K_{2}\left(m_{a}/T\right)}\,.$ (97)

  4. 4.

Fermion pair combination $b\bar{b}\rightarrow a$. A pair of fermion and anti-
fermion combines into a massive gauge boson

$\langle\sigma v\rangle_{b\bar{b}\rightarrow
a}^{T}=\frac{\sigma^{0}_{b\bar{b}\to
a}(m_{a}^{2})m_{a}\left(m_{a}^{2}-4m_{b}^{2}\right)K_{1}\left(m_{a}/T\right)}{8m_{b}^{4}TK_{2}^{2}\left(m_{b}/T\right)}\,,$
(98)

with $\sigma^{0}_{b\bar{b}\to a}$ representing the part without the delta
function in the cross-section $\sigma_{b\bar{b}\rightarrow a}(s)$.

In the above equations, $K_{1}$ and $K_{2}$ are the modified Bessel function
of the second kind and degrees one and two, respectively. The cross-sections
of all relevant processes discussed in this work are given in Appendix C.

Given a set of values
$\\{m_{\chi},M_{\gamma^{\prime}},g_{X},\delta,\epsilon\\}$ and the initial
value of $\eta=T/T_{h}$, one can calculate the full evolution of the hidden
sector particles as well as the change of the hidden sector temperature. The
initial value of $\eta=T/T_{h}$ will not modify the final abundance so much
but will alter the details of the dark particles evolution, as discussed in
[37, 64]. A complete analysis for the evolution of a $U(1)$ hidden sector
interacting feebly with the SM from the reheating stage, including the
determination of the initial value of $\eta$, is discussed in a companion
paper [65].

In the calculation, the expansion of non-relativistic thermally averaged
cross-section is useful, and is given in [61] as

$\langle\sigma
v\rangle=a^{(0)}+\frac{3}{2}a^{(1)}x^{-1}+\frac{15}{8}a^{(2)}x^{-2}+\cdots$
(99)

at low temperatures, where $x=m/T$, and $a^{(n)}$ indicates the $n$th
derivative of $\sigma v_{\rm lab}$ with respect to $\kappa=(s-4m^{2})/4m^{2}$
($m$ is the mass of the incoming particle) evaluated at $\kappa=0$. Similarly
for the thermally averaged decay width, one has the expansion at low
temperatures,

$\langle\Gamma\rangle_{a\rightarrow bc}=\Gamma_{a\rightarrow
bc}\left(1-\frac{3}{2}x^{-1}+\frac{15}{8}x^{-2}+\cdots\right)\,.$ (100)

##  4 Evolution of matter in visible and hidden sectors

###  4.1 $U(1)_{X}$ sub-GeV dark matter via freeze-in

In this section, we discuss the production of $U(1)_{X}$ sub-GeV millicharge
dark matter $\chi$ and dark photon $\gamma^{\prime}$ via freeze-in mechanism.
In some of the cases, the dark photon can also be stable and serve as a dark
matter candidate, although we will see its abundance is usually tiny. As
discussed in Section 2, millicharged dark matter can be produced via
Stueckelberg mass mixing, with a millicharge proportional to
$\epsilon=M_{2}/M_{1}$. In addition to the freeze-in production processes,
there are also quite strong interactions within the hidden sector like
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$, and
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}$ if the dark photon is heavier
than twice of $\chi$ mass. Since the couplings between the visible sector and
the hidden sector are feeble, the hidden sector temperature $T_{h}$ evolves
almost independently from the visible sector temperature $T$ (the temperature
of the Universe). As the freeze-in processes and the hidden sector
interactions (depending on the hidden sector temperature $T_{h}$) occur
simultaneously, we use the formalism discussed in Section 3 to calculate the
full evolution of the hidden sector particles as well as the hidden sector
temperature.

We will now explore six distinct scenarios upon the masses of dark matter
$\chi$ and the dark photon $\gamma^{\prime}$, each exhibiting markedly unique
characteristics in the evolution of the hidden sector. The common feature of
all these scenarios is that they have the same freeze-in production channel:

  * •

Dark fermion $\chi$ freeze-in production channels: SM fermion annihilation
$i\bar{i}\to\chi\bar{\chi}$, plasmon decay $\gamma^{*}\to\chi\bar{\chi}$. The
$Z$ boson decay channel $Z\to\chi\bar{\chi}$ is suppressed, since the
abundance of $Z$ in the early universe is Boltzmann suppressed at the
temperature when the dark fermion $\chi$ is majorly produced (around the mass
of $\chi$ in the sub-GeV region).

  * •

Dark photon $\gamma^{\prime}$ freeze-in production channels: For
$M_{\gamma^{\prime}}>2m_{e}$, the major production channel is the three-point
channels $i\bar{i}\to\gamma^{\prime}$ and four-point channels
$i\bar{i}\to\gamma\gamma^{\prime}$, $i\gamma\to i\gamma^{\prime}$,
$\bar{i}\gamma\to\bar{i}\gamma^{\prime}$. For $M_{\gamma^{\prime}}<2m_{e}$,
the production channels are $i\bar{i}\to\gamma\gamma^{\prime}$, $i\gamma\to
i\gamma^{\prime}$, $\bar{i}\gamma\to\bar{i}\gamma^{\prime}$. We emphasize here
again that in the model the couplings of the dark photon to neutrinos are
highly suppressed, c.f., Eq. (45), and thus either the decay of dark photon to
neutrinos or the freeze-in production via neutrinos combination can be safely
neglected.

We firstly discuss the case that the dark photon mass is greater than twice of
the electron mass, i.e., $M_{\gamma^{\prime}}>2m_{e}$. In this case, the dark
fermion $\chi$ is millicharged and serves as the dark matter candidate. The
dark photons, due to the mixing of $U(1)_{X}$ and the SM, will ultimately
undergo decay into electron and positron pairs, and are thus not dark matter
candidates. The dark photon also has to decay before BBN for a valid model.
Within this category we discuss three different cases:

##### Case 1: $m_{\chi}>M_{\gamma^{\prime}}>2m_{e}$, dark matter $\chi$

In this case hidden sector particles can interchange through the interaction
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ (dark matter
self-interactions such as $\chi\bar{\chi}\leftrightarrow\chi\bar{\chi}$ and
$\chi\chi\leftrightarrow\chi\chi$ will not change the number densities of dark
fermion and dark photon and we will not discuss at this point). When the
temperature drops down, the hidden sector encounters the dark fermion $\chi$
freeze-out. The dark photon will eventually decay into SM fermions through
$\gamma^{\prime}\to i\bar{i}$ (mostly to electron and positron pairs).

##### Case 2: $2m_{\chi}>M_{\gamma^{\prime}}>m_{\chi}>2m_{e}$, dark matter
$\chi$

In this case hidden sector particles interchange through the interaction
$\gamma^{\prime}\gamma^{\prime}\leftrightarrow\chi\bar{\chi}$. When the
temperature drops down, the hidden sector encounters the dark photon
$\gamma^{\prime}$ freeze-out. The leftover of the dark photon finally decay to
SM fermions through $\gamma^{\prime}\to i\bar{i}$.

##### Case 3: $M_{\gamma^{\prime}}>2m_{\chi}>2m_{e}$, dark matter $\chi$

In this case hidden sector particles interchange mostly through the three-
point interaction $\gamma^{\prime}\leftrightarrow\chi\bar{\chi}$. Because of
the rather strong interactions among the hidden sector, the dark photon will
soon decay to $\chi\bar{\chi}$ when the temperature drops down. Four-point
interactions can be neglected in this case. In addition, the decay channel of
dark photon to SM fermions can also be neglected since the decay of
$\gamma^{\prime}\to\chi\bar{\chi}$ is so much strong.

We then discuss the case that the dark photon mass is less than twice of the
electron mass, i.e., $M_{\gamma^{\prime}}<2m_{e}$. Again the dark fermion
$\chi$ is millicharged and serves as the dark matter candidate. Since the
decay channel $\gamma^{\prime}\to e\bar{e}$ is forbidden, the dark photon can
only decay to two neutrinos and to three photons which are usually quite rare,
c.f., Fig 1. Focusing on these two interactions, by tuning the $U(1)_{X}$
gauge coupling as well as the kinetic mixing parameter $\delta$ and mass
mixing parameter $\epsilon$, one can make the photon stable throughout the
entire age of the Universe, establishing the dark photon as a viable dark
matter candidate. However, although the dark photon’s lifetime is extended
beyond the age of the Universe, it can still undergo decay, even in minuscule
amounts. This decay contributes to the isotropic diffuse photon background
(IDPB), which can be precisely measured. The IDPB imposes an even more
stringent constraint on the couplings between the dark photon and SM fermions.
We incorporate this constraint into our analysis and discover that it further
diminishes the $\delta,\epsilon$ when compared to the constraints posed by
dark photon decay.

Within the category $M_{\gamma^{\prime}}<2m_{e}$ we also discuss three
different cases:666For the models with $M_{\gamma^{\prime}}<2m_{e}$, one
cannot make the dark photon unstable and decaying before BBN considering all
the relevant experimental constraints.

##### Case 4: $2m_{e}>m_{\chi}>M_{\gamma^{\prime}}$, dark matter
$\chi,\gamma^{\prime}$

The hidden sector evolves through the interaction
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$, and the freeze-
out of $\chi$ will occur when the temperature drops. The dark photon
$\gamma^{\prime}$ is stable and receives stringent constraints from both its
decay width as well as the IDPB.

##### Case 5: $2m_{e}>2m_{\chi}>M_{\gamma^{\prime}}>m_{\chi}$, dark matter
$\chi,\gamma^{\prime}$

The hidden sector evolves through the interaction
$\gamma^{\prime}\gamma^{\prime}\leftrightarrow\chi\bar{\chi}$, and the freeze-
out of $\gamma^{\prime}$ will occur when the temperature drops. The dark
photon $\gamma^{\prime}$ is stable so the leftover of $\gamma^{\prime}$ cannot
decay. The dark photon $\gamma^{\prime}$ again receives stringent constraints
from both its decay width as well as the IDPB.

##### Case 6: $2m_{e}>M_{\gamma^{\prime}}>2m_{\chi}$, dark matter $\chi$

For this case the hidden sector evolves via
$\gamma^{\prime}\leftrightarrow\chi\bar{\chi}$. The decay
$\gamma^{\prime}\to\chi\bar{\chi}$ is so strong given a relatively high
$U(1)_{X}$ gauge coupling $g_{X}$, when temperature drops all the dark photon
will decay into dark fermions. The evolution of the dark sector in this case
is almost identical to Case 3.

Upon all the above six distinct cases, the interactions among the dark
particles play very important roles in the hidden sector evolution. As the
dark particles $\chi$ and $\gamma^{\prime}$ are accumulating via the freeze-in
production, at some moment the hidden sector particles may reach thermal
equilibrium within the hidden sector, since the hidden sector (gauge) coupling
constant is much stronger than the freeze-in couplings. Depending on the
masses of $\chi$ and $\gamma^{\prime}$, the three-point interaction
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}$ and the four-point interactions
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ become active.
For Case 3 and Case 6, the three-point interaction
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}$ is allowed, which overwhelms
the four-point interaction
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$. Eventually, the
dark photon will entirely decay to $\chi\bar{\chi}$, leaving the dark fermion
as the dark matter. For Case 1 and Case 4, in which
$m_{\chi}>M_{\gamma^{\prime}}$, the dark freeze-out for $\chi$ will occur in
the hidden sector. As the interaction
$\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}$ goes out of equilibrium, the
dark freeze-out ended and the amount of the dark fermion $\chi$ freezes at the
final value. While for Case 2 and Case 5, in which
$2m_{\chi}>M_{\gamma^{\prime}}>m_{\chi}$, the dark freeze-out for
$\gamma^{\prime}$ will occur in the hidden sector. The dark freeze-out will
finish when the interaction $\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$
goes out of equilibrium, the leftover of the dark photon will either decay to
electron-positron pair if $M_{\gamma^{\prime}}>2m_{e}$ (Case 2) or will remain
as another component of dark matter if $M_{\gamma^{\prime}}<2m_{e}$ (Case 5)
and its lifetime exceeds the age of the Universe.

###  4.2 Constraints, benchmark models and phenomenology

First of all we discuss all the relevant experimental constraints for the
models we consider in this work. In the literature, people usually assume a
kinetic mixing between an additional $U(1)_{X}$ and the $U(1)_{\rm em}$ with a
mixing parameter $\bm{\delta}$, which is different from the kinetic mixing
parameter $\delta$ in our work indicating the kinetic mixing between
$U(1)_{X}$ and the hypercharge gauge field $U(1)_{Y}$. As experiments setting
constraints on the kinetic mixing parameter mostly from the dark photon
interaction with electron and positron pair, we will have a translated
relation of the current experimental constraints on our kinetic mixing
parameter $\delta$. The coupling of the dark photon to electron and positron
pair through the Lagrangian from Eqs. (7) and (8) is computed to be

$\displaystyle\mathcal{L}_{\gamma^{\prime}e\bar{e}}$
$\displaystyle=-\varepsilon
eJ^{\mu}A_{\mu}^{\prime}=-\frac{g_{2}g_{Y}}{\sqrt{g_{2}^{2}+g_{Y}^{2}}}\bm{\delta}\bar{e}\gamma^{\mu}eA_{\mu}^{\prime}\,.$
(101)

While for the case that $U(1)_{X}$ mixed with $U(1)_{Y}$ via kinetic terms and
the Stueckelberg mass mixing, the coupling of the dark photon to electron and
positron pair is given by Eq. (39)

$\mathcal{L}_{\gamma^{\prime}e\bar{e}}=-\frac{1}{2}\bar{e}\gamma^{\mu}\left(v_{e}^{\prime}-a_{e}^{\prime}\gamma^{5}\right)eA_{\mu}^{\prime}=-\frac{g_{2}^{2}g_{Y}^{\prime}}{g_{2}^{2}+g_{Y}^{\prime
2}}\left(\delta-\epsilon\right)\bar{e}\gamma^{\mu}eA_{\mu}^{\prime}\,.$ (102)

Comparing Eqs. (101) and (102), one finds the relation

$\displaystyle\left(\delta-\epsilon\right)=\frac{g_{Y}\left(g_{2}^{2}+g_{Y}^{\prime
2}\right)}{g_{2}g_{Y}^{\prime}\sqrt{g_{2}^{2}+g_{Y}^{2}}}\bm{\delta}\approx\frac{\sqrt{g_{2}^{2}+g_{Y}^{2}}}{g_{2}}\bm{\delta}\,,$
(103)

where $g_{Y}^{\prime}\approx g_{Y}$ for small $\delta$, $\epsilon$. Thus the
experimental constraint set on $\bm{\delta}$ (parameter for kinetic mixing
between $U(1)_{X}$ and the $U(1)_{\rm em}$ in the simplified model discussed
in Section 2.1) can be translated to a constraint on the parameter
$|\delta-\epsilon|$ for the model we consider in this work, discussed in
Section 2.4.1. We plot the current experimental constraints on the kinetic
mixing parameter $\delta$ minus mass mixing parameter $\epsilon$ in Fig. 7,
with all our benchmark models exhibited.

Figure 7: [color online] A display of current constraints (colored regions) on
the absolute value of the kinetic mixing parameter minus the mass mixing
parameter $|\delta-\epsilon|$, translated from the current constraints on the
simplified kinetic mixing model people usually referenced, reviewed in Section
2.1, using Eq. (103). The limits for $M_{\gamma^{\prime}}\gtrsim$ 1 MeV come
from supernova (SN1987A) [66], $(g-2)_{e}$ [67] and from various beam dump
experiments: E141 [68], E137 [69, 70, 71], $\nu$-Cal [72, 73]. The light
orange area in the upper right labeled “Others” includes some collider/fixed
target experiments: CMS [74], BaBar [75], LHCb [76], KLOE [77, 78, 79, 80] and
A1 [81]. For the dark photon mass in the 1 keV – 1 MeV window, the shaded
regions depict the existing constraints from solar lifetime (Sun-T), red
giants (RG), horizontal branches (HB) [83, 82, 84] and from isotropic diffuse
photon background (IDPB) [85, 86, 87]. All our benchmark models are also shown
in the plot.

As discussed in Section 2.4.1, when the dark photon is less massive than twice
of electron mass, dark photons can only decay to three photons, and to pairs
of neutrino and anti-neutrino, as shown in Fig. 1. Within the parameter values
we employ for our benchmark models, the decay of the dark photon to three
photons is dominant over the neutrino decay channel.

Photons from the decay of the dark photon (even its lifetime exceeds the age
of the Universe) will affect the observed $\gamma$-ray background and
consequently impose stringent constraints on the parameter space. Various
bounds on the decaying dark matter with a wide range of mass have been
explored in details in [88, 89, 90, 91, 92, 93], including 2-body and 3-body
decays. Specifically, constraints on dark photon dark matter in the keV-MeV
mass range were explored in [86], requiring that the $\gamma$-ray flux from
decaying dark photon does not exceed the total observed flux of $\gamma$-ray.
Model-independent constraints were established in [85] by demanding that the
photon flux from $\gamma^{\prime}\to 3\gamma$ decays must not surpass the
measurements of the diffuse $\gamma$-ray background:

$\frac{\mu\,\tau}{\mathrm{GeV\,s}}\lesssim 10^{27}\left(\frac{\omega}{\rm
GeV}\right)^{1.3}\left(\frac{\Omega h^{2}}{0.1}\right)\,,$ (104)

with $\mu,\,\tau,\,\omega$ representing mass, lifetime and the energy of the
decaying dark photon respectively. The lifetime $\tau$ can be calculated using
the decay width $\Gamma_{\gamma^{\prime}\to 3\gamma}$ given by Eq. (47), and
the corresponding constraint is plotted in Fig. 7 indicated by “IDPB”.

Figure 8: [color online] A display of current constraints on the millicharge
carried by the dark particle versus its mass. The constraints arise from: red
giants and white dwarfs (RGWD) [9], supernova (SN1897A) [16], SLAC milliQ
experiment [6], $N_{\rm eff}$ during BBN and CMB [9, 13], overclosure of the
Universe [4, 7] and multifarious collider experiments such as LEP, LHC, etc
[1, 7]. We also show the limit on the elastic scattering between baryons and
millicharged dark matter given by [15]. Our benchmark models are also shown in
the plot.

People also searched millicharge particles as dark matter candidates and
imposed constraints from laboratory, astrophysical and cosmological
experiments for a wide mass range of the millicharge particle. Bounds from
colliders including LEP, LHC, etc [1, 7], cover a wide mass range but are less
restrictive compared to other constraints. Astrophysical constraints involving
red giants, horizontal-branch stars, white dwarfs and supernova 1987A, are
caused by new energy-loss channels which will lead to observational
alterations to the standard stellar evolution process [3, 4, 5, 1, 94, 9, 16].
Additionally, particles carrying tiny electric charges have the possibility to
interact with the plasma and get thermally excited in the early universe. Both
dark photons and millicharge particles in the hidden sector may contribute to
$N_{\rm eff}$, leading to deviations from the standard value, and will thus
receive constraints from cosmological observations [95, 4, 96, 9, 13]. Another
important cosmological bound is from the overclosure of the Universe,
demanding the relic density of the millicharge dark matter to be less than the
critical density today $\left(\Omega<1\right)$ [4, 7]. Limits on the elastic
scattering cross-section between sub-GeV millicharge dark matter and baryons
are given by [15], following the treatment of previous works [97, 10, 11, 99,
98, 14] and using CMB temperature and polarization data [12] together with
Lyman-$\alpha$ flux power spectrum data [100]. The most relevant experimental
constraints on the millicharge carried by the dark particles are plotted in
Fig. 8, with all our benchmark models exhibited in the allowed region.

We mention in passing that the (stable) dark photon will not contribute
excessively to effective relativistic degree of freedom in the Universe in
Case 4 and Case 5 of our benchmark models. The dark photon contribution to
$\Delta N_{\rm eff}$ is given by

$\Delta N_{\rm
eff}=\frac{8}{7}g_{\gamma^{\prime}}\left(\frac{T_{h}}{T_{\gamma}}\right)^{4}\,.$
(105)

At low temperatures, $g_{\gamma^{\prime}}=3$ and ${T_{h}}/{T_{\gamma}}\ll 1$
and the above contribution is negligible.

Recall that the simplified model discussed in Section 2.1, which considers the
kinetic mixing between a (massless) $U(1)_{X}$ and the $U(1)_{\rm em}$. In the
mass eigenbasis the dark particle carries a millicharge
$g_{X}Q_{\chi}\delta/\sqrt{1-\delta^{2}}$, c.f., Eq. (11), where $g_{X}$ is
the $U(1)_{X}$ gauge coupling and $Q_{\chi}$ is the $U(1)_{X}$ charge carried
by the dark fermion.

In the literature, people usually takes the above expression for the
millicharge carried by the dark fermion even the $U(1)_{X}$ is not massless.
As we discussed in depth in Section 2, this treatment is inappropriate and
indeed, there would be no millicharge generated at all for the case of the
massless photon field mixing with an additional massive $U(1)$ via only
kinetic terms. Not mentioning considering $U(1)_{X}$ mixing directly with the
photon field, which is not a fundamental field before the electroweak symmetry
breaking, is not rigorous in the first place. Nevertheless, experiments can
still set bounds directly to the millicharge $Q_{\bm{\epsilon}}$ (if there is
any) carried by the dark particle, given by the expression

$Q_{\bm{\epsilon}}=g_{X}Q_{\chi}\delta/\sqrt{1-\delta^{2}}\equiv\bm{\epsilon}e\,.$
(106)

for the simplified model discussed in Section 2.1, and in the literature
bounds were set directly to the parameter $\bm{\epsilon}$.

Considering the mixing effects within the full electroweak theory, where the
$U(1)_{X}$ connects with $U(1)_{Y}$ through both kinetic mixing and mass
mixing, the dark matter carries a millicharge given by Eq. (49),

$Q_{\varepsilon}=\epsilon g_{X}Q_{\chi}\cos\theta\equiv\varepsilon e\,,$ (107)

giving rise to $\varepsilon=\epsilon g_{X}Q_{\chi}\cos\theta/e$. The
experimental constraints set on $\varepsilon$ can be then transformed into the
constraints on the mass mixing parameter $\epsilon$ in our models.

The current experimental constraints on the millicharge carried by the dark
particle with respect to the dark particle mass is plotted in Fig. 8. In the
model we discuss, the millicharge is proportional to the mass mixing parameter
$\epsilon$, and the constraints on $\epsilon$ can be obtained directly from
Eq. (107). We also mark our benchmark models on Fig. 8 indicating the
millicharge carried by the dark particles in the respective models.

Case | Model | $M_{\gamma^{\prime}}$ | $m_{\chi}$ | $g_{X}$ | $\delta$ | $\epsilon$ | $\varepsilon$ | $\Omega_{\chi}h^{2}$ | $\Omega_{\gamma^{\prime}}h^{2}$ | $\tau_{\gamma^{\prime}}$
---|---|---|---|---|---|---|---|---|---|---
1 | $m_{\chi}>M_{\gamma^{\prime}}>2m_{e}$ | $a$ | $20$ | $100$ | $0.0054$ | $1\times 10^{-13}$ | $1\times 10^{-10}$ | $1.57\times 10^{-12}$ | $0.120$ | $0$ | $0.616$
$b$ | $180$ | $250$ | $0.015$ | $1\times 10^{-12}$ | $1\times 10^{-9}$ | $4.36\times 10^{-11}$ | $0.120$ | $0$ | $6.82\times 10^{-4}$
2 | $2m_{\chi}>M_{\gamma^{\prime}}>m_{\chi}>2m_{e}$ | $c$ | $100$ | $60$ | $1.59\times 10^{-5}$ | $1\times 10^{-14}$ | $5\times 10^{-11}$ | $2.29\times 10^{-15}$ | $0.120$ | $0$ | $0.491$
3 | $M_{\gamma^{\prime}}>2m_{\chi}>2m_{e}$ | $d$ | $100$ | $10$ | $0.01$ | $1\times 10^{-14}$ | $5.6\times 10^{-13}$ | $1.62\times 10^{-14}$ | $0.120$ | $0$ | $2.48\times 10^{-19}$
4 | $2m_{e}>m_{\chi}>M_{\gamma^{\prime}}$ | $e$ | $0.09$ | $1$ | $0.20$ | $1\times 10^{-14}$ | $1.27\times 10^{-12}$ | $7.34\times 10^{-13}$ | $7.42\times 10^{-12}$ | $4.43\times 10^{-3}$ | $2.67\times 10^{30}$
5 | $2m_{e}>2m_{\chi}>M_{\gamma^{\prime}}>m_{\chi}$ | $f$ | $0.09$ | $0.06$ | $0.01$ | $1\times 10^{-14}$ | $1.27\times 10^{-12}$ | $3.67\times 10^{-14}$ | $5.94\times 10^{-3}$ | $2.58\times 10^{-9}$ | $2.67\times 10^{30}$
$g$ | $0.09$ | $0.06$ | $1.5\times 10^{-4}$ | $1\times 10^{-14}$ | $1.27\times 10^{-12}$ | $5.50\times 10^{-16}$ | $1.85\times 10^{-3}$ | $3.03\times 10^{-3}$ | $2.67\times 10^{30}$
6 | $2m_{e}>M_{\gamma^{\prime}}>2m_{\chi}$ | $h$ | $1$ | $0.05$ | $0.001$ | $4\times 10^{-13}$ | $1\times 10^{-11}$ | $2.9\times 10^{-14}$ | $0.120$ | $0$ | $6.20\times 10^{-18}$

Table 1: The benchmark models we consider in this work for six different types
of models. We compute the full evolution of all above benchmark models and
calculate the corresponding dark matter relic density for each model. One
cannot achieve a full occupation of the dark matter relic density for models
in the shaded blue region. The lifetimes (in the unit of seconds) of the dark
photon for each model are listed in the last column. $M_{\gamma^{\prime}}$ and
$m_{\chi}$ are in MeVs.

We consider the following benchmark models, shown in Table 1, to illustrate
our formalism of computing the hidden sector evolution. In the models we
consider, after the mixing, the dark fermion carried the millicharge
proportional to $\epsilon$, and $Z$ boson couples to the dark fermion with
strength proportional to $\delta$, c.f., Eqs. (50) and (51). These couplings
play important roles in calculating the evolution of the dark particles.
Specifically, the couplings of dark photon to SM fermions are proportional to
$|\delta-\epsilon|$, c.f., Eq. (52). Our investigation shows for each set of
the $\\{M_{\gamma^{\prime}},m_{\chi},g_{X}\\}$ value, the evolution of the
hidden sector temperature and the dark particles makes no difference for the
two cases $\delta\gg\epsilon$ and $\epsilon\gg\delta$, given the same value of
$|\delta-\epsilon|$. The same value of $|\delta-\epsilon|$ results in almost
identical shape of the plots for the two cases of $\delta\gg\epsilon$ and
$\epsilon\gg\delta$, although there is a separate constraint on $\epsilon$,
c.f., Fig. 8. In this paper we will not discuss the fine-tuned case that both
of $\delta$ and $\epsilon$ are large but their difference takes a tiny value,
i.e., $\delta\sim\epsilon\gg|\delta-\epsilon|$. For the following plots we
will only show the case of $\epsilon\gg\delta$ for a chosen value of
$|\delta-\epsilon|$, displayed in Fig. 7, which leads to a larger millicharge
carried by the dark fermion.

###  4.3 Evolution of hidden sector particles

In this subsection, we discuss in depth the evolution of the hidden sector
particles for eight benchmark models covering six distinct cases summarized in
the previous subsection. We calculate the coupled Boltzmann equations given by
Eqs. (84)-(86) and plot the evolution for the dark fermion $\chi$ and the dark
photon $\gamma^{\prime}$, as well as interaction rates for various hidden
sector interactions related to the dark freeze-out, versus Hubble. In the
calculation we choose an initial value of $\eta=10$ at the temperature
$T=10^{5}$ GeV. As discussed in [37, 64], the initial value of $\eta=T/T_{h}$
will not modify the final abundance so much but will alter the details of the
dark particles evolution. A comprehensive analysis for the evolution of a
$U(1)$ hidden sector interacting feebly with the SM from the reheating stage,
including the determination of the initial value of $\eta$, is discussed in a
companion paper [65].

Figure 9:  [Color online] The exhibition of the hidden sector evolution for
Case 1 (models a,b). For each model, the above figure shows the evolution of
the comoving number densities ($Y$) of the dark fermion $\chi$ and the dark
photon $\gamma^{\prime}$. The red horizontal dash line denotes the observed
dark matter relic density with dark matter mass equal to $m_{\chi}$. In the
lower figure we trace the interaction rates of three important processes
compared to the Hubble line. When the interaction rate is above the Hubble
line, the corresponding interaction takes place very rapidly and particles
participated in the interaction are in chemical equilibrium. The solid black
vertical line shows at this temperature the thermally averaged decay width of
the dark photon to SM fermions (mostly to the electron-positron pairs) going
above Hubble, indicating the decay of dark photon becomes rapid. The dash
black vertical line corresponds to the process of
$\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$ becoming inactive. The dash-
dotted black vertical line corresponds to the process of
$\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}$ becoming inactive.

##### Case 1 (dark matter $\chi$):

As we can see from Fig. 7, for the above 1 MeV region, the dark photon
receives the strongest constraint at the mass around 17 MeV, corresponding to
a constraint $|\delta-\epsilon|\lesssim 2\times 10^{-11}$. Using this lower
bound value, the dark photon decays not so rapidly and may spoil the BBN
constraint. Thus we choose the dark photon mass to be 20 MeV, a little off the
17 MeV point in the benchmark model a, corresponding to the lifetime of the
dark photon less than 1 second. When the dark photon mass is far from 20 MeV,
there would less restrictions on $|\delta-\epsilon|$, c.f., Fig. 7. In model
b, we take a larger value of $|\delta-\epsilon|$, resulting in a larger
coupling between the visible sector and the hidden sector. We show benchmark
models a,b for Case 1 in Fig. 9. In both models, the dark fermion can make up
the whole dark matter constituent.

For each benchmark model, we plot the evolution of the dark fermion $\chi$ and
the dark photon $\gamma^{\prime}$, and with the same temperature scales we
also plot the interaction rates of hidden sector particle evolutions versus
the Hubble line. When the interaction rate is above the Hubble scale at a
certain temperature, the corresponding interaction is quite intense and
particles involving in this interaction are in thermal equilibrium. For model
a: The black vertical dash line corresponds to the process of
$\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$ becoming inactive,
indicating the start of the dark freeze-out, represented by the number density
of $\chi$ about to drop; The black vertical dash-dotted line corresponds to
the process of $\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}$ becoming
inactive, indicating the finish of the dark freeze-out, represented by the
number density of $\chi$ tending to level off; The black vertical solid line
shows that at this temperature the thermally averaged decay width of the dark
photon to SM fermions (mostly to the electron-positron pairs) goes above
Hubble, indicating the decay of the dark photon becoming rapid, which
corresponds to a steep drop for the dark photon number density.

For model b, the constraint on $|\delta-\epsilon|$ is less restrictive and can
be taken to be a larger value compared to model a. Thus the freeze-in rate is
much larger and more dark particles are produced. In the plots, one can see
that even the purple line showing the dark photon decay becomes active in a
quite early stage, the robust freeze-in production still generates a massive
amount of both the dark fermion $\chi$ and the dark photon $\gamma^{\prime}$.
The evolutions of $\chi$ and $\gamma^{\prime}$ are determined by the combined
effects of freeze-in, hidden sector interactions and the dark photon decay to
SM particles. The dark photon finally diminishes to vanish before BBN.

Figure 10:  [Color online] The exhibition of the hidden sector evolution for
Case 2 (model c). The above figure shows the evolution of the comoving number
densities of the dark fermion $\chi$ and the dark photon $\gamma^{\prime}$,
with the red horizontal line denoting the observed dark matter relic density
with dark matter mass equal to $m_{\chi}$. The lower figures show the
interplay between the hidden sector interaction rates versus Hubble.

##### Case 2 (dark matter $\chi$):

For case 2 we study benchmark model c, and the full evolution is shown in Fig.
10. In this case, for a relatively large $g_{X}$, the hidden sector will have
strong freeze-out effect via the process
$\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$ and render an overproduction
of the dark matter $\chi$. The only way to achieve the observed dark matter
relic density is to weaken the hidden sector interaction
$\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$ with a rather small value of
$g_{X}$, such that less $\gamma^{\prime}$ will transform to the dark fermion
$\chi$, noticing that $\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$
dominates over $\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}$ as seen from
the lower figure in Fig. 10. The upper plot in Fig. 10 shows the evolution of
the hidden sector particles, whereas the lower plot exhibits the hidden sector
interaction rates are below the Hubble line, corresponding to the hidden
sector interactions being inactive. However, the effect of hidden sector
interactions play crucial roles generating the correct amount of dark matter
relic density. Later, when the dark photon decay began to dominate, shown by
the solid black vertical line, the number density of the dark photon
immediately drops down to zero. The abundance of the dark fermion, however,
gradually accumulates without reduction, and can finally reach the current
observed value of the dark matter relic density, given by suitable parameters
we choose.

Comparison of different calculations | $\Omega_{\chi}h^{2}$
---|---
All freeze-in processes included | 0.1195
Plasmon decay process excluded | 0.1195
Four-point $\gamma^{\prime}$ freeze-in processes excluded | 0.0643
Pure freeze-in for $\chi$ | $10^{-9}$
Table 2: A comparison of different calculations of the dark matter relic
density for the benchmark model c.

Since in this case the parameters are chosen in such a way that the hidden
sector interaction
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ never achieved
thermal equilibrium, it is a good occasion to show the contribution of the
four-point freeze-in processes of the dark photon which are usually ignored in
the literature, given by the last line in Eq. 85, as long as the three-point
freeze-in productions of the dark photon ($i\bar{i}\to\gamma^{\prime}$) are
present. In Table 2 we compare several different computations with or without
some important freeze-in processes. The first result corresponds to the
calculated relic density for the dark fermion $\chi$ including all relevant
processes given by Eqs. (84) and (85). The second result shows the $\chi$
relic density not including the contribution from the plasmon decay, which is
almost identical to the first result including all the processes, which agrees
with the previous findings that the production of the millicharge particle
from plasmon decay becomes important only when the millicharge particle is
lighter than the mass of the electron [43, 46]. The third result shows the
$\chi$ relic density without including all the four-point freeze-in processes
of the dark photon. We conclude: (1) although the interaction rates of
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ are always below
Hubble, the dark photon freeze-in production plays a crucial role in
generating the abundance of the dark fermion $\chi$, and (2) the four-point
freeze-in production of the dark photon is as important as the three-point
freeze-in production channels, and should not be neglected. The last result is
computed from the pure freeze-in production of $\chi$, ignoring the freeze-in
production of the dark photon and dropping out all hidden sector interactions,
which is a much smaller value compared with the previous one. This can be
understood as follows: the freeze-in production of $\chi$ depends on the
interaction $i\bar{i}\to\chi\bar{\chi}$, which always has an additional
$g_{X}^{2}$ suppression compared to the dark photon freeze-in production,
c.f., the rough estimation of the couplings in Eqs. (50)-(52). In model c we
choose a rather small $g_{X}\sim 10^{-5}$ and thus the pure freeze-in
production of $\chi$ is quite rare.

Figure 11:  [Color online] The exhibition of the hidden sector evolution for
Case 3 (model d). The top figure shows the evolution of the comoving number
densities of the dark fermion $\chi$ and the dark photon $\gamma^{\prime}$,
with the red horizontal line denoting the observed dark matter relic density
with dark matter mass equal to $m_{\chi}$. In the lower figure we trace the
hidden sector interaction rates compared to the Hubble line. The black
vertical solid line shows at this temperature the thermally averaged decay
width of $\gamma^{\prime}\to\chi\bar{\chi}$ becoming active, and the
corresponding interaction rate is extremely strong. The black vertical dash
line shows the process $\chi\bar{\chi}\to\gamma^{\prime}$ becoming active and
the black vertical dash-dotted line corresponds to the same process becoming
insignificant and all of the dark photon will eventually decay to dark
fermions.

##### Case 3 (dark matter $\chi$):

In this case the dark photon mass is larger than twice of the dark fermion
mass and thus a vigorous decay of the dark photon dominates all the
interaction channels. Especially, the decay of the dark photon to SM fermion
pairs $\gamma^{\prime}\to i\bar{i}$ can be completely neglected compared to
the process $\gamma^{\prime}\to\chi\bar{\chi}$ occurring at the same time. As
can be seen from Fig. 11, the abundance of the dark fermion $\chi$ always
increases before it reaches the final value. The first drop in the dark photon
abundance is due to $\gamma^{\prime}\to\chi\bar{\chi}$ becoming active,
indicated by the black vertical solid line. When the interaction rate of the
process $\chi\bar{\chi}\to\gamma^{\prime}$ goes above the Hubble line,
indicated by the black vertical dash line, the abundance of the dark photon
soon increases. The black vertical dash-dotted line corresponds to the process
of $\chi\bar{\chi}\to\gamma^{\prime}$ becoming insignificant and finally all
dark photon decay into $\chi\bar{\chi}$ and its abundance becomes zero.

Figure 12:  [Color online] The exhibition of the hidden sector evolution for
Case 4 (model e) and Case 5 (models f,g). For each model, the above figure
shows the evolution of the comoving number densities of the dark fermion
$\chi$ and the dark photon $\gamma^{\prime}$. The red horizontal line denoting
the observed dark matter relic density if the dark matter is totally made up
by the dark fermion $\chi$. whereas the blue horizontal line shows the
observed dark matter relic density if the dark matter is totally made up by
the dark photon $\gamma^{\prime}$. In the lower figures for each models, we
trace the interaction rates of two processes among the hidden sector compared
to the Hubble. The IDPB constraint is so strong that one can only take a much
smaller value for both $\delta$ and $\epsilon$ compared with other cases. As a
result, in Case 4 and Case 5 one cannot obtain a large enough value for the
dark matter relic abundance. In model g, we choose a certain $g_{X}$ value
such that the hidden sector interaction rate is about tangent to the Hubble
line, and we obtain the largest relic abundance for the dark matter
$\gamma^{\prime}$ and $\chi$.

##### Case 4 and Case 5 (dark matter $\chi,\gamma^{\prime}$):

In these two cases the dark photon mass is less than twice of the electron
mass and thus the dark photon can be stable throughout the history of the
Universe, as discussed in Section 2.4.1. As noted earlier, the IDPB constraint
is stronger than the experimental constraints on both the kinetic mixing and
the millicharge carried by the dark matter, and thus restrict the
$|\delta-\epsilon|$ value to be much smaller than Case 1,2,3 and 6. In these
two cases, both the dark fermion $\chi$ and the dark photon $\gamma^{\prime}$
are dark matter candidates. We conclude, for a single extra $U(1)$ mixed with
the SM, one cannot achieve a full occupation of the dark matter relic density
with the dark photon dark matter (at most about $5\%$). For a full occupation
of the dark photon dark matter, one needs at least two extra $U(1)$ hidden
sector and was studied in depth in [38].

In Case 4, we choose the benchmark model e such that they can produce the
highest abundance of the dark matter, corresponding to a rather large hidden
sector gauge coupling $g_{X}$. In this case one can see clearly a dark freeze-
out occurring right after the interaction rates of
$\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$ and
$\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}$ start to deviate. The dark
freeze-out ended at the interaction
$\chi\bar{\chi}\to\gamma^{\prime}\gamma^{\prime}$ become inactive, indicated
by the dash dotted line.

For Case 5 we choose two types of benchmark models: one type (model f) with
larger $U(1)_{X}$ gauge coupling and the hidden sector interactions are active
and the dark freeze-out for the dark photon takes place via the interaction
$\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$, where we see a quick drop
in $Y_{\gamma^{\prime}}$ around $T_{h}\sim 10^{-2}$ MeV. The dark freeze-out
ended when $\gamma^{\prime}\gamma^{\prime}\to\chi\bar{\chi}$ becoming inactive
indicated by the black vertical dash line. The $U(1)_{X}$ gauge coupling of
the other type (models g) is chosen to be small such that the hidden sector
interactions $\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$
never go beyond Hubble. In model g, we find the largest relic density for the
dark matter $\gamma^{\prime}$ and $\chi$. Thus for case that the SM extended
by a single $U(1)$, the $\mathcal{O}$(keV) dark photon dark matter can only
occupy at most $\sim 5\%$ of the total observed dark matter relic density.

As the dark fermion mass goes below the mass of the electron, the freeze-in
contribution from plasmon decay $\gamma^{*}\to\chi\bar{\chi}$ becomes
important [43, 46], especially when the dark fermion mass is around 50 keV.
However, the dark fermion number density cannot be computed separately, even
the interactions inside the hidden sector are ultraweak, as we already
discovered in the analysis for the benchmark model c. To calculate the final
relic abundance of the dark fermion, one needs to compute the full evolution
of the entire hidden sector including all of the hidden sector interactions.
The hidden sector interactions
$\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$ play crucial
roles in determining the final relic abundance of both $\chi$ and
$\gamma^{\prime}$.

In Table 3 we list the final result of the relic abundance of $\chi$ and
$\gamma^{\prime}$ with and without the inclusion of the plasmon decay. We can
see the contribution of the plasmon decay to the final relic densities of
$\chi$ and $\gamma^{\prime}$ are indeed not significant, because of the
interactions inside the hidden sector.

Benchmark models | Plasmon decay included | Plasmon contribution excluded
---|---|---
$\Omega_{\chi}h^{2}$ | $\Omega_{\gamma^{\prime}}h^{2}$ | $\Omega_{\chi}h^{2}$ | $\Omega_{\gamma^{\prime}}h^{2}$
Model e | $7.4193\times 10^{-12}$ | $4.4344\times 10^{-3}$ | $7.4191\times 10^{-12}$ | $4.4327\times 10^{-3}$
Model f | $5.9415\times 10^{-3}$ | $2.5842\times 10^{-9}$ | $5.9397\times 10^{-3}$ | $2.5807\times 10^{-9}$
Table 3: The results of dark matter relic densities including or excluding the
plasmon decay contribution for benchmark models e and f. In these two
benchmark models, both the dark fermion $\chi$ and the dark photon
$\gamma^{\prime}$ are dark matter candidates. Figure 13:  [Color online] The
exhibition of the hidden sector evolution for Case 6 (model h). The above
figure shows the evolution of the comoving number densities of the dark
fermion $\chi$ and the dark photon $\gamma^{\prime}$, with the red horizontal
line denoting the observed dark matter relic density with the dark matter
being the dark fermion $\chi$. The lower figure shows the interaction rates of
two processes among the hidden sector compared to the Hubble. In this case,
even the mass of the dark photon is taken to be less than twice of electron
mass, the massive decay of the dark photon to dark fermions plays the dominant
role.

##### Case 6 (dark matter $\chi$):

In this case, the intensive decay of the dark photon through channel
$\gamma^{\prime}\to\chi\bar{\chi}$ plays the dominant role. This case is
similar to Case 3: the abundance of the dark fermion $\chi$ always increases
before it reaches the final value. The drop in the dark photon abundance
around the black vertical dash line is due to the inverse interaction
$\chi\bar{\chi}\to\gamma^{\prime}$ becoming active. Finally all dark photon
decay into $\chi\bar{\chi}$ and its abundance becomes zero.

In this section, we investigate eight benchmark models (a-h) from six distinct
cases depending on the masses of the dark photon and the dark fermion, given
by Table 1. In models a,b,c,d,h, the dark matter candidate is the dark fermion
$\chi$; whereas in models e,f,g, the dark matter candidates are the dark
fermion $\chi$ and the dark photon $\gamma^{\prime}$. We have computed and
plotted the full evolution of the dark fermion $\chi$ and the dark photon
$\gamma^{\prime}$ as well as the hidden sector temperature, governed by the
coupled Boltzmann equations given in Section 3.2. We also keep track of the
hidden sector interaction rates versus Hubble to see the intensity of the
corresponding interaction. In models a,b,d,e,f,h we find the hidden sector
encountered apparent dark freeze-out.

We now succinctly outline the primary discoveries for the analysis of the
eight benchmark models encompassing six distinct cases, depending on the
masses of the dark photon and the dark fermion, given by Table 1:

  * •

For the extension of the SM with an extra $U(1)$ hidden sector, we find the
contribution from four-point freeze-in processes for the dark photon:
$i\bar{i}\to\gamma\gamma^{\prime}$, $i\gamma\to i\gamma^{\prime}$,
$\bar{i}\gamma\to\bar{i}\gamma^{\prime}$, which are usually ignored in the
literature, are as important as the three-point freeze-in processes
$i\bar{i}\to\gamma^{\prime}$, where $i$ denotes the SM fermions (see Table 2,
the third result). We also expect that for more general freeze-in dark matter
models, the four-point freeze-in production channels are also crucial and
should not be neglected, even the three-point freeze-in channels for
generating the same freeze-in particle are present.

  * •

If the hidden sector particles are produced via the freeze-in mechanism, as
long as the hidden sector has self-interactions involving various dark
particles, the self-interactions within the hidden sector (even ultraweak)
play significant roles in the evolution of the hidden sector (see Table 2, the
fourth result). One must take full consideration of the contribution from
hidden sector self-interactions at all times.

  * •

The freeze-in production of the dark fermion $\chi$ from the plasmon decay,
however, does not play an important role in determining the final relic
abundance of the $\chi$, as long as the hidden sector self-interactions (in
this case: $\chi\bar{\chi}\leftrightarrow\gamma^{\prime}\gamma^{\prime}$) are
present, see results from Table 2 and Table 3. Though for a comprehensive
analysis, it is advisable to incorporate the plasmon contribution into the
computation.

  * •

In the analysis of Case 4 and Case 5 (models e,f,g), we find
$\mathcal{O}$(keV) dark photon dark matter from a single $U(1)_{X}$ hidden
sector with feeble couplings to the SM, can only occupy at most $\sim 5\%$ of
the total observed dark matter relic density.

##  5 Conclusion

In this paper we conduct a comprehensive study on the sub-GeV millicharge dark
matter produced via the freeze-in mechanism. We discuss the SM extended by a
$U(1)_{X}$ hidden sector, and there exist tiny kinetic mixing and mass mixing
between the $U(1)_{X}$ and the hypercharge gauge field. The dark fermions
carrying the $U(1)_{X}$ quantum numbers are naturally dark matter candidates.
After the mass mixing of the $U(1)_{X}$ and the hypercharge gauge field, the
dark fermion carries a millicharge. We investigate sub-GeV millicharge dark
matter candidates produced via the freeze-in mechanism. We discuss in general
how such millicharge dark matter can be generated from a theoretical
perspective. We highlight that there are several inappropriate treatments and
misunderstandings in the literature regarding this subject.

The millicharge can be generated only in three ways:

  1. 1.

The dark particle carries a tiny amount of hypercharge as a prior, and one
such possibility is discussed in Section 2.4.2.

  2. 2.

Theory involving kinetic mixing between a massless $U(1)$ and the hypercharge
gauge field, and the generated millicharge is proportional to the kinetic
mixing parameter, as discussed in Section 2.4.3.

  3. 3.

The mass mixing between two massive $U(1)$ (one may produce a massless $U(1)$
in the final mass eigenbasis, see the discussion in Section 2.3), no matter
the kinetic mixing is present or not, and the millicharge generated is
proportional to the mass mixing parameter. In the electroweak theory, this can
be done by a Stueckelberg mass mixing between the $U(1)_{X}$ and the
hypercharge gauge field, as discussed in Section 2.4.1.

Inappropriate treatments and misunderstandings include:

  1. 1.

It is not appropriate to consider the kinetic mixing between the extra $U(1)$
and the photon field, especially from a theoretical perspective. The photon
field is a linear combination of the hypercharge gauge field and the neutral
$SU(2)$ gauge boson in the full electroweak theory. Thus the mixing between
the extra $U(1)$ and the hypercharge gauge field may not generate a direct
mixing between the extra $U(1)$ and the $U(1)_{\rm em}$ [28]. It is indeed
subtle to embed such effective term
$\mathcal{L}\sim\frac{\bm{\delta}}{2}F_{\mu\nu}^{\rm em}F_{X}^{\mu\nu}$ into a
UV-completed theory, as discussed in depth in Section 2.4.1.

  2. 2.

Even from the effective field theory point of view, the kinetic mixing between
the massless photon field and an extra massive $U(1)$ cannot generate a
millicharge, no matter how tiny the extra $U(1)$ mass is [28], summarised in
Section 2.2. Taking the mass of the massive $U(1)$ to be extremely small (but
still nonzero), cannot transit to the case of kinetic mixing between two
massless $U(1)$’s, which will generate a millicharge proportional to the
kinetic mixing parameter.

  3. 3.

The kinetic mixing does not play any role in generating a millicharge as long
as one of the two $U(1)$’s gets massive. In this case only the mass mixing,
e.g., Eq. (24) can generate a millicharge.

In the model exhibited in Fig. 2, both the dark fermion $\chi$ and the dark
photon $\gamma^{\prime}$ are produced via the freeze-in mechanism. After the
mass mixing, $\chi$ carries a millicharge in addition to the $U(1)_{X}$
quantum number. $\chi$ is always stable and thus is the dark matter candidate,
whereas the dark photon can be made stable when its mass is less than twice of
the electron mass. Within the hidden sector, $\chi$ and $\gamma^{\prime}$ may
simultaneously engage in strong or weak interactions during the freeze-in
production of $\chi$ and $\gamma^{\prime}$. The interactions within the hidden
sector alter the number densities of $\chi$ and $\gamma^{\prime}$ over time.
Hence, it is not possible to calculate the relic abundance of $\chi$ and
$\gamma^{\prime}$ solely through pure freeze-in production. With the formalism
established in [37] and reviewed in Section 3, we are capable of calculating
the complete evolution of the hidden sector particles as well as the hidden
sector temperature. This calculation formalism is general and can be applied
to models involving any hidden sectors with feeble connections to the visible
sector. In these models, the hidden sector particles are produced through the
freeze-in mechanism, and exhibit strong or weak (even ultraweak) self-
interactions within the hidden sector.

We perform a comprehensive analysis on eight benchmark models (a-h)
encompassing six distinct cases, depending on the masses of the dark photon
and the dark fermion, and determine the complete evolution of the $U(1)_{X}$
hidden sector. In some of the models with intensive hidden sector interactions
(models a,b,d,e,f,h), we observe apparent dark freeze-out; while in others
with weak or even ultraweak hidden sector interactions (models c,g), the dark
freeze-out is not apparent. The intensities of the hidden sector interactions
are characterized by the interaction rates, and are related to the hidden
sector temperature which is evidently different from the visible sector
temperature (the temperature of our observed Universe). The hidden sector
temperature can be tracked in the calculation that solving the coupled
Boltzmann equations given in Section 3.2. We have taken into account all of
the freeze-in processes as well as the hidden sector self-interactions
simultaneously. We plot the evolution of the hidden sector particles with an
emphasize on the dark particle interactions within the hidden sector.

The primary discoveries for the analysis of the eight benchmark models are
summarised at the end of Section 4.3. We tend to believe the following
findings are universal and applicable to broader freeze-in scenarios:

  1. 1.

For extensions of the SM, the four-point freeze-in processes are usually
ignored, especially when the three-point freeze-in production channels are
present for the same freeze-in particle. In our analysis we find that the
four-point freeze-in production channels are as important as the three-point
freeze-in channels, and should not be neglected.

  2. 2.

If the hidden sector possesses self-interactions (even ultraweak), these
interactions invariably occur simultaneously with the freeze-in production of
hidden sector particles. To accurately determine the relic abundance of the
freeze-in dark matter, it is essential to consider all pertinent freeze-in
processes, as well as all hidden sector interactions, and solve the coupled
Boltzmann equations while accounting for the temperature disparity between the
visible and hidden sectors.

Acknowledgments:

WZF acknowledges discussions with Xiaoyong Chu and Zuowei Liu. WZF is
supported by the National Natural Science Foundation of China under Grant No.
11935009.

##  Appendix A Kinetic mixing and mass mixing in the Stueckelberg extension

In this Appendix we derive the mass eigenbasis of neutral gauge bosons in the
Stueckelberg extension of the SM with kinetic mixing and mass mixing, which is
firstly done in [28]. We consider the Lagrangian given in Eqs. (3) and (4),
resulting in the mixing matrices Eq. (35). We consider the case that
$M_{2}\neq 0$ thus we have a tiny mass mixing parameterized by
$\epsilon=M_{2}/M_{1}$. In this case the physical photon in the mass
eigenbasis couples to hidden sector fermions with a strength proportional to
$\epsilon$, thus the hidden sector fermion is millicharged. A simultaneous
diagonalization of the kinetic matrix and mass matrix can be done in the
following steps. Firstly we rewrite the mass matrix in terms of $\epsilon$ as

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

A transformation

$K^{\prime}=KO^{\prime}=\begin{pmatrix}\frac{1}{\sqrt{1-\delta^{2}}}&0&0\\\
-\frac{\delta}{\sqrt{1-\delta^{2}}}&1&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}\cos\phi&-\sin\phi&0\\\
\sin\phi&\cos\phi&0\\\
0&0&1\end{pmatrix}=\begin{pmatrix}\frac{1-\delta\epsilon}{\sqrt{\left(1-\delta^{2}\right)\left(1-2\delta\epsilon+\epsilon^{2}\right)}}&\frac{-\epsilon}{\sqrt{1-2\delta\epsilon+\epsilon^{2}}}&0\\\
\frac{\epsilon-\delta}{\sqrt{\left(1-\delta^{2}\right)\left(1-2\delta\epsilon+\epsilon^{2}\right)}}&\frac{1}{\sqrt{1-2\delta\epsilon+\epsilon^{2}}}&0\\\
0&0&1\end{pmatrix}\,,$ (A.2)

with

$\phi=\arctan\left(\frac{\epsilon\sqrt{1-\delta^{2}}}{1-\delta\epsilon}\right)\,,$
(A.3)

diagonalizes the kinetic matrix and change the mass matrix into the form

$K^{\prime T}M_{St}^{2}K^{\prime}=\begin{pmatrix}\frac{v^{2}g_{Y}^{\prime
2}(\delta-\epsilon)^{2}+4M_{1}^{2}\left(1-2\delta\epsilon+\epsilon^{2}\right)}{4\left(1-\delta^{2}\right)}&\frac{v^{2}g_{Y}^{\prime
2}(\epsilon-\delta)}{4\sqrt{1-\delta^{2}}}&\frac{v^{2}\mathrm{g_{2}}g_{Y}^{\prime}(\delta-\epsilon)}{4\sqrt{\left(1-\delta^{2}\right)}}\\\
\frac{v^{2}g_{Y}^{\prime
2}(\epsilon-\delta)}{4\sqrt{1-\delta^{2}}}&\frac{v^{2}g_{Y}^{\prime
2}}{4}&-\frac{v^{2}\mathrm{g_{2}}g_{Y}^{\prime}}{4}\\\
\frac{v^{2}\mathrm{g_{2}}g_{Y}^{\prime}(\delta-\epsilon)}{4\sqrt{\left(1-\delta^{2}\right)}}&-\frac{v^{2}\mathrm{g_{2}}g_{Y}^{\prime}}{4}&\frac{v^{2}\mathrm{g_{2}^{2}}}{4}\end{pmatrix}\,,$
(A.4)

where we set $g_{Y}^{\prime}=g_{Y}/\sqrt{1-2\delta\epsilon+\epsilon^{2}}$. Now
the lower right $2\times 2$ matrix is of the same form with the SM
$\\{B,A_{3}\\}$ mass mixing matrix. The rotation matrix

$R_{1}^{\prime}=\begin{pmatrix}1&0&0\\\ 0&\cos\theta&-\sin\theta\\\
0&\sin\theta&\cos\theta\end{pmatrix}=\begin{pmatrix}1&0&0\\\
0&\frac{g_{2}}{\sqrt{g_{2}^{2}+g_{Y}^{\prime
2}}}&\frac{-g_{Y}^{\prime}}{\sqrt{g_{2}^{2}+g_{Y}^{\prime 2}}}\\\
0&\frac{g_{Y}^{\prime}}{\sqrt{g_{2}^{2}+g_{Y}^{\prime
2}}}&\frac{g_{2}}{\sqrt{g_{2}^{2}+g_{Y}^{\prime 2}}}\end{pmatrix}\,,$ (A.5)

will transform the mass matrix into

$\left(K^{\prime}R_{1}^{\prime}\right)^{T}M_{St}^{2}K^{\prime}R_{1}^{\prime}=\begin{pmatrix}\frac{v^{2}g_{Y}^{\prime
2}(\delta-\epsilon)^{2}+4M_{1}^{2}\left(1-2\delta\epsilon+\epsilon^{2}\right)}{4\left(1-\delta^{2}\right)}&0&\frac{1}{4}v^{2}g_{Y}^{\prime}\left(\delta-\epsilon\right)\sqrt{\frac{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}{1-\delta^{2}}}\\\ 0&0&0\\\
\frac{1}{4}v^{2}g_{Y}^{\prime}\left(\delta-\epsilon\right)\sqrt{\frac{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}{1-\delta^{2}}}&0&\frac{1}{4}v^{2}\left(\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}\right)\end{pmatrix}\,.$ (A.6)

A second rotation matrix

$R_{2}^{\prime}=\begin{pmatrix}\cos\psi&0&\sin\psi\\\ 0&1&0\\\
-\sin\psi&0&\cos\psi\end{pmatrix}$ (A.7)

will fully diagonalize the above matrix with

$\psi=\frac{1}{2}\arctan\frac{2v^{2}g_{Y}^{\prime}\left(\delta-\epsilon\right)\sqrt{\left(1-\delta^{2}\right)\left(\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}\right)}}{-4M_{1}^{2}\left(1-2\delta\epsilon+\epsilon^{2}\right)+v^{2}\left[\mathrm{g}_{2}^{2}\left(1-\delta^{2}\right)+g_{Y}^{\prime
2}\left(1-2\delta^{2}+2\delta\epsilon-\epsilon^{2}\right)\right]}\,.$ (A.8)

The full rotation matrix with mass mixing is then given by

$\displaystyle\mathcal{R}=KO^{\prime}R_{1}^{\prime}R_{2}^{\prime}$
$\displaystyle=\begin{pmatrix}\frac{1}{\sqrt{1-\delta^{2}}}&0&0\\\
-\frac{\delta}{\sqrt{1-\delta^{2}}}&1&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}\cos\phi&-\sin\phi&0\\\
\sin\phi&\cos\phi&0\\\ 0&0&1\end{pmatrix}\begin{pmatrix}1&0&0\\\
0&\cos\theta&-\sin\theta\\\
0&\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}\cos\psi&0&\sin\psi\\\
0&1&0\\\ -\sin\psi&0&\cos\psi\end{pmatrix}$
$\displaystyle=c\begin{pmatrix}\frac{\left(1-\delta\epsilon\right)\cos\psi}{\sqrt{1-\delta^{2}}}-\frac{g_{Y}^{\prime}\epsilon\sin\psi}{\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}&-\frac{g_{2}\epsilon}{\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}&\frac{g_{Y}^{\prime}\epsilon\cos\psi}{\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}+\frac{\left(1-\delta\epsilon\right)\sin\psi}{\sqrt{1-\delta^{2}}}\\\
\frac{(\epsilon-\delta)\cos\psi}{\sqrt{1-\delta^{2}}}+\frac{g_{Y}^{\prime}\sin\psi}{\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}&\frac{g_{2}}{\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}&-\frac{g_{Y}^{\prime}\cos\psi}{\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}+\frac{\left(\epsilon-\delta\right)\sin\psi}{\sqrt{1-\delta^{2}}}\\\
\frac{-g_{2}\sin\psi}{c\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}&\frac{g_{Y}^{\prime}}{c\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}&\frac{g_{2}\cos\psi}{c\sqrt{\mathrm{g}_{2}^{2}+g_{Y}^{\prime
2}}}\end{pmatrix}\,,$ (A.9)

where the constant $c=1/\sqrt{1-2\delta\epsilon+\epsilon^{2}}\approx 1$.

##  Appendix B Plasmon Effect

An electromagnetic wave propagating through the plasma consists of coherent
vibrations and encompasses both transverse and longitudinal components and
travels at speeds slower than that of light. The quantization of
electromagnetic waves within a plasma yields a spin-1 particle featuring two
transverse and one longitudinal spin polarizations. It is customary to refer
to all three polarization states as “plasmons” in the literature, underscoring
dispersion relations depend on the properties of the plasma. Thus in the early
universe, the photon acquires an effective mass
$m_{T}^{2}(k,T)=\omega_{T}^{2}(k,T)-k^{2}$,
$m_{L}^{2}(k,T)=\omega_{L}^{2}(k,T)-k^{2}$, for the transverse and
longitudinal modes respectively, where $k$ is the momentum of plasmon. The
massive photon, which we will refer to as plasmon in this paper, can decay
into light particle and anti-particle pair carrying electric charges. The dark
millicharged particles can then be also produced by the plasmon $\gamma^{*}$
via the decay $\gamma^{*}\rightarrow\chi\bar{\chi}$. In this appendix we
review briefly the plasmon effect and the effective plasmon frequency which is
used for calculating the decay of the plasmon. We follow the discussion of
[101], where the decay of plasmons to neutrinos was considered.

The typical electron velocity is defined by

$v_{\star}=\omega_{1}/\omega_{p}\,,$ (B.1)

where the first mode frequency and the plasma frequency are given by phase
space integrals

$\displaystyle\omega_{1}^{2}(T)$
$\displaystyle=\frac{4\alpha}{\pi}\int_{0}^{\infty}\mathrm{d}p\frac{p^{2}}{E}\left(\frac{5}{3}v^{2}-v^{4}\right)\big{[}n_{e}(E,T)+\bar{n}_{e}(E,T)\big{]}\,,$
(B.2) $\displaystyle\omega_{p}^{2}(T)$
$\displaystyle=\frac{4\alpha}{\pi}\int_{0}^{\infty}\mathrm{d}p\frac{p^{2}}{E}\left(1-\frac{v^{2}}{3}\right)\big{[}n_{e}(E,T)+\bar{n}_{e}(E,T)\big{]}\,,$
(B.3)

where $n_{e},\bar{n}_{e}$ are respectively the number densities of electron