Coscattering/coannihilation dark matter in a fraternal twin Higgs model

Dark matter candidates arise naturally in many models that address the hierarchy problem. In the fraternal twin Higgs model which could explain the absence of the new physics signals at the Large Hadron Collider (LHC), there are several viable dark matter candidates. In this paper we study the twin neutrino in the mass range ∼ 0.1–10 GeV as the dark matter. The thermal relic density is determined by the interplay of several annihilation and scattering processes between the twin neutrino, twin tau, and twin photon, depending on the order of the freeze-out temperatures of these processes. Besides the common coannihilation scenario where the relic density is controlled by the twin tau annihilation, it can realize the recently discovered coscattering phase if the scattering of the twin neutrino into the twin tau freezes out earlier than the twin tau annihilation. We also provide a method to calculate the thermal relic density in the intermediate regime where both coannihilation and coscattering processes contribute to the determination of the dark matter density. We show that the right amount of dark matter can be obtained in various scenarios in different regions of the parameter space. The current experimental constraints and future probes into the parameter space from direct detections, cosmological and astrophysical bounds, dark photon searches, and displaced decays at colliders, are discussed.


Introduction
The hierarchy problem and the dark matter are two main motivations for new physics near the electroweak (EW) scale. In the standard model (SM), the Higgs field receives large quadratically divergent contributions to its potential from the interactions with SM particles, in particular, the top quark and weak gauge bosons. For a natural EW symmetry breaking scale, new particles are expected to be close to the EW scale to cut off these quadratic contributions. On the other hand, a stable weakly interacting massive particle (WIMP) with a mass around the EW scale gives a right amount of thermal relic from the Hot Big Bang to account for the dark matter in the universe. It is called "WIMP miracle." Such a dark matter particle candidate also often appears naturally in models -1 -

JHEP09(2018)098
which address the hierarchy problem. The most popular and most studied examples are the supersymmetric (SUSY) extensions of SM. With a conserved R-parity, the lightest neutralino is stable and represents a good dark matter candidate. The models that can explain both the hierarchy problem and dark matter are particularly attractive because they provide a link between the two mysterious problems.
The new particles related to the hierarchy problem and the WIMPs have been extensively searched at various experiments. So far none of them has been discovered. The LHC has put very strong bounds on new colored particles that can cancel the SM top loop contribution to the Higgs mass. Except for some special cases, the bounds on the masses of new colored particles generically exceed 1 TeV. This would imply a quite severe tuning of the Higgs mass if the top loop is not canceled below 1 TeV. Direct searches of DM also put strong bounds on the scattering cross sections of the DM particle with nucleons. A big fraction of the expected region of the WIMP parameter space from typical SUSY models is excluded, though there are still surviving scenarios. These null experimental results have prompted people to wonder that the standard pictures such as SUSY might not be realized at the electroweak scale in nature. Alternative solutions to the hierarchy problem and DM where the interactions between new particles and SM particles are stealthier should be taken more seriously.
For the hierarchy problem, the "neutral naturalness" models gained increasing attentions in recent years. In these models, the top quark partners which regularize the top loop contribution to the Higgs mass do not carry SM color quantum numbers, and hence are not subject to the strong bounds from the LHC. The mirror twin Higgs model [1] is the first example and is probably the stealthiest one. The twin sector particles are SM singlets but charged under their own SU(3) × SU(2)[×U (1)] gauge group. They are related to the SM sector by a Z 2 symmetry. As a result, the mass terms of the Higgs fields of the SM and twin sectors exhibit an enhanced SU(4) symmetry. The 125 GeV Higgs boson arises as a pseudo-Nambu-Goldstone boson (PNGB) of the spontaneously broken SU(4) symmetry. The twin sector particles are difficult to produce at colliders because they do not couple to SM gauge fields. The main experimental constraints come from the mixing between the SM Higgs and the twin Higgs, which are rather weak. The model can still be natural without violating current experimental bounds.
The next question is whether the neutral naturalness models like the twin Higgs possess good dark matter candidates. In the fraternal version of the twin Higgs model [2], people have shown that there are several possible dark matter candidates [3][4][5][6]. The fraternal twin Higgs takes a minimal approach in addressing the hierarchy problem using the twin Higgs mechanism. In this model, the twin fermion sector only contains the twin partners of the third generation SM fermions, since only the top Yukawa coupling gives a large contribution to the Higgs mass that needs to be regularized below the TeV scale. The twin U(1) gauge boson can be absent or can have a mass without affecting the naturalness. The fraternal twin Higgs model can avoid potential cosmological problems of an exact mirror twin Higgs model which contains many light or massless particles in the twin sector. Refs. [3,4] showed that the twin tau can be a viable dark matter candidate. The correct relic density is obtained for a twin tau mass in the range of 50-150 GeV, depending on other model parameters. If twin hypercharge is gauged, then the preferred mass is lighter, in the 1-20 GeV range [4]. Another possibility is asymmetric dark matter from the twin baryon made of twin b-quarks, where the relic density is set by the baryon asymmetry in the twin sector [5,6].
In this paper, we explore a new scenario where the dark matter is the twin neutrino, in the mass range ∼ 0.1-10 GeV. In previous studies of twin tau dark matter, its stability is protected by the twin U(1) EM symmetry, which is assumed to be a good symmetry, either gauged or global. Here we consider that the twin U(1) EM is broken so that the twin photon acquires a mass to avoid potential cosmological problems. In this case, the twin tau and the twin neutrino can mix so the twin tau can decay to the twin neutrino if the twin tau is heavier. The twin neutrino, on the other hand, being the lightest twin fermion, can be stable due to the conservation of the twin lepton number or twin lepton parity. An interesting scenario is that if the twin photon, twin neutrino, and twin tau all have masses of the same order around a few GeV or below, the right amount of dark matter relic density from twin neutrinos can be obtained. The relic density is controlled by the coannihilation and recently discovered coscattering processes [7,8].
The coscattering phase is considered as the fourth exception in the calculation of thermal relic abundances in addition to the three classical cases enumerated in ref. [9]. It is closely related to the coannihilation case as both require another state with mass not far from the dark matter mass so that the partner state can play an important role during decoupling. The difference is that in the coannihilation phase the relic density is controlled by the freeze-out of the annihilation processes of these particles, while in coscattering phase the relic density is controlled by the freeze-out of the inelastic scattering of a dark matter particle into the partner state. Because of the energy threshold of the upward scattering, the coscattering process has a strong momentum dependence. This makes the relic density calculation quite complicated. The standard DM calculation tools such as micrOMEGAs [10], DarkSUSY [11], and MadDM [12] do not apply and one needs to solve the momentum-dependent Boltzmann equations. Also, because of the momentum dependence of the coscattering process, we find that there are parameter regions of mixed phase, i.e., the relic abundance is controlled partially by coannihilation and partially by coscattering. We investigate in detail the relevant parameter space and perform calculations of the relic abundance in different phases, including situations where it is controlled by coannihilation, by coscattering, or by both processes. The calculation in the mixed phase is more involved and we discuss a relatively simple method to obtain the DM abundance with good accuracies.
The twin neutrino DM does not couple to SM directly. Its couplings to matter through mixings of the photons or the Higgses between the SM sector and the twin sector are suppressed, so the direct detection experiments have limited sensitivities. Some of the main constraints come from indirect detections and searches of other associated particles. Its annihilation through twin photon is constrained by Cosmic Microwave Background (CMB), 21cm line absorption, and Fermi-LAT data. For associated particles, the light dark photon searches provide some important constraints and future probes.
The paper is organized as follows. In section 2 we first give a brief summary of the fraternal twin Higgs model and its possible DM candidates. Then we focus the discussion -3 -

JHEP09(2018)098
on the sector of our DM scenario, i.e., the twin neutrino as the DM, and its coannihilation/coscattering partners, the twin tau and the twin photon. In section 3 we enumerate the relevant processes and discuss their roles in controlling the DM abundance in different scenarios. Section 4 describes how to evaluate the DM relic density in different phases, including the coannihilation phase, the coscattering phase, and the mixed phase. Some details about the calculations are collected in the appendices. Our numerical results for some benchmark models are presented in section 5. In section 6 we discuss various experimental constraints and future probes of this DM scenario. The conclusions are drawn in section 7.

Fraternal twin Higgs and light DM
The twin Higgs model postulates a mirror (or twin) sector which is related to the SM sector by a Z 2 symmetry. The particles in the twin sector are completely neutral under the SM gauge group but charged under the twin SU(3) × SU(2) × U(1) gauge group. Due to the Z 2 symmetry, the Higgs fields of the SM sector and the twin sector exhibit an approximate U(4) (or O(8)) symmetry. The U(4) symmetry is spontaneously broken down to U(3) by the Higgs vacuum expectation values (VEVs). A phenomenological viable model requires that the twin Higgs VEV f to be much larger than the SM Higgs VEV v, f /v 3, so that the light uneaten pseudo-Nambu-Goldstone boson (pNGB) is an SM-like Higgs boson. (The other six Nambu-Goldstone bosons are eaten and become the longitudinal modes of the W, Z bosons of the SM sector and the twin sector.) This requires a small breaking of the Z 2 symmetry. The one-loop quadratically divergent contribution from the SM particles to the Higgs potential is cancelled by the twin sector particles, which are heavier than their SM counterparts by the factor of f /v. The model can be relatively natural for f /v ∼ 3 − 5.
The twin sector particles are not charged under the SM gauge group so it is difficult to look for them at colliders. However, if there is an exact mirror content of the SM sector and the couplings respect the Z 2 symmetry, there will be light particles (photon, electron, neutrinos) in the twin sector. They can cause cosmological problems by giving a too big contribution to N eff . 1 In addition, in general one expects a kinetic mixing between two U(1) gauge fields. If the twin photon is massless, its kinetic mixing with the SM photon is strongly constrained. On the other hand, these light particles have small couplings to the Higgs and hence play no important roles in the hierarchy problem. One can take a minimal approach to avoid these light particles by only requiring the Z 2 symmetry on the parts which are most relevant for the hierarchy problem. This is the fraternal twin Higgs (FTH) model proposed in ref. [2]. The twin sector of the FTH model can be summarized below.
• The twin SU(2) and SU(3) gauge couplings should be approximately equal to the SM SU(2) and SU(3) gauge couplings. The twin hypercharge does not need to be gauged. If it is gauged, its coupling can be different from the SM hypercharge coupling, as long as it is not too big to significantly affect the Higgs potential. Also, the twin photon can be massive by spontaneously breaking the U(1) gauge symmetry or simply writing down a Stueckelberg mass term.

JHEP09(2018)098
• There is a twin Higgs doublet. Together with the SM Higgs doublet there is an approximate U(4)-invariant potential. The twin Higgs doublet acquires a VEV f v, giving masses to the twin weak gauge bosons and twin fermions.
• The twin fermion sector contains the third generation fermions only. The twin top Yukawa coupling needs to be equal to the SM top Yukawa coupling to a very good approximation so that their contributions to the Higgs potential can cancel. The twin bottom and twin leptons are required for anomaly cancellation, but their Yukawa couplings do not need to be equal to the corresponding ones in the SM, as long as they are small enough to not generate a big contribution to the Higgs potential.
The collider phenomenology of the FTH model mainly relies on the mixing of the SM and twin Higgs fields. In typical range of the parameter space, one often expects displaced decays that constitute an interesting experimental signature. Here we focus on the DM. A natural candidate is the twin tau. Since its Yukawa coupling needs not to be related to the SM tau Yukawa coupling, its mass can be treated as a free parameter. It is found that a right amount of thermal relic abundance can be obtained for a twin tau mass in the range of 50-150 GeV if the twin hypercharge is not gauged [3,4]. It corresponds to a twin tau Yukawa coupling much larger than the SM tau Yukawa coupling. The requirement that the twin tau Yukawa coupling does not reintroduce the hierarchy problem puts a upper limit ∼ 200 GeV on the twin tau mass. A twin tau lighter than ∼ 50 GeV would generate an overabundance which overcloses the universe. The relic density can be greatly reduced if a light twin photon also exists, because it provides additional annihilation channels for the twin tau. If the twin photon coupling strength is similar to the SM photon coupling, the annihilation will be too efficient and it will be difficult to obtain enough DM. For a twin photon couplingα ∼ 0.03 α EM , a right amount of relic density can be obtained for a twin tau mass in the range of 1-20 GeV [4].
In the twin tau DM discussion, its stability is assumed to be protected by the twin U(1) EM symmetry. However, if twin U(1) EM is broken and the twin photon has a mass, the twin tau may be unstable and could decay to the twin neutrino if the twin neutrino is lighter. This is because that the twin tau and the twin neutrino can mix due to the twin U(1) EM breaking effect. On the other hand, if the twin lepton number (or parity) is still a good symmetry, the lightest state that carries the twin lepton number (parity) will be stable. In this paper we will assume that it is the twin neutrino and consider its possibility of being the DM.

Twin lepton mixings
If the twin U(1) EM (or equivalently twin hypercharge) is broken, it is possible to write down various Dirac and Majorana masses between the left-handed and right-handed twin tau and twin neutrino fields. For simplicity, we consider the case where the twin lepton number remains a good symmetry, which can be responsible for the stability of DM. This forbids the Majorana mass terms.

JHEP09(2018)098
The twin tau and twin neutrino receive the usual Dirac masses from the twin Higgs VEV, where the subscript B represents the twin sector fields, and H B = iσ 2 H * B transforms as (3, 2) −1/2 under the twin gauge group. The twin hypercharge breaking can be parameterized by a spurion field S which is a singlet under SU(3) B × SU(2) B but carries +1 twin hypercharge (also +1 twin electric charge). It can come from a VEV of a scalar field which breaks U(1) B spontaneously. The radial component is assumed to be heavier than the relevant particles (twin photon, tau, and neutrino) here and plays no role in the following discussion. Using the spurion we can write down the following additional lepton-number conserving mass terms, The mass matrix of the twin tau and twin neutrino is then given by The mass matrix can be diagonalized by the rotations where the mass eigenstates in the twin sector are labelled with a hat (ˆ). We assume that the off-diagonal masses |µ 1 |, |µ 2 | m τ B − m ν B so that the mixing angles θ 1 , θ 2 are small. This is reasonable given that µ 1 , µ 2 arise from higher dimensional operators and require an insertion of the twin hypercharge breaking VEV, which is assumed to be small for a light twin photon. For our analysis, to reduce the number of independent parameters, we further assume that one of the off-diagonal masses dominates, i.e., µ 1 µ 2 , so that we can ignore µ 2 . In this case, we obtain two Dirac mass eigenstatesτ andν which are labeled by their dominant components. The two mass eigenvalues are and the two mixing angles are given by in the small mixing angle limit. There is no qualitative difference in our result if µ 1 and µ 2 are comparable except that the two mixing angles become independent. We are interested in the region of parameter space where the twin tauτ , twin neutrinô ν, and twin photonγ have masses of the same order in the range ∼ 0.1 − 10 GeV, with mτ > mν > mγ. Compared with ref. [7], the twin neutrinoν plays the role of χ which is the DM,τ corresponds to ψ, the coannihilation/coscattering partner of the DM particle, andγ corresponds to the mediator φ. Following ref. [7], we define two dimensionless parameters, which are convenient for our discussion. The region of interest has r < 1 and 0 < ∆ 1.
The mass spectrum and mixing pattern would be more complicated if Majorana masses for the twin leptons are allowed. In addition to the standard Majorana mass for the righthanded twin neutrino, all other possible terms can arise from higher dimensional operators with insertions of the spurion field S (and the twin Higgs field H B ), filling the 4 × 4 mass matrix of (τ B,L , ν B,L , τ c B,R , ν c B,R ). There are four mass eigenstates and many more mixing angles. The stability of the lightest eigenstate can be protected by the twin lepton parity in this case. Because the twin photon couples off-diagonally to Majorana fermions, the analysis of annihilation and scattering needs to include all four fermion eigenstates, which becomes quite complicated. Nevertheless, one can expect that there are regions of parameter space where the correct relic abundance can be obtained through coannihilation and/or coscattering processes in a similar way to the case studied in this work.

Relevant processes for the thermal dark matter abundance
At high temperature, the SM sector and the twin sector stay in thermal equilibrium through the interactions due to the Higgs mixing and the kinetic mixing of the U(1) gauge fields. As the universe expands, the heavy species in the twin sector decouple from the thermal bath and only the light species including twin photon (γ), twin tau (τ ), and twin neutrino (ν) survive. These light species talk to SM mainly via the kinetic mixing term, −( /2)F µνF µν , between γ andγ. We assume that the kinetic mixing is big enough ( 10 −9 ) to keep the twin photon in thermal equilibrium by scattering off light SM leptons [7] during the DM freeze-out. The DM abundance is then controlled by several annihilation and scattering processes.
Annihilation processes: The coupling ofν toγ arises from mixing with the twin tau. As a result, for small mixings the usual annihilation process A for theν DM is suppressed by θ 4 1 , while the coannihilation process C A (where the subscript A stands for asymmetric) is suppressed by θ 2 1 . There is no mixing angle suppression for the coannihilation process C S (where the subscript S represents symmetric or sterile). On the other hand, the Boltzmann suppression goes the other way, the Boltzmann factors for the three processes are ∼ e −2mν /T , e −(mν +mτ )/T , and e −2mτ /T respectively. It is suppressed by θ 2 1 . Becauseτ is assumed to be heavier thanν, the initial states particleŝ ν andγ must carry enough momenta for this process to happen. The coscattering process therefore has a strong momentum dependence. Ignoring the momentum dependence for a moment, the Boltzmann factor can be estimated to be ∼ e −(mτ +mγ )/T .
In addition, there is also the decay procesŝ On-shell decay only occurs if mτ > mν + mγ (∆ > r). In this case the inverse decay (ID) plays a similar role as the coscattering process since both convertν toτ , but the rate is much larger. It turns out that if the (inverse) decay is open, the relic abundance is simply determined by the coannihilation because the inverse decay process decouples later. The majority of the parameter space we focus on has mτ < mν + mγ (∆ < r). In this case, the twin photon has to be off-shell then decays to SM fermions. It is further suppressed by 2 so it can be ignored during the freeze-out. It is however responsible for converting the remainingτ toν eventually after the freeze-out. If the mixing is large and/or ∆ is large so that θ 2 1 e −2mν /T > e −(mτ +mν )/T during freezeout, then the annihilation process A will dominate and we will have the usual WIMP scenario. However, for such a WIMP DM lighter than 10 GeV, this has been ruled out by the CMB constraint (see discussion in section 6). Therefore we focus on the opposite limit θ 2 1 e −2mν /T < e −(mτ +mν )/T , i.e., small mixing and small ∆. In this case we have C S > C A > A in terms of rates. The coscattering process S has the same θ 1 dependence as C A , but is less Boltzmann suppressed because mγ < mν, so the coscattering can keep τ andν in kinetic equilibrium after C A freezes out. In this simple-minded picture, the DM relic density is then determined by freeze-out of C S and S. Denoting their freeze-out temperatures by T C S and T S , then there are two main scenarios.
1. T C S > T S : this occurs if θ 2 1 e −(mτ +mγ )/T e −2mτ /T during freeze-out so that C S freezes out earlier. After that the total number ofν andτ in a comoving volume is fixed. The coscattering process only re-distributes the densities betweenν andτ , but eventually allτ 's will decay down toν's. The DM relic density is determined by C S . This is the coannihilation phase. It is schematically depicted in the left panel of figure 1.
2. T C S < T S : in the opposite limit, S freezes out before C S , and hence stops convertinĝ ν intoτ . On the other hand, C S is still active and will annihilate most of the leftoverτ 's. The relic density in this case is determined by the coscattering process S.
(Remember that A has frozen out earlier.) This is the coscattering phase discovered in ref. [7]. It is illustrated in the right panel of figure 1.
In the above discussion, we have associated each process with a single freeze-out temperature. This is a good approximation for the annihilation and coannihilation processes,

JHEP09(2018)098
Coa nnihila t ion Pha se Left: coannihilation phase where the DM relic density is dominantly determined by C S . Right: coscattering phase discussed in ref. [7] where DM relic density is determined by S.
but not for the coscattering process which has a strong momentum dependence. In the coscattering phase, the processsesνγ →νγ,νν →νν are suppressed by θ 4 1 and hence are expected to freeze out earlier and cannot re-equilibrate theν momentum. Consequently, different momentum modes in the coscattering process freeze out at different time, with low momentum modes freeze out earlier. If θ 2 1 e −(mτ +mγ )/T ∼ e −2mτ /T during freeze-out, we can have a situation that coscattering of the low momentum modes freezes out earlier than C S while the coscattering of the high momentum modes freezes out later than C S . This is illustrated in the left panel of figure 2. In this case, the relic density of low momentum modes is determined by the coscattering process and the relic density of high momentum modes is determined by the coannihilation process C S . We have a mixed coscattering/coannihilation phase where the relic density is determined by both S and C S .
If the masses of the twin photon and the twin neutrino are close, we have θ 2 1 e −(mτ +mγ )/T ∼ θ 2 1 e −(mτ +mν )/T < e −2mτ /T during freeze-out from our assumption. One expects that this belongs to the coscattering phase since T C S < T S . However, the rates of S and C A become comparable in this limit so we have T S ∼ T C A > T C S . Due to the strong momentum dependence of S, one can have a situation depicted in the right panel of figure 2. The coscattering of low momentum modes freezes out early, but their comoving density is still reduced by the coannihilation C A until C A freezes out. The relic density of high momentum modes is determined by S as in the coscattering phase. In this case we have another mixed coscattering/coannihilation phase where the relic density is determined together by C A and S. Finally, if the mixing is not very small so that the rates of C A and C S are not far apart, the freeze-out time of the coscattering process can even cut through that of both C A and C S , although it can only happen in some rare corner of the parameter space. The relic density is then determined by all three processes, with the low Figure 2. Schematic illustrations of the mixed coscattering/coannihilation phases. Left: the S/C S mixed phase: for low (high) momentum modes S freezes out earlier (later) than C S . Right: the C A /S mixed phase: for low (high) momentum modes C A freezes out earlier (later) than S. momentum modes controlled by C A , intermediate momentum modes governed by S, and high momentum modes determined by C S .
Due to the momentum dependence of the coscattering process, the relic density calculations of the coscattering and mixed phases are more complicated. We describe the calculation for each case in the next section.

Relic abundance calculations in various phases
In this section we describe the calculations of DM relic abundance in different phases.

Coannihilation
In the coannihilation phase, the coscattering process S decouples late enough to keepτ andν in chemical equilibrium even after all annihilation processes freeze out, so we have where n (n eq ) is the (equilibrium) number density. We can simply write down the Boltzmann equation for the total DM number density n tot (T ) = nτ (T ) + nν(T ) [9,18]: 2) It can be easily solved and is incorporated in the standard DM relic density calculation packages. In the parameter region that we are interested, the right-handed side is mostly dominated by the C S term.

Coscattering
The DM density calculation in the coscattering phase is more involved and was discussed in detail in ref. [7]. There the authors provided an approximate solution based on the integrated Boltzmann equation: However, as mentioned earlier, the kinetic equilibrium ofν will not be maintained during the freeze-out of the coscattering process and different momentum modes freeze out at different time. The simple estimate from eq. (4.3) is not always accurate. Here we reproduce the calculation from the unintegrated Boltzmann equation to keep track of the momentum dependence.
The unintegrated Boltzmann equation of the density distribution in the momentum space f (p, t) for the coscattering processν(p) +γ(k) →τ (p ) +γ(k ) is given by where the collision operator is defined as [7] C (4.5) In the above expression dΩ p = d 3 p/[(2π) 3 2E p ] is the Lorentz-invariant integration measure and |M| 2 is the squared amplitude averaged over initial and summed over final state quantum numbers.
In this phaseτ andγ can be assumed to be in thermal equilibrium with the thermal bath, fτ (γ) = f eq τ (γ) , from the processesτγ ↔τγ,ττ ↔γγ andγ interactions with SM fields. Using f eq τ (p , t)f eq γ (k , t) = f eq ν (p, t)f eq γ (k, t), the right-hand side of eq. (4.4) can be simplifed as where the reduced collision operator C(p, t) takes the form, with the Lorentz-invariant flux factor j(s) = 2E p 2E k |v p − v k |. The calculation of C from the coscattering is described in appendix A. The left-hand side of eq. (4.4) can be written as a single term by defining the comoving momentum q ≡ p a, then eq. (4.4) becomes a first order differential equation of the scale factor a for each comoving momentum q: Ha∂ a fν(q, a) = [f eq ν (q, a) − fν(q, a)] C(q, a), (4.9) -11 -

JHEP09(2018)098
where we have written the distribution f as a function of q and a, instead of p and t.
Taking the boundary condition fν(q, a 0 ) = f eq ν (q, a 0 ) at an early time a 0 , the solution is given by At early time and high temperature where C(q, a) H(a), the second term on the righthand side can be neglected andν density is given by the equilibrium density as expected. At late time and low temperature, C(q, a) H(a), theν comoving density stops changing. For each comoving momentum q, one can find a time a = a f (q) beyond which the exponent is small so that the exponential factor is approximately 1. Then the finalν density is roughly given by f eq ν (q, a f (q)). The a f (q) can be viewed as the freeze-out scale factor for coscattering of the comoving momentum q ofν. It is roughly determined by C(q, a f (q)) H(a f (q)).
As mentioned in the previous section, if mτ > mν + mγ (r < ∆), the inverse decaŷ ν +γ →τ plays a similar role as the coscattering process. Its contribution should be added to the collision operator, which is calculated in appendix B. It is larger than the coscattering contribution, as it requires less energy to produce the final state. In the parameter region that we consider, it always makes the second term in eq. (4.10) negligible before C S freezes out. Therefore, it goes back to the coannihilation phase if this on-shell decay is allowed.

Mixed phases
From the discussion in the previous section, mixed phases occur when T S ∼ T C S or T S ∼ T C A . We first consider the case T C A ∼ T S > T C S (C A /S mixed phase). This happens when mγ mν. Because the coannihilation process C A ,ν(p) +τ (k) →γ(p ) +γ(k ), is also important in this case, the contribution from C A to the collision operator, should be included in eq. (4.4) in addition to the coscattering contribution of eq. (4.5).
Since T C S is assumed to be small,τ stays in thermal equilibrium during the decoupling of S and C A . In contrast to S, there is no kinematic threshold in C A , so it has very weak momentum dependence. Hence we can treat it as a function of time/temperature only. In eq. (4.10) we can replace C by C S + C C A and conduct the computation in the same manner as in the coscattering phase. In this C A /S mixed phase, C C A > C S for low momentum modes so their freeze out temperature is determined by C A , while for high momentum modes C S > C C A and their contributions is given by the coscattering result, as illustrated in the right panel of figure 2 For the other mixed phase (S/C S ) where T C S ∼ T S < T C A , the calculation is more complicated. In this caseτ is no longer in thermal equilibrium with the thermal bath, fτ itself is unknown, hence cannot be set to equal f eq τ . As a result, eq. (4.6) no longer holds. Moreover, different from the usual coannihilation scenario,τ andν are not in chemical equilibrium. A complete solution requires solving the two coupled Boltzmann equations forτ andν in this case, which is numerically expensive. However, we can assume that -12 -

JHEP09(2018)098
τ is still in kinetic equilibrium with the SM sector (i.e., has the canonical distribution up to an unknown overall factor) due to the elastic scattering withγ. We can write fτ /f eq τ = Yτ /Y eq τ where Yτ = nτ /s is the comoving number density. The term in the square bracket of eq. (4.5) can be written as: This corresponds to replacing f eq ν (q, a) by (Yτ (a)/Y eq τ (a))f eq ν (q, a) in eq. (4.9). The solution eq. (4.10) is modified to [8]  Of course we do not know Yτ (a) in advance. One way to solve this problem is to employ the iterative method [8] which is described in appendix C. One starts with some initial guess of Yτ (a) to obtain the solution for eq. (4.13), then use that solution in the Boltzmann equation for nτ to obtain a new Yτ (a), and repeat the procedure until the result converges. We find that a good first guess is to simply take Yτ (a) to be the one obtained in the coannihilation calculation. In coannihilation,τ andν are in chemical equilibrium: where the superscript CA indicates that the result is obtained from the coanihillation-only estimation (eq. (4.2)). Notice that in this case the first term in eq. (4.13) is simply the contribution from coannihilation. If T S < T C S , S will still be active when C S freezes out, C(q, a )/(Ha ) will be large and the second term in eq. (4.13) will be suppressed. We obtain the correct coannihilation limit. On the other hand, if T S > T C S , the coannihilation C S is effective to keep Yτ (a)/Y eq τ (a) ≈ 1 and eq. (4.13) returns to the coscattering result in eq. (4.10). Using Yτ (a) from the coannihalation calculation in eq. (4.13) gives the correct results in both the coannihilation and coscattering limits. The expression interpolates between these two limits in the mixed phase and it turns out to be an excellent approximation to the correct relic density even without performing the iterations. (See appendix C.) For completeness, we can include the C A contribution in eq. (4.13), then this result also applies to the case if T S cuts through both T C A and T C S . Figure 3 shows the differential DM density as a function of the comoving momentum in different phases, all calculated from eq. (4.13) including all contributions. We will use this formula for numerical calculations of DM densities in the next section in all phases.

Numerical results
In this section we present the numerical calculations of the dark matter relic abundance for various parameter choices of the model. We fixê = 0.3 for the twin electromagnetic gauge  Figure 3. Benchmark results for relic density calculation in the coscattering (left), mixed (central), and coannihilation (right) phases. The red curves are the contributions from the first term of eq. (4.13) which corresponds to a pure coannihilation calculation. The blue solid curves are due to the coscattering contributions from the second term of eq. (4.13). The purple curves are the total contributions. We can see that in the coscattering phase and the coannihilation phase the total contribution is dominated by one term, while in the mixed phase both terms give comparable contributions. The blue dashed curves are calculated from the coscattering formula of eq. (4.4). Only in the coscattering phase the blue dashed curve approximates the correct result.

JHEP09(2018)098
coupling and calculate the DM density in units of the observed Ω obs h 2 = 0.12. All calculations are performed using eq. (4.13) including contributions from all relevant processes. In figure 4, we consider mν = 1 GeV and plot the DM density dependence on r = mγ/mν and the mixing angle θ 1 , for ∆ = (mτ − mν)/mν = 0.1 and 0.18. The contours are in log 10 (Ω DM /Ω obs ) and the "0" contour represents points which produce the observed DM density. The coscattering phase sits in the lower-right region, as for small θ 1 and large r the coscattering process S is suppressed and freezes out earlier. The upper-left region, on the other hand, belongs to the coannihilation phase. The green band separating them corresponds to the S/C S mixed phase. The boundaries of the green band are determined by the condition T S (q = p a = 0) = T C S and T S (q = 25) = T C S . The contribution to the relic density from modes with q > 25 is small and is ignored in the coscattering calculation. The orange vertical dashed line corresponds to T S (q = 0) = T C A . To the right of it C A becomes relevant and we enter the C A /S mixed phase. The DM relic density is slightly reduced by C A compared to a pure coscattering calculation. The red vertical dashed line indicates r = ∆. To the left of it the on-shell decay and inverse decayτ ↔νγ are open so this whole region is in the coannihilation phase.
The relic density in the coannhilation phase is mostly independent of θ 1 because it is mainly controlled by C S which hardly depends on θ 1 . Only at larger θ 1 values when C A and A become relevant the DM relic density shows some θ 1 dependence. The dependence on r of the relic density in the coannihilation phase is also mild, as it mainly affects the phase space of the coannihilation process. In the coscattering phase, the DM relic density increases as ∆ increases, which can be seen by comparing the two plots in figure 4. This is because a larger gap between mτ and mν requires a higher threshold momentum forγ to make S happen. Therefore the coscattering is suppressed, resulting in a larger relic density. For mν = 1 GeV, ∆ = 0.1, the observed DM density is produced in the coscattering phase, while for ∆ = 0.18 it moves to the mixed phase or coannihilation phase.
In figures 5 and 6 we show the results for mν = 100 MeV and 10 GeV. A larger DM mass will give a larger relic density if all other parameters are fixed. Consequently for a correct relic density we need a smaller (larger) ∆ for a larger (smaller) mν. The ∆ values are chosen to be 0.2 and 0.26 for the two plots with mν = 100 MeV, and 0.04 and 0.1 for the two plots with mν = 10 GeV. The behaviors of the contours are similar to the case of mν = 1 GeV.

Experimental constraints and tests
In this section we discuss the current experimental constraints and future experimental tests of this model. A comprehensive summary for this type of DM scenarios can be found in ref. [19].

Direct detection
The dark matterν interacts with SM particles through the Higgs portal or dark photon portal. The Higgs portal is suppressed by the Higgs mixing between the SM and twin sectors, and the twin neutrino Yukawa coupling. The dark photon portal is suppressed by the kinetic mixing , and also the mixing angle θ 1 . Most current DM direct detection experiments are based on heavy nuclei recoiling when the nuclei scatter with the DM particles. They become ineffective for light DM less than a few GeV, but could potentially constrain the parameter space of heavier twin neutrino region, as the Yukawa coupling also becomes larger for a heavier twin neutrino. For mν=10 GeV and f /v = 3, theν-nucleon cross section is of O(10 −47 ) cm 2 , which is dominated by the Higgs portal. Suchν DM is not yet constrained by recent liquid xenon DM detectors [20][21][22]. A conservative estimate according to ref. [23] gives an upper limit of ∼22(60) GeV on mν for f /v = 3(5). For even lighter DM, recent upgrades/proposals of detecting DM-electron scattering can largely increase the sensitivity for sub-GeV mass DM [24,25]. In our model, theν-e where the reference values of mixing parameters , θ 1 have been chosen close to the upper bounds to maximize the cross section. This is not yet constrained by recent electronscattering experiments, including SENSEI [26], Xenon10 [27], DarkSide-50 [28]. Future upgrades will be able to probe part of the parameter space with large mixings. It is also worth mentioning that there are also crystal experiments based on phonon signals coming from DM scattering off nuclei in the detector, such as CRESST-III [29]. Thanks to the low energy threshold (O(50 eV)), these experiments will also be sensitive to the sub-GeV DM mass region. For such low DM mass, the dark photon portal becomes important and can dominate over the Higgs portal interaction if the mixings are not too small. However, the current constraint still can not put any bounds on mν even for the f /v=3 case. Significant progress in the future could be helpful to constrain the parameter space for a lightν.

Indirect constraints induced from DM annihilation
Light DM is in general strongly constrained by indirect searches due to its high number density. WIMP models with annihilation cross section σv 10 −26 cm 3 /s and m DM 10 GeV have already been ruled out [30]. In our model the DM relic density is not determined by the DM annihilation process, but by the coannihilation and coscattering processes. The DM annihilation is dominated byνν →γγ → 4f , which is suppressed byê 4 θ 4 1 . An upper bound on the DM annihilation cross section gives a constraint on the combination of the parametersêθ 1 . Fermi-LAT data [31] has put an upper limit on the DM annihilation cross section for DM heavier than 6 GeV. A stronger constraint comes from CMB observables, which restrict the net energy deposited from DM annihilation into visible particles during the reionization era [32]. The constraints from the Fermi-LAT and the Planck data are plotted in figure 7. Note that the annihilation ofνν can produce 4e instead of 2e. This may modify the bounds derived from the 2e final state. The total energy injection is the same, and the 4e final state will result in more electrons but with lower energies. Ref. [33] performed a detailed study in comparing the constrains for cascade decays with different numbers of final sate particles. After convoluting with the energy dependence of the efficiency factor f eff [34], it is found that the effects due to multi-step decay is rather mild. The constraints on the 4e and 2e final states from the Planck data are roughly the same. On the other hand, a higher-step cascade tends to soften the spectrum and thus slightly weakens the constraint from the Fermi-LAT result [33]. The proposed ground-based CMB Stage-4 experiment [35] is expected to improve the constraint by a factor of 2 to 3 compared to Planck. The DM annihilation after the CMB era will also heat the intergalactic hydrogen gas and erase the absorption features of 21cm spectrum around z 17. The recent measurement by the EDGES experiment instead observed an even stronger absorption than the standard astrophysical expectation [36]. If one interprets this result as a constraint that the DM annihilation should not significantly reduce the absorption, the observed brightness temperature then suggests an even stronger bound than the one from CMB [37][38][39][40]. In figure 7 we also plot the most conservative constraint in terms ofê × θ 1 according to ref. [37], taking the efficiency factor to be 1. For mν = 1 GeV, the upper limit forêθ 1 is between 10 −2 and 10 −3 depending on other model parameters. 2 The constraint gets more stringent for lighter mν. 2 Theνν annihilation cross section depends on the model parameters as For a fixed mixing angle, the annihilation rate reaches its maximum around ∆ ∼ 0.5 for very small r, while for larger r, it decreases as ∆ increases.

Constraints induced by the light twin photon
In our scenarioγ is the lightest twin sector particle. An on-shellγ decays through its kinetic mixing with the SM photon and there is no invisible decay mode to the twin sector. In this case, it is well described by two parameters: the twin photon mass mγ = rmν and the kinematic mixing with the SM photon . Experimental constraints on dark photon have been extensively studied. Summeries of current status can be found in refs. [25,[41][42][43].
A lower bound on mγ comes from the effective number of neutrinos (N eff ) [44]. A lightγ can stay in thermal equilibrium with photons and electrons after the neutrinos decouple at T ∼ 2.3 MeV. The entropy transferred fromγ to the photon bath will change the neutrino-photon temperature ratio, and therefore modify N eff . Using the results in ref. [44] and the Planck data [30], we obtain a lower bound on mγ around 11 MeV. It may be further improved to ∼ 19 MeV by the future CMB-S4 experiment [35].
There are also many constraints on depending on mγ. The current upper bounds for mostly come from colliders, fixed target experiments and meson decay experiments, in searching for prompt decay products. (See ref. [43] for a summary and an extended reference list.) These experiments constrain 10 −3 in the mass range that we consider (except for a few narrow gaps at the meson resonances). The lower bounds for come from γ displaced decays from various beam dump experiments [45][46][47][48][49][50][51][52][53][54][55][56][57][58], and also from supernova SN1987A [59]. The Big Bang Nucleosynthesis (BBN) would also constrain the lifetime of γ. However, the decay ofγ is only suppressed by and the BBN does not introduce extra constraints for 10 −10 [59], which is required in this model to keep the DM sector in thermal contact with the SM. These bounds are summarized in figure 8.

Constraints induced fromτ decay
As we require mν > mγ, the twin photon only decays to SM particles, thusτ can only be pair produced in a lab via an off-shell twin photon or from the Higgs boson decay. The constraints fromτ pair production from off-shellγ are weaker compared to the ones fromγ visible decay modes described in the previous subsection. Moreover,ν pair production viâ γ * will be further suppressed by θ 4 1 , leaving h →νν/ττ to be the main production channel. At the LHC, theτ produced from h/γ * will be long-lived in general, if the two-body decayτ →νγ is forbidden (r > ∆), because the leading three-body decay is suppressed by θ 2 1 2ê2 e 2 . Assumingτ andν are Dirac fermions and taking mτ mν, the current upper bound of the Higgs invisible decay branching ratio, Br(h → invisible) < 24% [60,61], constrains mν to be 19 (52) GeV for f /v = 3 (5). HL-LHC is expected to improve the Higgs invisible branching ratio measurement to 6-8% [62], which would translate to a bound 11 (30) GeV for mν. Future e + e − colliders can probe BR(h → invisible) to the sub-percent level [63][64][65][66]. A 0.3% measurement can constrain mν down to 2 (6) GeV for f /v = 3 (5).
Theτ decay width can be expressed analytically in the small ∆ limit (1 r ∆): Γτ →νe + e − θ 2 1 ∆ 5 e 2ê2 2 60π 3 r 4 mν,   Figure 8. Constraints on the kinetic mixing parameter and the twin photon mass mγ. The green and cyan shaded regions are ruled out by lab experiments. The magenta shaded region is the constraint from SN1897A cooling. In the red shaded region is too small to keepγ in thermal equilibrium with the SM. We also plot 6 benchmark models which give the correct DM relic density from the numerical calculations in section 5. For small enough , models 1, 3, 4, and 6 are in the coscattering phase and model 2 is in the mixed phase. At large values of these model curves turn right because the three-body (inverse) decay rate becomes large and freezes out after the coannihilation process, driving the models into the coannihilation phase. Model 5 is in the coannihilation phase for all large enough to keepγ in thermal equilibrium. The ticks on each benchmark curve representτ lifetime, starting from τ (τ )=1 sec and increasing by 10 2 each tick below. The dashed parts of the curves are ruled out by the BBN constraint.
which strongly depends on ∆ and r. Numerically the proper decay length is given by where N f is the number of SM fermions that can appear in the final state. The dark sector could be probed by searching forτ displaced decays at the HL-LHC for cτ (τ ) ∼ O(1) m [67].

JHEP09(2018)098
only makes up a small fraction of Ω after it freezes out (Ωτ h 2 10 −3 ), also the fact that only a small fraction mτ −mν mτ = ∆ 1+∆ of energy would be injected. The strongest bound comes from the electromagnetic decay products and depends on the model parameters. Typically, the lifetime of the long-livedτ is only constrained by BBN to be shorter than ∼ 10 6 seconds [71,72]. The extra energy injection fromτ decays could distort the CMB blackbody spectrum which can potentially be captured by the proposed PIXIE mission [73]. It may give a slightly stronger bound than BBN on τ (τ ) [71], also around 10 6 sec for our typical benchmark models.
In figure 8, we also plot several benchmark points with fixed r, ∆, and θ 1 , that give rise to the correct DM relic density from our numerical results in section 5. The vertical parts indicates that the thermal relic density is mostly independent of the gauge kinematic mixing parameter , as long as it can keep the DM in thermal equilibrium with SM before freeze-out. For large values of , the three-body (inverse) decay rate gets larger, and can even freeze out later than the coannihilation process C S . It drives the benchmark models to the coannihilation phase even if it was in the coscattering or mixed phase for smaller values of . This occurs for all our benchmark models except for model 5 (red line) which is in the coannihilation phase for all . The transition to the coannihilation phase reduces the DM number density, hence it needs a larger DM mass to compensate the effect. This explains the turn of the curves at lager values. However, except for model 2 (orange line), the turns occur in the region which has been ruled out by other experiments. For smaller , theτ lifetime becomes longer, which is constrained by the BBN bound. The region that violates the BBN constraint is indicated by dashed lines.

Conclusions
The necessity of DM in the universe is one of the strongest evidences of new physics beyond the SM. Experimental searches in various fronts so far have not revealed the nature of the DM. For the most popular WIMP DM scenario, recent advancements in experiments have covered significant fractions of the allowed parameter space, even though there are still viable parameter space left. People have taken more seriously the possibility that DM resides in a more hidden sector, and hence has escaped our intensive experimental searches. Even in this case, it would be more satisfactory if it is part of a bigger story, rather than just arises in an isolated sector for no particular reason. In this paper, we consider DM coming from a particle in the twin sector of the fraternal twin Higgs model, which itself is motivated by the naturalness problem of the SM EW symmetry breaking and non-discovery of the colored top partners at the colliders. Although the relevant particles for the DM relic density in our study, i.e., the twin neutrino, twin tau, and twin photon, have little effect on the naturalness of the EW scale, they are an integral part of the full theory that solves the naturalness problem, just like the neutralinos in a supersymmetric SM.
To obtain the correct DM relic density, the interplay of the twin neutrino, twin tau, and twin photon is important. The DM relic density is determined by the order of the freeze-out temperatures of various annihilation and scattering processes. It is in the coannihilation phase if the twin tau annihilation freezes out earlier than the twin neutrino to -21 -JHEP09(2018)098 twin tau scattering. In the opposite limit, it realizes the recently discovered coscattering phase. There is also an intermediate regime where the DM relic density is determined by both coannihilation and coscattering processes due to the momentum dependence of the coscattering process. The calculation of the DM relic density in this mixed phase is more complicated and has not been done in the literature and we provide a reasonably simple way to evaluate it with very good accuracies.
There are many experimental constraints but none of them can cover the whole parameter space. Direct detection with nuclei recoiling can only constrain heavier DM, above a few GeV. The experiments based on electron scattering are not yet sensitive to this model. Future upgrades may be able to probe the region of the parameter space with large mixings between the twin tau and the twin neutrino. Indirect constraints from DM annihilation are more sensitive to smaller DM mass with large enough twin tau -twin neutrino mixings. Other constraints rely on the coannihilation/coscattering partners of the DM. Twin photon is subject to various dark photon constraints. Twin tau typically has a long lifetime so it is constrained by BBN and CMB. Its displaced decays may be searched at colliders with dedicated detectors or strategies. A big chunk of the parameter space still survives all the constraints. A complete coverage of the parameter space directly is not easy in the foreseeable future. An indirect test may come from the test of the whole fraternal twin Higgs model at a future high energy collider, if other heavier particles in the model can be produced.

A Calculation of the collision operator
The evaluation of the reduced collision operator C can be time consuming to achieve a high precision. It can be simplified by performing a partial analytic integration of the phase space, leaving a one-dimensional integral for the numerical calculation. A similar treatment for a case with a massless initial state was done in ref. [8]. Notice that since dΩ k and the final state integrated |M| 2 in eq. (4.8) are Lorentz invariant, C(p, t) can be evaluated in the rest frame ofν. The initial stateγ momentum k r in theν rest frame is a function of both p and k, obtained from a simple boost. In this frame, the density distribution ofγ is no longer spherically symmetric, but is still axially symmetric along the boost axis. In addition, the Mandelstam variable s is a function of | k r | only, so the integration over the angular variables can be easily performed and the result is When the decay D and the inverse decay ID [ν(p) +γ(k) ↔τ (p )] are kinematically allowed, the effect can also be included in the collision operator by an extra term C ID . The cross section of the inverse decay takes the form: where the flux j(s) = 4EνEγv. The unintegrated Boltzmann equation (4.6) now reads with the inverse decay contribution given by where kmax min = 1 2 p ∆(∆ + 2) − r 2 ± m 2 ν + p 2 (∆ 2 (∆ + 2) 2 + r 4 − 2 (∆ 2 + 2∆ + 2) r 2 ) .
(B.6) Notice that the kinematic threshold of the inverse decay is lower than that of the coscattering because it does not need to create a twin photon in the final state. For example, the ratio of theγ threshold momenta for p = 0 is k T (S) k T (ID) = ∆(∆ + 2)(∆ + 2r)(∆ + 2r + 2) |r 2 − ∆(∆ + 2)| , (B.7) which is larger than one for any ∆ and r. The interaction rate of S is hence exponentially suppressed compared to the rate of ID, besides that it has an extraê 2 suppression. Consequently ID decouples much later than S. To compare the decoupling time between the inverse decay and the coannihilation C S , we can simply compare C ID (p = 0, t) with H(t) when C S decouples. For p = 0, C ID is simplified to: We find that in the parameter space that we consider with θ 1 ≥ 10 −6 , the (inverse) decay is still active when C S decouples. As an example, for θ 1 = 10 −6 ,ê = 0.3, r = 0.24, ∆ = 0.26, mν = 0.1 GeV, ID(p = 0) decouples at x 48, much later than C S that decouples at x 30. The difference is even larger for lower values of ∆ and r. From these comparisons, we conclude that when r < ∆, the (inverse) decay will keepν andτ in chemical equilibrium and makes the relic density follow the coannihilation result.

JHEP09(2018)098
C The relic density calculation from iteration As the coscattering process S freezes out, the assumption of chemical equilibrium between τ andν (eq. (4.1)) breaks down. As a result, the combined Boltzmann equation forτ and ν eq. (4.2) no longer holds after S freezes out. In the mixed phase when T C S and T S are comparable, in principle one should solve the coupled Boltzmann equations for both Yν and Yτ . As the distribution ofτ remains canonical due toτγ scattering, we can integrate out the momentum dependence to obtain the integrated Boltzmann equation for nτ , nτ + 3Hnτ = − σv C S (n 2 τ − (n eq τ ) 2 ) − σv C A (nτ nν − n eq τ n eq ν ) − σv IS n eq γ nτ − n eq τ nν n eq ν , (C.1) where the last term is the contribution from the inverse coscattering processτγ →νγ, which does not have a threshold and has very weak pτ dependence. The combination of eq. (C.1) and eq. (4.13) will give the full solution of nν and nτ .
One way to solve the coupled equations is to use the iterative method [8]. With an initial guess of nν/n eq ν , we can obtain nτ /n eq τ from solving eq. (C.1). Plugging it into eq. (4.13) will return an improved value for nν/n eq ν . Repeating this process will eventually converges to the exact result.
We take the number densities obtained from the coannihilation calculation as our starting point. We argued that in this case eq. (4.13) gives the correct results both in the coannihilation limit and the coscattering limit without iterations. The only possibility of significant deviation is when their contributions are comparable, i.e., in the mixed phase. To examine the corrections from iterations, we perform numerical studies for the benchmark parameters, mν = 1 GeV, ∆ = 0.1, r = 0.5 and y = 0.3, with θ 1 = 10 −2 , 10 −3 , 10 −5 , which correspond to coannihilation, mixed, and coscattering phases respectively. We found that the correction of nν − n eq ν from a single iteration is always small. The results are shown in figure 9. The first plot represents a coannihilation-like phase. The correction due to the coscattering contribution from the first iteration is less than 0.1%. For smaller θ 1 the coscattering becomes more important one can see larger deviations in Yτ , due to the larger Yν from decoupling of S. In the coscattering phase as shown in the bottom plot, the correction of Yτ at late time can be significant. However, during S freeze-out (which occurs at x 10 for θ 1 = 10 −5 ),τ is still approximately in equilibrium, and the relic density is mostly given by nν. The correction to the total DM density is also small (O(10 −6 ) in this case). The largest correction from the iteration indeed happens in the mixed phase, although it is still quite small. For θ 1 10 −3 , the correction to Ωh 2 is 0.4% from the first iteration. The reason for such small corrections is that Yτ is not as sensitive to the early S decoupling as Yν is, therefore using Yτ from the coannihilation calculation in eq. (4.13) gives a very good approximation. Based on these results, in our numerical calculations we simply adopt this prescription without iterations.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.
-24 -JHEP09(2018)098 Figure 9. The ratios of Yτ from the first iteration and the initial input from the coannihilation calculation for 3 benchmarks differ by the mixing angle θ 1 described in the text. The top plot corresponds to a coannihilation phase, whereν andτ are in approximate chemical equilibrium during the freeze-out. The bottom plot corresponds to a coscattering phase, whereτ only deviates from the initial guess after S has already decoupled, thus it is unable to affect S freeze-out. The middle plot represents a mixed phase. In this case, the maximal deviation happens simultaneously as S decouples, resulting in a larger iterative correction to the DM density. However, the correction is still very small.