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)\times SU(4) \to 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.
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 [1] and WMAP data [2]. In a remarkable attempt to address the issue [3], 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). 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 [4][5][6], 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 [7,8]. Therefore, the theory naturally embeds two mass scales: one, the strong-coupling scale, is generated dynamically, while the other is related to the heavy quark mass. Hence a salient feature of such SQCDembedded 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 [9,10].
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.
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 [11,12]. 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) [13][14][15][16][17] 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 [18][19][20] within the MSSM. Here we consider an alternative scenario where the LSP relic abundance is small, while NGB-LSP interaction 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 [1]. 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 [3]. We consider the existence of a strongly coupled supersymmetric SU (N ) gauge sector having N f flavors of quark superfields denoted by Q i and Q 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 [4][5][6], 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 [22,23] 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 [3] 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 (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 5 having the superpotential (the relation to Eq.(2) is, T = M/Λ 0 , B = b/Λ 3 0 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.

Realization of the effective superpotential from
The low energy version of the supersymmetric SU (N ) gauge theory with N f = N + 1 (where N = 4) flavors is associated with mesons and baryons which are defined analogously to Eq.
(3), but due to the presence of an extra flavor (i = 1, 2, ...5), the baryons carry a free flavor and the meson matrix (T ij ) becomes correspondingly larger. Hence the baryonsB i andB i transform under the global group, as (N f , 1) and (1,N f ) respectively. Following Seiberg's [4,5] prescription, the system can then be represented by the superpotential, 5 Here the Q denote the SU (4) quark super-fields.
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 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 Comparing the above expression with Eq.(4), we can now identify B =B 5 ,B =B 5 and S = t. Hence the effective mass parameter involved in Eq.(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.(4) of N f = N (= 4) SQCD or equivalently to Eq.(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 [4][5][6] . Here we are interested in the specific point on the quantum moduli space (a la Eq.(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, 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 [7,8]. 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.(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 [7,8], 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 for σ 2 Λ eff Λ 0 .
Within the slow-roll approximation, the amplitude of curvature perturbation ∆ R , spectral index n s , and the tensor to scalar ratio r are given by where and η are the usual slow roll parameters. Assuming φ = M GUT = 2.86 × 10 16 GeV, then ∆ R = 2.2 × 10 −9 [1] implies Λ 0 4.3 × 10 17 GeV and Λ eff 1.8 × 10 15 GeV. When the number of e-folds is N e = 57, smooth hybrid inflation predicts n s 0.967 [24] and the very small value r 3 × 10 −6 [24] that are in good agreement with Planck 2016 data [1].
However, supergravity corrections to inflationary potential in Eq.(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 [24]. Particularly the value of n s ∼ 0.99 is in tension with observations [1]. To circumvent the problem one may use non-minimal Kähler potential as suggested in [10] 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 [9].

NGB as dark matter in SQCD
Once the inflation ends, the inflaton system slowly relaxes into its supersymmetric ground state specified by Eq. ( 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. (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 [11,12]: 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 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. (57)). 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. 10) and last terms in (1) R chiral symmetry, the middle term only respects the diagonal subgroup [25].
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 [3]; it can address the so-called µ-problem in the MSSM [26,27], 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 (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. (20)), is Higgsportal like: In terms of the physical mass eigenstates, the two Higgs doublets in MSSM can be written as follows [14][15][16][17]: 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 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.(22) as follows: where On the other hand, the couplings of the SM Higgs h to the vector fields and fermions are given by In view of Eq. (26), 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. (29)) and λ ∼ cos α(tan β − tan α) → 0 (see Eq. (28)), and hence is of limited interest. There is, however, a second solution: α β + π/2, so that tan β − cot α (see Fig.1) in which case the lightest CP-even scalar h will be SMlike and Eq. (27) closely resembles a Higgs-portal coupling [29][30][31]. 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.(27) These couplings are plotted in Fig.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. ( 27), 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 β.
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 [13,[18][19][20], 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 with mass m χ ) can be written as follows: where we assume m Gs > m χ . Here n G S and n χ 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 , χ} 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 pNBGs 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. (32) 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 →χχ ; 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. For example, consider the case where the LSP is dominated by the wino component. Then σv χχ→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.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 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 single-component Boltzmann equationṅ 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 [33] where ρ c = 1.05 × 10 −5 h 2 GeV c 2 cm −3 [32] is the critical density of the universe. Eq.(36) can be translated to T f with the freeze-out temperature T f . Next if one considers 6 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. (27). This situation is similar to the usual Higgs portal coupling with a singlet scalar, stabilized under a Z 2 symmetry [29][30][31] , with the important difference that the model we consider has two independent couplings, determined by λ and the angles α and β (cf. Eq. 28). 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. (38) are: The relic density Ω will be a function of the pNBG mass m Gs , the coupling λ and the angles α, β. In Fig. 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 [34].   and PANDA X [40] are shown with the expected sensitivity from XENONnT [41]. 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 [35,36]. In Fig. 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 [38], Xenon1T [39] and PANDA X [40] and future predictions from XENONnT [41] are also in the figure. Clearly, the figure shows that the degenerate case is excluded by 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 situaiont 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 in the game. which does not affect the direct detection cross section. We can allow a further smaller value of λ in this non-degenerate 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, therby 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 Fig.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 [1]. 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): 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. 20) 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; 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. 45 are coupled and must be solved numerically to find the freeze-out of the individual components. However, as it is shown in [38], for interacting multicomponent DM scenario, 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 , These approximate analytical results are in reasonably good agreement with the numerical solutions, as we show in Appendix III.
The DM phenomenology here crucially depends on the couplings λ and κ 1 , which we have varied freely for the scan. In Fig. 7   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 Fig 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 probablity 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 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 [37]. 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 [40] experimental limit.

Relic density for pNBG with small interactions to the MSSM sector
In this section we will consider some of the effects of the pNGB-LSP coupulings; 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 the evolution of the pNBG density is described bẏ 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 [30] 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. Figure 10: Feynman graph for Neutralino-pNBG interaction The Feynman graph responsible for the pNBG-LSP interactions is presented in Fig. 10.
The two vertex factors are named as α 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 [16]: 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. 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.(52) turns out to be a function of only Z 12 . In Fig. 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 Fig. 11.
The relic density for the case of 15 degenerate pNBGs as a function of pNBG DM mass is shown in Fig. 12, which takes into account the effects of the pNBG-LSP interaction. This  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.

Conclusions
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 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.
ACKNOWLEDGEMENT 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.

Appendix I: Dashen Formula:
Dashen formula for pNGB reads as When quark condensates like in form of ψψ = Λ 3 eff , we find for a=1, Correspondingly, (m 1 Similarly using other λ a 's for SU(4), we find (ζ = Λ/ χ ) The coefficients Z ij are the elements of the matrix which diagonalizes the neutralino mass matrix. In our set up we have considered the pNBG's interact with only the lighest neutralino by assuming rests are heavier than the pNBGs. The interaction lagrangian of χ 0χ0 h vertex is 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 well known 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 Fig.10 as obtained from Eq. (27). Using standard procedures we find Finally from Eq.(64) we obtain, 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 ) , For simplication purpose we have assumed C L = C R in our analysis. The Boltzman equations can be 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, 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 . Now we write the BEQ in terms of the newly defined quantities where the eqlibrium distribution has the form Y eq i (x) = 0.145 As we have already explored that the relic density of A-type DM having mass∼ m h 2 is negligible. So it freezes out at very later stage of the evolution of the universe. So during the freeze out of B and C types DM, it is quite trivial to write Y A Y eq A . Therefore, A-type DM does not have any role in the evolution of other two DM (B and C) candidates. We do further redefinitions and write the Boltzman equations as, where y eq i = 0.264M P √ g * µY eq i .
Once we obtain the freeze out temperature by solving the set of coupled equations (as in Eq.(71) numerically, we can compute the relic density for each of the DM species by