Matter density in closed/open universe

[fab13]

hello,

Friedmann equations lead to with cosmological constant =0 :

Omega_m -1= k/(H0^2 * R0^2) where H0 is Hubble constant and R0 the radius equal to 1 here.

If Omega_m > 1, we have necessairly k=1 but rho_m is determined by :

rho_m= rho_crit *(1+k/(H0^2 * R0^2)) (eq n°1)


I mean we cannot choose any value of rho_m, only one which is above (eq n°1) .

Moreover, if we take rho_crit= 5 protons mass and H0= 71 km/s/Mpc ( H0=71000/(3*10^(22)) s^(-1)) , rho_m takes a very big value with k=1.

Same case for Omega_m < 1, eq n°1 gives a negative value for rho_m with k=-1 and HO=71000(3*10^(22)) s^(-1).

From i know, k is an integer (-1,0,1), so, how can we choose k=1 and a free value for rho_m, and in the same time, verify the equation n°1 above ? the question is the same for k=-1.

Regards

 

[cepheid]

Hi fab13, welcome to PF!

Friedmann equations lead to with cosmological constant =0 :

Omega_m -1= k/(H0^2 * R0^2) where H0 is Hubble constant and R0 the radius equal to 1 here.

That seems about right. Using the version of the Friedman equation that I have, I get a result of:

$$\Omega_m - 1 = \frac{c^2}{H_0^2 \Re^2} \equiv \frac{c^2 \kappa}{H_0^2}$$​

where $\Re$ is the radius of curvature of the geometry of the universe at the present time, and $\kappa$ (kappa), the inverse of its square, is the curvature of the universe. I am assuming that this is the same as your result, except that you are using dimensionless versions of the curvature (k instead of kappa) and your radius of curvature is dimensionless (mine is measured in metres).

rho_m= rho_crit *(1+k/(H0^2 * R0^2)) (eq n°1)

I mean we cannot choose any value of rho_m, only one which is above (eq n°1) .

Yes eqn. 1 must be true, but that doesn't mean that the total density of matter in the universe is constrained to satisfy a certain value of k. If anything, it's the other way around. The total density of matter in the universe is something that is determined observationally (i.e. it is a fact of nature!). It, in turn, determines k. So, to see what different possibilities there are, we would set different values of $\rho_m$ and see how k responded. That is the beauty of the Friedman world models with $\Lambda = 0$. The total density of matter in the universe determines the spatial curvature. There is a one-to-one relationship between $\Omega_m$ and $\kappa$.

Moreover, if we take rho_crit= 5 protons mass and H0= 71 km/s/Mpc (

This doesn't make sense, for two reasons. The first reason is you can't choose the critical density and the Hubble constant independently of each other. They are dependent upon each other:

$$\rho_c = \frac{3H_0^2}{8 \pi G} = 1.88 \times 10^{-26} h^2\, \, \textrm{kg m}^{-3}$$​

where h is the dimensionless Hubble constant. That leads to the second reason why what you said doesn't make sense. The proton mass is a mass, but $\rho_m$ is a density. Even if you meant five proton masses per cubic metre, this value is still not quite right for h = 0.71, as shown here:

http://www.google.com/search?hl=en&...88e-26+kg*m^-3)+/+(proton+mass)&aq=f&oq=&aqi=

H0=71000/(3*10^(22)) s^(-1)) , rho_m takes a very big value with k=1.

I don't understand what you mean. Any value of $\rho_m$ that is greater than the critical density will result in k = 1. There must be something wrong with your version of the formula or with the calculations you did. Of course, $\kappa$ will certainly get larger and larger as $\rho_m$ does, but its normalized version, k, will always be unity.

so, how can we choose k=1 and a free value for rho_m, and in the same time, verify the equation n°1 above =? the question is the same for k=-1.

Again, any value of the matter density greater than the critical density will result in k = 1, and any value of the matter density less than the critical density will result in k = -1.

 

[fab13]

hi,

thanks for your help, it seems you have not understood my problem. I take the following equation :

$$\Omega_{m}-1=\frac{k}{H_{0}^{2} R_{0}^{2}}$$

with $$\Omega_{m}=\frac{\rho_{0}}{\rho_{crit}}$$

i have too : $$\rho_{crit}=\frac{3 H_{0}^{2}}{8 \pi G}$$

I assume $$\Omega_{m} > 1$$, so, k=1.

I can write : $$\rho_{0}=\rho_{crit} ( 1+\frac{1}{H_{0}^{2} R_{0}^{2}} )$$ (eq n°1)

The first problem is $$\rho_{0}$$ has a very big value if i take $$H_{0}=71 km/s/Mpc$$
and $$R_{0}=1$$. ( i get $$\rho_{0}=1.785*10^{35} * \rho_{crit}$$ ) .

But the main issue is that, for $$R_{0}$$ and $$H_{0}$$ given,there is only one value of $$\rho_{0}$$. I would like for example take $$\rho_{0}=2 * \rho_{crit}$$ with k=1 and verify in the same time eq n°1. I thought maybe we ca adjust $$R_{0}$$ to confirm this equation. We would have then in this case :

$$\rho_{0}=2 \rho_{crit}=\rho_{crit} ( 1+\frac{1}{H_{0}^{2} R_{0}^{2}} )$$

$$\Rightarrow$$

$$R_{0}=\frac{1}{H_{0}}$$

Regards

 

[marcus]

fab, you seem to be using conventional units. you have not set the speed of light c = 1, right?

In that case H₀ is a reciprocal time. It is often easier to work with if you treat it that way.

What do you get for 1/H₀ ? You should get somewhat over 13 billion years, or the equivalent number of seconds.

If you will tell me 1/H₀, in terms of some conventional time unit, like seconds (or years), I might be able to help sort things out for you. Put that into your formulas instead of the expression involving Megaparsecs and kilometers.

 

[marcus]

fab, if someone asked me to calculate 1/H₀ and I didn't know already what it was, I think I would type this into google

1 megaparsec/(71 km)

when you press return it doesn't do a search, it calculates with its built-in google calculator which knows various constants and units.
It will come back with something like 4.3 x 10¹⁷

That is the number of seconds that is the reciprocal of H₀, the socalled "Hubble time".

But it would be less lazy to type in this, and press return:

(1 megaparsec/(71 km)) seconds

Then it will come back with something like 1.377 x 10¹⁰ years.

That is the same as 13.77 billion years.

(the google calculator likes to express very long spans of time in terms of years instead of seconds)

 

[marcus]

Now in your posts you have several equations that suggest that you think that
1/H₀ is not a time, but is instead a distance.

That is the sort of thing that happens when cosmologists use special units in which c = 1.

The danger there is that you will copy one of their formulas which should have a c in it but where they have neglected to show the c, because they set it equal to one.
Then if you calculate with that formula, naively using ordinary metric units, you will get screwed up and numbers will come out way too big or too small by factors like 3x10⁸ (the speed of light in metric).

So for example you have an equation about the radius of curvature of the universe, which yu call R₀.
Don't confuse that with the Hubble Radius, which is not curvature radius at all, but quite different. The Hubble radius or Hubble distance is the distance from us of galaxies which are receding at exactly the rate c, the speed of light.
The Hubble distance = c/H₀.

In your equations R₀ is not the Hubble distance, it seems to be the radius of curvature of the universe, which is spatially closed in this case.

But now you write an equation R₀ = 1/H₀.

How can that be right!? You are equating a distance with a time. the radius of curvature is a distance, and the reciprocal Hubble is a time. In fact something like 13.7 billion years.

Something is fishy. The units don't match up.

 

[Chalnoth]

Let me see if I can clear things up a bit. Let's take a universe with some average density. Let's call that density $$\rho_m$$. This universe is also expanding at some rate at the current time. Let's call that rate $$H_0$$.

Bear in mind that there is nothing in physics that tells us what the relationship between those two quantities must be. We could have a very large density and a very small expansion rate (in which case we should expect a quick recollapse), or we might have a very low density and a very large expansion rate (in which case we should expect continual expansion).

So how do the Friedmann equations come in? They tell us how these two quantities change with time. If we define the critical density and density fractions as follows:

$$\rho_c = \frac{3 H_0^2}{8\pi G}$$
$$\Omega_m = \frac{\rho_m}{\rho_c}$$
$$\Omega_k = -\frac{kc^2}{H_0^2}$$

then we can write the first Friedmann equation as follows:

$$H^2(a) = H_0^2 \left(\frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2}\right)$$

In this equation, there are only two free parameters: the expansion rate ($$H_0$$) and the matter density ($$\Omega_m$$). In principle, either of these two parameters can take on any values at all (well, provided the matter density is positive). The remaining parameter, $$\Omega_k$$, shows how these other two affect the spatial curvature.

In particular, if it just so happens that the "critical density" determined by the expansion rate exactly equals the matter density, then we have zero spatial curvature. If it is expanding more slowly, then we have a closed universe that recollapses. If it is expanding more quickly, then we have an open universe that expands forever.

 

[cepheid]

In your equations R₀ is not the Hubble distance, it seems to be the radius of curvature of the universe, which is spatially closed in this case.

But now you write an equation R₀ = 1/H₀.

How can that be right!? You are equating a distance with a time. the radius of curvature is a distance, and the reciprocal Hubble is a time. In fact something like 13.7 billion years.

Something is fishy. The units don't match up.

Hey marcus,

I agree with you in general that the OP's problem is that he's not keeping proper track of units. However, I think that his R₀ really is measured in seconds, if you stick to SI units. (If you use a unit system where c=1, it doesn't matter). Bear with me here. I'm going to run something by you. Let me know if you agree. We know that the following equation I posted is definitely true:

$$\Omega_m - 1 = \frac{c^2}{H_0^2 \Re^2} \equiv \frac{c^2 \kappa}{H_0^2}$$​

where the symbols I'm using are defined in my first post (kappa is the spatial curvature). Now the question is, given the OP's form of the equation, how must R₀ be defined? Well, first of all, it seems pretty obvious that k = sgn(κ), where sgn() is the sign function. In other words, k = κ/|κ| if κ is non-zero, and k = 0 if κ is zero. Therefore, assuming κ is non-zero, and equating his version of the equation with mine, we have:

$$\Omega_m - 1 = \frac{c^2 \kappa}{H_0^2} = \frac{k}{H_0^2 R_0^2}$$

$$\frac{c^2 \kappa}{H_0^2} = \frac{\kappa}{|\kappa|H_0^2 R_0^2}$$

$$c^2 = \frac{1}{|\kappa|R_0^2} \Rightarrow R_0^2 = \frac{1}{|\kappa|c^2}$$​

So, I am asserting that this is how the OP's R₀ must be defined. Note to the OP: it is NOT dimensionless. From this formula, R₀² has units of (m^-2m²s^-2)^-1 = s². Therefore, when you say R₀ = 1, you must specify 1 WHAT? Setting R₀ = 1 s means that you are saying:

$$1 \, \, \textrm{s}^2 = \frac{1}{|\kappa|c^2}$$

$$c^2(1 \, \, \textrm{s}^2) = \frac{1}{|\kappa|} \equiv \Re^2$$

$$\Rightarrow \Re = c(1 \, \, \textrm{s})$$
​

Congratulations (to the OP). You picked a radius of curvature for the universe equal to c*1s which is ~ 300,000 km. Your radius of curvature is less than the distance between the Earth and the Moon. No wonder you need a ridiculous density of matter to pull that off (Ω_m ~ 10³⁵). Even if, instead of using SI units, you use a unit system in which c = 1 (which seems to be what your version of the formula was intended for), you can still see that setting R₀=1 is not reasonable, because 1 light-second is a very small radius of curvature for the universe (same result as before).

 

[fab13]

thanks for your answers, i begin to see more clearly. I posted here because I'm making in Matlab the simulation of the scale factor evolution in time, in the 3 cases that friedmann equations allow. (k=-1,0,1).
I managed to get the good curves ( see figure 1 in attachment) for :

$$\Lambda=\lambda*3H_{0}^{2}=0$$ and k=0 $$\Rightarrow \rho_{0}=\rho_{crit}$$ (see blue line)
$$\lambda=0.7$$ and $$\Omega_{m}=0.3$$ and k=0 $$\Rightarrow \rho_{0}=0.3* \rho_{crit}$$ (see red line)​

We can see on red line an accelerated expansion which is the expected result in this case.

Now, i'd like to get the recollapse case where $$\Omega_{m} > 1$$ and k=1 with $$\Lambda=0$$.
This case is represented on green line. I chose $$\Omega_{m}=2$$.You may notice there's not a recollapse of the universe with my code, the scale factor always increases, that's not normal. I use the MATLAB solver "ode45" to resolve the following differential system (it's based on Runge-Kutta algorithm):

eq number 1 :
$$\frac{R^{''}}{R}=-\frac{4\,\pi\,G}{3c^{2}}(\rho \, c^{2}+3\,p)+\frac{\Lambda}{3}$$

eq number 2 :
$$\bigg(\frac{R^{'}}{R}\bigg)^{2}=\frac{8 \pi G \rho}{3}+\frac{\Lambda}{3}-\frac{k c^{2}}{R^{2}}$$
​

You can see i added $$c^{2}$$ factor to $$\Omega_{k}$$. I decided to put the "meter" dimension to $$R(t)$$ according to your advice. So , $$\Omega_{k}=-\frac{k c^{2}}{H_{0}^{2} R_{0}^{2}}$$ is dimensionless like $$\Omega_{m}$$.

we can rewrite the equation number 1 in the form :

eq number 3 :
$$R^{''}=\bigg(\frac{-4 \pi G}{3c^{2}}(\rho c^{2}+3 p)+\frac{\Lambda}{3}\bigg) R^{'} \bigg(\frac{8 \pi G \rho}{3}+\frac{\Lambda}{3}-\frac{k c^{2}}{R^{2}}\bigg)^{-1/2}$$
​

with matter density :

eq number 4 :
$$\rho=\rho_{0} \bigg(\frac{R_{0}}{R(t)}\bigg)^{3}$$
​

In my source code, I decided to work with the following variable :

$$a(t)=\frac{R(t)}{R_{0}}$$​

wich gives :

eq number 5 :
$$a{''}=\bigg(\frac{-4 \pi G}{3c^{2}}(\frac{\rho_{0} c^{2}}{a^{3}}+3 p)+\frac{\Lambda}{3}\bigg) a^{'} \bigg(\frac{8 \pi G}{3} \frac{\rho_{0}}{a^{3}}+\frac{\Lambda}{3}-\frac{k c^{2}}{R_{0}^{2} a^{2}}\bigg)^{-1/2}$$
​

Finally, we have the following differential system :

y1'=a'=y2
y2'=a"

Here's the MATLAB code part for k=1 and $$\rho_{m}=2 \rho_{crit}$$ :

function friedmann

H0=71000/(3*10^(22));
cv=3*10^(8);

t_begin=380000*(3600*24*365); % time where uiverse become transparent
t_final=20*10^(9)*(3600*24*365); % 20 Gyr

init=[ 0.001 ; 10^(3/2)*H0 ]; % redshift z=1000=a_0 ; a'_0=10^(3/2)*H0[t3,y3]=ode45(@syseq3,[t_begin t_final],init);

figure(1);
t3=t3/(3600*24*365);
plot(t3,y3(:,1),'g');
end

function dydt3 = syseq3(t,y) 
H0=71000/(3*10^(22));
p=0;
k=1;
cv=3*10^(8);
G=6.673*10^(-11);
rho_crit=3*H0^(2)/(8*pi*G);
rho_0=2*rho_crit;                 % Omega_m=2
Lambda_red=0;                       
Lambda=Lambda_red*3*H0^(2); % Cosmological constant=0
R0=cv/H0;       % got by the relationship between Omega_k and Omega_m (Omega_m=2 here)

 dydt3=zeros(2,1);
 dydt3=[ y(2) ;                        
       (-(4*pi*G)/(3*cv^(2))*((rho_0/y(1)^(3))*cv^(2)+3*p)+Lambda/3)*(y(2))*((8*pi*G*rho_0/y(1)^(3))+Lambda/3-k*(cv/(R0*y(1)))^(2))^(-1/2);          
       ];
end

I'd like to get the orange line on "size_universe.jpg" (in attachment) with the famous "recollapse".
Do you know what's wrong in my modelisation ?

Sorry for this long post but i wanted to be as explicit as possible.

Regards.

 

[Chalnoth]

I haven't looked into your code in detail, but I suspect that the problem is a problem of units. The determination that k = 1, 0, or -1 is a particular choice of units, and may be messing you up. Basically, the reason why some theorists choose to set the curvature parameter in these terms is that if the curvature is postive, negative, or zero, you can always change your "reference time" so that k = 1, 0, or -1.

However, with the usual formulation of the Friedmann equations, we explicitly fix the scale factor to be equal to 1 at the current time, so this degree of freedom is taken up, and we have to allow k to take on any floating-point value.

Instead of working in this manner, though, I would recommend simply using the density fraction parameterization I give above, where you have the constraint:

$$\Omega_m + \Omega_k = 1$$.

You can, of course, use the definitions of $$\Omega_m$$ and $$\Omega_k$$ in order to determine what the value of $$k$$ must be for your choice of $$\rho_m$$, but that seems like unnecessary work to me.

 

[fab13]

For the case k=1, i determine automatically $$R_{0}$$ which appears only one time in my differential system :

$$R_{0}=\frac{c}{H_{0} (\Omega_{m}-1)^{1/2}}$$​

I have a strange result if i take $$\Omega_{m} \geq 3$$. Indeed, the scale factor ( in fact, $$a(t)={\frac{R(t)}{R_{0}}$$) increases initially and
stabilizes after. See on figure2.png in attachment. Actually, in the case of the recollapse, the problem may be "stiff". I tried others MATLAB solver for stiff problems like "ode15s", "ode23s", "ode23tb" but always the same behaviour, the scale factor follows an asymptote and doesn't decrease as it's expected. For $$1 \leq \Omega_{m} \leq 3$$, the scale factor always increases without asymptote. Why does this result occur from 3 ? a relationship with speed of light ?

Thanks.

 

[Chalnoth]

For the case k=1, i determine automatically $$R_{0}$$ which appears only one time in my differential system :

$$R_{0}=\frac{c}{H_{0} (\Omega_{m}-1)^{1/2}}$$​

I have a strange result if i take $$\Omega_{m} \geq 3$$. Indeed, the scale factor ( in fact, $$a(t)={\frac{R(t)}{R_{0}}$$) increases initially and
stabilizes after. See on figure2.png in attachment. Actually, in the case of the recollapse, the problem may be "stiff". I tried others MATLAB solver for stiff problems like "ode15s", "ode23s", "ode23tb" but always the same behaviour, the scale factor follows an asymptote and doesn't decrease as it's expected. For $$1 \leq \Omega_{m} \leq 3$$, the scale factor always increases without asymptote. Why does this result occur from 3 ? a relationship with speed of light ?

Thanks.

As I said, don't use the integer values for k. You have to fix other things in order to get sensible equations if you use k = 1, 0, -1 (also, k = 1 does not mean $$\Omega_k = -1$$). I suggest you use $$\Omega_k$$ instead. The value of $$\Omega_k$$ is exactly specified by the value you choose for $$\Omega_m$$, using the constraint I wrote previously.

 

[fab13]

hi,

i took your advice considering now in my differential system $$\Omega_{k}=-\frac{k c^2}{H_{0}^2 R_{0}^2}$$ with the constraint :

$$\Omega_{k}+\Omega_{m}=1$$​

i get so the following differential equation :

$$a{''}=\bigg(\frac{-4 \pi G}{3c^{2}}(\frac{\rho_{0} c^{2}}{a^{3}}+3 p)+\frac{\Lambda}{3}\bigg) a^{'} \bigg(\frac{8 \pi G}{3} \frac{\rho_{0}}{a^{3}}+\frac{\Lambda}{3}+\frac{\Omega_{k} H_{0}^{2} }{a^{2}}\bigg)^{-1/2}$$
​

here's the associate code part:

dydt3=[ y(2) ;                        
       (-((4*pi*G)/(3*cv^(2)))*((rho_0/y(1)^(3))*cv^(2)+3*p)+Lambda/3)*(y(2))*((8*pi*G*rho_0/y(1)^(3))+Lambda/3+(Omega_k*(H0/y(1))^(2)))^(-1/2);          
       ];

I have always got the same result for $$1 < \Omega_{m} < 3$$, no collapse. For $$\Omega_{m} \geq 3$$, the problem seems to be stiff, see the result on figure 3 in attachment ( procuced by "ode45" Matlab solver).

someone has already encountered this problem in differential system solving ?
is there a solution to this issue ?

Thanks

 

[Chalnoth]

What are you using for p?

Anyway, you should be able to do a similar thing as follows. First, using the parametrization I used above:

$$H^2(a) = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} \right)$$

Note that I no longer even need to worry about things like factors of G and pi and whatnot: everything is just in terms of density fraction and expansion rate, dramatically simplifying things.

Then, I can go about solving the differential equation by using the definition of the expansion:

$$\left( \frac{1}{a} \frac{da}{dt} \right)^2 = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} \right)$$

$$\left( \frac{da}{dt} \right)^2 = H_0^2 \left( \frac{\Omega_m}{a} + \Omega_k \right)$$

$$\frac{da}{dt} = \pm H_0 \sqrt{\frac{\Omega_m}{a} + \Omega_k}$$

Whether or not there is a recollapse depends entirely upon the argument under the square root. If $$\Omega_k < 0$$, then the argument under the square root will become negative for $$a > -\Omega_k/\Omega_m$$. Thus there is necessarily a maximum value for the scale factor a in this case.

A second essential feature of this formula is that it is symmetric: the expansion rate for a given density is the same whether that rate is positive or negative. This means that the recollapse is perfectly symmetric with respect to the prior expansion.

My suspicion is that when you're using this solver, you're running into some sort of radius of convergence issue where the value at $$\Omega_m = 1$$ is interfering with the code's ability to extrapolate beyond $$\Omega_m = 3$$. I bet that if you don't try to go beyond the point where the recollapse begins (because you know it's going to be symmetric with respect to the prior expansion), you won't run into any problems.

Edit:
Actually, the more I think about it, the more I think that you've done something incorrect with the pressure. The pressure p is a function of what sort of matter you have. For normal matter, it's zero (as far as cosmology is concerned).

Edit 2:
Just in case it helps, here's the acceleration equation assuming normal matter and cosmological constant only, using the parameterization above:

$$\dot{H} + H^2 = \frac{\ddot{a}}{a} = -\frac{H_0^2}{2}\left(\frac{\Omega_m}{a^3} - 2\Omega_\Lambda \right)$$

With the first Friedmann equation, including the cosmological constant, becoming:
$$H^2(a) = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} + \Omega_\Lambda \right)$$

 

[Chalnoth]

I figure there is one final thing I should mention. With this choice of variables,
$$H_0 \equiv H(a = 1)$$.

Just in case you weren't sure...

 

[fab13]

Hello,

I'm on the right track, the simulation gives coherent results. I think the problem of "no recollapse" is the change of sign in :

$$\frac{da}{dt} = \pm H_0 \sqrt{\frac{\Omega_m}{a} + \Omega_k}$$

Indeed, $$\frac{da}{dt}$$ is positive at the beginning and vanishes for :

$$a=-\frac{\Omega_{m}}{\Omega_{k}}$$

So, we've got to report this change of sign in the differential equation system. I used the "sign" command of Matlab in this way :

y1=a(t)=R(t)/R0
y2=a'(t)

dydt3=[ y(2) ;                        
       (-((4*pi*G)/(3*cv^(2)))*((rho_0/y(1)^(3))*cv^(2)+3*p)+Lambda/3)*y(2)/H0*sign((Omega_m/y(1)^(3))+Lambda_red+(Omega_k/y(1)^(2)))*(abs((Omega_m/y(1)^(3))+Lambda_red+(Omega_k/y(1)^(2))))^(-1/2);          
       ];

You ca see i put : " sign((Omega_m/y(1)^(3))+Lambda_red+(Omega_k/y(1)^(2)))" in front of the absolute value of "(Omega_m/y(1)^(3))+Lambda_red+(Omega_k/y(1)^(2))))^(-1/2)". Well, i respect the positive value for the square root.

Then, i get for $$\Omega_{m}=1.5$$, $$\Omega_{k}=-0.5$$ and $$\Omega_{\Lambda}=0$$ the figure number 4. I said the results are coherent because there's a break on the curve for a=3, which is expected because a' vanishes for $$a=-\frac{\Omega_{m}}{\Omega_{k}}=3$$ in this case. Now, i have an infinite expansion for $$\Omega_{m} > 1$$, no more problem with the previous post concerning factor 3.

From we can see on figure 4, i failed to tell the solver that there's a change of sign. Anyone could see a solution to this ? I thought "sign" MATLAB command would be the solution but after vanishing of a'(t), a(t) increases once again like we notice it on figure 4.

Thanks.

 

[fab13]

Hello,

i have reformulated the differential system through your equation :

$$H^2(a) = \bigg(\frac{1}{a} \frac{da}{dt}\bigg)^{2} = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} + \Omega_\Lambda \right)$$
​

which gives equation n°1 :

$$\left( \frac{da}{dt} \right)^2 = H_0^2 \left( \frac{\Omega_m}{a} + \Omega_k + \Omega_\Lambda a^{2} \right)$$
​

then we have :

$$a = \Omega_{m} \bigg( \frac{1}{H_{0}^{2}} \bigg(\frac{da}{dt}\bigg)^{2} - \Omega_{k} - \Omega_{\Lambda} a^{2} \bigg)^{-1}$$​

So, the differential equation can be written in this way :

$$a''=\bigg(\frac{-4 \pi G}{3c^{2}}(\frac{\rho_{0} c^{2}}{a^{3}}+3 p)+\frac{\Lambda}{3}\bigg) \Omega_{m} \bigg( \frac{1}{H_{0}^{2}} \bigg(\frac{da}{dt}\bigg)^{2} - \Omega_{k} - \Omega_{\Lambda} a^{2} \bigg)^{-1}$$​

Finally, i managed to get the curve representing the "recollapse" with $$\Omega_{m} > 1$$, $$a_{0} = 1$$, and $$a'_{0} = H_{0}$$ and $$\Omega_{\Lambda}=0$$ (see figure 5 in attachment) . The integration begins at $$t_{0}=13.7 Gyr$$. The results are validated by the value of the maximum of scale factor before the recollapse :

$$a_{max} = a ( - \frac{\Omega_{m}}{\Omega_{k}} )$$

$$t_{max} = - \frac{5}{-4} = 1.25$$ on the figure 5​

I would like to get the first part of this curve ( for $$t0 << 13.7 Gyr$$ ), ie from different initial conditions. I tried with these conditions :

$$a_{0} = 0.001$$, and $$a'_{0} =H_{0} (5*10^{3}-4)^{1/2} = 70.6824 H_{0}$$

respecting so equation n°1 for the value of $$a'_{0}$$. For this case, the Matlab solver begins integration at $$t_{0}=380.000 years$$, so a redshift z=1000.

I get the figure 6. There's well the recollapse design but the maximum is not equal to 1.25 (much more little, about 0.09 ) and the big crunch occurs fastier than on figure 5.

What's the problem, How it is done ? Why i don't get the same recollapse as in figure 5, i mean the same maximum and time "t_final" of big crunch ( $$a(t_{final})=0$$.

Thanks a lot.

 

[Chalnoth]

I don't think you've made proper use of the acceleration equation. What are you using for p?

Edit: Also, have you made certain that $$H_0 = H(a=1)$$ always?

 

[fab13]

for the last results, i took $$p = 0$$ but i think i should consider maybe the radiation pressure because i search to get scale factor evolution from a early age of the universe, where radiation was important.

in my code, the differential equation is expressed as follows ( see equation.png in attachment) :

$$a"=\bigg( \frac{-H_{0}^{2}}{2} ((\frac{\Omega_{m}}{a^{3}}) - 2 \Omega_{\Lambda}) \bigg) \Omega_{m} \bigg( \frac{1}{H_{0}^{2}} \bigg(\frac{da}{dt}\bigg)^{2} - \Omega_{k} - \Omega_{\Lambda} a^{2} \bigg)^{-1}$$

The problem is when i take $$a_{0} = 0.001$$ and $$a'_{0} = H_{0} (5*10^{3}-4)^{1/2} = 70.6824 H_{0}$$ ($$\Omega_{m} = 5$$), i am not sure if your condition , ie $$a_{now} = 1$$ and $$a'_{now} = H_{0}$$ is respected because i don't start the integration at $$t = t_{now}$$ but before, at $$t = 380000 years$$. I don't know how to constrain the equation n° 1 to this condition while taking initial conditions at a time younger.

 

[Chalnoth]

for the last results, i took $$p = 0$$ but i think i should consider maybe the radiation pressure because i search to get scale factor evolution from a early age of the universe, where radiation was important.

in my code, the differential equation is expressed as follows ( see equation.png in attachment) :

$$a"=\bigg( \frac{-H_{0}^{2}}{2} ((\frac{\Omega_{m}}{a^{3}}) - 2 \Omega_{\Lambda}) \bigg) \Omega_{m} \bigg( \frac{1}{H_{0}^{2}} \bigg(\frac{da}{dt}\bigg)^{2} - \Omega_{k} - \Omega_{\Lambda} a^{2} \bigg)^{-1}$$

The problem is when i take $$a_{0} = 0.001$$ and $$a'_{0} = H_{0} (5*10^{3}-4)^{1/2} = 70.6824 H_{0}$$ ($$\Omega_{m} = 5$$), i am not sure if your condition , ie $$a_{now} = 1$$ and $$a'_{now} = H_{0}$$ is respected because i don't start the integration at $$t = t_{now}$$ but before, at $$t = 380000 years$$. I don't know how to constrain the equation n° 1 to this condition while taking initial conditions at a time younger.

Maybe I see the problem.

$$a'_{0} = H_{0} (5*10^{3}-4)^{1/2} = 70.6824 H_{0}$$

Looks like you're missing a factor of a in this equation, as H = a'/a.

 

[fab13]

Hello,

the problem persists. Everything is OK for $$a_{0} = 1$$ and $$\big(\frac{da}{dt}\big)_{0} = H_{0}$$ (see https://www.physicsforums.com/attachment.php?attachmentid=22861&d=1262647111"

" Maybe I see the problem.

LaTeX Code: asingle-quote_{0} = H_{0} (5*10^{3}-4)^{1/2} = 70.6824 H_{0}

Looks like you're missing a factor of a in this equation, as H = a'/a. "​

If i take $$a = 0.001$$, $$\big(\frac{da}{dt}\big)_{0}$$ is given by the following formula :

$$\left( \frac{da}{dt} \right) = H_{0} \left( \frac{\Omega_m}{a} + \Omega_k + \Omega_\Lambda a^{2} \right)^{1/2}$$​

So, in my example, i have :

$$\Omega_{m} = 5$$, $$\Omega_{\Lambda} = 0$$ and $$\Omega_{k} = -4$$​

which gives:

$$\big(\frac{da}{dt}\big)_{0} = H_{0} (5*10^{3}-4)^{1/2} = 70.6824\,\,H_{0}$$​

and with this, the result of the plot is wrong , see https://www.physicsforums.com/attachment.php?attachmentid=22862&d=1262647111"
Indeed, the maximum of scale factor is equal to about 0.9 instead of $$a_{max} = - \frac{\Omega_{m}}{\Omega_{k}} = \frac{5}{4} = 1.25$$ here.

I put my Matlab code :

function friedmann

H0=71000/(3*10^(22));
cv=3*10^(8);

t_begin=380000*(3600*24*365);
t_final=40*10^(9)*(3600*24*365);
Omega_m=5;
Omega_red=0;
Omega_k=1-Omega_m-Omega_red;
a_0=0.001;          % value of a_0
a_prim_0=H0*(Omega_m/a_0+Omega_k)^(1/2);    %value of a'_0
init=[ a_0 ; a_prim_0 ];


[t3,y3]=ode45(@syseq3,[t_begin t_final],init,options);
t3=t3/(3600*24*365);   %time in year
plot(t3,y3(:,1),'b');

end

function dydt3 = syseq3(t,y) 

H0=71000/(3*10^(22));
cv=3*10^(8);
G=6.673*10^(-11);
cv=3*10^(8);
G=6.673*10^(-11);
rho_crit=3*H0^(2)/(8*pi*G);
Omega=0;
rho_0=5*rho_crit;
Omega_m=rho_0/rho_crit;
Omega_red=Omega/3*H0^2;
Omega_k=1-Omega_m-Omega_red;

 dydt3=zeros(2,1);
 dydt3=[ y(2) ;                        
 (-(H0^(2)/2)*(Omega_m/y(1)^(3)-2*Omega_red))*Omega_m*(1/(H0^(2))*y(2)^2-Omega_k-Omega_red*y(1)^(2))^(-1);                
       ];


end

As you can seen, with a redshift equal to 1000 ( because a_0 = 0.001), i begin the integration at t_begin ~ 380000 year.

I don't understand why the maximum is not the same than with $$a_{0} = 1$$ and $$a'_{0} = H0$$. If i take a value between "1 and 0.001" , "0.1" for example, the maximum "1.25" is nearly reached (see figure 7 in attachment, the max is equal to 1.24 ) .

I really don't know why this max is modified by other initial conditions than $$a_{0} = 1$$ and $$a'_{0} = H0$$. The equations seems to be correct, the code too. It comes maybe from the time at which i begin integration. This time is a function of the initial value of $$a_{0}$$. For example, i have for $$a_{0} = 0.001$$ a time corresponding to about 380000 years, when the universe became transparent.

however, The main goal is reached even if i have difficulties with earlier ages.

Thanks.

 

[Chalnoth]

If i take $$a = 0.001$$, $$\big(\frac{da}{dt}\big)_{0}$$ is given by the following formula :

$$\left( \frac{da}{dt} \right) = H_{0} \left( \frac{\Omega_m}{a} + \Omega_k + \Omega_\Lambda a^{2} \right)^{1/2}$$​

Oh, that's wrong. The correct formula is:
$$\left( \frac{1}{a}\frac{da}{dt} \right) = H_{0} \left( \frac{\Omega_m}{a} + \Omega_k + \Omega_\Lambda a^{2} \right)^{1/2}$$

 

[fab13]

i don't think, you wrote the following equation previously :

$$\left( \frac{1}{a} \frac{da}{dt} \right)^2 = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} \right)$$
​

 

[Chalnoth]

i don't think, you wrote the following equation previously :

$$\left( \frac{1}{a} \frac{da}{dt} \right)^2 = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} \right)$$
​

Oh, hmm, perhaps I am mistaken. I forgot that you'd already multiplied the a through. Okay.

I do think, however, that the problem is something very simple like this. You're still using the energy density rho in your equations, though. You don't have to. You might be able to fix the problem by changing completely to the much simpler formulation that only makes use of density fractions, and not using the constants G and c at all.

 

[fab13]

hello, the problem is solved, i had just to add a very small "relative " and "absolute " error as parameters for Matlab "ode" solver. I get the three models (open, closed and euclidean) as it's represented on figure in attachement.

Thanks for all.

 

[Chalnoth]

hello, the problem is solved, i had just to add a very small "relative " and "absolute " error as parameters for Matlab "ode" solver. I get the three models (open, closed and euclidean) as it's represented on figure in attachement.

Thanks for all.

Ah, excellent :) Sorry it was such a pain in the butt.