Dark matter as a remnant of SQCD inflation

We propose a strongly coupled supersymmetric gauge theory that can accommodate both the inflation (in the form of generalized hybrid inflation) and dark matter (DM). In this set-up, we identify the DM as the Goldstones associated with the breaking of a global symmetry (SU(4) × SU(4) → SU(4)) after inflation ends. Due to the non-abelian nature of this symmetry, the scenario provides with multiple DMs. We then construct a low energy theory which generates a Higgs portal like coupling of the DMs with Standard Model (SM), thus allowing them to thermally freeze out. While the scales involved in the inflation either have a dynamical origin or related to UV interpretation in terms of a heavy quark field in the supersymmetric QCD (SQCD) sector, the DM masses however are generated from explicit breaking of the chiral symmetry of the SQCD sector. We discuss DM phenomenology for both degenerate and non-degenerate cases, poised with DM-DM interactions and find allowed region of parameter space in terms of relic density and direct search constraints.


Introduction
Amidst the great success of the Standard Model (SM) of particle physics, there are several intriguing issues which indicate that the SM should be extended or supplemented by some other sector (including new fields and/or gauge symmetry) in order to provide a more complete description of nature. In particular, SM fails to accommodate a large share of energy density of the universe (25%), called dark matter (DM). On the other hand, the idea of primordial inflation serves as an elegant construction which can actually resolve some of the intricate problems (e.g horizon and flatness problems) of otherwise quite successful Big Bang cosmology. This inflationary hypothesis is further strengthened by its prediction on the primordial perturbation that leads to striking agreements with the observation of the cosmic microwave background radiation (CMBR) spectrum. However, SM alone can not be responsible for such a primordial inflation. In this work, we consider the existence of a different sector other than the SM, which can address primordial inflation and at the same time provides a suitable DM candidate. In [1], such an attempt to connect between primordial inflation and DM succesfully has been made.

JHEP10(2018)124
A successful model for inflation demands the existence a very flat potential for slowroll conditions to be satisfied during the inflationary epoch. Inclusion of supersymmetry protects the flatness of the scalar potential by non-renormalization theorem and is therefore a natural possibility. An inflation model embedded in a supersymmetric framework (for e.g., supersymmetric hybrid inflation models) usually contains one or more mass scales that are quite large compared to the electroweak scale, although smaller than the Planck scale as indicated by PLANCK [2] and WMAP data [3]. In a remarkable attempt to address the issue [4], it was shown that the inclusion of a hidden sector in the form of supersymmetric QCD (SQCD) can dynamically generate the scale of inflation, or relate it to the heavy quark mass of the electric theory (i.e. the theory at scales above the strong coupling scale). Inflation based on strongly coupled supersymmetric gauge theories has been studied in [5][6][7][8][9][10][11]. The properties of inflation in the SQCD framework are determined by the number of colors (N C ) and number of flavors (N f ) in the model. If one initially chooses N f = N C +1 [12][13][14], and then introduces a deformation of this pure SU(N C ) gauge theory by assuming the presence of one massive quark, then, upon integrating our this heavy quark, one obtains an effective superpotential of exactly the type used in smooth hybrid inflation [15,16]. Therefore, the theory naturally embeds two mass scales: one, the strongcoupling scale, is generated dynamically, while the other is related to the heavy quark mass. Hence a salient feature of such SQCD-embedded inflation model lies in the existence of a UV completion of the theory. Although in its original form the framework of smooth hybrid inflation embedded in supergravity is not consistent with all observational constraints, this can be corrected by considering a modified Kähler potential as shown in [17][18][19][20].
In this work, we demonstrate that DM candidates can arise from the same SQCD sector. We first note that at the end of the smooth hybrid inflation, one field from the inflation system gets a vacuum expectation value (VEV) and therefore breaks the associated global symmetry of the SQCD sector yielding a number of Nambu-Goldstone bosons (NGBs), whose interactions with the SM are suppressed by the spontaneous symmetry breaking scale, which is large, being related to the scale of inflation. Nambu-Goldstone bosons as DM has already been studied in some other contexts, see for example, [1,[21][22][23][24][25][26][27][28][29].
We then assume a small deformation of the UV theory by providing small masses (small compared to one heavy quark mass already present in the construction with N f = N c + 1 gauge theory) to the SQCD fermions; this generates the masses of NGBs through Dashen's formula [30,31]. Even with this deformation the NGBs are naturally stable and therefore serve as weakly interacting massive DM candidates; we will discuss this possibility in detail for various mass configurations. In the non-degenerate case, we show that DM-DM interactions play a crucial role in the thermal freeze-out, and therefore in relic density, and also impose the spin-independent (SI) direct search constraints from the XENON 1T.
We will assume that visible matter is included in a supersymmetric sector that is well described by the minimal supersymmetric Standard Model (MSSM) [32][33][34][35][36] at low energies. In this case the lightest supersymmetric particle (LSP) also serves as DM candidate in this framework. The case where this LSP contributes an important share of the relic density has been explored in various publications [37][38][39][40] within the MSSM. Here we consider an alternative scenario where the LSP relic abundance is small, while NGB-LSP interaction JHEP10(2018)124 plays a crucial role in surviving the direct search constraints, particularly for degenerate NGB DM scenario.
The paper is organized as follows. In section 2, we discuss the basic SQCD framework, which leads to the smooth hybrid inflation. We point out the prediction of such smooth hybrid inflation model in view of PLANCK result [2]. Then in section 3, NGBs are identified as DM and a strategy for introducing DM masses are discussed. We indicate here on the superpotential that would be responsible for generating DM interaction with the SM particles. Parameter space scan for relic density and direct search constraints on the model is elaborated in section 4 and we finally conclude in section 5.

Smooth hybrid inflation in SQCD
We start with a brief introduction to the SQCD framework that leads to a smooth hybrid inflation as was proposed in [4]. We consider the existence of a strongly coupled supersymmetric SU(N ) gauge sector having N f flavors of quark superfields denoted by Q i andQ i (i = 1, . . . , N f ) transforming as fundamental (N ) and anti-fundamental (N ) representation of the gauge group SU(N ) respectively. This theory also has a global symmetry: where the first U(1) is proportional to the baryon number and the second one is related to the anomaly-free R-symmetry. We are particularly interested in N f = N case, where in the electric (or UV) theory, the following gauge invariant (but unnormalized) operators can be constructed: Here a i correspond to the color indices and i j denote the flavor indices.
Classically, in absence of any superpotential, these invariant operators are required to satisfy the gauge and flavor-invariant constraint detM − bb = 0. As explained in [12][13][14], this constraint is modified by nonperturbative quantum contribution and becomes where Λ is a dynamically-generated scale. The corresponding quantum superpotential can be constructed by introducing a Lagrange multiplier field X, carrying R charge of 2 units, and is given by The necessity of introducing X follows from the fact that the expression of the quantum constraint does not carry any R charge and the superpotential W should have a R charge of 2 units. With N f = N = 2, it was shown in [41,42] that such a superpotential results into a low energy effective superpotential that is very much similar to the one responsible for supersymmetric hybrid inflation. However, since in this case the predictions of the supersymmetric hybrid inflation are not in accordance with the results of WMAP and PLANCK, we use [4] instead N f = N = 4, for which the corresponding effective superpotential still resembles the one in the smooth hybrid inflation scenario. In this case, the low energy JHEP10(2018)124 (or IR) theory below the strong coupling scale Λ 0 of the SU(N = 4) gauge theory, can be described in terms of meson fields T ij , baryon B and antibaryonB superfields fields 1 having the superpotential where Λ 0 is the strong coupling scale of the N f = N = 4 theory; then S can be identified with X The effective mass scale Λ eff can be interpreted in terms of holomorphic decoupling of one heavy flavor of quark (heavier than Λ 0 ) from a SU(N ) SQCD theory with N f = N + 1 flavors as we discuss below.
and the meson matrix (T ij ) becomes correspondingly larger. Hence the baryonsB i andB i transform under the global group, SU(N f ) × SU(N f ), as (N f , 1) and (1,N f ) respectively. Following Seiberg's [12,13] prescription, the system can then be represented by the superpotential, where Λ 0 is the strong coupling scale of SU(N = 4) SQCD. Note that the superpotential will have an R charge 2 if we assign only theT N f N f (= t) meson to carry R = 2. Next we introduce a tree level quark mass term in the superpotential, where we assume m Q m 1,2,3,4 . Considering first the case m Q >Λ 0 and m i=1,2,3,4 = 0 (later we will discuss the effect of having nonzero m i m Q ), the F -flatness conditions where we define the meson matrix T ij =T ij with i, j = 1, 2, 3, 4 and T 55 = t. Hence after integrating out the heavy N f th (the 5th) flavor of quark, we are left with the following effective superpotential for N f = N = 4 SQCD

JHEP10(2018)124
Comparing the above expression with eq. (2.4), we can now identify B =B 5 ,B =B 5 and S = t. Hence the effective mass parameter involved in eq. (2.4) is determined by the relation, Λ 2 eff = m Q Λ 0 and the Lagrange multiplier field turns out to be proportional to the N f th meson of the N f = N + 1 theory.
We now turn our attention to the superpotential in eq. (2.4) of N f = N (= 4) SQCD or equivalently to eq. (2.8) and discuss the vacua of the theory. Different points on the quantum moduli space associated with this N f = N (= 4) SQCD theory exhibits different patterns of the chiral symmetry breaking [12][13][14] . Here we are interested in the specific point on the quantum moduli space (a la eq. (2.4)) where B =B = 0, and T ij = (Λ 0 Λ eff ) 1/2 δ ij , which is the global vacuum of the theory. The corresponding chiral symmetry breaking pattern is then given by, (2.9) Hence along the direction B =B = 0 (the so-called meson branch of the theory), the superpotential reduces to with detT = χ 4 and S = −t; this superpotential is the same as the one used in the smooth hybrid inflationary scenario [15,16]. The SQCD construction of the superpotential serves as a UV completed theory and also the scales associated are generated dynamically. Below we discuss in brief the inflationary predictions derived from this superpotential which can constrain the scales involved, Λ 0 , Λ eff .

Inflationary predictions
The scalar potential obtained from eq. (2.10) is given by where the real, normalized fields are defined as φ = √ 2χ and σ = √ 2S (we use the same letters to denote the superfields and their scalar components). As was found in [15,16], the scalar potential has a local maximum at φ = 0 for any value of the inflaton σ and there are two symmetric valleys of minima denoted by φ = ±Λ 0 Λ eff /( √ 3σ). These valleys contain the global supersymmetric minimum which is consistent with our chosen point on the quantum moduli space given by detT = Λ 2 0 Λ 2 eff . At the end of inflation σ and φ will roll down to this (global) minimum During inflation σ 2 Λ eff Λ 0 and φ is stabilized at the local minimum φ = ±Λ 0 Λ eff /( √ 3σ). The inflationary potential along this valley is given by
However, supergravity corrections to inflationary potential in eq. (2.13) have important contributions. With minimal Kähler potential, previously obtained values of n s and r change to 0.99 and 1 × 10 −6 respectively [43]. Particularly the value of n s ∼ 0.99 is in tension with observations [2]. To circumvent the problem one may use non-minimal Kähler potential as suggested in [20] to bring back the value of n s within the desired range. On the other hand, to increase r to a detectable limit, further modification in Kähler potential may be required, as shown in [17,19].

NGB as dark matter in SQCD
Once the inflation ends, the inflaton system slowly relaxes into its supersymmetric ground state specified by eq. (2.12). This however spontaneously breaks the associated flavor symmetry from SU(4) L × SU(4) R to the diagonal SU(4) V subgroup. This would generate fifteen Nambu-Goldstone bosons (NBGs) those are the lightest excitations in the model. The Lagrangian so far considered has only the usual derivative couplings among the NGBs, which are suppressed by inverse powers of χ ; in particular, they have no interactions with the SM sector, and they are stable over cosmological time scale. In the following we will introduce additional interactions that will modify this picture.
At the global minimum we can write where G a , (a = 1, . . . , 15) denote the NGBs, and λ a are the generators of SU(4). This expression also identifies χ as the equivalent of the pion decay constant. Up to this point the NGBs are massless, which can be traced back to the choice of m i=1,2,3,4 = 0 in eq. (2.6). If we relax this assumption (while maintianing m i m Q ) the chiral symmetry is broken (explicitly) and, as a consequence, the NGBs acquire a mass that can be calculated using Dashen's formula [30,31]: where, as noted above, χ corresponds to the decay constant, the "+" subscript denotes the anticommutation, andQ a = 1 2 d 3 xψ † γ 5 λ a ψ are the SU(4) A axial charges with the quark state ψ = (Q 1 , Q 2 , Q 3 , Q 4 ) T ; λ a are the SU(4) generators normalized such that Tr[λ a λ b ] = δ ab 2 (a, b, . . . denote generator indices with a = 1, 2, . . . , 15). The details of the mass spectrum of these pseduo-NGBs (pNGB's) are provided in appendix I.
The pNGB spectrum is determined by the light-mass hierarchy. Among several possibilities, we will concentrate on the following two cases: • The simplest choice is to take m 1,2,3,4 = m in m diag (see eq. (2.6)). In this case all the fifteen pNGBs will be degenerate, having mass m 2 • A split spectrum can be generated if one assumes m 1 = m 2 = m 3 = m γ and m 4 m γ . Then we find three different sets of pNGB with different masses: (i) eight pNGBs will have mass m 2 We now turn to the discussion of the interactions of these pNGBs with the SM in this set-up. As noted in the introduction, we assume that at low energies the SM is contained in the MSSM. In this case, the S field (being neutral) serves as a mediator between the MSSM and SQCD sectors which then leads to the following modified superpotential where H u and H d are the two Higgs (superfield) doublets in MSSM, and κ 1,2 are phenomenological constants that can be taken positive. The terms within the curly brackets are phenomenologically motivated additions. Note that while the first (as in eq. (2.10)) and last terms in W Inf respect the full SU(N f ) L × SU(N f ) R × U(1) B × U(1) R chiral symmetry, the middle term only respects the diagonal subgroup [44]. The inclusion of such terms hence naturally lead to additional interactions of pNGBs. It is important to note that during inflation T is proportional to the unit matrix, so the term proportional to κ 1 vanishes and does not affect the smooth hybrid inflation scenario described before. Below we will see how incorporation of such terms can lead to the desired DM properties.
The term in W T proportional to κ 2 provides a connection between the SQCD sector and the minimal supersymmetric Standard model (MSSM). Inclusion of this term in the superpotential will have a significant effect on reheating after inflation [4]; it can address the so-called µ-problem in the MSSM [45][46][47], provided S acquires a small, ∼ O(TeV), vacuum expectation value. We will see that it also provides a useful annihilation channel for DM.
The linearity on S in these new contributions to the superpotential is motivated from R symmetry point of view. We note that any dimensionless coupling multiplying the first term in W Inf can be absorbed in a redefinition of Λ 0,eff ; in contrast, the couplings κ 1,2 JHEP10(2018)124 (assumed real) are physical and, as we show below, are constrained by observations such as the dark matter relic abundance and direct detection limits.
The scalar potential can now be obtained from V scalar = |∂W T /∂S| 2 + |∂W T /∂T | 2 . At the supersymmetric minimum, it is given by In obtaining this, we expanded T in powers of the NGBs: and χ is developed around its expectation value: Note that the first term in V s is of interest only when the pNGBs are not degenerate as otherwise it would not contribute to number changing process. In such a case with nondegenerate pNGBs, the heavier Gs can annihilate into the lighter ones. Hence in such a situation, κ 1 can also play a significant role in our DM phenomenology along with κ 1 κ 2 /2. In fact, we will show that the annihilation of the heavier ones to the lighter components will aid in freeze-out of the heavier component. This helps in evading the direct search bounds as the coupling κ 1 alone will not contribute to direct search cross section, thereby allowing a larger parameter space viable to our DM scenario. It is important to note here that even if we assume that the masses of the heavier pNGBs are very large, their annihilation cross-sections to the SM will be small enough for an early freeze-out, which leads to an unacceptably large relic abundance unless a large enough κ 1 allows them to annihilate to lighter ones. Note that the interactions among the G generated by Tr(∂ µ T † ∂ µ T ) are negligible since they are suppressed by powers of χ .
The interaction of the pNGBs with the MSSM sector (the last term in eq. (3.4)), is Higgs-portal like: In terms of the physical mass eigenstates, the two Higgs doublets in MSSM can be written as follows [33][34][35][36]: where h and H denote the light and heavy CP-even eigenstates respectively; H ± and A are the charged and CP-odd physical scalars respectively, and X 0,± are the would-be Goldstone bosons. As usual, h plays the role of the SM Higgs, and the vacuum expectation values of H 0 u , H 0 d (denoted by v u and v d respectively) are related by JHEP10(2018)124 The other mixing angle α appears as a result of the diagonalization of the CP-even Higgs mass-squared matrix (in the H 0 u − H 0 d basis) leading to the physical Higgses, h and H. The mixing angle α can be expressed in terms of β and the pseudoscalar A mass as We can now write the coupling of the pNBGs with the SM Higgs from eq. (3.6) as follows: On the other hand, the couplings of the SM Higgs h to the vector fields and fermions are given by In view of eq. (3.10), it can be noted that in the large pseudoscalar Higgs mass limit M A M Z , tan 2β tan 2α, that has two solutions. One possibility is α β, in which case couplings of the SM Higgs with W and Z vanish (see eq. (3.13)) and λ ∼ cos α(tan β − tan α) → 0 (see eq. (3.12)), and hence is of limited interest. There is, however, a second solution: α β + π/2, so that tan β − cot α (see figure 1) in which case the lightest CPeven scalar h will be SM-like and eq. (3.11) closely resembles a Higgs-portal coupling [48][49][50]. In the following we will continue to assume α = β + π/2. It is easy to show that with such a choice, λ , λ and λ are related with β as (using eq. (3.11) JHEP10(2018)124 Figure 2. Variation of λ and λ ( scaled by λ ) with tan β. We assume α = β + π/2 for both the cases.
These couplings are plotted in figure 2 as functions of tan β (note that v d is a monotonically decreasing function of tan β).

Relic density and direct search of pNBG DM
We now turn to the determination of the regions of parameter space where our DM candidates, the pNGB's, satisfy the relic density and direct detection constraints. The pNBGs interact with the SM through the Higgs portal as described in eq. (3.11), so the behavior of the pNBG sector is similar to that of a scalar singlet Higgs portal. In view of our consideration, α = β + π/2, the parameters for this sector are then effectively m Gs , λ, tan β. (4.1) The other factor which determines the relic density of pNGB's is their mass spectrum; we will show that depending on the choice of mass parameters employed in the Dashen formula, we can have several phenomenologically viable situations where one or more of the G contribute to the relic density. We hasten to note, however, that this model also contains an additional particle in the dark sector: the lightest supersymmetric particle (LSP), that we assume to be the lightest neutralino (χ 0 ), which is stable due to R parity conservation. The MSSM neutralino χ 0 has been studied as a DM candidate in various contexts [37][38][39][40], and can annihilate to SM and supersymmetric particles through many different but known channels, depending on the composition of Bino-Wino-Higgsino admixture. In this work we will concentrate on the case where the pNBGs dominate the DM abundance. We still must incorporate neutralino and MSSM phenomenology to some extent as the pNBGs interact with the neutralino DM through the Higgs portal coupling. This is because, there is a Wino-Higgsino-Higgs be completely neglected and the DM relic density is solely composed of pNBGs, and (ii) when the pNBG-LSP interaction is weak but non-vanishing.
In general, the coupled Boltzmann equations that determine the relic density for our two-component DM model (considering the presence of a single pNGB having mass m Gs and LSP with mass m χ 0 ) can be written as follows: where we assume m Gs > m χ 0 . Here n G S and n χ 0 denote the number density for the pNBG and LSP, respectively. Corresponding equilibrium distributions are given by where the index represents either of the DM species i = {G s , χ 0 } and ξ i denotes the degrees of freedom for the corresponding DM species, and we assume zero chemical potential for all particle species. In principle annihilations of the pNBGs occur to both SM and supersymmetric particles. However, assuming that the supersymmetric particles are heavier (as searches at LHC have not been able to find them yet), the dominant annihilation of the pN-BGs occurs to SM particles. After freeze out, relic density of the DM system is described by where the individual densities are determined by the freeze-out conditions of the respective DM components; which, in turn, are governed by the annihilation of these DM components to the SM, and as well as by the interactions amongst themselves. One can clearly see from eq. (4.2) that the Boltzmann equations for the two dark sector components, LSP and pNGBs, are coupled due to the presence of the terms containing σv G S G S →χ 0 χ 0 ; when this is of the same order as σv G S G S →SM the individual abundances will differ significantly from those obtained when σv G S G S →χ 0 χ 0 0. For example, consider the case where the LSP is dominated by the wino component. Then σv χ 0 χ 0 →SM is significant because of the large coupling of the wino to the Z, and from the co-annihilation channel involving the lightest chargino. Because of this large cross section we expect that the LSP relic density will be small: Ω χ 0 0.1 and the relic density will be composed almost solely of pNBGs: Ω G S h 2 ∼ Ωh 2 ∼ 0.1. In the following we will consider separately the cases where the LSP-pNBG interactions are negligible and when they are significant.
As stated earlier, our scenario allows for 15 pNBGs which can be degenerate. Hence their total contribution to Ω will be 15 times that of a single boson. However, the degeneracy of the pNBGs follow from making the simplifying assumption The simplest case we consider is that of completely degenerate pNGBs (which corresponds to m 1 = m 2 = m 3 = m 4 ); in this case the 15 G S contribute to DM relic density equally. In the absence of pNGB-neutralino interactions the individual abundance by the singlecomponent Boltzmann equatioṅ where the interactions between different pNBGs are also ignored because of the degeneracy. Therefore, total relic abundance can be obtained by adding all of the single component contributions. Now the relic density for a single component pNGB is given by [51] where ρ c = 1.05 × 10 −5 h 2 GeV c 2 cm −3 [52] is the critical density of the universe. Eq. (4.6) can be translated to where T f with the freeze-out temperature T f . Next if one considers 2 x f = 22, the total relic density turns out to be As mentioned earlier, the annihilation of pNBGs to SM states is controlled by two couplings λ and λ through eq. (3.11). This situation is similar to the usual Higgs portal coupling with a singlet scalar, stabilized under a Z 2 symmetry [48][49][50], with the important difference that the model we consider has two independent couplings, determined by λ and the angles α and β (cf. eq. (3.12)). In the limit where the psudoscalar mass is large, which we assume, the number of parameters is reduced by one because of the relation tan 2α = tan 2β. This results into the phenomenologically acceptable relation involving α and β as α = β + π/2.
The important cross sections contributing to eq. (4.8) are: JHEP10(2018)124  The relic density Ω will be a function of the pNBG mass m Gs , the coupling λ and the angles α, β. In figure 3 we evaluate Ω as a function of m Gs for 0.05 < λ < 1.0 and for tan β = 10 when α = β + π/2; the evaluation is obtained using MicrOmegas [53]. Note here, that the dominant annihilation of pNBG DM comes through λ . But the change in λ due to change in tan β is neatly balanced by the change in v d accompanying λ in all vertices making the relic density invariant under tan β. The results exhibit the usual resonant effect when m Gs ∼ m h /2 ∼ 62.5 GeV. We can also see that as λ increases Ω drops, a consequence of having larger cross sections. The allowed region in the m G S −λ plane by the DM relic density constraint is presented shown in figure 4 for tan β = 10 when all 15 pNBGs are degenerate. The allowed parameter space is similar to that of a Higgs portal scalar singlet DM. The difference is mainly due to having 15 particles: we require larger cross section, corresponding to a larger value of λ, to compensate for the factor of 15 in eq. (4.8).

JHEP10(2018)124
The non-observation of DM in direct search experiments imposes a very strong constraint on DM models, ruling out or severely constraining most of the simplest singlecomponent frameworks. It is therefore very important to study the constraints imposed on the pNBG parameter space by direct search data. Direct search reaction for the pNBG DM is mediated by Higgs boson in t-channel as in Higgs portal scalar singlet DM. Spinindependent direct search cross-section for pNBG DM is given by: where The subindex n refers to the nucleon (proton or neutron). We use default form factors for proton for calculating direct search cross-section: f p Tu = 0.0153 , f p T d = 0.0191 , f p Ts = 0.0447 [58,59] in micrOmegas. In figure 5, we show the spin-independent nucleon-pNBG DM cross-section for the chosen benchmark point, plotted as a function of DM mass m G S . The green points in this figure also meet the relic density constraint. Bounds from LUX 2016 [54], Xenon1T [55] and PANDA X [56] and future predictions from XENONnT [57] are also in the figure. Clearly, the figure shows that the degenerate case is excluded by JHEP10(2018)124 the direct search bound except in the Higgs resonance region. This is simply because for 15 degenerate DM particles, the values of λ required to satisfy relic density constraint correspond to a direct detection cross section large enough to be excluded by the data (except in the resonance region).

Non-degenerate pNBGs
We next consider the model when the pNGBs are not degenerate; we will see that in this case the allowed parameter space is considerably enlarged. For this it is sufficient to consider cases where there is a single lightest pNGB, and the simplest situation in which this occurs is when m 1 = m 2 = m 3 = m and m 4 = m. In this case the pNGB spectrum is type # of degenerate pNGBs mass and there are two cases: In the first case there is a single lightest pNGB (type C), while in the second there are 8 lightest states (type A).
Note that due to the presence of three types of pNGBs, the total relic density should be written as Ω T = n A Ω A + n B Ω B + n C Ω C , where n A/B/C are the number of degeneracies of the respective species. Compared to the case with all degenerate pNGBs, here the coupling κ 1 also comes into play along with λ. In case of all degenerate pNGBs, we have seen that the relic density satisfied region was ruled out by the direct detection cross-section limits excepting for Higgs resonance. This is because a large λ, as required to satisfy relic density, makes the direct detection cross section significantly higher than the experimental limits. Here we can make λ relatively free as κ 1 also enters into the game. which does not affect the direct detection cross section. We can allow a further smaller value of λ in this nondegenerate case, provided not all 15 pNGBs may not effectively contribute to the relic density. This can happen once we put mass of one type of pNGBs (out of A,B and C) near the resonance region ∼ m h /2 where the annihilation cross-section is large to make the corresponding relic density very small. In this case, the total relic density gets contribution from the two remaining types of pNGBs, thereby a smaller λ (compared to all degenerate case) can be chosen. Note that case II would be more promising compared to case I from this point of view as by putting m A m h /2, we effectively have remaining 7 pNGBs to contribute to relic. In figure 9, we demonstrate the masses, number of degeneracies and possible interactions among the pNGB DM candidates for the this case II. In this case m B (4m C + m h )/6 and m C > m h /2 and the relic density reads (since Ω A 0), with 0.1175 < Ω T < 0.1219 following PLANCK data [2]. In order to obtain Ω C,B we consider the set of coupled Boltzmann equations for the C, B i and A i number densities (the last included for completeness):

JHEP10(2018)124
In above we have ignored the interactions with the LSP. The numerical factors correspond to the number of final state particles each pNBG species can annihilate to. As is evident, a crucial role is played by the DM-DM contact interactions generated by (see eq. (3.4)) Since we assume that the A mass lies in the Higgs resonance region, σv G A G A →SM is very large, and produces very small relic density which we neglect in the following estimates;

JHEP10(2018)124
one can easily calculate that with m A ∼ m h /2, the A contributes less than 1% to the total relic if λ > 10 −3 . We will use this value of λ as a lower limit in our analysis.
Even if we neglect contributions from A type, pNBG, eqs. (4.15) are coupled and must be solved numerically to find the freeze-out of the individual components. However, as it is shown in [60], for interacting multicomponent DM scenario [61,62], annihilation of heavier components to lighter ones are crucial in determining the relics of only heavier components, while for lighter components it has mild effect. Hence we can derive an approximate analytic expressions for the individual relic densities by considering the annihilation of one pNGB kind to the lighter species. In this case we find , , In figure 8, we also show the relative contribution of relic density of one type of DM for a fixed κ 1 ; on the left (right) side, we choose κ 1 = 0.35 (κ 1 = 0.45). Contributions from different values of λ are shown in different colours. We note that Ω C yields the dominant contribution to Ω T , which occurs because C annihilation cross section is larger than for the B, by the contribution from the CC → BB process and also due to the larger degeneracy (6) of B component. We also see that with larger κ 1 , the regions of larger λ disappear (i.e. are inconsistent with the constraints).
Finally we show spin independent direct search cross section of C and B-types of DM in figure 9. A large region of parameter space is allowed by the LUX limit; this is because the DM-DM conversion allows different pNBG species meet the required relic density, without contributing to direct search cross sections. Clearly depending on how large one can choose κ 1 , the mass of the DM gets heavier to satisfy relic density and direct search constraints. Also due to the larger DM-DM conversion cross-section, the direct-detection probability for the C type of pNBG is smaller than for the B type.
We conclude the section by discussing a specific region of parameter space of the MSSM where the LSP is dominated by the wino/wino-Higgsino component, and also contributes negligibly to relic density. The LSP has four different contributions from the two Higgsinos  (H u,d ), wino (W ) and bino (B) [33][34][35][36]: where the Z 1j represent mixing angles. Now, since we are interested in the case where the DM density is dominated by the pNBGs, the relic density of neutralino (Ω χ 0 ) has to be very small. This is possible when the neutralino is generally dominated by the wino component or wino Higgsino components [63]. We tabulize two such examples (table 1) where we assume squarks, sleptons and gluinos of the order of 2 TeV, tan β = 5 and all trilinear couplings at zero, excepting A t = −1000 to yield correct Higgs mass. We find that the contribution of the LSP (neutralino) to the relic density is small and also the spin independent direct detection cross section σ SI χ 0 χ 0 →N N is well below the PANDA X [56] experimental limit.

Non-negligible pNBG-LSP interaction limit
In this section we will consider some of the effects of the pNGB-LSP couplings; for simplicity we will assume that the fifteen pNGBs are degenerate (it is straightforward to relax this assumption). As noted above, there are cases where the LSP receives a non-negligible contribution from the Higgsinos, in which case the LSP-pNGB interactions cannot be ignored, even though the DM relic density is still dominated by the pNBGs. In this case JHEP10(2018)124    Table 1. Relic density and corresponding SI direct detection cross section for wino/wino-Higgsino dominated neutralino with tan β = 5. The input parameters in terms Bino (M 1 ), Wino (M 2 ) and Higgsino (µ) masses and the output in terms of neutralino mass and mixing parameters are indicated. All the masses are in GeVs. Spin independent direct search cross section is in cm 2 . the evolution of the pNBG density is described bẏ

JHEP10(2018)124
where in the last term, we used n χ 0 = n eq χ 0 for the neutralinos since in the parameter region being considered they interact sufficiently strongly with the standard model to ensure they are in equilibrium; the decoupling of the LSP from the SM occurs much later than the decoupling of the pNBGs. We also assumed here that pNGBs are heavier than MSSM neutralino. Using then standard techniques [49], we find that the DM relic abundance is given approximately by where the presence of the second term in the dominator will be instrumental in accommodating the direct detection constraints and the numerical factor 15 to take care of fifteen degenerate pNGB species. The Feynman graph responsible for the pNBG-LSP interactions is presented in figure 10. The two vertex factors involved in the process are α 1 = 2λ v d and α 2 (elaborated below). This leads to the following annihilation cross section evaluated at threshold s = 4m 2 Gs : where m χ 0 is the neutralino mass and The vertex containing α 2 is generated by the Higgs-neutralino interaction [35]: where g 2 is the SU(2) L gauge coupling constant; as previously, we assumed α = β + π/2. The detailed calculation of cross-section is illustrated in appendix II. We perform a scan of the pNBG DM parameter space by varying DM mass and coupling of pNGB with SM (proportional to λ) with pNBG DM-Neutralino annihilations into account. For that we choose two benchmark points by fixing tanβ = 5 and using α = β + π 2 as mentioned in table 2. This is following from what we obtained from table 1, where neutralino has minimal relic density, but sizeable pNGB-neutralino interaction. The two interaction coefficients (λ and λ ) in this approximation turns out to be: Now as we are working with the wino or wino-Higgsino dominated neutralino, it is safe to consider Z 11 = 0. For simplification purpose we also assume Z 13 = Z 14 . Next we employ the constraint |Z 11 | 2 + |Z 12 | 2 + |Z 13 | 2 + |Z 14 | 2 = 1. In that case C in eq. (4.22) turns out to be a function of only Z 12 . In figure 11 a line plot has been presented to show the variation of C with Z 12 . For our analysis, we choose C = 0.2 and Z 12 = 0.88 (wino dominated LSP) which has been denoted by a red dot in figure 11.    Table 2. Benchmark points used to find out the relic density and SI direct detection cross section of pNGB DM using pNGB-LSP interaction. We mention physical masses of neutralino (χ 0 ) and pseudoscalar Higgs A and mixing pattern of neutralinos. All masses are in GeVs.
The relic density for the case of 15 degenerate pNBGs as a function of pNBG DM mass is shown in figure 12, which takes into account the effects of the pNBG-LSP interaction. This is scanned for λ = 0.005−0.25 (blue), 0.25−0.5 (green) and 0.5−2 (purple). The graph clearly shows flattening of relic density lines (which corresponds to fixed values of λ) for m G S ≥ m χ 0 , when the pNBG→LSP annihilation channel becomes kinematically allowed. The allowed parameter space for the pNBGs in mass-coupling plane is shown in figure 13.
In addition to predicting the expected relic density, the model should also comply with the constraints from direct detection experiments. The restrictions imposed by the LUX 2017 [54], XENON1T [55] and PANDA X [56] experiments are presented in figure 14 for the cases of 15 degenerate pNBGs with sizeable interaction with neutralinos. The results clearly indicate that in the presence of pNBG-LSP interactions the direct detection impose milder restrictions on parameter space compared to those cases where this coupling is negligible (compare figures 5 and 14). In particular, the fully degenerate pNGB case can now comply with the direct search constraints. It is again worth emphasizing the role played by the neutralino: it provides a new channel through which the pNBGs can annihilate (pNBG→LSP) that does not affect the direct detection cross section, this allows smaller values of λ (thus relaxing direct detection constraints), while keeping a large enough annihilation cross section, needed to meet the relic abundance requirements. As noted earlier, this occurs in the region where the LSP is wino/wino Higgsino dominated and in this reigon of parameter space Ω LSP Ω pNBG Ω.

Conclusions
DM as pNGB, arising out of breaking of a continuous symmetry has been studied in literature. Similar ideas have also been exploited to realize composite DM, which indeed appeal to a lot of astrophysical observations like non-cuspy halo profile etc. (see for example, [64][65][66][67][68][69][70][71][72]). Relating this type of DM to a consistent inflationary picture where the existence of pNGB is an artifact of the breaking of the continuous symmetry at the end of inflation, is the most interesting feature of our study. In this work, we have made use of the pNBGs, which are part of an SQCD framework in realizing early Universe inflation, as dark matter candidates. Due to the non-abelian nature of the chiral symmetry which was broken spontaneously at the end of inflation, a multiple such pNBGs as dark matter follow in the set-up. We have shown that depending on the explicit chiral symmetry breaking term, there could be different degree of degeneracy among the masses of these pNBGs. In addition, the presence of R-symmetry preserving supersymmetric Standard Model induces JHEP10(2018)124 in general another candidate of DM, called the neutralino (being LSP). Hence we end up having a multi-particle DM scenario, which eventually is considered to be dominated by the pNBGs only as far as the contribution to relic density is concerned. We then divide our analysis into two parts; one is when the interaction between LSP and pNBGs is completely neglected and the other one is with non-zero but small LSP-pNBGs interaction. We find that the case with all degenerate pNBGs can not lead to a successful situation consistent with the recent direct detection limit in the first case. On the other hand, the case with non-degenerate pNBGs without any effective contribution of the LSP toward relic density can be consistent with direct search bound. In case of small but non-zero LSP-pNBG interaction, we have found that this possibility alters our previous conclusions significantly.
For example, the case of fifteen degenrate pNBG DM now becomes a possibility. Therefore our model provides an interesting possibility of pNBG dark matter scenario.

Acknowledgments
SB would like to acknowledge DST-INSPIRE Faculty Award IFA-13 PH-57 at IIT Guwahati. AKS would like to acknowledge MHRD for a research fellowship.

A Dashen formula
Dashen formula for pNGB reads as

B pNBG annihilation to neutralino
In MSSM, after electroweak symmetry breaking, neutral gauginos and Higgsinos mix to yield four physical fields called neutralinos. In the mass basis, the neutralino can be written as a combination of wino, bino and two Higgssions. For example, the lightest one can be written as where the coefficients Z ij are the elements of the diagonalizing mass matrix and crucially control its interaction to other MSSM and SM particles. In our analysis, we have considered that the pNBG's can only annihilate to lightest neutralino by assuming the rest to be heavier than the pNBGs. The interaction lagrangian of χ 0χ0 h vertex is L = −ig 2χ0 (C L1 P L + C R1 P R )χ 0 h. Our aim is to calculate the cross section for G S G S →χ 0 χ 0 annihilation process that has been used in the estimation of pNGB DM relic density. It is a two body scattering process and the differential scattering cross section in centre of mass frame is given by ( p f is the final state momentum). The Feynman amplitude for the process is where α 2 = g 2 (C L P L + C R P R ) and α 1 = 2λ v d are two vertex factors in figure 10 as obtained from eq. (3.11). Using standard procedure, we find where C L = −Q * 11 sin α − S * 11 cos α, C R = −Q 11 sin α − S 11 cos α, Q 11 = 1 2 Z 13 (Z 12 − tan θ W Z 11 ) + Z 13 (Z 12 − tan θ W Z 11 ) , S 11 = 1 2 Z 14 (Z 12 − tan θ W Z 11 ) + Z 14 (Z 12 − tan θ W Z 11 ) , For simplification, we have assumed C L = C R in our analysis.

C Numerical estimate of coupled Boltzmann equations
Relic density allowed parameter space of non-degenerate multipartite DM components of the model (section 4.1.2 and section 4.2) has been obtained by using approximate analytic solution (eq. (4.15)). Here we will explicitly demonstrate the viability of such analytic solution to the exact numerical solution of the coupled Boltzmann equations (BEQ) that defines the freeze-out of such non-degenerate DMs. We will illustrate the case of nondegenerate pNGBs, with negligible interactions to LSP (section 4.1.2). Let us define Y = n i /s, where n i is the number density of i'th DM candidate and s is the entropy density of the universe. The BEQ is rewritten as a function of x = m/T , where m is the mass of DM particle and T is the temperature of the thermal bath. As we have three DM candidates of type A, B and C, we use instead a common variable x = µ/T where 1/µ = 1/m A + 1/m B + 1/m C and 1/x = T /µ = 1/x 1 + 1/x 2 + 1/x 3 . Assuming m C > m B > m A , the coupled BEQ then reads where the equilibrium distribution has the form Y eq i (x) = 0.145

JHEP10(2018)124
We have already explained that the relic density of A-type DM having mass ∼ m h 2 is negligible due to resonance enhancement of annihilation cross-section. Therefore, it freezes out much later than B and C type DMs. During the freeze out of B and C type DM, we can then safely write Y A Y eq A and ignore its contribution to the relic density. The coupled BEQ then effectively turns to dy C dx = − 1 x 2 ( σv CC→SM + 8 σv CC→A i A i )(y 2 C − y eq 2 C ) where y i = 0.264M P √ g * µY i , (C.4) y eq i = 0.264M P √ g * µY eq i . (C.5)

JHEP10(2018)124
Once we obtain the freeze out temperature by solving the set of coupled equations (as in eq. (C.3)) numerically, we can compute the relic density for each of the DM species by y i µ m i x ∞ indicates the value of y i evaluated at µ m i x ∞ , where x ∞ denotes a very large value of x after decoupling. For numerical analysis we have taken x = 500 which is a legitimate choice. We scan for κ 1 ∼ 0.1 − 0.3 and m C = 100 − 500 GeV to find out relic density allowed points for C type pNGB DM (Ω C h 2 < 0.1). We have shown it in terms of λ − m C in figure 15. Analytical solutions (eq. (4.17)) for relic density is plotted in the same graph for comparative purpose. We see for low values of κ : {0.1 − 0.2} the numerical solution (in blue) and approximated analytical solution (in red) falls on top of each other with very good agreement. With larger κ : {0.2 − 0.3}, the separation between numerical (in grey) and approximate analytical solution (in green) increases mildly within ∆m DM ∼ 10 GeV.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.