Undergrad A glitch in Jorrie’s Cosmo-Calculator?

[pbuk]

Common cosmological notation is $[something] = \frac{H\left(z\right)^2}{H_0^2} =: E\left(z\right)^2$.

Excellent, thank you, so adapting that slightly I can write:
$$ E(s)^2 = \Omega_{[m,0]} s^3 + \Omega_{[rad,0]} s^4 + \Omega_{[\Lambda,0]} s^{3(1+w)} + \Omega_{[K,0]} s^2 $$
[The last term in the equation above was corrected from $\Omega_{[K,0]} s^3$.]

[code lang=javascript title="Which I can implement for w = -1 as"]
const getESquaredAtS = (s) => {
const s2 = s * s;
return omegaM0 * s2 * s + omegaLambda0 + omegaRad0 * s2 * s2 + OmegaK * s2;
};
[/code]
.. and that saves a LOT of calculating square roots and then squaring, converting to different units and then back again.

 

[Jorrie]

Excellent, thank you, so adapting that slightly I can write:
$$ E(s)^2 = \Omega_{[m,0]} s^3 + \Omega_{[rad,0]} s^4 + \Omega_{[\Lambda,0]} s^{3(1+w)} + \Omega_{[K,0]} s^3 $$

I think the last term should be $\Omega_{[K,0]} s^2$, because curvature dilates by an inverse square law of the expansion factor.
I see in the code piece you have it right...

 

[Jorrie]

I am not sure which equations you mean. $\Omega = \Omega_r + \Omega_m + \Omega_\Lambda$ is always true, as is $\Omega + \Omega_k = 1$. If $\Omega = 1$, then $\Omega_k =0$, but if f $\Omega \neq 1$, then $\Omega_k \neq 0$.

Did you have other equations in mind?

No, I just wanted confirmation of my intuition, but I still need to think how to implement it in code.
Perhaps @pbuk has an implementation in mind.

 

[pbuk]

I see in the code piece you have it right...

Yes, $\LaTeX$ copy and paste error, thanks.

 

[George Jones]

No, I just wanted confirmation of my intuition, but I still need to think how to implement it in code.
Perhaps @pbuk has an implementation in mind.

I might be able to help, but I do not understand the notation. I have figured out "s" for "stretch". To start, I would like to ask about stuff in the following post.

It is calculated in the main html program, not in Calculate.js

line 356 var OmegaM = (Omega-Ynow * Ynow / Yinf / Yinf) * s_eq / (1+s_eq);

I'm not sure if this should have been "z_eq / (1+z_eq) "
Also take note that it is expressly rounded to 4 decimals for that calculation.
When calculated later in the tables, it is rounded to the requested number pf decimals.

Ps: The densities for the tables are calculated in
Calculation.js line 27: var Om = (Omega - Ol) * (s_eq + 1) / (s_eq + 2);
The (s_eq + 1) / (s_eq + 2) is a mistake that has been prviously pinted out and corrected in Lightcone8 by @pbuk. There is some further processing in lines 243 to 250.
These line numbers are only relevant to the LightCone7-2021-03-12 Version

Does OmegaM represent $\Omega_{m,0}$ or $\Omega_m \left(t\right)$?

What are Ynow and Yinf?

 

[pbuk]

Perhaps @pbuk has an implementation in mind.

I do: by eliminating unnecessary conversions I have got the calculations down to this:
[code lang=javascript]
const getParamsAtStretch = (s: number): LcdmModelVariables => {
const eSquared = getESquaredAtStretch(s);
const s2 = s * s;
const h = h0 * Math.sqrt(eSquared);
const omegaM = (omegaM0 * s2 * s) / eSquared;
const omegaLambda = omegaLambda0 / eSquared;
const omegaRad = (omegaRad0 * s2 * s2) / eSquared;
return {
h,
omegaM,
omegaLambda,
omegaRad,
temperature: cmbTemperature * s,
rhoCrit: rhoCrit0 * eSquared,
};
};
[/code]
... and this:
[code lang=javascript]
for (let i = integrationResults.length - 1; i >= 0; --i) {
const { s, t, dNow: dUnsafe, dPar } = integrationResults;

const params = model.getParamsAtStretch(s);
const hGy = params.h * model.kmsmpscToGyr;

// Force dNow to zero at zero redshift.
const dNow = s === 1 ? 0 : dUnsafe;

results.push({
z: s - 1,
a: 1 / s,
s,
t,
dNow,
d: dNow / s,
r: 1 / hGy,
dPar,
vGen: params.h / (s * model.h0),
vNow: dNow * model.h0Gy,
v: (dNow * hGy) / s,
...params,
});
}
[/code]

I only have one variable I don't quite understand - I've called it vGen, it's called XDpar in the LightCone7 code but displayed in the table as $V_{gen} / c$ and has the description

Is there anything snappier/more descriptive than vGen?

 

[pbuk]

Does OmegaM represent $\Omega_{m,0}$ or $\Omega_m \left(t\right)$?

What are Ynow and Yinf?

Yes, good questions, but the work is done moving away from that code now as in my last post. You can see the work in progress at https://github.com/light-cone-calc/light-cone-calc/blob/refactor/src/model.ts.

My next focus is on the calculation of the current density parameters, here is what we have now ($\Omega, \Omega_{[\Lambda,0]}$ and $s_{eq} = z_{eq} + 1$ are provided as parameters e.g. from Planck 2015).

  const rhoCrit0 = rhoConst * h0Seconds * h0Seconds;
  const omegaM0 = ((omega - omegaLambda0) * seq) / (seq + 1);
  const omegaRad0 = omegaM0 / seq;
  const OmegaK0 = 1 - omegaM0 - omegaRad0 - omegaLambda0;

 

[George Jones]

My next focus is on the calculation of the current density parameters, here is what we have now ($\Omega, \Omega_{[\Lambda,0]}$ and $s_{eq} = z_{eq} + 1$ are provided as parameters e.g. from Planck 2015).

$\Omega$ or $\Omega_0$? If $\Omega_0$, then this is equivalent to taking $\Omega_{r,0}$, $\Omega_{m,0}$, and $\Omega_{\Lambda,0}$ as input parameters.

  const rhoCrit0 = rhoConst * h0Seconds * h0Seconds;
  const omegaM0 = ((omega - omegaLambda0) * seq) / (seq + 1);
  const omegaRad0 = omegaM0 / seq;
  const OmegaK0 = 1 - omegaM0 - omegaRad0 - omegaLambda0;

At the risk of outlining what folks already know ... radiation-matter equilibrium happens when $\rho_r = \rho_m$, or, equivalently, when
$$\begin{align}
\Omega_{r,eq} &= \Omega_{m,eq} \\
\Omega_{r,0} s_{eq}^4 &= \Omega_{m,0} s_{eq}^3 \\
\Omega_{r,0} &= \Omega_{m,0} / s_{eq},
\end{align}$$
which is the second-last line in the above code.

To get the second line in the above code, use $\Omega = \Omega_r + \Omega_m + \Omega_\Lambda$ to substitute for $\Omega_{r,0}$.
$$\begin{align}
\Omega_{r,0} &= \Omega_{m,0} / s_{eq} \\
\Omega_0 - \Omega_{m,0} - \Omega_{\Lambda,0} &= \Omega_{m,0} / s_{eq}. \\
\end{align}$$

Solving for $\Omega_{m,0}$ gives the second line of the code.

 

[Jorrie]

Does OmegaM represent $\Omega_{m,0}$ or $\Omega_m \left(t\right)$?

What are Ynow and Yinf?

OmegaM represents $\Omega_{m,0}$. Later on OmegaMatterT represents $\Omega_m \left(t\right)$

Ynow and Yinf are inverse Hubble values, i.e. Hubble times now and at t -> infinity
These still come from the late Marcus's inputs.

 

[Jorrie]

I only have one variable I don't quite understand - I've called it vGen, it's called XDpar in the LightCone7 code but displayed in the table as $V_{gen} / c$ and has the description
View attachment 304033
Is there anything snappier/more descriptive than vGen?

AFAICR, the late Marcus and I were thinking of the worldlines of generic particles that presently cross our Hubble sphere. The variable XDpar appeared in the discussion as the ratio Vthen/Vnow and it was specifically included in the columns for purposes of plotting that interesting wordline mentioned. This is how Vgen[eric] came to be. It shows "dramatically" how a particle could enter the Hubble sphere early on and then be accelerated to cross it again, in the generic case, right now.
We used it to work out the crossing and inflection points.

 

[pbuk]

$\Omega$ or $\Omega_0$? If $\Omega_0$, then this is equivalent to taking $\Omega_{r,0}$, $\Omega_{m,0}$, and $\Omega_{\Lambda,0}$ as input parameters.

Good catch - yes this should be $\Omega_0$, I will change the variable name.

At the risk of outlining what folks already know
...
Solving for $\Omega_{m,0}$ gives the second line of the code.

Many thanks, it's good to have it confirmed.

 

[George Jones]

OmegaM represents $\Omega_{m,0}$. Later on OmegaMatterT represents $\Omega_m \left(t\right)$

Ynow and Yinf are inverse Hubble values, i.e. Hubble times now and at t -> infinity
These still come from the late Marcus's inputs.

Then, in

line 356 var OmegaM = (Omega-Ynow * Ynow / Yinf / Yinf) * s_eq / (1+s_eq);

$$\frac{Y_{now}^2}{Y_\infty ^2} = \frac{H_\infty ^2}{H_0^2} = \Omega_{r,0} s^4 + \Omega_{m,0} s^3 + \Omega_{\Lambda,0} + \Omega_{k,0} s^2$$
when $w = -1$ for dark energy.

If $R$ is used for the universe's scale factor, then $s := R_0 /R$. In universes close to our own, $R \rightarrow \infty$ as $t \rightarrow \infty$. Consequently, $s \rightarrow 0$ and $H\infty ^2 /H_0^2 \rightarrow \Omega_{\Lambda,0}$ as $t \rightarrow \infty$.

Thus, in universes close to our own,

line 356 var OmegaM = (Omega-Ynow * Ynow / Yinf / Yinf) * s_eq / (1+s_eq);

is equivalent to

const omegaM0 = ((omega - omegaLambda0) * seq) / (seq + 1);

The version with $\Omega_{\Lambda,0}$ should also work in universes that are not close to our own.

 

[pbuk]

I think I'm all done now, I've updated LightCone8 at https://light-cone-calc.github.io/, everything should work including unit conversions and all. One final thing: LightCone7 uses a CMB temperature of 2.75 K, I have been using 2.7548 K in LightCone8, hope that is OK.

The recession rate chart is exactly the same as LightCone7.

 

[pbuk]

The version with $\Omega_{\Lambda,0}$ should also work in universes that are not close to our own.

That's good as that's what I've used and the UI allows entering the parameters to suit different universes!

 

[Jorrie]

think I'm all done now, I've updated LightCone8 at https://light-cone-calc.github.io/, everything should work including unit conversions and all. One final thing: LightCone7 uses a CMB temperature of 2.75 K, I have been using 2.7548 K in LightCone8, hope that is OK.

Great Job! Now we have the UI and the tables agreeing on $\Omega_m$ and also $\Omega$ staying at precisely 1. 00... I have also checked that for $\Omega \ne 1$ the behavior seems to be as expected (and interesting!).

$\Omega$ had to start very much closer to 1 and diverged to the requested 1.1 at present, but then the cosmological constant will drive it back towards 1 (for $w=-1$ at least)

 

[JimJCW]

Minor corrections (red) in Calculation.js required to improve the issue:

var Ol = (Ynow / Yinf) * (Ynow / Yinf); // Lambda density parameter
var Om = (Omega - Ol) * s_eq / (s_eq + 1); // matter density parameter (corrected for radiation)
var Or = Om / s_eq; // Radiation density parameter (corrected for error; was /(s_eq+1)
var Ok = 1 - Om - Or - Ol; // curvature density par

rholNow = rhocritNow * Ol;
rhol = rholNow; // vacuum energy density remains constant
rhomNow = (rhoNow - rhol)*s_eq/(s_eq+1); // Matter energy density corrected for radiation
rhom = rhomNow *s*s*s; // Matter density at time T

I think these changes are incorrect. The following statements near Line 27 and Line 325 in Calculation.js,

correctly represent the following equations (note that s_eq = z_eq in LightCone7),

Ω_m = (Ω – Ω_Λ) (z_eq + 1) / (z_eq + 2)​

Ω_R = Ω_m / (z_eq + 1) .​

They are consistent with equations in Tutorial, part III:

Ω_m = (Ω - Ω_Λ) S_eq / (1 + S_eq)​

Ω_R = Ω_m / S_eq ,​

where S_eq = z_eq + 1.

However, the following statements near Line 261,

should be changed to,

rhomNow = (rhoNow - rhol) * (s_eq + 1) / (s_eq + 2);​

rhom = rhomNow *s*s*s;​

rhorNow = rhomNow / (s_eq + 1);​

I believe this is the source of the error discussed in Post #109. @pbuk

 

[pbuk]

I think these changes are incorrect.
...
should be changed to,

rhomNow = (rhoNow - rhol) * (s_eq + 1) / (s_eq + 2);​

rhom = rhomNow *s*s*s;​

rhorNow = rhomNow / (s_eq + 1);​

I believe this is the source of the error discussed in Post #109. @pbuk

Yes you are probably right, however LightCone8 now calculates everything differently as discussed above. By going back to square one on the Physics, aided by @George Jones cross-check to the equations used, I have eliminated the intermediate calculations involving the various densities as below:

https://github.com/cosmic-expansion...bc785935a20b65cda6beb725897/src/model.ts#L201

      const eSquared = getESquaredAtStretch(s);
      const s2 = s * s;
      const h = h0 * Math.sqrt(eSquared);
      const omegaM = (omegaM0 * s2 * s) / eSquared;
      const omegaLambda = omegaLambda0 / eSquared;
      const omegaRad = (omegaRad0 * s2 * s2) / eSquared;
      return {
        h,
        omegaM,
        omegaLambda,
        omegaRad,
        omega: omegaM + omegaLambda + omegaRad,
        temperature: temperature0 * s,
        rhoCrit: rhoCrit0 * eSquared,
      };

 

[George Jones]

aided by @George Jones cross-check to the equations used

@Jorrie and @pbuk have done tremendous work on LightCone; I haven't done that much. All the stuff I have done in this thread has been performed at coffee shops while "on vacation" at my in-laws' place (11 loud people in the house), and while 3300km from my books. Flying home today.

Sent from my Tim Hortons at Major MacK and Jane.

 

[pbuk]

@Jorrie and @pbuk have done tremendous work on LightCone; I haven't done that much.

That's exactly why your contribution has been valuable - it is always important to have eyes on the detail from outside the core team.

 

[Jorrie]

I have updated the Lightcone Tutorial sections to reflect the present status. Will the other members of the core team please comment, specifically on what we should have in Part III, the calculation section.

 

[JimJCW]

Some observations:

(1) The suggested correction to LightCone7 discussed in Post #136,

where, (s_eq = z_eq), is consistent with the following statements in LightCone8:

(2) It is noticed that whenever z_lower in LightCone8 is set to be 0, there is extra line printed. For example,

This is a little bit strange, but inconsequential.

@Jorrie, @pbuk

 

[pbuk]

Yes it is annoying: it is a consequence of ensuring that ranges spanning z = 0 always include z = 0 exactly (otherwise lines which are not smooth there look wrong). I'll add it to the bottom of the list

 

[JimJCW]

Yes it is annoying: it is a consequence of ensuring that ranges spanning z = 0 always include z = 0 exactly (otherwise lines which are not smooth there look wrong). I'll add it to the bottom of the list

I am sorry if I was a little pushy. I misunderstood your writing in Post #133:

“I think I'm all done now, . . .”​

 

[Jorrie]

Yes it is annoying: it is a consequence of ensuring that ranges spanning z = 0 always include z = 0 exactly (otherwise lines which are not smooth there look wrong). I'll add it to the bottom of the list

I have looked at all the code, tested as far as I could and I am very happy with it. Your model.ts code is a little outside of my capabilities as a coder, so I cannot really assist with the doubling of the z=0 line issue.

I'm happy to look after the UI and the Tutorial, so if anyone has any suggestion, we can discuss and implement if approved by the core team. The equations in Tutorial part III still represent the LCDM model accurately. Do you think we must put more in there to indicate how you have implemented the model?

 

[JimJCW]

Hi @pbuk and @Jorrie,

I am getting some puzzling results concerning output of LightCone8. Let’s use PLANCK 2015 Data and calculate OmegaM and OmegaR with three different methods:

1. LightCone8
2. Correct equations,

3. Wrong equations,

omegaM0 = (1 - OmegaL0) * zEq / (zEq + 1)​

omegaR0 = OmegaM0 / zEq​

The results are shown in the table below:

Is it possible that the LightCone8 result is calculated with the wrong equations?

 

[pbuk]

Is it possible that the LightCone8 result is calculated with the wrong equations?

No, I have looked into it and it is more subtle than that. The problem is that LightCone8 uses the same mapping to the UI as LightCone7 which does not pass the input value of $\Omega_\Lambda$ directly to the caculations, it actually passes
$$ \left ( \frac{\frac{978}{H_0}}{\frac{978}{H_0\sqrt{\Omega_\Lambda}}} \right )^2 !$$
This obviously introduces roundoff error, and because of the sqaring with a 53 bit mantissa I'd expect only about 8 decimal digits of accuracy - your calculation shows ~7.5. The fact that a similar result occurs when taking $z_{eq}$ instead of $s_{eq}$ is coincidence.

I can easily fix this by changing the inputs to the back-end to read directly from the input fields.

 

[pbuk]

I can easily fix this by changing the inputs to the back-end to read directly from the input fields.

Well I did that, which was very worthwhile but it didn't fix the problem, which was actually less subtle: the input field that is labelled $z_{eq}$ was actually

<input name="s_eq" />

so while it looked as though the conversion $s_{eq} = z_{eq} + 1$ had been done, it never was!

I have now fixed this, and tidied up the credits and a couple of other things, enjoy.

 

[JimJCW]

Well I did that, which was very worthwhile but it didn't fix the problem, which was actually less subtle: the input field that is labelled $z_{eq}$ was actually

<input name="s_eq" />

so while it looked as though the conversion $s_{eq} = z_{eq} + 1$ had been done, it never was!

I have now fixed this, and tidied up the credits and a couple of other things, enjoy.

Great! Now we have,

	z​	OmegaM​	OmegaL​	OmegaR​	OmegaT​
LightCone8	0​	0.3089083630000​	0.691​	0.0000916370107​	1​
Wrong Eqs	0​	0.3089083358054​	0.691​	0.0000916641946​	1​
Correct Eqs	0​	0.3089083629893​	0.691​	0.0000916370107​	1​

 

[pbuk]

Thanks, of course the UI only displays a maximum 9 digits, the results at full precision are:

• OmegaMatterT: 0.3089083629893239
• OmegaLambdaT: 0.691
• OmegaRadiationT: 0.0000916370106761566
• OmegaTotalT: 1

A useful feature for the UI would be display and/or download in CSV and/or JSON format of the full precision results.

 

[Jorrie]

Thanks, of course the UI only displays a maximum 9 digits, the results at full precision are:

----

A useful feature for the UI would be display and/or download in CSV and/or JSON format of the full precision results.

Could be nice to have for model comparison, but we know the cosmological inputs to no better than the 4th digit, with $H_0$ the main constraint. A more important UI change would be to have the latest (2018/2021) Planck data release as an option.

Also, historically, LightCone has concentrated on the $H_0$ and $\Omega_\Lambda$ as base inputs. This was partially because the late @marcus had a lot of threads explaining the usefulness of having $\Omega_\Lambda$ a parameter.

However, the Planck base group of sampled data does not include $\Omega_\Lambda$, but rather the two matter densities: baryonic and dark - see Table 2 on page 16 in the paper (partially shown below) - which they then statistically add together to get $\Omega_m$. So $\Omega_\Lambda$ is a derived parameter in their books.

In the past data releases, we took the LightCone inputs from the last row of the second group of derived inputs in Table 2. I suggest that we keep on doing that, except for giving priority to $\Omega_m$ rather than $\Omega_\Lambda$, the latter becoming a derived parameter.

Before we change the mapping between the UI and the calculation module, I think we should decide on this.

BTW, since the "glitch" has now been fixed, I think we should start a new thread on LightCone Improvement and close this one.