A multi-temperature universe can allow a sub-MeV dark photon dark matter

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


Introduction
Supergravity and strings models typically contain hidden sectors with gauge groups including U (1) gauge group factors. These hidden sectors with U (1) gauge groups can interact feebly with the visible sector and interact feebly or with normal strength with each other. The fields in the visible and hidden sectors in general will reside in different heat baths and the universe in this case will be a multi-temperature universe. The multi-temperature nature of the universe becomes a relevant issue if the observables in the visible sector are functions of the visible and the hidden sector heat baths. Such is the situation if dark matter (DM) resides in the hidden sector but interacts feebly with the visible sector. In this case an accurate computation of the relic density requires thermal averaging of cross sections and decay widths which depend on temperatures of both the visible and the hidden sector heat baths. In this work we develop a theoretical formalism which can correlate the evolution of temperatures of the hidden sector and of the visible sectors (for the specific case of two hidden sectors) in an accurate way.
The formalism noted above is used in the investigation of a dark photon and dark fermions of hidden sectors as possible DM candidates. There exists a considerable literature in the study of dark photons [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] (for review see [19][20][21]) to which the interested reader is directed. While axions and dark photons in the light to ultralight mass region (from keV to 10 −22 eV) have been investigated [11,13,[22][23][24][25][26], the sub-MeV dark photon mass range appears difficult to realize. The problem arises in part because with the visible sector interacting with a hidden sector via kinetic mixing, the twin constraints that the dark photon has a lifetime larger than the age of the universe, and also produce a sufficient amount of DM to populate the universe are difficult to satisfy. In addition to the relic density constraint there is also the constraint on dark photon lifetime which needs to be larger than the age of the universe as well as a constraint from BBN on the light degrees of freedom. For the case of one hidden sector, these constraints are difficult to satisfy. Specifically, ∆N eff at BBN time depends on the ratio (T hid /T γ ) 4 . The BBN temperature is typically ∼ 0.1 MeV and, as will be seen later, for the case of one hidden sector the ratio (T hid /T γ ) 4 ∼ 1 which gives a contribution to ∆N eff at BBN time in excess of the current experimental constraint which from the combined data from BBN, BAO and CMB [27] is ∆N eff < 0.214. On the other hand, for the case of two hidden sectors it is possible to satisfy all the current experimental constraints and for that reason we will focus on the two hidden sector model which is the minimal extension of one hidden sector model as discussed below.
In this study we show that a sub-MeV dark photon as DM can indeed be realized in a simple extension of the Standard Model (SM) where the hidden sector is constituted of two sectors X 1 and X 2 where the sector X 1 has kinetic mixing with the visible sector while the sector X 2 has kinetic and mass mixings only with sector X 1 , c.f., Fig. 1. We assume that the hidden sector X 1 has a dark fermion D and its gauge boson Z decays before the BBN while the hidden sector X 2 has only a dark photon γ which has a lifetime larger than the age of the Universe. This set up is theoretically more complex because here one has three heat baths and the computation of the relic density thus depends on three temperatures, i.e., the temperature T of the visible sector, the temperature T 1 of the hidden sector X 1 and the temperature T 2 of the hidden sector X 2 . In the following we develop a formalism that allows one to compute temperatures of all three heat baths in terms of one common temperature, which can be chosen to be T , T 1 or T 2 . In the analysis below it is found convenient to choose the reference temperature to be T 1 . The outline of the rest of the paper is as follows: in section 2 we discuss the particle physics model used for our multi-temperature universe and in sections 3 and 4 we write the coupled Boltzmann equations with three temperatures and derive the temperature evolution of T and T 2 relative to T 1 . Thermalization between the different sectors and dark freeze-out are explained in section 5 with a numerical analysis and a discussion of the astrophysical constraints. Conclusions are given in section 6. Further analytical results are given in Appendices A through D. Dark sector X 1 Dark sector X 2 Figure 1: The exhibition of the model we consider in the paper. The standard model has a direct coupling with hidden sector X 1 with strength proportional to δ 1 , whereas hidden sector X 2 only interacts directly with X 1 with strength proportional to δ 2 , and thus X 2 interacts with the standard model only indirectly.

A model for a multi-temperature universe
As mentioned in the introduction, supergravity and string models contain hidden sectors. These hidden sectors in general would have both abelian and non-abelian gauge groups and some of them could interact feebly with the visible sector while others may interact with each other as shown in Fig. 2. Thus, for example, in D-brane models one gets U (N ) gauge groups where U (N ) → SU (N ) × U (1). These extra U (1) factors in general can acquire kinetic mixing with the U (1) Y of the visible sector. Further, the gauge bosons of the extra U (1)'s can acquire mass via the Stueckelberg mechanism and also have Stueckelberg mass mixings with the hypercharge gauge boson. Additionally if there are several hidden sector U (1)'s they can have also gauge and Stueckelberg mass mixings among themselves. Interestingly, the possible existence of these hidden sectors can have significant effect on model building in the visible sector. As an example one phenomenon which is deeply affected by the existence of hidden sectors is dark matter which we discuss in further detail below.
In general the dark sectors with a gauge symmetry will contain gauge fields as well as matter, but, as noted, typically they will have feeble interactions with the SM particles and likely also with the inflaton. This means that these particles would not be thermally produced in the reheat period after inflation but would then acquire their relic density via annihilation and decay of the SM particles. Thus in general the temperatures of the visible and the hidden sectors will be different from the visible sector as well as from each other. This means that their relic densities will be governed by a set of coupled Boltzmann equations which depend on different temperatures, i.e., temperature of the visible sector and those for the hidden sectors. One of the central items in understanding of how to deal with such coupled systems with sectors involving different temperatures is to understand fully how the temperatures of the hidden sectors grow relative to the visible sector temperature. The formalism of how to correlate the hidden and the visible sector temperatures was worked out for the case of the visible sector interacting with one hidden sector in [5]. However, a general framework does not exist. Here we discuss the case where there are two hidden sectors X 1 and X 2 where the hidden sector X 1 interacts with the visible sector, while the hidden sector X 2 interacts only with the hidden sector X 1 as shown in Fig. 1. In this case two functions η −1 = ξ = T 1 /T and V (T ) H1(T1) H3(T3) H5(T5) · · · · · · Hn(Tn) Hidden sector H2(T2) H4(T4) H6(T6) Ho(To) Figure 2: A schematic diagram exhibiting the coupling of the visible sector with multiple dark sectors and of the dark sectors among themselves. The visible sector may have direct couplings with some of the dark sectors, or indirect couplings with others via interactions among the entire hidden sector. ζ = T 2 /T 1 enter in the coupled Boltzmann equations and we derive differential equations for their evolution. The above setup has a direct application in achieving a sub-MeV dark photon as dark matter as we show later in this work. However, a consistent analysis of the coupled dynamics of the visible sector and two hidden sectors is significantly more complex. This work develops the necessary machinery to do so and which can be extended to multiple hidden sectors.
We assume that the two sectors X 1 and X 2 have U (1) X 1 and U (1) X 2 gauge symmetries and that the field content of X 1 is (C µ , D) where C µ is the gauge field and D is a dark fermion, and the field content of X 2 is the gauge field D µ and there is no dark fermion in the sector X 2 . We invoke a kinetic mixing [28,29] between the hypercharge field B µ of SM and C µ and a kinetic mixing between C µ and D µ as well as a Stueckelberg mass growth [30][31][32][33][34] for all the gauge fields as well as a Stueckelberg mass mixing between the fields C µ and D µ . The extended part of the Lagrangian including both the kinetic and mass mixings is where L SM contains the SM terms and the kinetic part is given by and in the unitary gauge the mass Lagrangian is given by The D fermion is assumed charged under U (1) X 1 with interaction g XD γ µ DC µ . Canonical normalization of Eqs. (2.2) and (2.3) is carried out in Appendix A which gives the mass eigenstates γ , Z , Z, γ. The neutral current Lagrangian contained in L SM for the mass eigenstates γ , Z , Z, γ describing the couplings between the vector bosons γ , Z , Z, γ with the visible sector fermions is given by Here θ w is the weak angle and e is defined as (2.5) As seen from Eqs. (2.1)−(2.3), the framework of the model allows for the inclusion of both kinetic and mass mixing between the hidden and visible sectors. However, in the analysis presented in this work we will only take kinetic mixing between the visible and the hidden sectors, so that M 2 = 0. This is done because in this case the 4 × 4 neutral vector boson mass square matrix factors into a block diagonal form consisting of two 2 × 2 matrices as shown in Eq. (A.4). In this case we can carry out a set of eight GL(2, R) transformations to put the kinetic energy of the visible and the hidden sectors in a canonical form and at the same time to lowest order the mass matrix is also in a canonical form. This allows us to use perturbation theory around stable minima of the standard model mass matrix and therefore to deduce the couplings between the hidden sectors and the SM particles as given in Tables 2 and 3. As seen in Appendix A, the analysis is rather non-trivial and significantly more involved than for the case of one hidden sector.
With the inclusion of M 2 , the analysis becomes analytically intractable and an exhibition of results corresponding to those of Tables 2 and 3 is difficult. However, we note that even with M 2 = 0, we still have both kinetic and mass mixing in the hidden sector. Thus δ 2 takes account of kinetic mixing and M 3 takes account of mass mixing between the hidden sectors 1 and 2. So in summary in this model the visible sector has only kinetic mixing with the hidden sector 1 while the hidden sector 1 has both kinetic and mass mixing with hidden sector 2. It is seen that the mass mixing from M 3 does have significant effect on the model predictions, e.g., on the mass of the dark photon and on the relic density as seen in model point (d) in Table 1.
With M 2 = 0, the neutral current Lagrangian for coupling to the hidden sector fermion is given by where 3 Boltzmann equations for yields with three bath temperatures The relic densities of the dark photon and of the dark fermion arise in part from a freeze-in mechanism [35][36][37][38][39]. In general the visible sector and the dark sectors will have different temperatures [40][41][42][43][44][45] (a similar setup has been considered in Ref. [45] but with a different particle content, couplings and no explicit multi-temperature evolution. Also the DM candidate was not the dark photon as in our case). As mentioned above, we consider three different temperatures corresponding to the temperatures of the visible sector T and of the two hidden sectors, T 1 for X 1 and T 2 for X 2 . Defining the yield Y = n/s, where n is the number density and s is the entropy density, and the bath functions η and ζ so that T = ηT 1 and T 2 = ζT 1 , we write the Boltzmann equations for the yields as where H is the Hubble parameter given by with ρ = ρ v + ρ 1 + ρ 2 and s being the energy and entropy densities given by and the quantities J D , J Z and J γ are defined in terms of the collision terms as The collision terms C D , C Z and C γ are given by Eqs. (C.7)−(C.9) which allow us to write J D , J Z and J γ in terms of the yield as The thermally averaged cross sections are calculated using Eq. (C.11) for a single initial state temperature and with Eq. (C.12) in the case of two different initial state temperatures, while the thermal averaging of the width is given by Eq. (C.10). Note that the SM effective energy and entropy degrees of freedom (g v eff and h v eff ) are read from tabulated values [46,47], while those pertaining to the hidden sectors (g 1eff , g 2eff , h 1eff , h 2eff ) are calculated using Eqs. (C. 16), (C.17) and (C.18). Since the above Boltzmann equations are dependent on the parameters η and ζ, then one must consider the evolution of those parameters with temperature. We discuss the formalism in the next section.

Temperature evolution in the dark sectors versus in the visible sector
In this section we derive the evolution equations for temperatures T 1 and T 2 in the dark sectors and T in the visible sector, i.e., T /T 1 and T 2 /T 1 as a function of T . However, for numerical integration purposes it is found more convenient to use T 1 as the reference temperature. Thus we are interested in deriving the evolution equations for dη/dT 1 and dζ/dT 1 where we recall that η and ζ are defined so that To this end we look at the equations for the energy densities in the visible and the hidden sectors. In the analysis, we encounter the quantity dρ/dT 1 . The main difficulty in computing the quantity dρ/dT 1 is that ρ is constituted of three parts, ρ = ρ v + ρ 1 + ρ 2 which depends on three temperatures i.e., ρ v is controlled by T , ρ 1 is controlled by T 1 and ρ 2 is controlled by T 2 . Thus we need to express dρ v /dT 1 in terms of dρ v /dT and dρ 2 /dT 1 in terms of dρ 2 /dT 2 . Using the definitions of η and ζ we can write This means that a determination of dρ v /dT 1 and dρ 2 /dT 1 requires dη/dT 1 and dζ/dT 2 . Next, we derive the evolution equations for these quantities.
We note that ρ v , ρ 1 , ρ 2 satisfy the following evolution equations where j v , j 1 , j 2 are the corresponding sources. Instead of time we will use temperature so we will need to convert derivatives with respect to time to derivatives with respect to temperature. We note now that for any given temperatures T i the time derivative of temperature is given by As discussed above, we choose T 1 to be the reference temperature and the evolution equation for ρ v in this case can be written as From Eq. (4.5) we can deduce that In a similar fashion starting with the equation for dρ 2 /dt we can deduce Eqs. (4.6) and (4.7) are two coupled equations involving dρ v /dT 1 and dρ 2 /dT 1 which give the solution Using Eqs. (4.2) and (4.8) one can then obtain the relations Using the fact that j v +j 1 +j 2 = 0, eliminating j v in favor of j 1 and j 2 and inserting Eq. (4.9) in Eq. (4.10), one can further simplify Eq. (4.10) to cast dη/dT 1 and dζ/dT 1 in their final form as (4.14) The source terms j 1 and j 2 are given by One should not confuse the variable s in Eqs. (4.15) and (4.16) with the one in Eq. (4.20). The former is the entropy density while the latter corresponds to the Mandelstam variable.

Thermalization and dark freeze-out
In the analysis we make certain that the relic density of the dark relics is consistent with the Planck data [27], along with a 2σ corridor from theoretical calculations. Contribution to the relic density arise from γ and D while Z decays before BBN and is removed from the spectrum. In Table 1 we present four benchmarks which satisfy all the experimental constraints.
The relic density shown is that of γ while that of D is only O(10 −6 ) or less and thus negligible. We note here that in addition to the particle physics interactions generating the dark photon relic density, one could have in addition gravitational production [16][17][18]. However, such a production is highly dependent on the reheat temperature which is model dependent. From Eq. (46) of Ref. [16], for the case when the dark photon mass is ∼ 1 MeV, one finds that the gravitational production of dark photon would be suppressed when the reheat temperature H I < 10 11 GeV. So one may think of our model being valid for this restricted class of inflationary models.
Calculation of the relic density requires determining the yields by numerically solving the five stiff coupled equations, Eqs.  and γ as a function of the hidden sector temperature T 1 are shown in Fig. 3 for benchmarks (a) and (b). As the universe cools, the number densities of D, Z and γ increase gradually. At around T 1 = 100 GeV, the dark fermions D start to freeze-out, and the blue curve becomes flat around T 1 = 0.1 GeV. After the dark fermion decouples, the dynamics of Z and γ is affected mainly by SM particles freeze-in processes after T 1 = 0.1 GeV. However, since Z is unstable its density depletes to zero at T 1 ∼ 10 −4 GeV. The only particles that contribute to the relic density then are D (blue curve) and γ (yellow curve). As noted the analysis gives dark photon as the dominant component of DM. We note that the number changing processes in hidden sector 1 are driven by DD ←→ Z Z owing to the sizable value of the coupling g X . Since m Z m D , the reaction Z Z → DD shuts off early on as the temperature drops while the reverse reaction remains active. This causes a significant drop in Y D and as a consequence Y Z rises sharply as shown in Fig. 3. This is followed by a dramatic drop in Y Z due to the decay of Z to SM fermions.  Table 1.
The upper panel of Fig. 4 gives the evolution of ξ = T 1 /T as a function of T 1 which shows that ξ rises until it thermalizes with the visible sector, i.e., ξ ∼ 1. The middle panel of this figure gives the evolution ζ = T 2 /T 1 as a function of T 1 while the bottom panel of Fig. 4 gives the evolution of κ = T 2 /T as function of T . We note that X 2 does not thermalize with X 1 . This happens because the energy injection from X 1 into X 2 is not efficient enough. Consequently T 2 /T 1 which also has implications for ∆N eff as we explain later.
We note in passing that even though hidden sector 1 thermalizes with the visible sector, there is a distinction between how that happens for the case of the hidden sector versus the visible sector. In the presence of a coupling induced either by kinetic mixing or by mass mixing, the visible sector and the hidden sectors will eventually thermalize as long as the sectors do not thermally decouple according to the second law of thermodynamics. This is what happens in the top left panel of Fig. 4. However, we note that the speed with which a dark photon thermalizes is much slower relative to a visible sector particle such as a quark which has almost instantaneous thermalization with the photon background. Further, the thermalization will cease once the particles in the hidden sector fully decouple from the visible sector or from each other as seen in the right top panel of Fig. 4. This is the case for hidden sector 2.
In the left panel of Fig. 5 we show n σv and the thermally averaged Z decay width as a function of T 1 for benchmark (a). Also shown is the Hubble parameter H(T 1 ). As evident, while Z can enter into equilibrium with the visible sector for a period of time, the dark photon barely does so. We indicate by arrows the point at which the dark freeze-out of D and γ occurs. The dark photons decouple earlier followed by the dark fermions which is also evident in Fig. 3. We note that Γ Z overtakes H(T 1 ) at lower temperatures contributing to the depletion of Z number density. It is of interest to ask how thermal equilibrium of dark photons occurs once they are produced. Such an equilibrium can be achieved in our model by considering a massless complex scalar field φ in the second hidden sector with interactions with the dark photon of the type (κφ † ∂ µ φA γ µ + h.c.). Elastic scattering between γ and φ, γ φ → γ φ, can keep the dark photon in local thermal equilibrium. We follow the method of Ref. [49] to determine the temperature of kinetic decoupling of dark photons. For κ ∼ 10 −4 , we find that kinetic decoupling occurs at 10 keV which is much later than chemical decoupling which happens around 0.1 GeV. The right panel of Fig. 5 shows the temperature of kinetic decoupling. It is to be noted that the small value of κ has a minimal effect on the relic density of γ .
We note that the parameter space of Z and γ is constrained by experiments such as BaBar, CHARM and other beam-dump experiments as well as by astrophysical data from Supernova SN1987A and stellar cooling. Those limits become even stronger when the dark photon is assumed to be the dark matter particle. Thus, measurements of heating rates of the Galactic center cold gas clouds [50], the temperature of the diffuse X-ray background [51] as well as that of the intergalactic medium at the time of He ++ reionization [52][53][54] are affected by early γ → 3γ decays. Further constraints can be derived from energy injection during the dark ages [55] and spectral distortion of the CMB [54]. The presence of a long-lived sub-MeV particle species can contribute to the relativistic number of degrees of freedom ∆N eff during BBN and recombination [56]. All those constraints can exclude a sub-MeV dark photon down to a kinetic mixing coefficient O(10 −13 ).
In the model discussed here, the dark photon resides in a hidden sector X 2 that does not interact directly with the visible sector. Instead, the direct interaction is between the two hidden sectors X 1 and X 2 via kinetic and mass mixings. Since X 1 mixes kinetically with the visible sector, the interaction between X 2 and the visible sector becomes doubly suppressed and all coupling will be proportional to δ 1 (δ 2 − sin β). The quantity sin β is due to the mass mixing between the hidden sector and such a term can impart millicharges to the D fermions. However, the coupling between the photon and the D fermions is not only suppressed by δ 1 δ 2 sin β but also by the mass ratio as evident from the expression of c γ in Eq. (2.8). Therefore even for a modest value of sin β ∼ 10 −2 , the millicharges are very small and do not constitute a significant constraint on the model. This doubly suppressed coupling between the dark photon and the SM can alleviate the present constraints mainly from γ → 3γ as seen in Fig. 6, which still removes a part of the parameter space of our model. This constraint is derived from measurements of the intergalactic diffuse photon background. Another decay channel for the dark photon is to two neutrinos. This leads to a possible neutrino flux but experiments have not reached the required sensitivity to probe masses in the sub-MeV region. Experiments such as IceCube [57] have constrained only very heavy dark matter decays. The region of the parameter space which would produce a dark photon relic density within 2σ of the experimental value is shown in both panels of Fig. 6 and labeled 'Freeze-in'. For the case when a dominant component of the relic density arises from gravitational production, the parameter space of our model will be enlarged. The enlarged regions which are represented by the hatched area in Fig. 6. This area accommodates for a dark photon relic density down to ∼ 10 −4 .
It is argued in Ref. [51] that a dark photon with direct kinetic mixing with the SM can only give a subdominant contribution to the relic density and that such an observation can be dismissed if another production mechanism is in effect. The model discussed here presents exactly this counter argument required to produce a dominant dark photon dark matter. The main production mechanism for the dark photon in the current analysis is not via the freeze-in mechanism from the visible sector, iī → γ and iī → γ γ , because of the doubly suppressed coupling (see left panel of Fig. 7) but rather from interactions between the hidden sector particles. Thus, processes such as DD → γ γ have cross-sections proportional to g X δ 2 . As shown in Fig. 6, the sizes of g X and δ 2 are in the required ranges to produce a dark photon relic density which dominates that of D. A smaller value of g X will reduce the dark photon yield as expected (see right panel of Fig. 7). We note in passing that the gauge coupling g X is constrained by the main annihilation channel DD → Z Z → 4e from the Planck experiment [58,59] and for m D < 1 GeV, g X > 0.1 is excluded. However, this constraint does not exist for our model since the relic abundance of our D fermions is negligible. The effect of the forward process DD → γ γ can be clearly seen in Fig. 3 where the drop in the D fermions yield at a certain temperature is followed by a rise in the yield of γ . It is worth mentioning that the reverse process γ γ → DD works on reducing the number density of γ on the expense of D, but this process shuts off early on as shown in Fig. 7 allowing DD → γ γ to completely take over for lower temperatures. Figure 6: Exclusion limits from terrestrial and astrophysical experiments on dark photon which has kinetic mixings with the SM sector. Excluded regions are due to constraints from experiments which include electron and muon g − 2 [60], BaBar [61], CHARM [62,63], NA48 [64], E137 [65,66], NA64 [67,68], E141 [69] and ν-CAL [63,70,71]. The limits are obtained from darkcast [72]. The strongest constraints on a dark photon (mass less than 1 MeV) come from Supernova SN1987A (including a robustly excluded region and systematic uncertainties) [73], stellar cooling [74] and from decay to 3γ on cosmological timescales [51,75]. The islands in pink are constraints from BBN. The region where the freeze-in relic density is satisfied within 2σ of the experimental constraint is shown in different shades of blue corresponding to different choices of m D . The hatched area represents an enlargement to the 2σ region allowing a relic density as low as ∼ 10 where T γ = T . Using the ratio of the temperatures T 2 /T < 0.1 from Fig. 4, one finds that ∆N eff from dark photons is O(10 −4 ) which makes a negligible contribution to the SM N eff . With the inclusion of a complex scalar, the model has now five new bosonic degrees of freedom but this still gives a small contribution and does not violate the bound on ∆N eff . We note that the suppressed value of T 2 /T which arises from the non-thermalization of sectors 1 and 2 is due to the low density of dark fermions in sector 1. This can be seen by the yield for the dark fermion in Fig. 3. showing a diminishing Y γ due to a smaller g X for benchmark (d).

Conclusions
In this work we discussed the possibility that DM in the universe is constituted of sub-MeV dark photons which reside in the hidden sector. In this case a proper analysis of the relic density requires a solution to coupled Boltzmann equations which depend on multiple bath temperatures including the bath temperature for the visible sector and those for all the hidden sectors that are feebly coupled with the visible sector. In this work we discussed a model where the visible sector couples with two hidden sectors X 1 and X 2 and where the particles in the hidden sector consist of a dark fermion, a dark Z and a dark photon. The dark Z decays and disappears from the spectrum while the dark photon is a long lived relic. We show that the relic density of the dark photon depends critically on temperatures of both the visible and the hidden sectors. We present exclusion plots where a sub-MeV dark photon can exist consistent with all the current experimental constraints. We also show that the existence of a dark photon is consistent with the constraints on N eff from BBN. Thus a sub-MeV dark photon is a viable candidate for DM within a constrained parameter space of mass and kinetic couplings. The formalism developed here of correlated evolution of bath temperatures in the visible and hidden sectors may find application for a wider class of phenomena involving hidden sectors. A Canonical normalization of extended G SM × U (1) X 1 × U (1) X 2 Lagrangian with kinetic and Stueckelberg mass mixings Consider the Lagrangian C, B, A 3 ). This Lagrangian can be put in the canonical form by appropriate transformations on K E and M 2 . For the case when there is one dark sector it was done analytically in [33] and the basic reason which allows that to happen is that one of the eigenvalues of M 2 is zero corresponding to the photon which effectively reduces the analysis to two massive modes which can be handled analytically. In the present case since we have two hidden sectors, we have a 4×4 matrix, and while one of the eigenvalues corresponding to the photon is zero, one still has to deal with a cubic equation which, although possible to solve analytically, quickly becomes intractable in the presence of both kinetic and Stueckelberg mass mixing. However, because the kinetic mixings are typically small, it is possible to get accurate results by expanding couplings of the dark particles with the SM particles in powers of the kinetic mixings. In this case the relevant couplings can be recovered easily. However, such an expansion must occur around stable minima. This means that we must first diagonalize the SM mass squared matrix for the gauge bosons, and compute the kinetic mixing in this basis. We can then put the kinetic term in the canonical form. This step requires several GL(2, R) transformations because of several mixings of the hidden sector with the visible sector and the mixing of the two hidden sectors. After the kinetic energy is put in the canonical form, we must write the mass square matrix of the gauge bosons in the same basis which then undiagonalizes the said matrix. However, because of the smallness of the kinetic mixings, we can carry out a perturbation expansion of the mass square matrix where the zeroth order mass square matrix is diagonal and the perturbations are proportional to the kinetic mixings and are small. We make this analysis concrete in the formalism below.
Let us consider an orthogonal transformation V = RV (1) D is a diagonal matrix. In the V (1) basis the kinetic energy has the form K E = R T K E R. Next, let us make a transformation K such that V (1) = KV (n) , which could be a product of several sub-transformations, such that the kinetic energy is in the canonical form, i.e., In our case we will have n = 8 (as discussed below). In the V (8) basis, while the kinetic energy is in the canonical form, the mass matrix M 2 = K T M 2 D K is not. However, as explained above since the kinetic mixings are small we can expand M 2 around δ 1 = 0 = δ 2 so that Now since the kinetic mixing is supposed to be small, K differs from a unit matrix only by a small amount and thus ∆M 2 is small relative to M 2 D and one may carry out perturbation expansion in ∆M 2 to arrive at the kinetic and mass mixing effects in the physical processes.
To compute ∆M 2 we need K defined by Eq. (A.2). The computation of K is significantly more complicated than for the case of one hidden sector. Below we give its computation in some detail.
While the procedure outlined above is general we will discuss the specific case where M 2 is block diagonal so that where the upper right 2×2 matrix is for the hidden sector and the lower left 2×2 matrix is for the case of the standard model in the basis with θ w being the weak mixing angle. Here R β diagonalizes the hidden sector mass squared matrix while R w diagonalizes the standard model mass squared matrix, and the diagonalization gives Since the mass square matrix is now diagonal, it is a good starting point to diagonalize and normalize the kinetic energy matrix. This is a bit non-trivial and requires several steps which we outline below. In the basis (D, C, B, A 3 ), the kinetic Lagrangian is given by where we use an abbreviated notation so that B 2 = B µν B µν , BC = B µν C µν , etc. Next we write L KE in the basis (D (1) , C (1) , B (0) , A 3 ) in which the mass square matrix of the gauge bosons is diagonal. Thus, using allows us to write L KE as The diagonal kinetic terms for D (1) and C (1) have the form To normalize them to unity we make a transformation from the basis V (1) T = (D (1) , C (1) , B (0) , A 3 ) to V (2) T = (D (2) , C (2) , B (0) , A 3 ) so that After the transformation, L KE in the V (2) basis has the form We now note that there is a C (2) D (2) mixing term in Eq. (A.13) which can be removed by a GL(2, R) transformation. We do this by going from the basis V (2) 3 ) so that (A. 16) In the V (3) basis, L KE becomes We note that while there are no kinetic mixing terms between C (3) and D (3) , there are kinetic mixing terms between them and the fields B (0) and A 3 . The mixing term between D (3) and B (0) can be removed by the transformation 3 ), s δ 3 and c δ 3 are defined similar to Eq. (A.16) and δ 3 is defined by In the V (4) basis the Lagrangian takes the form A mixing term between D (4) and A In the V (5) basis, L KE takes the form Next we look at the kinetic mixing of (C (3) , B (1) ). This mixing can be removed by the transformation (A. 25) Here V (6) T = (D (5) , C (4) , B (2) , A 3 ) and where s δ 5 and c δ 5 are defined as in Eq. (A. 16). After the transformation the kinetic energy Lagrangian in the V (5) basis has the form 3 B (2) .
Next we examine the kinetic mixing of the fields (C (4) , A 3 ). This mixing term can be eliminated by the transformation where V (7) T = (D (5) , C (5) , B (2) , A 3 ) and where δ 6 is defined by and s δ 6 and c δ 6 are defined as usual.
After the transformation the kinetic energy Lagrangian in the V (7) basis has the form We are now left with the last kinetic mixing term involving the fields B (2) and A 3 . To eliminate this mixing we make the final transformation 3 ). In the basis V (8) the kinetic energy for all the gauge fields is in the canonical form so that The free Lagrangian in the V (8) basis is then As discussed in the beginning of this section we now make the expansion of Eq. (A.3). Below we exhibit K and ∆M 2 in the limit δ 1 , δ 2 1. In this case s δ 1 ∼ δ 1 , c δ 1 ∼ 1, etc. and K and ∆M 2 have the following form The interactions relevant for our computation arise from ∆M 2 and the relation  where z = m Z /m Z and γ = m γ /m Z . The relevant couplings with the visible sector are summarized in Tables 2 and 3. Thus Table 2 gives the couplings of Z and Z to the visible sector fermions ff and Table 3. gives the coupling of γ to ff .
The triple gauge boson couplings of γ , Z , Z are given by are Couplings in the limit of large δ 2 : Some of the processes such as the lifetime of the dark photon require only that the product δ 1 δ 2 be small which could be achieved by δ 1 being small while δ 2 is O(1) size. Thus we list below the vector and axial-vector couplings in the limit of small δ 1 and β while δ 2 is not necessarily small (B.14) (B.15) The couplings with the D fermions become , (B.16) 19) and the triple gauge boson couplings take the form (B.20) (B.22)

C Deduction of three temperature Boltzmann equations
Next we give a deduction of the Boltzmann equations for the case of three heat baths. We will use T 1 as the reference temperature. Let us consider a generic number density n i . In this case n i R 3 is conserved during the expansion if there is no injection and we have d(n i R 3 ) dt = 0 where R is the scale factor, while in the presence of injection one has where C i represent the integrated collision terms. Next, if S = sR 3 is the total entropy, it is conserved which implies that Using Eqs. (C.1) and (C.2), and the fact that n i = sY i , one finds We can convert this equation to one that uses temperature T 1 rather than time which gives where dρ/dT 1 is given by The Boltzmann equations for Y D , Y Z , Y γ may now be written as where C D , C Z and C γ are given by In Eqs. (C.7), (C.8) and (C.9) one encounters thermally averaged decay width and thermally averaged cross sections. The thermally averaged decay width is given by 10) and the thermally averaged cross-section is given by (C.11) K 1 and K 2 are the modified Bessel functions of the second kind and degrees one and two, respectively. For the case when the annihilating particles have different masses m 1 and m 2 and are at different temperatures T 1 and T 2 , the thermally averaged cross-section becomes where and where (C.14) Note that in the limit T 1 → T 2 , I(s) → allows us to recover Eq. (C.11) using Eq. (C.12). The equilibrium yield of species i is given by The hidden sectors degrees of freedom is given by (C. 18) In Eq. (C.18), V = Z , γ and we take g γ = g Z = 3 and g D = 4.

D Dark photon and dark fermion scattering cross sections and Z decay width
Here s is the Mandelstam variable which gives the square of the total energy in the CM system. Further, the notation used above is as follows: f = e, µ, τ : T 3f = −1/2 and where .

Processes
, (D. 16) and c A = a f , c V = v f for V = Z and c A = a f , c V = v f for V = γ .

Process: Z → ff , νν
The decay width of Z to SM fermions is given by