#  MS-TP-22-15 A tower of hidden sectors: a general treatment and physics
implications

Amin<EMAIL_ADDRESS>and Pran
<EMAIL_ADDRESS>

aInstitut für Theoretische Physik, Westfälische Wilhelms-Universität Münster,
Wilhelm-Klemm-Straße 9, 48149 Münster, Germany
bDepartment of Physics, Northeastern University,
111 Forsyth Street, Boston, MA 02115-5000, USA

###### Abstract

An analysis of a tower of hidden sectors coupled to each other, with one of
these hidden sectors coupled to the visible sector, is given and the
implications of such couplings on physics in the visible sector are
investigated. Thus the analysis considers $n$ number of hidden sectors where
the visible sector couples only to hidden sector 1, while the latter couples
also to hidden sector 2, and the hidden sector 2 couples to hidden sector 3
and so on. A set of successively feeble couplings of the hidden sectors to the
visible sector are generated in such a set up. In general each of these
sectors live in a different heat bath. We develop a closed form set of coupled
Boltzmann equations for the correlated evolution of the temperatures and
number densities of each of the heat baths. We then apply the formalism to a
simplified model with scalar portals between the different sectors.
Predictions related to dark matter direct detection experiments and future CMB
probes of dark radiation are made.

###### Contents

  1. 1 Introduction
  2. 2 A model with $n$ hidden sectors
  3. 3 Theory application: The Stueckelberg extension of the Standard Model
    1. 3.1 Generation of very feeble couplings from kinetic mixings via a chain of hidden sectors
    2. 3.2 Neutral currents of the coupled visible and hidden sectors
  4. 4 Numerical application: The visible sector coupled to three hidden sectors
    1. 4.1 A simple model
    2. 4.2 Evolution equations
    3. 4.3 Numerical analysis
  5. 5 Implications for observables
    1. 5.1 Thermal effect of hidden sectors on direct detection of dark matter
    2. 5.2 Dark radiation from hidden sectors
  6. 6 Conclusion
  7. A Mass and temperature dependence of $g_{\rm eff}$ and $h_{\rm eff}$
  8. B Collision and source terms

##  1 Introduction

Most modern theories of particle physics based on supergravity, strings and
branes contain hidden sectors which are neutral under the Standard Model (SM)
gauge group but may possess their own gauge groups. However, even though the
hidden sectors are neutral under the SM gauge group, they may possess feeble
interactions with the visible sector via a variety of portals such as the
Higgs portal, kinetic and Stueckelberg mass mixings between the hidden sector
$U(1)$ gauge fields and the visible sector $U(1)_{Y}$ hypercharge gauge field,
or via higher dimensional operators which are dually charged under the gauge
groups of the Standard Model and of the hidden sectors. The visible and the
hidden sectors will in general reside in different heat baths. In prior
literature, it is often assumed (see, e.g. [1, 2]) that the visible and the
hidden sectors temperature evolution are governed by entropy conservation
carried out separately for the visible and the hidden sectors. However, such
an assumption is untenable if there is a coupling between the hidden sectors
and the visible sector. Thus a proper treatment of the temperatures of the
visible and the hidden sectors implies that one perform their evolution
consistent with conservation of the total entropy and not the entropy
conservation for each sector separately. An analysis of such a formalism was
given in ref. [3] for the case of one hidden sector and for the case of two
hidden sectors in ref. [4]. Further, an application of ref. [3] was made for
the explanation of EDGES anomaly [5] in ref. [6] (for related works see [7, 8,
9, 10]). It is to be noted that the extension to two hidden sectors brings in
new physics. Thus with one hidden sector, it is not possible to have dark
photon as dark matter because it is not possible to satisfy the twin
constraints that the dark photon simultaneously have a lifetime larger than
the lifetime of the universe and at the same give us significant amount of
dark matter consistent with the Planck data [11, 12]. We expect that further
extensions of the formalism involving more hidden sectors would also bring in
new physics which may have implications for particle physics and cosmology.

Thus in this work we extend the analysis to the case of $n$ number of hidden
sectors which are linked to the visible sector in a sequential manner.
Specifically, if we label the visible and the hidden sectors as
$S_{a},a=0,1,\cdots,n$ where $S_{0}$ stands for the visible sector and
$S_{a},a=1,\cdots,n$ stand for the hidden sectors, then we assume that the
sectors have only next neighbor couplings. That is to say that a sector
$S_{i}$ has couplings only with $S_{i-1}$ and $S_{i+1}$. In this setup, we
investigate the correlated evolution of temperatures in each sector relative
to a given reference temperature. While such a reference temperature can be
chosen arbitrarily, in this analysis we assume the temperature of the visible
sector to be the reference temperature or the ‘clock’. Specifically, we obtain
a set of $n$ coupled Boltzmann equations for the temperature evolution
functions $\xi_{a}=T_{a}/T,a=1,2,\cdots,n$, where $T$ is the temperature of
the visible sector and $T_{a}$ ($a=1,\cdots,n$) is the temperature of hidden
sector $a$. An analysis is also given of how very feeble couplings can be
obtained by a sequential couplings of the hidden sectors to the visible
sector. For some early works related to the evolution of hidden and visible
sectors see, e.g., [13, 14].

The outline of the rest of the paper is as follows: In section 2, we discuss
the model with $n$ hidden sectors, and in section 3 we apply this
generalization to $n$ hidden sectors to the Stueckelberg extension of the SM.
In section 4 we discuss the special case of three hidden sectors for a simple
model consisting of a dark Dirac fermion and dark scalar particles. In section
5 we study the thermal effect of hidden sectors on experimentally measured
observables such as the spin-independent proton-DM scattering cross-section
and on dark radiation given by $\Delta N_{\rm eff}$ for extra light relics.
Conclusions are given in section 6. A detailed derivation of the evolution
equations for visible and hidden sector temperatures is given in Appendix A
and the collision and source terms used in the Boltzmann equations are listed
in Appendix B.

##  2 A model with $n$ hidden sectors

The model consists of $n$ number of hidden sectors linked in a chain with the
first hidden sector connected to the visible sector as shown in Fig. 1.

Figure 1: The visible sector linked to a chain of $n$ number of hidden
sectors.

The above is a generalization of the cases considered previously for $n=1$ and
$n=2$. As mentioned earlier, we choose the temperature of the visible sector
as the clock. The set of Boltzmann equations for the energy density
$\rho_{\alpha}$ including energy exchange between the sectors is
$d\rho_{\alpha}/dt+3H(p_{\alpha}+\rho_{\alpha})=j_{\alpha}$
($\alpha=0,1,2,\cdots,n$). For the range of masses for the dark sector
particles we consider in this analysis, it is safe to assume radiation
domination333Early matter domination is possible, e.g., in the presence of
heavy long-lived particles and rapidly oscillating fields (such as the
inflaton field). In string theory, oscillating moduli fields displaced from
the minimum of their potential can dominate the energy density in the early
universe. None of these situations are relevant in our analysis. The light
spectrum of our dark species justify the radiation domination assumption.
down to $T\sim 10^{-5}$ GeV and so $p_{\alpha}=1/3\rho_{\alpha}$. Thus, we get

$\displaystyle\frac{d\rho_{\alpha}}{dt}+4H\rho_{\alpha}$
$\displaystyle=j_{\alpha},~{}~{}\alpha=0,1,2,\dots,n\,,$ (2.1)

where $\alpha=0$ refers to the visible sector and where $j_{\alpha}$ encodes
in it all the possible processes exchanging energy between neighboring
sectors. As discussed above, in order to properly describe the thermal
evolution of the universe with visible and hidden sectors in different heat
baths, we need to to have evolution equations for the $d\xi_{\alpha}/dT$.
However, to compute $d\xi_{\alpha}/dT,\alpha=1,2,\cdots,n$ we first derive
equations for $d\rho_{\alpha}/dT~{}(\alpha=1\cdots n)$ in terms of
$d\rho_{v}/dT$. We proceed by considering the Boltzmann equation of the total
energy density $\rho$, i.e., $d\rho/dt+4H\rho=0$. Knowing that
$d\rho/dt=(d\rho/dT)(dT/dt)$, we have

$\displaystyle\frac{dT}{dt}$ $\displaystyle=-\frac{4H\rho}{d\rho/dT}.$ (2.2)

Next we look at the Boltzmann equations for the individual energy densities
$d\rho_{\alpha}/dt+4H\rho_{\alpha}=j_{\alpha}$ which after using Eq. (2.2)
will immediately give us

$\displaystyle\frac{d\rho_{\alpha}}{dT}=\frac{(4H\rho_{\alpha}-j_{\alpha})}{4H\rho}\frac{d\rho}{dT}.$
(2.3)

Next we define

$\displaystyle\sigma_{v}$ $\displaystyle=\rho-\rho_{v},$ (2.4)
$\displaystyle\sigma_{\alpha}$
$\displaystyle=\rho-\rho_{\alpha}=\rho_{v}+\sum_{\beta\neq\alpha}\rho_{\beta}\quad(\alpha,\beta=1\cdots
n),$

where the total energy density is

$\rho=\rho_{v}+\sum_{\alpha=1}^{n}\rho_{\alpha}\,.$ (2.5)

To reduce clutter, we define
$K_{\alpha}\equiv(4H\rho_{\alpha}-j_{\alpha})/4H\rho$ and so Eq. (2.3) becomes

$\displaystyle\frac{d\rho_{\alpha}}{dT}$
$\displaystyle=K_{\alpha}\frac{d\rho}{dT}=K_{\alpha}\left(\frac{d\rho_{v}}{dT}+\frac{d\rho_{\alpha}}{dT}+\sum_{\beta\neq\alpha}\frac{d\rho_{\beta}}{dT}\right).$
(2.6)

Rearranging, we get

$\displaystyle\frac{d\rho_{\alpha}}{dT}-C_{\alpha}\sum_{\beta\neq\alpha}\frac{d\rho_{\beta}}{dT}=C_{\alpha}\,\frac{d\rho_{v}}{dT},$
(2.7)

where we have defined

$C_{\alpha}\equiv\frac{K_{\alpha}}{1-K_{\alpha}}=\frac{4H\rho_{\alpha}-j_{\alpha}}{4H\rho-4H\rho_{\alpha}+j_{\alpha}}=\frac{4H\rho_{\alpha}-j_{\alpha}}{4H\sigma_{\alpha}+j_{\alpha}}.$
(2.8)

Further, we may write out the coupled equation for $d\rho_{\alpha}/dT$ as
follows

$\displaystyle d\rho_{1}/dT-C_{1}d\rho_{2}/dT-C_{1}d\rho_{3}/dT\cdots-
C_{1}d\rho_{n}/dT=C_{1}d\rho_{v}/dT,$ $\displaystyle-
C_{2}d\rho_{1}/dT+d\rho_{2}/dT-C_{2}d\rho_{3}/dT\cdots-
C_{2}d\rho_{n}/dT=C_{2}d\rho_{v}/dT,$ $\displaystyle\cdots\cdots$
$\displaystyle-C_{n}d\rho_{1}/dT-C_{n}d\rho_{2}/dT-
C_{n}d\rho_{3}/dT\cdots+d\rho_{n}/dT=C_{n}d\rho_{v}/dT.$ (2.9)

We can write this in a matrix form so that

$\displaystyle\left(\begin{matrix}1&-C_{1}&-C_{1}&\cdots&-C_{1}&-C_{1}\\\
-C_{2}&1&-C_{2}&\cdots&-C_{2}&-C_{2}\\\
-C_{3}&-C_{3}&1&\cdots&-C_{3}&-C_{3}\\\
\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\\
-C_{n-1}&-C_{n-1}&-C_{n-1}&\cdots&1&-C_{n-1}\\\
-C_{n}&-C_{n}&-C_{n}&\cdots&-C_{n}&1\end{matrix}\right)\left(\begin{matrix}d\rho_{1}/dT\\\
d\rho_{2}/dT\\\ d\rho_{3}/dT\\\ \cdots\\\ d\rho_{n-1}/dT\\\
d\rho_{n}/dT\end{matrix}\right)=\left(\begin{matrix}C_{1}\\\ C_{2}\\\ C_{3}\\\
\cdots\\\ C_{n-1}\\\ C_{n}\end{matrix}\right)d\rho_{v}/dT.$ (2.10)

This equation can be written in a compact form so that

$\displaystyle\sum_{\beta=1}^{n}{\cal
C}_{\alpha\beta}\frac{d\rho_{\beta}}{dT}=C_{\alpha}\frac{d\rho_{v}}{dT},$
(2.11)

where the matrix ${\cal C}$ is the square matrix on the left hand side of Eq.
(2.10). Thus $d\rho_{\alpha}/dT$ can be written as

$\displaystyle\frac{d\rho_{\alpha}}{dT}=\sum_{j=1}^{n}({\cal
C}^{-1})_{\alpha\beta}C_{\beta}\frac{d\rho_{v}}{dT}.$ (2.12)

Note that one may also write

$\displaystyle\rho_{\alpha}^{\prime}$
$\displaystyle=\frac{d\rho_{\alpha}}{dT}=P_{\alpha}+Q_{\alpha}\xi_{\alpha}^{\prime},~{}~{}~{}\alpha=1,2,\cdots,n-1,n,$
(2.13)

where $\xi^{\prime}_{\alpha}\equiv d\xi_{\alpha}/dT$ and where $P_{\alpha}$
and $Q_{\alpha}$ are defined so that

$\displaystyle
P_{\alpha}=\xi_{\alpha}\frac{dP_{\alpha}}{dT_{\alpha}},~{}~{}Q_{\alpha}=T\frac{d\rho_{\alpha}}{dT_{\alpha}}.$
(2.14)

Thus we have an equation for $d\xi_{\alpha}/dT$ which takes the form

$\displaystyle\frac{d\xi_{\alpha}}{dT}=-\frac{P_{\alpha}}{Q_{\alpha}}+\sum_{\beta=1}^{n}(C^{-1})_{\alpha\beta}C_{\beta}\frac{\rho_{v}^{\prime}}{Q_{\alpha}},~{}~{}~{}~{}\alpha=1,2,\cdots,n-1,n.$
(2.15)

Eqs. (2.15) give us $n$ differential equations for the evolution functions
$d\xi_{\alpha}/dT$. These have to be solved along with the Boltzmann equations
governing the number density evolution of the hidden sector particles. This
will allow us to determine the relic densities of all stable species and
describe the thermal evolution of this coupled system.

For the special case of three hidden sectors, the matrix ${\cal{C}}$ is given
by

$\displaystyle{\cal C}=\left(\begin{matrix}1&-C_{1}&-C_{1}\\\
-C_{2}&1&-C_{2}\\\ -C_{3}&-C_{3}&1\end{matrix}\right).$ (2.16)

Therefore we can write Eq. (2.12) in a matrix form so that

$\displaystyle\left(\begin{matrix}d\rho_{1}/dT\\\ d\rho_{2}/dT\\\
d\rho_{3}/dT\end{matrix}\right)=\frac{1}{D}\left(\begin{matrix}(1-C_{2}C_{3})&C_{1}(1+C_{3})&C_{1}(1+C_{2})\\\
C_{2}(1+C_{3})&(1-C_{1}C_{3})&C_{2}(1+C_{1})\\\
C_{3}(1+C_{2})&C_{3}(1+C_{1})&(1-C_{1}C_{2})\end{matrix}\right)\left(\begin{matrix}C_{1}\\\
C_{2}\\\ C_{3}\end{matrix}\right)d\rho_{v}/dT,$ (2.17)

where $D$ is given by

$\displaystyle D$
$\displaystyle=1-C_{1}C_{2}-C_{1}C_{3}-C_{2}C_{3}-2C_{1}C_{2}C_{3}.$ (2.18)

Since $\rho_{\alpha}^{\prime}=P_{\alpha}+Q_{\alpha}\xi_{\alpha}^{\prime}$, we
have

$\displaystyle\xi_{\alpha}^{\prime}\equiv\frac{d\xi_{\alpha}}{dT}=-\frac{P_{\alpha}}{Q_{\alpha}}+\frac{1}{Q_{\alpha}}\rho^{\prime}_{\alpha}.$
(2.19)

Therefore the individual expressions can be easily obtained as

$\displaystyle\frac{d\xi_{1}}{dT}$
$\displaystyle=-\frac{P_{1}}{Q_{1}}+\frac{1}{Q_{1}D}\left[C_{1}+C_{1}C_{2}+C_{1}C_{3}+C_{1}C_{2}C_{3})\right]d\rho_{v}/dT,$
(2.20) $\displaystyle\frac{d\xi_{2}}{dT}$
$\displaystyle=-\frac{P_{2}}{Q_{2}}+\frac{1}{Q_{2}D}\left[C_{2}+C_{1}C_{2}+C_{2}C_{3}+C_{1}C_{2}C_{3})\right]d\rho_{v}/dT,$
(2.21) $\displaystyle\frac{d\xi_{3}}{dT}$
$\displaystyle=-\frac{P_{3}}{Q_{3}}+\frac{1}{Q_{3}D}\left[C_{3}+C_{2}C_{3}+C_{1}C_{3}+C_{1}C_{2}C_{3})\right]d\rho_{v}/dT.$
(2.22)

These equations give the evolution of the temperature ratios
$\xi_{\alpha}=T_{\alpha}/T,\alpha=1,2,3$.

Before proceeding further, we note here that a chain of $n$ hidden sectors
analogous to the one of Fig. 1 appears in well known scenarios of physics
beyond the Standard Model such as in moose/quiver gauge theories (see, e.g.,
[15, 16, 17, 18]). Thus an analysis of the type above may find application in
a broader class of particle physics models beyond the Standard Model.

##  3 Theory application: The Stueckelberg extension of the Standard Model

The Stueckelberg extension of the SM with an extra hidden $U(1)_{X}$ gauge
group offers the possibility of a hidden sector feebly coupled to the visible
sector. The gauged $U(1)_{X}$ gives rise to a dark photon which, under certain
conditions, can be a dark matter candidate [19, 4]. In this section we give
the analytic expressions for a Stueckelberg extension of the SM to $n$ hidden
sectors which the above formalism can easily be applied to.

###  3.1 Generation of very feeble couplings from kinetic mixings via a chain
of hidden sectors

As discussed in section 2, we consider the following type of couplings: the
visible sector couples to the hidden sector 1 $(H_{1})$, the hidden sector 1
couples to hidden sector 2 $(H_{2})$ and so on, up to hidden sector ($n-1$)
which couples to the hidden sector $n$ $(H_{n})$. The couplings between the
sectors could be via kinetic mixing [20], via Stueckelberg mass mixing [21],
via both kinetic mixing and Stueckelberg mass mixing [22] and contain matter
in the hidden sectors [21, 23] (for recent works based on Stueckelberg
extensions of the Standard Model see also [24, 25]). Specifically we assume
that each of the $n$ number of hidden sectors possess a $U(1)_{X_{I}}$ gauge
invariance with a gauge field $A_{I}^{\mu},I=1,\cdots,n$. Including the gauge
field $B_{\mu}$ for the hypercharge $U(1)_{Y}$, we have $n+1$ $U(1)$ gauge
fields which we collectively denote by $V_{i}^{\mu},i=1,\cdots,n+1$ where we
define the elements of $V_{i}$ as follows:

$\displaystyle(V)=(V_{1}=B,V_{2}=A_{1},\cdots,V_{n+1}=A_{n}),$ (3.1)

where $B_{\mu}$ is the hypercharge field of the visible sector and
$A_{I\mu},~{}I=1,\cdots,n$ are the gauge fields for the hidden sector
$U(1)_{X_{I}},~{}I=1,\cdots,n$. We assume the gauge fields have kinetic
mixings and Stueckelberg mass mixings and such mixings to be sequential and
thus the Lagrangian for the gauge fields is assumed to be of the form

$\displaystyle\mathcal{L}_{\rm km}=$
$\displaystyle-\frac{1}{4}\left[\sum_{i=1}^{n+1}G_{i\mu\nu}G_{i}^{\mu\nu}+2\sum_{i=1}^{n}\delta_{i}G_{i\mu\nu}G_{i+1}^{\mu\nu}\right]$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{n}(M_{i}V_{i\mu}+\mu_{i}V_{i+1\mu}+\partial_{\mu}\sigma_{i})^{2}.$
(3.2)

Note that there are $n$ axion fields $\sigma_{i},~{}i=1,\cdots,n$ which are
absorbed by the $n$ gauge fields of the hidden sectors to become massive.

We note that the fields in the Lagrangian above are not in their canonical
form. Further, this Lagrangian must be coupled with the Standard Model Higgs
field in the electroweak symmetry breaking. After the electroweak symmetry
breaking, one must make the appropriate transformations to put the Lagrangian
in the canonical form. We give now a brief discussion of this procedure. In
the analysis below we will include the third component of the gauge fields for
$SU(2)_{L}({\cal{A}}_{3})$ of the SM since that fields along with the
hypercharge gauge field enter in giving mass to the $W$ and $Z$ bosons in the
electroweak symmetry breaking. The transformation that puts the kinetic energy
in the canonical form is given by

$\displaystyle{\cal{V}}_{\mu}=K{\cal{V}}^{\prime}_{\mu}\,,$ (3.3)

where our notation is the following444Note that while $V_{i\mu}$ has $i$
taking values $i=1,\dots,n+1$, ${\cal{V}}_{a}$ has $a$ taking values
$a=0,1,\cdots,n+1$.

$\displaystyle{\cal{V}}_{\mu}^{T}$
$\displaystyle=(\\{{\cal{V}}_{a}\\})=({\cal{A}}_{3\mu},V_{1\mu},V_{2\mu},V_{3\mu},\cdots,V_{n\mu},V_{n+1,\mu}),a=0,\cdots,n+1,$
$\displaystyle{\cal{V}}_{\mu}^{{}^{\prime}T}$
$\displaystyle=(\\{{\cal{V}}^{\prime}_{a}\\})=({\cal{A}}_{3\mu},V^{\prime}_{1\mu},V^{\prime}_{2\mu},V^{\prime}_{3\mu},\cdots,V^{\prime}_{n\mu},V^{\prime}_{n+1,\mu}),a=0,\cdots,n+1,$

and where $K$ is an $(n+2)\times(n+2)$ dimensional matrix defined as

$\displaystyle K$
$\displaystyle=\left(\begin{matrix}1&0&0&0&0&\cdot&\cdot&0\\\
0&1&-s_{1}&s_{1}s_{2}&-s_{1}s_{2}s_{3}&\cdot&\cdot&(-1)^{n}s_{1}s_{2}\cdots
s_{n}\\\
0&0&c_{1}&-c_{1}s_{2}&c_{1}s_{2}s_{3}&\cdot&\cdot&(-1)^{(n-1)}c_{1}s_{2}\cdots
s_{n}\\\ 0&0&0&c_{2}&-c_{2}s_{3}&\cdot&\cdot&(-1)^{(n-2)}c_{2}s_{3}\cdots
s_{n}\\\ \cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
0&0&0&0&0&0&c_{(n-1)}&-c_{(n-1)}s_{n}\\\
0&0&0&0&0&0&0&c_{n}\end{matrix}\right),$ (3.4)

where in general $c_{k}$ and $s_{k}$ (for $k=1,\dots,n$) are given by

$\displaystyle c_{k}=$
$\displaystyle\frac{\sqrt{1-\sum_{i=1}^{k-1}}\delta_{i}^{2}}{\sqrt{1-\sum_{i=1}^{k}}\delta_{i}^{2}},~{}~{}s_{k}=\frac{\delta_{k}}{\sqrt{1-\sum_{i=1}^{k}}\delta_{i}^{2}}.$
(3.5)

Note that $c_{k},s_{k}$ satisfy the relation $c_{k}^{2}-s_{k}^{2}=1$. In the
transformed coordinates the kinetic energy part of the Lagrangian takes the
form
$\mathcal{L}=-\frac{1}{4}{\cal{V}}^{{}^{\prime}T}_{\mu\nu}{\cal{V}}^{{}^{\prime}\mu\nu}.$
Next, we focus on the mass terms. In addition to the Stueckelberg mass terms
involving the visible and the hidden sectors, we have mass terms for
$A_{0},{{\cal{A}}_{3}}$ from the Higgs mechanism which have the form

$\displaystyle\mathcal{L}^{m}_{\text{SM}}=-\frac{1}{2}\left(\frac{1}{4}g_{Y}^{2}v^{2}B^{\mu}B_{\mu}+\frac{1}{4}g_{2}^{2}v^{2}{\cal{A}}^{\mu}_{3}{\cal{A}}_{3\mu}+\frac{1}{2}g_{1}g_{2}v^{2}B^{\mu}{\cal{A}}^{\mu}_{3}\right),$
(3.6)

where $v=(\sqrt{2}G_{F})^{-1/2}$ and where $G_{F}$ is the Fermi constant. Thus
we write the total mass term for the gauge fields so that

$\displaystyle\mathcal{L}_{m}=-\frac{1}{2}{\cal{V}}_{\mu}^{T}M^{2}{\cal{V}}^{\mu},$
(3.7)

where $M^{2}$ includes the mass terms both from the Higgs mechanism and from
the Stueckelberg mass growth mechanism. After the transformation that
diagonalizes the kinetic energy, the mass term is given by

$\displaystyle\mathcal{L}_{m}=-\frac{1}{2}{\cal{V}}_{\mu}^{\prime T}M^{\prime
2}{\cal{V}}^{\prime\mu},$ (3.8)

where ${\cal{V}}^{{}^{\prime}\mu}$ is defined in Eq. (3.3) and $M^{\prime
2}=K^{T}M^{2}K$. Since $M^{\prime 2}$ is a symmetric matrix it can be
diagonalized by an orthogonal transformation $R$ so that $R^{T}R=I$. Thus the
transformation ${\cal{V}}^{\prime}=R{\cal{E}}$, or ${\cal{V}}=KR{\cal{E}}$
diagonalizes the mass matrix.

$\displaystyle M_{D}^{2}$ $\displaystyle={\mathcal{O}}^{T}M^{2}{\mathcal{O}},$
$\displaystyle{\cal{V}}$ $\displaystyle={\mathcal{O}}{\cal{E}},$ (3.9)

where ${\cal{O}}=KR$ and both the kinetic and the mass terms are in a
canonically diagonal basis and we have

$\displaystyle\mathcal{L}=-\frac{1}{4}{{\cal{E}}}^{T}_{\mu\nu}{{\cal{E}}}^{\mu\nu}-\frac{1}{2}{{\cal{E}}}^{T}_{\mu}M_{D}^{2}{\cal{E}}^{\mu}.$
(3.10)

We note that the subcase $n=2$ was discussed in [26]. We arrange the vector
mass eigenstates ${\cal{E}}_{\mu}$ so that
${\cal{E}}^{T}=(Z,A_{\gamma},A_{\gamma_{2}},\dots,A_{\gamma_{3}},\cdots,A_{\gamma_{n}}$).
Thus the mass eigenstates consist of the $Z$ boson, the photon, and $n$ number
of dark photons in the hidden sector.

###  3.2 Neutral currents of the coupled visible and hidden sectors

Since the mixings involve the neutral gauge bosons of the visible and the
hidden sectors, we begin by displaying the interactions of these gauge bosons
in their respective sectors. Thus the neutral current interactions of the
visible and the hidden sectors together are

$\displaystyle\mathcal{L}_{\rm
NC}=g_{2}J_{2}^{\mu}{\cal{A}}_{3\mu}+g_{Y}J_{Y}^{\mu}B_{\mu}+\sum_{i=1}^{n}g_{Xi}J^{\mu}_{Xi}A_{i\mu}.$
(3.11)

In the equation above, matter charged under a hidden sector gauge field
$A_{i\mu}$ will reside in the source term $J^{\mu}_{Xi}$. Such matter may
consist of dark fermions and dark complex scalar fields charged under the
$U(1)_{Xi}$. Thus if we have Dirac fermions $D_{Xi}$ charged under $U(1)_{Xi}$
we will have $J^{\mu}_{Xi}=\bar{D}_{Xi}\gamma^{\mu}D_{Xi}$.

To write the neutral currents in a compact notation in the ${\cal{V}}$ basis,
we define

$\displaystyle({\cal{V}}_{a}):$
$\displaystyle({\cal{A}}_{3\mu},B_{\mu},A_{1\mu},\cdots,A_{n\mu}),$
$\displaystyle({\cal{G}}_{a}):$
$\displaystyle{\cal{G}}_{0}=g_{2},~{}{\cal{G}}_{1}=g_{Y},~{}{\cal{G}}_{3}=g_{X_{1}},\cdots,{\cal{G}}_{n+1}=g_{X_{n}}~{}a=0,\cdots,n+1,$
$\displaystyle({\cal{J}}_{a}):$
$\displaystyle{\cal{J}}^{\mu}_{0}=J^{\mu}_{2},{\cal{J}}^{\mu}_{1}=J^{\mu}_{Y},\cdots,{\cal{J}}^{\mu}_{i+1}=J^{\mu}_{X_{i}},~{}a=0,\cdots,n+1$
$\displaystyle({\cal{E}}_{a}):$
$\displaystyle{\cal{E}}^{\mu}_{0}=Z^{\mu},~{}{\cal{E}}^{\mu}_{1}=A^{\mu}_{\gamma},~{}{\cal{E}}^{\mu}_{2}=A^{\mu}_{\gamma_{2}},\cdots,,~{}{\cal{E}}^{\mu}_{n+1}=A^{\mu}_{\tilde{\gamma}_{n}},a=0,\cdots,n+1,$
(3.12)

where $Z^{\mu}$ is the $Z$-boson field, $A^{\mu}_{\gamma}$ is the photon
field, and $A^{\mu}_{\tilde{\gamma}_{i}},\cdots,A^{\mu}_{\tilde{\gamma}_{n}}$
are the fields for the dark photons. We now rewrite the neutral current
interaction so that

$\displaystyle\mathcal{L}_{\rm
NC}=\sum_{a=0}^{n+1}{\cal{G}}_{a}{\cal{J}}^{\mu}_{a}{\cal{V}}_{a}.$ (3.13)

Next we write $\mathcal{L}_{\rm NC}$ in terms of the mass eigenstates
${\cal{E}}$ so that

$\displaystyle\mathcal{L}_{\rm
NC}=\sum_{a=0}^{n+1}{\cal{G}}_{a}{\cal{J}}^{\mu}_{a}{\cal{O}}_{ab}{\cal{E}}_{b\mu}.$
(3.14)

We can now display the coupling of the $Z$ boson and the photon and for the
dark photon in an explicit manner. Thus we can decompose the couplings into
several parts so that

$\displaystyle\mathcal{L}_{\rm NC}=\mathcal{L}^{\rm NC}_{\rm
SM}+\mathcal{L}_{Z\gamma
D}+\mathcal{L}_{\text{SM}\tilde{\gamma}}+\mathcal{L}_{D\tilde{\gamma}},$
(3.15)

where

$\displaystyle\mathcal{L}^{\rm NC}_{\rm SM}$
$\displaystyle=\left(g_{2}J_{2}^{\mu}{\cal{O}}_{00}+g_{Y}J^{\mu}_{Y}{\cal{O}}_{10}\right)Z_{\mu}+\left(g_{2}J_{2}^{\mu}{\cal{O}}_{01}+g_{Y}J^{\mu}_{Y}{\cal{O}}_{11}\right)A^{\gamma}_{\mu},$
(3.16) $\displaystyle\mathcal{L}_{Z{\gamma}D}$
$\displaystyle=\sum_{a=2}^{n+1}g_{Xa}J^{\mu}_{Xa}\mathcal{O}_{a0}Z_{\mu}+\sum_{a=2}^{n+1}g_{Xa}J^{\mu}_{Xa}\mathcal{O}_{a1}A^{\gamma}_{\mu},$
$\displaystyle\mathcal{L}_{\text{SM}\tilde{\gamma}}$
$\displaystyle=\sum_{i=1}^{n}\left(g_{2}J_{2}^{\mu}{\cal{O}}_{0(i+1)}+g_{Y}J_{Y}{\cal{O}}_{1(i+1)}\right)A^{\tilde{\gamma}_{i}}_{\mu},$
$\displaystyle\mathcal{L}_{D\tilde{\gamma}}$
$\displaystyle=\sum_{i,j=1}^{n}g_{X_{i}}J_{X_{i}}\mathcal{O}_{(j+1)(i+1)}A^{\tilde{\gamma}_{j}}_{\mu}.$

Here $\mathcal{L}^{\rm NC}_{\rm SM}$ is the modified neutral current of the
Standard Model, $\mathcal{L}_{Z{\gamma}D}$ is the coupling of the SM gauge
bosons $Z,A_{\gamma}$ with the dark sector,
$\mathcal{L}_{\text{SM}\tilde{\gamma}}$ is the coupling of the dark photons
with the SM sector and $\mathcal{L}_{D\tilde{\gamma}}$ is the coupling of the
dark photons with the dark sources. The various parts of the couplings of Eqs.
(3.16) have the following meaning. As stated $\mathcal{L}_{\rm SM}$ is the
modified neutral current of the Standard Model. $\mathcal{L}_{Z{\gamma}D}$ and
$\mathcal{L}_{\text{SM}\tilde{\gamma}}$ enter in the freeze-in production of
the hidden sector particles, $\mathcal{L}_{\text{SM}\tilde{\gamma}}$ enters in
the decay of the dark photons, and $\mathcal{L}_{D\tilde{\gamma}}$ controls
freeze-out in the dark sector, while $\mathcal{L}_{Z{\gamma}D}$ enters in the
analysis of kinetic equilibrium of the dark sources after chemical decoupling.

##  4 Numerical application: The visible sector coupled to three hidden
sectors

###  4.1 A simple model

In this section we will use the general expression of Eq. (2.15) to determine
the temperature evolution of a chain of hidden sectors coupled to the visible
sector as described by Fig. 1. The hidden sectors will eventually get
populated by dark particle species and so the evolution of the temperature is
coupled with the evolution of the particle number densities. In refs. [3, 4]
the cases of one and two hidden sectors have been considered in the context of
a Stueckelberg $U(1)$ extension of the SM. Here, we add a third hidden sector
but consider a simpler model to illustrate the effect of three hidden sectors
on dark matter phenomenology and on cosmology. Further, we assume that
particles in the different sectors communicate via a scalar portal which will
become evident next.

Our set up of the matter content of the three hidden sectors is as follows:
The first hidden sector has a Dirac fermion $D$ (which serves as one of the
possible dark matter candidates) and a pseudo-scalar $\phi_{1}$ (mediator).
The field $\phi_{2}$ is a scalar while $\phi_{3}$ is a pseudo-scalar and are
not particularly associated with any specific hidden sector. The model
presented below is phenomenological and its UV completion remains to be
explored. The Lagrangian of the model we use for the rest of the analysis is

$\displaystyle\mathcal{L}$ $\displaystyle=\mathcal{L}_{\rm
SM}+\Delta\mathcal{L},$ (4.1)

where $\mathcal{L}_{\rm SM}$ is the SM Lagrangian and $\Delta\mathcal{L}$ is
the Lagrangian for the hidden sector including portal interactions between the
visible and the hidden sectors. We note here in passing that we are
considering the hidden sector to have no gauged $U(1)$ symmetry and for that
reason there is no kinetic energy for the hidden sector gauge field or its
interactions included in $\Delta\mathcal{L}$ below. Thus $\Delta\mathcal{L}$
is given by

$\displaystyle\Delta\mathcal{L}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{3}(\partial_{\mu}\phi_{i}\partial^{\mu}\phi_{i}-m_{i}^{2}\phi^{2}_{i})+\bar{D}(i\gamma^{\mu}\partial_{\mu}-m_{D})D$
$\displaystyle~{}~{}~{}-y_{q}\phi_{1}\bar{q}\gamma^{5}q-y_{\ell}\phi_{1}\bar{\ell}\gamma^{5}\ell-
y_{D}\phi_{1}\bar{D}\gamma^{5}D$
$\displaystyle~{}~{}~{}-\frac{\kappa_{1}}{3}\phi_{1}^{2}\phi_{2}-\frac{\kappa_{2}}{3}\phi_{2}\phi_{3}^{2}-\frac{1}{4}\sum_{i=1}^{3}\lambda_{i}\phi_{i}^{4}-\frac{1}{4}\sum_{i}\sum_{j>i}\lambda_{ij}\phi_{i}^{2}\phi_{j}^{2}\,,$
(4.2)

where the dark fermion $D$ resides in the hidden sector 1 and couples only to
the field $\phi_{1}$ in the hidden sector 1. The fields $\phi_{1}$ and
$\phi_{3}$ as pseudoscalars and $\phi_{2}$ as CP-even scalar are chosen purely
on phenomenological grounds which we use to illustrate the main features of
the visible sector coupling to three hidden sectors. For simplicity, we assume
that $\phi_{1}$ only couples to the first generation of quarks and leptons via
the Yukawa couplings $y_{q}$ and $y_{\ell}$ and to the $D$ fermions via
$y_{D}$. Thus in this model, the field $\phi_{1}$ is the portal connecting the
visible and hidden sectors. Further, for simplicity, we will assume equality
of the portal field coupling to the Standard Model quarks and leptons, i.e.,
$y_{q}=y_{\ell}=y_{f}$. The size of the Yukawas could be very small. Thus, for
example, a higher dimensional operator such as
$\frac{H^{\dagger}H}{\Lambda^{2}}\bar{q}q$ after electroweak symmetry breaking
would give $y_{q}\sim 10^{-9}$ for $\Lambda={\cal O}(10^{6})$ GeV. A string
scale such as this can arise in Type IIB string theory. The quartic couplings
$\lambda_{i}$ do not play a role in the temperature and number density
evolution while we assume the couplings $\lambda_{ij}$ are small compared to
$\kappa_{1}$ and $\kappa_{2}$.

###  4.2 Evolution equations

In the set up described above, only the first hidden sector possesses direct
couplings to SM (the visible sector) through the Yukawa couplings $y_{q}$ and
$y_{\ell}$ while the second hidden sector connects to the first and to the
third hidden sectors. We assume that the visible sector is held at a
temperature $T$ and the three hidden sectors at temperatures $T_{1}$, $T_{2}$
and $T_{3}$ where one prescribes the initial values of the temperatures
$T_{i}\ll T$ ($i=1,2,3)$. The Boltzmann equations for the number and total
energy densities are given by

$\displaystyle\frac{dn_{r}}{dt}+3Hn_{r}=\mathcal{Q}_{r}\,,$ (4.3)
$\displaystyle\frac{d\rho}{dt}+3H(\rho+p)=0,$ (4.4)

with $r\in\\{D,\phi_{1},\phi_{2},\phi_{3}\\}$. In the above, $\rho$ is the
total energy density, $p$ is the total pressure, $H$ is the Hubble parameter
and $\mathcal{Q}_{r}$ are the collision terms which include all possible
particle number-changing processes. They are listed in Appendix B.

After defining the co-moving number density or yield, $Y=n/\mathbb{s}$, where
$\mathbb{s}$ is the total entropy density and using Eqs. (4.3) and (4.4), we
write

$\displaystyle\frac{dY_{r}}{dT}=$
$\displaystyle-\mathbb{s}\frac{d\rho/dT}{4H\rho}\mathcal{Q}_{r},$ (4.5)

where $H$ is related to the energy density by the Friedmann equation

$H^{2}=\frac{8\pi G_{N}}{3}\rho.$ (4.6)

One can define the individual energy and entropy densities in terms of
effective number of degrees of freedom so that

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

We note that $g_{\rm eff}$ and $h_{\rm eff}$ of the hidden sectors are defined
as integrals which depend on the mass of the relevant particle and
temperature. In appendix A, we exhibit the $m/T$ dependence of the integrals.
In the above, we need to compute $T_{i}$ so that $T_{i}=\xi_{i}T$. Thus we
need to write out the evolution equations for $\xi_{i}$ in an explicit form.
In this case we have three evolution functions
$\xi_{1}=T_{1}/T,\xi_{2}=T_{2}/T$ and $\xi_{3}=T_{3}/T$ and using Eq. (2.15)
we find that they obey the following equations

$\displaystyle\frac{d\xi_{1}}{dT}$
$\displaystyle=-\frac{P_{1}}{Q_{1}}+\frac{1}{Q_{1}D}\left[C_{1}+C_{1}C_{2}+C_{1}C_{3}+C_{1}C_{2}C_{3})\right]d\rho_{v}/dT,$
(4.9) $\displaystyle\frac{d\xi_{2}}{dT}$
$\displaystyle=-\frac{P_{2}}{Q_{2}}+\frac{1}{Q_{2}D}\left[C_{2}+C_{1}C_{2}+C_{2}C_{3}+C_{1}C_{2}C_{3})\right]d\rho_{v}/dT,$
(4.10) $\displaystyle\frac{d\xi_{3}}{dT}$
$\displaystyle=-\frac{P_{3}}{Q_{3}}+\frac{1}{Q_{3}D}\left[C_{3}+C_{2}C_{3}+C_{1}C_{3}+C_{1}C_{2}C_{3})\right]d\rho_{v}/dT,$
(4.11)

where $D$ is given by

$\displaystyle D$
$\displaystyle=1-C_{1}C_{2}-C_{1}C_{3}-C_{2}C_{3}-2C_{1}C_{2}C_{3}.$ (4.12)

For the case of one hidden sector ($n=1$), we have $C_{2}=C_{3}=0$, $D=1$ and
$\xi_{2}$ and $\xi_{3}$ are absent. In this case Eq. (4.9) reads

$\displaystyle\frac{d\xi_{1}}{dT}=-\frac{P_{1}}{Q_{1}}+C_{0}\frac{\rho_{v}^{\prime}}{Q_{1}},$
(4.13)

which agrees with the result of ref. [3]. Next we look at the $n=2$ case (two
hidden sectors). Here we set $C_{3}=0$ and $D=1-C_{1}C_{2}$, so that

$\displaystyle\frac{d\xi_{1}}{dT}$
$\displaystyle=-\frac{P_{1}}{Q_{1}}+\frac{C_{1}+C_{2}C_{1}}{1-C_{1}C_{2}}\frac{\rho_{v}^{\prime}}{Q_{1}}.$
(4.14)

Similarly we have

$\displaystyle\frac{d\xi_{2}}{dT}$
$\displaystyle=-\frac{P_{2}}{Q_{2}}+\frac{C_{1}C_{2}+C_{2}}{1-C_{1}C_{2}}\frac{\rho_{v}^{\prime}}{Q_{2}}.$
(4.15)

These agree with the analysis of refs. [4, 19]. We note here that
$C_{1},C_{2},C_{3}$ are given by

$\displaystyle C_{1}$
$\displaystyle=\frac{4H\rho_{1}-j_{1}}{4H\sigma_{1}+j_{v}},~{}~{}C_{2}=\frac{4H\rho_{2}-j_{2}}{4H\sigma_{2}+j_{2}},~{}~{}C_{3}=\frac{4H\rho_{3}-j_{3}}{4H\sigma_{3}+j_{3}},$
(4.16)

where the source terms admit the relation $j_{v}+j_{1}+j_{2}+j_{3}=0$. The
source terms $j_{1}$, $j_{2}$ and $j_{3}$ are listed in Appendix B. Carrying
out some further simplifications to Eqs. (4.9)$-$(4.11), we can write them as

$\displaystyle\frac{d\xi_{1}}{dT}$
$\displaystyle=-\frac{\xi_{1}}{T}+\left(\frac{4H\rho_{1}-j_{1}}{4H\rho_{v}+j_{1}+j_{2}+j_{3}}\right)\frac{d\rho_{v}/dT}{T\frac{d\rho_{1}}{dT_{1}}},$
(4.17) $\displaystyle\frac{d\xi_{2}}{dT}$
$\displaystyle=-\frac{\xi_{2}}{T}+\left(\frac{4H\rho_{2}-j_{2}}{4H\rho_{v}+j_{1}+j_{2}+j_{3}}\right)\frac{d\rho_{v}/dT}{T\frac{d\rho_{2}}{dT_{2}}},$
(4.18) $\displaystyle\frac{d\xi_{3}}{dT}$
$\displaystyle=-\frac{\xi_{3}}{T}+\left(\frac{4H\rho_{3}-j_{3}}{4H\rho_{v}+j_{1}+j_{2}+j_{3}}\right)\frac{d\rho_{v}/dT}{T\frac{d\rho_{3}}{dT_{3}}}.$
(4.19)

The factor $\frac{d\rho/dT}{4H\rho}$ appearing in Eq. (4.5) involves a
derivative of the total energy density with respect to the visible sector
temperature. We differentiate each of the energy densities with respect to $T$
and remembering that $\xi_{\alpha}=\xi_{\alpha}(T)$, then Eq. (4.5) can be
cast in the form

$\displaystyle\frac{dY_{r}}{dT}=$
$\displaystyle-\frac{\mathbb{s}}{H}\Big{(}\frac{d\rho_{v}/dT}{4\rho_{v}-j_{v}/H}\Big{)}\mathcal{Q}_{r}\,.$
(4.20)

It is to be noted from Eq. (4.20) that the source term $j_{v}$ now appears in
the Boltzmann equations for the particle yields. The formalism given above
provides us with a set of coupled Boltzmann equations in temperature and yield
which need to be numerically solved to get a full description of the thermal
evolution of the hidden sectors. We do this next.

###  4.3 Numerical analysis

In our analysis here we consider the case of self-interacting DM, where we
take $m_{\phi_{1}}/m_{D}\ll 1$ which will allow us to consider the most recent
results on the DM-nucleon spin-independent (SI) cross-section from PandaX-II
experiment [27] as an experimental constraint. In the following we give an
analysis of three type of mass hierarchies among the matter fields of the
hidden sectors which are: Case 1: $m_{\phi_{1}}>2m_{\phi_{2}}$,
$m_{\phi_{2}}=m_{\phi_{3}}$; Case 2: $m_{\phi_{2}}>2m_{\phi_{1}}$,
$m_{\phi_{2}}>2m_{\phi_{3}}$; Case 3: $m_{\phi_{2}}>2m_{\phi_{1}}$,
$m_{\phi_{3}}=0$. For each of these cases we compute the thermal evolution
given by $\xi_{i}=T_{i}/T$ for the three hidden sectors, and compute the
yields. Specifically the analysis indicates the fraction of the relic density
of the dark matter contributed by each of the dark sectors. We discuss each of
these cases in detail below.

### Case 1: $m_{\phi_{1}}>2m_{\phi_{2}}$ and $m_{\phi_{2}}=m_{\phi_{3}}$

In Fig. 2 we exhibit the particle yields (left panel) and the temperature
ratios $\xi_{i}$ (right panel) as a function of the visible sector
temperature. As the temperature $T$ drops, the hidden sectors are gradually
populated by the dark particle species. The process happens sequentially where
the first hidden sector gets populated first whereas the second and third
hidden sectors see an increase in the yield much later. This is because the
first hidden sector is the only one with direct couplings to the visible
sector whereas the other two hidden sectors rely on the first hidden sector to
get populated.

Figure 2: The case of three initially cold dark sectors. Left panel: evolution
of the yields of the dark fermion $D$ and of the spin zero particles of the
dark sector versus the visible sector temperature. Right panel: evolution of
the temperature ratios $\xi_{i}$. The input parameters are: $m_{D}=1.0$ GeV,
$m_{\phi_{1}}=50$ MeV, $m_{\phi_{2}}=m_{\phi_{3}}=1$ MeV, $y_{f}=1.3\times
10^{-9}$, $y_{D}=0.053$ and $\kappa_{1}=\kappa_{2}=10^{-6}$. Total relic
density $\Omega h^{2}=0.114$, where $\Omega h^{2}_{D}=0.072$, $\Omega
h^{2}_{\phi_{2}}=0.0175$ and $\Omega h^{2}_{\phi_{3}}=0.0275$. We note that
the kink in $\xi_{i}$ around $T=10^{-4}$ GeV is due to the variation in the
degrees of freedom. Specifically, this is due to neutrino decoupling and
$e^{+}e^{-}$ annihilation. Similar kinks appear in Figs. 4, 5 and 6.

Further, the energy densities of the hidden sectors depend on the temperature
evolution of the respective sectors. This is illustrated in the right panel of
Fig. 2. Here we notice how $\xi_{1}$ rises almost immediately starting from
large $T$ while $\xi_{2}$ and $\xi_{3}$ remain at their lowest values before
picking up at $T\sim 30$ GeV. Thus the process of thermalization begins with
the first hidden sector which then permeates to the second and the third
hidden sectors. However, while $\xi_{1}$ reaches 1 the other hidden sectors
remain out of thermal equilibrium with $\xi_{2},\xi_{3}\ll 1$ and
$\xi_{2}\neq\xi_{3}$. For the benchmark considered in Fig. 2, the relic
density is shared among $D$, $\phi_{2}$ and $\phi_{3}$ ($\phi_{1}$ decays back
to the SM particles). Note that $\phi_{2}$ can still have a variety of off-
shell decays:
$\phi_{2}\to\phi_{1}^{*}\phi_{1}^{*}\to\bar{q}q\bar{q}q,\ell^{+}\ell^{-}\ell^{+}\ell^{-},\bar{q}q\ell^{+}\ell^{-},\bar{q}q\nu\nu,\ell^{+}\ell^{-}\nu\nu,\nu\nu\nu\nu$.
However, if the mass of $\phi_{2}$ is constrained so that $m_{\phi_{2}}\simeq
1$ MeV, all the decays except for the final state $\nu\nu\nu\nu$ will be
forbidden. In this case the decay is extremely suppressed since
$\Gamma_{\phi_{2}}\propto\kappa_{1}^{2}y^{4}_{f}$,

$\displaystyle\Gamma_{\phi_{2}}$
$\displaystyle=\frac{\kappa_{1}^{2}y^{4}_{f}}{16384\pi^{5}m_{\phi_{2}}}\int_{0}^{m_{\phi_{2}}^{2}}dq_{1}^{2}\int_{0}^{(m_{\phi_{2}}-q_{1})^{2}}dq_{2}^{2}$
$\displaystyle\hskip
113.81102pt\times\frac{\lambda^{1/2}(q_{1}^{2},q_{2}^{2},m_{\phi_{2}}^{2})q_{1}^{2}q_{2}^{2}}{[(q_{1}^{2}-m_{\phi_{1}}^{2})^{2}+m_{\phi_{1}}^{2}\Gamma_{\phi_{1}}^{2}][(q_{2}^{2}-m_{\phi_{1}}^{2})^{2}+m_{\phi_{1}}^{2}\Gamma_{\phi_{1}}^{2}]}.$
(4.21)

Thus we have $\tau_{\phi_{2}}\sim\mathcal{O}(10^{25})$ yrs for the input
mentioned in the caption of Fig. 2.

Figure 3: Plots of the energy exchange rate between different sectors versus
the visible sector temperature $T$. The input for the left panel (Case 1) are
same as in Fig. 2 while the input for the right panel (Case 2) are the same as
in Fig. 5.

As is apparent from the preceding discussion, the thermalization between the
first hidden sector and the visible sector is mainly driven by the processes
$\phi_{1}\leftrightarrow i\,\bar{i}$, where $i\,\bar{i}$ refers to the SM
particles. To illustrate this further, we plot in Fig. 3 the energy exchange
rate for different processes as well as the energy density dilution term due
to the Hubble expansion, $4H\rho_{1}$. In the left panel of Fig. 3, we
consider case 1 which corresponds to the analysis of Fig. 2. We see that all
processes are subdominant to the dilution term except for $i\,\bar{i}\to
D\bar{D}$ and $\phi_{1}\leftrightarrow i\,\bar{i}$. Of particular interest are
the processes $\phi_{1}\leftrightarrow i\,\bar{i}$ (green and brown curves)
which become comparable in rate near $10^{-2}$ GeV which is the region where
the two sectors thermalize as seen in the right panel of Fig. 2. The rates
slightly overtake $4H\rho_{1}(T_{1})$ for a short period of time before dying
out due to Boltzmann suppression of the number density. The right panel of
Fig. 3 pertains to case 2 whose discussion we leave for the next section.

Since the visible and the hidden sectors constitute a chain, a change in one
coupling will propagate to all sectors. In Fig. 4 we consider a smaller
$y_{f}$ value, i.e. a weaker coupling between the visible and first hidden
sector. As expected, the relic density of $D$ becomes smaller and due to a
smaller abundance of dark species in the first hidden sector, the yields of
$\phi_{2}$ and $\phi_{3}$ become negligible (the yield of $\phi_{3}$ does not
appear in the left panel as it is extremely small). Furthermore, the first
hidden sector no longer thermalizes with the visible sector as evident from
the right panel of Fig. 4. We also note the change in the evolution of
$\xi_{2}$ and $\xi_{3}$.

Figure 4: Same input as in Fig. 2 but for $y_{f}=1\times 10^{-10}$. Relic
density $\Omega h^{2}_{D}=3.6\times 10^{-2}$ and $\Omega
h^{2}_{\phi_{2}}=\Omega h^{2}_{\phi_{3}}\sim 0$.

### Case 2: $m_{\phi_{2}}>2m_{\phi_{1}}$ and $m_{\phi_{2}}>2m_{\phi_{3}}$

For this mass hierarchy, both decay processes $\phi_{2}\to\phi_{1}\phi_{1}$
and $\phi_{2}\to\phi_{3}\phi_{3}$ are kinematically allowed. We solve the
coupled Boltzmann equations for this case and plot in Fig. 5 the yields of the
dark species as well as $\xi_{i}$ against the visible sector temperature.

Figure 5: The case of three initially cold dark sectors. Left panel: evolution
of the yields of the dark fermion $D$ and and of the spin zero particles of
the dark sector versus the visible sector temperature. Right panel: evolution
of the temperature ratios $\xi_{i}$. The input parameters are: $m_{D}=1.0$
GeV, $m_{\phi_{1}}=50$ MeV, $m_{\phi_{2}}=105$ MeV, $m_{\phi_{3}}=45$ MeV,
$y_{f}=1\times 10^{-9}$, $y_{D}=0.05$, $\kappa_{1}=10^{-7}$ and
$\kappa_{2}=6\times 10^{-6}$. Total relic density is $\Omega h^{2}=0.1198$,
where $\Omega h^{2}_{D}=0.102$ and $\Omega h^{2}_{\phi_{3}}=0.0178$.

As one would expect, $\phi_{1}$ as well as $\phi_{2}$ completely decay leaving
$D$ and $\phi_{3}$ as the two possible DM components. For the benchmark
considered in Fig. 5, the energy exchange rate due to the processes
$\phi_{2}\to\phi_{3}\phi_{3}$ and $\phi_{2}\to\phi_{1}\phi_{1}$ are comparable
and reach values larger than the energy density dilution due to the Hubble
expansion (see the magenta and orange curves in the right panel of Fig. 3)
which results in a rapid thermalization between the three hidden sectors as
seen in the right panel of Fig. 5. Notice that all three sectors remain very
much in close thermal contact while also reaching thermal equilibrium with the
visible sector for a period of time.

### Case 3: $m_{\phi_{2}}>2m_{\phi_{1}}$ and $m_{\phi_{3}}=0$

With the same setup as before we assume now that particle $\phi_{3}$ is
massless. The evolution of the yields of the dark particles as well as the
temperature ratios $\xi_{i}$ are shown in Fig. 6. The DM relic density is
satisfied as seen from the left panel.

Figure 6: The case of three initially cold dark sectors. Left panel: evolution
of the yields of the dark fermion $D$ and the scalars of the dark sector
versus the visible sector temperature. Right panel: evolution of the
temperature ratios $\xi_{i}$. The input parameters are: $m_{D}=5.0$ GeV,
$m_{\phi_{1}}=10$ MeV, $m_{\phi_{2}}=105$ MeV, $m_{\phi_{3}}=0$,
$y_{f}=6.7\times 10^{-10}$, $y_{D}=0.17$ and $\kappa_{1}=\kappa_{2}=10^{-7}$.
Total relic density is $\Omega h^{2}=\Omega h^{2}_{D}=0.123$. As discussed in
subsection 5.2, this case leads to extra relativisitc degrees of freedom at
BBN time so that $\Delta N_{\rm eff}^{\rm BBN}=0.323$.

Even though $\phi_{3}$ resides in the third hidden sector, the coupling
$y_{f}$ between the first hidden sector and the visible has an important
effect on $\xi_{3}$ owing to the coupled nature of the system. Smaller values
of $y_{f}$ will lead to smaller $\xi_{3}$ and hence a lower contribution to
$\Delta N_{\rm eff}$. We study this in more detail in the next section.

##  5 Implications for observables

In this section we discuss the physical consequences of the formalism given
above. In subsection 5.1 we discuss the thermal effects of hidden sectors on
the proton-DM scattering cross sections. In subsection 5.2 we discuss the
contribution to dark radiation from hidden sectors.

###  5.1 Thermal effect of hidden sectors on direct detection of dark matter

In this subsection we discuss the effect of a multi-hidden sector model on the
determination of the relic density and on the spin-independent DM-proton cross
section. The first observation to be made is regarding the effect of such a
coupled system on the DM relic density. Here we consider case 1 but with
$\kappa_{1},\kappa_{2}\sim 10^{-9}$ so that the contributions from $\phi_{2}$
and $\phi_{3}$ to the relic density are negligible. We solve the coupled
Boltzmann equations for two cases: (A) $\xi_{i_{0}}=1$ and (B) $\xi_{i_{0}}\ll
1$. We note that cases A and B result in different DM relic density such that
$(\Omega h^{2})_{A}<(\Omega h^{2})_{B}$. The difference between the two can
reach up to 50% in some cases, especially for $y_{f}\lesssim 10^{-8}$. To see
why this is the case, one needs to examine the main annihilation process
responsible for DM depletion. For $m_{\phi_{1}}\ll m_{D}$, the DM relic
density is dominantly controlled by the annihilation process
$D\bar{D}\to\phi_{1}\phi_{1}$ whose thermal average depends on the temperature
$T_{1}$. In Fig. 7, we plot $n_{D}\langle\sigma v\rangle$ for
$D\bar{D}\to\phi_{1}\phi_{1}$ versus the visible sector temperature for cases
A and B. Also shown (in dashed line) is the Hubble parameter as well as the
process $D\bar{D}\to i\,\bar{i}$ for comparison. The left panel of Fig. 7,
which corresponds to a light mediator ($m_{\phi_{1}}=1$ MeV) and small
$y_{f}$, shows that $D\bar{D}\to\phi_{1}\phi_{1}$ is far more dominant than
$D\bar{D}\to i\bar{i}$ which can be explained by the relative sizes of the
couplings $y_{f}$ versus $y_{D}$. Furthermore, the rate of annihilation to
$\phi_{1}\phi_{1}$ for case B (blue curve) overtakes the expansion rate
earlier than for case A (red curve) and consequently decoupling also happens
earlier. Thus for case A, the annihilation process continues for a longer time
leading to a smaller relic density. This is in contrast to the case of a
heavier mediator ($m_{\phi_{1}}=1$ GeV) and a larger $y_{f}$ as shown in the
right panel of Fig. 7 where one can see that $D\bar{D}\to\phi_{1}\phi_{1}$
decouples at the same instant for both cases A and B and thus the final yield
is insensitive to whether the two sectors were initially at the same
temperature or not.

Figure 7: Display of $n_{D}\langle\sigma v\rangle$ for the processes
$D\bar{D}\to\phi_{1}\phi_{1}$ and $D\bar{D}\to i\bar{i}$ versus $T$ for two
cases: $\xi_{i_{0}}=1$ and $\xi_{i_{0}}\neq 1$. Input parameters for the left
panel are: $m_{D}=9.0$ GeV, $m_{\phi_{1}}=10$ MeV,
$m_{\phi_{2}}=m_{\phi_{3}}=1$ MeV, $y_{f}=4.5\times 10^{-10}$, $y_{D}=0.13$
and $\kappa_{1}=\kappa_{2}=5\times 10^{-8}$. Input parameters for the right
panel are: $m_{D}=5.0$ GeV, $m_{\phi_{1}}=1$ GeV,
$m_{\phi_{2}}=m_{\phi_{3}}=1$ MeV, $y_{f}=9.6\times 10^{-6}$, $y_{D}=0.18$ and
$\kappa_{1}=\kappa_{2}=1\times 10^{-5}$.

To illustrate more concretely the effect of hidden sectors on observables in
the visible sector, we use the latest constraints on the spin-independent (SI)
DM-proton scattering cross-section from the PandaX-II collaboration. The
collaboration presented their limits based on a secluded dark sector model of
ref. [28] which considers a thermal DM distribution in the early universe.
Using those limits, we recast the PandaX-II results to our model parameters
and display the limits in the left panel of Fig. 8 (solid lines). The analysis
which corresponds to the solid lines is for case A where we assume that the
sectors have the same initial temperature. The dashed lines correspond to case
B where the hidden sectors are initially much colder than the visible sector.
We note that for 1 MeV and 10 MeV mediators the limits have moved up, i.e.,
have become more relaxed while for a 1 GeV mediator there is virtually no
effect.

As discussed earlier, sectors at different temperatures lead to a larger relic
density and so to reduce the value of the latter, one increases the coupling
$y_{D}$ which leads to a larger SI cross-section. Thus the limits move up as
shown in the left panel of Fig. 8 (dashed lines). However, this seems not to
be the case for a 1 GeV mediator. Also as discussed earlier, the main
annihilation process decouples at the same time for cases A and B which
results in no change in the relic density. This is due to the fact that in
this case a larger $y_{f}$ is required for those limits ($y_{f}\gtrsim
10^{-7}$) which takes us closer to a pure freeze-out case as the one used by
PandaX-II to produce those limits. Thus as the coupling $y_{f}$ increases, it
gradually phases out the freeze-in mechanism [29] and we begin to agree with
the limits presented in ref. [27]. This can be clearly seen in the right panel
of Fig. 8 where the yield of $D$ (red curve) reaches close to its equilibrium
value (dashed brown curve) and then undergoes a freeze-out in a manner
different from what appears for the cases of smaller $y_{f}$ (see Figs. 2 and
4$-$6).

Figure 8: Left panel: A recast of the exclusion limits of the spin-independent
proton-DM cross section, $\sigma_{pD}$, given by the PandaX-II Collaboration
[27] for our model with two sets of initial hidden sector constraints, i.e.,
case A with $\xi_{i0}=1$ (solid lines) and case B with $\xi_{i0}<<1$ (dashed
lines) as a function of $m_{D}$ for three different values of the pseudo-
scalar $\phi_{1}$ mass, i.e., $m_{\phi_{1}}$=(1 MeV (red lines), 10 MeV (blue
lines) and 1 GeV (green lines)). The coupling $y_{D}$ for each point is chosen
to give the correct DM relic density. Right panel: a plot of the yield of the
dark particle species with the input: $m_{D}=48$ GeV, $m_{\phi_{1}}=1$ GeV,
$m_{\phi_{2}}=m_{\phi_{3}}=1$ MeV, $y_{f}=10^{-7}$, $y_{D}=0.55$ and
$\kappa_{1}=\kappa_{2}=1\times 10^{-5}$.

###  5.2 Dark radiation from hidden sectors

Hidden sectors in general contain extra relativistic degrees of freedom which
will contribute to the so-called dark radiation which is constrained by BBN
and CMB data. Thus, for example, for case 3, we have a massless spin zero
field $\phi_{3}$ which will remain in the spectrum and contribute to the
relativistic degrees of freedom in the universe as dark radiation. Such a
contribution over the existing neutrino species is determined by $\Delta
N_{\rm eff}$, so that

$\Delta N_{\rm eff}=\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\frac{\rho_{\rm
DS}}{\rho_{\gamma}},$ (5.1)

where the energy density of the dark sector is $\rho_{\rm
DS}=\rho_{D}(T_{1})+\rho_{\phi_{3}}(T_{3})$ and $\rho_{\gamma}$ is the photon
energy density. Based on the latest estimate of the effective number of
neutrino species, $N_{\rm eff}$, from BBN [30]

$N_{\rm eff}=2.88\pm 0.27,$ (5.2)

where one can deduce a $2\sigma$ upper bound on $\Delta N_{\rm eff}$, so that
$\Delta N_{\rm eff}<0.42$. If $\phi_{3}$ is produced efficiently and reaches
an equilibrium distribution, then $\rho_{\phi_{3}}\propto T_{3}^{4}$ and from
Eq. (5.1) we have $\Delta N_{\rm eff}\sim(T_{3}/T)^{4}=\xi_{3}^{4}$. So of
particular importance to the determination of $\Delta N_{\rm eff}$ is the
evolution of $\xi_{3}$ which is tightly coupled to the rest of the evolution
parameters. In Fig. 6 we display the evolution of the yields and $\xi_{i}$ for
the case of a massless $\phi_{3}$. For the set of assumed couplings, we have
$\Delta N_{\rm eff}^{\rm BBN}=0.323$ consistent with the BBN bound.

As noted, the constraint given in Eq. (5.2) comes from BBN while other
estimates of $N_{\rm eff}$ rely on the power spectrum data and CMB lensing
from Planck and Baryon Acoustic Oscillations (BAO) [11, 31]. Further, local
measurements of the Hubble parameter $H_{0}$ [32, 33, 34] show a significant
deviation from estimates of $H_{0}$ at earlier times (see ref. [35] for a
review of tensions in cosmology and ways to alleviate them and ref. [36] for a
cosmologically consistent particle physics model aimed at alleviating the
Hubble tension). A fit of the cosmological parameters to the CMB, BAO and
local $H_{0}$ data points to a larger $\Delta N_{\rm eff}$ at the surface of
last-scattering. The corresponding ranges for $\Delta N_{\rm eff}$ based on
these analyses are shown in Fig. 9 as red and blue bands. The future CMB-S4
experiment [37, 38] is expected to reach a much better sensitivity in
determining $N_{\rm eff}$ which will put stronger constraints on models
predicting extra light relics. The projected sensitivity of CMB-S4 is shown as
the dashed blue line in Fig. 9.

The curves in Fig. 9 show the variation in $\Delta N_{\rm eff}$ as a function
of the coupling $y_{f}$ for two sets of $\kappa_{1},\kappa_{2}$. Since the
hidden sectors form a chain with the visible sector then a change in $y_{f}$
will eventually propagate to the third hidden sector thus affecting the energy
density of $\phi_{3}$ which is a function of the temperature $T_{3}$ (here
enters the dependence on $\xi_{3}$). Larger $y_{f}$ leads to larger
contributions to $\Delta N_{\rm eff}$ until it saturates for $y_{f}\gtrsim
10^{-9}$ for the particular benchmark under consideration. Furthermore, larger
values of $\kappa_{1}$ and $\kappa_{2}$ lead to a stronger coupling between
the hidden sectors which eventually leads to an enhanced contribution to
$\Delta N_{\rm eff}$ as evident when comparing the two curves in Fig. 9. The
color gradient of the curves shows the change in the relic density as $y_{f}$
varies. One can see that in the range of interest $\Omega h^{2}\leq 0.12$.
Based on our model, CMB-S4 can probe feeble couplings between the visible and
hidden sectors down to $\mathcal{O}(10^{-10})$. Deviations from $N_{\rm
eff}^{\rm SM}$ which are $\mathcal{O}(10^{-2})$ can point to relativistic
relics that reside in hidden sectors at temperatures much lower than the
visible sector temperature.

Figure 9: A plot of $\Delta N_{\rm eff}$ for the case of massless $\phi_{3}$
as a function of $y_{f}$ for sets of values of $\kappa_{1}$ and $\kappa_{2}$.
The color gradient of the curves reflects the DM relic density. Other input
parameters are same as in Fig. 6.

##  6 Conclusion

In a variety of supergravity, string and brane models one encounters hidden
sectors which although neutral under the Standard Model gauge group can
interact feebly with the visible sector via a variety of portals. In this work
we have extended previous works where the visible sector was coupled to one or
two hidden sectors to the case where an arbitrary number of hidden sectors are
sequentially coupled to the visible sector. Thus the visible sector couples
only to hidden sector 1 and the latter couples to hidden sector 2, and so on,
i.e., in general the $i$-th hidden sector couples to hidden sectors ($i\pm
1$). In general the visible and the hidden sectors will be in different heat
baths, each with their own temperature. Of prime interest is to obtain the
evolution equations for the ratio of temperatures in different heat baths
after one has chosen a given temperature to serve as the clock. In the work
here we use the temperature $T$ of the visible sector as the clock and obtain
the evolution equations for the ratios $T_{i}/T,i=1,2,\cdots,n$. A detailed
analysis for the case of the visible sector coupled to three hidden sectors is
given with an application to a simple model with scalar portals. We solve a
set of coupled Boltzmann equations for particle yields and the temperature
ratios $\xi_{i}$ to trace the evolution of the dark species’ population as
well as the temperatures in each of the hidden sectors. We’ve shown that the
presence of a chain of hidden sectors affects the DM relic density which comes
out to be larger than one would get assuming all sectors are at the same
temperature. A direct consequence of this observation is a modification of the
recent PandaX-II limits on the SI DM-proton scattering cross sections for the
case of self-interacting DM. Also we have seen that the hidden sectors
thermalize with each other once the decay channels across the sectors are
active. Further, assuming that one of the spin zero particles is massless,
contributions to $\Delta N_{\rm eff}$ become an important constraint from
cosmology. Weaker couplings between the different sectors lead to a smaller
dark species’ energy density and temperature ratios $\xi_{i}$. We expect the
future CMB-S4 experiment to probe $\Delta N_{\rm eff}$ down to
$\mathcal{O}(10^{-2})$ which would constraint our model in the very feeble
coupling regime. We note that applications of the formalism developed here to
other portals and to larger values of $n$ beyond $n=3$ and possibly to explore
the region of large $n$ should be of interest for future investigations.
Finally, as noted in section 2, a chain of $n$ hidden sectors appears in a
broad class of beyond the Standard Model physics scenarios such as
moose/quiver theories. Thus an analysis of the type discussed above may find
application in this broad class of models.

## Acknowledgments

The research of AA was supported by the BMBF under contract 05P21PMCAA and by
the DFG through the Research Training Network 2149 “Strong and Weak
Interactions – from Hadrons to Dark matter”, while the research of PN was
supported in part by the United States NSF Grant PHY-1913328.

##  Appendix A Mass and temperature dependence of $g_{\rm eff}$ and $h_{\rm
eff}$

The total effective number of energy density and entropy degrees of freedom in
a sector containing a Dirac fermion and a boson are given by

$\displaystyle g_{\rm eff}$ $\displaystyle=g^{b}_{\rm
eff}+\frac{7}{8}g^{f}_{\rm eff},~{}~{}~{}\text{and}~{}~{}~{}h_{\rm
eff}=h^{b}_{\rm eff}+\frac{7}{8}h^{f}_{\rm eff},$ (A.1)

where the superscript $b~{}(f)$ indicates bosonic (fermionic) degrees of
freedom where at temperature $T$, $g^{b}_{\rm eff},h^{b}_{\rm eff}$ and
$g^{f}_{\rm eff},h^{f}_{\rm eff}$ are given by [39]

$\displaystyle g^{b}_{\rm eff}$
$\displaystyle=\frac{15d_{b}}{\pi^{4}}\int_{x_{b}}^{\infty}\frac{\sqrt{x^{2}-x_{b}^{2}}}{e^{x}-1}x^{2}dx,~{}~{}~{}\text{and}~{}~{}~{}h^{b}_{\rm
eff}=\frac{15d_{b}}{4\pi^{4}}\int_{x_{b}}^{\infty}\frac{\sqrt{x^{2}-x_{b}^{2}}}{e^{x}-1}(4x^{2}-x_{b}^{2})dx,$
(A.2) $\displaystyle g^{f}_{\rm eff}$
$\displaystyle=\frac{15d_{f}}{\pi^{4}}\int_{x_{f}}^{\infty}\frac{\sqrt{x^{2}-x_{f}^{2}}}{e^{x}+1}x^{2}dx,~{}~{}~{}\text{and}~{}~{}~{}h^{f}_{\rm
eff}=\frac{15d_{f}}{4\pi^{4}}\int_{x_{f}}^{\infty}\frac{\sqrt{x^{2}-x_{f}^{2}}}{e^{x}+1}(4x^{2}-x_{f}^{2})dx.$

where $d_{b}~{}(d_{f})$ are the intrinsic degrees of freedom for a bosonic
(fermionic) field, so that for a massive spin 1, $d_{b}=3$ and for a Dirac
fermion $d_{f}=4$. Further, $x_{b}$ and $x_{f}$ are defined so that
$x_{b}=m_{b}/T$ and $x_{f}=m_{f}/T$. The limit $x_{b}\to 0$ gives $g^{b}_{\rm
eff}=h^{b}_{\rm eff}\to d_{b}$ and the limit $x_{f}\to 0$ gives $g^{f}_{\rm
eff}=h^{f}_{\rm eff}\to d_{f}$.

##  Appendix B Collision and source terms

In this appendix we list the collision terms appearing in Eq. (4.20) and the
source terms appearing in Eqs. (4.17)$-$(4.19) for particles $D$, $\phi_{1}$,
$\phi_{2}$ and $\phi_{3}$. The respective collision terms are given by

$\displaystyle\mathcal{Q}_{D}=$ $\displaystyle\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to f\bar{f}}(T)(Y_{D}^{\rm
eq}(T)^{2}-Y_{D}^{2})-\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to\phi_{1}\phi_{1}}(T_{1})\left(Y^{2}_{D}-Y_{D}^{\rm
eq}(T_{1})^{2}\frac{Y^{2}_{\phi_{1}}}{Y^{\rm
eq}_{\phi_{1}}(T_{1})^{2}}\right),$ (B.1)
$\displaystyle\mathcal{Q}_{\phi_{1}}=$ $\displaystyle\langle\sigma
v\rangle_{\phi_{1}\phi_{1}\to f\bar{f}}(T)(Y_{\phi_{1}}^{\rm
eq}(T)^{2}-Y_{\phi_{1}}^{2})+\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to\phi_{1}\phi_{1}}(T_{1})\left(Y^{2}_{D}-Y_{D}^{\rm
eq}(T_{1})^{2}\frac{Y^{2}_{\phi_{1}}}{Y^{\rm
eq}_{\phi_{1}}(T_{1})^{2}}\right)$ $\displaystyle+\langle\sigma
v\rangle_{f\bar{f}\to\phi_{1}}(T)\left(Y_{q}^{\rm eq}(T)^{2}-Y_{q}^{\rm
eq}(T)^{2}\frac{Y_{\phi_{1}}}{Y^{\rm
eq}_{\phi_{1}}(T_{1})}\right)-\langle\sigma
v\rangle_{\phi_{1}\phi_{1}\to\phi_{2}\phi_{2}}(T_{1})\left(Y^{2}_{\phi_{1}}-Y_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{Y^{2}_{\phi_{2}}}{Y^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)$ $\displaystyle+\langle\sigma
v\rangle_{\phi_{3}\phi_{3}\to\phi_{1}\phi_{1}}(T_{3})\left(Y^{2}_{\phi_{3}}-Y_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{Y^{2}_{\phi_{1}}}{Y^{\rm
eq}_{\phi_{1}}(T_{1})^{2}}\right)-\langle\sigma
v\rangle_{\phi_{1}\phi_{1}\to\phi_{2}}(T_{1})\left(Y_{\phi_{1}}^{2}-Y_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{Y_{\phi_{2}}}{Y^{\rm eq}_{\phi_{2}}(T_{2})}\right),$ (B.2)
$\displaystyle\mathcal{Q}_{\phi_{2}}=$ $\displaystyle\langle\sigma
v\rangle_{\phi_{1}\phi_{1}\to\phi_{2}\phi_{2}}(T_{1})\left(Y^{2}_{\phi_{1}}-Y_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{Y^{2}_{\phi_{2}}}{Y^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)+\langle\sigma
v\rangle_{\phi_{1}\phi_{1}\to\phi_{2}}(T_{1})\left(Y_{\phi_{1}}^{2}-Y_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{Y_{\phi_{2}}}{Y^{\rm eq}_{\phi_{2}}(T_{2})}\right)$
$\displaystyle+\langle\sigma
v\rangle_{\phi_{3}\phi_{3}\to\phi_{2}\phi_{2}}(T_{3})\left(Y^{2}_{\phi_{3}}-Y_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{Y^{2}_{\phi_{2}}}{Y^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)+\langle\sigma
v\rangle_{\phi_{3}\phi_{3}\to\phi_{2}}(T_{3})\left(Y_{\phi_{3}}^{2}-Y_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{Y_{\phi_{2}}}{Y^{\rm eq}_{\phi_{2}}(T_{2})}\right),$ (B.3)
$\displaystyle\mathcal{Q}_{\phi_{3}}=$ $\displaystyle-\langle\sigma
v\rangle_{\phi_{3}\phi_{3}\to\phi_{2}\phi_{2}}(T_{3})\left(Y^{2}_{\phi_{3}}-Y_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{Y^{2}_{\phi_{2}}}{Y^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)-\langle\sigma
v\rangle_{\phi_{3}\phi_{3}\to\phi_{1}\phi_{1}}(T_{3})\left(Y^{2}_{\phi_{3}}-Y_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{Y^{2}_{\phi_{1}}}{Y^{\rm
eq}_{\phi_{1}}(T_{1})^{2}}\right)$ $\displaystyle-\langle\sigma
v\rangle_{\phi_{3}\phi_{3}\to\phi_{2}}(T_{3})\left(Y_{\phi_{3}}^{2}-Y_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{Y_{\phi_{2}}}{Y^{\rm
eq}_{\phi_{2}}(T_{2})}\right)\Bigg{]}.$ (B.4)

The source terms are given by

$\displaystyle j_{1}$ $\displaystyle=2n_{f}^{\rm
eq}(T)^{2}\left(1-\frac{n_{D}^{2}}{n^{\rm
eq}_{D}(T_{1})^{2}}\right)J(f\bar{f}\to D\bar{D})+2n_{q}^{\rm
eq}(T)^{2}\left(1-\frac{n_{\phi_{1}}^{2}}{n^{\rm
eq}_{\phi_{1}}(T_{1})^{2}}\right)J(f\bar{f}\to\phi_{1}\phi_{1})$
$\displaystyle+n_{q}^{\rm
eq}(T)^{2}J(f\bar{f}\to\phi_{1})-n_{\phi_{1}}J(\phi_{1}\to
f\bar{f})+2\left(n_{\phi_{2}}-n_{\phi_{2}}^{\rm
eq}\frac{n_{\phi_{1}}^{2}}{n^{\rm
eq}_{\phi_{1}}(T_{1})^{2}}\right)J(\phi_{2}\to\phi_{1}\phi_{1})$
$\displaystyle-2\left(n_{\phi_{1}}^{2}-n_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{n_{\phi_{2}}^{2}}{n^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)J(\phi_{1}\phi_{1}\to\phi_{2}\phi_{2})$
$\displaystyle-2\left(n_{\phi_{1}}^{2}-n_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{n_{\phi_{3}}^{2}}{n^{\rm
eq}_{\phi_{3}}(T_{3})^{2}}\right)J(\phi_{1}\phi_{1}\to\phi_{3}\phi_{3}),$
(B.5) $\displaystyle j_{2}$
$\displaystyle=-\left(n_{\phi_{2}}-n_{\phi_{2}}^{\rm
eq}\frac{n_{\phi_{1}}^{2}}{n^{\rm
eq}_{\phi_{1}}(T_{1})^{2}}\right)J(\phi_{2}\to\phi_{1}\phi_{1})+2\left(n_{\phi_{1}}^{2}-n_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{n_{\phi_{2}}^{2}}{n^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)J(\phi_{1}\phi_{1}\to\phi_{2}\phi_{2})$
$\displaystyle+2\left(n_{\phi_{3}}^{2}-n_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{n_{\phi_{2}}^{2}}{n^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)J(\phi_{3}\phi_{3}\to\phi_{2}\phi_{2})-\left(n_{\phi_{2}}-n_{\phi_{2}}^{\rm
eq}\frac{n_{\phi_{3}}^{2}}{n^{\rm
eq}_{\phi_{3}}(T_{3})^{2}}\right)J(\phi_{2}\to\phi_{3}\phi_{3}),$ (B.6)
$\displaystyle j_{3}$
$\displaystyle=-2\left(n_{\phi_{3}}^{2}-n_{\phi_{3}}^{\rm
eq}(T_{3})^{2}\frac{n_{\phi_{2}}^{2}}{n^{\rm
eq}_{\phi_{2}}(T_{2})^{2}}\right)J(\phi_{3}\phi_{3}\to\phi_{2}\phi_{2})+2\left(n_{\phi_{2}}-n_{\phi_{2}}^{\rm
eq}\frac{n_{\phi_{3}}^{2}}{n^{\rm
eq}_{\phi_{3}}(T_{3})^{2}}\right)J(\phi_{2}\to\phi_{3}\phi_{3})$
$\displaystyle+2\left(n_{\phi_{1}}^{2}-n_{\phi_{1}}^{\rm
eq}(T_{1})^{2}\frac{n_{\phi_{3}}^{2}}{n^{\rm
eq}_{\phi_{3}}(T_{3})^{2}}\right)J(\phi_{1}\phi_{1}\to\phi_{3}\phi_{3}).$
(B.7)

In the above we have

$\displaystyle n_{f}^{\rm eq}(T)^{2}J(f\bar{f}\to XX)(T)$
$\displaystyle=\frac{T}{32\pi^{4}}\int_{4m^{2}_{D}}^{\infty}ds~{}\sigma_{XX\to
f\bar{f}}s(s-s_{0})K_{2}(\sqrt{s}/T),$ (B.8) $\displaystyle n_{q}^{\rm
eq}(T)^{2}J(f\bar{f}\to\phi_{1})(T)$
$\displaystyle=\frac{T}{32\pi^{4}}\int_{s_{0}}^{\infty}ds~{}\sigma_{f\bar{f}\to\phi_{1}}s(s-s_{0})K_{2}(\sqrt{s}/T),$
(B.9) $\displaystyle n_{X}(T_{i})^{2}J(XX\to YY)(T_{i})$
$\displaystyle=\frac{n_{X}(T_{i})^{2}}{8m^{4}_{X}T_{i}K_{2}^{2}(m_{X}/T_{i})}\int_{s_{0}}^{\infty}ds~{}\sigma_{XX\to
YY}s(s-s_{0})K_{2}(\sqrt{s}/T_{i}),$ (B.10) $\displaystyle n_{X}J(X\to
YY)(T_{i})$ $\displaystyle=n_{X}m_{X}\Gamma_{X\to YY}\,.$ (B.11)

## References

  * [1] L. Ackerman, M. R. Buckley, S. M. Carroll and M. Kamionkowski, Phys. Rev. D 79, 023519 (2009) doi:10.1103/PhysRevD.79.023519 [arXiv:0810.5126 [hep-ph]].
  * [2] L. Heurtier, Y. Mambrini and M. Pierre, Phys. Rev. D 99, no.9, 095014 (2019) doi:10.1103/PhysRevD.99.095014 [arXiv:1902.04584 [hep-ph]].
  * [3] A. Aboubrahim, W. Z. Feng, P. Nath and Z. Y. Wang, Phys. Rev. D 103, no.7, 075014 (2021) doi:10.1103/PhysRevD.103.075014 [arXiv:2008.00529 [hep-ph]].
  * [4] A. Aboubrahim, W. Z. Feng, P. Nath and Z. Y. Wang, JHEP 06, 086 (2021) doi:10.1007/JHEP06(2021)086 [arXiv:2103.15769 [hep-ph]].
  * [5] J. D. Bowman, A. E. E. Rogers, R. A. Monsalve, T. J. Mozdzen and N. Mahesh, Nature 555, no.7694, 67-70 (2018) doi:10.1038/nature25792 [arXiv:1810.05912 [astro-ph.CO]].
  * [6] A. Aboubrahim, P. Nath and Z. Y. Wang, JHEP 12, 148 (2021) doi:10.1007/JHEP12(2021)148 [arXiv:2108.05819 [hep-ph]].
  * [7] J. B. Muñoz and A. Loeb, Nature 557, no.7707, 684 (2018) doi:10.1038/s41586-018-0151-x [arXiv:1802.10094 [astro-ph.CO]].
  * [8] J. B. Muñoz, C. Dvorkin and A. Loeb, Phys. Rev. Lett. 121, no.12, 121301 (2018) doi:10.1103/PhysRevLett.121.121301 [arXiv:1804.01092 [astro-ph.CO]].
  * [9] H. Liu, N. J. Outmezguine, D. Redigolo and T. Volansky, Phys. Rev. D 100, no.12, 123011 (2019) doi:10.1103/PhysRevD.100.123011 [arXiv:1908.06986 [hep-ph]].
  * [10] A. Berlin, D. Hooper, G. Krnjaic and S. D. McDermott, Phys. Rev. Lett. 121, no.1, 011102 (2018) doi:10.1103/PhysRevLett.121.011102 [arXiv:1803.02804 [hep-ph]].
  * [11] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].
  * [12] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].
  * [13] R. Foot and S. Vagnozzi, Phys. Rev. D 91, 023512 (2015) doi:10.1103/PhysRevD.91.023512 [arXiv:1409.7174 [hep-ph]].
  * [14] R. Foot and S. Vagnozzi, JCAP 07, 013 (2016) doi:10.1088/1475-7516/2016/07/013 [arXiv:1602.02467 [astro-ph.CO]].
  * [15] I. Rothstein and W. Skiba, Phys. Rev. D 65, 065002 (2002) doi:10.1103/PhysRevD.65.065002 [arXiv:hep-th/0109175 [hep-th]].
  * [16] M. R. Douglas and G. W. Moore, [arXiv:hep-th/9603167 [hep-th]].
  * [17] N. Arkani-Hamed, A. G. Cohen and H. Georgi, Phys. Rev. Lett. 86, 4757-4761 (2001) doi:10.1103/PhysRevLett.86.4757 [arXiv:hep-th/0104005 [hep-th]].
  * [18] C. T. Hill, S. Pokorski and J. Wang, Phys. Rev. D 64, 105005 (2001) doi:10.1103/PhysRevD.64.105005 [arXiv:hep-th/0104035 [hep-th]].
  * [19] A. Aboubrahim, W. Z. Feng, P. Nath and Z. Y. Wang, [arXiv:2106.06494 [hep-ph]].
  * [20] B. Holdom, Phys. Lett. B 166, 196-198 (1986) doi:10.1016/0370-2693(86)91377-8
  * [21] B. Kors and P. Nath, Phys. Lett. B 586, 366-372 (2004) doi:10.1016/j.physletb.2004.02.051 [arXiv:hep-ph/0402047 [hep-ph]].
  * [22] D. Feldman, Z. Liu and P. Nath, Phys. Rev. D 75, 115001 (2007) doi:10.1103/PhysRevD.75.115001 [arXiv:hep-ph/0702123 [hep-ph]].
  * [23] K. Cheung and T. C. Yuan, JHEP 03, 120 (2007) doi:10.1088/1126-6708/2007/03/120 [arXiv:hep-ph/0701107 [hep-ph]].
  * [24] A. Aboubrahim, W. Z. Feng and P. Nath, JHEP 02, 118 (2020) doi:10.1007/JHEP02(2020)118 [arXiv:1910.14092 [hep-ph]].
  * [25] M. Du, Z. Liu and P. Nath, [arXiv:2204.09024 [hep-ph]].
  * [26] A. Aboubrahim, T. Ibrahim, M. Klasen and P. Nath, Eur. Phys. J. C 81, no.8, 680 (2021) doi:10.1140/epjc/s10052-021-09483-0 [arXiv:2012.10795 [hep-ph]].
  * [27] J. Yang et al. [PandaX-II], Sci. China Phys. Mech. Astron. 64, no.11, 11 (2021) doi:10.1007/s11433-021-1740-2 [arXiv:2104.14724 [hep-ex]].
  * [28] R. Huo, M. Kaplinghat, Z. Pan and H. B. Yu, Phys. Lett. B 783, 76-81 (2018) doi:10.1016/j.physletb.2018.06.024 [arXiv:1709.09717 [hep-ph]].
  * [29] L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 03, 080 (2010) doi:10.1007/JHEP03(2010)080 [arXiv:0911.1120 [hep-ph]].
  * [30] C. Pitrou, A. Coc, J. P. Uzan and E. Vangioni, Phys. Rept. 754, 1-66 (2018) doi:10.1016/j.physrep.2018.04.005 [arXiv:1801.08023 [astro-ph.CO]].
  * [31] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A1 (2020) doi:10.1051/0004-6361/201833880 [arXiv:1807.06205 [astro-ph.CO]].
  * [32] A. G. Riess, W. Yuan, L. M. Macri, D. Scolnic, D. Brout, S. Casertano, D. O. Jones, Y. Murakami, L. Breuval and T. G. Brink, et al. [arXiv:2112.04510 [astro-ph.CO]].
  * [33] A. G. Riess, S. Casertano, W. Yuan, L. Macri, J. Anderson, J. W. MacKenty, J. Bradley Bowers, K. I. Clubb, A. V. Filippenko and D. O. Jones, et al. Astrophys. J. 855, no.2, 136 (2018) doi:10.3847/1538-4357/aaadb7 [arXiv:1801.01120 [astro-ph.SR]].
  * [34] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri and D. Scolnic, Astrophys. J. 876, no.1, 85 (2019) doi:10.3847/1538-4357/ab1422 [arXiv:1903.07603 [astro-ph.CO]].
  * [35] E. Abdalla, G. Franco Abellán, A. Aboubrahim, A. Agnello, O. Akarsu, Y. Akrami, G. Alestas, D. Aloni, L. Amendola and L. A. Anchordoqui, et al. JHEAp 34, 49-211 (2022) doi:10.1016/j.jheap.2022.04.002 [arXiv:2203.06142 [astro-ph.CO]].
  * [36] A. Aboubrahim, M. Klasen and P. Nath, JCAP 04, no.04, 042 (2022) doi:10.1088/1475-7516/2022/04/042 [arXiv:2202.04453 [astro-ph.CO]].
  * [37] K. N. Abazajian et al. [CMB-S4], [arXiv:1610.02743 [astro-ph.CO]].
  * [38] K. Abazajian et al. [CMB-S4], [arXiv:2203.08024 [astro-ph.CO]].
  * [39] M. Hindmarsh and O. Philipsen, Phys. Rev. D 71, 087302 (2005) doi:10.1103/PhysRevD.71.087302 [arXiv:hep-ph/0501232 [hep-ph]].