Composite 2HDM with singlets: a viable dark matter scenario

We study the non-minimal composite Higgs model with global symmetry $\text{SO}(7)$ broken to $\text{SO}(5)\times \text{SO}(2)$. The model results in a composite Two-Higgs doublet model (2HDM) equipped with two extra singlets, the lightest of which can be a viable dark matter candidate. The model is able to reproduce the correct dark matter relic density both via the usual thermal freeze-out and through late time decay of the heavier singlet. In the case of thermal freeze-out, it is possible to evade current experimental constraints even with the minimum fine tuning allowed by electroweak precision tests.


Introduction
The spectacular success of the Standard Model (SM) in describing particle physics phenomena, culminated with the discovery of the Higgs boson [1,2], still leaves us with several open problems. Two of the most pressing questions to address are: what is the dynamics protecting the electroweak (EW) scale from large ultraviolet (UV) corrections? What is the nature of the dark matter (DM) of the universe?
The Composite Higgs (CH) paradigm [3][4][5][6] provides a very appealing framework to answer both questions at once, see e.g. the reviews [7][8][9]. In CH models, a new strongly-coupled sector symmetric under a global symmetry G is assumed to exist above the electroweak scale. The Higgs boson arises as a pseudo-Nambu-Goldstone boson (pNGB) of the spontaneous breaking G → H at a scale f . The Higgs boson mass is naturally light, originating from explicit breaking of G, in a similar fashion to what happens for the QCD pion. In most Composite Higgs constructions other particles in addition to the Higgs doublet arise from the symmetry breaking, depending on the specific breaking pattern. In some cases, one of these may be stable and thus possibly have the right properties to account for the non-baryonic DM component of the universe. Such a scenario is one of the few in which the DM candidate is naturally at the same mass scale as the Higgs boson since the DM explanation is tightly linked with the solution of the hierarchy problem (the other notable example being neutralinos in TeV-scale supersymmetry).
We construct and study a CH model based on SO(7) → SO(5) × SO (2). The pNGB field content of this theory consists of two Higgs doublets and two real scalars, the lightest of which is stable and is our dark matter candidate. As will be described in more detail below and in the rest of the paper, the presence of the second doublet will be important to relax the constraints from electroweak precision tests (EWPTs) and DD, while the second singlet can offer an interesting alternative mechanism for DM production in the early Universe. Another desirable feature of this coset is the absence of a Wess-Zumino-Witten anomaly. This model is an extension of the composite Two-Higgs doublet model (2HDM) based on the coset SO(6) → SO(4) × SO (2), studied in detail in refs. [12,[26][27][28]. Several features of our setup are shared with this smaller coset, which however does not include the two singlets pNGB, the lightest of which is our DM candidate.
In this paper, we build the low-energy effective theory of the pNGBs using the tools of Callan-Coleman-Wess-Zumino (CCWZ) construction [29,30] and naive dimensional analysis (NDA) [31,32], paying particular attention to the spectrum and interactions of the pNGBs.
A viable DM candidate needs to satisfy a variety of phenomenological constraints. Most notably, it should reproduce the correct relic abundance and it should not be excluded by direct detection experiments. This implies severe constraints on the model parameters and, in the context of CH models, typically requires f much larger than the electroweak scale, and then relatively large fine tuning. Instead, we find that a viable DM candidate consistent with all phenomenological constraints can be achieved in this model, without paying the price of an excessive fine tuning on the symmetry breaking scale f with respect to the one dictated by EWPTs. This is largely due to the contribution of the second Higgs doublet, which helps both to partially compensate the SM Higgs contribution to direct detection of DM, as well as to relax the EWPT constraints.
Furthermore, we find that DM production may be different from the usual thermal freeze-out mechanism, and proceed non-thermally through decays of the heavier singlet pNGB. This feature is possible because of the richer structure of the model, and (to the best of our knowledge) it is novel in DM models within the CH paradigm.
The present paper is organized as follows: in section 2, we introduce and construct the effective Lagrangian of the model, while in section 3 we describe the spectrum and the interactions of the pNGBs; in section 4, we study the thermal DM candidate of the model, how it can achieve the correct relic abundance, and the corresponding constraints from LHC searches, direct detection (DD) and indirect detection (ID); in section 5, the non-thermal production of DM by decays of the heavier singlet is studied; finally, we conclude in section 6. The appendices contain technical supplementary material, such as the group generators (appendix A), the constraints from EWPTs (appendix B), the detailed expressions of the effective couplings of the NGBs interactions (appendix C), and the calculation of the DM relic density (appendix D).

Effective Lagrangian construction
In this section, we construct the low-energy effective theory, valid below the compositeness energy scale Λ, of pNGBs based on the coset SO(7)/SO(5) × SO(2). In the following three subsections, we present the pNGB fields of this theory, the details of the partial compositeness mechanism employed and the radiatively generated pNGB potential.

Coset and the pseudo-NGB
We consider a new strongly coupled sector lying at an energy scale Λ = m * ∼ (few) TeV and assume that it respects a global symmetry G = SO (7), spontaneously broken to a subgroup H = SO(5) × SO(2) at a scale f ∼ m * /g * by a condensate of the strong dynamics, where by g * we indicate a typical strong coupling of the composite sector. This spontaneous symmetry-breaking pattern produces a set of ten NGBs transforming as a (5, 2) of H.
The global symmetry G of the strong sector is explicitly broken by the SM gauging and the interactions which generate the Yukawa couplings. This breaking induces a potential for the pNGB. The EW gauge group G EW = SU(2) L × U(1) Y is assumed to be embedded in a subgroup H , which in general does not coincide with the subgroup individuated by the vacuum of the theory, H. This well known mechanism of vacuum misalignment between these two subgroups of G, shown schematically in fig. 1, is responsible for the spontaneous breaking of G EW .
In general, we can then consider two basis of generators: one, {T θ }, related to the breaking G → H, and a second one, {T }, related to G → H . The groups H and H are misaligned by an angle θ, which in general is a vector, since more than one field can acquire a vacuum expectation value (VEV). The corresponding vacua are related by a rotation matrix r θ and, if we assume that the generators are normalized as Tr{T A T B } = δ AB , we have: We then introduce the Goldstone matrices in the two basis as: where we defined Π ≡ Π IT I and Π θ ≡ Π IT I θ , withT I ,T I θ being the broken generators in the two basis and Π I the respective pNGB fields.
The rotation r θ can be obtained by considering the Goldstone matrix in the gauge (non-rotated) basis {T }, and setting the NGBs at the corresponding VEVs, i.e.: Our choice for the generator basis is described as follows: H H Figure 1. The effect of gauging H ⊂ G is a breaking of the global symmetry due to vacuum misalignment.
The generators of SO(5) ×SO (2) , whose expressions can be found in appendix A, are then block-diagonal in the basis we adopted. We indicate with the first letters of the alphabet the indices of a generic SO(7) transformation, a, b, c, · · · = 1, . . . , 7.
Because of its block-diagonal form, instead, an SO(5) ×SO(2) transformation will haveā,b, . . . indices, whereā = {i, µ}, with i and µ being SO(5) five-plet and SO(2) doublet indices, respectively. In the following, we do not distinguish between upper or lower indices, identifying however the first and second indices as row and column ones, respectively. As already stated, the generators in the physical vacuum (G/H) basis are related to these by T θ = r θ T r −1 θ ; the virtue of this approach is that we can expand U θ in the fields to extract the interactions, while keeping the exact, trigonometric, expression for the parameters related to the vacuum.
In order to identify the accidental symmetries of the theory and the quantum numbers of the pNGBs, it is useful to start by considering the limit of no misalignment, i.e. r θ = 1. In this limit, and with the generator basis specified above, the NGB matrix takes the simple form: The two NGB five-plets of SO(5) can be decomposed under representations of the custodial symmetry SO(4) c ⊂ SO(5) as (5, 2) = 2 × 4 + 2 × 1: The two 4 describe the two Higgs doublets of the theory: where h is the SM-like Higgs, G i are the would-be longitudinal polarizations of the EW gauge bosons, and H 0 and A 0 are the CP -even and -odd components of the neutral scalar, respectively. The lightest of the two singlets, η, will be the DM candidate.
As discussed in detail in ref. [26], it is useful to introduce a discrete transformation of the pNGBs, namely: . This symmetry is important to protect the second Higgs doublet from taking a sizeable VEV, thereby spontaneously violating custodial symmetry beyond the allowed limit. As we will see, the interactions of the strong sector with SM fermions in general break this symmetry. Another useful parity is: under which (η, κ) → −(η, κ), which can thus stabilize the singlets. In the following, we impose P 7 as a symmetry of the theory. An appealing property of the coset under consideration is that no Wess-Zumino-Witten term [33,34] is generated since the fifth de Rham cohomology group of SO(7)/SO(5) × SO(2) vanishes [35,36] 1 . This implies that the P 7 parity remains unbroken by the strong dynamics to all orders in the chiral expansion. The parities of the NGBs and their representation under SO(4) are summarized in table 1.
Lagrangian terms which break explicitly the global symmetry G could in general also break these parities. In the following we show how, in general, C 2 is indeed broken by the interaction of the top quark with the composite sector and by the potential, while P 7 can remain a good symmetry without being broken at any level.
Compatibly with P 7 and CP , the general misalignment of the vacuum with respect to the gauged subgroup H can be described by two angles θ 1 = h /f and Field SO(4) C 2 P 7 Table 1. Representations and quantum numbers of the NGB fields under SO(4) , C 2 and P 7 . φ 1 can be identified as the SM Higgs doublet if C 2 is a symmetry of the NGB potential or more in general if φ 2 does not take a VEV.

CCWZ Lagrangian
The leading operator describing the low-energy NGB effective theory is, in CCWZ language [29,30]: This Lagrangian contains the kinetic terms for the Nambu-Goldstone fields, mass terms for the SM EW gauge bosons and their interactions with pNGBs, as well as an infinite series of two-derivative interactions among NGBs. The mass of the W boson is given by: It is then convenient to define: so that, defining m 2 W = g 2 v 2 /4, we get: ξ ≡ v 2 f 2 = sin 2 θ 1 + sin 2 θ 2 , tan β = sin θ 2 sin θ 1 , (2.14) with v = 246 GeV being the SM VEV.
On the other hand, the prediction for the Z mass is: leading to a tree-level positive contribution to theT parameter: This is due to the fact that H 0 explicitly breaks SO(4) c to SO (2). The different CP structure with respect to [26] implies that we have custodial breaking even if we assume CP to be preserved. While O(1) values of β are disfavored by EWPTs, as they would require very high fine tuning, values β 0.1 are allowed. Such a contribution might even help to improve the fit to electroweak precision observables (see appendix B for more details). As shown in section 3.1, small values of β are obtained naturally in this model. Interestingly enough, the positive contribution tô T given by eq. (2.16) can help to relax the usual EWPT limits on ξ. We find that for β ≈ 0.1, a fine tuning up to ξ ≈ 0.08 is compatible with both EW and Higgs data (see appendix B).

Partial compositeness
To couple the SM fermions to the Higgs field and generate their masses, we resort to the partial compositeness paradigm: the basic idea is that quarks are linearly coupled to fermionic operators O L,R belonging to the strong sector [37]. In the following, we assume for simplicity that the operators coupled to the top quarks transform in the fundamental representation of SO(7), although other choices are possible (see e.g. [17] for other representations of SO(7)).
As usual in Composite Higgs models, the group G has to be enlarged to correctly reproduce the SM quantum numbers: to this purpose we consider SO(7)×U(1) X , where the charge X is X = 2/3 for the top quark. The hypercharge is then identified with Y = T 3 R + X. The 7 decomposes under SO(5) ×SO(2) and SU(2) L ×U(1) Y as: (2.17) We see from this decomposition that the right-handed quark t R can be coupled to the 1 of SO(5) and the singlet in the 5, while the left-handed doublet q L can only couple with the 5. We consider the following Lagrangian for the top quark: where a = 1, . . . , 7 is an SO(7) index, α = 1, 2 is the flavor index of the quark doublet and Y L,R are the spurions. The SM fermions are assumed to be even under both C 2 and P 7 . We promote the couplings Y L,R to fields (the spurions) whose transformations under G are dictated by the ones of the operators O L,R . Compatibly with P 7 , and rotating away unphysical components with the elementary U(2) el L and U(1) el R symmetries under which the spurions and the quark fields transform, the most general VEVs for the spurions are: with y L and y R real, and the fifth component of Y R set to zero by the P 7 parity. It is evident that the VEV of Y R breaks C 2 unless θ t = 0. It is important to distinguish between two types of symmetries: the spurionic ones are symmetries of the strong sector which are unbroken before the spurions acquire a VEV; the residual ones are symmetries at the electroweak scale, which remain unbroken even after the spurions have acquired a VEV. In the following, we assume that spurionic symmetries are respected by the theory. In order to build the partial compositeness Lagrangian, we "dress" the spurions with the NGB matrix and define: This definition is consistent with the standard one, i.e.Ȳ = U † Y: this can be easily checked by going to the basis of the VEVs, where U θ = 0 and U = r θ . In general, the dressing procedure has the effect to take an object transforming with an index a of G into a new object transforming with an indexā of H. It is then understood that whenever a barred quantity appears, barred indices are implicit.
In order to write the low-energy effective Lagrangian obtained from integrating out the composite sector, we follow the standard procedure, detailed for example in ref. [26]: we use the dressed spurions and SM fields to write operators invariant under H. This will also assure their invariance under the full G. The pNGB dependence will be included in the dressed spurions.
In the aligned limit, θ 1,2 = 0, the dressed spurions transform as a (5, 1) ⊕ (1, 2) under SO(5)×SO(2), with components given by: with i = 1, . . . , 5, while µ = 1, 2 is the index associated to SO (2). The effective Lagrangian can then be constructed by combining them with δ ij , δ µν and µν . The latter possibility, however, violates the C 2 -spurionic and will not be considered. In addition, only left-right combinations have to be considered because of chirality. Finally, the two invariants which can be constructed with δ symbols are not independent, due to the singlet one can obtain by combining two 7.
The leading order operator generating the top mass is thus given by: where the coefficient m * /g 2 * comes from NDA and is such that c t is a coefficient expected to be O(1). From this effective Lagrangian one obtains both the top mass and the top-NGBs interactions. Expanding around β = 0, the top quark Yukawa coupling is given by: Note that the factor in parenthesis approaches β for θ t → π/2. This suppression can be compensated by a slightly larger value of c t or of y L y R /g * . We can consider an analogue Lagrangian for the compositeness of the other quarks, but since Y b Y t , the other contributions are expected to be subleading. In particular, given that the choice has no major effect on the potential, we take θ b = 0.
A crucial point is that eq. (2.22) leads to an interaction of the type ihA 0t γ 5 t after the spurions acquire a VEV: this interaction explicitly violates C 2 (in particular, the one which is broken is the C 2 -residual in the language introduced before), implying that in general also the second doublet takes a VEV.
Furthermore, it also turns out that η and κ have opposite CP -parities: we assume that η is even and κ is odd.
In order to avoid sizable flavor changing neutral currents (FCNCs), the embedding of all the other SM fermions should be fixed carefully. Choosing the fundamental representation for the fermions embedding leads in principle to two independent strong sector invariants; however, the spurionic C 2 forbids one of these. Nevertheless, to avoid FCNCs, different families should have the same embedding, including the same choice of θ u = θ c = θ t , θ d = θ s = θ b . Even so, it is well known that in Composite Higgs models operators from the strong dynamics at the energy scale Λ can induce potentially dangerous flavor violating effects both in the quark and in the lepton sector. However, a more detailed discussion of the flavor phenomenology of this model is beyond the scope of this work.

Pseudo-NGB potential
The pNGB potential is generated from the explicit breaking of the Goldstone symmetry due to the gauging of the EW subgroup of G and to the mixing between SM fermions and the composite sector. Using naive dimensional analysis, the radiativelygenerated potential can be written schematically as (see e.g. ref. [32]): whereV is a dimension-less function of the NGBs, L counts the number of loops at which each term is generated, and µ G and µ F count the required insertions of the gauge and fermionic spurions, respectively. The construction of the different terms in the potential by building invariants from the spurions follows closely the discussion presented in [26]. In the following we shortly describe only the main parts.

Gauge contributions
One source of explicit breaking of the global symmetry of the strong sector are interactions between the SM gauge bosons and the pNGBs. It is convenient to introduce a set of spurions: transforming under g ∈ G as G X → gG X g † . They can be dressed with NGBs as: Their componentsḠ X A = Tr Ḡ X T A transform as the following multiplets of SO(5) × SO(2): 27) associated to {T L , T R , T 5 }, {T 1 ,T 2 } and T 2 , respectively. We can thus organize the components ofḠ X as: with I = 1, . . . , 10 being an index in the adjoint, while i = 1, . . . , 5 and µ = 1, 2 being the indices associated to SO(5) and SO(2), respectively. The set of independent invariants with two spurion insertions, compatible with C 2 and P 7 , is: The gauge contribution to the NGB potential is then given by: where c (i) g,g are O(1) coefficients.

Fermionic contribution
The main source of explicit breaking of the Goldstone symmetry is due to the coupling of the composite sector with elementary quarks, and in particular with the top. The relevant Lagrangian was already introduced in eq. (2.18). The first step to build the possible invariants which can enter in the potential is to construct combinations of spurions which are invariant under the elementary gauge symmetry of eq. (2.18):∆āb The independent set of invariants which can be obtained at O(y 2 ) are: At O(y 4 ), the non-vanishing invariants are: where the indices have to be interpreted as already indicated. While the operators indicated with (1) are generated at one loop, all the other ones are generated at two loops [26], and are thus accompanied by a further factor of g 2 * /(4π) 2 . The general form of the scalar potential was given in eq. (2.24); for the fermionic case, it can be expressed as: where I (i) (n L ,n R ) is an invariant formed with n L,R powers of∆ L,R , and c (i) (n L ,n R ) are O(1) coefficients. Since the fermions in the loop generating the potential are colored, there is a factor N c accounting for the number of colors; in the following, we take which is the reason for the denominator in the previous estimate for the potential. Since we assumed there is no further CP breaking coming from the effective Lagrangian, we set c (0,2) to 0, since the associated invariant contains CP breaking terms.
It turns out that c (1,1) and c (2) (0,2) are the most relevant coefficients for numerical estimates. In our numerical scans, we take for simplicity all the other coefficients (namely c (2) (2,0) and the ones coming from the gauge invariants) equal to 1, since they do not play a relevant role. We generically denote by c i the O(1) coefficients, and define three possible ranges of variation of these coefficients, depending on how close they are to unity: • strictly natural coefficients: 0.2 ≤ |c i | ≤ 5; • loosely natural coefficients: 0.1 ≤ |c i | ≤ 10; • unnatural coefficients: |c i | < 0.1 or |c i | > 10.

NGB dynamics
In this section, we summarize the main properties of the pNBGs, such as their vacuum structure, their spectrum and interactions.

Vacuum structure
By setting to zero the pNGB fields in the potential V tot = V gauge + V fermion , we can find the minimum for the misalignment angles θ 1 and θ 2 . In practice, we impose that the minimum is found for the required benchmark values of ξ and β by solving for two of the free coefficients. Specifically, we solve for c (1) (1,0) and c (1) (1,1) and check that the solution lies within the desired naturalness range. An approximate expression for ξ, obtained at leading order in y L,R /g * 1 is: A tuning among the coefficients in the numerator must be imposed in order to reproduce the desired misalignment, the amount of which is of order ∆ ∼ ξ −1 . As already discussed in ref. [26], a hierarchy θ 2 θ 1 , i.e. β 1, is instead naturally obtained in this model. By minimizing the potential, we get approximately: .

(3.2)
A strong suppression is automatically obtained for θ t ∼ π/2. Furthermore, for g y R cos 2θ t , one can approximate tan β ∼ y 2 L /g 2 * tan 2θ t which shows clearly a suppression if y L g * . Another interesting region we will study is close to θ t ≈ π/4. In this case the hypercharge term in the denominator cannot be neglected, but values of tan β ∼ 0.1 are still naturally obtained.

Spectrum
Due to the smallness of β, we can perform a power expansion in the expression for the pNGB masses we obtain from the potential. A mixing between h and H 0 is present in general, and can be diagonalized via a rotation by an angle α ≈ β. Once the conditions fixing ξ and β have been imposed, the physical Higgs mass at leading order in ξ and β is given approximately by: where we omitted contributions from gauge or two-loop coefficients. In the second line we substituted the expression for the top-Yukawa (cf. eq. (2.23)). A small value of g * , i.e. light top partners, helps to avoid a further tuning in order to obtain the correct Higgs mass. For this reason in the numerical analysis we fix g * = 3. In practice, we impose the measured value m h ≈ 125 GeV, by solving the (exact) m h expression for the coefficient c (2,0) . Once ξ, β, and m h have been fixed, the masses of the other pNGBs, H 0 , A 0 , H ± , η and κ as functions of the remaining coefficients, to the leading order in ξ, are: Due to our choice of coefficients the gauge contribution cancels at the first order in ξ, but it is present at the next to leading one. As we can see, H 0 , A 0 and H ± are almost degenerate in mass. We then assume that π/4 ≤ θ t ≤ π/2 and c (1,1) > 0; with this choice, all the mass parameters are positive. The sperum is shown in fig. 2 for ξ = 0.061 and β = 0.1, as well as different values of c (0,2) and c (5) (1,1) . When θ t ≈ π/2, also the first order in ξ can play an important role for m η ; if all the coefficients are O(1), we can approximate it as: (0,2) and c (1,1) in the strictly natural range, |c i | ∈ [0.2, 5].
θ t π/4 and with η and κ very close in mass, of O(1 TeV); in this case, if κ has a long enough lifetime, it can freeze-out in the early universe and then decays into η, which is the stable DM relic, giving rise to non-thermal DM production. These two scenarios are studied in sections 4 and 5, respectively.

pNGB interactions
In order to study the pNGBs interactions, it is convenient to consider separately those coming from CCWZ and the ones coming from partial compositeness and the potential.

Interactions from CCWZ
The CCWZ Lagrangian in eq. (2.11) contains pNGB interactions with the SM EW gauge fields, as well as derivative self-interactions: where L kin contains the pNGB kinetic terms, k der = 2ξ/3, and Φ 1,2 are the two pNGB five-plets of physical fields introduced in eq. (2.6). We neglected all couplings which break custodial symmetry and which become negligible once the limits from EWPTs are taken into account, and omitted other interactions with two gauge bosons and two pNGBs, less relevant for the phenomenology discussed in the following. It is worth noticing that interactions with three pNGBs and two derivatives, such as η 2 h, are absent from the Lagrangian above. This might seem surprising at first, since such interactions have been long known to be present in similar scenarios, and their relevance has been often stressed (see e.g. refs. [10][11][12]17]). Such interactions usually arise from Π 4 terms, once the Higgs(es) takes a VEV. However, since we employ the description of pNGB fields from the misaligned vacuum (as discussed in section 2.1) no field takes a VEV and therefore these terms are not generated. Another way to easily understand their absence is to set to zero the gauge and Yukawa couplings. In this case, the global symmetry is exact and all vacua are degenerate, so that the SO(5) and SO(5) subgroups are physically equivalent. In this limit, the two-derivative NGBs interactions start at O(Π 4 ) in both vacua (and are of the form specified above). Now, switching on the gauge and Yukawa couplings selects SO(5) as the true vacuum; however, since derivative interactions do not depend on these couplings and since in our descriptions fields do not take a VEV, the derivative interactions are not affected and therefore cubic ones are not generated.
The connection with the description most commonly employed in the literature (i.e. describing fields from the gauge vacuum SO(5) and allowing then the Higgs to take a VEV) can be easily obtained via a non-linear field redefinition ( [12,26]). For example, in the limit of θ 2 = 0 this is given by (h andη are the physical fields in the gauge description): Such a transformation generates cubic derivative interactions from the kinetic terms, as well as non-derivative interactions from the pNGBs mass terms. The net effect of these is to keep physical observables invariant under such transformations.

Interactions with fermions and from the potential
Let us now list other phenomenologically relevant pNGB interactions, in particular those with SM fermions and self-interactions from the potential: The expressions of the effective couplings are reported in appendix C. We require them to be always less than 4π for perturbative reasons: this usually forces c (0,1) to be smaller than 1.

Thermal dark matter scenario
The first region of interest is the one for θ t π/2, where η is the lightest pNGB and the dark matter abundance is generated via a thermal freeze-out. The fact that θ t ≈ π/2 has several implications: from fig. 2, we see that η is much lighter than the other resonances, so that the effect of other pNGBs can be largely neglected for the freeze-out computation; a small value of β is natural and does not require a further tuning or unnaturally small value of the potential coefficients, see eq. (3.2); finally, since cos θ t is small, the expression for the top mass, eq. (2.23), has a mild suppression which must be compensated by having either a large c t or y L y R g * ; we choose the second option for naturalness reasons, and take y L = 2, y R = 3 and g * = 3. This choice is consistent with perturbativity and with the phenomenological requirement of having a less composite left-handed top rather than the right-handed one. Other choices are possible and we do not expect them to change the results qualitatively.
We also explored other regions of the parameter space: in particular, the region in the range θ t ≈ π/4 could be potentially interesting because of coannihilations with other pNGBs (cf. fig. 2); however, we verified that in this region it is not possible to reproduce the observed relic density via the standard freeze-out mechanism. Instead, in this region a non-thermal DM production mechanism can take place, whose discussion is deferred to section 5.

Relic density
The main contributions to the relic abundance are given by DM annihilations into SM EW gauge bosons, Higgs and top quark. In our computations, we also included subleading contributions, and all the details can be found in appendix D. The relic density profile as a function of the DM mass is shown in fig. 3. The darker (lighter) region is obtained by letting the coefficients vary inside the strictly (loosely) natural range.
It is useful to remember that the relic density is inversely proportional to the integrated thermally-averaged annihilation cross section: Ωh 2 ∝ ∼ 1/ σv . This is the reason why a plateau appears for m η m h , where the cross section is dominated by annihilations into SM gauge bosons: since the latter do not depend on the c i 's, the annihilation cross section is always bounded from below. The situation changes at larger masses, where new annihilation channels open up: in principle, we expect the relic density to decrease with increasing mass. However, it is possible that a cancellation in the main contributions to the effective cross section occurs: this is precisely the case in the region m η ≈ 400 GeV of fig. 3, where η is sufficiently heavy so that the exchange of H 0 compensates the exchange of h in s-channel, the two contributions having opposite signs.
Of all the effective couplings, the one that plays the most important role is g ηh , describing the interactions between two h's and one η. It enters different processes with different signs, so it is non-trivial to describe its role analytically. The results highly depend on the parameters ξ and β and, guided by EWPTs, we decide to focus our attention on ξ = 0.061 (corresponding to f = 1 TeV) and β = 0.1. From fig. 3, we see that there are two good mass regimes which give the correct relic density, Ωh 2 = 0.1198 [38]: m η ≈ m h /2 and m η ≈ 150 GeV; a third one, at m η ≈ 400 GeV, also reproduces the correct relic density if the coefficients are allowed to vary inside the loosely natural range. However, as we will show, this high mass range is already excluded by direct detection.

LHC searches
In the region of parameter space where m η ≤ 62.5 GeV, the Higgs can decay into two DM particles, h → ηη. Experimental constraints on invisible Higgs decays have been obtained by ATLAS and CMS [39,40] for various possible assumptions on the Higgs couplings to SM particles. Since the couplings of the Higgs to quarks and gauge bosons in our model are different from the SM ones, the widths of the decays into SM particles have to be appropriately rescaled. However, no major departure from the SM result is expected for the values of ξ and β allowed by the electroweak precision test. We thus take as an experimental bound BR inv < 19% at 95% CL [39]. The HL-LHC will reach a 95% CL exclusion sensitivity of 1.9%, while future electron-positron colliders would be able to reach sub-percent precision (see, e.g., ref. [41] for a recent review of Higgs boson measurements at future colliders).
In our model, the invisible Higgs decay width is given by: which depends on the coefficients c i 's via the effective coupling g ηh . The corresponding prediction for the invisible branching ratio of the Higgs, obtained by letting the c i 's vary inside the strictly (loosely) natural range are shown as dark (light) yellow regions in fig. 4. As can be seen, the model's predictions are always lower than the current experimental bound. In the same figure, the small region where also the correct relic density is obtained is shown in blue. It is interesting to notice that HL-LHC will be able to test most of this parameter space via invisible Higgs decays.  The missing energy trace could also be produced by direct production of η particles or by the decay of other massive scalars (for instance H 0 , which is also linearly coupled to η). In this case, we need to look at specific tags: in our model, the most relevant ones will be an energetic jet, monojet signature (MJ), or two well-separated jets, Vector Boson Fusion (VBF) signature. In both cases, it is important to describe the effective coupling of gluons to the massive scalars, eq. (3.9b), since this gives the main contribution to the production cross section. We implemented the model in FeynRules [42,43] and generated simulated events with MadGraph5 [44]. Even tough h is always lighter than H 0 in the parameter space we are considering, diagrams involving H 0 still need to be considered in the MJ and VBF processes, since these can give sizable contributions. We obtain that for masses above 50 GeV both monojet and VBF do not put any constraint on the model, being always at least an order of magnitude below the experimental limits ( [45,46]), for any reasonable values of ξ, β and the c i .
Overall, DM searches at LHC do not put important constraints on the parameter space of our model.

Direct detection
Direct detection experiments usually put strong constraints on DM models, including the one studied here. In our model, the DM candidate η can interact with quarks either via the contact interaction generated by partial compositeness, proportional to the coupling g q (see eq. (3.9a)), or through an exchange of h and H 0 . A convenient way to evaluate the spin-independent DM-nucleon cross section is to parametrize the interaction Lagrangian as: L where: As already stressed, in order to avoid FCNCs, one must have a t = a c = a u and a b = a s = a d . Because of the different signs between the first term and the parenthesis, it is their interplay to determine the allowed parameter space. In particular, since g q is negative (see eq. (C.2a)), negative values of the parenthesis will be favored, leading to a partial cancellation of the two terms in eq. (4.3). Notice that the value of a q actually depends on the coefficients c i , so that a cancellation in the scattering amplitude can occur, allowing to evade the DD constraint. From the effective Lagrangian above we derive the spin-independent DM-nucleon cross section [47]. At present the strongest constrain on it comes from the XENON1T experiment [48]. As can be seen in fig. 5, this casts important constraints on all the three regions in m η where the relic density is realized in our model, in particular excluding completely the large-mass one.
In our plots we also show the prospects for the future exclusion bounds coming from the XENONnT experiment [49]. We observe that the majority of the parameter space, for the values of ξ and β used, will be tested.

Indirect detection
While being somewhat beyond the scope of the paper, we briefly consider the constraints from indirect detection as well. We mainly focus on limits from dwarf spheroidal galaxies (dSphs) given by Fermi-LAT, and reported in ref. [50]; the relevant one for our model is given by DM-annihilation into bb. In the region m η ≈ m h /2, only the direct process ηη → bb has to be considered; for higher DM masses, instead, also the intermediate productions of W, Z, h (and possibly other NGBs) are important: as a conservative estimate, we assume that these intermediate states completely decay into bb. We also took into account possible branching ratios, but the result is only slightly modified by this correction.
As already stated in ref. [17], also anti-protons bounds from AMS-02 are worth exploring; however, since the size of the systematic uncertainties is still unclear, we limit ourselves with constraints from dSphs, leaving a comprehensive treatment of ID in this model to future work.

Discussion
We summarize the main results for f = 1 TeV (ξ = 0.061) and β = 0.1 in fig. 5. As anticipated, three regions are possibly interesting: m η ≈ m h /2, m η ≈ 150 GeV and m η ≈ 400 GeV. The orange (purple) hatched area is excluded by DD (ID), while in the blue one the correct DM abundance at 3σ, Ωh 2 = 0.1198 ± 0.0036 [38], is reproduced. In these plots, for better readability, we only show the regions obtained by letting the coefficients vary within the loosely natural range. In the gray region the coefficients are beyond this range (unnatural). Current indirect detection limits from dSphs are close to parameter space where the correct relic density is reproduced, but unable to exclude any portions of it. A future stronger bound should be able to test this model.
The correct relic density is reproduced for masses just below and just above the on-shell Higgs production threshold of 62.5 GeV. Since this regime is good due to the Higgs resonance, the allowed mass range is very narrow and given that there is no symmetry or dynamical argument to expect such a value for m η , this would represent a further tuning in the model parameters. The experimental results coming from the Higgs invisible BR and ID do not exclude any parts of the region where the correct relic density is obtained, while DD sets an upper bound on g ηh . Even if ID does not exclude any points in the parameter space at the current state, the cross section in our model is only an O(1) factor below the experimental constraint, so that upgraded searches are expected to either find a positive result or exclude this region of the parameter space. It is also worth mentioning that HL-LHC will test most of the region below the Higgs pole via a precise measurement of the invisible branching ratio of the Higgs; however, in the near future DD and ID seem more promising directions.
Intermediate mass range (m η ≈ 150 GeV) At masses larger than around 100 GeV, the most relevant constraints come from direct detection and the requirements of naturalness of the coefficients. Limits from LHC experiments are not relevant for this mass regime and are not expected to be able to probe it in the near future. EWPTs and Higgs couplings constraints are safe because of the values of ξ and β we considered. We are intentionally overestimating the bound from ID, and yet this search is not putting constraints on the model. Upper and lower bounds on m η are set by DD results (cf. fig. 5). If we were to limit to strictly natural coefficients, then masses below 135 GeV would be excluded. The feature at m η ≈ 180 GeV is because of a cancellation in the cross section for ηη → tt, due to the different sign between the effective couplings g ηh and g ηH 0 (cf. eq. (D.1c)).
This seems to be the most promising mass range, because ID is pretty weak and DD leaves a significant region of the parameter space available. Also, such m η values are naturally obtained for θ t π/2. Decreasing the value of ξ has the effect of enlarging the allowed mass range, so that it is interesting to investigate how the latter varies with varying fine tuning: this can be seen in fig. 6. We observe that DM phenomenology in this model allows for low values of f , up to f ≈ 750 GeV for strictly natural coefficients and even below 600 GeV for loosely natural ones. Indeed, the relevant lower bound on f in our setup is the one due to EWPTs and Higgs coupling measurements. This should be compared to other similar non-minimal composite DM models in the literature, for which significantly larger values of f were found to be necessary (see e.g. [12,16]). The reason why we can have viable DM phenomenology with a fine tuning which is lower than in other models is two-fold: on the one hand, the contribution tô T from the second doublet is positive, so that a smaller fine tuning is allowed; in addition, because of the richness of the model, it is possible to have a cancellation in the amplitude relevant for DD, resulting in a viable region of the parameter space compatible with EWPTs.
What happens varying β is less trivial: by increasing β, we increase the effect of subleading terms and the coefficient dependence, and this can in principle affect the results; however, we verified that this is not the case for the values of β allowed by EWPTs.

Large mass range (m η ≈ 400 GeV)
The third region of interest is for m η ≈ 400 GeV. As already stated, its existence is due to a cancellation in the main contributions to the relic density between the terms with the exchanges of h and H 0 . The correct relic abundance can only be reproduced in the loosely natural range of the c i coefficients. Our benchmark point is already excluded by current DD constraints, and in order to evade these limits one would need ξ 0.01. For this reason, we do not study this region further.

Non-thermal dark matter production
With the larger number of pseudo-Goldstones present in this model, another possible production mechanism for the DM candidate η, other than the thermal freeze-out discussed above, could be via the decay of heavier pNGBs.
It turns out that the interesting case is when κ and η are close enough in mass to allow for a sufficiently long decay of κ into η. From fig. 2, we see that this happens at roughly θ t ≈ π/4.
In this scenario, the η relic density receives two contributions. The first is due to the standard thermal DM freeze-out. The second mechanism is due to the freeze-out of κ, which is sufficiently long-lived, followed by its decay into η. Since a single η is produced per each decay, the κ number-abundance is completely converted into η. The total η relic density is thus given by (see e.g. ref. [51]): It is important to notice that since H 0 , A 0 , H + are all lighter than η (cf. fig. 2), they also have to be included as final states of η-and κ-annihilations. On the other hand, as in standard coannihilations, processes which simply convert η in κ (and viceversa), are not relevant for the determination of the relic density. We do not report the formulas for all the annihilation channels, but these can be easily computed from the effective Lagrangian given in eq. (3.6) and eqs. (3.9a) to (3.9c). As shown below, in the parameter space of interest the relevant decay channel is only κ → ηbb. The decay width for the process κ → ηqq is: The κ lifetime has to be compared to the age of the universe at the time of the freeze-out of η, which is given by: with g * ,EU ≈ 100 being the effective number of relativistic species in the early universe and x F ≈ 25. Given that Γ κ→ηqq ∝ m 2 q , if ∆m κ,η ≡ m κ − m η > 2m t , the decay into tt is so quick that κ always decays before η freezes-out. This can be avoided by either having η and κ close in mass or by taking g ηκq very small; we choose the first option as it requires less fine tuning. On the other hand, if ∆m κ,η < 2m t , the two contributions to the relic density are of the same order, i.e. Ω η h 2 ≈ Ω κ h 2 . Compared to the thermal case, a higher fine tuning on ξ is needed in order to evade DD constraints. The phenomenology has substantially changed from the thermal case because, by requiring that the two singlets are close in mass, we have put ourselves in the portion of parameter space near θ t = π/4, where both singlets are very heavy.
Analogously to what we did in fig. 5, we show in fig. 7 the results in the m η -g ηh plane, for ξ = 0.01, β = 0.2, y L = 1 = y R , g * = 3 and loosely natural coefficients. While DD excludes a portion of the parameter space, current bounds from ID are ineffective. In general, a fine tuning on the masses is needed, since in all the plane 20 GeV ∆m κ,η 50 GeV.
As one can see from eq. (3.4d), this range of masses for m η roughly corresponds to θ t ≈ π/4, as anticipated. While large mass splittings tend to favour a fast decay for κ, non-thermal effects are always possible for small mass splittings, although a larger and larger unnaturalness of the coefficients is required (corresponding to a larger and larger fine tuning for ∆m κ,η /m η ).
The non-thermal production mechanism represents an intriguing feature of this model; we think this is one of the most peculiar and interesting aspects of the model. The greater level of complexity with respect to the minimal case has been traded for a richer spectrum of NGBs which can play an active role in DM phenomenology.

Conclusions
In this paper we carried out a detailed construction of the CH model based on the symmetry breaking pattern SO(7) → SO(5) × SO (2). The lightest pNGB, η, is electrically neutral, stable and is a potential DM candidate.
We studied the DM phenomenology by requiring correct relic abundance and evading the constraints coming from invisible Higgs decay, direct detection and indirect detection, and found large portions of parameter space satisfying all of them (see section 4). In particular, we identified a viable region of DM mass around m η ≈ 130 ÷ 160 GeV (see fig. 6), which is realized for a symmetry breaking scale as low as the minimum required by the ElectroWeak Precision Tests (f 0.8 TeV). This is mostly due to a cancellation in the couplings of ηη to qq between Higgs-exchange and H 0 -exchange in the s-channel, which allows to enhance the relic abundance and deplete the direct detection cross section. This feature is peculiar of the low-energy 2HDM-like structure of the model.
Another important aspect of the model is that the extra pNGB κ may freeze-out in the early universe and subsequently decay to η, thus providing an extra (nonthermal) contribution to the DM density (see section 5). To the best of our knowledge, no other model studied so far for DM in the CH framework provides such a possibility.
The main results of this paper may be summarized as follows: • the CH model based on the symmetry breaking pattern SO(7) → SO(5)×SO(2) delivers a viable DM candidate, consistent with the current phenomenological constraints; • the correct amount of DM can be achieved with relatively low amount of fine tuning on the symmetry breaking scale f 0.8 TeV; • it is possible to produce DM also non-thermally via late-time decays of the heavier pNGB.
A more exhaustive study of the indirect detection constraints (e.g. using anti-protons data), of the detection prospects at future colliders, as well as a detailed analysis of the UV completion of the theory (by including, for instance, the top partners) is beyond the scope of the present paper and left for future work.

B Higgs couplings fit and EWPTs
To obtain the constraints on the parameters of our model from Higgs measurements we use the recent global Higgs coupling analysis published by ATLAS with 80 fb −1 of luminosity [52]. Restricting the fit to the deviations relevant in the model one has: For the EWPT, we use the updated combined limits on S = 4s 2 WŜ /α and T = T /α from GFitter [53]. The dependence of theŜ andT parameters [54] on the model's parameters is: where we use Λ = g * f and g * = 3. In fig. 8 (left) we show the exclusion limits from the combination of Higgs and EWPT data. Solid (dashed) lines correspond to a value of θ t = π/2 (π/4) while we fix θ b = 0. The reason why in the combined fit for values β ≈ 0.1 the limit on ξ relaxes up to almost the maximum value allowed by Higgs data is due to the fact that a small but non-zero β induces a positive contribution toT , which helps to relax the EWPT bound, as is shown in the right panel for the specific benchmark point used throughout the paper.

C Expression of the effective couplings
In order to describe the effective NGB interaction couplings in terms of the parameters of the potential, it is convenient to introduce a new set of coefficients. Since each invariant generating the potential brings some power of the fundamental couplings, it is useful to define: (1,0) y 2

Lc
(2) g g 2 In the following, we provide the list of all the effective couplings used in the computations. The couplings have been listed separately depending on the particles involved in the interaction.

C.1 Interactions between NGBs and gauge bosons
The first set of couplings is generated by CCWZ and involves the massive gauge bosons: As mentioned in the text, in order to properly discuss LHC phenomenology the couplings of the NGB to the gluons generated by loop of top-quarks has to be included: where f (τ X ) is the usual function appearing for gluon effective couplings (see for instance eq. (1.198) of [55]). Notice that since this vertex is generated by a quark-loop the resulting coupling for H 0 and η will be suppressed by a factor ξ with respect to the one to the SM Higgs and so it is not expected to play a relevant role, although it has been included in our simulations.

C.2 Interactions between NGBs and fermions
Interactions with quarks are generated by the partial compositeness Lagrangian and depend on the specific embedding of the fermions: with α q = 1 (−1) for quarks with charge 2/3 (−1/3).

D Relic density
The main contributions to the relic density are given by η annihilations into h, W , Z, t and b; below the W, Z threshold, however, also the channels ηη → W W * , ZZ * have to be included. With the vertices given in appendix C, the thermally-averaged cross sections are: with α V = 1 (1/2) for W (Z). Finally, also the process ηη → V V * can play an important role below the W and Z bosons production threshold. The thermally-averaged cross section for this process is (in this case, the exchange of H 0 in the s-channel is completely negligible): where N (f ) c is the number of colors of the final state f , and: with: ε V ≡ m V mη and ζ f ≡ (m f 1 + m f 2 )/(2m η ), f 1 and f 2 being the final states of V * decay. Obviously, also the coefficients τ (V ) and χ (V ) depend on the final states.
The relic density is computed from the effective cross section as: σv eff 1 pb , (D.5) with: The same formulas also apply for the case of the non-thermal contributions of section 5, where also κ freezes-out before decaying into η.