Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matter density in closed/open universe

  1. Dec 31, 2009 #1

    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.

  2. jcsd
  3. Dec 31, 2009 #2


    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    Hi fab13, welcome to PF!

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

    [tex] \Omega_m - 1 = \frac{c^2}{H_0^2 \Re^2} \equiv \frac{c^2 \kappa}{H_0^2} [/tex]​

    where [itex] \Re [/itex] is the radius of curvature of the geometry of the universe at the present time, and [itex] \kappa [/itex] (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).

    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 [itex] \rho_m [/itex] and see how k responded. That is the beauty of the Friedman world models with [itex] \Lambda = 0 [/itex]. The total density of matter in the universe determines the spatial curvature. There is a one-to-one relationship between [itex] \Omega_m [/itex] and [itex] \kappa [/itex]. :wink:

    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:

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

    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 [itex] \rho_m [/itex] 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:


    I don't understand what you mean. Any value of [itex] \rho_m [/itex] 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, [itex] \kappa [/itex] will certainly get larger and larger as [itex] \rho_m [/itex] does, but its normalized version, k, will always be unity.

    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.
  4. Dec 31, 2009 #3

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

    [tex]\Omega_{m}-1=\frac{k}{H_{0}^{2} R_{0}^{2}}[/tex]

    with [tex]\Omega_{m}=\frac{\rho_{0}}{\rho_{crit}}[/tex]

    i have too : [tex]\rho_{crit}=\frac{3 H_{0}^{2}}{8 \pi G}[/tex]

    I assume [tex]\Omega_{m} > 1[/tex], so, k=1.

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

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

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

    [tex]\rho_{0}=2 \rho_{crit}=\rho_{crit} ( 1+\frac{1}{H_{0}^{2} R_{0}^{2}} )[/tex]

    [tex] \Rightarrow[/tex]

    [tex] R_{0}=\frac{1}{H_{0}} [/tex]

    Last edited by a moderator: Jan 1, 2010
  5. Dec 31, 2009 #4


    User Avatar
    Science Advisor
    Gold Member
    Dearly Missed

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

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

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

    If you will tell me 1/H0, 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.
  6. Jan 1, 2010 #5


    User Avatar
    Science Advisor
    Gold Member
    Dearly Missed

    fab, if someone asked me to calculate 1/H0 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 1017

    That is the number of seconds that is the reciprocal of H0, 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 1010 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)
  7. Jan 1, 2010 #6


    User Avatar
    Science Advisor
    Gold Member
    Dearly Missed

    Now in your posts you have several equations that suggest that you think that
    1/H0 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 3x108 (the speed of light in metric).

    So for example you have an equation about the radius of curvature of the universe, which yu call R0.
    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/H0.

    In your equations R0 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 R0 = 1/H0.

    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.
  8. Jan 1, 2010 #7


    User Avatar
    Science Advisor

    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 [tex]\rho_m[/tex]. This universe is also expanding at some rate at the current time. Let's call that rate [tex]H_0[/tex].

    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:

    [tex]\rho_c = \frac{3 H_0^2}{8\pi G}[/tex]
    [tex]\Omega_m = \frac{\rho_m}{\rho_c}[/tex]
    [tex]\Omega_k = -\frac{kc^2}{H_0^2}[/tex]

    then we can write the first Friedmann equation as follows:

    [tex]H^2(a) = H_0^2 \left(\frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2}\right)[/tex]

    In this equation, there are only two free parameters: the expansion rate ([tex]H_0[/tex]) and the matter density ([tex]\Omega_m[/tex]). In principle, either of these two parameters can take on any values at all (well, provided the matter density is positive). The remaining parameter, [tex]\Omega_k[/tex], 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.
  9. Jan 1, 2010 #8


    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    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 R0 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:

    [tex] \Omega_m - 1 = \frac{c^2}{H_0^2 \Re^2} \equiv \frac{c^2 \kappa}{H_0^2}[/tex]​

    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 R0 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:

    [tex] \Omega_m - 1 = \frac{c^2 \kappa}{H_0^2} = \frac{k}{H_0^2 R_0^2}[/tex]

    [tex] \frac{c^2 \kappa}{H_0^2} = \frac{\kappa}{|\kappa|H_0^2 R_0^2}[/tex]

    [tex] c^2 = \frac{1}{|\kappa|R_0^2} \Rightarrow R_0^2 = \frac{1}{|\kappa|c^2}[/tex]​

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

    [tex] 1 \, \, \textrm{s}^2 = \frac{1}{|\kappa|c^2} [/tex]

    [tex] c^2(1 \, \, \textrm{s}^2) = \frac{1}{|\kappa|} \equiv \Re^2 [/tex]

    [tex] \Rightarrow \Re = c(1 \, \, \textrm{s}) [/tex]

    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 ~ 1035). 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 R0=1 is not reasonable, because 1 light-second is a very small radius of curvature for the universe (same result as before).
  10. Jan 1, 2010 #9
    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 :

    [tex]\Lambda=\lambda*3H_{0}^{2}=0[/tex] and k=0 [tex]\Rightarrow \rho_{0}=\rho_{crit}[/tex] (see blue line)
    [tex]\lambda=0.7[/tex] and [tex]\Omega_{m}=0.3[/tex] and k=0 [tex]\Rightarrow \rho_{0}=0.3* \rho_{crit}[/tex] (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 [tex] \Omega_{m} > 1[/tex] and k=1 with [tex]\Lambda=0[/tex].
    This case is represented on green line. I chose [tex]\Omega_{m}=2[/tex].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 [tex] c^{2}[/tex] factor to [tex]\Omega_{k}[/tex]. I decided to put the "meter" dimension to [tex]R(t)[/tex] according to your advice. So , [tex]\Omega_{k}=-\frac{k c^{2}}{H_{0}^{2} R_{0}^{2}}[/tex] is dimensionless like [tex]\Omega_{m}[/tex].

    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 :


    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 :


    Here's the matlab code part for k=1 and [tex]\rho_{m}=2 \rho_{crit}[/tex] :

    Code (Text):

    function friedmann


    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);


    function dydt3 = syseq3(t,y)
    rho_0=2*rho_crit;                 % Omega_m=2
    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=[ y(2) ;                        
    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.


    Attached Files:

    Last edited: Jan 1, 2010
  11. Jan 1, 2010 #10


    User Avatar
    Science Advisor

    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:

    [tex]\Omega_m + \Omega_k = 1[/tex].

    You can, of course, use the definitions of [tex]\Omega_m[/tex] and [tex]\Omega_k[/tex] in order to determine what the value of [tex]k[/tex] must be for your choice of [tex]\rho_m[/tex], but that seems like unnecessary work to me.
  12. Jan 1, 2010 #11
    For the case k=1, i determine automatically [tex]R_{0}[/tex] which appears only one time in my differential system :

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

    I have a strange result if i take [tex]\Omega_{m} \geq 3[/tex]. Indeed, the scale factor ( in fact, [tex]a(t)={\frac{R(t)}{R_{0}}[/tex]) 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 [tex] 1 \leq \Omega_{m} \leq 3[/tex], the scale factor always increases without asymptote. Why does this result occur from 3 ? a relationship with speed of light ?


    Attached Files:

  13. Jan 1, 2010 #12


    User Avatar
    Science Advisor

    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 [tex]\Omega_k = -1[/tex]). I suggest you use [tex]\Omega_k[/tex] instead. The value of [tex]\Omega_k[/tex] is exactly specified by the value you choose for [tex]\Omega_m[/tex], using the constraint I wrote previously.
  14. Jan 2, 2010 #13

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

    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:

    Code (Text):

    dydt3=[ y(2) ;                        

    I have always got the same result for [tex]1 < \Omega_{m} < 3[/tex], no collapse. For [tex] \Omega_{m} \geq 3 [/tex], 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 ?


    Attached Files:

  15. Jan 3, 2010 #14


    User Avatar
    Science Advisor

    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:

    [tex]H^2(a) = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} \right)[/tex]

    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:

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

    [tex]\left( \frac{da}{dt} \right)^2 = H_0^2 \left( \frac{\Omega_m}{a} + \Omega_k \right)[/tex]

    [tex]\frac{da}{dt} = \pm H_0 \sqrt{\frac{\Omega_m}{a} + \Omega_k}[/tex]

    Whether or not there is a recollapse depends entirely upon the argument under the square root. If [tex]\Omega_k < 0[/tex], then the argument under the square root will become negative for [tex]a > -\Omega_k/\Omega_m[/tex]. 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 [tex]\Omega_m = 1[/tex] is interfering with the code's ability to extrapolate beyond [tex]\Omega_m = 3[/tex]. 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.

    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:

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

    With the first Friedmann equation, including the cosmological constant, becoming:
    [tex]H^2(a) = H_0^2 \left( \frac{\Omega_m}{a^3} + \frac{\Omega_k}{a^2} + \Omega_\Lambda \right)[/tex]
    Last edited: Jan 3, 2010
  16. Jan 3, 2010 #15


    User Avatar
    Science Advisor

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

    Just in case you weren't sure...
  17. Jan 3, 2010 #16

    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, [tex] \frac{da}{dt} [/tex] is positive at the beginning and vanishes for :


    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 :


    Code (Text):

    dydt3=[ y(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 [tex]\Omega_{m}=1.5[/tex], [tex]\Omega_{k}=-0.5[/tex] and [tex]\Omega_{\Lambda}=0[/tex] 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 [tex]a=-\frac{\Omega_{m}}{\Omega_{k}}=3[/tex] in this case. Now, i have an infinite expansion for [tex]\Omega_{m} > 1[/tex], 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.


    Attached Files:

    Last edited: Jan 3, 2010
  18. Jan 4, 2010 #17

    i have reformulated the differential system through your equation :

    [tex]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 :

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

    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} [/tex] ​

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

    [tex] a_{max} = a ( - \frac{\Omega_{m}}{\Omega_{k}} )[/tex]

    [tex] t_{max} = - \frac{5}{-4} = 1.25 [/tex] on the figure 5 ​

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

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

    respecting so equation n°1 for the value of [tex] a'_{0} [/tex]. For this case, the Matlab solver begins integration at [tex] t_{0}=380.000 years [/tex], 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 ( [tex] a(t_{final})=0 [/tex].

    Thanks a lot.

    Attached Files:

    Last edited: Jan 4, 2010
  19. Jan 4, 2010 #18


    User Avatar
    Science Advisor

    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 [tex]H_0 = H(a=1)[/tex] always?
  20. Jan 4, 2010 #19
    for the last results, i took [tex] p = 0 [/tex] 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) :

    [tex]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} [/tex]

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

    Attached Files:

    Last edited: Jan 5, 2010
  21. Jan 5, 2010 #20


    User Avatar
    Science Advisor

    Maybe I see the problem.

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

    Looks like you're missing a factor of a in this equation, as H = a'/a.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook