Solving Einsteins equations for the FLRW metric with ideal quantum gases for matter and radiation, assuming that pressure is proportional to the energy density $p=w\rho$ and deriving $w$ from statistical mechanics for the different particles boils down to the following function for the Hubble parameter:


Where $H_0\approx 68 \frac{km}{s\cdot Mpc}$ is the Hubble parameter value today, $\Omega_i=\frac{\rho_i}{\rho_{crit}}$ is the energy density ratio today for the different contributions ($M$ for matter, $R$ for radiation and $\Lambda$ for vacuum), and the critical energy density is given by $\rho_{crit}=\frac{3H_0^2c^2}{8\pi G}\approx \left(9.18\cdot 10^{-27}kg/m^3\right)c^2$. The individual powers of $a$ are all obtained from the general expression $\rho\sim a^{-3(w+1)}$ where $w$ are taken from equations of state in stat mech.

Now, apparently, since we can parametrize time evolution by $0<a<1$, the above equation provides us with time evolution for the energy density ratios. One could write:

$$\Omega_M(a)=\frac{\Omega_M}{a^3}\\\Omega_R(a)=\frac{\Omega_R}{a^4}\\\Omega_\Lambda(a)=\Omega_\Lambda$$ Then todays value is recovered by taking $\Omega_i(a=1)$. Now, lets do an experiment! Assume for todays values $\Omega_M\approx 0.3$ and $\Omega_\Lambda\approx 0.7$ (that would make $\Omega_R\approx 0$ today, since we know that these contributions have to sum to $1$). If we then consider the time evolution of these energy ratios, it appears that the radiation contribution stays zero at all times if it is zero today! But that cannot possibly be true, since we know that there was a time where radiation dominated and we have a CMB temperature of $2.73K$ remaining today! So how can the above equation still be valid and properly describe the time evolution of all density ratios? Maybe it has to be modified somehow? Did I miss something? Could someone clarify please? Thanks for any suggestion!

  • $\begingroup$ I don’t see how $\rho_i \sim a^{-3(w + 1)}$ and $\Omega_i(a) = \Omega_i a^{-3(w + 1)}$ can simultaneously be true. $\endgroup$
    – Socob
    Commented Jul 24, 2022 at 23:32

2 Answers 2


The only error is in taking $\Omega_\mathrm{R} \approx 0$ to be exactly the same as $\Omega_\mathrm{R} = 0$. The radiation density is not exactly $0$ today; it's just small enough to no longer matter.

However, as you go back in time ($a \to 0$), $\Omega_\mathrm{R}(a)$ grows faster than either $\Omega_\mathrm{M}(a)$ or $\Omega_\Lambda(a)$. Thus there is some $a$ for which $\Omega_\mathrm{R}(a)$ is larger than the other two. This occurs at very roughly $z = 3000$. Working backward, this corresponds to $a = 1/3000$, which means $\Omega_\mathrm{M}$ must be $3000$ times greater than $\Omega_\mathrm{R}$ today.

One way to get at the magnitude of $\Omega_\mathrm{R}$ is to look at the CMB as a blackbody photon gas, which has energy density per unit volume $$ u = \frac{8\pi^5(k_\mathrm{B}T)^4}{15c^3h^3}. $$ Plugging in $T = 2.73\ \mathrm{K}$ tells us $u/c^2 = 4.6\times10^{-31}\ \mathrm{kg/m^3}$ today, or $5\times10^{-5}$ times the value of $\rho_\mathrm{crit}$ you cite.


If you wanted your densities as functions of the scalefactor it would have been easy because there is an analytical solution, but if you need it by time it will be a little more tricky because you have to do some numerics. First you have to integrate the proper time by scalefactor if you want to take care of all Omegas:

$$t(a)=\frac{\int_0^{a'} \frac{1}{a \sqrt{\frac{\text{$\Omega $R}}{a^4}+\frac{\text{$\Omega $M}}{a^3}+\frac{\text{$\Omega $K}}{a^2}+\Omega \Lambda }} \, da}{\text{H0}}$$

and solve numerically for the scalefactor by time (Plot 1). Then you can solve for the density parameters. To get the percentage of an Omega, you have to divide it through the critical density of the corresponding time (Plot 2 and 3):

enter image description here

Sorry for not putting everything into Latex due lack of time.

For the derivation of the equation see the course Alan Guth @ MIT: The early universe (Youtube)

Edit: I put in the real values according to Planck 2013 Edit 2: minor corrections and generalizations

  • $\begingroup$ Wow, this is great! Also, I wonder how you put these black boxes around regions in Mathematica? $\endgroup$
    – Kagaratsch
    Commented Feb 10, 2015 at 23:53
  • 1
    $\begingroup$ Menu: "Palettes" -> "Classroom Assistant", then go to "Writing and Formatting", there you find "Frame". When you have your region selected (or select all) you can add all kind of frames. By the way, maybe the numerical values of the Ωs are not up to date, I did not take much care about the value of ΩK (the curvature density). Better check yourself if it is consistent with the Planck 2013 data. $\endgroup$
    – Yukterez
    Commented Feb 11, 2015 at 0:25
  • 1
    $\begingroup$ By the way, if you do more complex modells than FLRW/ΛCDM, such as phantom dark energy or dynamic equation of state etc you should add some NIntegrate[equation, {x, 0, X}, Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 100, Method -> "GaussKronrodRule"}, MaxRecursion -> 100] into your code. $\endgroup$
    – Yukterez
    Commented Feb 11, 2015 at 3:23
  • $\begingroup$ The plots look fine, so I must be missing something here: If $\omega_i$ are supposed to be the time-dependent density parameters, why do the formulae $\omega_i[\tau\_] {{}:}= \Omega_i \alpha[\tau]^{-n} / \rho_c [HF[\alpha[\tau]]]$ work? I would’ve expected $\omega_i[\tau\_] {{}:}= \Omega_i \alpha[\tau]^{-n} \rho_c[H0] / \rho_c [HF[\alpha[\tau]]]$ to get $0 \le \omega_i \le 1$ (how do the units even work out?). $\endgroup$
    – Socob
    Commented Jul 24, 2022 at 23:43
  • $\begingroup$ Can someone explain what is going on one the case of $\alpha$ in In[13]? $\endgroup$ Commented Jan 19 at 20:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.