1 Introduction

Multicomponent Dark Matter in Supersymmetric Hidden Sector Extensions

Daniel Feldman111e-mail: djfeld@umich.edu, Zuowei Liu222e-mail: liu@max2.physics.sunysb.edu, Pran Nath333e-mail: nath@neu.edu, and Gregory Peim444e-mail: peim.g@neu.edu footnotetext: Preprint Numbers: MCTP-10-15, YITP-SB-10-08, NUB-3266

Michigan Center for Theoretical Physics, Ann Arbor, Michigan 48104, USA

C.N. Yang Institute for Theoretical Physics, Stony Brook, New York, 11794, USA

Department of Physics, Northeastern University, Boston, Massachusetts 02115, USA

Abstract

Most analyses of dark matter within supersymmetry assume the entire cold dark matter arising only from weakly interacting neutralinos. We study a new class of models consisting of hidden sector extensions of the minimal supersymmetric standard model that includes several stable particles, both fermionic and bosonic, which can be interpreted as constituents of dark matter. In one such class of models, dark matter is made up of both a Majorana dark matter particle, i.e., a neutralino, and a Dirac fermion with the current relic density of dark matter as given by WMAP being composed of the relic density of the two species. These models can explain the PAMELA positron data and are consistent with the antiproton flux data, as well as the photon data from FERMI-LAT. Further, it is shown that such models can also simultaneously produce spin-independent cross sections which can be probed in CDMS-II, XENON-100 and other ongoing dark matter experiments. The implications of the models at the LHC and at the next linear collider (NLC) are also briefly discussed.

## 1 Introduction

Recently several particle physics models have been constructed that connect the standard model (SM) to hidden sectors and lead to massive narrow vector boson resonances as well as other signatures which can be detected at colliders [1, 2, 3]. The connection to the hidden sector arises via mass mixings and kinetic mixings[1, 2, 3, 4, 5, 6] and via higher dimensional operators. Models with the above forms of communication between the sectors also have important implications for dark matter [7, 3, 6] (for a review see [8, 9]). In this work we show that multicomponent dark matter can arise from extensions of the minimal supersymmetric standard model (MSSM) with Abelian hidden sectors which include hidden sector matter. Our motivation stems in part from the results of several dark matter experiments that have recently appeared. Thus the PAMELA Collaboration [10] has observed a positron excess improving previous results from HEAT and AMS experiments [11]. One possible explanation of such an excess is via the annihilation of dark matter in the galaxy[12]. Additionally, recent data from CDMS-II hints at the possibility of dark matter events above the background, and this will be explored further by the upgraded XENON experiment [13, 14].

For a thermal relic, the PAMELA data and CDMS-II data taken together at face value do raise a theoretical puzzle if indeed both signals arise from the annihilation of cold dark matter. Thus most models which aim to explain the PAMELA positron excess do not give a significant number of dark matter events in the direct detection experiments currently operating. Conversely, models which can give a detectable signal in direct detection experiments typically do not explain the PAMELA data without the use of enormous so-called boost factors. As we will show here, this can be circumvented in models where the dark matter has several components. Thus, motivated in part by the recent cosmic anomalies we develop supersymmetric models which contain minimally a hidden Abelian sector broken at the sub-TeV scale where the mass generation of the hidden states involves nontrivial mixings with the field content of the electroweak sector of the minimal supersymmetric extension of the standard model leading to dark matter which can have several components which can be both bosonic and fermionic.

More specifically, in this work we go beyond the simple theoretical construction that thermal dark matter compatible with WMAP observations is composed of a single fundamental particle. There is no overriding principle that requires such a restriction, and nonbaryonic dark matter (DM) may indeed be constituted of several components, so in general one has , where refers to the various species of dark particles that can contribute to the total nonbaryonic . In fact we already know that neutrinos do contribute to dark matter although their contribution is relatively small. Thus we propose here a new class of multicomponent cold dark matter models in Abelian extensions of MSSM which can simultaneously provide an explanation of the PAMELA and WMAP data through a Breit-Wigner enhancement [12], while producing detectable signals for the direct searches for dark matter with CDMS/XENON and other dark matter experiments.

A simultaneous satisfaction of the PAMELA positron excess and the satisfaction of WMAP relic density constraints can also occur if there is a nonthermal mechanism for the annihilation of dark matter with a wino lightest (R parity odd) supersymmetric particle (LSP) [15, 16, 17, 18, 8, 9]. However, a detectable spin-independent cross section in such a nonthermal framework does require that a pure wino is supplemented by a suitable admixture of Higgsino content as in the analysis of [19] and in [20], the later for a thermal relic. We remark that multiple factors and its influence on dark matter have very recently been studied [20, 21]. We also remark, some other works have recently looked at dark matter with more than 1 component [22]. The models proposed and analyzed here are very different from these.

The outline of the rest of the paper is as follows: In Sec.(2) we give a detailed description of the two models one of which is based on a extension of the MSSM where is a hidden sector gauge group with Dirac fermions in the hidden sector. This model allows for dark matter consisting of Dirac, Majorana, and spin zero particles. The second model is based on a extension of MSSM, where is a gauged leptophilic symmetry and , as before, is the hidden sector gauge group which also contains Dirac particles in the hidden sector. This model too has Dirac, Majorana, and spin zero particles as possible dark matter. In both cases we will primarily focus on the possibility that dark matter consists of Dirac and Majorana particles, and we will not discuss in detail the possibility of dark matter with bosonic degrees of freedom. In Sec.(3) we discuss the relic densities in the two component models. In Sec.(4) we give an analysis of the positron, antiproton, and photon fluxes in the two models. In Sec.(5) we give an analysis of event rates for the proposed models for CDMS-II and for XENON-100. We give the analysis within the framework of supergravity grand unified models [23, 24] defined by the parameters , and sign with nonuniversalities (NUSUGRA) defined by in the gaugino sector so that gaugino masses at the grand unified theory (GUT) scale are given by () (see, e.g., [25] and references therein). We also discuss the possible new physics one might observe at the LHC (for a recent review see also [9]) and elsewhere for these models. Conclusions are given in Sec.(7).

## 2 Multicomponent Hidden Sector Models

### 2.1 Multicomponent U(1)x model

A extension of the minimal supersymmetric standard model involves the coupling of a Stueckelberg chiral multiplet to vector supermultiplets , where is a real scalar and is an axionic pseudoscalar. Here is the vector multiplet which is neutral with respect to the SM gauge group with components and is the vector multiplet with components , where the components are written in the Wess-Zumino gauge. The chiral multiplet transforms under both and and acts as the connector sector between the visible and the hidden sectors. The total Lagrangian of the system is given by

 L=LMSSM+LU(1)X+LSt (1)

where is the kinetic energy piece for the vector multiplet and is the supersymmetric Stueckelberg mixing between the and the vector multiplets so that [1, 7] (see also [26, 27, 20])

 LSt=∫d2θd2¯θ (M1X+M2B+S+¯S)2 , (2)

where and are mass parameters. The Lagrangian of Eq.(1) is invariant under the and gauge transformations, i.e., under

 δXX=ζX+¯ζX ,δXS=−M1ζX,  δYB=ζY+¯ζY ,δYS=−M2ζY, (3)

where is an infinitesimal transformation chiral superfield. In component form we have for the Stueckelberg sector with

 LSt = −12(M1Xμ+M2Bμ+∂μσ)2−12(∂μρ)2−iχSσμ∂μ¯χS+2|FS|2 (4) +ρ(M1DX+M2DB)+¯χS(M1¯λX+M2¯λB)+χS(M1λX+M2λB) .

In addition, one may include a supersymmetric kinetic mixing term between the and gauge fields [7] leading to , where

 LU(1)X+LKM = −14XμνXμν−iλXσμ∂μ¯λX+12D2X (5) −δ2XμνBμν−iδ(λXσμ∂μ¯λB+λBσμ∂μ¯λX)+δDBDX .

One can also add additional terms as in [7]. Both Stueckelberg and kinetic mixings of the gauge fields and are constrained by the electroweak data[2]. As a consequence of the mixings, the extra gauge boson of the hidden sector couples with the standard model fermions and can become visible at colliders. The Lagrangian for matter interacting with the gauge fields is given by

 Lmatt = ∫d2θd2¯θ∑i[¯Φie2gYQYB+2gXQXXΦi+¯Φhid,ie2gYQYB+2gXQXXΦhid,i] . (6)

where the visible sector chiral superfields are denoted by (quarks, squarks, leptons, sleptons, Higgs, and Higgsinos of the MSSM) and the hidden sector chiral superfields are denoted by . In the above, is the hypercharge normalized so that . As mentioned already, the SM matter fields do not carry any charge under the hidden gauge group and vice versa, i.e. and . The minimal matter content of the hidden sector consists of a left chiral multiplet and a charge conjugate so that and have opposite charges and form an anomaly-free combination. A mass for the Dirac field arises from an additional term in the superpotential , where is composed of and . The scalar fields acquire soft masses of size from spontaneous breaking of supersymmetry by gravity mediation, and in addition acquire a mass from the term in the superpotential so that

 m2ϕ=m20+M2ψ=m2ϕ′. (7)

After spontaneous breaking of the electroweak symmetry there would be mixing between the vector fields , where is the third component of the field , . After diagonalization can be expressed in the terms of the mass eigenstates as follows:

 Vi=OijEi,  i,j=1−3,   E=(Z′,Z,γ). (8)

The neutral vector mass squared matrix is of the form given in Ref. [1] of [6]. Further, the chiral fermions in the multiplet together with the MSSM gauginos and Higgsinos will form a neutralino mass matrix whose eigenstates are six neutralino states , , where we assume that the set is the regular set of neutralinos and are the two additional neutralinos that arise in the extension. From the components and that appear in Eq.(4), we can form two Majorana fields and as follows:

 ΛX=(λXα¯λ˙αX), ψS=(χα,S¯χ˙αS). (9)

These components combine with the MSSM gauginos and Higgsinos to form a neutralino mass matrix whose eigenstates are the six neutralinos , . Thus and can be expanded as linear combination of , i.e.,

 ΛX=R1aχa,  a=1−6, ψS=R2aχa,  a=1−6 (10)

where is the unitary matrix which diagonalizes the neutralino mass matrix. Further the CP even Higgs sector is extended by the additional state [1]. The results outlined here give the following types of interactions:

1. There are interactions of the Dirac fermion in the hidden sector with the standard model particles via ,, interactions. Thus, the Dirac dark matter can annihilate into standard model particles via exchange of ,, in the early universe and in the galaxy. Depending on which of the two, Dirac or Majorana, is the heavier one may have Dirac particles annihilating into Majoranas or the Majorana particles annihilating into Dirac fermions in the galaxy:

 ¯ψψ→χχ  or  χχ→¯ψψ . (11)
2. In addition to the above we have fermion-neutralino-sfermion couplings in the hidden sector as given by Eq.(6). Thus interactions of the type , etc. can produce decays such as if they are kinematically allowed.

3. The scalar field is CP even and mixes with the MSSM Higgs fields. Through these mixings has couplings to the SM fermions and through these couplings it can decay into the SM fermions.

It is instructive to list all the new particles in this model as summarized below:

 New particles of the U(1)X model spin 0: ρ,ϕ,ϕ′, spin 12:ψ,χ05,χ06, spin 1:Z′. (12)

We assume that the lightest -parity odd particle (LSP) is the least massive neutralino () and resides in the visible sector and thus the masses of are larger than the LSP mass, and consequently are unstable and decay into SM particles and . The bosons and are unstable and decay into SM fermion pairs with the decay of the going dominantly through the process or if . The remaining three particles are all milli charged and, consequently, at least one of them is stable. If we assume , at least is always stable and the other two may or may not be stable. These along with the LSP give rise to various possible candidates for dark matter. Thus, depending on the relative masses of the Majorana, Dirac, and spin 0 particles there are three possibilities for the constituents of dark matter as outlined below.

Two component dark matter: Majorana + Dirac.– This model arises as follows: consider the case where . In this case the decays , will occur and will be unstable. Thus is stable and so is under the assumption of R parity conservation. Consequently, we will have two dark matter particles; namely, one a Majorana which is the LSP in the visible sector and the other a Dirac in the hidden sector. The Majorana and Dirac particles once created will annihilate as follows:

 ψ+¯ψ→Z,Z′,γ→SM+SM′, (13)
 χ+χ→(s:Z′,Z,h,H,A,ρ),(t/u:~fa,χi,χ±k)→SM+SM′. (14)

where and refer to and or channel exchanges. In addition to Eq.(14) there are coannihilation processes which contribute to the relic density. Since both and are stable, the total relic density of dark matter will be the sum of the relic densities for the two, the sum being constrained by the WMAP data. These constraints are discussed further in Sec.(3).

Three component dark matter: Dirac and two spin 0 particles.– Suppose the mass of is larger than the sum of the masses of the Dirac plus the scalar , i.e., . In this case the decay will occur and, consequently, is unstable. On the other hand, , and are stable since they cannot decay into anything else. Thus, here we have three dark matter particles: one Dirac, and the other two spin 0. Processes that lead to the annihilations of these particles are those in Eq.(13) for , and also for and , they are similar to those in Eq.(13), i.e., . In this three component dark matter model all the components reside in the hidden sector and thus their couplings to the standard model particles are extra weak. Consequently, they will have very small spin-independent cross sections in direct detection experiments. For this reason, this class of models is less preferred compared to the two component model.

Four component dark matter: Majorana, Dirac, and two spin 0 particles.– Finally, we consider the case when either of the following two situations occur: (i) , (ii) . In these cases all four particles, one Majorana, one Dirac, and two spin 0 particles, are stable and thus are possible dark matter candidates. These particles will annihilate to the SM particles as in Eq.(13), Eq.(14) and for and via processes in the three component dark matter model as described above. This model is in many ways similar to the two component model and like the two component model this model too should lead to detectable signals in experiments for the direct detection of dark matter.

### 2.2 Multicomponent Leptophilic U(1)X×U(1)C model

We discuss now another model which contains two additional Abelian vector bosons where one of the extra bosons is leptophilic. Leptophilic s have a long history [28] and have been revisited [29] over the recent past in the context of dark matter. Here we will consider a model where the as before is in the hidden sector, and is a leptophilic symmetry. As in the model, we also assume that the hidden sector has a pair of Dirac fermions and which are charged under but are neutral under the standard model gauge group and under . Regarding we assume it to be , i.e., a difference of family-lepton numbers, which is anomaly free, and can be gauged. The corresponding gauge field couples only to families and nothing else. The total Lagrangian in this case is

 L=LMSSM+LU(1)2+LSt, (15)

where is the kinetic energy for the and multiplets and for we assume the following form:

 LSt=∫d2θd2¯θ (M1C+M′2X+M′3B+S+¯S)2 +∫d2θd2¯θ (M′1C+M2X+M′′3B+S′+¯S′)2, (16)

where is the vector multiplet with components and and are the and multiplets as discussed before. The gauge transformations under , , and are

 δCC=ζC+¯ζC ,δCS=−M1ζC ,δCS′=−M′1ζC δXX=ζX+¯ζX ,δXS=−M′2ζX ,δXS′=−M2ζX , δYB=ζY+¯ζY ,δYS=−M′3ζY ,δYS′=−M′′3ζY , (17)

where , etc. are the infinitesimal transformation chiral superfields. The quantities and are the mass parameters. In the vector boson sector assumes the form

 LSt=−12(M1Cμ+M′2Xμ+M′3Bμ+∂μσ)2−12(M′1Cμ+M2Xμ+M′′3Bμ+∂μσ′)2. (18)

The mass matrix in the vector boson sector in the basis is given by

 ⎛⎜ ⎜ ⎜ ⎜ ⎜⎝M21+M′21M1M′2+M′1M2M1M′3+M′1M′′30M1M′2+M′1M2M22+M′22M′2M′3+M2M′′30M1M′3+M′1M′′3M′2M′3+M2M′′3M′23+M′′23+M2Y−MYMW00−MYMWM2W⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (19)

where is the boson mass and , and where is the weak angle. The dynamics of the model of Eq.(19) is rather involved. We will focus, therefore, on a simpler version of this more general case where we neglect the mixings with , i.e., we set . Inclusion of these coupling in the analysis would not drastically change the analysis or the conclusions of this work as long as we keep the mixing parameters very small. After neglecting the mixings with , the mass matrix is block diagonal and so we can diagonalize the top left hand corner mass matrix independent of the standard model sector. We are interested in the limit of small mixing between and and thus consider111Note these mass terms are different than those considered in Sec. 2.1

 M′1,M′2 ≪ M1,M2. (20)

In the above approximation the eigenvalues of this mass matrix are

 M2Z′≃M22+M′22−ΔM2,   M2Z′′≃M21+M′21+ΔM2, ΔM2≃(M1M′2+M′1M2)2(M21+M′21−M22−M′22). (21)

The corresponding mass eigenstates are and , where

 Cμ=cosθXZ′′μ−sinθXZ′μ, Xμ=sinθXZ′′μ+cosθXZ′μ, tanθX≃M1M′2+M′1M2M21+M′21−M22−M′22. (22)

Because of Eq.(20) . In the above, the Dirac fermions in the hidden sector have no couplings with the photon and are electrically neutral. However, by a small mixing of with in  Eq.(18), we can generate a milli charge for the Dirac particles in the hidden sector consistent with all electroweak data.

We discuss now the gaugino/chiral fermions in the extra sectors which arise from the superfields . From the gaugino components , , and from the chiral fermion components in the extra sectors , one can construct four component Majorana spinors two of which are exhibited in Eq.(9) and the remaining two are given by

 ΛC=(λCα¯λ˙αC),ψS′=(χα,S′¯χ˙αS′). (23)

The neutralino mass matrix in the model takes a block diagonal form

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣U(1)X×U(1)C04×4sector04×4MSSMsector⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦8×8 . (24)

Thus, the Stueckelberg mass generation produces a mass matrix in the hidden gaugino/chiral fermion sector which is decoupled from the neutralino mass matrix in the visible sector. Specifically in the 4 component notation the gaugino/chiral fermion mass matrix in the sector is given by

 LmassU(1)X×U(1)C=−⎛⎜ ⎜ ⎜ ⎜⎝¯ψS¯ψS′¯ΛC¯ΛX⎞⎟ ⎟ ⎟ ⎟⎠T⎛⎜ ⎜ ⎜ ⎜⎝00M1M′200M′1M2M1M′100M′2M200⎞⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜⎝ψSψS′ΛCΛX⎞⎟ ⎟ ⎟⎠. (25)

In the diagonalized basis we can label the extra neutralinos by . Since the hidden sector and the neutralinos of the visible sector are decoupled, the diagonalization of the neutralinos in the visible sector, i.e., of , is not affected. Further, as for the case of the model, it is instructive to list all the new particles in this model as summarized below:

 New particles of U(1)C×U(1)Xmodel spin 0: ρ,ρ′,ϕ,ϕ′, spin 12:ψ,χ05,χ06,χ07,χ08, spin 1:Z′,Z′′. (26)

We discuss now the stability of the new particles in this model. As before we assume that the mass of (and of ) is larger than the mass of . Thus will be stable since it cannot decay into anything. If kinematically allowed the fields and can decay only via the process as in the model. Of the remaining fields obviously and are unstable as they decay into as well as into depending on the mass of . As already noted, a small milli charge can develop for the hidden sector matter via small couplings of the and fields. The phenomenology of such models will be very similar to the one we are discussing here.

The extra neutralinos of Eq.(26) can also be all unstable. Thus couples with leptons-sleptons ( etc.) via coupling of the type , etc. and after diagonalization of the gaugino/chiral fermion mass matrix all the , will have coupling with leptons-sleptons of the type indicated. Further, two of the have roughly a mass of size while the remaining two have roughly a mass of size . Thus, if , which is what is assumed in this work, all the neutralinos of the hidden sector will be unstable and decay into final states of the type , etc. Regarding the field , there is an interaction of type

 M1gCρ(~f∗QfC~f),  f=e,μ . (27)

With this interaction will decay as follows: provided this process is kinematically allowed which we assume is the case. A similar situation occurs for the case of . Additionally, if there is a mixing with in the Stueckelberg sector then, as in the analysis of the model, the fields and will mix with the Higgs sector and can have decays of the type , , etc  . Thus, in the end we are left with a similar set of possibilities for dark matter as in the model, i.e., (i) a two component model with and , (ii) a three component model with , , , and (iii) a four component model with , , , and . However, as in the case we will focus on the two component model consisting of Dirac and Majorana dark particles.

We assume and that the annihilation of dark matter occurs close to the pole for reasons that will become apparent shortly. As a consequence, the annihilation of dark matter in the early universe and in the galaxy is controlled by the pole and the effect of the pole on the analysis is essentially negligible. The basic interaction of and of with matter is given by

 Lint=gXQX¯ψγμψXμ+gCQfC¯fγμfCμ (28)

where runs over and families and where . In the mass diagonal basis the interaction of Eq.(28) assumes the form

 Lint=(gXQX¯ψγμψcosθX−gCQfC¯fγμfsinθX)Z′μ +(gXQX¯ψγμψsinθX+gCQfC¯fγμfcosθX)Z′′μ. (29)

The interaction of Eq.(29) leads to the annihilation of into and via the poles for which we assume a Breit-Wigner form. Thus, the annihilation cross section takes the form

 σψ¯ψ→f¯f = aψ∣∣(s−M2Z′+iΓZ′MZ′)−1−(s−M2Z′′+iΓZ′′MZ′′)−1∣∣2, (30) aψ = βf(gXgCQXQfCsin(2θX))264πsβψ[s2(1+13β2fβ2ψ)+4M2ψ(s−2m2f)+4m2f(s+2M2ψ)],

where . The relevant partial decay widths are given by

 Γ(Z′→f¯f) = (gCQfCsinθX)2MZ′12π,  f=e,μ, (32) Γ(Z′→ψ¯ψ) = (33)

and similarly for the partial decay widths of the with and in Eq.(32) and in Eq.(33).

A constraint on comes from the contribution of the and to [30]. Their exchange gives

 Δ(gμ−2)=g2Cm2μ24π2[sin2θXM2Z′+cos2θXM2Z′′]. (34)

Using the current error [30] of in the determination of and assuming is small, one finds the following constraint on :

 αC ≲ 0.001(MZ′′300 GeV)2 , (35)

where . We note that if the mixing angle is small, the decay width of () and of will be narrow while the decay width of () and of will be of normal size. However, when the decay width into will also be small due to the kinematic suppression factor . In this case we will have the total width of the to be rather narrow. Thus for annihilation near the Breit-Wigner pole we will have a large enhancement of due to the narrowness of the [12]. It was shown in the analysis of Feldman-Liu-Nath in [12] that near the Breit-Wigner pole such annihilations allow one to fit the relic density as well as allow an enhancement of in the galaxy. We note that while decay width is very small this is not necessarily the case for which can decay into with normal strength. Thus neglecting the contribution of which is small due to the suppression, one finds the total width of to be . We will see in Sec.(4) that the needed in the analysis of the relic density is relatively small compared to normal electroweak coupling and, consequently, the width of though significantly larger than the width is still relatively small compared to what one might expect for a in a GUT model and certainly much smaller than the width for a arising as a Kaluza-Klein excitation in the compactification of an extra dimension [31, 32]. Finally, the annihilation of the Dirac particles in the early universe goes by the processes

 ψ¯ψ→Z′,Z′′→e+e−,μ+μ−,νe¯νe,νμ¯νμ, (36)

which is to be contrasted with the processes Eq.(13) in the model.

## 3 Relic Density in a Two Component Model

Here we discuss the relic density in models with two components. A general analysis requires solving the Boltzmann equations in a Friedmann-Robertson-Walker universe [33, 34], and includes coannihilations [35] and an accurate integration over pole regions. As in the MSSM alone, one will generally encounter the and Higgs poles [36] and these need to be treated with care. The number changing processes include

 ψ¯ψ↔SM SM′,  ψ¯ψ↔χχ,  χχ↔SM SM′. (37)

Note that the process is not allowed since connect only to and , neither of which can connect to the standard model particles. For the simplest two component model with dark matter particles with the assumption that the only relevant processes in the annihilation of are final states. Since is heavier than its freeze-out occurs earlier (at a higher ) than for . Thus, the Boltzmann equations for (which includes fermions and antifermions) and for for the and for the two component models are given by

 dnψdt=−3Hnψ−12⟨σv⟩ψ¯ψ(n2ψ−n2ψ,eq), (38)
 dnχdt=−3Hnχ−⟨σv⟩χχ(n2χ−n2χ,eq)+12⟨σv⟩ψ¯ψ→χχ(n2ψ−n2ψ,eq). (39)

Here refers to , and stands for . For the spin averaged cross section for the Dirac case, the extra factor of is to account for the fact that we are dealing with a Dirac fermion. The number densities are and are their values at equilibrium, i.e., , where and . Since the two dark matter particles are sub-TeV in mass, they will freeze-out at temperatures that are not drastically different. One can solve the Boltzmann equation for with the appropriate boundary conditions to compute the freeze-out temperature and the relic density of at the current temperatures. To compute the freeze-out temperature for the particles , one uses solutions for as computed from the Boltzmann equation for as input in the Boltzmann equation for keeping in mind that