{'fraction_non_alphanumeric': 0.03197743789998915, 'fraction_numerical': 0.014101312506779477, 'mean_word_length': 4.175836514709185, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Charged Dark Matter in Supersymmetric Twin Higgs models

Supersymmetric Twin Higgs models ameliorate the fine-tuning of the electroweak scale originating from the heavy scalar top partners required by the non-discovery of them at the Large Hadron Collider. If the Lightest Supersymmetric Particle resides in the twin sector, it may play the role of dark matter even if it is charged under twin gauge interactions. We show that the twin stau is a viable candidate for charged dark matter, even if the twin electromagnetic gauge symmetry is unbroken, with thermal relic abundance that naturally matches the observed dark matter abundance. A wide parameter space satisfies all the experimental constraints including those on dark matter self-interactions. Twin stau dark matter can be observed in future direct detection experiments such as LUX-ZEPLIN. The stau has a mass in the range of 300-500 GeV, and in the minimal scenario, has a decay length long enough to be observed as a disappearing track or a long-lived particle at the Large Hadron Collider.


Introduction
The electroweak (EW) scale is unnaturally small in comparison with other energy scales of the fundamental law such as the Planck scale. Such a hierarchy may be explained by supersymmetry (SUSY) with the masses of the superpartners of the Standard Model (SM) particles around the EW scale [1][2][3][4]. No such particles have been found so far, and their masses should be above the TeV scale [5,6]. Because of the large masses of superpartners, in the minimal embedding of the SM into a supersymmetric theory called the minimal supersymmetric SM (MSSM), the EW scale is obtained by fine-tuning of the parameters of the theory.
The Twin Higgs (TH) mechanism [7] can relax the required degree of fine-tuning in supersymmetric theories [8][9][10][11][12][13] and provide a better understanding on the origin of the EW scale. The mechanism is based on a Z 2 symmetry and predicts twin partners of the SM particles as well as those of the superpartners of the SM particles. As in the MSSM, the Lightest Supersymmetric Particle (LSP) is a dark matter (DM) candidate. The LSP may be the twin superpartner of a SM particle. Ref. [14] investigates the possibility that the twin bino-like LSP is DM and studies its phenomenology such as DM direct detection.
In this paper, we point out that the LSP may be a twin charged particle, such as a twin stau, a twin stop, and a twin gluino. Although they are charged, they do not have a long-range force nor strong interactions with SM particles, and are good DM candidates. If the LSP is twin-colored, assuming that the abundance is determined by the freeze-out of its annihilation, its mass is required to be above a few TeV. 1 Since the Higgsino mass must be also above a few TeV in this case, the supersymmetric Higgs mass parameter is large and the small EW scale requires fine-tuning. We thus focus on the twin stau LSP, whose mass may be just a few 100 GeV for a successful freeze-out scenario.
DM in the framework of Twin Higgs has been intensively studied in recent years. The role of DM may be played by twin neutrons, twin neutral atoms, or twin electrons [15][16][17][18][19][20][21][22]. Many more DM candidates have been proposed in various versions of Fraternal Twin Higgs model [23], e.g., twin taus, neutrinos, mesons, or bottom baryons [15,16,[24][25][26][27][28][29][30][31][32]. In models of twin charged DM proposed so far, such as twin electrons or twin taus, the twin electromagnetic gauge symmetry is either broken or not present at all; otherwise long-range interactions mediated by twin electromagnetism would exclude these DM candidates. A novel feature of twin stau DM is its long-range self-interaction that is compatible with astrophysical data. This is possible due to soft supersymmetry breaking contributions to its mass O(100) GeV.
This paper is organized as follows. In Sec. 2, we review supersymmetric Twin Higgs models and explain how the twin stau can be the LSP. In Sec. 3, experimental and observational constraints on the twin stau LSP are formulated. The lifetime of the MSSM stau, which is important for collider constraints and the computation of the thermal twin stau abundance, is computed in Sec. 4. In Sec. 5, we put together the constraints and show the viable parameter space and future prospects. The compatibility of the scenario with various solutions to the dark radiation problem in the Twin Higgs mechanism is discussed in Sec. 6. In some case, the LSP mass may be correlated with the abundance of dark radiation. Finally, Sec. 7 is devoted to the conclusions.

Supersymmetric Twin Higgs models
In this section, we briefly review SUSY Twin Higgs models and discuss the possibility of the twin stau LSP.

Higgs potential and tuning
Two ingredients are needed in successful SUSY TH models with a natural EW scale: • The potential of the SM and twin Higgses must be approximately SU(4) symmetric, so that the SM Higgs becomes a pseudo Nambu-Goldstone boson associated with the spontaneous breaking of the SU(4) symmetry by the twin Higgs, and as a result, the correction to the SM Higgs mass parameter is suppressed. The quadratic terms are automatically SU(4) symmetric because of the Z 2 symmetry. The minimal quartic terms given by the D term of the EW and twin EW gauge interactions, on the other hand, are not SU(4) symmetric. Therefore, an additional large SU(4) invariant quartic term is necessary.
• The SM Higgs mass of 125 GeV must be obtained without large stop masses that lead to large tuning of the Higgs mass parameter.
In the original SUSY Twin Higgs model [8,9], an F-term of a singlet chiral field S generates an SU(4) invariant quartic term via a superpotential W ∼ S(H u H d + H u H d ). The EW D-term generates an SU(4)-breaking quartic coupling which contributes to the Higgs mass. These quartic terms depend on the ratio of the vacuum expectation value (VEVs) of the up-type Higgs to that of the down-type Higgs, tan β (assumed to be the same for the twin sector). If one takes a large tan β to increase the tree-level Higgs mass from the D-term, so that the required stop mass is small, the SU (4)-invariant F-term becomes small. In this case, the Twin Higgs mechanism does not work well. On the other hand, taking small tan β allows for a large SU (4) symmetric quartic term, at the price of a small tree-level Higgs mass. This requires a large stop mass. As a result, the EW scale is fine-tuned at least at the level of 1% [10]. An alternative was proposed in [11], leading to fine-tuning better than a few percent and, in some of the parameter space, as good as ten percents. The MSSM and its twin copy are supplemented by an abelian U(1) X gauge symmetry generating an SU (4)-invariant D term. The U(1) X is spontaneously broken, but if the masses of the symmetry breaking fields are dominated by the SUSY breaking soft mass, the D-term potential does not decouple. In this model the quartic coupling is maximised in the limit of large tan β, and the Higgs mass of 125 GeV is easily obtained without large stop masses. In order to obtain a large SU(4) invariant quartic, the U(1) X gauge coupling is required to be large and consequently its Landau pole appears below the Planck scale. However, the U(1) X gauge symmetry may be embedded into a non-abelian group, so that the Landau pole is above the Planck scale or the theory is asymptotically free [12,13]. The model is then fully perturbative below the Planck scale.

The lighest supersymmetric particle
If R-parity is conserved, the LSP is absolutely stable. In principle, the LSP may be a superpartner of a twin particle rather than that of a SM particle. The twin bino LSP has been already investigated in [14]. In this paper, we consider the twin stau LSP.
Assuming that soft SUSY breaking terms in the MSSM and the twin sector are the same due to the Z 2 symmetry, the physical masses of the MSSM and the corresponding twin states differ only due to the larger vacuum expectation value (VEV) of the twin Higgs, v , compared to that of the SM Higgs, v. For tan β > 1, both the right-handed and left-handed (twin) stau obtain a positive mass squared from the D-term coupling with the (twin) Higgs. As a consequence, the pure right-handed and left-handed twin staus are heavier than their MSSM counterparts. Moreover, for the pure left-handed twin stau the twin sneutrino is even lighter and play the role of the LSP.
Nevertheless, the twin stau can be the LSP because of the mixing between the lefthanded and right-handed states. The off-diagonal entry of the squared mass matrix of MSSM (twin) staus is proportional to m τ (τ ) (A τ − µ tan β), where µ is the Higgsino mass parameter and A τ the soft SUSY breaking tau trilinear coupling, taken to be the same for the MSSM and twin sector. Hence, the mixing effects are bigger in the twin stau sector, favouring a twin stau LSP. The twin stau LSP is either mostly right-handed or a highlymixed left-right state. It can never be strongly dominated by the left-handed component because in such a case twin sneutrino would be the LSP.

Experimental constraints
In this section, we discuss the experimental and observational constraints on the twin stau LSP. The resulting constraints on the parameter space of the theory are presented in Sec. 5.

Thermal abundance
We will focus on a mass spectrum such that the LSP is the twin stau and the twin Higgsino is not too heavy in order to obtain the EW scale without fine-tuning. We also allow the bino to be light, since a light bino has important phenomenological implications. The wino can also affect the relic density but its effect is minor in comparison with the bino except for the parameter space where the LSP has an O(1) left-handed twin stau component. We thus mainly take a large wino mass for simplicity, but comment on the impact of a light wino. All the other twin and MSSM particles (squarks, first and second generation sleptons, and gluino) do not affect the relic density computation as long as they are not degenerate with the stau so that coannihilation with them is ineffective. We mostly consider the case where the twin and the MSSM tau Yukawa couplings are equal, but also comment on the case with different couplings.
The dominant twin stau annihilation channels resemble those of the MSSM stau [33][34][35] with the MSSM states replaced by the corresponding twin states. The twin stau annihilation channels depend mainly on the size of the left-right mixing. For the righthanded dominated twin stau the dominant annihilation channel is into twin photons unless the twin bino is light in which case annihilation into twin taus may become dominant. Finally, in the presence of a larger left-handed twin stau component of the LSP, the annihilation into twin W 's becomes relevant.
We compute the relic density with our modified version of Micromegas [36][37][38]. In our setup Micromegas automatically takes into account the coannihilation within particles in the twin sector. Coannihilation with the MSSM sector becomes relevant for the relic density computation if there are states with masses relatively close to the twin stau, ∆ i = (m i − mτ )/mτ 0.05. In the simplified case in which all superpartners except for MSSM and twin staus are decoupled, only the MSSM lightest stau and the twin sneutrino may satisfy this condition.
In order to capture the effect of the coannihilation with the MSSM stau, we use the approximation that the annihilation cross section σ of the twin stau and the MSSM stau are the same, which is justified by the approximate Z 2 symmetry. We also neglect mixed annihilation cross-section between the MSSM stau and the twin stau. 2 We then compute the effective annihilation cross section following [39], which is given by where x f = mτ /T f ∼ 25, with T f the freeze-out temperature. Assuming that the effects of coannihilation on the freeze-out temperature are subdominant, the relic density can be simply computed by the scaling of the abundance without the coannihilation Ωh 2 by the cross sections, Ω coann h 2 = Ωh 2 σ/σ eff . Within this simplified scheme, the effects of the coannihilation can be controlled by the masses of the spectrum alone. The most interesting case is when the thermal relic abundance of the twin stau LSP computed above, Ω th h 2 , fits the observed density of dark matter Ω obs h 2 ≈ 0.12 [40]. We refer to this case as thermal DM. Nevertheless, in analysing the constraint from and sensitivity of DM direct detection, we also consider two other possible scenarios in which the abundance of the twin stau LSP, Ωτ , may be different from its thermal abundance: • Non-thermal DM with Ωτ = Ω obs , which may be achieved in a non-standard cosmological history with non-thermal DM production (for Ω th < Ω obs ) or entropy production (for Ω th > Ω obs ).
• Multi-component DM with Ωτ = Ω th < Ω obs . In this scenario it is assumed that apart from the twin stau LSP there exist an additional component of DM which ensures that the total relic density fits the observed value.
See [41] for the up-to-date estimation of non-thermal DM abundance with the finiteness of the thermalization time scale properly taken into account. The non-standard cosmological history can solve the dark radiation problem in the Twin Higgs mechanism as we discuss in Sec. 6.

Twin stau self-interactions
The twin stau possesses a twin electromagnetic charge, and hence DM in our scenario has long-range self-interactions. This has interesting implications on galaxy formation. In particular, the strongest constraint comes from the triaxial structure of galaxy halos. The self-interactions can create a more isotropic dark matter velocity distribution in two ways. On one hand, interactions with large momentum exchange may immediately reduce anisotropy significantly. On the other hand, the cumulative effect of many soft interactions makes the velocity distribution more isotropic. The resulting constraints have been analysed in Ref. [42] for a Dirac fermion charged under a dark U(1) gauge symmetry, using measurements of the non-zero ellipticity of the gravitational potential of NGC720. Those results can be straightforwardly applied to the case of twin stau DM since the selfscattering cross-section of the twin stau in the non-relativistic limit is the same as that for a Dirac fermion. Assuming that the twin and the SM electromagnetic gauge couplings are equal, the ellipticity constraint sets a lower bound on the twin stau mass of about 210 GeV. As we will see, the twin stau mass leading to the correct thermal abundance of DM is above this bound. Notice, however, that this bound should be taken with caution.
Previously, the same data have been analysed in Ref. [43] resulting in a much stronger lower bound on the twin stau mass of about 1 TeV. However, Ref. [42] has made several improvements in extracting the constraint with respect to Ref. [43]. In particular, Ref. [43] neglects the saturation effect in the evolution of the ellipticity, approximating the change of the ellipticity linearly. Moreover, the constraint should be taken at larger radius with lower DM density, further suppressing the bound. A detailed analysis of the differences between these constraints can be found in [42]. We will use the more recent bound, that however should be taken with caution. Indeed, there are many challenges to the ellipticity measurements on DM self-interactions.
For smaller dark matter abundance, the characteristic timescale for a particle to change its kinetic energy by an O(1) factor, interpreted as the timescale on which the velocity completely randomizes, becomes longer. Therefore, for Ωτ h 2 0.12, the lower bound on the twin stau mass from self-interactions is even below 210 GeV.

Direct detection
A significant portion of the parameter space of the model can be probed via directdetection experiments. The interaction of the twin stau DM with nucleons is mediated by the Higgs portal and kinetic mixing. 3

Higgs portal
The twin stau couples to the twin Higgs, which mixes with the SM Higgs. The effective Higgs-twin stau coupling reads where s x and c x stands for sin x and cos x respectively, tan β = v u /v d is the ratio between the VEVs acquired by H u and H d , and α is the mixing angle necessary to diagonalize the neutral CP-even Higgs mass matrix. The masses of the twin tau, W, and Z are denoted by m τ , m W and m Z respectively, while g is the electroweak coupling constant. The second line is obtained in the decoupling limit by using cos(β − α) = 0. The twin stau mixing angle θτ diagonalizes the twin stau mass matrix and therefore relates the mass eigenstates τ 1 andτ 2 toτ L andτ R via an orthogonal transformation The twin stau thus interacts with nucleons via Higgs exchange, with a cross section where µτ N is the nucleon-twin stau reduced mass, mτ is the twin stau mass, and c N is the twin stau-nucleon effective coupling. The sums run over q = u, d, s and Q = c, b, t and the coefficients c q,Q are given by for up-and down-type quarks, respectively.

Kinetic mixing
The twin photon may have a kinetic mixing with the SM photon, The twin stau interacts with nucleons via the kinetic mixing, and XENON1T [44] gives a constraint [45] < 10 −10 mτ 100 GeV In the next section, we will consider a case where the bino and the twin bino have a mass mixing. For such a case, the strong upper bound on the kinetic mixing requires a structure of the theory such that the mass mixing is sizable while the kinetic mixing is suppressed.

Collider
The main collider signature of this scenario is the existence of relatively light staus in the range of several hundred GeV. The experimental constraints on the stau mass significantly depends on the lifetime of the stau. LHC constraints are absent for promptly decaying staus, because the mass splitting between the stau and the twin stau LSP is small (few tens of GeV at most), resulting in soft taus in the final states. Moreover, the lightest stau in this scenario is dominated by the right-handed component and the LHC is not yet sensitive to direct production of right-handed staus even in the limit of massless LSP [46]. We should also note that in the present scenario the decays of the lightest stau are threebody, as shown in Fig. 1, in contrast to the MSSM-motivated LHC stau searches assuming two-body decays of the stau, so the taus in the final states are even softer. LHC constraints on the wino and Higgsino that decay into staus may also be relevant, but they are not constraining as long as staus decay promptly. The wino decay proceeds as follows: the wino decays into a tau/neutrino and a left-handed slepton, which subsequently decays into Z/W/h and a right-handed stau, as shown in Figs. 2a and 2b, with the righthanded stau decaying into a tau, a twin tau, and a twin stau, see Fig. 1. The search that exploits taus and missing energy does not constrain the parameter region where the LSP mass is above 100 GeV [47], which is the region of parameter space relevant for twin stau DM. The search using Z/W/h and missing energy could give a stronger constraint, but the lower bound on the wino mass is 600 GeV even with a massless LSP, and the constraint is absent if the LSP mass is above 300 GeV [47,48]. In our setup, the left-handed slepton is boosted, so the Z/W/h tends to be aligned with the tau from the right-handed stau decay, and leptons from Z/W/h may be vetoed by the isolation criteria imposed in the search, weakening the bound. If the first two generation sleptons are also lighter than the wino, the lower bound on the wino mass becomes stronger and is 1300 GeV in the limit of massless LSP. This limit, however, is absent for the LSP mass above 800 GeV [47].
The Higgsino decays into a tau and a left-handed slepton, or a tau and a right-handed slepton, see Fig. 2c. Both give the same signal as the wino discussed above, except that the decay into a tau and a right-handed lepton does not involve Z/W/h. Because of the smaller cross section, the constraint on the Higgsino mass is even weaker. Even with the first two generation of sleptons lighter than the third generation, the bound on the Higgsino mass is expected to be weak since the Higgsino dominantly decays into a tau and a stau.
The left-handed stau is required to be light so that the twin stau becomes the LSP via left-right mixing. The left-handed stau decays into Z/W/h and a right-handed tau, which subsequently decays into a tau, a twin tau, and a twin stau. The topology of the event is similar to that of the wino decay, but the constraint on the left-handed stau mass is absent because of the smaller cross section.
Finally, for sufficiently long-lived staus the LHC constraints on the stau mass and the masses of other MSSM particles that decay into stau become more important. Since the lifetime of the stau affects not only the collider phenomenology of the model but also the relic density computation, we discuss it in some more detail in the next section.

The lifetime of stau
In this section, we discuss the lifetime of the stau, which affects the phenomenology in several ways.
• The stau may be observed as a long-lived charged particle at collider experiments.
The LHC searches for long-lived charged particles such as stau are the most sensitive if the lifetime is longer than about 10 ns (that corresponds to the decay length of O(1) m) [49]. In such a case the lower bound on the stau mass is about 430 GeV for the right-handed stau [50]. Most of the parameter region with successful thermal freeze-out is then excluded, as we will show in section 5. Furthermore, since the stau is also produced from the decay of heavier MSSM particles, their mass is also bounded from below. For example, the stop, gluino, wino, Higgsino, and left-handed slepton masses must be above about 1600, 2200, 1300, 1200, and 500 GeV, respectively.
For the stau decay length between about 10 cm and 1 m LHC searches for disappearing charged tracks may become sensitive [51,52]. The upper limit on the production cross-section at the LHC with √ s = 13 TeV can be as small as about 10 fb for the decay length O(1) m [51] which corresponds to the production cross-section of the 200 GeV right-handed staus. 4 Thus, the disappearing track searches are not yet sensitive to our scenario. To sum up, the LHC constraints are relevant only if the stau decay length is above O(1) m and put lower bound on the masses of MSSM particles.
• The decay product of the stau contains a tau. If the lifetime is longer than 10 5 s, the abundance of the stau at the time of decay must be much smaller than the dark matter abundance, so that Big-Bang Nucleosynthesis (BBN) is not disturbed [53].
Since the annihilation rate of the stau is comparable to that of the twin stau, the stau abundance before it decays is of the same order as the twin stau abundance. As a consequence the twin stau cannot be DM. We found, however, that the lifetime of stau is generically much shorter than 10 5 s and BBN does not constrain our scenario, unless the stau is nearly degenerate with the twin stau.
• If the decay of the stau occurs after the freeze-out of the twin stau annihilation around T mτ /25, the abundance of the twin stau is enhanced. Since the stau abundance before it decays is of the same order as the twin stau freeze-out abundance, the final dark matter abundance is larger by a factor of about 2 with respect to the freeze-out value. However, the coannihilation between stau and twin stau is no longer important if the stau decay is ineffective. This occurs if the lifetime is longer than We first discuss the case in which the stau lifetime is determined by the Higgsino-twin Higgsino mass mixing and significantly constrains the parameter space of the model. We then discuss possible bino-twin bino mass mixing and show that the lifetime can be short enough that none of the above issues arises. In both cases, the mediators of the decay include the bino. We assume that the lighter stau is dominantly the right-handed one. When the lighter stau contains a significant fraction of the left-handed twin stau, the decay rate becomes smaller by a factor of up to 16. The suppression, however, can be compensated by the change of parameters by O(1) factors.

Higgsino mass mixing
Since the Higgs and the twin Higgs couple with each other via a quartic term, we expect that the Higgsino also couples with the twin Higgsino. As a result the neutral Higgsino has a mass mixing with the neutral twin Higgsino. The magnitude of the mixing depends on the UV model. In order to be concrete, we consider a D-term Twin Higgs model where the mixing is given by Here the mixing H is generated by the exchange of a heavy gaugino and is as small as vv m soft /v 2 S , where v S is the symmetry breaking scale of the gauge symmetry whose D term is responsible for the SU (4) invariant Higgs quartic potential and m soft is a soft SUSY breaking mass of the symmetry breaking field. For natural choices of parameters we found that typically H ∼ O(0.01 − 0.1)v . A detailed computation is given in Appendix A.
The main stau decay channel isτ →τ † τ τ . (Another decay channelτ →τ τ τ † is p-wave suppressed.) In the limit of large tan β and µ, M 1 g v , we find the decay length where M 1 is the (twin) bino mass parameter. Here we neglect the tau and twin tau masses, but we include them in our numerical results. The decay rate is suppressed for large |M 1 | and |µ|. In order for the decay length of the stau to be below 1 m via the Higgsino mixing, M 1 and µ must not be much above the stau mass. The requirement that the stau decays before freeze-out can be achieved with an interplay between the different parameters involved in the decay rate: M 1 , µ, tan β, m 3L , m 3R , and A τ . The parameter M 1 can be as small as the stau mass. However, this unavoidably leads to twin stau annihilation via twin bino exchange which leads to a larger mass for the twin stau LSP consistent with Ωh 2 ≈ 0.12. The parameters µ, tan β, m 3L , m 3R and A τ are all related and affect the decay width with competing effects. The µ parameter is bounded to be below ∼ 1.5 TeV by the need to solve the naturalness problem. In principle, small µ implies an enhanced decay rate. However, the smaller µ the smaller is the mixing between the left-and right-handed twin staus. This makes the twin stau heavier than the MSSM counterpart (if mostly right handed) or than the twin sneutrino (if mostly left handed) and the LSP is not the twin stau. This can be compensated by requiring large tan β or a Z 2 breaking in the tau Yukawa couplings. We will show the allowed parameter space in section 5.

Bino mass mixing
The bino and twin-bino may have mass mixing, A stau can then decay into a tau, a twin stau, and a twin tau via bino exchange. The effective operator describing the decay is whereτ is the right-handed tau and we use the approximation M 1 mτ . The lifetime of the stau is Γ −1 = 2 ns mτ 300 GeV and may be easily shorter than 10 ns, so that the stau decays before the freeze-out and the collider constraints from long-lived particle searches are evaded.
While the mass mixing parameter M must be large enough to make the stau shortlived enough, the photon-twin photon kinetic mixing must be suppressed to avoid too large a twin stau-nucleon scattering cross section, cfr section 3.3.2. This can be achieved by a supersymmetry breaking field Z with a coupling and VEV, where M * is a cut-off scale. Although Z ∼ M * may be avoided by a discrete symmetry of Z, non-zero VEV of Z may be induced by supergravity effect and quantum corrections. We investigate those effect for an explicit example in Appendix B and find that the VEV is sufficiently small and the kinetic mixing is much below the constraint from DM detection experiments shown in Eq. (7).

Results
In this section we assemble together all the bounds described above and show representative plots. We start by considering the case of a decoupled (twin) bino and wino, with  a Higgsino at ∼ 1.5 TeV and tan β = 10. Fig. 3 shows the contours of Ω th h 2 = 0.12 (blue) and the twin stau masses (black) in the twin stau soft masses parameter space. The correct thermal relic density is achieved for twin stau masses between about 300 and 400 GeV. In particular, in the region where the twin stau is mostly right-handed, coannihilation with the MSSM stau plays an important role in decreasing the effective cross section. As a consequence the thermal density is larger and Ω th h 2 = 0.12 is obtained for smaller values of the SM mass. It is remarkable that the whole region preferred by the thermal abundance of the twin stau is consistent with astrophysical constraints on twin stau self-interactions which exclude masses below about 200 GeV. Part of the parameter space is excluded simply due to the mass spectrum. The purple region leads to tachyons. In the red areas the twin stau is not the LSP; in the lower right region the twin sneutrino becomes the LSP, while in the upper left region the stau is lighter than the twin stau.
The yellow region is excluded by Xenon1T [54], while the green region shows the potential reach of LUX-ZEPLIN (LZ) [55]. The dashed yellow (green) curve shows the direct detection bound from Xenon1T (potential reach of LZ) assuming thermal abundance. The yellow (green) solid curve, on the other hand, assumes that Ωh 2 = 0.12 in the whole parameter space. This includes the case in which DM is produced non-thermally and the case in which the thermal abundance is too large but is reduced to the observed one because of entropy production. This plot shows that, in the scenario in which the thermal abundance of the LSP explains all of the dark matter, most of the parameter space is allowed by Xenon1T but most of it could be probed by LZ. In the thermal dark matter scenario with under-abundant LSP, which happens for twin stau masses smaller than ∼ 300 GeV, Xenon1T excludes significant portions of parameter space while LZ could probe most of it. Note, however, that LZ will not be sensitive enough to probe some region of parameter space corresponding to the LSP dominated by the right-handed twin stau. This region tanβ Ω th h 2 = 0.12  surrounds a blind spot in direct detection corresponding to a vanishing twin stau-Higgs coupling due to cancellation between the first and the last term in the squared bracket in eq. (2). However, in this region the decayτ →τ † τ τ is kinematically closed and the stau lifetime is generically too long to be consistent with the BBN. 5 We should note, however, that for the parameters chosen in Fig. 3, in the absence of bino mass mixing, the decay length would be much above 1 m because of the large µ and M 1 . As a consequence, the stau would decay after the freeze-out, which can enhance the relic density. When the MSSM and twin staus are nearly degenerate with each other the contours of Ω th h 2 remain almost unchanged while when their fractional mass difference is more than 1/x f , the enhancement factor is about 2 and requires a smaller stau mass. Furthermore, mτ < 430 GeV is excluded by the LHC searches for charged long-lived particles, as explained in section 4. Therefore, thermal dark matter case would be excluded. In order to avoid these conclusions it is enough to introduce the bino-twin bino mixing parameter M ∼ O(10 −1 ), see eq. (13), since a typical mass splitting between twin stau and stau is O(10) GeV.
In Fig. 4, we show contours of Ω th h 2 = 0.12 (blue) and mτ (black) in the m 3R vs µ (left) and m 3R vs M 1 (right) plane. We fix M 2 = 10 TeV, v /v = 3 and m 3L = m 3R + 100 GeV. In the left panel, in order to keep the stau mixing parameter fixed as µ varies, we take tan β = 15/µ TeV, while M 1 = 1 TeV. The left panel shows that decreasing µ has no impact on the relic abundance of the twin stau unless the twin Higgsino is almost degenerate with the LSP, for which co-annihilation with the twin Higgsino becomes effective and the twin stau mass reproducing the observed dark matter abundance can be as large as about 450 GeV. For smaller µ, the twin Higgsino becomes the LSP. The red curves show the contours of the decay length dτ = 0.1, 1, and 10 m via the Higgsino mixing discussed in Sec. 4. A large portion of the parameter space with µ 600 GeV has dτ below 1 m. In such a region bino mass mixing is not necessary. We also note that the magnitude of the Higgsino mixing can be larger than our benchmark value by an O(1) factor, which can decrease the decay length by an O(10) factor and further expands the viable parameter space without resorting to bino mass mixing. For larger µ one has to still rely on bino mass mixing, although M can be somewhat smaller, down to about 10 −2 , due to smaller M 1 than in Fig. 3. In the right panel of Fig. 4 we show the dependence on M 1 for µ = 1 TeV. The presence of a light bino enhances the twin stau annihilation cross-section into twin taus. As a result, for a sub-TeV bino Ω th h 2 = 0.12 is obtained for a larger twin stau mass up to almost 500 GeV. In order to have the stau decay length below 1 m without introducing bino mixing, M 1 600 GeV is preferred for µ = 1 TeV. LZ will probe both the thermal and non-thermal scenario. In the non-thermal scenario, the region where the twin stau mass is above about 800 GeV will not be tested by LZ, but it is generically in conflict with the BBN. To suppress the stau decay length without introducing bino mass mixing, the bino mass is preferred to be below 800 GeV.
In Fig. 5 we show how the stau parameter space opens up for µ = M 1 = 700 GeV assuming M 2 = 10 (1) TeV in the left (right) panel. In order to have large enough left-right stau mixing, we compensate the smaller value of µ by a large tan β = 30. A large portion of the parameter space has a stau decay length below 1 m and is allowed by Xenon1T. Moreover, in the region with the stau decay length above 1 m the mass of the twin stau playing the role of thermal DM is above 430 GeV, even after considering the enhancement of the relic abundance due to stau decays after the freeze-out. As a consequence, this region is allowed by the LHC constraints on long-lived charged particles. However, this requires also taking µ and M 2 above about 1.2 and 1.3 TeV respectively (increasing µ does not affect the relic abundance of the twin stau, and doing so for M 2 also does not in the parameter region with a stau decay length above 1 m), otherwise the production rate of long-lived staus from Higgsino and wino decays would exceed the rate allowed by the LHC searches. From the plots in Fig. 5 we can also see how a light wino impacts the relic abundance of the twin stau LSP. As expected, the contour of Ω th h 2 = 0.12 does not change with a lighter wino as long as the LSP is dominated by the right-handed twin stau. However, in the vicinity of the twin sneutrino LSP region the left-handed stau component is nonnegligible and the 1 TeV wino increases the twin stau mass corresponding to Ω th h 2 = 0.12 up to about 550 GeV, compared to about 450 GeV in the case of decoupled wino. This also helps to avoid the Xenon1T constraint in that region and only a small part of the Ω th h 2 = 0.12 contour is excluded.
We have been assuming so far strict Z 2 symmetry in the Yukawa couplings. Let us now discuss what happens if Z 2 breaking in the tau Yukawa coupling is present. This breaking does not spoil the TH solution to the hierarchy but can help to solve the dark radiation problem of TH models [17,18]. A larger twin tau Yukawa coupling increases the left-right mixing of twin staus. Moreover, the larger mixing increases the mass splitting between the twin stau LSP and stau, which suppresses the stau lifetime. One consequence is that smaller values of tan β are now viable. The impact of larger y τ is shown in Fig. 6 for tan β = 5. For y τ /y τ above about two, the twin stau can be the LSP. However, in order to have the correct thermal relic abundance and evade the BBN bounds, the ratio should be above four. Furthermore, the stau decay length can be sufficiently small for most of the parameter space with thermal DM without introducing bino mass mixing.
An important feature of the scenario with an enhanced twin tau Yukawa coupling is that a larger twin stau mass leads to Ω th h 2 ≈ 0.12, consequence of the increased annihilation cross-section into SM-like Higgses and twin tops. The twin stau mass may be up to about 600 GeV, above which the twin neutralino becomes the LSP. However, there is an upper bound on y τ /y τ (of about six in the thermal scenario) set by Xenon1T constraints. This is because the direct-detection cross-section is enhanced by the larger twin tau Yukawa coupling. Nevertheless, even after taking into account the Xenon1T constraint the twin stau thermal DM can exceed 500 GeV. The whole thermal scenario will be probed by LZ also in this case.
The range of allowed y τ /y τ depends on tan β. For smaller (larger) values of tan β the lower bound on y τ /y τ gets stronger (weaker) while the upper bound from Xenon1T gets weaker (stronger). The stau decay length in the parameter space with thermal twin stau DM varies by orders of magnitude even without introducing bino-twin bino mixing. Some of the parameter space will be probed by disappearing track searches at the LHC but promptly decaying staus are possible for large enough y τ /y τ .

Dark radiation and dark matter abundance
It is well known that Twin Higgs models in their simplest realizations suffer from too large amount of dark radiation. In particular, in Mirror Twin Higgs models ∆N eff ≈ 5.6 [58]. This is order of magnitude above the upper bound from Planck of about 0.3 at 95% C.L. [40]. Possible solutions to this problem are: • late entropy production which dilutes dark radiation [58,59], • Z 2 breaking in the Yukawa couplings of light fermions [17,18], which reduces the entropy of the mirror sector by the Boltzman suppression, • Fraternal Twin Higgs [23] where the first and second generation mirror fermions are removed (the mirror photon, however, may be still in the spectrum).
In this section, we discuss the compatibility of the parameter space of the twin stau LSP DM with these solutions. In Figs. 3, 4, and 5, we assume Z 2 symmetric yukawa, so only late time entropy production and Fraternal Twin Higgs are applicable. In the Fraternal Twin Higgs solution, the thermal abundance is maintained, so the prediction shown by the blue lines is applicable.
With entropy production, the twin stau abundance is affected. If the production of the twin stau is negligible during entropy injection, which is guaranteed when the particle responsible for it is lighter than the twin stau, the twin stau abundance is simply diluted. As a consequence, the abundance consistent with the observed value would require a larger mass of the twin stau DM. The thermal abundance is roughly proportional to m 2 LSP . Assuming that ∆N eff is diluted by a factor D 4/3 , the DM abundance is reduced by a factor D and we can link ∆N eff with the twin stau LSP mass that gives the observed DM abundance after taking into account entropy production: where m th is the twin stau mass that gives thermal abundance equal to the observed value. Here we consider the limit where the heating of the twin sector is negligible during the entropy production. If not, ∆N eff is even larger. We see that in order to satisfy the Planck constraint on ∆N eff , the twin stau DM mass would be around 1 TeV in the case of decoupled bino and even heavier for lighter bino. Interestingly, ∆N eff will be tested with much better precision and the bound on ∆N eff is expected to be set by CMB-S4 experiments [60] around 0.06 at 95% C.L. [61]. If CMB-S4 does not find signals of dark radiation this would set the lower bound on the twin stau mass of almost 2 TeV. This would result in an even heavier Higgsino, leading to fine-tuning of the EW scale worse than 1%.
If the particle responsible for entropy production is heavier than the twin stau, twin staus may be produced during entropy production. The stau mass may be smaller than that predicted from Eq. (15). It is even possible that the stau mass is smaller than the prediction of the thermal DM scenario shown by blue lines.
In Fig. 6, since we assume Z 2 breaking yukawas, all of the solutions are applicable, and the case with full mirror sector and without entropy production would be the minimal one. Then the thermal abundance is maintained and the prediction shown by the blue lines is applicable. One can in addition introduce entropy production, for which case the stau mass may be different from that prediction. Because of the Z 2 breaking in the yukawa, Eq. (15) is no longer applicable.

Conclusions
We have investigated the possibility that DM in the Universe is a charged particle in supersymmetric Twin Higgs models. We found that the twin stau is a perfectly viable candidate for such a DM particle. In order to ensure that the twin stau is the lightest supersymmetric particle and hence stable, it must be a mixture of the right and left components. We found that the thermal relic abundance of the twin stau LSP matching the observed DM density implies that the twin stau mass is generically in the range of 300-400 GeV, but it can even exceed 500 GeV if some of (twin) neutralinos are not much heavier than the LSP. Intriguingly, the preferred range of the twin stau mass is above a lower bound of about 200 GeV from twin stau self-interactions arising from the non-zero twin electromagnetic charge of the twin stau.
Signatures of twin stau DM are very different from standard supersymmetric DM candidates. The scattering cross-section off nuclei for the twin stau is naturally suppressed because the twin stau interacts with nuclei only via the mixing of the twin Higgs with SMlike Higgs. As a result, the latest constraints from direct detection experiments are satisfied in most of the parameter space in contrast to standard supersymmetric scenarios in which fine-tuning of the parameters is necessary (coannihilation [33,34] and blind-spots [62][63][64] require fine-tuning in the parameters relevant for the DM phenomenology, while a pure Higgsino or wino requires fine-tuning of the EW scale.). Nevertheless, the Higgs-twin Higgs mixing cannot be arbitrarily small in natural Twin Higgs models so a signal of twin stau dark matter is expected at the LZ experiment. The main distinctive feature of this scenario are self-interactions of the LSP which is not present in standard supersymmetric models and future astrophysical data may differentiate between these classes of DM.
In this scenario the mass splitting between the twin stau LSP and the MSSM stau is generically small, below few tens of GeV. This may lead to long lifetime of the stau which, in turn, may have important implications for cosmology and collider searches for staus. The lifetime of the stau is model-dependent and can be short enough to avoid all cosmological and collider constraints. Nevertheless, in some of the parameter space one may expect signals of long-lived staus at the LHC, e.g., in disappearing track searches.
The gauge symmetry is broken by VEVs of S andS with charges 1/2 and −1/2, respectively. The VEVs may be induced by a superpotential where Y is a chiral multiplet. We assume the same soft masses of S andS so that S = S ≡ v S (= O(M )) and the D-term potential of U (1) X does not give large Higgs masses. The mass of the X gauge boson is g X v S . Y obtains a tadpole term ∼ λm soft v 2 S Y via supersymmetry breaking. It also obtains a positive mass squared = 2λ 2 v 2 S |Y | 2 . From these two terms, λv Y ∼ m soft .
The mass terms of relevant fermions are where mX is the SUSY breaking Majorana mass ofX. After integrating outX,S, and S, the mixing between the MSSM and twin Higgsinos is given by To have a non-decoupling D-term potential of the SM and twin Higgses from U (1) X , we need m soft ∼ m X = g X v S . The parameter H in Eq. (9) is then estimated as where we used λv Y ∼ m soft and neglected mX . The electroweak precision measurements requires that v S > 4 TeV [12]. In our numerical calculations of the stau decay length we take, for concreteness, g X = 2, v S = 5 TeV and λv Y = m soft = m X = 10 TeV which leads to H ≈ 20 GeV(v /v)/3 but we should keep in mind that H can vary by O(1) factor depending on the soft masses of the U (1) X symmetry breaking sector.

B Generation of bino-twin bino mixing
In this appendix, we analyze a setup with a significant bino-twin bino mass mixing but with a negligible kinetic mixing. We introduce a gauge singlet chiral multiple Z and the following interaction, This interaction preserves the following approximate Z 2 × Z 2 symmetry, where V is a hypercharge gauge multiplet. We may give a non-zero F term to Z while fixing it near the origin via, for example, the O'Raifeartaigh model [65], where X and Y are chiral multiplets. The Z 2 symmetry is promoted to Z 4 symmetry with Z → −Z, X → iX and Y → −iY . A dimensionful parameter Λ is a spurion of Z 4 symmetry breaking. We may obtain it dynamically by embedding this into the IYIT dynamical supersymmetry breaking model [66,67], where m ∼ Λ is also generated dynamically. Z obtains an F term of Λ 2 while is fixed at the origin by the quantum correction from X and Y , Z obtains a mass term ∼ k 4 Λ 4 /(16π 2 m 2 ) and a tadpole term ∼ Λ 2 m 3/2 Z, so obtains a non-zero VEV ∼ 16π 2 m 2 m 3/2 /(k 4 Λ 2 ). For k ∼ 1 and m ∼ Λ, the resultant kinetic mixing is much smaller than the experimental upper bound. Explicit breaking of the symmetry also induces a non-zero field value of Z. Indeed, the action on the hypercharge gauge field is charge-conjugation symmetry, which is explicitly broken by the electroweak gauge interaction. Quantum corrections to D terms of V and V gives Zm 4 soft /(256π 4 M * ). There also exists a tree-level contribution from the Higgs VEVs ∼ Zv 2 v 2 /M * . The VEV of Z induced by these tadpole terms is sufficiently small.