Pionic dark matter

We study a phenomenological model where the lightest dark matter (DM) particles are the pseudo-Goldstone excitations associated with a spontaneously broken symmetry, and transforming linearly with respect to an unbroken group \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ {{\mathcal{H}}_{\mathrm{DM}}} $\end{document}. For definiteness we take \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ {{\mathcal{H}}_{\mathrm{DM}}} $\end{document} = SU(N) and assume the Goldstone particles are bosons; in parallel with QCD, we refer to these particles as dark-matter pions. This scenario is in contrast to the common assumption that DM fields transform linearly under the full symmetry of the model. We illustrate the formalism by treating in detail the case of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ {{\mathcal{H}}_{\mathrm{DM}}} $\end{document} = SU(2), in particular we calculate all the interactions relevant for the Boltzmann equations, which we solve numerically; we also derive approximate analytic solutions and show their consistency with the numerical results. We then compare the results with the constraints derived from the cold DM and direct detection experiments and derive the corresponding restrictions on the model parameters. We also briefly comment on constraints from indirect detection of DM.


Introduction
Dark matter (DM) is the most promising hypothesis proposed to explain astrophysical and cosmological observations related to the motion of stars in galaxies [1], the motion of galaxies in clusters [2][3][4], structure formation [5] and the inhomogeneities in the CMBR [6,7].Having not direct experimental information about this component of the universe the theoretical efforts to understand DM have been couched within realistic extensions of the Standard model (SM) [8][9][10][11][12][13][14], or have taken a purely phenomenological approach [15][16][17][18][19][20], in which case simplicity has been used as a guide and constraint.
In this publication we will investigate a phenomenological model for DM based on general assumptions concerning the dark sector, explicitly, we will assume that the lightest particles in that sector are the pseudo-Goldstone bosons resulting from a broken symmetry [21].Operationally this implies that the lightest particles (that we take as scalars for simplicity) transform non-linearly under a continuous symmetry group, a situation similar to the one occurring in low energy hadron physics.Accordingly, we will refer to them as dark matter pions (DMP) (we emphasize however, that these are quite distinct form the pions in the hadronic sector, in particular they do not have direct couplings to the standard model (SM) W ± and photon (in this we fundamentally differ from the assumptions made in [22]).This approach is in contrast with most phenomenological approaches where the dark-sector fields are assumed to transform under a discrete symmetry, or linearly under a continuous one [15][16][17][18].
In the following we will study this type of DM model based on the nonlinear realization of a spontaneously broken symmetry group G DM .However, given the difficulties of hot dark matter gas in dealing with structure formation [23], we will also assume that the Goldstone bosons receive their masses through an explicit breaking of the original symmetry.We also require that all SM particles are singlets under the dark-sector symmetries and that the dark particles are singlets under the SM local symmetries.
The interaction between these two sectors (SM and DM) is presumably effected by the exchange of some heavy mediators whose nature we do not need to specify, but only assume are much heavier than the typical scales in either sector.Therefore the typical interactions are of the form where O DM , O SM are operators invariant under the internal symmetries of the corresponding sector, but they need not be Lorentz invariant (though, of course, L DM−SM must be).The details of these interaacitons will be elaborated blow.This paper is organized as follows: in the next section we describe the formalism behind our model, and construct the Lagrangian we will use in our calculations.In section 3 we calculate the SM-DM interactions that we then use in sections 4 and 5 to derive the relic abundance of this type of dark matter.These results are compared with the experimental constraints in section 6 with our brief conclusions are presented in section 7. A few details are relegated to the two appendices.

Nonlinear realization of G DM
Models where the symmetry is nonlinearly realized have been extensively studied (see, e.g.[24]); here we summarize some of the results for completeness.We assume there is a subgroup H DM ⊂ G DM under which the vacuum is invariant and, following [24], we denote the generators of H DM by V i and the remaining generators of G DM by T a .Then the fields can be chosen as {π, ψ} with the following properties: • Under H DM they transform linearly: π → D(h)π, ψ → D(h)ψ for h ∈ H DM ; where D and D are some matrix representations of H DM .
• Under a general g ∈ G DM π → ξ(π, g) , ψ → D e u.V ψ ; u = u(π, g) , where D is the same representation as above, and ξ and u are defined by ge π.T = e ξ.T e u.V . (2.2) Note that the transformation of π depends only on g and π, and is non-linear; while that of ψ depends on g, ψ and π.Because of their transformation properties the π are massless and correspond to the Goldstone bosons generated under the spontaneous breaking G DM → H DM , and accordingly the number of these fields equals that of the broken generators T a .We will refer to the π as the "dark-matter pions" (DMP) or dark pions.
To be specific we concentrate on the familiar case [25,26] of a unitary chiral theory where G DM = SU (N ) × SU (N ) and H DM = SU (N ), the diagonal subgroup.In this case the above general formalism is realized by introducing a unitary field Σ and transforms as where Σ = exp(iπ.T/f ) and f is a mass scale associated with the spontaneous breaking of the symmetry.The diagonal subgroup corresponds to the choice R = L.
As it is well known [26,27], the leading fully chirally invariant operator is Expanding (2.4) in terms of the π we find that this Lagrangian describes a series of massless particles1 which are difficult (though not impossible [29]) to reconcile with structure formation.We will therefore also include an explicit breaking of the G DM symmetry that generate a mass for these excitations; for the chiral model this corresponds to a term of the form This term is invariant under the diagonal (unbroken) subgroup H DM .
In order to construct the DM-SM interactions of the form (1.1) we need the list of the lowest-dimensional SM gauge-invariant (tough not necessarily Lorentz invariant) operators.These are easily listed; for dimension ≤ 2 we have where φ denotes the SM scalar doublet and B the hypercharge gauge field.The dimension 3 operators (that we will not use here) are φ † D µ φ and ψγ µ ψ , where ψ and ψ are any two fermion fields carrying the same gauge group representation (e.g. e R and τ R ); higher dimensional operators are similarly constructed.Then, the simplest DM-SM coupling is clearly where v = φ ∼ 174 GeV.
The coupling Σ to B µν is less straightforward since there are no G DM -invariant operators that can be constructed out of Σ and its derivatives and which transforms as the (0, 1) + (1, 0) representation of the Lorentz group2 .Noting however, that (2.5) is invariant only under the diagonal subgroup H DM , we will only require the Σ − B coupling to have the same property, and in this case, For our choices of G DM and H DM the Lagrangian for our model is obtained from (2.4, 2.5, 2.7, 2.8); explicitly, where, as before, (2.10) In parallel with the usual strong-interaction pions, we will call f the DMP decay constant.The T a are the broken Hermitian generators normalized by and obeying (with a, b, . . .= 1, 2, . . ., N 2 − 1).In the Cartan basis with root generators T ±α and Cartan generators T i we have [30] [ where N α,β = 0 if α + β is not a root.We could also add another φ − π coupling by replacing To lowest order this coupling is of the form |φ| 2 π 2 and its effects have been studied extensively [31].Given our interest in studying the effects of the new interactions listed in (2.9) we will neglect λ V in the following.Writing Σ = exp(iσ) and using the Lagrangian can be written (in a Hermitian basis) where and h is the Higgs field; in unitary (SM) gauge φ T = (v + h/ √ 2)(0, 1).In the Cartan basis, and for the case of N = 2 (that we will develop later as a specific illustrative case): where π o is associated with the SU (2) Cartan generator, and π ± = π ±α , where α is the single root in this group.

Conserved currents
The Lagrangian (2.9) is invariant under the global transformations which give rise to a set of conserved Noetherian currents (2.21) Ignoring the interactions with the SM the canonical momentum are ℘ a = g ab πa in terms of which the charges (again ignoring the SM interactions) become and (ignoring possible sigma terms and other anomalies [32]) satisfy the algebra as expected.
The number of commuting conserved charges equals the rank of the group, which, in a Cartan basis, can be conveniently chosen as those associated with the π i : Assuming that these relations do not exhibit commutator anomalies [32] the charges Q i will be conserved; in particular this property will be reflected in the Boltzmann equations.It follows from the expression for Q i that the π i carry no charge, while π ±α carry opposite i-charges when α i = 0.

Parameters of the model
The model we consider has then 4 parameters: the DMP mass M , the DMP decay constant f , the coupling constant of the DMP to the Higgs λ h , and λ V , the coupling constant of the DMP to the hypercharge vector field B (from which follow the coupling to the Z boson and the photon).
In the calculations below we will take λ V coupling to be real with magnitude λ V = 0.63 . (2.25) We will see later that as far as the Boltzmann equations are concerned, any change in λ V can be absorbed in a redefinition of the other parameters (cf. the end of Sec.5), so this choice does not represent a loss of generality and is made for computational ease only.It is worth noting that according to naive dimensional analysis (NDA) [33] its value is λ V ∼ g /(4π) 2 0.0023, where g is the U (1) Y gauge coupling constant in the Standard Model.
For the rest of the parameters we impose just some loose constraints.We require that in order to ensure the model remain perturbative 3 .We will see later that all the experimental constraints on the model also have simple scaling dependence on the couplings λ h (see Sec.6.1),so this constraint will also not restrict the generality of our results.Since we assume that the DMP are the pseudo-Goldstone bosons of some underlying theory and are generated by the breaking of G DM to H DM at some scale Λ, consistency of the resulting chiral model requires [27] 4πf M . (2.27) For large values of N the left hand side is expected to be suppressed by a factor of 1/ √ N [34], which we do not include because we will restrict ourselves to low values of N .
Another constraint can be derived by requiring loop corrections not to dominate over the tree-level terms.In particular this should hold for the radiative corrections generated by the term proportional to λ V in (2.16), which includes vertices of the form (λ V /f n+2 )Z µν ∂ µ π ∂ ν π π n .Two such vertices will generate loop corrections to the ∂ µ π ∂ µ π π k /f k vertex of the first term in (2.16): where we have assumed that all the terms in (2.16) that explicitly violate G DM are associated with the scale M , which we have used as an UV cutoff.We require (2.28) not to be larger than the tree-level contribution, which implies (since L can be arbitrarily large) (2.29)

DMP interactions
In this section we calculate the cross sections for the processes that dominate the Boltzmann equations that describe possible equilibration between the dark and SM sectors, and within the dark sector.The relevant interactions (2.16) separate into those that involve only DMP, and those that involve DMP and the SM scalar φ or the vector boson B. We also derive the reactions relevant for direct detection of the DMP.In all the calculations below we only consider 2 → 2 processes and will use the Cartan basis for the DMP.

DMP → SM interactions
There are two kinds of reactions: Processes with only SM particles in the final state.These are of the form for which the interaction terms in (2.16) are and the processes are shown in Fig. 1.The cross sections for these processes are: where We neglected the Higgs width in the expression for σ(ππ → hh) since it is never resonant (resonance occurs at s ∼ m 2 h while the reaction occurs only if s > 4m 2 h ) and current data [35] suggests Γ h Γ (SM) h 4 MeV and m h = 125 GeV so that Γ (SM) h /m h 3.2 × 10 −6 .For the W , Z and t reactions we can also ignore Γ h in σ SM (defined in eq.3.4); the same is true for the other reactions if M > m h /2.
Processes involving DMP in the final state.These correspond to ππ ↔ πZ/γ for which the Lagrangian is given by So there are 3 types of reactions (the first present only for SU (N ), N > 2): σ where s V the number of spin degrees of freedom: s Z = 3, s γ = 2, and In the center of momentum (CM) frame K V = |k| = |l| denotes the magnitude of the V 3-momentum, and P = |p| = |q| the magnitude of the 3-momentum of the pions not paired with the vector boson:

Direct-detection reaction
The most important process that can contribute to the scattering of the DMP off heavy nuclei (relevant for direct DM detection [36][37][38]) is πψ → πψ, where ψ is SM fermion, and occurs through a t-channel h exchange.The averaged amplitude-squared is so that, in the CM frame, the corresponding cross section for this process is given by where P denotes the momentum of the incoming particles in the CM frame.When M, m h P, m f this cross section is approximated by At low momentum transfer the effective interaction obtained from integrating the Higgs using (3.2) and the Standard Model

Pure DMP scattering
Finally, we obtain the cross sections responsible for equilibrium within the DMP sector, ππ → ππ, Fig. 3.The lowest-order terms (taking M real) in (2.16) are where and we have dropped terms that vanish on shell and will no contribute to the S-matrix.In terms of DMP defined in the Cartan basis we have the following reactions:

Decays of SM particles to DMP
Limits on the DMP parameters can be derived either from collider reactions or from potential deviations from SM decays.Reactions of the form f f → ππ, where f is a SM fermion, or W fusion reactions W W → ππ, would mimic neutrino production at colliders.The limits, however, are very weak since these processes proceed through a virtual h and so the amplitude will be proportional to small Yukawa coupling, or, for the case of heavy initial quarks, suppressed distribution functions.
The main limits are then derived form the two leading decay processes, Fig. 4, namely, h → ππ and Z → πππ, to which we now turn.h → ππ decay Using (2.16) and choosing a Hermitian π basis we find that the width is given by Recent data [35] favors a Higgs decay close to the SM prediction of ∼ 4 MeV and a mass m h ∼ 125 GeV; this requires M > m h /2, or M < m h /2 and Γ hππ < 4 MeV, hence the constraint we use is In the numerical solutions to Boltzmann equations for DMP for the SU (2) case (discussed below), we consider DMP masses in the interval 50 GeV ≤ M ≤ 2000 GeV so the h → ππ constraint plays an important role only for comparatively small values of M .
Z → πππ decay The calculation is straightforward; using again a Hermitian DMP basis we find where s w = sin θ w , while E, K denote the usual Elliptic functions, and For H DM = SU (N ) and our normalization conventions (2.11,2.12) the summation involving the structure constants is given by Using the uncertainty in the invisible width of the Z, Γ(Z) inv we have limit which implies where The function Q is monotonic; it vanishes as r → 3 and approaches 0.75 as r → ∞.Taking N = 2, and λ V = 0.63, the most conservative limit (corresponding to taking When λ V = 0.063, this limit becomes f > 23.87 GeV.
In the numerical analysis, we choose to work with DMP mass ≥ 50 GeV and therefore the constraint from Z → πππ is of no importance.

Thermal history of DMP
We now turn to the derivation of the relic abundance of DMP.We follow the standard treatment (see e.g.[39]) and will consider only 2 → 2 processes.

Boltzmann equations
The change in the number density of particle of type a due to collisions and the expansion of the universe is given by ṅa where dΠ denotes the phase-space volume and g is the number of internal degrees of freedom.The amplitude-squared |A| 2 for the a + b → c + d process is understood to be averaged over initial and final states, and to include symmetry factors for identical particles in the final states.The functions f are the particle phase-space distribution functions; the corresponding particle number density is We will assume that interactions are such that kinetic equilibrium is maintained [40]; we will also assume that particles densities are sufficiently small to ignore the effects of quantum statistics.In this case the energy dependence in the distribution functions is given by the Boltzmann factor: f = ζ exp(−E/T ).Since we are interested in the epoch when the DMP first decouple, all distribution functions will have the same temperature T ; this will continue after decoupling provided no mass thresholds are crossed, or phase transitions occur.
The equilibrium distributions for a particle of mass m is given by where z is the fugacity in equilibrium.For the SM z SM = 1 to very good accuracy [41]; for the DMP, however, we will allow non-zero chemical potentials.Using the definition in (2.24) and the discussion below it, it follows that where µ a denotes the chemical potential for particle a associated with charge Q i so that z = 1 for those particles with non-zero conserved charges, as defined in Sec.2.1.
Substituting these definitions in the expression for C and using the standard definition of the scattering cross section σ we find In the definition of s o we used the condition (contained in the cross section) that s should be large enough to create c and d.
For the pure DMP scattering processes that appear in the Boltzmann equations the averaged cross sections can be evaluated in closed form.We obtain, for example with similar expressions for the other relevant processes; in deriving this we used (3.18) and (3.19).For the relevant initial states (π i π i or π α π −α ) we have z i = 1 = z α z −α so that in all cases of interest (see below) we can replace z πa z π b → 1.Also u is defined in (3.19), while B is defined as and µ is given in (3.15).In deriving the above result we used With the above preliminaries we can now find the relevant collision terms C a (c.fEq.(4.1)) for the cases a = π i and a = π α that we abbreviate as C i and C α respectively.We will assume that all SM particles remain in equilibrium, so that n SM = n (eq) SM .The tables of the relevant reactions (which do not cancel in C i,α ) are where V represents Z or γ, β = −α, and a summation over j and β is assumed.Now, using (4.6) and noting that (4.5) implies and similarly for the equilibrium densities, we find and where the contributions coming from π α V → π i π α (V = Z, γ) and π α π i → V π α cancel, as do those from π α π β → V π α+β and π α V → π β π α−β .We have also defined, using (4.6), and similarly for π α π −α → SM .

Contributions from SM→DMP decays
The effects of Higgs decays into DMP, when kinematically allowed, can be included in the Boltzmann equation in two equivalent ways.We can include them in the total h width: and use this expression in the cross sections involving Higgs exchange.Or, alternatively, we can exclude these effects from the Higgs propagators (see e.g.[42]) : and include them in suitable additions C (decay) i,α to the collision terms; explicitly (see Appendix A) where h counts the number of produced π i : If we assume that the recently observed particle at the LHC [35] is the SM Higgs, it's very small total width ensures that the effects from Higgs decay to DMP are negligible.We have checked that for realistic DMP masses the contribution of Z → πππ decays in the Boltzmann equations (see Appendix A) are also negligible.

Solving the Boltzmann equations for the SU (2) case
The simplest non-trivial group is H DM = SU (2), which we consider as an illustrative example of the formalism; the same approach can be used for any N , though with the calculations become increasingly cumbersome.For N = 2 there is a single conserved charge and 3 DMP states that we label as o, ±, with the first associated with the Cartan generator.
As usually we find it convenient to rewrite the Boltzmann equations (BE) (4.1, 4.13, 4.14) by defining where T denotes the photon temperature and s the entropy density: here k runs over all particles, T k is the temperature of particle k and g k its number of internal degrees of freedom, and r k = 1 (7/8) when k is a bosons (fermion).We will also make use of Friedman's equation, In the following we will take T k for all SM particles (assuming T is above that of the e + e − annihilation epoch), so that g s (T ) = g(T ); we use the expression for g(T ) in Ref. [43].The explicit form of the equilibrium distribution is where z r is the fugacity for particle r and g r the number of internal degrees of freedom.
We will also consider model parameters where the SM and DM sectors are in equilibrium for temperatures T > T f , such that T f > M , so that the region of interest is x > 1 and the DMP will not contribute 5 to the effective number of relativistic degrees of freedom g(T ) = g SM (T ).
In terms of Y the Boltzmann equations take the form where the collision terms are and where we used Y (eq) + Y (eq) − = Y (eq) o

2
, and also and defined (5.9) For the SU (2) case there is a single non-trivial chemical potential (4.5) and an associated conserved charge q = Y − − Y + . (5.10) Using q, the two independent Boltzmann equations become From (4.6) we find that where σ Z/γ are given in (3.8) and P, K V are defined in (3.9).The σv are plotted in Fig. 5 for a representative parameter space point.The SM cross section is almost x-independent (corresponding to a predominance of s-wave scattering), while the γ/Z cross section is proportional to 1/x, indicating a predominance of p-wave scattering.It is interesting to note that the DMP→DMP cross section has an unusual 1/ √ x behavior for large x that results from all particles having the same mass and the amplitude being non-zero and finite at threshold, which for this model is a consequence of the chiral couplings of the DMPs.One can see, that σv ππ→SM is much smaller than σv ππ→πV or σv ππ→ππ for the particular choice of parameters.The relevance of σv ππ→SM can be understood by referring to Fig. 6 where we compare σv ππ→SM and σv ππ→πV at the decoupling temperature (the point at which the DMP particle density begins to deviate significantly from its equilibrium value -see Sec.5.1) for points that satisfy the cold-dark matter (CDM) relic-abundance constraint (see Eq.(6.1) below).Cross sections for a representative set of parameters, (M, f, λ h , λ V ) = (1000 GeV, 950 GeV, 0.01, 0.63), for which the model satisfies the cold-dark matter and directdetection constraints.Top curve: 10 7 x σv ππ→πV ; middle curve: 10 12 σv ππ→SM ; bottom curve: 10 8 √ x σv ππ→ππ .The prefactors are chosen to fit the curves into the same graph and to illustrate the leading x behavior.All the cross sections are in GeV −2 .
To obtain the particle densities and their freeze out temperatures it is necessary to solve a system of coupled linear differential equations for {Y o , Y + } given by (5.11).The boundary conditions are determined by requiring that at low x the DM sector is in equilibrium with the SM: (5.13) Note that (5.11) and (5.13) imply that both the equations and initial conditions are invariant under Y + ↔ Y − and q ↔ −q.
For the following it is useful to note that σv πoπo→SM depends on λ h only in the combination λ h /f 2 , while σv π + π − →πoV depends on λ V only as λ V /f 3 .This implies that we can take M, f and λ h as independent parameters, fixing λ V at some convenient value as in (2.25); any other value of λ V can be obtained by appropriate rescaling of f and λ h .

Zero charge solutions
When q = 0 all DMP will have the same initial equilibrium distribution, the relevant solutions to the BE then correspond to Y o,± = Y ; substituting this (and q = 0) in (5.11) we find where we drop the o, ± subindices.Approximate solutions to this equation are readily obtained.We find that to good accuracy (see Fig. 5) the cross sections have an s and p wave behaviors for x > 10: where σ SM,V are approximately x-independent.
Finally the decoupling 'temperature' x f is obtained from the condition ∆(x f ) = cY (eq) (x f ), where c is a numerical constant.This gives where a is defined in (5.4) and ϑ SM , ϑ V in (5.16); this result is better suited for the case ϑ V ϑ SM than the one presented in [39].We will follow this reference and choose c(c+2) = 1 or, c 0.414.In calculating the relic abundance it is important to remember that Y ∞ refers to each DMP species, so that the total abundance will be proportional to 3Y ∞ .
An alternative definition of x f can be derived by assuming Y is close to Y (eq) and casting (5.14) in the form so x f can be defined as the point where Γ/H = 1.A plot of Γ/H for representative values of the parameters, and a comparison with the previous definition of x f is given in Fig. 7.This also illustrates that x f in general is large enough for the approximations (5.15) to be valid.In Fig. 8 we compare the relic abundance derived numerically with the one obtained from (5.18), showing that, at least in this instance, the latter is reasonably accurate.From this figure one can also see that the decoupling point inferred from the numerical solutions equals the analytically obtained values within 10%.

Behavior for small values of |q|
We now turn to the case where q is small but non-vanishing.In this case it is convenient to define in terms of which Eqs. (5.11) become where while the initial conditions (5.13) correspond to Now Y t,d are even in q, and assuming they are analytical in q it follows that they depend on q 2 ; at q = 0, we have Y t = 3Y and Y d = 0. Taking a derivative of (5.21) with respect to q 2 and evaluating at q = 0 gives where . Initially, (5.25) Now, a differential equation of the form has solution In particular, if v(x) > 0 for all x, and Z i > 0, then Z(x) > 0 for x > x i .Applying this to Z = (∂Y t,d /∂q 2 ) q=0 , that have initial values ∼ 1/Y (eq) o (x i ) > 0, we find that The relic abundance is obtained from the expression [39] Ω DM h 2 = 2.7711 since Y t (q = 0) > Y t (q = 0) (at least for small q and with the other parameters fixed), it follows that If Ω DM (f, M, λ h , λ V ; q = 0) < Ω CDM for some parameters {f, M, λ h , λ V }, then there will be a non-zero q such that Ω DM (f, M, λ h , λ V ; q) = Ω CDM .That is, if the predicted abundance falls below the observations when q = 0, one can always "make-up" the difference by introducing an appropriate q (at least when the difference is small).It follows that the the region in parameter space that can satisfy the CDM constraints is determined by A non-zero value of q does not, of course, affect the direct-detection probability.We illustrate Boltzmann equation solutions for small q in Fig. 9.In general, there is a small range of |q| ∼ 10 −12 − 10 −13 for which differences among the Y + , Y − and Y 0 abundances and between these and their equilibrium values are easily distinguished (it these cases the freeze-out temperatures for all three DMP components are very close).For smaller values, the effect of q is negligible, while for larger values the effects of q dominate the relic abundance and we find that 6 Experimental limits on model parameters

Constraints from the cold dark matter (CDM) relic density measurements
In this section we will obtain the numerical solution to the Boltzmann equations for the case q = 0, when 6 and find the region of parameter space that meets the relic-abundance constraint [44] As noted at the end of Sec.5 the solutions will depend on 3 independent parameters that we choose as M , f and λ h ; without loss of generality, we fix λ V to the value (2.25).We scan the 3-dimensional parameter space (M, f, λ h ) in the ranges 50 GeV ≤ M ≤ 2 TeV, 50 GeV ≤ f ≤ 1.5 TeV, 10 −4 ≤ |λ h | ≤ 1 for points allowed by (6.1); we also impose the constraint (2.27) and the one derived from h → ππ decay, which is open in the low M region (cf.Sec.3.4);note that in this region of parameter space the decay Z → πππ is kinematically forbidden, so that the restriction (3.26) does not apply.The q = 0 case is included by considering only the upper inequalities (see (5.31)).In the next section we consider the constraints direct-detection results from XENON100 and XENON1T experiments [38].In particular, using Y o + Y + + Y − |q| for q 10 −12 (cf. the end of Sec.5.2) we fin that (5.31) satisfies ( 6 In Fig. 10 we plot the relic abundance Ω DM h 2 and low-temperature distribution Y ∞ as functions of M .In Fig. 11 we show the region in the M − f plane allowed by the CDM constraint (6.1) as well as the region allowed by q = 0. Note, from the bottom panel of this figure, that Y ∞ cannot be assumed to be M independent as usually assumed in many models.In Fig. 11 we present the region in the M − f plane allowed by the CDM constraint We see from that figure that Ω DM increases with λ h : and the region of sufficiently small (large) λ h corresponds to an under (over)-abundance of DM.This is in contrast to models where the leading coupling to the DM fields is through the Higgs-portal interaction [31].We trace this difference to the presence of the ππ → Zπ interaction: comparing Fig. 6 and Fig. 11 we see that the region where the relic abundance is small (but still allowed by the data) corresponds to small values of λ h and also to σv ππ→SM (x = x f ) > σv ππ→πV (x = x f ); while large values of λ x correspond to the larger allowed values of the relic abundance and to σv ππ→SM (x = x f ) < σv ππ→πV (x = x f ).
The q = 0 allowed region in Fig. 11 can be approximated analytically by 39.65 We now use this result to extend the CDM limits with reasonable accuracy to the whole region of parameter space of interest.To do that note first that the s-wave contribution to σv ππ→SM is generated by the ππ → hh contribution (cf.Eq.(3.3)) so that in (5.15) σ SM ∼ (λ h M/f 2 ) 2 where the factor (|λ h |/f 2 ) 2 comes from the vertices, while the factor of M 2 is needed to get the right units (the other mass scales can be ignored for M > m h /2).
Similarly σ V ∼ (λ V M 2 /f 3 ) 2 where the factor (|λ V |/f 3 ) 2 comes from the vertices, while the factor of M 4 is needed to get the right units.Using this in (5.18) and (5.29) we find that up to a weak logarithmic dependence the parameters, 1/(h 2 Ω DM ) will depend on a linear combination of (λ h M/f 2 ) 2 and (λ V M 2 /f 3 ) 2 .Comparing then Fig. 6 and Fig. 11 we find that the upper limit in (6.3) corresponds to pa-rameters where σ SM dominates and where the upper limit in (6.1) is saturated; while the lower limit in (6.3) corresponds to parameters where σ V dominates and where the lower limit in (6.1) is saturated.Using this in conjunction with (6.3) we find that the CDM constrain reduces to 4.04 where δ q0 vanishes when q = 0 so that there is no upper limit in (6.4) in this case.

Direct detection constraints
The direct detection experiments probe the elastic scattering of DM particles off different kinds of materials [36][37][38].For the present model the leading interaction is the πN → πN scattering of DMP off the material's nucleons N (Fig. 12) through a t-channel Higgs exchange.The corresponding hard process was discussed in Sec.3.2where we show that the DMP-quark scattering cross section (3.13) is proportional to . Direct detection process.
The parton-level interaction is converted to the nucleon level by using effective nucleon f q N (N = p, n) couplings defined as [45] N |m q ψq ψ q |N = f q N M N , ( where M N is the nucleon mass and f p u = 0.0160, f p d = 0.0193, f p s = 0.0410, for the proton; f n u = 0.0108, f n d = 0.0284, f n s = 0.0409 for the neutron; while for the heavy quarks the f N q are generated by gluon exchange with the nucleon and are given by Then, DMP scattering with a nucleon composed of Z protons and A − Z neutrons is [45]  XENON1T is projected to exclude all points above the lower (red) solid line and would correspond to the constraint λ h /f 2 < 10 −6.5 .and the sum is over all quarks.The α q are effective couplings of DMP with the q-quarks, L = − 1 2 α q ψ q ψqππ that can be read off (3.14): Using microOMEGAs [45] we evaluate numerically the DMP-nucleon scattering cross section for direct detection and then compare these results to the XENON100 and XENON1T bounds.The results are presented in Fig. 13.As indicated above, if M is fixed the cross section depends only on λ h /f 2 and, in fact, the XENON bounds give rather simple expression for the constraints on this ratio: XENON100 : f 2 /λ h > 10 5.5 , XENON1T : f 2 /λ h > 10 6.5 . (6.9) The corresponding restrictions on the M − f plane over the CDM constrain are presented in Fig. 14.

Combined constraints on DMP model
The parameters in the model are constrained by the relations (2.29), (3.21), (6.4), and (6.9) that we collect here for convenience: Figure 14.Left: region in the M − f plane allowed by the CDM constraint and allowed (green) or disallowed (red) by the XENON100 data (6.9);black points are disallowed by (3.21).Right: same for the predicted XENON1T exclusion region in red and allowed in blue.We took q = 0, λ V = 0.63 and Higgs decay : where f, M are in GeV, and we used the XENON100 limit.The δ q,0 factor indicates that the corresponding limit disappears when non-zero values of q are allowed.
The resulting allowed regions in parameter space are given in Fig. 15 for our benchmark value of λ V = 0.63 as well as for the smaller natural value λ V = 0.0023 derived by NDA (see Sec.2.2).As can be seen from this figure if λ h 0 current data excludes DMP masses below ∼ 100 GeV while XENON1T would push this limit above 1 TeV.These limits do not apply when λ h 0; in this case low values (< 100 GeV) for M and f are allowed; in this case a non-zero value of q can always be found that meets all constraints (see Eq. 6.2).

Conclusions
We have studied a phenomenological model, where dark matter particles are pseudo-Goldstone bosons associated with the spontaneous breaking G DM → H DM ; we refer to these particles as dark matter "pions".The self-couplings and the couplings to the SM for such pionic DM differ from those of conventional scalars due to their chiral nature.We have illustrated the formalism for the case G DM = SU (2) × SU (2), H DM = SU (2) for which we have calculated all possible interactions and solved the Boltzmann equations to study the thermal history of when q = 0 for λ V = 0.0023.The various bands correspond to λ h = {0, 0.5, 1, 1.5, 2, 2.5, 3} from bottom to top, respectively; the darker regions correspond to those allowed by XENON1T.Top right panel: same for λ V = 0.63.Bottom panels: same as the top panels when q = 0. such pionic dark matter.We have also derived approximate analytic solutions and shown that they are consistent with the numerical calculations.
Our model of pionic dark matter satisfy relic abundance and direct detection constraint in a large region of parameter space.When the coupling to the Higgs is not too small the DMP mass M is required to lie above ∼ 100 GeV, and this lower limit will increase to ∼ 2 TeV if XENON1T does not detect a signal, since the absence of direct detection corresponds to relatively large values of f 2 /λ h .For each value of M the DMP decay constant f is moderately constrained to a range of values which is ∼ 200 GeV wide.
Collider signature of such dark matters at LHC is hard to see.The channel to study is essentially jets with missing energy [46], which is similar to many other dark matter model signatures [47].This requires a careful analysis to see if the existing bound in such channels put further constraints on the DMP parameter space, which lies beyond the scope of this paper.We will consider this in a future publication.
The DM couples to the SM via Z, γ and h, therefore it does not distinguish between fermion flavors.In particular there is no mechanism for suppressing the effects of the π at XENON experiments and enhancing them at DAMA/LIBRA [36].
As in QCD, there will presumably be baryons in this model (corresponding to solitons in the chiral theory, stabilized by higher derivative terms such as the Skyrme term [28]), but though they are SM singlets, they carry DM baryon number, so they do not couple singly to the SM, and they do not look like RH neutrinos.
A Effects on the Boltzmann equations of the SM particle decays to DMP.
The decay of the SM particles to the DMP require modification of the Boltzmann equation collision term by adding two terms C h and C Z corresponding to the h → ππ and Z → πππ decays.For the first, where the prefactor of N (i) h corresponds to the number π i produced, and we approximated (1 + f π ) 1. Since Γ does not depend on E h = p 2 + m 2 h , and using f h = e −E h /T (we assume a vanishing Higgs chemical potential), it follows with κ i defined in (3.4), Γ(h → ππ) is given in (3.20), and where we used (4.4).In complete analogy, the corresponding contribution from Γ → πππ is where Γ(Z → πππ) is given in (3.22).Note that for this decay the final state has a single π i (and a π ±α pair) so the prefactor corresponding to N (i)

B Kinetics of pure DMP
Using expressions from Sec.4 and Sec.5, and Eq.(3.18) the Boltzmann equations for pure DMP scattering are where dτ = ξ dx with and the last factor is explicitly given in (4.8).We solve these equations in two special cases • Suppose Y i = Y j = Y C for all i, j and Y α = Y β = Y R for all α, β; then with solutions where N is a constant and that are real only for 3δ 2 ≤ 1; note that the τ -dependent solutions interpolate between them.Only the second set has a region (|δ| ≤ 1/2) where they are all positive, so these correspond to the steady-state solutions.

which are presented in Fig. 2 .Figure 2 .
Figure 2. DMP scattering with Z and γ.

Figure 7 .Figure 8 .
Figure 7. Plot of Γ/H for the same parameters as in Fig.5.We also include the values of x f obtained from the condition ∆ = cY (eq) for c = 0.414, 0.732, 1 (left, center and right heavy dots on the dashed line, respectively).The freeze-out condition Γ = H corresponds to x f 31.3 which coincides almost exactly with the c = 1 value.

Figure 10 .
Figure 10.Ω DM h 2 (top left) and Y ∞ (top right) dependence on the DMP mass M for all values of f, λ h in the region scanned, and when q = 0 and λ V = 0.63.Red points: DM over-abundance (Ω DM h 2 > 0.13); blue points: region allowed by the CDM constraint (6.1); green points: DM underabundance (Ω DM h 2 < 0.094), which are allowed for appropriately chosen non zero q.The CDM-allowed region for Y ∞ is amplified in the bottom panel in order to better see the dependence on M .

Figure 13 .
Figure 13.Direct detection constraints from XENON experiments.XENON100 excludes all points above the solid line in purple at the top, which corresponds to the constraint λ h /f 2 < 10 −5.5 .XENON1T is projected to exclude all points above the lower (red) solid line and would correspond to the constraint λ h /f 2 < 10 −6.5 .

Figure 15 .
Figure 15.Top left panel: region in the f − M plane allowed by the combined constraints (6.10) when q = 0 for λ V = 0.0023.The various bands correspond to λ h = {0, 0.5, 1, 1.5, 2, 2.5, 3} from bottom to top, respectively; the darker regions correspond to those allowed by XENON1T.Top right panel: same for λ V = 0.63.Bottom panels: same as the top panels when q = 0.