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 $\sim$ 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.


I. 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 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. 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].
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. [8]. 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 [9], DarkSUSY [10], and MadDM [11] 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 was not discussed in the original coscattering paper [7].
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 Sec. II we first give a brief summary of the fraternal twin Higgs model and its possible DM candidates. Then we focus the discussion 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 Sec. III we enumerate the relevant processes and discuss their roles in controlling the DM abundance in different scenarios. Sec. IV 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 Sec. V. In Sec. VI we discuss various experimental constraints and future probes of this DM scenario. The conclusions are drawn in Sec. VII. 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.

II. FRATERNAL TWIN HIGGS AND LIGHT DM
• The twin SU (2) and SU (3) [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. 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 where 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 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 right- 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.
(Co)scattering process:νγ 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 Coannihilation Phase 2. T C S < T S : In the opposite limit, S freezes out before C S , and hence stops convertingν 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 Fig. 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, 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

DM stops annililation
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. 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 Fig. 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 Fig. 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 momentum modes controlled by C A , intermediate momentum modes governed by S, and high momentum modes determined by 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.

IV. RELIC ABUNDANCE CALCULATIONS IN VARIOUS PHASES
In this section we describe the calculations of DM relic abundance in different phases.

A. 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 ) [8,17]: 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.

B. 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: nν + 3Hnν = − σv S (nν − n eq ν )n eq γ .
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. (11) 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 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. (12) can be simplifed 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. (12) can be written as a single term by defining the comoving momentum q ≡ p a, then Eq. (12) 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), 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. (18) negligible before C S freezes out. Therefore, it goes back to the coannihilation phase if this on-shell decay is allowed.

C. 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. (12) in addition to the coscattering contribution of Eq. (13). 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. (18) 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 Fig. 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. (14) 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τ 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. (13) can be written as: fτ (p , t)fγ(k , t) − fν(p, t)fγ(k, t) = f eq τ (p , t) This corresponds to replacing f eq ν (q, a) by (Yτ (a)/Y eq τ (a))f eq ν (q, a) in Eq. (17). The solution Eq. (18) is modified to Of course we do not know Yτ (a) in advance. One way to solve this problem is to employ the iterative method [18] which is described in Appendix C. One starts with some initial guess of Yτ (a) to obtain the solution for Eq. (21), 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. (10)). Notice that in this case the first term in Eq. (21) 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. (21) 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. For completeness, we can include the C A contribution in Eq. (21), then this result also applies to the case if T S cuts through both T C A and T C S . Fig. 3 shows the differential DM density as a function of the comoving momentum in different phases, all calculated from Eq. (21) including all contributions. We will use this formula for numerical calculations of DM densities in the next section in all phases. 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 Fig. 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.

VI. 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].
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 [25], Xenon10 [26] and DarkSide-50 [27]. Future upgrades will be able to probe part of the parameter space with large mixings.

B. 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 [28]. 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 [29] 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 [30]. The constraints from the Fermi-LAT and the Planck data are plotted in Fig. 7. The proposed ground-based CMB Stage-4 experiment [31] 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 [32]. 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 [33][34][35][36]. In Fig. 7 we also plot the most conservative constraint in terms ofê × θ 1 according to Ref. [33], 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ν.

C. 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. [24,[37][38][39].
A lower bound on mγ comes from the effective number of neutrinos (N eff ) [40]. 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 neutrinophoton temperature ratio, and therefore modify N eff . Using the results in Ref. [40] and the Planck data [28], we obtain a lower bound on mγ around 11 MeV. It may be further improved to ∼ 19 MeV by the future CMB-S4 experiment [31].

D. 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γ 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,  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 .
Theτ decay width can be expressed analytically in the small ∆ limit (1 r ∆): 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 [63].
If τ (τ ) O(1) second, the decays ofτ thermal relic will inject energy during the Big Bang Nucleosynthesis (BBN) era and the recombination era. However, the constraint is weakened by the fact thatτ 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 [67,68]. The extra energy injection fromτ decays could distort the CMB blackbody spectrum which can potentially be captured by the proposed PIXIE mission [69]. It may give a slightly stronger bound than BBN on τ (τ ) [67], also around 10 6 sec for our typical benchmark models.
In Fig 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 Sec. V. 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. 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 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. final state integrated |M| 2 in Eq. (16) 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 where k rt = 1 2mν (m 2 τ − m 2 ν )[(mτ + 2mγ) 2 − m 2 ν ] is the threshold momentum ofγ for the upward scattering.
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)| , 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: C ID (0, t) = 1 4π 2 π 2m 2 ν |M ID | 2 kf eq γ (k, t).
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.
FIG . 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.
One way to solve the coupled equations is to use the iterative method [18]. With an initial guess of nν/n eq ν , we can obtain nτ /n eq τ from solving Eq. (C1). Plugging it into Eq. (21) 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. (21) 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 Fig. 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. (21) gives a very good approximation. Based on these results, in our numerical calculations we simply adopt this prescription without iterations.