FRW Models

This notebook numerically integrates the modified FRW equation 



with



to compute a specific homogeneous isotropic cosmological model for the inputs (1) the Hubble constant (Subscript[H, 0 ]), and the three Ω's for radiation, matter, and cosmological constant. Note bene: The program assumes these parameters are such that the universe started with a big bang, which we know may not always be the case.  But it appears to be for our universe.

Clearing the variables used:

Clear all the variables that will be used in the calculation:

Clear[Ωr,Ωm,ΩΛ,Ωk,h,Th, x,y,a0,t0,Ωcrit,yend,ymax]

Parameters of a FRW model:

Four parameters specify a FRW cosmological model. First, there are the three Ω 's :

Ωr=.1*10^(-5)

Ωm=.3

ΩΛ=.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])=67.3(km/sec)/Mpc,

h=.7

Th=9.788/h

We define (Subscript[Ω, k]) for curvature. It is related to the other Ω's by:

Ωk:=1-Ωr-Ωm-ΩΛ

For the parameters above

Ωk

Integrating the FRW equation in dimensionless variables:

Its convenient to rewrite the FRW equation in terms of  dimensionless variables as above.  Here it is convenient to write x for (OverTilde[t]), and y for (OverTilde[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)( -Ωr/y^2 -Ωm/y-ΩΛ y^2)

That is, the FRW equation becomes:

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

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

yend=4

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

plotomc:=Plot[Ωk/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) (([(Subscript[Ω, k])( (Subscript[Ω, r]), (Subscript[Ω, m]),(Subscript[Ω, v]))-2U(w; (Subscript[Ω, r]), (Subscript[Ω, m]), (Subscript[Ω, v]))])^((-1/2)))dw                                           (*)

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

That's the answer.

For values of (Subscript[Ω, k]) 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 Ωcrit.

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

Ωcrit:=If[ΩΛ>.000001,-crit[[1]], 0]

Ωcrit

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[Ω, k])/2=U(tp). This is relevant only for closed models where (Subscript[Ω, k])/2 < (Subscript[U, max]).  (If the program fails to identify the correct 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[Ωk==2U[tp], tp]

listtp=tp/.listtprules

listpostp=If[Ωk<Ωcrit,Select[listtp,#>0 &],{}]

listpostp=Sort[listpostp]

If[Ωk<Ωcrit, 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]],2x[ymax]-x[yb[z]]]

yp[z_]:=If[Ωk<Ωcrit, yb[z],yu[z]]

xp[z_]:=If[Ωk<Ωcrit, 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[Ωk])^(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 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[Ωk<Ωcrit, agebigcrunch=th xp[1]]

If[Ωk<Ωcrit, maxradius=a0 ymax]