Question about quantities rho_m and rho_r

[akshay-ghalsasi]

Hello,

I am a little confused about the definitions of rho_m and rho_r in background.c and perturbations.c

I understand that for components behaving like matter and radiation, they need to be included, however for the scalar field with a given potential, can have behavior distinct from radiation or matter. Yet the contribution from pressure of scalar field is included in rho_r and an equivalent amount is subtracted fro rho_m. See attached screenshot from background.c

However this change to rho_p and rho_m does not exist for the fluid module. See attached screenshot below.

Why is there a contribution to rho_r and rho_m in one case but not the other? Should the fluid module have that correction to rho_r and rho_m as well? Any help is appreciated. Thanks in advance!

[schoeneberg]

Hey @akshay-ghalsasi ,

I think the rho_m and rho_r contributions are not super important for the rest of the code, and are mostly used to compute "useful" quantities, like the matter-radiation equality time. I think you are right that adding the SCF contributions to rho_m and rho_r is very dangerous in general. Indeed, for many potentials this leads to non-sensical rho_m(z) and rho_r(z), as I have experienced in some of my previous works. I would love to ask the person who originally implemented it. As for your question, it's better not to add the fld to rho_m/rho_r, as you get a very strange splitting for something that is not a lot like radiation OR matter. E.g. if you use the typical split as 3p and rho-3p for dark energy, you get -3rho and 4 rho as the respective contributions, which can make the total radiation energy negative, and completely overestimates the matter density. As such, in principle you ONLY want to use this splitting for things close to radiation or matter. I think the SCF part could be removed in principle.

[ThomasTram]

Hi both

As @schoeneberg says, the quantities rho_r/Omega_r and rho_m/ Omega_m in background_functions() are not so important for the rest of the code. The most important cases are:

1. Checking that a=1e-14 is early enough for the Universe to be radiation dominated which is an assumption made when setting initial conditions for the background.
2. Computing the number of effective relativistic d.o.f.s (Neff) which is used for interpolating the BBN prediction for the Helium abundance, YHe.

If the scalar field is important at early times, i.e. rho_scf/rho_tot is non-negligible, keeping this contribution is important. For instance, the scalar-field might behave as matter at early times and dominate the Universe for at while, meaning that the initial condition for the background ODE equations would be wrong. Similarly, if the scalar field is important around BBN, its contribution to Neff should be included so that the assumptions about the expansion-rate during BBN are approximately correct. (Incidentally, Neff is computed at the initial time and not the BBN time for this purpose which should be fixed at some point.)

In both case 1 and 2, if the densities are oscillating, it should really be the oscillation-averaged values we want, and not just whatever the value is at the given time. Although this is not ideal, removing them would not be a solution. (And note that the current implementation already works as intended for non-oscillatory models.)

Concerning the fluid implementation, the rho_r and rho_m contributions should in fact be added if the fluid contribution is important at early times.

For perturbations, there is an additional quantity, delta_m which is used to compute the power spectrum. A scalar field is pretty general and in some cases it could be Dark Matter. In that case we need to include the scalar field in the matter power spectrum. On the other hand, if the scalar field is like quintessence, we probably should not include it in the matter power spectrum.

A part from delta_m, rho_m and rho_r are used in multiple locations for setting initial conditions in the perturbations module. The basic idea is exactly the same as in the background: we need the contributions from the scalar field, otherwise the initial conditions will not be accurate. To be concrete, we need to set initial conditions for theta_g, the velocity-divergence of the photon distribution. It depends on frac_nu which is the fraction of neutrinos, but more precisely it is the radiation energy density in free-streaming species divided by the radiation energy density in "non-free-streaming" species at that time. (At this time and scale, photons are tightly coupled which has been checked earlier.)

Now, the fluid is not included here, since an earlier check already checks that it is sub-dominant at early times. The scalar field is also not included, and I am afraid that we forgot to check this case. Thus, one must be very careful if one wants to use a scalar field or fluid that is relevant at early times.

Cheers, Thomas

[akshay-ghalsasi]

Hi @schoeneberg @ThomasTram . Thanks very much for the reply. It clarified quite a bit for me.

However I do have a question regarding the issue @schoeneberg pointed out. I have a fluid that behaves like matter (w = 0) in the early universe and kination (w = 1) later. It crosses w = 1/3 in between. It can (may or may not) dominate the energy budget of the universe before recombination or matter radiation equality. In such a case if the fluid dominates and my w > 1/3 and I do rho_m = rho_m - 3 p_fld I will end up with negative rho_m. Surely that is not physical?

I can put an if condition i.e. if w_kin < 1/3: rho_m = rho_m - 3 p_fld

but that seems rather abrupt and wont necessarily solve my problem depending on how much the fluid dominates over rho_m.

Is there a physical reason for doing this substitution of rho_r = rho_r + 3 p_fld rho_m = rho_m - 3 p_fld

For instance if I try to to derive initial conditions for my perturbations following Ma and Bertschinger, will I get the same result for my initial conditions as doing the substitutions above?

Once again thanks for the taking the time to reply.

Akshay

[schoeneberg]

Dear @akshay-ghalsasi : Indeed, what you describe seems like a case where sadly things are not as trivial to do. The rho_m ~= rho - 3p, rho_r ~= 3p setting is only an approximation that can be used for simple cases that are very much alike radiation or matter (e.g. massive neutrinos) in order to avoid having to derive the full complicated initial conditions (in background and perturbations) and the impact on BBN. Thanks to @ThomasTram for pointing out the other locations where these quantities are used. For an arbitrary fluid that dominates early dynamics, the assumptions made in deriving the initial conditions in the background and perturbations module fail, and you might have a strong impact on BBN, potentially causing different element abundances. As such, if you really have a model with w=0 in the early universe and w=1 later, you will have to consider when exactly the transition occurs. If it occurs for z>>10^9, then you will have to look at the impact of BBN and re-derive even the background initial conditions. If it occurs for z << 10^9, then probably BBN and the ICs of the background are unaffected, but you might have to re-derive the perturbation ICs which are set as late as slightly before the MR-equality (z~3000). If you are anyway at lower redshifts, then there's not much to worry about, and you can safely ignore the impact. In the other cases, it seems unavoidable to me to re-derive at least some of the ICs. (i.e. the substitutions that you wrote above would not be equivalent to deriving perturbations following Ma&Bertschinger).

Hope that helps!

[akshay-ghalsasi]

Thanks @schoeneberg ,

It briefly dominates for z << 10^9 but above 3000. I will just rederive the initial conditions and implement

Thanks, Akshay