# Cold Dark Matter Production In Early Universe

Guided by:

## Abstract

For many decades, astronomical observations have hinted at the presence of non-luminous matter, the Dark Matter (DM). From galactic to cosmological scales, the presence of dark matter is felt through its gravitational effects. There is no suitable candidate for DM within the framework of the standard model (SM) of particle physics; neither it finds an explanation in astrophysical compact objects.

Recent cosmic microwave background raditaion (CMBR) observations infer that DM accounts for approximately 85% of matter in the universe and about a quarter of its energy density. We assume that DM consists of weakly interacting massive particles. This report consists of a review of different cold dark matter production mechanisms which can produce the observed relic density.

## Abbreviations

DM | Dark Matter |

CMBR | Cosmic Microwave Background Radiation |

SM | Standard Model of particle physics |

## Evidence/Background Of Dark Matter

Before we start discussing dark matter, let us first see why there is a need for dark matter in the first place. As already pointed out, there have been many observations that do not co-incide with our understanding of the Universe. Let us look at a few of them in order of their length scales. (Joseph Silk, 2005)

## The Galactic Scale

The most convincing and direct evidence of Dark matter on galactic scale comes from the observation of the rotation curves of disc galaxies, which is the graph of circular velocities of stars and gas as a function of thier distance from galactic centre. In newtonian dynamics, the circular velocity is expected to be:

$\frac{mv^2}{r} = \frac{GMm}{r^2}$

or $v = \sqrt{\frac{GM(r)}{r}}$

where as usual $M(r)=4\pi\int\rho(r)r^2dr$and $\rho(r)$is the mass density profile and should be falling as $\frac{1}{\sqrt{r}}$.

Thus we estimate the mass inside a radial distance by: $M(r) = \frac{v^2R}{G}$.

This is a very powerful result which is capable of telling us important information about mass distributions in galaxies, provided we have spherical symmetry. The observed rotation curve for a disc galaxy is shown in the figure below.

We see that v(r) is approximately constant, which implies the existence of unseen mass with

$M(r) \propto r and \rho(r) \propto \frac{1}{r^2}$. Rotation curve measurements of low surface brightness galaxies indicate an extremely high mass-to-light ratio, meaning that stars and luminous gas contribute only very little to the overall mass balance of LSB.

## The Galactic Clusters Scale

A particularly important evidence for the existence of dark matter is the discrepency between the luminous mass and the total mass of a galaxy cluster. The total mass of a galaxy cluster can be measured in three different ways :

1. From the scatter in radial velocities of the galaxies within clusters. Zwicky (1933) noticed that the outer members of the Coma cluster were moving far too quickly to be merely tracing the gravitational potential of the visible cluster mass. The only way the observed velocities of the cluster members could be reconciled with the Virial theorem was to postulate that the cluster also contained another large, but unseen, mass component: dark matter.

2. From X-rays emitted by hot gas in clusters.

3. From gravitational lensing.

We now know that most of the atoms in the clusters of galaxies are not seen in observations with visible light. Because these clusters generate enormously deep potential wells, it is easy for hydrogen gas to leak out and fill the whole volume of the cluster. These atoms acquire large velocities and emit X-rays. The X-ray spectrography of the clusters, however, does not remove the mystery. It can only account for at the most 20% of the mass of the cluster and cannot explain the origin of the deep potential well. For this, we must postulate that the clusters are also filled with a new invisible, weakly interacting form of matter.

Bullet Clusters provide yet another strong evidence for the existence of dark matter. Bullet clusters are those clusters which have recently collided. In the case of bullet clusters, lensing map exhibits a large amount of matter that has passed straight through the gas clouds, and appears undisturbed after the collision. This shows that the dark matter does not have to necessarily track luminous matter and that it does not interact strongly with gas or itself; this means it is effectively collisionless.

## Big Bang Nucleosynthesis

A critical prediction of the Hot Big Bang cosmology is that the protons and the neutrons were fused in the primordial reballs to create the light elements as it cooled to temperatures of the order of an MeV. The resulting elemental abundances depend only on the nuclear reaction rates and the baryon-to-photon ratio ( $\eta$) at that time. Laboratory and theoretically calculated nuclear reaction rates can therefore be used to derive the primordial abundances of the elements as a function of $\eta$.

The abundances of D and the He isotopes point very consistently to a baryon density of $\Omega_B = 0.04$, where $\Omega = \frac{\rho_m}{\rho_c}$, $\rho_m$ is present mass density and $\rho_c$ is critical density. However, recent surveys from large-scale galaxy structures indicates $\Omega_m$ = 0:29. The independent measurements of $\Omega_B$ from Big Bang Nucleosynthesis and $\Omega_m$ from large scale structures together provide evidence for the existence of dark matter. Since all luminous matter essentially consists of baryons, $\Omega_m$ = 0:29 and $\Omega_m$ = 0:04 together imply that the remaining $\Omega$(leftover) = 0:25 must be dark matter. Thus we are now able to quantify the amount of dark matter that must be present. (Edward Kolb, 2018)

## Cosmic Microwave Background

The CMB is a form of electromagnetic radiation which is a remnant from the early stage of universe, and is also known as relic radiation. Very precise observations have been carried out to measure the CMB spectrum. It is known to be highly isotropic and only shows anisotropy at subdegree levels. An extraordinary wealth of information lies in the cosmic microwave background anisotropy.

Each fluctuation rang in this way with the frequency dependent on its size. The dark matter, which does not interact with photons, remained unaffected by this ringing effect. Areas into which matter had fallen were hotter than the surroundings. Areas from which matter had streamed out, by contrast, were cooler. The CMB measurements have now verified that Baryons account for about four percent of the critical density. Further, they suggest that the matter density is some ten times higher than this, implying the existence of non-baryonic dark matter and dark energy. (Scott Dodelson, 2002)

## Identifying Dark Matter

## Can dark baryons fit the bill?

Dark baryons are objects that do not shine. One class of dark baryons is MACHO or Massive Astrophysical Compact Halo Object, such as black holes, neutron stars, black dwarfs. The recent black hole image is a concrete proof that galaxies contain a super-massive black hole. Can these MACHOs be our unknown dark matter? The answer is NO. Although MACHO's have been detected through gravitational microlensing observations, they are insuffcient for explaining all of the galaxy's dark mass. Thus gravitational microlensing observations have mostly ruled out the possibility of MACHO's as the explanation of dark matter. BBN also limits the cosmological density of baryons to a fraction of the matter density, and it is hard to imagine that the missing baryons could pop out late.

Within the standard model of particle physics, neutrino appears to be a promising DM candidate. A relic thermal bath of neutrinos, produced alongside the cosmic microwave background, could yield the required mass density if the neutrino mass was 10eV. However, numerical simulations have shown the top-down structure formation in a neutrino dominated universe, where superclusters form first. Hence neutrinos cannot fit the bill. (Stefano Profumo, 2015)

## Can it be modified gravity?

An even better argument is given by the observed "smoothness" of the CMB sky. It is almost impossible to explain the timely formation of non-linear structures in the universe unless baryons falls into pre existing gravitational potential wells seeded by dark matter. In truth, certain Tensor-Vector-Scalar or TeVeS theories do predict large enough inhomogeneities for structures to form in time. However, such models vastly overpredict baryon accoustic oscillations and the matter power spectrum they predict is typically in contradiction with data.

## Weakly Interacting Massive Particles (WIMPs)

Weakly interacting massive particles are hypothetical particles that are thought to constitute dark matter. A WIMP is an elementary particle which interacts via gravity and any other force, potentially not part of the standard model itself, which is as weak or weaker than the weak nuclear force, but also non-vanishing in its strength. A WIMP must have been produced thermally in the early universe, similarly to the particles of standard model according to standard cosmology, and usually will constitute dark matter. The mass of a WIMP is expected to be in the range of 100GeV. Because supersymmetric extensions of the standard model of physics predict a new particle with these properties, the apparent coincidence is known as "WIMP MIRACLE". (Griest 2006)

The identity of the WIMP remains a mystery. The earliest idea was that the WIMP was a heavy fourth-generation Dirac or Majorana neutrino. Since then, numerous other candidates have been put forth. Of all the candidates, perhaps the most well-motivated and certainly the most theoretically well developed WIMP candidate is the lightest super-symmetric particle(LSP). Vast experiments are now being conducted to detect WIMPs. The most promising techniques involve direct-detection in low-background laboratory detectors and indirect detection through observation of energetic neutrinos from annihilations of WIMP that have accumulated in the sun and/or the earth.

## Thermal Production Of WIMPs In Early Universe

## Thermal Relics

Thermal decoupling consists of the process where a particle species at high temperatures, participates in a reaction that keeps it in "thermal equilibrium" in the statistical mechanical sense. As the universe cools, the particles' number density decreases, and eventually the rate $\Gamma$ for the reaction becomes smaller than the expansion rate of the universe at that temperature. $\Gamma$ << H(T) indicates that the reaction keeping the particle species in the universe is too slow: The particle species is no longer in thermal equilibrium. This is called thermal decoupling, and the species is said to have frozen out of thermal decoupling.

For much of its early history, most of the constituents of the universe were in thermal equilibrium. However, there have been a number of notable departures from thermal equilibrium. The departures from equilibrium have led to important relics - the light elements, the cosmic microwave background, relic WIMPs etc. Once a species totally decouples from the plasma, its evolution is very simple: particle number density decreasing as $R^{-3}$ and particle momenta decreasing as $R^{-1}$. The evolution of the phase space distribution of a species which is in LTE or is completely decoupled is simple. It is the evolution of particle distributions that is challenging. The rough criterion for a particle species to be either coupled or decoupled involves the comparison of interaction rates with the expansion rate

$\Gamma$> H (Coupled)

$\Gamma$ < H (Decoupled)

The relic density of any species can be calculated using the Boltzmann equation.

## Boltzmann Equation

One of the most powerful tools in cosmology is the Boltzmann equation. It gives us the power to calculate the relic density of species, which were thermally produced in the early universe. If we consider species to be stable, only annihilation and inverse annihilation processes ( $\psi \bar{\psi} \leftrightarrow \chi \bar{\chi}$) can change the number of $\psi$ and $\bar{\psi}$ in a comoving volume. The Boltzmann equation in this case can be written as:

$\frac{dn_\psi}{dt} + 3Hn_\psi = -\langle\sigma v\rangle(n^2 - n_{eq}^2)$

where $n_\psi$ is actual number density of $\psi s$.

$\sigma$ is the cross-section of the reaction that keeps the species in equilibrium.

$n_{eq}$ is the equilibrium number density of $\psi s$.

The $3H{n}_{\psi}$ term accounts for the dilution effect of the expansion of the universe. $<\sigma v>$accounts for the interactions that change the number of $\psi s$ present. In the absence of interaction, ${n}_{\psi}$ decreases with expansion of the universe and goes as ${n}_{\psi}\propto {R}^{-3}$.

It is usually useful to scale out the effect of the universe by considering the evolution of the number of particles in a comoving volume. This greatly simplifies the equation and is done by using the entropy density s (which is considered to be constant). Defining a new variable

$Y=\frac{{n}_{\psi}}{s}$ and using the conservation of entropy per comoving volume, it follows that

$\frac{d{n}_{\psi}}{dt}+3H{n}_{\psi}=s\frac{dY}{dt}$

Furthermore, since the interaction term will usually depend explicitly on temperature rather than time, it is useful to introduce a new independent variable $x=m/T$. The Boltzmann equation in terms of our new defined variables is

$\frac{dY}{dx}=\frac{-x<\sigma v>s(Y^2-Y_{eq}^2)}{ H(m)}$

where $H(m)=1.67g_\ast^\frac12\frac{m^2}{m_{pl}}$.

It is important to realize the physics behind the equation. If $\Gamma >H$, the RHS term will dominate. Interaction being the dominant term, there will occur enough reactions to keep the species in equilibrium. The actual number density will then follow equilibrium number density, which is given by Boltzmann distribution function. When the interaction term becomes less than the Hubble parameter, LHS dominates. The actual number density no more follows the equilibrium number density and the species fall out of equilibrium, the number density per co-moving volume becomes constant.

## Cosmological Abundnace Of WIMP

The Boltzmann equation for the evolution of the abundance of a species is a particular form of Riccati equation, for which there are no general, closed form solutions. We therefore computationally solve the equation to calculate the relic density of WIMPs. The recent observations infer dark matter density to be $\Omega {h}^{2}=0.120$. (Perdereau et al 2016). For WIMPS of mass $\simeq 100GeV$, we can solve the boltzmann equation to find the reaction cross-section that will produce the observed relic density.

Solving Boltzmann equation numerically is a notoriously difficult thing to do because in the interesting region, the Boltzmann equation depends on the small difference between two big numbers. Rushing naively into solving the equation can lead to numerical instabilities and grumpy code. We therefore make a substitution: $y=\lambda Y$. The Boltzmann equation becomes:

$\frac{dy}{dx}=-{x}^{-2}({y}^{2}-{y}_{eq}^{2})$

where $y_{eq}= 0.192\frac{ g}{\sqrt{g*}}\frac{ m}{\sqrt{8\pi G}}\langle\sigma v\rangle x^{1.5}e^{-x}$

The mathematica code which solves our first order non-homogeneous differential equation for different values of m and $<\sigma v>$ is given in the Appendix at the end of this report The results are shown in the graph below. At small values of x, Y follows Yeq and the graph is exponentially decreasing. At larger values of x, Y becomes constant indicating the freeze out of species.

*In calculatidng the relic density of WIMPs, the value of all the constants have been taken in GeV units considering c=1 and $\hbar$=1. The Dark Matter is considered as a Majorana fermion with degrees of freedom g=2.*

Mass (GeV) | Reaction Cross-section ( $c{m}^{3}{s}^{-1}$) |

1 | $4.8\times {10}^{-26}$ |

100 | $2\times {10}^{-26}$ |

1000 | $2.12\times {10}^{-26}$ |

## Analytical Approximation

The Boltzmann equation can be solved analytically using certain approximations. In order to do so, we must first parameterize the temperature dependence of the annihilation cross-section. Considering annihilation cross-section has velocity dependence $\langle\sigma v\rangle \propto v^t$ where t = 0 for s-wave annihilation and t = 2 corresponds to p-wave annihilation. Since $\langle v\rangle \propto T^{\frac{1}{2}}$, $\langle\sigma\;v\rangle\propto T^n$, where n = 0 is for s-wave annihilation. Therefore we have

$\langle\sigma v\rangle=\sigma_0(\frac Tm)^n=\sigma_0x^{-n}$

With this parameterization, the Boltzmann equation for the abudance of $\psi s$ become

$\frac{dY}{dx} = -\lambda x^{-n-2}(Y^2 - Y_{eq}^2)$

where, $\lambda = [\frac{x <\sigma v>s}{H(m)}] = 0.264(\frac{g_{*s}}{g_*^{\frac{1}{2}}})m m_{pl}\sigma_0$

To analytically solve it further, we consider the departure from equilibrium $\Delta = Y - Y_{eq}$. The differemtial equation becomes

$\frac{d(\Delta + Y_{eq})}{dx} = -\lambda x^{-n-2}(Y-Y_{eq})(Y+Y_{eq})$

$\frac{d\Delta}{dx}=-\lambda x^{-n-2}\Delta(\Delta+Y_{eq}+Y_{eq})-Y'_{eq}$

$\frac{d\Delta}{dx}=-\lambda x^{-n-2}\Delta(\Delta+2Y_{eq})-Y'_{eq}$

Consider $x_f$ the point of freeze out of species. For 1 < x < $x_f$ , Y tracks Yeq very closely, and both $\Delta$ and $|\Delta'|$ are small. An appropriate solution is obtained by setting $\Delta'$ = 0.

$-\lambda x^{-n-2} \Delta (\Delta +2Y_{eq})-Y'_{eq} = 0$

$-\lambda x^{-n-2} \Delta (\Delta +2Y_{eq})=Y'_{eq}$

$\Delta = \frac{x^{n+2}Y'_{eq}}{2\lambda Y_{eq}}$

Since Y tracks $Y_{eq}$ closely, we can approximate $Y'_{eq}$ = $Y_{eq}$. Our equation becomes:

$\Delta = \frac{x^{n+2}}{2\lambda}$

At late times x >> $x_f$, Y tracks $Y_{eq}$very poorly : $\Delta \simeq Y >> Y_{eq}$, and terms involving $Y'_{eq}$ and $Y_{eq}$can be safely neglected.

$\Delta' = -\lambda x^{-n-2}\Delta (\Delta +2Y_{eq})$

$\Delta' = -\lambda x^{-n-2}\Delta^2$

To calculate the relic density of species, we must first calculate Y at time when the species fall out of equilibrium. Thus we integrate the above equation from x = $x_f$to x = $\infty$. The above equation is an example of Cauchy-Euler differential equation, and can be solved by setting z = $\frac{1}{\Delta}$.

$\frac{dz}{dx} = -\frac{1}{\Delta^2}\frac{d\Delta}{dx}$

$\frac{dz}{dx} = -\lambda x^{-n-2}$

$\int dz = \lambda \int_{x_f}^\infty x^{-n-2}dx$

$z=\lambda(\frac{x^{-n-1}}{-n-1})$

$Y_\infty = \Delta_\infty = \frac{n+1}{\lambda}x_f^{n+1}$

Now our job is to determine $x_f$. x = $x_f$is the time when Y ceases to track $Y_{eq}$. At this moment $\Delta$becomes of the order $Y_{eq}$and we can assume $\Delta(x_f)$= c $Y_{eq}$, where c is a constant. The early time solution becomes

$\Delta(x_f) = \frac{x_f^{n+2}Y'_{eq}}{\lambda(2Y_{eq}+cY_{eq})}$

$\Delta(x_f) \simeq \frac{x_f^{n+2}}{\lambda(2+c)}$

Substituting in previous equation and solving we get

$x_f = ln[(2+c)\lambda ac]-(n+\frac{1}{2})ln[ln[(2+c)\lambda ac]]$

Choosing c(c + 2) = n(n + 1) gives best fit to the numerical results for the final abundance $Y_\infty$. With this choice we obtain

$x_f = ln[0.038(n+1)\frac{g}{\sqrt{g_*}}m_{pl}m\sigma_0]-(n+\frac{1}{2})ln[ln[0.038(n+1)\frac{g}{\sqrt{g_*}}m_{pl}m\sigma_0]]$

## Freeze Out of Dark Matter

We will use the above result to calculate the temperature at which our m = 100GeV candidate froze out. Taking

m=100

mpl = $1.2211\times 10^{19}$

sv = $1.7\times 10^{-9}$

we find that $x_f$= 19.62. Numerical calculations in previous section gives $x_f\simeq 20$. This is a very precise result which gives

$T_f=\frac{m}{x_f}$

$T_f=5GeV = 5.8025\times 10^{13} K$

Once we know the temperature of the particle, we can estimate its velocity. Using $v=\sqrt{\frac{8k_bT}{\pi m}}$, we find that the velocity of dark matter particles at the time of freeze out should be equal to

$v=1.07\times10^{8}=0.3c$

Boltzmann equation can be numerically solved to calculate the relic abundance of WIMPs today. We use the same trick that we used for symmetric dark matter to scale out the numerical instabilities. Multiplying equation with, $\lambda$λ, we get

## Asymmetric Dark Matter

Up till now, we considered the case of the universe that exhibits dark matter symmetry i.e. there are equal number of dark matter particles and anti-particles. We know today that our universe exhibits baryon asymmetry. The ordinary matter in our universe is almost entirely made up of baryons, with anti-baryons contributing only a tiny fraction. It is thus natural to consider asymmetric dark matter, for which particles and anti-particles are not identical.

For symmetric dark matter, we have produced observed relic density considering WIMP annihilation cross-section is of weak size. In this section we generalize this calculation to the case of asymmetric dark matter. We make the following assumptions:

1. WIMPs can annihilate only with their anti-particles.

2. The number of particles minus the number of anti-particles in a co-moving volume is conserved.

3. Dark matter asymmetry is created well before dark matter annihilation reactions freeze out. (Iminniyaz et al 2011).

## Boltzmann Equation For Asymmetric Dark Matter Abundance

Consider a dark matter particle denoted by $\chi$ that is not self-conjugate, i.e , the anti-particle $\chi$ is not equal to $\bar{\chi}$. The relic abundance of $\chi$ and $\bar{\chi}$ are determined by solving the boltzmann equation. Under the assumption that only $\chi$ $\bar{\chi}$ pairs can annihilate into standard model particles, while $\chi$ $\chi$ and $\bar{\chi}$ $\bar{\chi}$ cannot, the relevant Boltzmann equation are,

$\frac{dn_\chi}{dt} +3Hn_\chi= -\langle\sigma v\rangle(n_\chi n_{\bar{\chi}} - n_{eq,\chi}n_{eq,\bar{\chi}})$

$\frac{dn_{\bar{\chi}}}{dt} +3Hn_{\bar{\chi}}= -\langle\sigma v\rangle(n_\chi n_{\bar{\chi}} - n_{eq,\chi}n_{eq,\bar{\chi}})$

We assume that the dark matter particles were already non-relativistic at decoupling. The equilibrium number densities are given by

$n_{\chi,eq} = g_\chi(\frac{m_\chi T}{2\pi})^\frac{3}{2}e^{-\frac{(m_\chi + \mu_\chi)}{T}}$

$n_{\bar\chi,eq} = g_\chi(\frac{m_\chi T}{2\pi})^\frac{3}{2}e^{-\frac{(m_\chi - \mu_\chi)}{T}}$

where $m_\chi$is the mass of the WIMP. $\mu_\chi$ is the chemical potential of the particles, and $g_\chi$ is the number of internal degrees of freedom. We have used the fact that $\mu_{\bar\chi}=-\mu_{\chi}$ in the equilibrium. We follow the standard picture of dark matter evolution and rewrite Boltzmann equation in terms of dimensionless quantities, $Y_\chi = \frac{n_\chi}{s}$ and $x=\frac{m}{T}$.

$\frac{dY_\chi}{dx} = -\frac{\lambda<\sigma v>}{x^2}(Y_\chi Y_{\bar\chi}-Y_{\chi,eq}Y_{\bar\chi,eq})$

$\frac{dY_{\bar\chi}}{dx} = -\frac{\lambda<\sigma v>}{x^2}(Y_\chi Y_{\bar\chi}-Y_{\bar\chi,eq}Y_{\bar\chi,eq})$

where $\lambda=\frac{4\pi}{\sqrt{90}}m_\chi m_{pl}\sqrt{g_*}$ . Subtracting second from first, we obtain

$\frac{dY_\chi}{dx} - \frac{dY_{\bar\chi}}{dx}=0$

This implies

$Y_\chi - Y_{\bar\chi} = C$

where C is a constant. This follows from our assumption that $\chi$ and $\bar\chi$ only annihilate with each other. Our Boltzmann equations become

$\frac{dY_\chi}{dx} = -\frac{\lambda<\sigma v>}{x^2}(Y_\chi Y_{\bar\chi}-CY_{\chi}-P)$

$\frac{dY_{\bar\chi}}{dx} = -\frac{\lambda<\sigma v>}{x^2}(Y_\chi Y_{\bar\chi}+CY_{\bar\chi}-P)$

where $P=Y_{\chi,eq}Y_{\bar\chi,eq}=(0.145\frac{g_\chi}{g_*})^2x^3e^{-2x}$

*We have considered the difference in co-moving densities of particles and

antiparticles to be of order of $10^{11}$. Annihilation from only s-wave is taken into

consideration.*

Boltzmann equation can be numerically solved to calculate the relic abundance of WIMPs today. We use the same trick that we used for symmetric dark matter to scale out the numerical instabilities. Multiplying equation with, $\lambda$λ, we get

$\frac{dy}{dx} = -\frac{<\sigma v>}{x^2}(y^2 - \lambda Cy-P )$

where $y=\lambda Y$and $P=\frac{g_*}{90}(4\pi m_\chi m_{pl}g_* 0.145\frac{g_\chi}{g_*})^2x^3e^{-2x}$

We take help from mathematica to solve our differential equation numerically. The code can be found in the appendix at the end of this report. The results are shown in the graph below. Unlike symmetric DM, the relic abundance of asymmetric DM depends on the initial asymmetry.

## Dark Matter Coannihilation

The cases of symmetric and asymmetric dark matter have been discussed. One possibility is that there was a whole sector of particles $\chi_1, \chi_2, ... \chi_k$ , where k can tend to infinity. The relic particle is the lightest of this set and also our DM candidate. In this case the relic abundance of the lightest particle is determined not only by its annihilation cross section, but also by the annihilation of the heavier particles, which will later decay into the lightest. We call this the case of coannihilation.(Griest & Seckel 1991).

## Assumptions::

(1) We assume a whole sector of particles $\chi_1, \chi_2, ... \chi_k$. $\chi_1$ is the lightest particle, and our DM candidate, and all the other particles eventually decay into $\chi_1$.

(2) All the particles were in thermal equilibrium with the standard model particles.

(3) We consider all possible annihilations between particles, that kept them in equilirium.

(4) The particles masses are nearly degenerate.

## Boltzmann Equation For Dark Matter Coannihilation

Consider the annihilation of N dark matter particles $\chi_i(i=1 ... N)$ with masses $m_i$ and internal degrees of freedom $g_i$. The evolution of the number density $n_i$ of particle i is (Gondolo & Edsjö 1999).

$\frac{dn_i}{dt} = -3Hn_i - \sum_{j=1}^N\langle\sigma v\rangle_{ij}(n_in_j - n_{i,eq}n_{j,eq}) -\sum_{j not equal i}\{\Gamma_{ij}(n_i-n_{i,eq})-\Gamma_{ji}(n_j-n_{j,eq})\}$

The first term on the right hand side is the dilution due to the expansion of the universe, H is the Hubble parameter.

The second term describes $\chi_i\chi_j$ annihilations.

The third term accounts for $\chi_i$ deacys.

The final abundance of the lightest particle is simply described by the sum of the density of all supersymmetric particles,

$n=\sum_{i=1}^N n_i$

Summing equations of all $n_i$, we get evolution equation for n

$\frac{dn}{dt} = -3Hn-\sum_{i,j=1}^N \langle\sigma v\rangle_{ij}(n_in_j-n_{i,eq}n_{j,eq})$

where the decay terms cancel in the sum. The $\chi_i$ distributions remain in thermal equilibrium, and in particular their ratios are equal to the equilibrium values,

$\frac{n_i}{n}\simeq \frac{n_{i,eq}}{n_{eq}}$

We then get

$\frac{dn}{dt} = -3Hn-\langle\sigma_{eff}v\rangle(n^2-n_{eq}^2)$

where $\langle\sigma_{eff}v\rangle=\sum_{ij}(\sigma v)_{ij}\frac{n_{i,eq}}{n_{eq}}\frac{n_{j,eq}}{n_{eq}}$.

This equation is the same as the Boltzmann equation for a single WIMP, except the WIMP annihilation rate is replaced by the total effective annihilation rate.

## A Simple Case Of Two Particles Coannihilation

Let us assume that the dark sector comprises two particle species $\chi_1$ and $\chi_2$. $\chi_1$ is the lighter particle whose relic density is the topic of interest, and $\chi_2$ eventually decays into $\chi_1$. The mass difference beween $\chi_1$ and $\chi_2$ is assumed to be small, $\Delta m << m_1$, where $m_1$ is the mass of species $\chi_1$. We evaluate the final Boltzmann equation for this case and study the effect of $\langle\sigma_{eff}v\rangle$ on the relic abundace of $\chi_1$.

The first particle Boltzmann equation can be written as

$\frac{dn_{\chi_1}}{dt}=-3Hn_{\chi_1}-\langle\sigma v\rangle_{12}(n_{\chi_1}n_{\chi_2}-n_{\chi_1,eq}n_{\chi_2,eq})-\langle\sigma v\rangle_{11}(n_{\chi_1}^2-n_{\chi_1,eq}^2)+\Gamma_{21}(n_{\chi_2}-n_{\chi_2,eq})$

For the second particle we have

$\frac{dn_{\chi_2}}{dt}=-3Hn_{\chi_2}-\langle\sigma v\rangle_{12}(n_{\chi_1}n_{\chi_2}-n_{\chi_1,eq}n_{\chi_2,eq})-\langle\sigma v\rangle_{22}(n_{\chi_2}^2-n_{\chi_2,eq}^2)-\Gamma_{21}(n_{\chi_2}-n_{\chi_2,eq})$

Adding the above two equations, we get

$\frac{dn}{dt}=-3Hn-2\langle\sigma v\rangle_{12}(n_{\chi_1}n_{\chi_2}-n_{\chi_1,eq}n_{\chi_2,eq})-\langle\sigma v\rangle_{11}(n_{\chi_1}^2-n_{\chi_1,eq}^2)-\langle\sigma v\rangle_{22}(n_{\chi_2}^2-n_{\chi_2,eq}^2)$

where $n=n_1+n_2$.

We have $\langle\sigma_{eff}\rangle = \sum_{i,j}^N\sigma_{ij}r_ir_J$

Here $r_i=\frac{n_i^{eq}}{n^{eq}} = \frac{g_i(1+\Delta_i)^{\frac{3}{2}}exp(-x\Delta_i)}{g_{eff}}$ and $\Delta_i=\frac{m_i-m_1}{m_1}$

$\langle\sigma_{eff}\rangle$ for the case of two particles hidden sector becomes

$\sigma_{eff}=2\sigma_{12}\frac{g_1g_2}{g^2_{eff}}(1+\Delta_1)^\frac{3}{2}(1+\Delta_2)^\frac{3}{2}(-x(\Delta_1+\Delta_2))+\sigma_{11}\frac{g_1^2}{g_{eff}^2}(1+\Delta_1)^2exp(-2x\Delta_1)+\sigma_{22}\frac{g_2^2}{g_{eff}^2}(1+\Delta_2)^2exp(-2x\Delta_2)$

where $\Delta_1=\frac{m_1-m_1}{m_1}=0$ and $\Delta_2=\frac{m_2-m_1}{m_1}=\frac{\Delta m}{m_1}$

Thus we obtain

$\sigma_{eff}=2\sigma_{12}\frac{g_1g_2}{g^2_{eff}}(1+\frac{\Delta m}{m_1})^\frac{3}{2}exp(-x\frac{\Delta m}{m_1})+\sigma_{11}\frac{g_1^2}{g_{eff}^2}+\sigma_{22}\frac{g_2^2}{g_{eff}^2}(1+\frac{\Delta m}{m_1})^2exp(-2x\frac{\Delta m}{m_1})$

Two cases are possible: first when $\sigma_{11}$ is greater than $\sigma_{22}$ and the second when $\sigma_{22}$ is greater than $\sigma_{11}$. In each case we calculate $\sigma_{eff}$. Due to the prescence of exponential in $\sigma_{22}$ term, $\sigma_{eff}$ in the second case becomes larger than the first one. Since $\sigma_{eff}$ increaes in second case, less relic abundance is expected.

## Dark Matter Cannibalization

A new type of dark matter can be considered, where the number changing processes keeps it in equilibrium. During this era the dark matter cannibalizes its rest mass to keep warm. The self-interactions of the dark matter are assumed to be strong at the time when dark matter becomes non-relativistic, but the couplings to ordinary matter are assumed feeble, the dark matter will have a temperature T' that is different from photon temperature T.

We consider the freeze out via N-to-N' reactions. As DM only annihilate within the dark sector, it is also refeered to as 'dark freeze-out'. This scenario can give rise to both large DM self-interactions and the proper DM relic abundance. We consider for now $4\rightarrow 2$ annihilating processes. The Boltzmann equation which describes the evolution of the DM number density n'=n'(T') is (Bernal and Xiaoyong 2016):

$\frac{dn'}{dt}+3H(T)n' = -\langle\sigma v^3\rangle_{4-2}[n'^4-n'^2n'^2_{eq}]$

We will make a distinction between thermal and chemical equilibrium. When we speak of thermal equilibrium, what we mean is that the dark matter has a well defined temperature, and the occupation number of momenum states is determined only in terms of chemical potential and temperature. Thermal equilibrium can be maintained by elastic collisions between pairs of dark matter particles, which presumably are much faster than $4\rightarrow 2$ processes. in contrast, chemical equilibrium means that number changing processes are also fast, which assumes that the chemical potential is zero, and consequently the occupation number can be written in terms of temperature only. (Carlson et al 1992)

## Numerical Solution Of Boltzmann Equation

The evolution of DM number density is given by the boltzmann equation:

$\frac{dn'}{dt}+3H(T)n' = -\langle\sigma v^3\rangle_{4-2}[n'^4-n'^2n'^2_{eq}]$

which in case where the SM energy density dominates the expansion of the universe can be rewritten in terms of the SM temperature T as

$xH(x)\frac{dY}{dx}=-s(x)^3Y^2\langle\sigma v^3\rangle_{4-2}[Y^2-Y_{eq}^2]$

where $x=\frac{m}{T}$ and $Y=\frac{n'}{s(x)}$.

Assuming different temperatures in both sectors, we solve the boltzmann equation for temperature ratio $\frac{T'}{T}=\frac{1}{4}$. We observe that the observed relic abundance can be obtained as long as we assume a temperaure ratio T'/T less than one. Since HS (Hidden Sector) remains colder than SM, the equilibrium DM abundnace is smaller than in the case where both sectors thermalize. During the freeze out, it is thus not necessary to deplete as much DM as before, so a smaller cross-section leads to the observed relic abundance.

We get

$Y_{eq}(x)=\frac{2\pi^2}{45(2\pi)^{\frac{3}{2}}}\frac{(mT')^{\frac{3}{2}}}{T^3}\frac{g}{g_{*s}}exp(-m/T')$

$Y_{eq}(x)=\frac{2\pi^2}{45(2\pi)^{\frac{3}{2}}}(\frac{m}{T})^{\frac{3}{2}}(\frac{T'}{T})^{\frac{3}{2}}\frac{g}{g_{*s}}exp(-\frac{m}{T}\frac{T}{T'})$

$Y_{eq}(x)=0.018x^{\frac{3}{2}}\frac{g}{g_{*s}}exp(-4x)$

Once again we take the help of Mathematica to solve the equation numerically for different values of interaction rate. The code which is given in the appendix produces the graph given below. The reaction cross-section rate that produces the correct relic abundance is not very high, and can sustain the DM galaxy halos present today.

For the temperature ratio $\frac{T'}{T}=\frac{1}{4}$, the correct relic abundance is obtained for $<\sigma v^3>=2 GeV^{-4}$.

## Freeze-In Production Of Dark Matter

There is an alternate mechanism for dark matter production, "freeze-in". We consider a long lived particle $\chi$, having interactions with the bath that are so feeble that $\chi$ is thermally decoupled from plasma. At very high temperatures we assume a negligible initial $\chi$ abundance. As the universe evolves, $\chi$ particles are produced from collisions or decay of bath particles, for instance by $\sigma \rightarrow\chi \chi$, where $\sigma$ is a particle in the visible sector heat bath. The freeze-in yield is active until the number density of $\sigma$ particles becomes boltzmann suppressed. The comoving number density of DM particles $\chi$ then becomes a constant and the DM abundnace freezes-in. (Hall et al 2010).

Freeze-out and Freez-in share a crucial common aspect: the final relic abundance, given the relevant particle masses and couplings, can be computed solely from an initial state of bath particles that are in thermal equilibrium. A complete understanding of the history of universe is not required in both the cases. Freeze-in can be viewed as the opposite of freeze-out.

## Freeze-In Calculations

The freeze-in process is dominated by decays or inverse-decays of bath particles to the FIMP (Feebly Interacting Massive Particles). Processes suvch as $2 \rightarrow 2$ and other scatterings give sub-dominant contribution to the frozen-in abundance. If the initial DM number density is small compared to the equilibrium number density, the backward reaction term in the Boltzmann equation can be neglected, and the evolution of $n_\chi$ is given by

$\frac{dn_\chi}{dt}+3Hn_\chi=2\Gamma_{\sigma \rightarrow \chi \chi}\frac{K_1(m_\sigma/T)}{K_2(m_\sigma/T)}n_\sigma^{eq}$

where K's are modified Bessel functions of second kind, $\Gamma_{\sigma \rightarrow \chi \chi}$ is the decy width and $n_\sigma^{eq}$ is the equilibrium number density of $\sigma$. Defining then $Y=n_\chi/s$ and $x=m_\sigma/T$ and assuming again that the number of relativistic degrees of freedom remain constant during DM production, the Boltzmann equation can be rewritten as

$\frac{x}{Y_\sigma^{eq}}\frac{dY}{dx}=2\frac{\Gamma_{\sigma \rightarrow \chi \chi}}{H}\frac{K_1(x)}{K_2(x)}$ (Bernal et al 2017)

where $Y_\sigma^{eq}=0.145\frac{g}{g_{*s}}x^{\frac{3}{2}}e^{-x}$ and $\Gamma=n_\chi\langle\sigma v\rangle$.

In the radiation dominated era, H is given by: $H^2=\frac{4\pi^3}{45m_{pl}^2}g_{eff}T^4$. The equation becomes

$\frac{dY}{dx}=0.368\frac{g}{g_*^{\frac{1}{2}}}m_\sigma m_{pl}\langle\sigma v\rangle x^{-\frac{1}{2}}e^{-x}Y\frac{K_1(x)}{K_2(x)}$

Once again we take the help of mathematica to solve the equation.

The result of our calculation is shown in the graph below:

Freeze-in unlike freeze-out, depends on the the intial condition. A change in the initial conditions will reflect directly on the relic abundance of the FIMPs. If the interaction rate is too high, the species can re-annihilate before freezing-in, and we obtain once again the observed relic abundance. The freeze-in production mechanism offers a variety of cases that could have produced the relic abundance. The one we have discussed is the simplest one.

## Conclusion

There is considerable evidence to assume that DM exists, but what actually it is remains an open mystery. Despite decades of search for DM, its interaction with ordinary matter, beyond gravity, remains unknown. Lots of DM models have been devised to generate observed abundance, few of which have been reviewed in this report. Different DM models are consistent with different DM properties, and our understanding of the early universe. These DM properties are what we look for in the collider experiments.

While WIMP remains our favourite DM model, a courtesy of WIMP miracle, there is no denying that DM could have been produced by other mechanisms. Future experiments might give us an insight into the correct production mechanism. Till then we hope for another miracle to show up.

## APPENDIX

## MATHEMATICA CODES:

## 1. Thermal Production Of WIMPs

$m=100$

$mpl=1.2211*{10}^{19}$

$sv=1.75*{10}^{-9}$

$sol1=NDSolve\left[\right\{y\text{'}\left[x\right]==-{x}^{-2}*\left(y\right[x{]}^{2}-(0.192*0.2*mpl*m*sv*{x}^{1.5}*{e}^{-x}{)}^{2}),$

$y\left[1\right]==0.192*0.2*mpl*m*sv*{e}^{-1}\},y,\{x,1,20\left\}\right]$

$bc=Evaluate\left[y\right[20]/.sol1]$

$sol2=NDSolve\left[\right\{y\text{'}\left[x\right]==-{x}^{-2}*\left(y\right[x{]}^{2}-(0.192*0.2*mpl*m*sv*{x}^{1.5}*{e}^{-x}{)}^{2}),$

$y\left[20\right]==bc\},y,\{x,20,10000\left\}\right]$

$LogLogPlot\left[Evaluate\right[y\left[x\right]/.sol1],\{x,1,20\},PlotRange->\{\{1,10000\},\{1,{10}^{12}\}\left\}\right]]$

$LogLogPlot\left[Evaluate\right[y\left[x\right]/.sol2],\{x,20,10000\},PlotRange->\{\{1,10000\},\{1,{10}^{12}\}\left\}\right]$

$LogLogPlot[0.192*0.2*m*mpl*sv*{x}^{1.5}*{e}^{-x},\{x,1,1000\},PlotRange->\{\{1,10000\},\{1,{10}^{12}\}\},PlotStyle->RGBColor[0.93,0.4,0.33\left]\right]$

$Show[\%,\%\%,\%\%\%]$

## 2. Asymmetric Dark Matter

$m=100$

$m_{pl}=2.44\times10^{18}$

$sv=1.5\times10^{-10}$

$k=10^{-11}$

$g=2$

$g_s=90$

$lam=4*\pi*m*m_{pl}*sqrt[gs]/sqrt[90]$

$sol1 = NDSolve[\{y'[x]==-sv*x^{-2}*(y[x]^2-(k/lam)*y[x]-(0.145*4*\pi*m*m_{pl}*(g/gs))^2*x^3*e^{-2x}),y[1]==0.145*\pi*m*m_{pl}*(g/gs)*1^{1.5}*e^{-1}\},y,\{x,1,20\}]$

$sol3 = NDSolve[\{y'[x]==-sv*x^{-2}*(y[x]^2+(k/lam)*y[x]-(0.145*4*\pi*m*m_{pl}*(g/gs))^2*x^3*e^{-2x}),y[1]==0.145*\pi*m*m_{pl}*(g/gs)*1^{1.5}*e^{-1.1}-(lam*k)\},y,\{x,1,20\}]$

$cd=Evaluate[y[20]/.sol3]$

$sol4 = NDSolve[\{y'[x]==-sv*x^{-2}*(y[x]^2+(k/lam)*y[x]-(0.145*4*\pi*m*m_{pl}*(g/gs))^2*x^3*e^{-2x}),y[20]==cd\},y,\{x,20,10000\}]$

$LogLogPlot[Evaluate[y[x]/.sol1]/lam,\{x,1,20\},PlotRange->\{\{1,10000\},\{1,10^{-12}\}]$

$LogLogPlot[Evaluate[y[x]/.sol2]/lam,\{x,20,10000\},PlotRange->\{\{1,10000\},\{1,10^{-12}\}]$

$LogLogPlot[Evaluate[y[x]/.sol3]/lam,\{x,1,20\},PlotRange->\{\{1,10000\},\{1,10^{-12}\},PlotStyle->Dashed]$