Velocisaurus: EFT for Peculiar Velocities

• Velocisaurus is a fast Python tool for computing effective-field-theory predictions of redshift-space peculiar velocity statistics and density-velocity correlations.
• It uses a generating-function formulation with modified IR resummation to accurately capture BAO features and extend standard EFT treatments.
• Validated with N-body simulations and Zeldovich predictions, the code achieves sub-percent precision in growth rate measurements for future surveys.

velocisaurus is a fast Python code released with “Effective Theories of Redshift-Space Galaxy Peculiar Velocities” for computing effective-field-theory predictions for peculiar-velocity statistics in large-scale structure. Its scientific role is to extend the EFT of large-scale structure, already widely used for redshift-space galaxy clustering, to the observables measured in peculiar-velocity surveys and closely related redshift-space distortion analyses. The code targets 2-point, density-weighted pairwise velocity statistics in redshift space, especially the mean pairwise velocity and pairwise velocity dispersion / second moment, together with the redshift-space galaxy power spectrum that generates them. The public release is associated with the repository https://github.com/sfschen/velocisaurus (Chen et al., 31 Jul 2025).

1. Scientific scope and target observables

The code is designed around the family of redshift-space spectra

$$P^s_{LL'}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ u_{s,n}^L(s_1) u_{s,n}^{L'}(s_2) \right\rangle_{s=s_1-s_2},$$

or, equivalently, in real-space variables,

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$

The especially useful combinations are the density-weighted pairwise moments

$$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$

with

$$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$

Accordingly, velocisaurus is aimed at the mean pairwise velocity and pairwise velocity dispersion / second moment, jointly with the redshift-space galaxy power spectrum (Chen et al., 31 Jul 2025).

A central conceptual feature is that these velocity statistics are modeled consistently with redshift-space galaxy densities using the same EFT parameter set, up to a finite extension in the generating-function language. This makes the code directly relevant for joint analyses of galaxy clustering and peculiar velocities. The same framework is explicitly connected to upcoming peculiar-velocity surveys and to future applications to the kinetic Sunyaev–Zeldovich effect, for which the paper notes that kSZ velocity statistics “have yet to be predicted within the EFT.” A plausible implication is that velocisaurus is intended not merely as a calculator for isolated velocity observables, but as a common inference layer for density–velocity analyses.

2. Generating-function formulation

The redshift-space power spectrum is treated as the generating function of the pairwise moments,

$$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$

with moments obtained by differentiation,

$$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$

Conversely,

$$P_{s}(k;\lambda) = \sum_{n=0}^\infty \frac{(ik\mu)^n \lambda^n }{n!} \tilde{\Xi}^{(n)}(k).$$

This relation is central to the code’s scientific role because it permits redshift-space density and velocity observables to be predicted within one EFT framework, with matched bias, counterterm, stochastic, and Fingers-of-God-like parameters (Chen et al., 31 Jul 2025).

This generating-function structure is not a purely formal convenience. It organizes the connection between galaxy clustering and pairwise-velocity observables and underlies the code’s claim that $P$, $v$, and $\sigma^2$ can be fit jointly with common EFT coefficients. In that sense, the code’s defining contribution is not only that it computes peculiar-velocity statistics, but that it does so in a way that preserves consistency with the density sector.

3. Lagrangian and Eulerian EFT realizations

The theoretical implementation is at next-to-leading, 1-loop order in both Lagrangian and Eulerian formulations of the EFT of large-scale structure (Chen et al., 31 Jul 2025). In Lagrangian perturbation theory, tracer motion is written in terms of displacements $$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$0, with redshift-space displacement

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$1

and biased tracers use the Lagrangian bias expansion

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$2

The redshift-space power spectrum is then

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$3

In Eulerian perturbation theory, the code uses the standard one-loop Eulerian bias basis,

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$4

and computes real-space moments $$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$5, then builds the redshift-space moments through the moment expansion. At 1-loop, real-space moments are nonzero up to $$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$6, so, for example,

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$7

The EFT ingredients are explicit. For redshift-space counterterms, the generating function receives

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$8

with only five independent counterterms needed in the generating function at 1-loop, and therefore the same five controlling all pairwise moments up to that order in real and redshift space. Stochastic terms are

$$P^s_{LL'}(k) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot\Delta u}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2))\ u_{n}^L(q_1) u_{n}^{L'}(q_2) \right\rangle_{q=q_2-q_1}.$$9

and the leading beyond-1-loop FoG-like pieces included are

$$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$0

The repeated emphasis is that the same EFT coefficients can be used to fit $$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$1, $$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$2, and $$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$3 jointly.

4. Infrared resummation and BAO response

What most distinguishes velocisaurus scientifically is its treatment of infrared resummation for baryon acoustic oscillations in velocity statistics (Chen et al., 31 Jul 2025). For galaxy densities, long-wavelength bulk flows smear the BAO mainly through the familiar exponential damping $$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$4. For peculiar velocities, however, the observables directly probe those bulk flows, so there are additional enhanced terms beyond pure damping. The BAO-displacement scale is

$$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$5

and in redshift space

$$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$6

The modified resummation for the wiggle part of the first pairwise moment is

$$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$7

not just a damped $$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$8. At 1-loop,

$$\tilde{\Xi}_s^{(n)}(k) = \int d^3 s\ e^{ik\cdot s} \left\langle (1 + \delta_{g,s}(s_1))(1+\delta_{g,s}(s_2))\ \Delta u_{n}^n \right\rangle_{s=s_1-s_2},$$9

Similarly, for the second moment,

$$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$0

The authors describe this as a recursive implementation in powers of the velocity in both Eulerian and Lagrangian approaches: the resummed power spectrum remains the generator, and derivatives with respect to $$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$1 generate the connected pieces order by order. In LPT they write

$$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$2

and then

$$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$3

This is the Lagrangian analogue of the recursive derivative scheme. A common misconception is that velocity BAO can be modeled by applying the usual exponential-only damping used for galaxy densities. The paper’s analysis indicates that this can fail badly, especially for $$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$4, because peculiar-velocity observables respond directly to long-wavelength bulk flows.

5. Implementation profile and likelihood use

The paper is more explicit about the theory backend than about software engineering (Chen et al., 31 Jul 2025). It does not give a detailed software manual, API, or benchmark table for velocisaurus; the explicit claims are that it is a fast Python code and that the authors “make our likelihood publicly available as part of velocisaurus.” For that reason, velocisaurus should be understood as a public Python implementation of the paper’s EFT predictions for peculiar-velocity observables and associated likelihood machinery, rather than as a general-purpose simulation or reconstruction package.

The implementation is connected to existing numerical methods. LPT calculations build on methods previously developed for velocileptors, especially FFTLog/Hankel-transform evaluation of the angular integrals. EPT predictions use the previously implemented one-loop real-space pairwise moments from velocileptors and extend them to redshift space with the new resummation. This suggests a division of labor in which velocisaurus specializes the peculiar-velocity and redshift-space extensions while reusing established perturbative infrastructure where appropriate.

Practical use is illustrated by the fitting analyses in the paper. The authors fit $$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$5, $$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$6, and $$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$7 jointly with common EFT parameters, use wide priors in some tests, vary $$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$8, and run Gaussian-likelihood analyses with Cobaya, analytically marginalizing linear parameters. The methodological significance is that the same generating-function-based EFT description underlies both prediction and inference.

6. Validation, fitting reach, and scientific relevance

Validation proceeds in two stages (Chen et al., 31 Jul 2025). First, the theory is checked in a Zeldovich universe, where the recursive IR-resummation scheme can be compared to exact Zeldovich predictions. In that setting, the modified resummation captures the BAO phase and amplitude in $$v_s(k) = \tilde{\Xi}^{(1)}_s = 2 P^s_{01}, \qquad \sigma^2_s(k) = \tilde{\Xi}^{(2)}_s = 2P^s_{02} - 2 P_{11}.$$9 and $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$0, whereas an exponential-only damping scheme can fail badly, especially for $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$1.

Second, the full EFT predictions are compared to fully nonlinear N-body simulations, specifically 25 AbacusSummit boxes of volume $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$2, using CompaSO halo catalogs for halos with $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$3 at $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$4 and $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$5. Agreement is shown for the first three redshift-space pairwise moments $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$6, $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$7, $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$8, including multipoles. The quoted fitting reach is roughly $$P_s(k ; \lambda) = \int d^3 q\ e^{ik\cdot q} \left\langle e^{ik\cdot(\lambda \Delta u)}\ (1 + \delta_{g}(q_1))(1+\delta_{g}(q_2)) \right\rangle,$$9 at $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$0 and $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$1 at $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$2, with some higher multipoles restricted to slightly smaller $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$3. The BAO wiggles isolated from the data are also reproduced well by both Eulerian and Lagrangian predictions.

In fits where all cosmological parameters except the growth rate are fixed, and $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$4 is allowed to vary jointly with the EFT nuisance parameters, the model recovers the true growth rate to better than percent level: the reported constraints are $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$5 and $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$6 on $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$7 at $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$8 and $$\tilde{\Xi}_s^{(n)}(k) = \frac{1}{(ik\mu)^n}\left( \frac{d^n P_s(k;\lambda)}{d \lambda^n} \right)_{\lambda=1}.$$9, with best-fit values deviating by at most $$P_{s}(k;\lambda) = \sum_{n=0}^\infty \frac{(ik\mu)^n \lambda^n }{n!} \tilde{\Xi}^{(n)}(k).$$0 from the truth. This matters because the code is intended for upcoming peculiar-velocity surveys whose statistical precision will demand sub-percent to percent-level modeling accuracy; the authors explicitly state that this is “well beyond the statistical precision required for upcoming peculiar velocity surveys.”

Taken together, these features explain the code’s significance. Standard EFT galaxy-clustering codes focus on density and redshift-space distortion observables, whereas velocisaurus targets peculiar velocity statistics themselves in a framework where densities and velocities are derived consistently from the same redshift-space generating function and the same EFT parameter set. Its distinctive methodological contribution is the demonstrated need for modified BAO IR resummation in peculiar-velocity statistics, rather than merely the standard exponential damping used for galaxy densities. That makes it directly relevant for low-redshift peculiar-velocity surveys, joint density–velocity analyses, RSD/BAO inference, and likely future EFT treatments of kSZ observables.