Cosmology of the Twin Higgs without explicit $\mathbb{Z}_2$ breaking

The cosmology of the Twin Higgs requires the breaking of the $\mathbb{Z}_2$ symmetry, but it is still an open question whether this breaking needs to be explicit. In this paper, we study how the Mirror Twin Higgs could be modified to be compatible with current cosmological constraints without explicit $\mathbb{Z}_2$ breaking. We first present a simple toy model that can realize baryogenesis without explicit $\mathbb{Z}_2$ breaking or reaching temperatures that would lead to domain walls. The model can also either solve the $N_{\text{eff}}$ problem and bring the abundance of mirror atoms to an allowed level or provide the correct dark matter abundance. We then present another simple model that leads to mirror neutron dark matter and thus acceptable dark matter self-interactions. We also include in appendix a series of results on energy exchange between different sectors that might prove useful for other cosmological problems.


Introduction
The Twin Higgs [1,2] attempts to solve the little hierarchy problem by introducing partners that are neutral under the Standard Model (SM) gauge groups and is the prime example of Neutral Naturalness. Its simplest version is the Mirror Twin Higgs. In this model, a copy of every SM field is introduced. The principal difference is that these mirror partners are instead charged under new gauge groups that reflect those of the Standard Model. The Higgs doublet and its mirror partner can then be combined to write a potential that respects an approximate global SU (4) symmetry. Its spontaneous breaking to SU (3) results in seven (pseudo)-Goldstone bosons. Three are eaten by the massive SM gauge bosons and three by their partners. The remaining one corresponds to the experimentally observed Higgs boson. Its mass is protected at one loop by a Z 2 interchange symmetry which ensures that the leading correction to the potential respects SU (4) and hence does not contribute to the Higgs mass directly. The latter is then effectively protected by neutral partners. This symmetry imposes an equality between the Yukawa and gauge couplings of the two sectors. Without unacceptable tuning, the mirror partners are typically a factor of a few heavier than their SM equivalents. It is a well established fact that the Z 2 symmetry must be broken for the Twin Higgs to be compatible with the Higgs signal strengths measurements. During the early days of the model and often still to this day, this was done by introducing explicit soft Z 2 breaking. The possibility of doing without this explicit Z 2 breaking and having the symmetry only broken spontaneously is both aesthetically appealing and likely to facilitate UV completions. Refs. [3][4][5][6][7] demonstrated that, at least as far as collider constraints are concerned, it is possible to do so and in addition the amount of tuning required decreases.
At the same time, the Twin Higgs model can be consistent with cosmology, which is already the subject of a considerable literature [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. 1 In a similar fashion to the Higgs signal strengths, the cosmology of the Twin Higgs also requires Z 2 breaking. It was amply demonstrated that viable cosmology can be obtained via explicit Z 2 breaking. However, one question that is still unanswered in the literature is whether cosmology requires this breaking to be explicit. Indeed, all previous works either included some sort of explicit Z 2 breaking, be it some gauge or Yukawa couplings being different, some particles being absent in one sector or even more esoteric possibilities, or never presented a full model that could address all known issues while providing an adequate dark matter candidate.
The fact however is that the cosmology of the Twin Higgs with explicit Z 2 breaking already presents major challenges. With only spontaneous breaking, these challenges are exacerbated, as respecting this symmetry imposes additional constraints and reduces our set of tools to address them.
The first such challenge is baryogenesis. The Z 2 symmetry being only broken spontaneously almost unavoidably leads to domain walls, which may overclose the Universe. As long as inflation lasts long enough, their density can thankfully be brought to acceptable levels. In addition, domain walls will not be reintroduced during reheating as long as the reheating temperature does not reach the Z 2 restoration scale. However, most standard baryogenesis mechanisms take place at temperatures higher or not too far from the expected Z 2 restoration temperature. One obvious way to solve this apparent conflict is for baryogenesis to take place below the Z 2 restoration scale. However, this could be challenging in general for mechanisms like electroweak baryogenesis [37][38][39] or leptogenesis [40]. This is even more difficult if one wishes for dark hadrons to represent dark matter via some realization of Asymmetric Dark Matter [41][42][43].
The second challenge is the contribution of the mirror photon and mirror neutrinos to the effective number of relativistic degrees of freedom N eff , which is severely constrained by both the Cosmic Microwave Background (CMB) and Big Bang Nucleosynthesis (BBN) [44,45]. With explicit Z 2 breaking, such particles can be removed or made heavier, but this is non-trivial when Z 2 is not broken explicitly. Note however that some existing solutions to the N eff problem only need some expectation values to differ between the two sectors and could in principle be accommodated without explicit Z 2 breaking (see e.g. Refs. [8][9][10]12]).
The third challenge is that, even if it could explain the observed dark matter abundance via Asymmetric Dark Matter, the standard Mirror Twin Higgs would lead to dark matter in the form of mirror atoms. The problem with this scenario is that these would display self-interactions similar to normal atoms. If dark atoms were to represent the entirety of dark matter, their self-interactions would be ruled out by orders of magnitude or would require tuning at an unacceptable level [46][47][48]. It is then crucial to be able to modify the model such that the dark matter takes a more acceptable form, such as mirror neutrons.
With this context in mind, the goal of the present paper is to study the feasibility of constructing cosmologically viable Twin Higgs Models without explicit Z 2 breaking. The construction of a full model is rather ambitious and we will instead limit ourselves to studying whether it is possible to individually solve the three challenges mentioned above. More specifically, we will study whether it is possible to realize baryogenesis without reintroducing domain walls, whether the same process can also generate the correct dark matter abundance and/or solve the N eff problem and whether dark matter can be converted into an acceptable form.
The end result will be that it is indeed possible to overcome these challenges. This will be demonstrated by presenting two different models. The first one includes two Majorana fermions. The heaviest one is assumed to dominate the energy content of the Universe at early times. It then decays and produces a net amount of both baryons and mirror baryons, as well as some amount of the lighter Majorana fermion. As the Universe expands, the lighter Majorana fermion comes to dominate the energy abundance. Because of kinematical reasons, it then decays mainly to the Standard Model sector, thus reheating that sector. The model can provide the correct matter abundance while maintaining temperatures that are low enough not to reintroduce domain walls. It can also either solve the N eff problem and reduce the abundance of dark atoms to an acceptable level or provide the correct dark matter abundance.
The second model solves the remaining problem of dark matter self-interactions. A set of vector quarks is introduced in each sector. These mix with their respective up quarks via Yukawa interactions with the Higgs. This mixing has the effect of adding to the mass of the up quark of a given sector a contribution proportional to the vev of the Higgs of that sector cubed. This results in the mass of the mirror up quark increasing faster than the mass of the mirror down as the vev of the mirror Higgs increases. As such, the mirror proton can be made considerably heavier than the mirror neutron. Dark matter then consists of mirror neutrons and the abundance of mirror atoms can be brought to negligible levels. The smallness of the mass of the up quark ensures that the required amount of mixing is small enough to be comfortably below any current experimental constraints.
The article is organized as follows. The first model is introduced, its mechanism explained, its constraints discussed and its parameter space studied. The second model is then introduced, its constraints discussed, its parameter space studied and alternative models presented. An appendix presents some useful results on energy exchange between different sectors for cosmological evolution. Additional appendices discuss the decay asymmetry, scattering asymmetries, the evolution equations, the Higgs signal strengths and the computation of the dark atom abundance.
2 Addressing baryogenesis, dark matter abundance and N eff We begin this paper by introducing a toy model which can potentially address baryogenesis, dark matter abundance and N eff while maintaining temperatures below the Z 2 restoration scale. The model serves as a proof of principle and it goes without saying that variations are possible. This section contains a description of the model, an explanation of the mechanisms involved, a discussion of the different constraints and some summary scans of parameter space. To avoid obscuring the discussion with technicalities, all mathematical details concerning the cosmological evolution are relegated to Appendices A, B, C and D.

Model summary
The field content of the model is as follows. First, a complete copy of the Standard Model is introduced. Fields from the SM sector are labelled with an A and those of the mirror sector with a B. 2 The Higgs doublets are labelled as H M and obtain expectation values with v B larger than v A by a factor of a few to satisfy Higgs signal strength requirements. How these vevs are acquired is irrelevant to the present discussion, but can be done via spontaneous breaking [3][4][5][6][7]. In addition, several fields without SM equivalents are introduced. Every field labelled by χ is a left-handed Weyl spinor and those labelled by φ are complex scalars. The fields are where we used the notation The Lagrangian containing the interactions relevant to us can be separated into two parts. The first one involves the fermions and can be written as This can be rewritten in terms of Majorana spinors N i as It is easy to verify that this Lagrangian allows for both baryon number and CP violation. The second part of the Lagrangian is responsible for providing different masses to φ A and φ B without explicit Z 2 breaking. This can be done in several ways. First, the Lagrangian could contain the term (2.5) Replacing the Higgs doublets by their expectation values will affect the masses of φ A and φ B differently. Alternatively, new scalar fields could be introduced and play a similar role to H M . This can be done for example by introducing another Higgs doublet and its partner or by introducing a real scalar and its partner. In these cases, the different vevs can be obtained again by spontaneous breaking of the Z 2 symmetry. Since there are so many possibilities and since this is sufficient, we will simply work with the effective Lagrangian (2.6)

Description of the mechanism
We now proceed to describe how this model solves the issues it is designed to address. We refer to Fig. 1 for illustration of the evolution of different relevant quantities in a given benchmark. The parameters are taken as λ 323 = 0.05, λ 431 = 0.05, λ 432 = 0.0005 e iπ/4 . (2.7) All unspecified λ 3ij and λ 4ij are set to zero. The initial density of N 2 is set to 10 7 GeV 3 and its temperature to zero. All other initial densities are set to zero. The mass of φ A is chosen to reproduce the correct baryon abundance. This benchmark also leads to an r T = T B /T A of 0.408, which satisfies the bounds on N eff as will be discussed in Sec. 2.3.  The abundance of dark baryons is Ω B = 6.02 × 10 −4 . Even assuming all dark baryons are mirror atoms, this is still considerably below experimental bounds, as will be discussed in Sec. 2.3. The value of v B /v A satisfies the Higgs signal strengths, which are discussed in Appendix E. The temperature of sector M is labelled as T M , its energy density as ρ M , its entropy density as s M , its net baryon density as ∆B M and This benchmark is not special and a summary exploration of the parameter space will be performed in Sec. 2.4. 3 The whole process can be separated into three qualitative phases. Initially, the energy content of the Universe is dominated by the heavier Majorana fermion N 2 . This situation could easily be realized through the decay of the inflaton if its coupling to N 2 is considerably stronger than its other couplings. Early on, N 2 starts to decay. It can decay to the A sector mainly via three channels: to three quarks, to three antiquarks or to N 1 , a quark and an antiquark. Similar decays to the mirror sector are also present. Because of the presence of a third decay channel in each sector and in conjunction with the Nanopoulos-Weinberg theorem [49], N 2 can present an asymmetry in its decay to baryons and antibaryons and similarly for its decay to the mirror sector. In practice, this comes from the interference of the diagrams of Fig. 7. The asymmetry in the B sector can be adjusted by changing the ratio m φ B /m φ A . Once the N 2 are mostly decayed, the Universe is populated with particles from the A and B sectors as well as some N 1 . Baryon asymmetries are also present in both sectors.
As time passes, the expansion of the Universe dissolves the energy densities of the different particles. As the A and B sectors are radiation dominated, their energy densities scale as a −4 , where a is the scale factor. Since N 1 is mostly non-relativistic, its energy density instead scales as a −3 and soon comes to dominate the energy abundance.
Finally, the N 1 population starts to decay. In principle, N 1 could decay to particles of either sector. However, there exists a sizable region of parameter space, which includes the benchmark, where the decay to the B sector is strongly suppressed because of kinematics. In the benchmark, N 1 can decay to the A sector as an off-shell top, a bottom and a strange. It however cannot decay to a mirror top, a mirror bottom and a mirror strange or even two mirror bottoms, a mirror strange and a mirror W as both the mirror top and mirror W are too heavy to be produced on-shell. This results in the N 1 population transferring its energy almost exclusively to the A sector and thus a relative reheating of that sector. This constitutes the main mechanism through which the N eff problem is solved. This is also why the N i were assumed to couple mainly to up-type quarks of the third generation, as having the decay of N 1 only being possible to one sector is easy to accomplish thanks to the large mass of the top quark. The decay of N 1 does not generate any sizable asymmetry and in fact partially dissolves the asymmetries by injecting entropy.
The end result of this mechanism is an A sector with a net population of baryons and a B sector with a much smaller net population of mirror baryons. This both explains baryogenesis and satisifies the bounds on dark matter self-interactions associated with the dark atoms. The fact that the mirror sector is much cooler also ensures that the N eff constraints are satisfied. The temperature of each sector is also maintained at all times considerably below the electroweak scale.

N eff
The number of effective relativistic degrees of freedom N eff is measured by Planck to be 2.99 ± 0.17 [45]. During BBN and the creation of the cosmic microwave background, the numbers of relativistic degrees of freedom in both sectors are about the same, which puts a limit on of ∼ 0.44 at 95% CL, where we used the fact that the SM value of N eff is 3.046.

Fraction of dark atoms
Without additional model building, dark matter would take the form of mirror atoms. However, a too large fraction of dark atoms X DA is excluded by limits on dark matter selfinteractions. Ref. [50] claims that this fraction can still be as high as about 10%, though the amount of uncertainty on this number is rather unclear. In addition, Ref. [14] claims that the limit on X DA might be brought to the few percent level in the not-so-distant future. When relevant, we will present contours of X DA and emphasize that the region above 10% is disfavoured.

Big Bang Nucleosynthesis
If N 1 is sufficiently long-lived, it will disturb BBN by injecting energetic hadrons and modify the observed abundances of primordial elements. Unfortunately, the cosmology of the model is rather exotic and no study of the decay of metastable particles during BBN perfectly mimics it. As such, we will simply ask that the lifetime of N 1 be below 0.1 s, which is the typical bound (see for example Refs. [51][52][53]). Considering that BBN limits are generally not strongly dependent on parameters such as the mass of the metastable particle and its branching ratio to hadronic channels, a more advanced treatment is not expected to change this constraint much.

Higgs signal strengths
The limits on the Higgs signal strengths are applied using the results of Appendix E.

Direct collider searches
The only new coloured particle in the model is the colour-triplet scalar φ A . Its pair production at the 13 TeV LHC in the mass range we consider ( 5 TeV) is negligible ( 1 event). The fermions N 1 and N 2 are gauge singlets and do not need to have any significant couplings that involve pairs of light quarks, so their direct production is irrelevant too. In part of the parameter space, N 1 can be produced in top quark decays, but the branching fraction is highly suppressed by the mass of φ A , the small couplings and phase space. Part of the B sector particles can be produced in Higgs decays (and escape the detectors invisibly), but the resulting effect on the visible branching fractions of the Higgs is too small to be seen in the current datasets (see Appendix E).

Parameter space and comments
We now provide some summary scans and comment on various properties of the model. All other λ couplings are set to zero. The initial density of N 2 is set to 10 7 GeV 3 and its temperature to zero. All other initial densities are set to zero.   a sufficient amount of dark matter, a sufficiently low ∆N eff or a sufficiently low dark atom abundance. These plots use the same parameters as Fig. 2, except for different values of v B /v A and m N 1 . The region excluded by BBN is outside the plot to the right.
As can be seen, there are regions of parameter space that can individually provide the correct matter abundance, the correct dark matter abundance or a sufficiently low ∆N eff . At low m φ A , the amount of normal matter is small because wash-out effects erase most of the asymmetry. At high m φ A , the amount of matter is also small because the decay asymmetry is suppressed. The amount of matter is therefore optimized for an intermediary value of m φ A . The same discussion applies to the B sector. The only difference is that the wash-out effects are less important because they involve the mirror top, which is very heavy and thus suppresses wash-out. The contribution to N eff simply decreases as m φ A and m φ B become larger, as this increases the lifetime of N 1 and decreases the efficiency of              processes that destroy N 1 around the time N 2 decays.
In addition, there are regions that can meet several of these requirements at the same time. Two of them are especially interesting. First, the blue/green region provides a sufficient amount of matter and dark matter. It however does not provide a sufficiently low ∆N eff and the dark matter self-interactions are too large. These two issues could however be addressed by additional model buildings. Second, the blue/purple/yellow region provides a sufficient amount of matter and a sufficiently low ∆N eff . It also has the benefit of leading to an amount of mirror atoms sufficiently low to pass the dark matter self-interactions bounds. This region is especially interesting, as it only requires an additional source of dark matter for it to provide a complete valid cosmology. Assuming the mechanism responsible for the production of this extra dark matter is over when N 2 starts to decay, our mechanism wouldn't be affected much by the small amount of dark matter that would lead to the current density. Otherwise, the consequences of the production of this extra dark matter are too model dependent to make a general statement. More extensive scans of the parameters of Eq. (2.10) and initial conditions did not reveal any region that could at the same time provide the correct matter and dark matter abundance with a sufficiently low ∆N eff . If such a region exists, it is most likely unnatural or requires a large reheating temperature. Different branching ratios of N 2 to the A and B sectors have relatively little impact on the final temperature ratio. This is because the A and B sectors acquire sufficiently high temperatures to reach thermal equilibrium via Higgs boson exchange. The initial temperature difference is simply erased. Also, this thermal equilibration between the A and B sectors makes it necessary for N 1 to decay after the two sectors have decoupled. However, the amount of time required for N 1 to dominate the energy abundance is typically much larger and generally controls the required lifetime of N 1 .
The process N 1 N 1 →b M b M is responsible for the destruction of a large fraction of the N 1 . This is why relatively so few of them are present in Fig. 1a immediately after the decay of N 2 , despite its relatively large branching ratio to N 1 b MbM .
The peak that can be seen in Fig. 1d at 10 −5 s corresponds to the B sector QCD phase transition quickly followed by the A sector QCD phase transition. Despite its rather striking nature, it affects relatively little the final results as the A and B sectors have decoupled by this point. If the decoupling of the two sectors had taken place between the two QCD phase transitions, this could have contributed to a partial solution of the N eff problem, which was explored in Ref. [8].
Arguably the most crucial question concerning baryogenesis is whether it can be done while maintaining temperatures low enough not to reintroduce domain walls. A rigorous answer to that question however requires the knowledge of the full potential and not simply v A and v B as we have only provided. As such, this question cannot be fully answered here. However, it can be seen in Fig. 2d that the reheating temperatures can be comfortably below the electroweak scale. As long as N 1 is sufficiently light, this reheating temperature is mostly independent of m N 1 and v B /v A . Barring any esoteric model building, the Z 2 symmetry restoration temperature of a given Twin Higgs model should be far higher than such temperatures and domain walls should not be a problem. We also mention that these temperatures are considerably above the lower bound for reheating from BBN [54].

Mirror neutrons as dark matter
As explained before, dark matter cannot realistically take the form of dark atoms in the Mirror Twin Higgs. In this section, we discuss a model without explicit Z 2 breaking in which dark matter consists of mirror neutrons. The model is summarized and the constraints discussed. Two alternative models are then presented.

Model summary
One easy albeit not necessarily obvious way to make the mirror proton heavier than the mirror neutron is via the inclusion of vector quarks. Introduce the vector fermions The part of the Lagrangian that controls the up-type quark masses is whereH M = iσ 2 (H M ) * and where we considered only the first generation. The mass of the lightest eigenstateû M 1 of sector M is then given approximately by A similar mixing could take place for the down quark, but we will assume it to be negligible. The main point of this mechanism is the correction to mûM 1 that goes as (v M ) 3 . Assuming this term is negligible for the down quark, the presence of this correction ensures that the mass of the mirror up quark increases more rapidly than the mirror down quark as v B increases. Barring any experimental constraints, this is sufficient to make the mirror up quark heavier than the mirror down quark and results in the mirror proton being heavier than the mirror neutron, which decreases the abundance of mirror atoms.
As a more technical aside, the presence of the (v M ) 3 term is especially interesting. Ref. [22] studied the mirror neutron as a dark matter candidate in the Mirror Twin Two Higgs Doublet Model (MT2HDM). One of the major challenges was that increasing the mirror vevs could certainly increase the splitting between the mass of the mirror proton and mirror neutron, but it also reduces the mirror Fermi constant. This reduction has the effect of making processes that convert mirror protons to mirror neutrons freeze-out earlier.
These two effects partially cancel each other, either requiring tan β A to be small or forcing certain parameters to be closer to their experimental limits. This dependence on (v M ) 3 ensures that this cancellation is much weaker, avoiding the main issue of the MT2HDM.

Constraints
N eff , Higgs signal strengths and fraction of dark atoms The constraints on N eff and the Higgs signal strengths are applied as in Sec. 2.3. The fraction of dark atoms is computed using the procedure of Appendix F.

Precision measurements
The mixing of chiral quarks with vector quarks also affects electroweak precision measurements. The masses of the up-like quarks coming from Eq. (3.2) are is the positively (negatively) charged part of Q M . This can be diagonalized by performing the basis change whereû M i are the mass eigenstates of sector M ordered from lightest to heaviest. We consider three types of precision measurements.
First, the S and T parameters can be computed using the results of Refs. [55][56][57][58] (see also Refs. [59,60] for their use in relation to vector fermions) where N c is the number of colours, s W (c W ) is the sin (cos) of the weak mixing angle, The new physics contributions to the oblique parameters are then Second, the weak nuclear charges of 133 Cs and 204 Tl from atomic parity violation are computed using the results of Ref. [61]. This gives a contribution from new physics of where Z and N are respectively the number of protons and neutrons in an element. Third, mixing of the chiral up with vector quarks leads to violation of the unitarity of the CKM matrix. The first row is the most precisely measured and the sum of the absolute values of its elements squared becomes where

Direct collider searches
Searches for vector partners of the light quarks have been performed by ATLAS and CMS in the 8 TeV, 20 fb −1 dataset [62,63] and excluded pair production of such quarks up to masses of 845 GeV or lower, depending on branching fractions. Even though dedicated searches for such signatures have not yet been done on the full currently available dataset, it is reasonable to assume that vector quarks with masses ∼ 1.5 TeV and higher are still consistent with the data, given that recent dedicated searches for vector partners of the heavy quarks, whose decays include b jets (which is an easier signature) set limits only up to ∼ 1.6 TeV [64,65]. Additionally, scenarios with large mixing have significant cross sections for single and pair production of vector quarks via electroweak processes, which can be constrained by various LHC measurements [66].  Fig. 6 shows contours of X DA for other values of r T . As can be seen, the fraction of dark atoms can easily be brought to extremely low levels. This can be done while leading to contributions to experimental measurements well below any current limits. In addition, mixing of the chiral up quark with vector quarks could in principle contribute enough to its mass that it might require some amount of tuning for it to remain light. As such, we can define a measure of tuning as

Parameter space and comments
where The tuning is then given by t = 1/∆. As can be seen in Fig. 5e, all constraints can be satisfied without t needing to be small. Do note that a sufficiently large splitting between the masses of the mirror down and mirror up could eventually lead to the spin-3/2 baryon d B d B d B being lighter than the mirror neutron. This would reintroduce the dark atoms problem. A naive estimate in combination with the lattice results of Ref. [67] reveals that this takes place at values of Y U = Y Q and v B /v A much larger than those required to obtain a sufficiently low X DA .
In simple terms, the mechanism works so well because the mass of the up quark is so small that it can be considerably modified without introducing much mixing.

Alternative models
In this section, we describe two alternative models to obtain mirror neutrons as dark matter candidates. We only present the models and leave detailed studies of their constraints for future work.
Both models are inspired by Ref. [22], which showed that mirror atoms could be brought to acceptable abundances in the MT2HDM with explicit Z 2 breaking. The idea of the paper was to introduce two Higgs doublets H A 1 and H A 2 and their partners H B 1 and H B 2 . By assumption, H A 2 provides mass to the up-type quarks and H A 1 to the down-type quarks. The masses of the quarks of the A and B sectors obey the following relation where tan β M = H M 2 / H M 1 . As such, taking a sufficiently large tan β B / tan β A leads to a mirror proton heavier than the mirror neutron and should decrease the abundance of dark atoms. The main challenge however is that increasing the mirror vevs decreases the mirror Fermi constant and makes the processes that convert mirror protons to mirror neutrons freeze-out earlier. It is then necessary to go to relatively low tan β A or be willing to accept a mass of the up quark closer to its experimental upper limit. In the end, the model is compatible with current bounds and does not require additional tuning besides the one necessary to pass the Higgs signal strengths requirements.
The correct structure of vevs was obtained in Ref. [22] by including soft masses that explicitly broke the Z 2 symmetry. The idea of the models of this section is to obtain a similar vevs structure without any explicit Z 2 breaking. In the first model, a pair of new real scalars S A and S B are introduced. The following potential can then be introduced (3.16) Assuming α > 0, the Z 2 symmetry will be broken spontaneously by the first line. At tree level, only one of S A or S B will get a vev and we can assume it to be S B . The second line of Eq. (3.16) then effectively acts as soft Z 2 breaking masses that can be adjusted to reproduce the results of the MT2HDM with explicit Z 2 breaking. The second model is based on Ref. [3]. The following potential is introduced First, assume B µ is zero. If α i is positive and the other negative, the i Higgs spontaneously breaks the Z 2 symmetry by obtaining a vev in only one sector, which can be taken to be the B sector. The other Higgs obtains a vev that maintains the Z 2 symmetry. Once the B µ term is turned on, the Z 2 breaking is transmitted from the broken to the unbroken Higgs sector. It was shown that such a vev structure can pass the Higgs signal strengths requirements. There are then two standard behaviors: (1) α 1 > 0 and α 2 < 0 : tan β A > 1, tan β B < 1, (2) α 1 < 0 and α 2 > 0 : The first possibility is the exact opposite of what is required. The second possibility however leads to a tan β B / tan β A that can be considerably larger than one and at the same time a low tan β A . All the tools necessary to obtain a sufficiently low abundance of dark atoms are then present. The main drawback is that the model leads to a low tan β A , which can complicate UV completions.

Conclusion
The Twin Higgs attempts to solve the little hierarchy problem by introducing a mirror copy of the Standard Model related by a Z 2 symmetry. Because of the measurements of the Higgs signal strengths and cosmology, the Z 2 symmetry must however be broken. The possibility of only breaking this symmetry spontaneously is certainly aesthetically appealing. It was already demonstrated that this can be done for the Higgs signal strengths, but it remained an open question as to whether this could be done for cosmology. As such, the goal of this paper was to determine whether it is possible to create Twin Higgs models in which the Z 2 symmetry is only broken spontaneously that can successfully lead to baryogenesis, provide the correct dark matter abundance, solve the N eff problem and provide a viable dark matter candidate. We found that it is indeed possible to create models that address the above issues. To demonstrate this, we built and studied two of them. In the first model, a pair of Majorana fermions is introduced. In the early Universe, the heaviest Majorana fermion dominates the energy abundance. It then decays, producing a net amount of baryons and mirror baryons as well as some amount of the lighter Majorana fermion. The latter eventually comes to dominate the energy abundance of the Universe. Because of the masses of the particles involved, the lighter Majorana fermion then decays almost exclusively to the Standard Model sector thus reheating it. This model can provide the correct matter abundance without reaching temperatures that would reintroduce domain walls. It can also either solve the N eff problem and generate an acceptably low abundance of dark atoms or provide the correct dark matter abundance.
The second model attempts to convert the dark matter to a form compatible with limits on dark matter self-interactions. This is done by introducing vector quarks that mix with the up quark of their respective sector via Yukawa interactions involving the Higgs. This contributes to the mass of the up quark of a given sector a term proportional to the vev of the Higgs of that sector cubed. This can easily make the mirror up heavier than the mirror down and thus result in a mirror proton heavier than the mirror neutron. The dark matter then takes the form of mirror neutrons and the amount of mirror atoms can be brought to negligible levels. All considered experimental constraints can easily be satisfied and the model can be combined with the first one without adverse side effects.
As a closing word, the models presented in this paper indeed show that the challenges associated to the cosmology of the Twin Higgs without explicit Z 2 breaking can be solved individually and sometimes multiple at a time. However, whether there exists a simple model that can solve all of these problems at the same time is still an open question.

A Thermal averages and energy exchange rates
In this appendix, we present the computations for the thermally averaged cross sections and energy exchange rates. This is done by expanding the work of Ref. [68], from which we reuse the notation.

A.1 General approach to energy exchange
We focus on 2 → 2 processes of the form where i, j, m and n are a set of particles not necessarily of distinct species. We will refer to a generic particle from this set by a lower case Greek letter. The mass and number of internal degrees of freedom of particle α are labelled respectively as m α and g α . In a fixed 'laboratory' frame, the momentum of particle α is labelled as p α , its energy as E α and its three-momentum as p α . A convenient and complete basis for these energies is All thermal averages we will be concerned with are of the form where T α is the temperature of particles α, F (s, E + , E − ) is a generic function, s = (p i +p j ) 2 , v ij is the Møller velocity given by and f α the Maxwell-Boltzmann distribution for particle α Since v ij can be expressed as a function of E + , E − and s, the derived result will still be generic. The inclusion of the v ij factor is simply more convenient. The denominator is trivially given by where n eq α is the equilibrium number density of particle α at temperature T α where K n is the modified Bessel function of the second kind of order n. To simplify the treatment of the numerator, introduce the notation This can be used to rewrite f i f j in the more convenient form Considering that the only non-trivial angular dependence of the differential element is on the angle between the momenta of particles i and j, it can be rewritten as where the equality is as far as integration is concerned. In terms of these variables, the region of integration is given by where p ij is the norm of the center-of-mass (CM) three-dimensional momentum of particle i or j and is given by With this change of variables, the numerator becomes (A.14) Finally, the thermal average is given by .

A.2 E − computation
Consider a given collision ij → mn. In addition to the 'laboratory' frame, one can define a center-of-mass frame. Its three-velocity with respect to the 'laboratory' is labelled as v CM and has norm v CM . Conversely, the three-velocity of the 'laboratory' in the CM frame is labelled v lab and has norm v lab = v CM . Quantities in the CM frame are labelled with a CM subscript. As long as the coordinate systems are properly aligned, the following holds 16) where p ± = p i ± p j and The quantity E − is then related to its CM value by a simple Lorentz transformation where p ± = p m ± p n . The first term of Eq. (A.18) is easily evaluated in terms of standard 2 → 2 kinematics and gives The second term can be evaluated as follows. First, decompose (p m ) CM as where p mn = |(p m ) CM | and θ m is the angle between (p m ) CM and (p i ) CM . The three-vector (p ⊥ m ) CM is the component of (p m ) CM perpedicular to (p i ) CM . In E − , it leads to a term proportional to cos of an azimuthal angle. In the thermal averages, this term vanishes once integrated over that angle as long as axial symmetry is respected. We will ignore (p ⊥ m ) CM from now on. Then, we have v lab · (p − ) CM = 2v CM p mn cos θ lab cos θ m , where θ lab is the angle between v lab and (p i ) CM and we used the fact that (p m ) CM = −(p n ) CM . The quantity cos θ lab is given by The three-vector (p − ) CM is related to its 'laboratory' value by a Lorentz transformation taking the dot product of Eq. (A.23) and p + leads to Assembling everything finally leads to the main result of this section 4

A.3 Results for thermal averages and energy exchange rates
With the results of the previous two sections, it is a trivial matter to obtain the thermally averaged cross sections and energy exchange rates. It suffices to use Eq. (A.15) and then perform the integral over E − , which can easily be done analytically. The results are where the indices on v ij are now implicit, s min = max{(m i + m j ) 2 , (m m + m n ) 2 }, where t is the standard Mandelstam variable and With these results, the computation of any exchange rate is trivial. 5 In a given 2 → 2 process, it suffices to use Eqs. (A.27) and Eq. (A.2) to know exactly the rate at which a specific incoming or outgoing particle gains or loses energy. Knowing the rate at which particles of a given type either gain or lose energy is then trivial. In the limit of T i = T j = T , 5 It is of course understood that σvE − the results of Eq. (A.27) reduce to

B Decay asymmetry
In this section, we compute the asymmetry between the decay of N 2 to baryons and antibaryons. Similar albeit partial results can be found in Refs. [69][70][71]. In an effort to make the result applicable to more generic models, we do the computation for the toy Lagrangian The down-type quarks d and d are assumed distinct. Two assumptions are made: d and d are massless and φ is heavy. These assumptions are made to simplify the calculations, but are not crucial to the mechanism. The mass of u is labelled m u and is not neglected. The fermion N 2 is assumed heavier than both N 1 and u. The leading order diagram for the decay of N 2 to three quarks is shown in Fig. 7a and the next-to-leading order diagram in Fig. 7b. The interference term of these diagrams leads to an asymmetry in the decay to baryons and antibaryons ∆Γ N 2 where x and y are integration variables. Then, also definẽ Figure 7: (a) Tree-level decay of N 2 to three quarks. (b) NLO correction to that decay.
Finally, define When B > A, these functions are given by (B.5) When A > B, they are instead given by (B.6) -26 -With all this, we get the asymmetry (B.7)

C Scattering asymmetries
In a similar fashion to decays, scattering processes of the form N iq → qq and N i q →qq can present an asymmetry in their cross sections. This is due to variations of the diagrams of Fig. 7. We maintain the notation and assumptions of Sec. B. There are then two possibilities. First, there is the asymmetry Obviously, only terms where i = j contribute. In practice, Eq. (C.1) means that N 2 can always present an asymmetry in this scattering, but N 1 can only for a sufficiently large center-of-mass energy. Second, there is also the asymmetry More concretely, these functions are given by (C.5) Because of the kinematics, only N 2 can present an asymmetry and only for N 1 lighter than u and sufficiently low center-of-mass energy. The channel N iū → dd does not present an asymmetry at this order of perturbation.

D Evolution equations
In this section, we present the evolution equations that are used to compute the relic densities in Sec. 2. The relevant processes are first introduced and some important properties are then discussed. To simplify the treatment, we will work with the Lagrangian D.1 Decay: The decay width is The average energy fraction of N 1 in the centre-of-mass frame is obtained by computing the expectation value of m 2 23 = (p d + pd) 2 , which gives It is then a basic exercise in kinematics to compute the average energy fraction of N 1 in the centre-of-mass frame, which gives The decay width is In some regions of parameter space, the decay of N 1 to three quarks is forbidden. The particle N 1 is then forced to go through a four-body decay where u M is off-shell. This decay width is computed numerically, including the width of u M and the mass of d M . The numerical result is also used for m N 1 not too far removed from m u M , as the narrow width approximation is not necessarily a good approximation when these two masses are close.
The cross section is The cross section is The cross section is The cross section is The cross section is The cross section between an A sector fermion and a B sector fermion that is not necessarily its partner is 19) and (D.20) The cross section is where N c p is the number of colours of particle p and

D.10 Combined evolution equations
We now combine the different results together. The number density of particle p is labelled as n p , its energy density as ρ p and its pressure as P p . 6 The differences between the baryon number densities are labelled as ∆B M . The equations are A few comments are in order: These appear in the equations because the baryon asymmetries redistribute themselves amongst the different quarks of a given sector.
• We defined Γ T via 24) In practice, this is simply the decay rate corrected by the fact that boosted particles decay more slowly.
• In equilibrium at temperature T , massive particles respect These equations do not hold out-of-equilibrium, but ratios such as ρ/n remain unchanged. For massive particles, this means that we can obtain their temperature from the ratio of their energy and number densities. At leading order, we have This is the result expected from the equipartition theorem.
• Assume a process ij → mn and its inverse mn → ij. At a given temperature T , the definition of thermal equilibrium implies σv T,T ij→mn n eq i (T )n eq j (T ) = σv T,T mn→ij n eq m (T )n eq n (T ). (D. 27) This equation was used to simplify the evolution equations. Generalization to energy transfers and decays is trivial.
• The Hubble constant is given by where ρ tot = ρ N 1 + ρ N 2 + ρ A + ρ B is the total energy density.
• The quantity g * M corresponds to the effective number of relativistic degrees of freedom in sector M and is computed following standard procedure. A prime represents a derivative with respect to T M .
• The QCD and mirror QCD coupling constants are assumed to unify at high enough scale and are run at one loop order.
• In the dρ M /dt equations, the first term is simply 3H(ρ M + P M ), where we took into account that P M is not exactly ρ M /3 when particles are not fully relativistic. The correction factor takes values in the range [3/4, 1]. This is interesting as, when g * goes to infinity, the energy density scales as the number density. From Eq. (D.26), this means that the temperature remains constant during that time. This is why a loss of degrees of freedom in one sector results in the temperature of that sector rising with respect to the other sector.
• The inverse decay of N 2 to N 1 and quarks or mirror quarks is neglected. Its treatment is complicated, but Boltzmann suppression ensures that it is negligible.
• In some regions of parameter space, the decay u M → N id Md M is allowed. In the B sector, the large mass of the mirror top and Boltzmann suppression render this effect negligible. In the A sector, the decay width to this channel is typically sufficiently small that it would only come into play once the top density is negligible. As such, this effect is neglected.
• At sufficiently high temperatures, certain processes like Z A Z B → Z A Z B can contribute significantly to energy exchange between the A and B sectors. However • The evolution equations are not very stiff in the regions of parameter space studied in this paper. The only exception is for annihilation/scattering between fermions of both sectors via Higgs exchange. At very high temperatures, this process can take place at a rate too high to easily manage numerically. Thankfully, this also means that the two sectors have extremely close temperatures. As such, the problem can be circumvented by treating the A and B sectors as a single population. This is done when the rate at which a sector can exchange energy with the other sector via Higgs exchange is much larger than the rate it receives energy from other sources.

E Higgs signal strengths constraints
In this section, we discuss how the bounds on the Higgs couplings are applied. We follow the procedure of Ref. [22] which is based on the κ formalism [72], albeit the present situation is considerably simpler. Assume a production mechanism i with cross section σ i or decay process i with width Γ i . The parameter κ i is defined such that where σ SM i and Γ SM i are the corresponding SM quantities. For the Mirror Twin Higgs, all κ i are equal at leading order and given by In addition, the Higgs can also decay to mirror particles that escape the detector unseen.
The decay width to a pair of mirror fermions f B is where N c is the number of mirror colours of f B . The decay width to a pair of mirror gluons is given by and α B S is the mirror strong coupling constant. Because of the constraints on v B /v A , decays to mirror massive gauge bosons require both gauge bosons to be off-shell and can therefore be neglected. All other decays to mirror particles are also negligible.
With the above results, limits on the ratio v B /v A can be obtained using the searches of Ref. [73] by ATLAS and Ref. [74] by CMS. These are the most up-to-date available global fits of the Higgs signal strengths and provide all the information necessary to perform a χ 2 fit within the κ formalism. The branching ratio to invisible is not directly constrained by these searches, but indirectly via the reduction of the signal strengths of visible channels. In the allowed range of v B /v A , this branching ratio is far below current constraints (see for example Refs. [75,76]), a fact that was already noted for the Twin MSSM in Ref. [77]. A simple χ 2 fit is performed combining the results of the two experiments and assuming no correlations between them. The results are shown in Fig. 8

F Computation of the mirror atom abundance
In this section, we explain how the mirror atom abundance is computed under the assumption that the mirror proton is heavier than the mirror neutron. We follow the procedure of Ref. [22], from which we summarize the most important elements and to which we refer for more details. Three quantities first need to be computed.
• The mirror QCD scale is computed by requesting that the strong coupling constants of both sectors unify at high enough scale using one loop beta functions.  [73] and [74]. The horizontal lines correspond to the 95% and 99% confidence level limits on χ 2 .
• The binding energy of mirror deuteron B D B is computed via with B 1 ≈ 0.033 and B 2 ≈ −0.011. This equation is obtained by a fit of the lattice QCD results of Refs. [78][79][80][81][82][83] which compute the binding energy of deuteron for different pion masses and the appropriate rescaling. Note that these results contain large uncertainties and as such we will limit all computations to simple approximations.
• The difference between the masses of the mirror proton and neutron is given by where α EM is the fine structure constant, C 1 ≈ 0.86, C 2 ≈ 0.54 and C 0 is fixed to reproduce the equivalent SM value of m A pn . This result also comes from lattice QCD and is obtained from Fig. 3 of Ref. [84] or alternatively Table 2.
With these three quantities, the ratio of mirror proton and mirror neutron abundances can be computed following Refs. [14,85,86]. At high temperatures, collisions with electrons and neutrinos maintain the protons and neutrons in equilibrium. This proceeds at a rate The mirror deuterium bottleneck is crossed when the mirror sector reaches the tem- is 0.08 MeV [86]. This occurs at t DB B ≈ 0.301g −1/2 m Pl r 2 T /(T DB B B ) 2 , meaning that n p B /n n B further decreased by a factor of f 2 ≈ exp(−Γ p B n B e B ν B e t DB B ). At the onset of deuterium formation, the proton to neutron abundance ratio is then (n p B /n n B ) DB B ≈ f 1 × f 2 . Almost all mirror protons are quickly absorbed into Helium-4 nuclei. If the splitting between the mirror up and mirror down is not too extreme, Helium-4 can safely be assumed to be stable and almost all neutrons end up in this isotope. If the splitting is very large, the abundance of mirror protons is bound to be far below experimental constraints and its exact abundance is irrelevant to us. As such, we perform the computation assuming mirror Helium-4 to be stable and mention that the results might not be accurate for extremely low abundances. The final fraction of dark atoms X DA is then