CDM Phase-Space Statistic Kernel Calculator in the Zel'dovich Approximation

Author: Svetlin V. Tassev

This repository accompanies the paper Dominated-Convergence Failure in Cosmological Perturbation Theory and a Numerical Foundation for BBGKY+ZA. It contains the C++ pipeline that evaluates the Zel'dovich phase-space two-point kernel $\mathcal{P}(\boldsymbol k,\boldsymbol l,\boldsymbol m)$ of the paper's Section sec:corrZA, the Python helpers that prepare its cosmology input, and the Python figure scripts that produced the numerical-results panels of Section sec:Pnumresults and Appendix app:asympt.

The kernel $\mathcal{P}$ is the obstacle that has so far kept the BBGKY+ZA programme from being implemented at the cumulant level. This code provides a stable, arbitrary-precision implementation of the angular algebra and an FFTLog-based implementation of the radial Hankel transforms that together yield $\mathcal{P}$ to $\mathcal{O}(1\%)$ over the displayed $(k,p,\mu_{lm})$ range of the paper.

Repository Layout

Path	Role
makefile	Builds the nine maintained executables (seven FINALP_exp lane + two FINALP cross-check lane)
*.cpp, *.h, *.hpp	C++ sources of the angular, radial, merge and contraction stages
Compute_*.sh, Merge_*.sh, Run_*.sh	Shell wrappers for the maintained pipeline stages
Python/	Cosmology-input generation and HDF5 comparison helpers (Python)
analyze/CLM_check/	Python scripts that produce the figures of the paper from the HDF5 outputs
data/	Runtime input/cache directory, kept in git via .gitkeep
DOCUMENTATION.md	Detailed code, file-format, and code-to-paper reference
FIGURES.md	Provenance of every paper figure that this repository produces
LICENSE	GNU GPL v3 license text

The executables are intended to be run from the repository root so that the relative runtime directory data/ resolves correctly.

Paper-to-Code Variable Map

The code stores quantities under names that differ from those in the paper by historical convention. The mapping is:

Paper symbol (tex)	Code symbol / file	Notes
$\mathcal{P}(\boldsymbol k,\boldsymbol l,\boldsymbol m)$	mathcalP	the kernel of eq. (Pcurly), eqs. (FINALP_exp,FINALP)
$\mathcal{C}_{LM}(s,t,\mu_{lm})$	mathcalC_LM2D (HEALPix-sampled)	cosmology-independent angular polynomial, eq. (CLM_FINAL)
$\mathcal{C}_{LM,ab}(\mu_{lm})$	mathcalC_coeff_LM2D (coefficient table)	$(X,Z)$-substituted coefficients, eq. (CLM_ab_def)
$\mathcal{G}_{LM}(k,p,\mu_{lm})$	mathcalG (merged files)	radial Hankel transform piece, eq. (Glmn); stored scaled by $(2\pi/k)^{3/2}$
$\mathcal{H}_M(q,p,\mu_{lm})$	inner integrand of compute_mathcalG_noSumN.cpp	eq. (Hlmn)
$\zeta(q)$, $\gamma(q)$	zeta, gamma (in read_cosmo.cpp, data/zeta_gamma.txt)	cosmological inputs, derived from $P_L(k)$
$\sigma_v^2$	sigmav2	velocity-dispersion scalar
$P_Z(k,D)$	pk_zeldovich, pk_zeldovich_standard (in power_spectra.h5)	Zel'dovich density power spectrum, Appendix app:PZ
$P_L(k)$	pk (in cosmology_input.h5)	linear matter power spectrum
$X$ (paper, eq. Psub_def)	$P$ (code symbol, P_powers axis)	substituted angular variable
$Z$ (paper, eq. Psub_def)	$Q$ (code symbol, Q_powers axis)	substituted angular variable
$\mathcal{P}_*(p,k)$	reference slice $(\mu_{lm}{=}-1,X{=}0,Z{=}0)$ of mathcalP_PQ_coefficients	eq. (Pstar_def)

The translation P_code = X_paper, Q_code = Z_paper is centralized in init_pq_linear_substitution() inside C_LM_PQ.cpp. The HDF5 datasets P_powers[i] and Q_powers[j] carry the integer exponents stored along the third and fourth axes of mathcalC_coeff_LM2D and mathcalP_PQ_coefficients, so the reconstruction reads

$$\mathcal{P}(\boldsymbol k,\boldsymbol l,\boldsymbol m) \;=\;\sum_{i,j}\;\text{coefficient}[i,j]\;\;X^{\,P_{\text{powers}}[i]}\;Z^{\,Q_{\text{powers}}[j]} ,$$

with $X=P_{\text{code}}$ and $Z=Q_{\text{code}}$.

One subtlety to keep in mind. The code's substitution (eq. Psub_def) $s = (1-\mu_{lm})(1-P) + 2\mu_{lm}Q$, $t = \tfrac{1}{2}(\mu_{lm}-1)(1-P) + Q$ has both $s$ and $t$ carrying $P$ with a coefficient proportional to $(1-\mu_{lm})$ ($s_P=\mu_{lm}-1$, $t_P=\tfrac{1}{2}(1-\mu_{lm})$), so every power $P^a$ pulls a factor $(1-\mu_{lm})^a$ into the stored coefficient when $t^x s^y$ is expanded. Eq. (Pexp_XZ) of the paper writes the same kernel as $\mathcal{P}=\sum_{a,b} c_{ab}(\mu_{lm},p,k)\,(1-\mu_{lm})^a X^a Z^b$, keeping the $(1-\mu_{lm})^a$ factor outside $c_{ab}$. The stored coefficient is therefore the paper's $c_{ab}$ already multiplied by the $(1-\mu_{lm})^a$ prefactor at $a=P_{\text{powers}}[i]$, $b=Q_{\text{powers}}[j]$; to recover the paper-bare $c_{ab}$ at a non-collinear $\mu_{lm}$, divide by $(1-\mu_{lm})^a$.

The paper's truncation of the $(L,M)$ sums of eq. (FINALP_exp) is exposed at run time via the environment variables MATHCAL_MAX_L and MATHCAL_MAX_M, with maintained defaults 30/100. The numerical tiering of Appendix app:numerics — $(10,50)$, $(30,100)$, $(50,120)$ across three $p$ tiers — is realised by re-running the maintained pipeline at each tier and concatenating the results across $p$.

Pipeline

The kernel $\mathcal{P}$ is built in four stages, all keyed by a discrete index i_LM $\in[0,499]$ that maps to a logarithmic grid in $p\equiv lm$:

$$p(i_{LM})\;=\;10^{(i_{LM}/499)\,(\log_{10} 16 - \log_{10} 10^{-8}) + \log_{10} 10^{-8}}\,.$$

The four stages and their outputs are:

1. Angular stage. Either ComputeMathcalC_LM2D (paper's HEALPix-sampled $\mathcal{C}_{LM}(s,t,\mu_{lm})$ via eq. (SLM1)) or ComputeMathcalC_LM_PQ (paper's $(X,Z)$-substituted $\mathcal{C}_{LM,ab}(\mu_{lm})$ via eq. (CLM_ab_def)). Cosmology-independent. Produces one HDF5 file per mu_lm grid point on the linear grid $\mu_{lm}=\cos(\pi\,\mathrm{NUM}/100)$ with $\mathrm{NUM}\in[0,100]$.
2. Radial stage. ComputeMathcalG_noSumN evaluates the Hankel transform in eq. (Glmn) at bias = -0.9 and bias = +0.9 using FFTLog. Produces one raw HDF5 file per i_LM per bias.
3. Merge stage. MergeMathcalGs stitches the two biases onto a single $k$ grid using the matched-derivative criterion described in Appendix app:numerics. Produces one merged HDF5 file per i_LM.
4. Contraction stage. ComputeMathcalP_noSumN (HEALPix lane) or ComputeMathcalP_QP (coefficient lane) contracts the angular and merged radial files into $\mathcal{P}$. The figures of the paper are produced from the coefficient-lane output mathcalP_PQ_i.h5.

The figures of the paper come from the coefficient-lane (PQ) outputs mathcalP_PQ_<i_LM>.h5. The HEALPix-lane outputs P_<i_LM>.h5 are retained as a fully resummed reference and are cross-checked against the PQ reconstruction by Compare_healpix_PQ.sh (see also analyze/CLM_check/plot_mathcalP_healpix_vs_PQ.py).

ExportPowerSpectra is optional and writes a 1D diagnostic file power_spectra.h5 tied to Appendix app:PZ; it is not consumed by any of the four stages above.

A second, independent representation of $\mathcal{P}$ (eq. (FINALP) of the paper, the fully separated $\mathcal{E}_{LMN}$ form) is also implemented as a cross-check. The lane has two stages:

1. ComputeMathcalE (driven by Compute_MathcalE.sh): the radial stage that sums over the explicit-$N$ index inside the Hankel transform. Writes one mathcalE_LMN_scaled_by_2pi_over_k_to_3_halves_LM_<LM>_MaxL_10_MaxM_100_MaxN_100.h5 file per LM=l*m value from a hard-coded list of discrete LM points ({1e-10, ..., 3.5} in compute_mathcalE.cpp::main()).
2. ComputeMathcalP_E_LMN (driven by Compute_MathcalP_E_LMN.sh): a single-pixel HEALPix contraction that pairs one mathcalE_LMN_*.h5 file with one mathcalC_LM2D_<NUM>_<DEN>.h5 slice and prints the contracted $\mathcal{P}$ to stdout at a hard-coded test point $(l,m,I_{\rm pixel})=(4,4,N_{\rm pixels}-1)$.

This is the FINALP-vs-FINALP_exp cross-check that Section sec:Pnumerics of the paper refers to. The production figures use the FINALP_exp ($\mathcal{G}_{LM}$) lane because it converges to a fixed precision faster on the author's hardware. To run the cross-check at a different $(l,m,I_{\rm pixel})$ or a different LM list, edit the hard-coded values in compute_mathcalE.cpp::main() and compute_mathcalP_E_LMN.cpp::main().

Requirements

• g++ and GNU make (the source warns that clang was not reliable for this code path)
• HDF5 with C++ headers and libraries
• FFTW3 (with OpenMP support; fftw3l if USE_LONG_DOUBLE=1)
• FLINT / ARB / ACB (for the arbitrary-precision angular and FLINT-path radial branches)
• GSL
• OpenMP
• HEALPix C++ (only required by the HEALPix-sampled angular stage and the HEALPix-lane contraction)
• Python 3 with numpy, scipy, matplotlib, h5py, pyccl (for the cosmology generator and the figure scripts)

The figure-production Python scripts in analyze/CLM_check/ additionally use matplotlib's LaTeX backend (text.usetex = True) at run time; install a working texlive if you want to reproduce the rendered panels exactly.

Python Environment

A working, least-fragile setup is a dedicated conda environment from conda-forge:

conda create -n vlasov-2pt -c conda-forge python=3.12 pyccl numpy scipy matplotlib h5py
conda init zsh
exec zsh
conda activate vlasov-2pt

If activation is not available in the current shell, use:

eval "$(conda shell.zsh hook)"
conda activate vlasov-2pt

For one-off runs, skip activation and prefix commands with conda run -n vlasov-2pt.

Quick Start (coefficient lane, used by the paper figures)

From the repository root:

# 1. write data/cosmology_input.h5 from pyccl
conda run -n vlasov-2pt python3 Python/calculate_cosmology_input.py --no-show

# 2. build all nine maintained executables
make

mkdir -p output/mathcalC_coeff output/mathcalG output/mathcalP_PQ

# 3. angular stage (coefficient lane)
MATHCALC_COEFF_DIR=output/mathcalC_coeff ./Compute_MathcalC_LM_PQ.sh

# 4. radial stage (writes two raw biases per i_LM into MATHCALG_DIR)
RADIAL_OUTPUT_DIR=output/mathcalG ./Compute_MathcalRadial.sh

# 5. merge the two bias branches per i_LM in place
MATHCALG_DIR=output/mathcalG ./Merge_MathcalGs.sh

# 6. contract into mathcalP_PQ_<i_LM>.h5
MATHCALC_COEFF_DIR=output/mathcalC_coeff \
MATHCALG_DIR=output/mathcalG \
OUTPUT_DIR=output/mathcalP_PQ \
./Compute_MathcalP_QP.sh

The HEALPix lane uses the same data preparation, radial, and merge stages, and substitutes Compute_MathcalC_LM2D.sh and Compute_MathcalP.sh (writing P_<i_LM>.h5) for stages 3 and 6.

By default the four stages use the maintained truncation (MATHCAL_MAX_L, MATHCAL_MAX_M) = (30, 100). The paper's tiered production uses (10, 50), (30, 100), and (50, 120) across three $p$-tiers; set those values via environment variables before each tier to switch to the corresponding custom lane. Outputs at non-default cutoffs gain a MaxL_<L>_MaxM_<M> filename suffix so the maintained and custom lanes do not collide.

The shared-cache wrappers Run_mathcalP_PQ_shared_RESULTS.sh, Run_mathcalP_PQ_shared_RESULTS_range.sh, and Run_mathcalP_Healpix_shared_RESULTS.sh codify the production layout used to produce the figures of the paper. They default the shared runtime/output root to /run/media/user/SEAGATE/RESULTS; export OUTPUT_ROOT=<your path> before invocation, or edit the OUTPUT_ROOT="${OUTPUT_ROOT:-...}" line at the top of each script, to run them on a different host.

Reproducing the Paper Figures

The provenance of every figure that this repository produces is documented in FIGURES.md. In short, the figures of Section sec:Pnumresults are produced by Python scripts in analyze/CLM_check/ from mathcalP_PQ_<i_LM>.h5 files; the appendix asymptotic figures (Figs. fig:pz_asympt_late_D_low_k and fig:pz_asympt_late_D_high_k) and the angular reference figures (Figs. fig:L=1,fig:L=3,fig:L=10) are produced externally (Mathematica) and not from this repository's Python code.

Pipeline Notes

• All stage wrappers are resumable. Each one skips outputs that already exist; set the corresponding FORCE_RECOMPUTE_* environment variable (see DOCUMENTATION.md for the full list) to overwrite.
• Compute_MathcalC_LM2D.sh and Compute_MathcalC_LM_PQ.sh accept MATHCALC_PARALLEL_JOBS=<N> and use per-output lock directories so that two processes do not compute the same NUM/DEN file at once.
• Compute_MathcalRadial.sh accepts MATHCALG_PARALLEL_JOBS=<N> and launches separate bias runs in isolated worker directories.
• Merge_MathcalGs.sh preserves the raw bias-split inputs by default; set DELETE_RAW_MATHCALG=1 to delete them after a successful merge. If a merge attempt fails for a given i_LM the wrapper retries with thresholds enlarged by 1.5x up to three times before skipping that i_LM.
• Set TARGET_I_LM=<i> to make Compute_MathcalRadial.sh, Merge_MathcalGs.sh, Compute_MathcalP.sh, and Compute_MathcalP_QP.sh work on one i_LM value instead of their default ranges.
• Set FAST_MATHCALP=1 to keep the same convergence/error calculation while switching the mathcalP contraction to a faster plain double accumulator (the default uses Neumaier-compensated long double).
• Compute_MathcalP.sh assumes the full NUM = 0..100 angular grid exists.

Full details of all wrapper environment variables, HDF5 file schemas, the merge logic, and the convergence flag are in DOCUMENTATION.md.

Citation

If you use this code, please cite the accompanying paper.

License

This software is distributed under the terms of the GNU General Public License version 3 or later. See LICENSE for the full text.

---

Note: the documentation files in this repository (README.md, DOCUMENTATION.md, FIGURES.md) were generated automatically by inspection of the source code and the accompanying paper. They are intended to be code-first: where the documentation and the source disagree, the source is authoritative.