Theoretical foundations of the reduced relativistic gas in the cosmological perturbed context

TL;DR

This paper develops a theoretical framework for the Reduced Relativistic Gas (RRG) in cosmological perturbation theory, validating its use as a perfect fluid approximation for warm dark matter and other cosmological components.

Contribution

It derives the Einstein-Boltzmann equations for RRG and analyzes the validity of the perfect fluid approximation in cosmological scales.

Findings

Perfect fluid approximation valid for keV-mass particles at large scales

Derived initial conditions for RRG in early universe

Studied potential evolution in a toy RRG universe

Abstract

The Reduced Relativistic Gas (RRG) is a simplified version of the ideal relativistic gas, which assumes that all particles have the same momentum magnitude. Although this is a very idealized situation, the resulting model preserves the phenomenology of Maxwell-Boltzmann distribution and, in some situations, can be described as a perfect fluid, without introducing large errors in both cosmological background and first-order perturbations. The perfect fluid description of RRG model was already used to study the warmness of dark matter, massive neutrinos and interaction of baryons and photons before recombination, showing very good agreement with previous works based on the full Einstein-Boltzmann system of equations. In order to understand these results and construct a more general and formal framework for RRG, we develop a theoretical description of first-order cosmological perturbations of RRG, based on a distribution function which encodes the simplifying assumption that all particles have the same momentum magnitude. The full set of Einstein-Boltzmann equations for RRG distribution are derived and quantities beyond the perfect fluid approximation are studied. Using RRG to describe warm dark matter, we show that, for particles with $m \sim \text{keV}$, the perfect fluid approximation is valid on scales $k < 10\, \text{h}/\text{Mpc}$, for most of the universe evolution. We also determine initial conditions for RRG in the early universe and study the evolution of potential in a toy model of universe composed only by RRG.

Tables

Table 1. Table 1 : Six values for the parameter b 2 superscript 𝑏 2 b^{2} in a scenario of WDM and their corresponding values of mass and thermal velocity.

$\#$	Ref. [11]	This work
$b^2$	$m_{\omega }\text{ [keV]}$	$m\text{ [keV]}$	$v_{thℛ,0}\text{ [m/s]}$
${10}^{-16}$	11.68	40.7	3.00
${10}^{-15}$	4.65	12.87	9.48
${10}^{-14}$	1.85	4.07	30.00
${10}^{-13}$	0.73	1.28	94.86
${10}^{-12}$	0.29	0.40	300.00
${10}^{-11}$	0.11	0.12	948.68

Equations

$$f_{\mathcal{R}}\left(p,p_{\mathcal{R}}\right)\equiv Cp\delta \left(p-p_{\mathcal{R}}\right),$$

$$n_{\mathcal{R}}=\int \frac{d^{3}p}{{\left(2\pi \right)}^3}{f}_{\mathcal{R}}\Rightarrow C=\frac{2{\pi }^2{n}_{\mathcal{R}}}{p_{\mathcal{R}}^3}\text{.}$$

$${\rho }_{\mathcal{R}}=\int \frac{d^{3}p}{{\left(2\pi \right)}^3}E\left(p\right)f_{\mathcal{R}}=n_{\mathcal{R}}E\left(p_{\mathcal{R}}\right)$$

$$P_{\mathcal{R}}=\int \frac{d^{3}p}{{\left(2\pi \right)}^3}\frac{1}{3}\frac{p^2}{E\left(p\right)}{f}_{\mathcal{R}}=\frac{n_{\mathcal{R}}}{3}\frac{p_{\mathcal{R}}^2}{E\left(p_{\mathcal{R}}\right)}.$$

$$P_{\mathcal{R}}=\frac{{\rho }_{\mathcal{R}}}{3}\left[1-{\left(\frac{{\rho }_{d\mathcal{R}}}{{\rho }_{\mathcal{R}}}\right)}^2\right]\text{,}$$

$${\rho }_{\mathcal{R}}={\rho }_{d\mathcal{R}}\left[\frac{3}{2}\left(\frac{P_{\mathcal{R}}}{{\rho }_{d\mathcal{R}}}\right)+\sqrt{\frac{9}{4}{\left(\frac{P_{\mathcal{R}}}{{\rho }_{d\mathcal{R}}}\right)}^2+1}\right]\text{.}$$

$$P_{\text{MB}}=nT\text{ and }{\rho }_{\text{MB}}=3nT+nm\frac{K_1\left(m/T\right)}{K_2\left(m/T\right)}\text{,}$$

$${\rho }_{\text{MB}}={\rho }_{d\text{MB}}\left[3\left(\frac{P_{\text{MB}}}{{\rho }_{d\text{MB}}}\right)+\frac{K_1\left({\rho }_{d\text{MB}}/P_{\text{MB}}\right)}{K_2\left({\rho }_{d\text{MB}}/P_{\text{MB}}\right)}\right]\text{.}$$

$$v_{\text{thMB}}=\sqrt{1-{\left[\frac{K_1\left({\rho }_{d\text{MB}}/P_{\text{MB}}\right)}{K_2\left({\rho }_{d\text{MB}}/P_{\text{MB}}\right)}+3\frac{P_{\text{MB}}}{{\rho }_{d\text{MB}}}\right]}^{-2}}\text{.}$$

$$v_{\text{th}\mathcal{R}}=\sqrt{1-{\left(\frac{{\rho }_{d\mathcal{R}}}{{\rho }_{\mathcal{R}}}\right)}^2}=\sqrt{1-4{\left(3\frac{P_{\mathcal{R}}}{{\rho }_{d\mathcal{R}}}+\sqrt{9{\left(\frac{P_{\mathcal{R}}}{{\rho }_{d\mathcal{R}}}\right)}^2+4}\right)}^{-2}}.$$

$$F_{\rho }\equiv \frac{{\rho }_{\mathcal{R}}-{\rho }_{\text{MB}}}{{\rho }_{\text{MB}}}\text{ and }{F}_{v_{\text{th}}}\equiv \frac{v_{\text{th}\mathcal{R}}-v_{\text{thMB}}}{v_{\text{thMB}}}\text{.}$$

$$\frac{1}{3}\frac{dV}{V}=-\frac{U_{\mathcal{R}}d{U}_{\mathcal{R}}}{U_{\mathcal{R}}^2-U_{d\mathcal{R}}^2},$$

$${\rho }_{\mathcal{R}}\left(a\right)={\rho }_{d\mathcal{R},\text{ref}}{\left(\frac{a_{\text{ref}}}{a}\right)}^3\sqrt{1+b^2{\left(\frac{a_{\text{ref}}}{a}\right)}^2}$$

$${\omega }_{\mathcal{R}}\left(a\right)=\frac{1}{3}{v}_{\text{th}\mathcal{R}}^2=\frac{1}{3}\left[\frac{b^2}{{\left(a/a_{\text{ref}}\right)}^2+b^2}\right],$$

$$b=\frac{v_{\text{th}\mathcal{R},\text{ref}}}{\sqrt{1-v_{\text{th}\mathcal{R},\text{ref}}^2}}\text{.}$$

$$H^2=\frac{8\pi G}{3}\left({\rho }_{\mathcal{R}}+{\rho }_{\text{NR}}+{\rho }_{\text{UR}}+{\rho }_{\Lambda }\right)-\frac{\kappa }{a^2}\text{,}$$

$${\dot{u}}^2=4\frac{8\pi G}{3}{u}^2{\rho }_{\mathcal{R}}\left(u\right)=4H_0^2{\Omega }_{d\mathcal{R},\text{ref}}\sqrt{u+b^2}\text{,}$$

$$a\left(t\right)=a_{\text{ref}}\sqrt{{\left(\frac{3H_0\sqrt{{\Omega }_{d\mathcal{R},\text{ref}}}}{2}t+b^{\frac{3}{2}}\right)}^{\frac{4}{3}}-b^2\text{.}}$$

$$b^2={\left[\frac{K_1\left(x_0\right)}{K_2\left(x_0\right)}+\frac{3}{x_0}\right]}^2-1\text{,}$$

$$T^{{}_{\text{MB}}}\propto\left\{\begin{array}[]{lll}D/a&\text{for}&T>m\\ Q/a^{2}&\text{for}&T\leq m\end{array}\right.\text{,}$$

$$T^{{}_{\text{MB}}}\left(a\right)=\left\{\begin{array}[]{lll}T_{\gamma 0}\left(\frac{a_{0}}{a}\right)&\text{for}&T>m\\ T_{\gamma 0}\left(\frac{T_{\gamma 0}}{m}\right)\left(\frac{a_{0}}{a}\right)^{2}&\text{for}&T\leq m\end{array}\right.\text{.}$$

$$x_0=\frac{m}{T_0^{{}_{\text{MB}}}}={\left(\frac{m}{T_{\gamma 0}}\right)}^2\text{.}$$

$$b^2\simeq 3{\left(\frac{T_{\gamma 0}}{m}\right)}^2\text{ or }m\simeq \sqrt{3}\left(\frac{T_{\gamma 0}}{b}\right)\simeq \frac{4.07}{10^{7}b}\text{ keV.}$$

$$m_{\omega }=\frac{4.65}{10^6{b}^{4/5}}\text{ keV.}$$

$$d{s}^2=-\left(1+2\Psi \right)d{t}^2+a^2\left(1+2\Phi \right){\delta }_{ij}d{x}^{i}d{x}^j\text{,}$$

$$P^{\mu }=\left[E\left(1-\Psi \right),p\frac{\left(1-\Phi \right)}{a}{\hat{p}}^i\right]\text{.}$$

$$p_{\mathcal{R}}\left(\vec{x},\hat{p},t\right)={\bar{p}}_{\mathcal{R}}\left(t\right)+\delta p_{\mathcal{R}}\left(\vec{x},\hat{p},t\right),$$

$$f_{\mathcal{R}}\left(\vec{x},p,\hat{p},t\right)=Cp\delta (p-p_{\mathcal{R}})\text{,}$$

$$f_{\mathcal{R}}\simeq {\bar{f}}_{\mathcal{R}}+{\frac{\partial f_{\mathcal{R}}}{\partial p_{\mathcal{R}}}|}_{p_{\mathcal{R}}={\bar{p}}_{\mathcal{R}}}\delta p_{\mathcal{R}}={\bar{f}}_{\mathcal{R}}+{\mathcal{F}}_{\mathcal{R}}\text{,}$$

$${\mathcal{F}}_{\mathcal{R}}\left(\vec{x},p,\hat{p},t\right)\equiv -p\frac{\partial {\bar{f}}_{\mathcal{R}}}{\partial p}\mathcal{M}\left(\vec{x},\hat{p},t\right)\equiv -\frac{p}{{\bar{p}}_{\mathcal{R}}}\frac{\partial {\bar{f}}_{\mathcal{R}}}{\partial p}\delta p_{\mathcal{R}}$$

Full text

1 Introduction

The ideal relativistic gas is a good description of a real gas when interactions between the constituent particles can be neglected. In the classical regime, when average occupation numbers are much smaller than unity, this model is described by the Maxwell-Boltzmann distribution function 222Frequently, the Maxwell-Boltzmann distribution function is called Jöuttner-Synge distribution function [1] due to F. Jüttner’s work in 1911 [2] and J.L. Synge’s in 1957 [3].. Despite the simplicity of the model, the resulting equations in cosmological context can not be solved analytically, both in background and first-order perturbations. With a further simplifying condition, that all particles of the relativistic gas have the same momentum magnitude, the relativistic gas and the cosmological equations become simpler, at the cost of an error smaller than $2.5\%$ in the density evolution with respect to Maxwell-Boltzmann description [4]. This model is called Reduced Relativistic Gas (RRG), and was first introduced in 1966 by Sakharov for Cosmic Microwave Background applications [5].

The RRG was recently revived in Ref. [4], where background analytical solutions for scenarios with radiation and cosmological constant were studied. It was also used to provide a simpler fluid description of Warm Dark Matter (WDM), which, based on linear power spectrum data, produced stringent constraints on the thermal velocities of the particles without the necessity of solving the full set Einstein-Boltzmann equations [6]. Moreover, analytical solutions for flat and non-flat background in the presence of other two or three components (relativistic and non-relativistic matter and cosmological constant) were found in Ref. [7]. Other uses of RRG include description of barion-photon interaction before recombination [8], massive neutrinos [9], studies of relativistic matter in anisotropic universes [10] and a model independent description of warm dark matter [11], where it was found that the transfer function for RRG model reproduces the full Einstein-Boltzmann calculation for thermal relics [12] with $1\%$ precision.

Given these interesting applications of RRG model in cosmology, it is valuable to study the model in a more general theoretical framework. By studying the RRG beyond the perfect fluid description we can have a better understanding of its approximations and more effects such as dissipation and collisional terms can be introduced. In this paper, we implement the main idea of RRG, i.e., that all particles have same momentum magnitude, in distribution function, both for background and first-order perturbations, and derive the corresponding Einstein-Boltzmann equations. Moreover, we also identify the conditions which reproduce the fluid description of RRG.

The paper is organized as follows. In Sec. 2, we present the concepts of RRG, review its background evolution and estimate the mass of WDM particles as a function of RRG warmness parameter. In Sec. 3, we derive the first-oder Einstein-Boltzmann system of equations for RRG. In Sec. 4, we study the warm approximation of RRG, which can be used to study WDM in model independent manner, find the initial conditions for RRG perturbations in the early universe and study a toy model of a universe composed solely by RRG. We make our final comments in Sec. 5.

2 Background: concepts and basic equations

The RRG is based on the premise which all particles have the same magnitude for the momentum. Thus, one can reasonably assume its distribution function is proportional to a Dirac delta function in the momentums i.e.

[TABLE]

where $p_{\mathcal{R}}$ sets the magnitude for the momentums and $C$ is a dimensionless constant given by

[TABLE]

where $n_{\mathcal{R}}$ is the numerical particle density of RRG. Besides, the energy density $\rho_{\mathcal{R}}$ and pressure $P_{\mathcal{R}}$ for RRG are given by

[TABLE]

and

[TABLE]

Combining these equations and using the dispersion relation $E^{2}-p_{\mathcal{R}}^{2}=m^{2}$, it is possible to obtain the equation of state of RRG

[TABLE]

where $\rho_{d\mathcal{R}}\equiv n_{\mathcal{R}}m$ is the rest energy density. In the non-relativistic (NR) and ultra-relativistic (UR) limits, we have $\rho_{\mathcal{R}}\simeq\rho_{d\mathcal{R}}\Rightarrow$ $P_{\mathcal{R}}\simeq 0$ and $\rho_{\mathcal{R}}\gg\rho_{d\mathcal{R}}$ that $P_{\mathcal{R}}\simeq\rho_{\mathcal{R}}/3$, respectively. It is also possible to isolate $\rho_{\mathcal{R}}$ in the last equation by obtaining

[TABLE]

The first to deduce the equation of state of RRG from an assumption about its distribution function were S. C. dos Reis and I. L. Shapiro [10].

As we know, the standard description for a classical relativistic ideal gas of massive particles is given by Maxwell-Boltzmann statistics (MB), whose pressure $P_{\text{MB}}$ and energy density $\rho_{\text{MB}}$ are expressed as [2, 1]

[TABLE]

where $n\equiv N/V$ is the numerical particle density, $T$ is the temperature and $K_{1}$ and $K_{2}$ are modified Bessel functions of the second kind. Combining these two equations, we can write

[TABLE]

From the thermal velocity definition $v_{\text{th}}\equiv\sqrt{(\left\langle E\right\rangle^{2}-m^{2})/\left\langle E\right\rangle^{2}}$ [1], we obtain for Maxwell-Boltzmann gas the following expression:

[TABLE]

On the other hand, due to the peculiar RRG distribution function, its thermal velocity is given by:

[TABLE]

Considering that both gases have the same mass and numerical particle density, so that $\rho_{d\mathcal{R}}=\rho_{d\text{MB}}=\rho_{d}$, we can compare them from the functions:

[TABLE]

In Fig. 1 we show the behavior of these two functions in relation to $P/\rho_{d}$. From this figure, we found that the maximum relative deviation is about $2.5\%$ for the energy density (as was already obtained in Ref. [4]) and about $1.6\%$ for the thermal velocity. We also verified that the relative difference between Maxwell-Boltzmann statistic and RRG model reaches its maximum in a relativistic region where the fraction $P/\rho_{d}$ is close to $1/2$. In the NR ($P/\rho_{d}\simeq 0$) and UR ($P/\rho_{d}\gg 1$) regions, this difference becomes completely negligible.

2.1 RRG in the cosmological context

Now let us consider the evolution of RRG density in Friedmann-Lemaître-Robertson-Walker metric, where $V\propto a^{3}$. Replacing Eq. (2.5) in conservation law $dU=-PdV$, we obtain

[TABLE]

where $U_{\mathcal{R}}=$ $V\rho_{\mathcal{R}}$ is the total energy and $U_{d\mathcal{R}}=V\rho_{d\mathcal{R}}=\text{const.}$ represents the rest energy density. Integrating Eq. (2.12), the evolution of $\rho_{\mathcal{R}}$ with scale factor can be written in the form [4, 11]

[TABLE]

where $b=\sqrt{U^{2}_{\mathcal{R},\text{ref}}/U^{2}_{d\mathcal{R}}-1}$ is a dimensionless parameter which quantifies the relativistic degree of RRG and is defined at some chosen $a_{\text{ref}}$.

Moreover, the function $\omega\equiv P/\rho$ for RRG is given by

[TABLE]

where for $a=a_{\text{ref}}$,

[TABLE]

Observe that $\omega_{\mathcal{R}}\simeq v_{\text{th}\mathcal{R}}\simeq 0$ for $a\gg ba_{\text{ref}}$ (NR limit) and for $a\ll ba_{\text{ref}}$ (UR limit) we have $v_{\text{th}\mathcal{R}}^{2}\simeq 1$ and $\omega_{\mathcal{R}}\simeq 1/3$. Therefore, $b$ is a relative parameter, since this measure is established by comparing the value of the scale factor $a$ with the product of $ba_{\text{ref}}$. Thus, the value of the parameter $b$ that quantify the relativistic degree of the fluid depends on the specific choice made for $a_{\text{ref}}$.

For a better interpretation of these relations, see Fig. 2. In the figure to the left, we show the dependence between the energy density and the scale factor for NR (continuous black curve), UR (continuous red curve) and RRG (dashed blue curve) components, assuming $b=1$. Note that for $a\ll ba_{\text{ref}}$ and $a\gg ba_{\text{ref}}$ RRG exhibits a typically behavior of UR and NR components, respectively. The $\rho_{\mathcal{R}}\left(a\right)$ interpolates between a typical dependence of a NR component to a UR one at $a\sim ba_{\text{ref}}$. The same conclusion is obtained analyzing the dependence on the scale factor of the function $\omega_{\mathcal{R}}$ in Fig. 2 to the right. Because of a smooth interpolation between these two regimes, RRG can be used to describe the period when the matter content of the universe is in the transition from the radiation to NR matter [10].

Friedmann equation for a universe with curvature $\kappa$, cosmological constant $\Lambda$, RRG, NR and UR components can be expressed as follows:

[TABLE]

where $H\equiv\dot{a}/a$ is the Hubble parameter and "$\cdot$" represent the derivative with respect to time. The simplicity of the equation of state for RRG makes it possible to obtain analytical solutions of Friedmann equation with one or multiple components including RRG (see Refs. [4, 7, 9] for details).

For instance, considering only RRG and introducing a new variable $u=a^{2}/a_{\text{ref}}^{2}$ in Eq. (2.16) we get

[TABLE]

where $\Omega_{d\mathcal{R},\text{ref}}\equiv\rho_{d\mathcal{R},\text{ref}}/\rho_{\text{crit},0}$ with $\rho_{\text{crit},0}=3H_{0}^{2}/8\pi G$. The subscript "[math]" indicates the current value of the quantities. Therefore, integrating Eq. (2.17) and assuming that $a\left(t\simeq 0\right)\simeq 0$, we obtain

[TABLE]

Observe that for $a\gg ba_{\text{ref}}$ (NR limit) the standard NR behavior $a\left(t\right)\propto t^{2/3}$ is obtained. On the other hand, assuming $a\ll ba_{\text{ref}}$ (UR limit) leads to typical UR behavior $a\left(t\right)\propto\sqrt{t}$.

Therefore, the dependence of the scale factor with the time given by Eq. (2.18), interpolates between the typical solutions of a universe dominated by UR matter to NR ones at $a/a_{\text{ref}}\sim b$. This statement is illustrated in Fig. 3, where we plot the Eq. (2.18) for different values of $b$.

2.2 The relationship between the parameter b and the mass m

In the description via RRG, the relativistic degree is defined by the ratio $a/(a_{\text{ref}}b)$. On the other hand, in the description via Maxwell-Boltzmann statistics the relativistic degree is determined by the ratio of the mass of the particle to its temperature $m/T$. To obtain a consistent relationship between these two descriptions it is reasonable to assume that at some specific time $t_{\text{ref}}$ both are equal and then independently evolve. Thus, at $t_{\text{ref}}$, the two descriptions have the same energy density and numerical density of particles, i.e., $\rho_{\mathcal{R}}\left(a_{\text{ref}}\right)=\rho_{\text{MB}}\left(m/T_{\text{ref}}\right)$ and $n_{\mathcal{R}}\left(a_{\text{ref}}\right)=n_{\text{MB}}\left(m/T_{\text{ref}}\right)$. Moreover, since $a_{\text{ref}}$ is an arbitrary scale factor, we adopt for convenience $a_{\text{ref}}=a_{0}$ which represents the current scale factor.

In this case, combining the Eqs. (2.13) and (2.8) at $t=t_{\text{ref}}=t_{0}$, we can obtain the following relation:

[TABLE]

where $x_{0}\equiv\rho_{d\text{MB}}/P_{\text{MB}}=m/T_{0}^{\text{MB}}$ and $T_{0}^{\text{MB}}$ is the current temperature of the relativistic Maxwell-Boltzmann gas. Therefore, it is necessary to determine the $T_{0}^{\text{MB}}$ to obtain a unique relationship between $b$ and $m$. The difficulty here is the $T_{0}^{\text{MB}}$ generally depends on the cosmological model adopted. Nevertheless, from reasonable assumptions, we can get a relatively general estimate for it.

To obtain the $T_{0}^{\text{MB}}$ in a simple and approximate way we assume that the gas of relativistic particles undergoes an abrupt transition from UR to NR as the universe expands. In the UR case, one has known that $T_{\text{UR}}^{\text{MB}}\propto a^{-1}$. On the other hand, for NR case, it can be shown from the conservation equation that $T_{\text{NR}}^{\text{MB}}\propto a^{-2}$.

The next step is to assume that the abrupt transition from UR to NR occurs at $T_{\text{tran}}^{\text{MB}}=m$ such that

[TABLE]

where $D$ and $Q$ are constants to be determined. In addition, we suppose that at some point in the early universe Maxwell-Boltzmann gas was in equilibrium with the photons. Thus, neglecting any extra contribution of photons due to the pair annihilation processes333The most common example of such processes occurs shortly before nucleosynthesis when the annihilation of electron-positron pairs causes a slight increase in photon temperature with respect to neutrino temperature., the adiabatic expansion ensures $T_{\gamma}=T_{\text{UR}}^{{}_{\text{MB}}}$. Thus, since $T_{\gamma}=T_{\gamma 0}\left(a_{0}/a\right)$ and $T_{\gamma}=T_{\text{UR}}^{{}_{\text{MB}}}$, we have $D=T_{\gamma 0}a_{0}$. Setting $a\left(T_{\text{tran}}^{\text{MB}}\right)=a_{\text{tran}}$, we have $T_{\text{UR}}^{{}_{\text{MB}}}\left(a_{\text{tran}}\right)=T_{\text{NR}}^{{}_{\text{MB}}}\left(a_{\text{tran}}\right)$, so that $Q=T_{\gamma 0}a_{0}a_{\text{tran}}$ where $a_{\text{tran}}/a_{0}=T_{\gamma 0}/T_{\text{tran}}^{\text{MB}}=T_{\gamma 0}/m$. Therefore,

[TABLE]

According to the Eq. (2.21), the temperature $T^{{}_{\text{MB}}}(a)$ calculated today is $T_{0}^{{}_{\text{MB}}}=(T_{\gamma 0})^{2}/m$, then the parameter $x_{0}$ is given by:

[TABLE]

Substituting this last expression into Eq. (2.19), we can obtain an expression for $b$ as a function of mass $m$, where $T_{\gamma 0}=2.3486\times 10^{-7}$keV [13].

In the specific case of WDM, it is natural to expect $m\gg T_{\gamma 0}$. Thus, adopting the RRG to describe the WDM, we can use the asymptotic expansions for the Bessel functions to approximate the Eq. (2.19) obtaining the expression $b^{2}\simeq 3/x_{0}$. Finally, using Eq. (2.22) in the latter result, we conclude that

[TABLE]

Note that due to the simplifications used this result should be regarded as an approximate relation between $m$ and $b$.

A more complex procedure involving the transfer function of thermal relics [12] can be used to obtain an equation analogous to Eq. (2.23). From a numerical method using RRG and thermal relics, W. S. Hipólito-Ricaldi et al. [11] obtained the following fit-formula:

[TABLE]

In Fig. 4 we plot the mass as a function of the parameter $b^{2}$ given by Eq. (2.23) and Eq. (2.24). Analyzing this figure, we verify that the expression obtained in our work and that obtained in Ref. [11] are equivalent in order of magnitude. In addition, we observed the largest discrepancy occurs for smaller values of $b$ when we consider a range between $10^{-16}$ and $10^{-6}$. These same conclusions can be obtained from Table 1 where we assign six values for $b^{2}$ and show the equivalent value of mass and thermal velocity. It is worth noting that constraints from high redshift Lyman-$\alpha$ forest data impose a lower bound $m_{\text{WDM}}\geq 3.3$ keV for particles of WDM in the form of early decoupled thermal relics [14]. Thus, using Eq. (2.23) and Eq. (2.24) we find that the values of $b^{2}$ corresponding to this mass limit are $b^{2}\lesssim 1.52\times 10^{-14}$ and $b^{2}\lesssim 2.35\times 10^{-15}$, respectively. As shown in Ref. [11], the RRG reproduces the relative standard WDM transfer function with an accuracy of $1\%$ for $b^{2}\lesssim 10^{-14}$.

3 Linear perturbation level: distribution function and dynamic equations

In the linear level, we restrict the analyses to scalar perturbations in a spatially flat universe. The scalar metric in the Newtonian gauge is given by

[TABLE]

where the perturbation of the metric $\Psi(x,t)$ corresponds to the Newtonian potential and $\Phi(x,t)$ is associated with the perturbation of the spatial curvature. We define the quadri-momentum as $P^{\nu}\equiv dx^{\nu}/d\lambda$ and the energy by $E\equiv\sqrt{p^{2}+m^{2}}$, where $m$ and $p^{2}\equiv g_{ij}P^{i}P^{j}$ are the mass and proper momentum of the particle, respectively. Using $g_{\mu\nu}P^{\mu}P^{\nu}=-m^{2}$, it is possible to show that [15]:

[TABLE]

3.1 Distribution function and energy-momentum tensor of RRG

In the linear perturbation level, we relax the basic assumption of RRG by allowing a small variation around the momentum $\bar{p}_{\mathcal{R}}\left(t\right)$ of the particles, so the new momentum will be

[TABLE]

with $\bar{p}_{\mathcal{R}}\left(t\right)\gg\delta p_{\mathcal{R}}\left(\vec{x},\widehat{p},t\right)$. From here on, our bar will be used to identify background quantities. Observe the momentum $p_{\mathcal{R}}$ is a parameter of distribution function (such as the temperature for photons) which in the background is only time dependent. On the other hand, at a linear perturbation level, we consider the most general dependence for $p_{\mathcal{R}}$, namely, a dependence on $\vec{x}$ (inhomogeneities), $\widehat{p}$ (anisotropies) and obviously on time $t$. By definition of RRG all the particles have the same magnitude of momentum hence we do not consider $p_{\mathcal{R}}$ depending on the magnitude of the moment $p$.

Thus, the RRG distribution function can be rewritten as

[TABLE]

where $p_{\mathcal{R}}$ can be interpreted as a constraint parameter of the momentum variable $p$. In linear perturbation level, Eq. (3.4) can be expanded as

[TABLE]

where $\bar{f}_{\mathcal{R}}\left(p,t\right)\equiv Cp\delta\left(p-\bar{p}_{\mathcal{R}}\right)$ is the background RRG distribution function and

[TABLE]

is the first-order RRG distribution function. Note that, $p_{\mathcal{R}}$ is for RRG distribution function which the temperature $T$ is for the photon distribution function. In this sense, $\bar{p}_{\mathcal{R}}$, $\delta p_{\mathcal{R}}$ and $\mathcal{M}$ play the role of $\bar{T}$, $\Delta T$ and $\Theta\equiv\Delta T/\bar{T}$, respectively.

Once the RRG distribution function is established, we can find its energy-momentum tensor from the general expression [15, 16],

[TABLE]

To proceed with the calculations, we performed a multipolar expansion of $\mathcal{M}\left(\vec{x},\widehat{p},t\right)$ in Fourier space,

[TABLE]

where $P_{l}\left(\mu\right)$ are the Legendre polynomials, $\mu\equiv\widehat{k}^{i}\widehat{p}_{i}=\widehat{k}\cdot\widehat{p}$ and $\mathcal{M}_{l}$ (with $l=0,1,2,\cdots$) are the multipoles of the momentum contrast, given by the expression

[TABLE]

Using the metric (3.1), the momentum $P^{\mu}$ and the distribution function, we can show that

[TABLE]

where

[TABLE]

and

[TABLE]

are the energy density and pressure perturbations, respectively. The irrotational velocity perturbation $v_{\mathcal{R}}$ in the Fourier space is given by the expression:

[TABLE]

The $\pi_{\mathcal{R}}{}^{i}{}_{j}$ is the traceless longitudinal part of the energy-momentum tensor, defined as

[TABLE]

and quantified by the variable $\sigma_{\mathcal{R}}$

[TABLE]

where $\widehat{k}^{i}\equiv k^{i}/\sqrt{k_{j}k^{j}}$ is the unit wavevector which arises when we take the Fourier transform.

Combining Eq. (3.11) with Eq. (3.12), the pressure perturbation is given by

[TABLE]

This equation can also be obtained perturbing Eq. (2.5) with respect to the momentum by taking into account the role of $\rho_{d\mathcal{R}}$ without referring to any particular gauge. Thus, it can be identified as the sound speed squared in the rest frame of RRG:

[TABLE]

Moreover, one can check $c^{2}_{s}$ is equal to the adiabatic sound speed squared, $c^{2}_{a}=\dot{P_{\mathcal{R}}}/\dot{\rho_{\mathcal{R}}}$, hence RRG induces no intrinsic entropy perturbations.

Applying UR limit $\bar{\rho}_{d\mathcal{R}}^{2}/\bar{\rho}_{\mathcal{R}}^{2}\ll 1$ in the previous equations, we obtain

[TABLE]

and

[TABLE]

Note that these expressions have the same functional form as those obtained for photons and massless neutrinos. On the other hand, applying NR limit $\bar{\rho}_{d\mathcal{R}}^{2}/\bar{\rho}_{\mathcal{R}}^{2}\rightarrow 1$, we obtain similar expressions for those of cold dark matter and baryons:

[TABLE]

and

[TABLE]

At NR limit, we can not neglect $\sqrt{3\omega_{\mathcal{R}}}$, because the density perturbations themselves induce velocities through the continuity equation [15]. Therefore, in the linear treatment of RRG at NR limit, we maintain the $\sqrt{3\omega_{\mathcal{R}}}$ terms and neglect the $\omega_{\mathcal{R}}$ ones, which are equivalent to terms of first and second order in $\bar{p}_{\mathcal{R}}/\bar{E}$, respectively.

3.2 Einstein and conservation equations

Using the metric (3.1) and assuming a universe constituted solely (or dominated) by RRG whose energy-momentum tensor is expressed by Eq. (3.10), the perturbed Einstein equation in Fourier space leads to the following expressions:

[TABLE]

and

[TABLE]

In this same context, the conservation equation of the energy-momentum tensor leads to the well-known energy and momentum equation for fluids with the correspond RRG quantities:

[TABLE]

and

[TABLE]

where $\delta_{\mathcal{R}}=\delta\rho_{\mathcal{R}}/\bar{\rho}_{\mathcal{R}}$ is the density contrast. As expected, in the absence of interactions, we obtain equations with functional forms identical to those of cold dark matter (or baryons) and massless neutrinos (or photons), when we apply NR and UR limits, respectively.

It is important to emphasize that up to now, our description of RRG in the perturbed cosmological context does not constitute a closed system of the equation since there are six variables ($\Phi,\Psi,\delta\rho_{\mathcal{R}},v_{\mathcal{R}},\sigma_{\mathcal{R}},\delta P_{\mathcal{R}}$) and only five independent equations namely Eqs. (3.22), (3.23), (3.24), (3.25) and (3.16). If for some physical reason the anisotropic stress term $\sigma_{\mathcal{R}}$ can be neglected or determined as a function of other variables (e.g. Ref. [17]) this set of equations becomes a closed set. Otherwise, it is necessary to use the Boltzmann equation to complete the description. In the following section, we study RRG in this context.

3.3 Boltzmann equation

The Boltzmann equation in the differential form is given by

[TABLE]

The right-hand side represents the collision terms $C[f]$ which are usually complicated functions of the distribution function. In the Newtonian gauge, we can express the complete derivative $df/dt$ in terms of partial derivatives, so

[TABLE]

For the zero order of perturbation, the previous equation becomes:

[TABLE]

where we establish $C[f]=0$ since we are interested only in the collisionless Boltzmann equation. In the specific case of RRG, we can manipulate the previous expression such that

[TABLE]

Hence, integrating both sides with respect $p$, we conclude that

[TABLE]

This if implies the RRG distribution function satisfies the zero-order Boltzmann equation $\bar{p}_{\mathcal{R}}\propto a^{-1}$. This is expected since we can show from the geodesic equation that the momentum $p$ of any particle is proportional to $a^{-1}$ in a Friedmann universe.

For the first-order of perturbation, we can express the left-hand side of Eq. (3.29) as

[TABLE]

Using Eq. (3.6), the relation $d\bar{p}_{\mathcal{R}}/dt=-\bar{p}_{\mathcal{R}}H$ and properties of the Dirac delta function we can rewrite the time derivative in terms of $p$. In the absence of collision terms, we obtain the following expression for the Boltzmann equation in the Fourier space:

[TABLE]

This expression is directly related to the equation of energy-momentum tensor conservation since by multiplying Eq. (3.34) by $Ed^{3}p/(2\pi)^{3}$ and integrating it, we obtain Eq. (3.26). Similarly, multiplying Eq. (3.34) by $p\mu d^{3}p/(2\pi)^{3}$ and integrating it, we obtain Eq. (3.27). It is worth noting that the first term between parentheses in Eq. (3.34) can not be neglected, since $\bar{f}_{\mathcal{R}}=Cp\delta\left(p-\bar{p}_{\mathcal{R}}\right)$ and therefore it must be considered at the time of integration.

To get new equations necessary to have a closed system of equations, we must go to higher orders of $\mu$. We know the highest multipolar moments are always subdominant compared to their predecessors, and they decay rapidly when the particles become NR. Thus, it is possible to choose a much smaller $l_{\text{max}}$ for the massive particles than for the approximately massless ones [18]. If the interest is to apply RRG to describe slightly relativistic particles, where the quadrupole and higher terms can be neglected, there is no need to use Boltzmann hierarchical equations. However, if we wish to apply RRG to describe relativistic or UR particles, where the term of quadrupole and higher are relevant, it is necessary to use Boltzmann hierarchical equations to complete the system of equations.

The Boltzmann hierarchical equations are obtained multiplying Eq. (3.34) by $p\mu^{m}d^{3}p/(2\pi)^{3}$ where $m=0,1,2,\cdots$ and integrating it. Proceeding in this way, it is possible to obtain an infinite set of hierarchical equations:

[TABLE]

This set of equations governs the evolution of RRG distribution in the phase space. Note that in UR limit, where $\bar{\rho}_{d\mathcal{R}}^{2}/\bar{\rho}_{\mathcal{R}}^{2}\ll 1$, the Boltzmann hierarchical equations for RRG have the same functional form as the massless neutrinos and the collisionless photons (see, e.g., Ref. [16]).

The practical use of this infinite set of hierarchical equations occurs through truncation in some maximal multipolar order $l_{\text{max}}$. A simple but imprecise method is to define $\mathcal{M}_{l}=0$ for $l>l_{\text{max}}$. A better truncation scheme is based on the extrapolation of the behavior of $\mathcal{M}_{l}$ to $l=l_{\text{max}}+1$, presented in Ref. [18].

The deduction of RRG distribution function makes it possible to study its dynamics through the Boltzmann equation, which in turn allows the introduction of collision terms in a natural way. For example, possible applications of RRG include self-interaction and dissipative dark matter, which can solve some small scale discrepancies between predictions of $\Lambda$CDM model and small scale observations [19, 14, 20, 21, 22, 23, 24, 25].

3.4 Gauge invariant equations

So far we are working with scalar perturbations in a Newtonian gauge. However, there are other different gauges each one with its own advantages. Thus, the ability to switch between different coordinate systems is useful. To facilitate this process we will rewrite all relevant equations using gauge invariant variables.

In general, the metric scalar perturbation can be written in terms of four quantities [16]

[TABLE]

where $\nabla^{2}\equiv\delta^{lm}\partial_{l}\partial_{m}$ and $\eta$ is the conformal time. In addition, with these four quantities, we can construct two invariant gauge variables,

[TABLE]

In a similar way, we can construct invariant gauge variables for density, velocity, anisotropic stress and pressure perturbations, so that we get:

[TABLE]

and

[TABLE]

respectively. Note that, the Newtonian gauge is realized choosing $w=\mu=0$. Thus, the equations obtained in previous sections are easily converted into gauge invariant equations performing the change $X_{N}\rightarrow X^{\left(gi\right)}$.

Therefore, the complete set of gauge invariant equations for RRG, takes the following form:

Einstein equations:

[TABLE]

and

[TABLE]

Conservation equations:

[TABLE]

and

[TABLE]

Boltzmann hierarchical equation:

[TABLE]

where prime indicates derivative with respect to conformal time $\eta$, $\mathcal{H}\equiv aH$ and $\delta_{\mathcal{R}}^{\left(gi\right)}\equiv\delta\rho_{\mathcal{R}}^{\left(gi\right)}/\bar{\rho}_{\mathcal{R}}$.

4 Applications and discussions

So far we have derived general equations for the cosmological evolution of RRG first-order perturbations. In the following, we consider an application of RRG as a model independent description of WDM and perform a pedagogical application.

We study the warm approximation of RRG, which can be used to study WDM and, based on current limits for DM mass, argue that, even in the case of very small thermal velocities, RRG is relativistic in the early universe. Then, we derive the initial conditions for RRG perturbations at this epoch. We also study the cosmological evolution of the gravitational potential in a toy model where the universe is solely composed by RRG.

4.1 RRG in the warm approximation

One of the most interesting applications of RRG model is for the description of WDM, which, for instance, can be keV sterile neutrinos [26]. If dark matter is cold, its simplest representation is a pressureless perfect fluid and the only relevant variables are density and fluid velocity. However, if dark matter has a small but non-negligible thermal velocity, it is important to study the impact of pressure and anisotropic stress. Such studies were conducted using Maxwell-Boltzmann distribution function, e.g., [27, 28]. Since the error of RRG model with respect to Maxwell-Boltzmann description for non-relativistic thermal velocities is essentially negligible, a warm fluid can also be studied in this framework.

Let us consider that the total energy density is given by a small increment above the rest energy density

[TABLE]

Using Eqs. (2.10) and (2.14) we can expand the fluid quantities, Eqs. (3.11), (3.12), (3.13) and (3.15) in terms of $x$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In order to describe a cold fluid we need to retain $x^{0}$ and $x^{1/2}$ orders, which are associated with density and fluid velocity, respectively. Beyond $x^{1/2}$ order, pressure and anisotropic stress are introduced with $x^{1}$ terms and energy density receives a correction of the same order. This is the first correction to the cold fluid description, which we call warm fluid, i.e., all terms up to $x^{1}$ order.

We can relate $x$ with the warmness parameter $b$. Using Eq. (2.14) we get

[TABLE]

The sound speed and $\omega_{\mathcal{R}}$ can also be written as a function of $x$:

[TABLE]

[TABLE]

Then, up to $x^{1}$ order we have

[TABLE]

Once warmness is introduced, in principle, we need to consider pressure, anisotropic stress and all the higher order multipoles in Boltzmann equations for RRG. Up to order $x^{1}$, these equations can be written as

[TABLE]

[TABLE]

[TABLE]

Note that we can not take $c_{s}\rightarrow 0$ because this corresponds to neglect the fluid velocity already at cold approximation.

However, the usual argument used for relativistic particles to neglect multipoles with $l\geq 2$ on scales larger then horizon, $k\eta\ll 1$, can be rephrased for warm RRG in terms of scales larger than the sound horizon $c_{s}k\eta\ll 1$, which is much smaller then horizon in the warm approximation. Therefore, pressure becomes important only on scales smaller then the sound horizon of the warm fluid, $c_{s}k\eta>1$ and, as usual, higher multipoles will be important on even smaller scales. This explains the fact that the perfect fluid approximation of RRG performed in [11] can reproduce with very good accuracy the transfer function for relics computed with the full set of Einstein-Boltzmann equations [12]. For instance, considering $b^{2}\sim 10^{-14}$ ($a_{\rm ref}=1$ assumed), which corresponds to particles with $m\sim\text{keV}$, pressure becomes important for scales roughly $k>10^{3}$ $\text{h}/\text{Mpc}$ at $z=0$. Assuming matter domination ($a\propto\eta^{2}$) from now until matter-radiation equality, $a_{\rm eq}=10^{-4}$, the scale at which pressure and higher multipoles become important increases to $k>10$ $\text{h}/\text{Mpc}$ until matter-radiation equality. During radiation dominated era ($a\propto\eta$), this scale remains unchanged as long as the warm approximation is valid.

Naturally, at even earlier times, the warm approximation breaks down and higher multipoles have to be considered. From Eq. (2.19) in the non-relativistic limit, we have $b^{2}\simeq 3(T_{0}/m)$, then, considering only the non-relativistic evolution of temperature, $T\propto a^{-2}$, the temperature and mass will be roughly the same at $a\sim b$, which can also be seen in Fig. 2. Hence, keV particles, which corresponds to $b^{2}=10^{-14}$, become relativistic at $a\sim 10^{-7}$ (curiously at this time the temperature of photons is roughly the same). In this case, the warm approximation breaks down very deep inside the radiation dominated era. Therefore, even WDM models that are very non-relativistic today are relativistic in the early universe. In the following, we consider this scenario in order to determine the initial conditions of RRG perturbations.

4.2 Initial conditions assuming adiabatic primordial modes

Let us consider adiabatic initial conditions for RRG perturbations in the early universe. As we saw, if we use RRG to describe particles with $m\simeq$ keV, they will be relativistic for $a<10^{-7}$. Typically, one is interested to determine the evolution of cosmological perturbations after neutrino decoupling, which corresponds to universe with $T<1$ MeV, or $a>10^{-10}$. In this scenario, we can safely consider RRG as a relativistic component. According to what was discussed in the previous section, this would not be case only for $b^{2}\lesssim 10^{-20}$, from Eq. (2.23) this corresponds to particles with $m\gtrsim 407$ keV. For simplicity, hereafter we drop the use of the superscript "$\left(gi\right)$".

In the described situation, RRG and neutrinos (see, e.g., Ref. [16]) will have analogous initial conditions. Therefore, in a super-horizon regime, we have

[TABLE]

and

[TABLE]

where $\zeta$ is the comoving curvature perturbation. In addition, the initial value for $\Phi$ and $\Psi$ is given by

[TABLE]

with

[TABLE]

where $\sum\bar{\rho}_{j}$ represents the sum of all non-interacting relativistic components and $\bar{\rho}_{\text{tot}}\simeq\bar{\rho}_{\gamma}+\sum\bar{\rho}_{j}$ in an era dominated by relativistic components.

For example, if one wants to describe standard model neutrinos with RRG in a universe with cold dark matter, then we have

[TABLE]

where $\bar{\rho}_{\text{tot}}\simeq\bar{\rho}_{\gamma}+\bar{\rho}_{\mathcal{R}}$ and $N_{\nu}^{\text{eff}}\equiv N_{\nu}g_{\nu}=3.0395$ [29] is the effective number of families. On the other hand, for the description of WDM with RRG we have

[TABLE]

where $\bar{\rho}_{\text{tot}}\simeq\bar{\rho}_{\gamma}+\bar{\rho}_{\nu}+\bar{\rho}_{\mathcal{R}}$ and $\bar{\rho}_{\mathcal{R}}/\bar{\rho}_{\gamma}$ is the initial energy density fraction of WDM relative to photons. In this case, the quantity $\bar{\rho}_{\mathcal{R}}/\bar{\rho}_{\gamma}$ contributes as a new relativistic species which can be determined by observational data.

4.3 Potential evolution

As discussed previously, RRG interpolates between a typical behavior of a UR component to a typically NR component (see Fig. 2). In this way, the RRG can be used to approximately describe the background where the relevant constituents are only UR and NR matter [4, 10]. In a similar context, we analyze the evolution of the gravitational potential $\Phi\left(k,a\right)$ in the linear level of perturbation, assuming a hypothetical universe constituted only by RRG with negligible anisotropic stress.

Neglecting the anisotropic stress ($\sigma_{\mathcal{R}}\simeq 0$), Eq. (3.44) leads to $\Phi=-\Psi$. Thus, combining Eqs. (3.42) and (3.45), we obtain a closed-form equation for the potential,

[TABLE]

where $\mathcal{H}$ can be obtained by manipulating Eqs. (2.16) and (2.13).

Applying the limit $a/a_{\text{ref}}\ll b$ in Eq. (4.19), we obtain

[TABLE]

whose solution is [15]

[TABLE]

This solution shows that in a super-horizon regime $k\eta\ll 1$, the potential $\Phi$ is constant, equal to its primordial value $\Phi_{\text{p}}$. However, when the modes enter the horizon $k\eta>1$, $\Phi$ amplitude rapidly decays with $\left(k\eta\right)^{-2}$.

Otherwise, applying the limit $a/a_{\text{ref}}\gg b$ in Eq. (4.19), we obtain

[TABLE]

whose general solution is $\Phi=C_{1}+C_{2}\left(k\eta\right)^{-5}$. Neglecting the decay mode, this solution shows that the gravitational potential is constant when $a/a_{\text{ref}}\gg b$. Therefore, we find that Eq. (4.20) and Eq. (4.22), as well as their solutions, are functionally the same as the expressions obtained for the gravitational potential $\Phi$ in an era dominated by the radiation and dust, respectively.

In order to obtain a complete description of the evolution of the potential, we define $a_{\text{ref}}=a_{0}\equiv 1$ and perform a change of variable from $\eta$ to $a$ in Eq. (4.19). Thus, we obtain

[TABLE]

where $\mathcal{H}$ can now be expressed as

[TABLE]

and the initial conditions given by $\Phi\left(a\rightarrow 0\right)=\Phi_{\text{p}}$ and $d\left[\Phi\left(a\rightarrow 0\right)\right]/da=0$ were assumed. In Fig. 5, we show the complete evolution of potential $\Phi\left(a\right)$ through numerical solutions of Eq. (4.23). This plot corroborates what we discussed previously, since, for $a\ll b$ (UR behavior), we have a constant potential and one which decays proportional to $\left(k\eta\right)^{-2}$, for the super-horizon ($k=10^{-3}$ Mpc*-1*) and sub-horizon ($k=1$ Mpc*-1*) modes, respectively. It is also evident from Fig. 5 that at $a\gg b$ (NR behavior), the potential is constant in $k$. However, when it crossing from region $a<b$ to $a>b$, the super-horizons modes ($k=10^{-3}$ Mpc*-1*) decrease by $1/10$, so that we have $\Phi\left(a\ll b\right)=(9/10)\Phi\left(a\gg b\right)$. In a realistic universe model, this decrease occurs in the equivalence era [15, 16]. For the hypothetical model considered here, this decrease occurs around $a=b$ (see dotted curves in Fig. 5), where RRG passes from an UR to a NR regime.

It is important to emphasize that, although this model well reproduces the predicted evolution of potential in realistic universes444In the absence of neutrinos., there are no physical reasons for neglecting the quadrupole term (related to anisotropic stress) in the sub-horizon regime since RRG is a non-interacting relativistic gas. Thus, only the results obtained for the super-horizon regime are physically consistent since in this regime we can neglect the quadrupole term. In realistic universe models, the quadrupole term of the photons component can by neglected while they are tightly coupled with the baryons.

5 Final comments

In this work we established a theoretical framework for cosmological linear perturbations of RRG described by a distribution function which encodes the main hypothesis of RRG, i.e., that all particles have the same momentum magnitude. We derived the Einstein-Boltzmann system of equations for the model, explicitly verifying that the non-relativistic and ultra-relativistic limits are recovered.

We also studied RRG in the warm approximation, showing that particles with $m\sim$ keV can be described as a perfect fluid on scales $k<10$ h/Mpc from now until, at least, $a_{\rm eq}$. For this mass scale, the warm approximation breaks up around $a\sim 10^{-7}$ and the full system of Einstein-Boltzmann equations have to be used. Our results corroborate previous ones [4, 11, 9] which state that RRG is a good approximation to Maxwell-Boltzmann statistics and, therefore, a sufficiently reliable approach to understand and study the behavior of WDM or neutrinos in a model-independent way, without the necessity of solving the full Einstein-Boltzmann system of equations.

Initial conditions for relativistic RRG in the early universe were also determined and we showed that the evolution of the gravitational potential through radiation and matter dominated eras can be described by a model with RRG only, recovering the general behavior of realistic models.

It is also interesting to note that the formal description of RRG in terms of a distribution function allows one to include other effects in the Boltzmann equation, such as collisional terms that can model self-interacting dark matter. These models became popular as a possible solution to the Core/Cusp and the Missing Satellites problems [30].

The simple equation of state of RRG motivates the analytical study of relativistic components in cosmological perturbative regime. In particular, it is interesting to look for analytical solutions of density perturbations in a model with WDM, radiation and $\Lambda$, where WDM is modeled by RRG. This kind of study involving generalizations of Mészáros equation [31] is being carried out by the authors.

Acknowledgments

G. Pordeus da Silva thanks CNPq-Brazil (141165/2017-0), government of the state of Paraíba-Brazil and CAPES-Brazil for the financial support.