FRW Models

This notebook numerically integrates the Friedman equation (18.78) to find a homogeneous isotropic cosmological model. The inputs are the four parameters that specify a FRW model —the Hubble constant (Subscript[H, 0 ])and the three Ω's for radiation, matter, and vacuum energy. The program assumes these parameters are such that the universe started with a big bang [cf Figure 18.10].

Clearing the variables used:

First,  clear all the variables that will be used in the calculation:

Clear[omr, omm, omv, omc, h, Th,  x, y, a0, t0, omcrit, yend, 
      ymax]

Parameters of a FRW model:

Four parameters specify a FRW cosmological model. First, there are the three Ω 's which are the ratios of the present density to critical density for radiation, matter, and vacuum energy. These dimensionless numbers are denoted by omr, omm, and oml respectively. These are specified in the following three statements:

omr = 0.00008

omm =  .3

omv =  .7

The last parameter is the Hubble constant (Subscript[H, 0]), or equivalently the Hubble time (Subscript[T, 0])=1/(Subscript[H, 0]) which we denote here by Th.   (Subscript[T, 0]) has the dimensions of time.  A billion years (Gyr)  is a convenient unit of time for cosmology and (Subscript[T, 0])=9.788 (h^(-1)) Gyr,  where  h = (Subscript[H, 0])/[100 (km/s)/Mpc]. Thus, for (Subscript[H, 0])=72(km/sec)/Mpc,

h =  .72

Th = 9.788/h

It is also convenient to define an (Subscript[Ω, c]) for curvature as in the text. Here we denote it by omc and its related to the other Ω's by:

omc := 1 - omr - omm - omv

For the parameters above

omc

Integrating the FRW equation in dimensionless variables:

Its convenient to rewrite the FRW equation in terms of  dimensionless variables as in(18.78). In the text we used  \(TraditionalForm\`t\&~\)=(Subscript[H, 0])t=t/(Subscript[T, h]), and \(TraditionalForm\`\(a\&~\)\)=a(t)/a((Subscript[t, 0])).  (The value of \(TraditionalForm\`a\&~\) at the present time (Subscript[t, 0]) is therefore always 1.)   Here its's notationally convnient to write x for \(TraditionalForm\`t\&~\), and y for \(TraditionalForm\`a\&~\). The FRW equation is then the same as that of a non-relativistic particle moving in a potential U(y) where y is the particle's position  and x is the time.

U[y_] := (1/2) ( -omr/y^2  - omm/y - omv y^2)

That is, the FRW equation (18.78) becomes:

       (((dy/dx))^2)= (Subscript[Ω, c])((Subscript[Ω, r]), (Subscript[Ω, m]),  (Subscript[Ω, v]))-2 U(y; (Subscript[Ω, r ]), (Subscript[Ω, m]), (Subscript[Ω, v]))  .

We plot this U(y) out to a value y=yend, also showing a horizontal line at the valueof (Subscript[Ω, c])/2. (You might have to change yend and edit  PlotRange in the Show  statement below to get aninformative plot.)

yend = 3

plotu := 
      Plot[U[y], {y, 0, yend}, AxesLabel -> {"y", "U"}, 
        DisplayFunction -> Identity]

plotomc := 
      Plot[omc/2,  {y, 0. , yend}, DisplayFunction -> Identity]

Show[plotu, plotomc, PlotRange -> {{0, yend}, {0, -3}}, 
      DisplayFunction -> $DisplayFunction]

To solve for a big bang FRW model (one that starts at a=0 or y=0)  is equivalent to carrying out the integral:

       x(y)=(∫_0^y) \(\([\(Ω\_c\)(\ Ω\_r, \ Ω\_m, Ω\_v) - 
                    2 \( U(w; \ Ω\_r, \ Ω\_m, \ Ω\_v)\)]\)\^\((\(-1\)/2)\)\)dw                                           (*)

x[y_] := 
      NIntegrate[(omc - 2  U[w])^(-1/2),  {w, 0, y}]

That's the answer. We will plot some specific cases in the next subsection.

For values of (Subscript[Ω, c]) above the maximum of the potential the expansion is unbounded in time;  for values that are below the expansion is bounded and the universe recollapses.  The value just at the maximum divides these two behaviors. We denote that value by omcrit.

crit :=  FindMinimum[-2 U[yy], {yy,   .5}]

omcrit := If[omv >  .000001, -crit[[1]],  0]

omcrit

For bounded cases the integration in (*) can extend only to the first turning point where the denominator vanishes.  This represents the expansion of a universe up to its maximum scale factor ymax.  Since the integral is singular at the turning point it is necessary to stop it just a bit before by small fraction ϵ .

eps =  .0001

The following list of statements looks for the turning point which is the smallest positive value of y where (Subscript[Ω, c])/2=U(tp). This is relevant only for closed models where (Subscript[Ω, c])/2 < (Subscript[U, max]).  (If the program fails to identify the corret turning point, its possible that you might have to insert a  statement ymax=  to give it the correct value and then run the program again.)

listtprules = NSolve[omc == 2  U[tp],  tp]

listtp = tp /. listtprules

listpostp = If[omc < omcrit, Select[listtp, # > 0  &], {}]

listpostp = Sort[listpostp]

If[omc < omcrit,  ymax := listpostp[[1]] (1 - eps), 
      ymax :=  yend]

ymax

In the unbounded cases the integration can extend to the limit of the plot,  yend.

Plotting the Results

A parameter z is introduced which runs from 0 to 1 along the curve to be plotted, either from the big bang to the big crunch for bounded models or from the big bang to yend for unbounded models.

ymax

yu[z_] := yend z

xu[z_] := x[yu[z]]

yb[z_] := 
      If[z <  .5, ymax (2 z), ymax (2 (1 - z))]

xb[z_] := If[z <  .5, x[yb[z]], 2  x[ymax] - x[yb[z]]]

yp[z_] := If[omc < omcrit,  yb[z], yu[z]]

xp[z_] := If[omc < omcrit,  xb[z], xu[z]]

plotxy := 
      ParametricPlot[{xp[z], yp[z]}, {z, 0, 1}, 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> Identity]

Show[plotxy, 
      Graphics[{Dashing[{ .01}],  Line[{{x[1], 0}, {x[1], 1}}]}], 
      DisplayFunction -> $DisplayFunction]

The dashed vertical line is a the time of today.

Rescaling to get a(t) as a function of t:

The Hubble time (Subscript[T, 0]) (Th) can be used to rescale the  dimensionless measure of the scale factor y and the dimensionless measure of the time x  that were convenient for computation to give the usual scale factor a(t) as follows:

t[z_] := Th xp[z]

a0 := (Th/Abs[omc])^(1/2)

a0 // N

a[z_] := a0 yp[z]

We can then plot a(t)  vs  t .   (Its the same plot as above but on rescaled axes.)

plotat := 
      ParametricPlot[{t[z], a[z]}, {z, 0, 1}, 
        AxesLabel -> {"t (Gyr)", "a(t)"}, 
        DisplayFunction -> Identity]

Show[plotat, 
      Graphics[{Dashing[{ .01}],  Line[{{Th x[1], 0}, {Th x[1], a0}}]}], 
      DisplayFunction -> $DisplayFunction]

The vertical dotted line is the present time. The present age (in Gyr)  can then be evaluated, and, if the evolution is bounded, the total age (in Gyr) and the maximum size (in billions of light years.)

agepresent = Th x[1]

If[omc < omcrit,  agebigcrunch = th xp[1]]

If[omc < omcrit,  maxradius = a0 ymax]