Plot Structures for Friedman-Robertson-Walker metric, k=+1

Before anything else, we know that we will need some special routines that Maple uses for plotting solutions of

ordinary differential equations which are not always in memory; therefore, we must load them in:

> with(plots):

A few numbers that will be useful later on are the following:

cMpc is the speed of light in Megaparsecs per eon (an eon is a billion years)

gMpc is the factor to multiply by to change 1/(Megaparsecs)^2 to grams/cm^3

> cMpc:=306.937; gMpc:= 1/(7.0554*10^20);

We may now begin by stating the values of the constants that will be used here;

those values are explained and/or justified in a different place. We do this one with constants so that

the lifecylce of the universe is periodic, or, if you prefer, of finite length in time.

We have need for constants, A, B, k, and Lambda; we denote the value for B^2 by B2 and the value for

Lambda/3 by Lambda3. A is in Mpc, B2 in Mpc squared and Lambda in inverse Mpc squared.

> A:=10772; B2:=24240; k:=1; Lambda3:=1/8226^2;

We may then write down the fourth-order polynomial appropriate for this problem, where I use the

symbol x for the usual arc-time variable, eta:

> R1:=A*R(x)+B2-k*(R(x))^2-Lambda3*(R(x))^4;

and then we may use that to describe the differential equations that need to be solved. However, because

of the fact that the principal equation has the derivative squared, there are actually two branches for the

square root involved, and we must deal with each of them separately:

> sys1:=diff(R(x),x)=R1^(0.5), diff(t(x),x)=R(x);

> sys2:=diff(R(x),x)=-R1^(0.5), diff(t(x),x)=R(x);

To actually solve the differential equations we must append initial conditions, which are of course different for the

two cases. In the first case, when the function R=R(x) has positive derivative, the universe is beginning, so we

begin with R(0)=0 and t(0)=0:

> rw1:=dsolve({sys1,R(0)=0,t(0)=0},[R(x),t(x)],type=numeric):

For the second half, when the derivative is negative, we want to join the two halves at the correct point.

Therefore, we need to find out exactly where that point is. First we find the place where the derivative vanishes,

since that is the turn-around point, i.e., the point at which they will be joined. The value for R at that point is

fairly easy to come by; it is simply the appropriate (real, positive) root of the polynomial:

> r1:=subs(R(x)=r,R1);

We now plot this function to get a feel for its shape, knowing that we want it to begin at R=0 and go until the polynomial

has its first, real root. We change the number 7000 in the plot command below two or three times until we get it right.

> plot(r1,r=0..7000);

Then we come along and get a slightly more accurate value for the location of that root.

> rot:=fsolve(r1,r);

> Rmax:= rot[2];

To determine the corresponding value for eta, i.e., for x, we go ahead and create the data structure for a plot, and

look in detail at its entries. We do this with the command just below, but, first, with a semi-colon at the end of

the command. The semi-colon, and the fact that it has been given a name, i.e., RW1, causes Maple to print out

a list of all 500 of its computations for the values of (x,R(x)). One then reads forward to the place where the

square root vanishes, and then goes imaginary. At that point the list of points in RW1 reads [FAIL,FAIL] from

then on---since it doesn't plot complex numbers---and so we can pick off the appropriate value of eta = x,

which turns out to be 2.122. After looking at that printout, and picking off that value, we then change the

semi-colon to a colon, so that the printout is ``retracted.'' The numpoints comment is simply to put in enough

points to have a reasonable accuracy. As well, looking at the one for t=t(x), we see that the corresponding

value for t is 6539.497, which is translated, just below, in billions of years, namely 21.3 eons.

> etamax:=2.122; tmax:=6539.497/cMpc;

> RW1:=odeplot(rw1,[x,R(x)],0..2.2,color=red,numpoints=500):

> RWt1:=odeplot(rw1,[x,t(x)],0..2.2,color=green,numpoints=500):

> display([RW1,RWt1], title=`R(eta) and t(eta), first half`);

To find the density at maximum expansion, we ignore the radiation density, since we know it has been getting

smaller faster, and calculate the matter density at that time from its R-dependence: the first number is in

inverse Megaparsecs squared, the second is in grams/cubic centimeter.

> rhomax:=(3/evalf(8*Pi))*A/Rmax^3; Rhomax:=rhomax*gMpc;

The above graph gives us an idea what this first half looks like, with the red curve being R, and the green one t.

Yet a different alternative is to actually plot R versus t:

> RWrt1:=odeplot(rw1,[t(x),R(x)],0..2.2,numpoints=500,title=`R=R(t)`):

> display(RWrt1,title=`R=R(t), first half`);

We are also asked for various other quantities that relate to the first half of the age of the universe:

The first one of these is the time when the matter and radiation densities were the same, i.e., when

A/R^3 = B2/R^4, which means that R was equal to B2/A. The answer for R is given in Megaparsecs.

Looking that value up in the tables created above give us the corresponding values of eta=x and of time.

The first answer for time is also in Megaparsecs, but it is then changed to years by dividing by the speed of

light and multiplying by 1 billion:

> Requal:=evalf(B2/A); etaequal:=.01197; Tequal:=.01312; tequal:=Tequal*10^9/cMpc;

We also want to know about NOW. An input in my approach is that the current ``radius'' of the universe is

given by R0. That input allows us to consult the tables above and determine the values of eta and time that

correspond to now; that last ones being, of course, the current ``age'' of the universe, given first in Megaparsecs

and then in billions of years, i.e., eons.

> R0:=4500; eta0:=1.40; T0:=2341; t0:=T0/cMpc;

We also want to know the current value of the Hubble constant, in km/sec/Megaparsec

the current rate of increase of the size of the universe, i.e., dR/dt, in Megaparsecs/eon,

do note that this number is greater than the speed of light measured locally,

and the current value for q0, which is the second derivative, appropriately normalized:

q0 = -(d^2R/dt^2)(R/H0^2), which is dimensionless:

> H0:=sqrt(evalf(A/R0^3+B2/R0^4 -k/R0^2-Lambda3))*(2.9979*10^5);

> Rdot:=H0*R0*cMpc/(2.9979*10^5);

> qnum:=A/(2*R)+B2/R^2+Lambda3*R^2;

> qnum0:=evalf(subs(R=R0,qnum));

> qdenom:=r1/r^2;

> qdenom0:=evalf(subs(r=R0,qdenom));

> q0:=qnum0/qdenom0;

We also want to know about the source that is observed on Earth today with a redshift of Z=5. This implies

immediately that it was emitted when R=R(t) was R0/(Z+1) = 750 Mpc. Again looking at the tables, this

gives us a value for eta of 0.506 and for t of 134.32 Mpc, which is 0.438 billions of years ago.

> RZ:=R0/6; etaZ:=0.506; TZ:=134.32; tZ:=TZ/cMpc;

To determine the current metric distance of that source from us, we multiply the current value of R(t),

namely R0, by the difference in eta between the emission event and the reception event, the result being

given in Megaparsecs. We see that this source is now 89% of the distance R0 away from us.

> LZ:=R0*(eta0-etaZ); fract:=LZ/R0;

Lastly we are asked about the particle horizon. There we pretend that light could now be received that was

originally emitted at (or almost at) the time of the Big Bang, so that it has had an arc-time, eta, to travel of the

current value of eta, called eta0 above. The current metric distance to a galaxy at that horizon is then given by

the following calculation, where we see that it is 140% of the distance R0 away from us, at this time.

> LPart:=R0*(eta0-0); fractPart:=LPart/R0;

Now, we need to arrange for the "rest" of the graph, the other half that is. We need the opposite sign for the

square root, as given above in sys2, and to begin at x = x0=2.1222, R at its appropriate value, and also t:

> rw2:=dsolve({sys2,R(2.1222)=6575,t(2.1222)=6539.5},[R(x),t(x)],type=numeric):

As before, we are unsure just how far to plot, so we first begin at the left-hand end and guess a number about twice

that value, and then look at the output, to determine when R goes back to 0, and therefore the "big" crash has come,

and we don't have any need to go any further; also as before, after we obtain that information, we tell Maple NOT to \

print it out for us in the final version by changing the semi-colon to a colon. We find that x=4.243 is approximately

where R goes negative, so that this is the end. The corresponding time is 13086.16.

> RW2:=odeplot(rw2,[x,R(x)],2.1222..4.25,color=red,numpoints=500):

> RWt2:=odeplot(rw2,[x,t(x)],2.1222..4.25,color=green,numpoints=500):

> RWrt2:=odeplot(rw2,[t(x),R(x)],2.1222..4.25,color=blue,numpoints=500):

> display({RW1,RW2},title=`complete R(eta)`);

> display({RWt1,RWt2},title=`complete t(eta)`);

> display({RWrt1,RWrt2},title=`complete R=R(t)`);

>