Loading...

Summer Research Fellowship Programme of India's Science Academies

Cold Dark Matter Production In Early Universe

Badal Bhalla

Indian Institute Of Technology Jodhpur, Rajasthan

Guided by:

Arindam Chatterjee

Physics and Applied Mathematics Unit, Indian Statistical Institute, Kolkata 700108.

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.

Keywords: Dark Matter, Cosmology, Standard Model Of Particle Physics, Boltzmann Equation, CMB

Abbreviations

 DMDark 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 fi rst 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:

mv2r=GMmr2\frac{mv^2}{r} = \frac{GMm}{r^2}

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

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

Thus we estimate the mass inside a radial distance by: M(r)=v2RGM(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 fi gure below.

rotation1.jpg
    Rotation curve of Disc Galaxy (Milky Way). This figure is adopted from the site: www.ircamera.as.arizona.edu/

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

    M(r)r and  ρ(r)1r2M(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 ΩB=0.04\Omega_B = 0.04, ​​where  Ω=ρmρc \Omega = \frac{\rho_m}{\rho_c}​, ρm\rho_m is present mass density and ρc\rho_c is critical density. However, recent surveys from large-scale galaxy structures indicates Ωm\Omega_m = 0:29. The independent measurements of ΩB\Omega_B from Big Bang Nucleosynthesis and Ωm\Omega_m from large scale structures together provide evidence for the existence of dark matter. Since all luminous matter essentially consists of baryons, Ωm\Omega_m = 0:29 and Ωm\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.

    It so happens that the CMB isotropy can only be explained if dark matter is taken into account when matter and radiation are decoupled. Before the neutral hydrogen formed, the matter was distributed almost uniformly in space - although small variations occured in the densities of both normal and dark matter owing to quantum mechanical fluctuations. Gravity pulled the normal and dark matter in toward the center of each fluctuation. While the dark matter continued to move inward, the normal matter fell in only until the pressure of photons dragged it back (matter drag radiation with itself through electro-magnetic interactions), causing it to flow outward until the gravitational pressure overcame the photon pressure, and the matter began to fall in once more.

    CMB2018_Planck_4672.jpg
      CMB Temperature Fluctuaions (​​Akrami et al 2019​​) 

      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 verifi ed 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?

      Are we really sure that the explanation to these observations is dark matter? Can it possibly be modifi ed gravity after all? The bullet cluster system's X-ray versus gravitational potential map shows that the baryons do not sit where most of the gravitationally interacting matter, as traced by weak gravitational lensing. It is naturally very hard to explain observations of this type with modi fied 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 R3R^{-3} and particle momenta decreasing as R1R^{-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:

       dnψdt+3Hnψ=σv(n2neq2) \frac{dn_\psi}{dt} + 3Hn_\psi = -\langle\sigma v\rangle(n^2 - n_{eq}^2)

      where nψn_\psi is actual number density of ψs\psi s.

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

      neqn_{eq} is the equilibrium number density of ψs\psi s.

      The 3Hnψ term accounts for the dilution effect of the expansion of the universe. <σv>accounts for the interactions that change the number of ψs present. In the absence of interaction, nψ decreases with expansion of the universe and goes as nψ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=nψs and using the conservation of entropy per comoving volume, it follows that

      dnψdt+3Hnψ=sdYdt

      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

      dYdx=x&lt;σv&gt;s(Y2Yeq2)H(m)\frac{dY}{dx}=\frac{-x&lt;\sigma v&gt;s(Y^2-Y_{eq}^2)}{ H(m)}

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

      It is important to realize the physics behind the equation. If Γ>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 Ωh2=0.120.​​​ (​Perdereau et al 2016​​). For WIMPS of mass 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=λY. The Boltzmann equation becomes:

      dydx=-x-2(y2-yeq2)

      where yeq=0.192ggm8πGσvx1.5exy_{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 <σ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.*

      DM_1.png
        Variation of number density per co-moving volume with x.
        Reaction cross-section corresponding to different values of m
         Mass (GeV) Reaction Cross-section ( cm3s-1)  
        1  4.8×10-26
         100 2×10-26
         1000 2.12×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 σvvt\langle\sigma v\rangle \propto v^t where t = 0 for s-wave annihilation and t = 2 corresponds to p-wave annihilation. Since  vT12\langle v\rangle \propto T^{\frac{1}{2}}, σ&ThickSpace;vTn\langle\sigma\;v\rangle\propto T^n, where n = 0 is for s-wave annihilation. Therefore we have

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

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

        dYdx=λxn2(Y2Yeq2)\frac{dY}{dx} = -\lambda x^{-n-2}(Y^2 - Y_{eq}^2)

        where, λ=[x&lt;σv&gt;sH(m)]=0.264(gsg12)mmplσ0\lambda = [\frac{x &lt;\sigma v&gt;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 Δ=YYeq\Delta = Y - Y_{eq}. The differemtial equation becomes

        d(Δ+Yeq)dx=λxn2(YYeq)(Y+Yeq) \frac{d(\Delta + Y_{eq})}{dx} = -\lambda x^{-n-2}(Y-Y_{eq})(Y+Y_{eq}) 

        dΔdx=λxn2Δ(Δ+Yeq+Yeq)Yeq\frac{d\Delta}{dx}=-\lambda x^{-n-2}\Delta(\Delta+Y_{eq}+Y_{eq})-Y&#x27;_{eq}

        dΔdx=λxn2Δ(Δ+2Yeq)Yeq\frac{d\Delta}{dx}=-\lambda x^{-n-2}\Delta(\Delta+2Y_{eq})-Y&#x27;_{eq}

        Consider xfx_f the point of freeze out of species. For 1 < x < xfx_f , Y tracks Yeq very closely, and both Δ\Delta and Δ|\Delta&#x27;| are small. An appropriate solution is obtained by setting Δ\Delta&#x27; = 0.

        λxn2Δ(Δ+2Yeq)Yeq=0-\lambda x^{-n-2} \Delta (\Delta +2Y_{eq})-Y&#x27;_{eq} = 0

        λxn2Δ(Δ+2Yeq)=Yeq-\lambda x^{-n-2} \Delta (\Delta +2Y_{eq})=Y&#x27;_{eq}

        Δ=xn+2Yeq2λYeq\Delta = \frac{x^{n+2}Y&#x27;_{eq}}{2\lambda Y_{eq}}

        Since Y tracks YeqY_{eq} closely, we can approximate YeqY&#x27;_{eq} = YeqY_{eq}. Our equation becomes:

        Δ=xn+22λ\Delta = \frac{x^{n+2}}{2\lambda}

        At late times x >> xfx_f, Y tracks YeqY_{eq}very poorly : ΔY&gt;&gt;Yeq\Delta \simeq Y &gt;&gt; Y_{eq}, and terms involving YeqY&#x27;_{eq} and YeqY_{eq}can be safely neglected.

        Δ=λxn2Δ(Δ+2Yeq)\Delta&#x27; = -\lambda x^{-n-2}\Delta (\Delta +2Y_{eq})

        Δ=λxn2Δ2\Delta&#x27; = -\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 = xfx_fto x = \infty. The above equation is an example of Cauchy-Euler differential equation, and can be solved by setting z = 1Δ\frac{1}{\Delta}.

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

        dzdx=λxn2\frac{dz}{dx} = -\lambda x^{-n-2}

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

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

        Y=Δ=n+1λxfn+1Y_\infty = \Delta_\infty = \frac{n+1}{\lambda}x_f^{n+1}

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

        Δ(xf)=xfn+2Yeqλ(2Yeq+cYeq)\Delta(x_f) = \frac{x_f^{n+2}Y&#x27;_{eq}}{\lambda(2Y_{eq}+cY_{eq})}

        Δ(xf)xfn+2λ(2+c)\Delta(x_f) \simeq \frac{x_f^{n+2}}{\lambda(2+c)}

        Substituting in previous equation and solving we get

        xf=ln[(2+c)λac](n+12)ln[ln[(2+c)λac]]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 fi nal abundance YY_\infty. With this choice we obtain

        xf=ln[0.038(n+1)ggmplmσ0](n+12)ln[ln[0.038(n+1)ggmplmσ0]]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×10191.2211\times 10^{19}
        sv = 1.7×1091.7\times 10^{-9}​​

        we fi nd that xfx_f= 19.62. Numerical calculations in previous section gives xf20x_f\simeq 20. This is a very precise result which gives

        Tf=mxfT_f=\frac{m}{x_f}

        Tf=5GeV=5.8025×1013KT_f=5GeV = 5.8025\times 10^{13} K

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

        v=1.07×108=0.3cv=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,

        dnχdt+3Hnχ=σv(nχnχˉneq,χneq,χˉ) \frac{dn_\chi}{dt} +3Hn_\chi= -\langle\sigma v\rangle(n_\chi n_{\bar{\chi}} - n_{eq,\chi}n_{eq,\bar{\chi}}) 

        dnχˉdt+3Hnχˉ=σv(nχnχˉneq,χneq,χˉ) \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χ,eq=gχ(mχT2π)32e(mχ+μχ)Tn_{\chi,eq} = g_\chi(\frac{m_\chi T}{2\pi})^\frac{3}{2}e^{-\frac{(m_\chi + \mu_\chi)}{T}}

        nχˉ,eq=gχ(mχT2π)32e(mχμχ)Tn_{\bar\chi,eq} = g_\chi(\frac{m_\chi T}{2\pi})^\frac{3}{2}e^{-\frac{(m_\chi - \mu_\chi)}{T}}

        where mχm_\chiis the mass of the WIMP. μχ\mu_\chi is the chemical potential of the particles, and gχ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χ=nχsY_\chi = \frac{n_\chi}{s} and x=mTx=\frac{m}{T}.

        dYχdx=λ&lt;σv&gt;x2(YχYχˉYχ,eqYχˉ,eq)\frac{dY_\chi}{dx} = -\frac{\lambda&lt;\sigma v&gt;}{x^2}(Y_\chi Y_{\bar\chi}-Y_{\chi,eq}Y_{\bar\chi,eq})

        dYχˉdx=λ&lt;σv&gt;x2(YχYχˉYχˉ,eqYχˉ,eq)\frac{dY_{\bar\chi}}{dx} = -\frac{\lambda&lt;\sigma v&gt;}{x^2}(Y_\chi Y_{\bar\chi}-Y_{\bar\chi,eq}Y_{\bar\chi,eq})

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

        dYχdxdYχˉdx=0\frac{dY_\chi}{dx} - \frac{dY_{\bar\chi}}{dx}=0

        This implies

        YχYχˉ=CY_\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

        dYχdx=λ&lt;σv&gt;x2(YχYχˉCYχP)\frac{dY_\chi}{dx} = -\frac{\lambda&lt;\sigma v&gt;}{x^2}(Y_\chi Y_{\bar\chi}-CY_{\chi}-P)

        dYχˉdx=λ&lt;σv&gt;x2(YχYχˉ+CYχˉP)\frac{dY_{\bar\chi}}{dx} = -\frac{\lambda&lt;\sigma v&gt;}{x^2}(Y_\chi Y_{\bar\chi}+CY_{\bar\chi}-P)

        where P=Yχ,eqYχˉ,eq=(0.145gχg)2x3e2xP=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 101110^{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

        dydx=&lt;σv&gt;x2(y2λCyP)\frac{dy}{dx} = -\frac{&lt;\sigma v&gt;}{x^2}(y^2 - \lambda Cy-P )

        where y=λYy=\lambda Yand P=g90(4πmχmplg0.145gχg)2x3e2xP=\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.

        ADM.png
          Variation of number density per co-moving volume with x for asymmetric dark
          matter. It is evident that in case of asymmetric dark matter, freeze out occurs at  xfx_f = 22. 
          ADM2.png
            Difference in number densities of χ\chi and χˉ\bar{\chi} after freeze out. Dashed line depicts the
             evolution of number density of χˉ\bar{\chi}.

            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 χ1,χ2,... χk\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 χ1,χ2,... χk\chi_1, \chi_2, ...  \chi_k. χ1\chi_1 is the lightest particle, and our DM candidate, and all the other particles eventually decay into χ1\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 χi(i=1 ... N)\chi_i(i=1  ...  N) with masses mim_i and internal degrees of freedom gig_i. The evolution of the number density nin_i of particle i is ​​(​Gondolo & Edsjö 1999​​).

            dnidt=3Hnij=1Nσvij(ninjni,eqnj,eq)j not equal i{Γij(nini,eq)Γji(njnj,eq)}\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 χiχj\chi_i\chi_j annihilations.

            The third term accounts for χi\chi_i deacys.

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

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

            Summing equations of all nin_i, we get evolution equation for n

            dndt=3Hni,j=1Nσvij(ninjni,eqnj,eq)\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 χi\chi_i distributions remain in thermal equilibrium, and in particular their ratios are equal to the equilibrium values,

            ninni,eqneq\frac{n_i}{n}\simeq \frac{n_{i,eq}}{n_{eq}}

            We then get

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

            where σeffv=ij(σv)ijni,eqneqnj,eqneq\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 χ1 \chi_1 and χ2\chi_2. ​ χ1 \chi_1 is the lighter particle whose relic density is the topic of interest, and χ2\chi_2 eventually decays into χ1 \chi_1. The mass difference beween χ1\chi_1 and χ2\chi_2 is assumed to be small, Δm&lt;&lt;m1\Delta m &lt;&lt; m_1, where m1m_1 is the mass of species χ1 \chi_1. We evaluate the final Boltzmann equation for this case and study the effect of σeffv\langle\sigma_{eff}v\rangle on the relic abundace of χ1 \chi_1.

            The first particle Boltzmann equation can be written as

            dnχ1dt=3Hnχ1σv12(nχ1nχ2nχ1,eqnχ2,eq)σv11(nχ12nχ1,eq2)+Γ21(nχ2nχ2,eq)\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

            dnχ2dt=3Hnχ2σv12(nχ1nχ2nχ1,eqnχ2,eq)σv22(nχ22nχ2,eq2)Γ21(nχ2nχ2,eq)\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​

            dndt=3Hn2σv12(nχ1nχ2nχ1,eqnχ2,eq)σv11(nχ12nχ1,eq2)σv22(nχ22nχ2,eq2)\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=n1+n2n=n_1+n_2.​

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

            Here ri=nieqneq=gi(1+Δi)32exp(xΔi)geffr_i=\frac{n_i^{eq}}{n^{eq}} = \frac{g_i(1+\Delta_i)^{\frac{3}{2}}exp(-x\Delta_i)}{g_{eff}} and Δi=mim1m1\Delta_i=\frac{m_i-m_1}{m_1}

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

            σeff=2σ12g1g2geff2(1+Δ1)32(1+Δ2)32(x(Δ1+Δ2))+σ11g12geff2(1+Δ1)2exp(2xΔ1)+σ22g22geff2(1+Δ2)2exp(2xΔ2)\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 Δ1=m1m1m1=0\Delta_1=\frac{m_1-m_1}{m_1}=0 and Δ2=m2m1m1=Δmm1\Delta_2=\frac{m_2-m_1}{m_1}=\frac{\Delta m}{m_1}

            Thus we obtain

            σeff=2σ12g1g2geff2(1+Δmm1)32exp(xΔmm1)+σ11g12geff2+σ22g22geff2(1+Δmm1)2exp(2xΔmm1)\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 σ11\sigma_{11} is greater than σ22\sigma_{22} and the second when σ22\sigma_{22} is greater than σ11\sigma_{11}. In each case we calculate σeff\sigma_{eff}. Due to the prescence of exponential in σ22\sigma_{22} term, σeff\sigma_{eff} in the second case becomes larger than the first one. Since σeff\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 424\rightarrow 2 annihilating processes. The Boltzmann equation which describes the evolution of the DM number density n'=n'(T') is ​(​Bernal and Xiaoyong 2016​​):

            dndt+3H(T)n=σv342[n4n2neq2]\frac{dn&#x27;}{dt}+3H(T)n&#x27; = -\langle\sigma v^3\rangle_{4-2}[n&#x27;^4-n&#x27;^2n&#x27;^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 424\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:

            dndt+3H(T)n=σv342[n4n2neq2]\frac{dn&#x27;}{dt}+3H(T)n&#x27; = -\langle\sigma v^3\rangle_{4-2}[n&#x27;^4-n&#x27;^2n&#x27;^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)dYdx=s(x)3Y2σv342[Y2Yeq2]xH(x)\frac{dY}{dx}=-s(x)^3Y^2\langle\sigma v^3\rangle_{4-2}[Y^2-Y_{eq}^2]

            where x=mTx=\frac{m}{T} and Y=ns(x)Y=\frac{n&#x27;}{s(x)}​.

            Assuming different temperatures in both sectors, we solve the boltzmann equation for temperature ratio TT=14\frac{T&#x27;}{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

            Yeq(x)=2π245(2π)32(mT)32T3ggsexp(m/T)Y_{eq}(x)=\frac{2\pi^2}{45(2\pi)^{\frac{3}{2}}}\frac{(mT&#x27;)^{\frac{3}{2}}}{T^3}\frac{g}{g_{*s}}exp(-m/T&#x27;)

            Yeq(x)=2π245(2π)32(mT)32(TT)32ggsexp(mTTT)Y_{eq}(x)=\frac{2\pi^2}{45(2\pi)^{\frac{3}{2}}}(\frac{m}{T})^{\frac{3}{2}}(\frac{T&#x27;}{T})^{\frac{3}{2}}\frac{g}{g_{*s}}exp(-\frac{m}{T}\frac{T}{T&#x27;})

            Yeq(x)=0.018x32ggsexp(4x)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 TT=14\frac{T&#x27;}{T}=\frac{1}{4}, the correct relic abundance is obtained for &lt;σv3&gt;=2 GeV4&lt;\sigma v^3&gt;=2  GeV^{-4}.

            Cannibilization.png
              Evolution of actual number density in cannibilization regime.

              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 222 \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χn_\chi is given by

              dnχdt+3Hnχ=2ΓσχχK1(mσ/T)K2(mσ/T)nσeq\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σeqn_\sigma^{eq} is the equilibrium number density of σ\sigma. Defining then Y=nχ/sY=n_\chi/s and x=mσ/Tx=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

              xYσeqdYdx=2ΓσχχHK1(x)K2(x)\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σeq=0.145ggsx32exY_\sigma^{eq}=0.145\frac{g}{g_{*s}}x^{\frac{3}{2}}e^{-x} and Γ=nχσv\Gamma=n_\chi\langle\sigma v\rangle.

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

              dYdx=0.368gg12mσmplσvx12exYK1(x)K2(x)\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.png
                LogLog plot for the evoluion of Y in freeze-in mechanism.

                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*1019

                sv=1.75*10-9

                sol1=NDSolve[{y'[x]==-x-2*(y[x]2-(0.192*0.2*mpl*m*sv*x1.5*e-x)2),

                y[1]==0.192*0.2*mpl*m*sv*e-1},y,{x,1,20}]

                bc=Evaluate[y[20]/.sol1]

                sol2=NDSolve[{y'[x]==-x-2*(y[x]2-(0.192*0.2*mpl*m*sv*x1.5*e-x)2),

                y[20]==bc},y,{x,20,10000}]

                LogLogPlot[Evaluate[y[x]/.sol1],{x,1,20},PlotRange->{{1,10000},{1,1012}}]]

                LogLogPlot[Evaluate[y[x]/.sol2],{x,20,10000},PlotRange->{{1,10000},{1,1012}}]

                LogLogPlot[0.192*0.2*m*mpl*sv*x1.5*e-x,{x,1,1000},PlotRange->{{1,10000},{1,1012}},PlotStyle->RGBColor[0.93,0.4,0.33]]

                Show[%,%%,%%%]

                2. Asymmetric Dark Matter

                m=100m=100

                mpl=2.44×1018m_{pl}=2.44\times10^{18}

                sv=1.5×1010sv=1.5\times10^{-10}

                k=1011k=10^{-11}

                g=2g=2

                gs=90g_s=90

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

                sol1=NDSolve[{y[x]==svx2(y[x]2(k/lam)y[x](0.1454πmmpl(g/gs))2x3e2x),y[1]==0.145πmmpl(g/gs)11.5e1},y,{x,1,20}]sol1 = NDSolve[\{y&#x27;[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\}]

                bc=Evaluate[y[20]/.sol1]\displaystyle bc=Evaluate[y[20]/.sol1]

                sol3=NDSolve[{y[x]==svx2(y[x]2+(k/lam)y[x](0.1454πmmpl(g/gs))2x3e2x),y[1]==0.145πmmpl(g/gs)11.5e1.1(lamk)},y,{x,1,20}]sol3 = NDSolve[\{y&#x27;[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]cd=Evaluate[y[20]/.sol3]

                sol4=NDSolve[{y[x]==svx2(y[x]2+(k/lam)y[x](0.1454πmmpl(g/gs))2x3e2x),y[20]==cd},y,{x,20,10000}]sol4 = NDSolve[\{y&#x27;[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&gt;{{1,10000},{1,1012}]LogLogPlot[Evaluate[y[x]/.sol1]/lam,\{x,1,20\},PlotRange-&gt;\{\{1,10000\},\{1,10^{-12}\}]

                LogLogPlot[Evaluate[y[x]/.sol2]/lam,{x,20,10000},PlotRange&gt;{{1,10000},{1,1012}]LogLogPlot[Evaluate[y[x]/.sol2]/lam,\{x,20,10000\},PlotRange-&gt;\{\{1,10000\},\{1,10^{-12}\}]

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