Big Bang initial conditions and self-interacting hidden dark matter

A variety of supergravity and string models involve hidden sectors where the hidden sectors may couple feebly with the visible sectors via a variety of portals. While the coupling of the hidden sector to the visible sector is feeble its coupling to the inflaton is largely unknown. It could couple feebly or with the same strength as the visible sector which would result in either a cold or a hot hidden sector at the end of reheating. These two possibilities could lead to significantly different outcomes for observables. We investigate the thermal evolution of the two sectors in a cosmologically consistent hidden sector dark matter model where the hidden sector and the visible sector are thermally coupled and their thermal evolution occurs without the assumption of separate entropy conservation for each sector. Within this framework we analyze several phenomena to illustrate their dependence on the initial conditions. These include the allowed parameter space of models, dark matter relic density, proton-dark matter cross section, effective massless neutrino species at BBN time, self-interacting dark matter cross-section, where self-interaction occurs via exchange of dark photon, and Sommerfeld enhancement. Finally fits to the velocity dependence of dark matter cross sections from galaxy scales to the scale of galaxy clusters is given. The analysis indicates significant effects of the initial conditions on the observables listed above. The analysis is carried out within the framework where dark matter is constituted of dark fermions and the mediation between the visible and the hidden sector occurs via the exchange of dark photons. The techniques discussed here may have applications for a wider class of hidden sector models using different mediations between the visible and the hidden sectors to explore the impact of Big Bang initial conditions on observable physics.


Introduction
Hidden sectors appear in most modern models of particle physics beyond the standard model and have become increasingly relevant in analyses of particle physics phenomena. Success of precision electroweak physics tell us that the hidden sector couplings to the standard model must be feeble. But what about the coupling of the hidden sector to the inflaton? If the coupling of the hidden sector to the inflaton is also feeble relative to the coupling of the standard model the population of the hidden sector particles would be negligible and their temperature would be much colder than of the standard model particles. On the other extreme the hidden sector and the visible sectors could couple democratically, i.e.,with equal strength, to the inflaton and thus be essentially in thermal equilibrium at the end of reheating. These cases represent two extreme possibilities with a variety of other possibilities in between. Because of the interactions between the visible and the hidden sectors, the two sectors are thermally coupled and thus their evolution is constrained by the initial condition on the hidden sector at the end of inflation which can be codified by the ratio ξ 0 ≡ T 0 h /T 0 , where T 0 h is the temperature of the hidden sector and T 0 is the temperature of the visible sector initially after reheating. It is thus of relevance to ask the influence of the initial conditions on physical observables at low energy. In this work we study the effect of ξ 0 on a variety of physical observables, i.e., on the relic density of dark matter, on the proton-DM scattering cross-sections, on the number of massless degrees of freedom at BBN, and on DM self-interaction cross-sections. For DM self-interaction cross-section, we further analyze the effect of ξ 0 on its velocity dependence and on Sommerfeld enhancement and analyze the effect of ξ 0 on fits to the galactic dark matter cross sections from the scale of dwarf galaxies to the scale of galaxy clusters. The portal we utilize in the analysis consists of a hidden sector with a U (1) X gauge invariance with kinetic mixing [1] and Stueckelberg mass growth of the U (1) X gauge boson [2]. By numerically solving the Schrödinger equation, we are able to achieve a comprehensive understanding of the dark matter self-interacting cross-section in this model. While our analysis is done in a specific choice of the portal one may expect similar effects discussed in this work using other portals connecting visible and the hidden sectors.
The outline of the rest of the note is as follows: Section 2 gives a summary of the thermal evolution of coupled visible and hidden sectors while section 3 discusses a specific model with the visible sector coupled to one hidden sector where the coupling arises via kinetic mixing along with Stueckelberg mass generation for the hidden sector gauge boson. Section 4 discusses the effect of ξ 0 on dark freeze-out and on the relic density of dark matter. Here we also discuss the dependence of ∆N eff at BBN time on ξ 0 and further the effect of ξ 0 on the allowed parameter space and on spin-independent proton-DM cross section. In section 5 we discuss the effect of ξ 0 on Sommerfeld enhancement for the self-interacting dark matter cross section. In section 6 we discuss the effect of ξ 0 on fits to the galaxy data on dark matter cross sections from low relative velocities to high relative velocities which encompass scales from dwarf galaxies to galaxy clusters. In section 7 we give a comparison of our exact analysis of the thermal evolution of hidden sector vs the visible sector with the thermal evolution arising under the approximation that the entropies of the visible and the hidden sectors are separately conserved. Conclusions are given in section 8. Additional details related to the analysis are given in sections 9-10.
2 Cosmologically consistent evolution of coupled visible and hidden sectors in DM analysis.
As mentioned in the previous section, most models of particle physics based on physics beyond the standard model contain hidden sectors which may be feebly coupled to the visible sector. In this case the thermal evolution of each is interdependent on the other. Thus the approximation typically made that the entropy of the visible and the hidden sectors are separately conserved is invalid and it is only the sum of them together that is conserved. Further, the hidden sector by itself may consist of several sectors some of which may interact directly with the visible sector while others indirectly via their interactions with other hidden sectors which couple with the visible sector. First, in this case the Hubble expansion is affected by the hidden sectors via their energy densities so that where ρ v is the energy density of the visible sector and ρ i the energy density of the i-th hidden sector where ρ ′ s have temperature dependence so that (2.2) and the total entropy density of the visible and hidden sectors is given by Here g v,h eff and h v,h eff are the energy and entropy degrees of freedom and are temperature dependent. A full expression for them for the specific model we will consider is given in section 3. In [3] an analysis was given where the visible sector (V) at temperature T is coupled to the hidden sector H 1 at temperature T 1 , the hidden sector H 1 is coupled to the hidden sector H 2 at temperature T 2 , and so on and finally that the hidden sector H n−1 at temperature T n−1 is coupled to the hidden sector H n at temperature T n . In the analysis of [3] radiation dominance was assumed. Here we extend the analysis to include radiation and matter. In this case the energy densities for various sectors obey the following set of coupled Boltzmann equations dρ α dt + 3H(ρ α + p α ) = j α , α = 0, 1, 2, . . . , n , (2.4) Here ρ α and p α are the energy and momentum densities for the sector α, and where α = 0 refers to the visible sector, and α = 1, 2, · · · , n to the hidden sectors, and where j α encodes in it all the possible processes exchanging energy between neighboring sectors. We note now that the total energy density ρ = n α=0 ρ α in an expanding universe satisfies the equation where p = n α=0 ρ α is the total pressure density. In the analysis it is convenient to introduce the functions ζ = 3 4 (1 + p ρ ) and ζ α = 3 4 (1 + pα ρα ), where ζ α = 1 for radiation dominance and ζ α = 3 4 for matter dominance. More generally ζ and ζ α are temperature dependent and this dependence is taken into account in the evolution equations. Thus ρ α and ρ satisfy the evolution equations We use the visible sector temperature T as the clock as thus wish to write the evolution equations Eqs. (2.6) and (2.7) in terms of temperature T rather than time. This is accomplished using the relation Next decomposing ρ so that ρ = ρ v + n i=1 ρ i , one finds that dρ i /dT can be written as (2.10) Note that one may also write (2.14) where ρ ′ v = dρ v /dT . Eqs. (2.14) give us a set of n differential equations for the evolution functions dξ i /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 case of the visible sector coupled to one hidden sector we have The source term j h is discussed in section 9. With this notation specific to the case of the visible sector and one hidden sector we have the following equation for ξ which governs the temperature evolution of the hidden sector relative to that of the visible sector We note that g v eff and h v eff are pre-calculated and we use tabulated results from micrOMEGAs [4]. As noted already g h eff and h h eff for the hidden sector that enter Eq. (2.2) and Eq. (2.3) are temperature dependent [5,6] and their explicit expressions are given in Eq.(3.6).

The model coupling visible and hidden sectors
There are a variety of portals that allow communication between the visible and the hidden sectors. These include the Higgs field portal [7], kinetic mixing of two gauge fields [1], Stueckelberg mass mixing [2,8], kinetic and Stueckelberg mass mixing [9], Higgs-Stueckelberg portal [10], as well as other possibilities such as higher dimensional operators. In this work we focus on kinetic mixing along with the mass growth for the hidden sector gauge field via the Stueckelberg mechanism. Thus for analysis in this work we consider a specific model for dark matter which is an extension of the standard model with an SU (3) × SU (2) × U (1) Y × U (1) X gauge invariance where the U (1) X gauge field has kinetic mixing with the visible sector U (1) Y gauge field [1] and in general a Stueckelberg mass mixings [2,9,11,12]. We assume that the U (1) X hidden sector has a dark fermion D which interacts with the U (1) X gauge field. Thus the extended SU (3) × SU (2) × U (1) Y × U (1) X Lagrangian consisting of the SM part L SM and the extended part L ext is given by Here B µ is the gauge field for the U (1) Y , C µ is the gauge field of U (1) X , σ is an axionic field which gives mass to C µ and is absorbed in the unitary gauge and D is the dark fermion where Q X is the U (1) X charge of D and g X is gauge coupling of U (1) X . Further, δ is the kinetic mixing parameter between the field strengths of C µ and B µ , and M 1 and M 2 are the Stueckelberg mass parameters. A non-vanishing M 2 will lead to a milli-charge for the dark fermion D and we assume neutrality of dark matter and thus set M 2 = 0 in the analysis 1 . The spontaneous breaking of the SU (2)×U (1) Y electroweak symmetry along with the Stueckelberg mass growth gives rise to mixing among the three gauge fields C µ , B µ , A µ 3 where A µ 3 is the third component of the SU (2) L gauge field A µ a (a = 1, 2, 3) of the standard model. The mixings give rise to a 3 × 3 mass square matrix which can be diagonalized by the three Euler angles (θ, ϕ, ψ) which are given in Eq. (3.5). The diagonalization gives the following mass eigenstates: the Z boson, a massive dark photon γ ′ and the massless photon γ. The Lagrangian governing the interaction of dark matter, which is the dark fermion in our case, is given by The interaction of Eq. 3.2 involves two massive gauge bosons (Z, γ ′ ). For the case when the kinetic mixing is small one has g D γ ′ ≃ g X Q X and ϵ D Z = O(δ 3 ) which is negligible. We note here that models with the vector boson as the mediator between the hidden sector and the visible sector have been considered in several previous works . Axions and dark photons in the light to ultralight mass region have also been investigated [45,47,[57][58][59][60][61], and dark photons have been used in explaining astro-physical phenomena including galactic γ-rays [62,63] and PAMELA positron excess [64][65][66][67][68].
We further note that the dark photon in this model even when very light and kinematically disallowed to decay into e + e − will eventually decay via the modes γ ′ → νν and γ ′ → 3γ and not contribute to dark matter density unless its lifetime is larger than the lifetime of the universe and even in that case only if it has non-negligible relic density, which is not the case we consider. Thus the dark fermion will be the only constituent of dark matter.
In addition to the above, in the frame with the canonically diagonalized kinetic energy and mass square matrices in the gauge sector the standard model fermions (i.e., quarks and leptons) have feeble interactions with the dark photon given by where g 2 is the SU (2) L gauge coupling constant, f stands for the standard model fermions and angle θ is defined in Eq. (3.5). The vector and axial vector couplings of the dark photon with the SM fermions f are given by Here s δ = sinh δ, T 3f is the third component of isospin, and Q f is the electric charge for the fermion f . The angles θ and ψ which along with ϕ are the three Euler angles with diagonalize the 3 × 3 gauge boson mass square matrix involving the fields C µ , B µ , A µ 3 are defined as [9] tan ϕ = −s δ , where Further, to compute ζ h = 3 4 (1 + p h ρ h ), we need ρ h and p h , which are given by Here g γ ′ = 3 and g D = 4. For total ζ, we need energy and pressure densities for both the visible and the hidden sectors and we use the relation Before proceeding further we note that the extension to include both the kinetic and mass mixings in Eqs.(3.2)-(3.5) is straightforward as has been discussed in [9] and we exhibit them here for easy reference to guide the discussion in section 7. With inclusion of both the kinetic mixing parameter δ and the mass mixing parameter ϵ = M 2 /M 1 , the neutral current Lagrangian in the hidden sector takes the form There is also a corresponding modification of the neutral current in the visible sector as discussed in [9]. The constraints on δ and ϵ arising from the fits to the electroweak data are mild and one finds that |ϵ − δ| can be as large as 0.05 consistent with the same level of χ 2 fits to the electroweak data as the standard model [9].

Effect of ξ 0 on dark freeze-out and on relic density
In the analysis here we will discuss the effect of ξ 0 on the dark freeze-out which generates the relic density of D. Computationally the quantities of interest for this purpose are the yields for the dark fermion Y D and for the dark photon Y γ ′ , where the yield is defined so that Y = n/s where n is the number density and s is the entropy density. Thus for Y D and Y γ ′ the evolution equations are given by Here ⟨σv⟩ DD→iī is the annihilation cross-section of DD into standard model particles which are denoted by iī, ⟨σv⟩ DD→γ ′ γ ′ is their annihilation into dark photon while ⟨σv⟩ iī→γ ′ gives the annihilation of standard model particles into a dark photon, and n D and n γ ′ are the number densities of the D fermion and the dark photon γ ′ . In the above the cross-section for the process DD → γ ′ γ ′ is given by while the rest of the cross-sections are given in appendix of [13]. Here s, t, u are Mandelstam variables where s + t + u = 2m 2 D + 2m 2 γ ′ , R 11 and R 21 are matrix elements of R which diagonalizes the mass and kinetic energy matrices of Eq.(3.1) as given in [9]. Further, we note that in addition to DD → γ ′ γ ′ we also have γ ′ γ ′ → DD which enters in the yield equations when kinematically allowed and is related to DD → γ ′ γ ′ so that In the above equation, the thermally averaged cross section and decay widths are given by The equilibrium yield of the i-th particle is given by where s is the entropy density. In Eqs.(4.5), (4.6) and (4.7), K 1 and K 2 are the modified Bessel functions of the second kind and of degrees one and of degree two. Further, cross sections for the processes DD → iī, where i,ī are the standard model particles, can be found in Appendix D of [69]. As noted above in this model the dark photon is unstable and decays and does not contribute to the relic density and the entire DM relic density arises from the dark fermion where at current times the relic density Ω D h 2 is given by where s 0 is the current entropy density, Y 0 D which is Y D at current times can be gotten using Eqs. (2.15), (4.1) and (4.2), ρ c is the critical energy density needed to close the universe, and h is defined so that H 0 = 100h km s −1 Mpc −1 , where H 0 is the Hubble parameter today.
The procedure for solving the evolution equations involves simultaneous analysis of coupled equations Eq.  Table 1 gives 6 model points used in this paper all of which are consistent with the current experimental constraints [70] including those from a variety of experiments, i.e., BaBar [71,72], HPS [73], LHCb [74], Belle-2 [75], SHiP [76], SeaQuest [77,78] and NA62 [78], CHARM [78], νCal [78][79][80], E137 [81], E141 [82], NA64 [83], NA48 [84]. For a sub-MeV dark photon mass stringent constraints arise from Supernova SN1987A [85] and from BBN, stellar cooling [86] and from the decay to 3γ on cosmological timescales [30,87] An analysis of these constraints in limiting the parameter space have been analyzed in [69,70]. The parameter space chosen in the current analysis is consistent with these constrains.  We note here that the mass of the dark photon is in the sub MeV region and is long lived with its most dominant decay mode being γ ′ → 3γ. For kinetic mixing the decay width for the mode γ ′ → 3γ is given by [30,87,88] (4.9) where α = e 2 4π , ϵ γ γ ′ is the kinetic mixing parameter of coupling between dark photon γ ′ and photon γ given by ϵ γ γ ′ = g Y √ 1 +ε 2 R 21 as defined in [9], m γ ′ is the dark photon mass and m e is the electron mass. The dark photon lifetimes for the different model points are given in Table (1). Here we find that the dark photon lifetimes are smaller than the age of the universe and thus there is no contribution of the dark photon to the relic density and consequently no constraint on the allowed parameter space regarding the relic density constraint. We note, however, that even if the dark photon was long lived with a lifetime greater than the lifetime of the universe, its contribution to the relic density would be negligible. A recent analysis [69] in accord with the analysis of [30] shows that with one hidden sector it is not possible to get both a long lived dark photon which can contribute to the relic density and simultaneously achieve a significant amount of dark matter relic density. To do that one needs at least a two hidden sector model [69] in which a dark photon as dark matter can produce a non-negligible amount of dark matter.
10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 1.  We discuss now the dependence of the dark matter freeze-out and of the relic density on the initial conditions. In the top left panel of Fig. 1 we exhibit the dependence of the dark freeze-out and specifically the decoupling of the dark photon and the dark fermion on ξ 0 where we consider the cases: ξ 0 = 0.01 and ξ 0 = 1. The top right panel is the zoom in of the top right in the region of the freeze-out. From Fig. 1 we see that the process Model   Fig. 2 we exhibit the dependence of the relic density on ξ 0 for the six model points of Table 1. Here we find that the relic density can change up to 40% as ξ 0 varies in the range (0, 1).

Dependence of ∆N eff at BBN on ξ 0
One of the predictions of beyond the standard model physics is N eff , the number of effective relativistic degrees of freedom at BBN. For the standard model N eff = 3.046. The current experimental constraint on N eff is summarized in Fig. 39 of the Planck Collaboration [89] which shows the spread in N eff . Thus the Planck Collaboration gives N eff = 2.99 ± 0.17 while the joint BBN analysis of deuterium/helium abundance and the Planck CMB data gives N eff = 3.41 ± 0.45. Here we will use the conservative constraint on ∆N eff = N exp eff − N sm eff so that ∆N eff ≤ 0.25. In the model under discussion, the dark fermion D and dark photon γ ′ will contribute to the effective neutrino number. Such contribution is given by where g h eff can computed from Eq.(3.6) and Eq.(4.10) is to be evaluated at the BBN temperature T BBN = 1 MeV (for related works see, e.g., [90,91]). In Table 2, ∆N eff is computed for the six model points of Table 1 for ξ 0 = 0.01 and ξ 0 = 1 while the right panel of Fig.  2 exhibits ∆N eff (BBN) for the six model points for ξ 0 in the range (0 − 1). The analysis in general indicates that hidden sectors which start off cooler than the standard model at the end of reheating contribute a smaller amount to ∆N eff than those which are relatively hotter at the end of reheating. Further, the analysis indicates that models where ξ ≃ 0 could accommodate more massless degrees of freedom allowing for the possibility of building a wider class of models with more hidden sector particles which may still be consistent with the ∆N eff constraint at BBN time.  4.3 Effect of ξ 0 on the allowed parameter space and on spinindependent proton-DM cross section Next, we investigate the influence of ξ 0 on the allowed parameter space consistent for a chosen range of relic density. To this end we constrain the relic density to lie in the range 0.012 ≤ Ωh 2 ≤ 0.12 and m D to lie in the range of 5 GeV to 10 TeV. Specifically we explore the allowed region for the two cases: ξ 0 = 0.01 and ξ 0 = 1. The result of our analysis is exhibited in the left panel of Fig. 3 which gives a scatter plot of the allowed models in δ vs m γ ′ where those with color blue correspond to ξ 0 = 0.01 and those with color red correspond to ξ 0 = 1.
One of the interesting result that emerges is that for ξ 0 = 1 most of the models lie in the range 10 −10 < δ < 10 −5 while for ξ 0 = 0.01 the allowed range is 10 −9 < δ < 10 −4 . Thus, the analysis shows that the initial choice of ξ 0 significantly impacts the model's allowed parameter space. ξ 0 also has significant effect on the proton-DM scattering cross section in the direct detection experiments for dark matter. Specifically we consider the spin independent proton-DM cross section σ SI:p−DM . Here we use the micrOMEGAs [4] to find the spin independent cross section. In the right panel of Fig. 3 we exhibit σ SI:p−DM for the six model points of Table 1 and their dependence on ξ 0 in the range (0.01-1) is indicated by the small vertical lines for each of the model points. The numerical values of the σ SI:p−DM for ξ 0 = 0.01 and ξ 0 = 1 are exhibited in Table 3 for the models of Table 1. Here one finds that the variation of the cross-section can be as large as 40%. Thus some of the models that are eliminated for the ξ 0 = 1 case would still be viable for the case ξ 0 = 0.01. Model  Table 3: Table of spin-independent proton-DM cross section σ SI:p−DM for the model points of Table 1 for ξ 0 = 0.01 and ξ 0 = 1.0.
5 Self-interacting dark matter, Sommerfeld enhancement, and dependence on ξ 0 The self-interacting dark matter cross sections arises from the processess DD → DD, DD → DD, andDD →DD via the exchange of a dark photon. The Lagrangian of Eq. (3.1) leads to a Yukawa potential between the D-fermions due to the dark photon exchange in the non-relativistic limit so that where the plus sign is for DD → DD andDD →DD and the minus sign is for DD → DD.
In some of the regions of parameters (i.e. m γ ′ m D ≤ (g X ) 2 4π ), tree-level scattering or the Born approximation is no longer valid and one has contributions from higher order dark photon  exchanges as shown in Fig. 4 which contribute to scattering. In this case we need to numerically solve the Schrodinger equation to find the accurate scattering cross sections. The radial equation one needs to solve is given by where p is the particle momentum and V (r) is the potential. The substitution x = pr and R p,l = N pΦ l (x)/x gives [94] The non-perturbative effect arising from the repeated exchange of the mediator is often encoded in Sommerfeld enhancement and has been discussed in several previous works (see, e.g., [24,36,[95][96][97][98] and the references therein). Thus including non-perturbative effects the annihilation cross section times the velocity v (where v is the relative velocity in the CM system) for the cross section σ ab for the elastic scattering process a + b → a + b may be written as where (σ 0 ab v) is the tree level cross section and S E is the Sommerfeld enhancement. As noted in the present context the contribution to the Sommerfeld enhancement arises from multiple exchanges of the dark photon γ ′ . The solution of the differential equation Eq.(5.3) has the form: where δ l is the phase shift for the l−th partial wave. We write the Sommerfeld enhancement of l-th partial wave cross-section for the case of the Yukawa potential so that where [94], S E l = (1 · 3 · · · (2l + 1)/A) 2 . Using Eq. (5.5), we get Taking x larger than 30 gives a good enough approximation to the exact solution. Typically an attractive potential leads to Sommerfeld enhancement of cross section at low collision velocities, but one may also have Sommerfeld suppression for a repulsive potential. In the where S E0 is for the Yukawa potential by solving numerically (red), S C0 is for the Coulomb potential (black), and S H0 is for the Hulthen potential as a function of g X for the case when m D = 250 GeV, m γ ′ = 45 MeV, and v = 10 km/s. Right panel: Plot of S-wave Sommerfeld suppression for the case of a repulsive potential where we use the same symbols S E0 , S H0 , S C0 for suppression as for enhancement to avoid a proliferation of notation.
left panel of Fig. 5 we exhibit Sommerfeld enhancement for the case of a negative Yukawa potential. Here we see that Sommerfeld enhancement can be very significant and further the enhancement shows oscillatory behavior with g X . To check the accuracy of our numerical analysis and to explain the oscillatory behavior we compare our result with those from the Hulthen potential as an approximation to the Yukawa potential for which one can obtain a good analytic approximation for the S-wave. The Hulthen potential is given by [99,100] It is known that Hulthen potential is a very good approximation to Yukawa potential both at short and at long distances. With it one can find an analytic solution for the S-wave and thus find a good analytic approximation to the S-wave Sommerfeld enhancement [101,102]: From Eq.(5.9), valid for the attractive potential case, it is obvious that the oscillation is due to the existence of the cosine term. For the Coulomb potential the Sommerfeld enhancement for the S-wave is given by where plus is for repulsive potential and minus is for attractive potential. The left panel of Fig.5 gives a comparison of the S-wave Sommerfeld effect for three different potentials: Yukawa, Hulthen and Coulomb. The analysis shows that Hulthen potential gives a good approximation to the Yukawa potential and also explains the deep oscillations as a function of g X . For the case of a repulsive potential (α negative), the analysis is very different. A comparison of the numerical analysis using Yukawa potential and the analytic solution using Hulthen potential for the case of a repulsive potential is given in the right panel of Fig. 5.
Here again one finds that the numerical analysis and the Hulthen potential result fully agree.
Having checked the numerical accuracy of our analysis in Fig. 5 we next investigate the effect of the Big Bang initial conditions on Sommerfeld enhancement. In Fig. 6 we compare S-wave Sommerfeld enhancement for the cases ξ 0 = 0.01 and ξ 0 = 1 for the case of an attractive Yukawa potential. The left panel of Fig. 6 shows Sommerfeld enhancement vs v and here one finds that ξ 0 = 1 (red) gives an enhancement which is larger than for the case ξ 0 = 0.01 (blue). In the analysis we keep the relic density fixed at ∼ 0.12 for ξ 0 = 0.01 and ξ 0 =1 by allowing g X to vary. The right panel of Fig. 6 displays Sommerfeld enhancement as a function of m γ ′ /m D and here one finds that the oscillation peaks for the case ξ 0 = 1 (red) are significantly larger than those for the case ξ 0 = 0.01 (blue). A similar analysis for a repulsive Yukawa potential is carried out in Fig. 7. However, in this case we have Sommerfeld suppression rather than an enhancement where the Sommerfeld suppression is vs v for the left panel and vs m γ ′ /m D for the right panel and the red curve is for ξ 0 = 1.0 and the blue curve for ξ 0 = 0.01. For both cases the Sommerfeld suppression is significantly larger for ξ 0 = 1 relative to ξ 0 = 0.01.  6 Effect of ξ 0 on fit to galaxy data Several analyses of galaxy data indicate that dark matter is collisional at the scale of dwarf galaxies and appears collision-less at the scale of galaxy clusters [33,38,103]. Thus, for dwarf galaxies one finds collisional velocity ⟨v⟩ of dark matter in the range 10-100 km/s and 1 cm 2 /g< σ/m < 50 cm 2 /g [33,103] where σ is the cross section and m is the mass of DM particle. For midsize galaxies such as the low surface brightness galaxies (LSB) and the Milky Way one finds ⟨v⟩ in the range 80-200 km/s and 0.5 cm 2 /g< σ/m < 5 cm 2 /g. The galaxy clusters exhibit ⟨v⟩ > 1000 km/s. Here it is estimated that the σ/m is maximally 1 cm 2 /g [33,103,104] and could be as low as 0.065 cm 2 /g< σ/m < 1 cm 2 /g [38,105,106].
As is well known one interesting possibility to account for the velocity dependence of the DM cross sections is that DM is self-interacting by Spergel and Steinhardt [107] and there is considerable follow up work on this idea [22, 25-28, 36, 108-115]. An analysis of fit to the data within the dark photon model was previously done in [11] (see also [116,117]).
Here we study the dependence of the fits on ξ 0 . Further, here the analysis goes beyond the Born approximation used in [11] taking into account non-perturbative effects encoded in the Sommerfeld enhancement with also inclusion of identical particle exchange effects. Since the dark matter is constituted of dark Dirac fermions consisting of D andD constituents, we will have processes of the type DD → DD, DD → DD andDD →DD. Thus the total cross section σ DM is given by where the factor of 1/2 arises due to identical nature of particles.
To numerically calculate the cross section, we start with Eq.(5.3) and use the method of [27].
Here in the computation of DM cross sections, we need to calculate the phase shifts (δ ℓ ) for DD → DD separately from the phase shifts (δ ′ ℓ ) for DD → DD while the phase shifts for the processDD →DD will be the same as for the process DD → DD. Including all contributions, i.e., from DD → DD, DD → DD andDD → D¯D, and taking account of the identical nature of particles in DD → DD andDD → D¯D scattering we find where f ℓ = e iδ ℓ sin δ ℓ /k and f ′ ℓ = e iδ ′ ℓ sin δ ′ ℓ /k. The details leading to Eq.(6.2) are given in section 10. The result of our numerical analysis to fit the galaxy data on σv/m D in the range of velocities from 10 km/s to 10 4 km/s is given in Fig. 8 which exhibits the dependence of the fits on ξ 0 in the range ξ 0 = 0.01 to ξ 0 = 1. In the analysis we allow g X to vary to keep the relic density fixed at Ω D h 2 ∼ 0.12 as ξ 0 varies between 0.01 and 1. The analysis shows that the variation of σv/m with ξ 0 is significant and can sometimes be as large as O(1) (see Model (f)) in Fig. 8. We note that the plots include Sommerfeld enhancement effects but these effects are relatively small. The reason for it is that the Sommerfeld enhancement strongly depends on g X as can be seen from the left panel of Fig. 5. However, in the analysis of galaxy fits of Fig. 8, we find that g X is relatively small which suppresses the Sommerfeld enhancement.
Our result here is consistent with a similar observation on Sommerfeld enhancement in the work of [36] (see also [24]). In Fig. 8 the Born approximation results are also plotted for comparison with the exact solutions. Further, we note that more fine tuned fits to the galaxy data can be gotten by adjustment of the model parameters such that resonances appear in some of the low lying partial waves, e.g., S, P and D waves. This is exhibited in Fig. 9 where in the left panel we see enhancements in the S and the P waves appear to simulate the oscillations in the data at ⟨v⟩ ∼ 10 2 km/s and at ⟨v⟩ ∼ 10 3 km/s. On the right panel of Fig. 9, we plot the cross section contributed from each partial wave separately. It is clear that the peak at ⟨v⟩ ∼ 10 2 km/s is largely due to the S-wave while the one at ⟨v⟩ ∼ 10 3 km/s has a large contribution from l = 5 although sum of all partial waves up to l = 5 enter in the fit given on the left panel.     Figure 8: A fit to the galaxy data taken from [33] which studies the dependence of σv/m D on ξ 0 in the range (0.01 − 1) for the six models of Table 1. Here solid lines are for ξ 0 = 1 and the dashed line for ξ 0 = 0.01 exhibiting the dependence of σv/m D on ξ 0 . The fits (in red) are done using the full analysis by numerically integrating the Schrodinger equation including identical particle effects as well as Sommerfeld enhancement. For comparison we also exhibit the tree-level QFT cross section shown by black curves which does not consider the effect of identical scattering. In the analysis we allow g X to vary but keep the relic density fixed at ∼ 0.12 as ξ 0 varies.  Figure 9: Left: Exhibition of peaks in fits to the galaxy data with specific parameters. Right panel: Here it is shown how the peaks in the left panel at ⟨v⟩ ∼ 10 2 km/s and at ⟨v⟩ ∼ 10 3 km/s arise from successive additional of higher waves. Thus the peak at ⟨v⟩ ∼ 10 2 km/s arises mainly from S and P contributions, while the one at ⟨v⟩ ∼ 10 3 km/s arises from contributions from up to l = 5.
Besides the total cross section, transfer cross section [22,23,[25][26][27][28]118] (σ T ) which have been used in the simulations of long range interactions [108,111,119]. Further, viscosity cross sections (σ V ) are also widely considered in analysis of SIDM [27,120,121]. They are defined so that In terms of partial waves σ T and σ V are given by Details of their computation in terms of partial waves are given in section 10. Fig. 10 shows that σ tot , σ T , σ V differ significantly from each other. We also exhibit the tree level cross section for comparison.  Figure 10: A comparison of σv/m D vs v for different cross sections which include the total cross section σ total , the transverse cross section σ T , the viscosity cross section σ V as defined by Eq. (6.3) and the tree level QFT cross section with and without identical particle effects for model (f) of Table 1. The analysis is done for ξ 0 = 1.

On the validity of entropy conservation approximation
In several previous works (see, e.g., [29]) an assumption of entropy conservation per comoving volume separately for the visible and the hidden sectors is made to relate ξ(T ) at different temperatures. The above implies that the ratio s h /s v is unchanged at different temperatures where s v and s h are the entropy densities for the visible and the hidden sectors where Specifically it is assumed that the following relation between the temperatures T 0 and T holds Noting that T h = ξ(T )T and T 0h = ξ 0 T 0 where ξ 0 ≡ ξ(T 0 ) we can write the above equation as follows Note that the left hand side is a highly non-linear function of ξ(T ) since for our model In Fig. 11 we give a comparison of the exact analysis of the evolution of the ξ(T ) vs the evolution given by the approximation of entropy conservation in co-moving volume for the visible and the hidden sectors separately. The analysis shows that as ξ 0 gets progressively smaller deviations of the approximate solutions gets progressively worse and especially in the freeze-in region where ξ 0 = 0.001, the deviations of the approximate from the exact is huge for temperatures in the visible sector below 10 5 GeV. More importantly for any choice of ξ 0 in the range (0, 1) which includes both the freeze-out and the freeze-in regions, the prediction of ξ 0 for the approximation is always inaccurate at the BBN temperature of ∼ 1 MeV. The right panel gives a plot of ξ as a function of the visible sector temperature for different values of δ for the case when ξ 0 = 0.001. Here one finds that the approximation (dashed line) gives reasonably accurate result for the case when δ = 0, i.e., there is no kinetic mixing but it gives highly inaccurate results for the case when δ in non-vanishing, even as small as δ ∼ 10 −10 . Thus, our conclusion, is that the entropy conservation approximation in thermal evolution is not suitable for any precision analysis. In the analysis done so far we assumed M 2 = 0. For generality we consider now the case where we include the mass mixing parameter ϵ along with kinetic mixing δ. Thus we discuss again the thermal evolution when there are both kinetic mixing and mass mixing present where we use the relations given by Eq.(3.9) and Eq.(3.10) and related relations given in [9]. In Fig. 12 we investigate the effect of including ϵ along δ on the evolution of ξ(T ). The left panel is for the case ϵ = 0.9δ with δ = 4 × 10 −8 and as expected the evolution for different ξ 0 shows a pattern similar to the left panel of Fig. 11. The right panel of Fig. 12 shows that evolutions with different ϵ follow a similar path at high temperatures but begin to diverge at T ∼ 10 GeV. This divergence results in signficantly different values of ξ(T ) at BBN temperature. As expected we find that since δ and ϵ together control the thermal evolution, there is a significant difference in the pattern of evolution relative to the right hand panel of Fig. 11. However, for both Fig. 11 and Fig. 12 one finds that the predictions for ξ(T ) given by the entropy approximation equation Eq.(7.2) shown by dashed curves differs significantly from the result of exact analysis using Eq.(2.15) over wide regions and specifically at BBN temperature.

Conclusion
Hidden sectors are ubiquitous in models of extra dimensions, in extended supergravity and in strings and appear in a variety of beyond the standard model constructions such as in moose/quiver gauge theories (see, e.g., [122][123][124][125]). While the hidden sectors are neutral under the standard model gauge group they can couple feebly with the standard model. However, the couplings of the hidden sector to the inflaton could vary over a wide range. Thus on one extreme the hidden sector coupling to the inflaton could be negligible relative to the coupling of the standard model. In this case at the end of reheating there would be essentially no production of the hidden sector particles, except via gravitational production, and the hidden sector would likely be colder than the standard model. On the other extreme the hidden sector and the visible sectors could couple democratically to the inflaton and in this case the hidden sector and the visible sectors would be in thermal equilibrium at the end of reheating. These two extremes would have significantly different thermal evolution and would result in significant differences in their predictions of the physical observables. In this work we have investigated these effects in the context of a specific hidden sector model which arises from a U (1) X extension of the standard model gauge group. The contents of the hidden sector consists of a dark fermion which has gauge interactions with the U (1) X gauge field. The communication between the hidden sector and the visible sector arises from kinetic mixing between the U (1) X and U (1) Y gauge fields where the U (1) X gauge field acquires mass via the Stueckelberg mechanism. In view of the asymmetric coupling of the visible and the hidden sectors to the inflation field, the temperature of the hidden sector T 0 h and of the visible sector T 0 will in general be unequal at the end of reheating. Thus the ratio of the two, i.e., ξ 0 = (T 0 h /T 0 ) enters in the thermal evolution of the hidden and the visible sectors and affects phenomena at low energy.
The analysis of the work provides a cosmologically consistent framework in that it involves a synchronous evolution of the coupled hidden and visible sectors with only the total entropy of the two sectors being conserved. In the above framework we investigate a number of phenomena and their dependence on ξ 0 . These include dark freeze-out, relic density and the extra number of relativistic degrees of freedom at the BBN time, and the proton-DM cross-section. Further, we investigate the effects of ξ 0 on the self-interaction cross section and on Sommerfeld enhancement. The model is then used in fitting self-interacting dark matter cross sections from galaxy scales to the scale of galaxy clusters. Here we find that fits to data show a significant variation sometime as much as O(1) for ξ 0 in the range (0, 1). Thus the analysis indicates that inclusion of hidden sectors which appear in a variety of models of particle physics beyond the standard model, the initial constraints on the hidden sector at the end of reheating and specifically on ξ 0 could have significant influence on observables and thus their inclusion will be relevant for accurate description of physical phenomena. While our analysis is done for the case of one portal, the general techniques discussed here would be valid for a broader class of models. Finally, we show that the approximation often made in the thermal evolution of visible and hidden sectors by assuming entropy conservation for each of the sectors separately gives widely inaccurate results even for the case for very feeble interactions such as, for example, with the kinetic mixing parameter as low as δ = 10 −10 . Such an analysis is thus a poor approximation to the exact analysis we carry out for the thermal evolution of the visible and the hidden sectors in a synchronous manner using Eq.(2.15). For generality we also consider the case with the mass mixing parameter ϵ which reaches a similar conclusion. Thus an accurate thermal evolution is essential for the current and future precision analyses in cosmology while analyzing physics involving hidden sectors.

Appendix A: Source functions
The source term j h that appears in Eq.(2.15) is defined by:  10 Appendix B: Self-interacting dark matter cross sections 10.1 DD → DD Let us first consider DD → DD scattering. Here the wave function for scattering of a plane wave scattering from a central potential is given by ψ(⃗ r) ∼ [e ikz + f (θ)e ikr /r]. In this case the scattering amplitude f (θ) has an expansion in terms of the partial wave amplitudes f l = e iδ l sin δ l k where δ l is the phase shift for the l-th partial wave so that f (θ) has the expansion f (θ) = ∞ l=0 (2l + 1)f l P l (cos θ) and σ DD tot have the expansion σ DD total = |f (θ)| 2 dΩ = (4π) l (2l + 1)|f l | 2 . can be expanded in terms of partial waves using the relations for l = l ′ − 2 ,

DD → DD,DD →DD
Here the scattering involves identical particles which are fermions so the overall wave function for the particles must be anti-symmetric. This can happen in two ways: (i) spin antisymmetric and space symmetric or (ii) spin symmetric and space anti-symmetric. Now the two spin particles can have a total spin 1 (triplet state) or total spin is zero (singlet state). For the triplet state the space wave function must be anti-symmetric in θ → π − θ and for the singlet state the space wave function must be symmetric. Thus we have for σ DD (θ) the expression Further, σ DD tot is given by where the front factor of 1/2 is to take account of the identical nature of the scattering particles. The partial wave analysis of σ DD tot gives σ DD tot = 2π l (2l + 1) |f l | 2 + 2 1 − 1 2 (−1) l |f ′ l | 2 . (10.10) Here f ′ l = e iδ ′ l sin δ ′ l since the potential governing DD → DD scattering is different from the one that governs DD → DD scattering. Similar calculations are done for transfer cross section and for viscosity cross section. Further, σDD tot = σ DD tot .