Update velocities to 2LPT

[fjaviersanchez]

Currently the code uses linear velocities. @andreufont suggested that for Ly-alpha, taking into account 2LPT contributions to velocities might be important (not so much to the QSOs themselves as per @damonge). This can be a cool addition for the future.

[cramirezpe]

We are trying to include 2LPT velocities in CoLoRe right now. But we are unsure about how to implement it.

The relevant equations here are (from https://arxiv.org/abs/astro-ph/0112551): This means that the same potentials used for computing the displacement are the same that we should use to modify velocities. However, we are unsure on how to proceed, since LPT displacement fields and velocities are computed in different places of the code:

• The displacement fields computed through Lagrangian Perturbation Theory are computed inside the function compute_physical_density_field (called in main.c:69, defined in density.c:1160). The potentials are computed in:
  • For 1LPT in density.c:436.
  • For 2LPT in density.c:751.
• Velocities for catalog objects are computed inside the function srcs_set_cartesian (called in main.c:84, defined in srcs.c:285). Then there are multiple calls (srcs_set_cartesian_single -> get_rvel). This last function, defined in srcs.c:120 is the one that uses potentials to compute displacements.
• Velocities for skewers are computed inside the function get_beam_properties (called in main.c:125, defined in beaming.c:293). Then there are multiple calls (srcs_get_beam_properties_single -> interpolate_from_grid -> get_element). This last function is defined in beaming.c:31, here at lines 69, velocities are computed; in a similar way as the previous one.

So, putting it simple:

• LPT potentials are computed in main.c:69.
• Velocities for objects are computed in main.c:84.
• Velocities for skewers are computed in mainc.c:125. I think the best option is to store velocity in the r direction in a new array, then use it (if defined) in the next two steps. This means in a 4096 box: 4096^3*4bytes/1024^3=256Gb extra on memory; does it sound reasonable?

With the previous approach I am not including velocities not in the radial direction, there is a recent pull request including transverse velocities. If we want to also compute this, memory requirements will be x3(?).

I can start implementing it if this sounds ok, and will only need help reviewing it before merging.

[andreufont]

@damonge , @fjaviersanchez - any thought about the plan above? Should we have a call at some point? We could also invite Marc Manera and Santi Avila, who will join the project at IFAE.

[fjaviersanchez]

I think the plan sounds good. Just a couple of things that I think can be useful and may help with the memory usage: 1) You can set an option to allow the user run with the Gaussian velocities, so that way, even if you select LPT or 2LPT you don't need additional memory. 2) You don't need to stick to the original grid to store the radial velocity in memory, you can probably do it in a different grid size (I guess that it may slow down the code though). A call sounds good to me :).

[cramirezpe]

Thank you Javi for the comments. We could arrange a call for the first week of October (Santi will already be at IFAE by then). I will send an email so we can find a time that works for everyone.