Control variates from Eulerian and Lagrangian perturbation theory: Application to the bispectrum

TL;DR

The paper tackles variance reduction in matter-bispectrum measurements from N-body simulations by constructing perturbative control variates from Eulerian and Lagrangian theory. It reveals that pure Eulerian control variates suffer exponential decorrelation due to large-scale displacements, while Lagrangian approaches (ZA and 2LPT) preserve correlations but have analytic limitations, motivating the shifted control variate that combines tractability with strong correlation. The shifted control variate, along with, in some regimes, 2LPT, achieves substantial variance reductions—up to factors of around $10^4$ for some configurations—and enables sub-2% precision for the matter bispectrum in a single $1\,\mathrm{Gpc}/h$^3 box, paving the way for accurate bispectrum emulators and extensions to higher-order statistics. The work suggests broad future directions, including applying shifted control variates to biased tracers and redshift-space bispectra, and extending the framework to the full $n$-point hierarchy.

Abstract

Control variates have seen recent interest as a powerful technique to reduce the variance of summary statistics measured from costly cosmological $N$-body simulations. Of particular interest are the class of control variates which are analytically calculable, such as the recently introduced 'Zeldovich control variates' for the power spectrum of matter and biased tracers. In this work we present the construction of perturbative control variates in Eulerian and Lagrangian perturbation theory, and adopt the matter bispectrum as a case study. Eulerian control variates are analytically tractable for all $n$-point functions, but we show that their correlation with the $N$-body $n$-point function decays at a rate proportional to the sum-of-squared wavenumbers, hampering their utility. We show that the Zeldovich approximation, while possessing an analytically calculable bispectrum, is less correlated at low-$k$ than its Eulerian counterpart. We introduce an alternative -- the 'shifted control variate' -- which can be constructed to have the correct tree-level $n$-point function, is Zeldovich-resummed, and in principle has an analytically tractable bispectrum. We find that applying this shifted control variate to the $z=0.5$ matter bispectrum is equivalent to averaging over $10^4$ simulations for the lowest-$k$ triangles considered. With a single $V=1({\rm Gpc}/h)^3$ $N$-body simulation, for a binning scheme with $N\approx 1400$ triangles from $k_{\rm min} = 0.04 h {\rm Mpc}^{-1}$ to $k_{\rm \max} = 0.47 h {\rm Mpc}^{-1}$, we obtain sub-2% precision for every triangle configuration measured. This work enables the development of accurate bispectrum emulators -- a probe of cosmology well-suited to simulation-based modeling -- and lays the theoretical groundwork to extend control variates for the entire $n$-point hierarchy.

Figures (15)

• Figure 1: Visualization of the Eulerian and Lagrangian perturbative control variate fields considered in this work, as a function of smoothing scale. The left panels show the $N$-body density field centered around the largest overdensity of the simulation at $z=0.5$, projected across 20 Mpc/$h$. When filtered on a smoothing scale with $R_{\rm smooth} > \Sigma \sim 4.5 {\rm Mpc}/h$ (at which EPT decorrelates) the $N$-body, Eulerian and Lagrangian fields are all visually correlated. However, smoothing at a smaller scale reveals the rapid decorrelation of EPT, while the Lagrangian fields maintain a large degree of similarity to the $N$-body distribution. The color scale is chosen to be symmetric around $\pm 3 \sigma_{\rm lin}(R)$, where $\sigma_{\rm lin}(R)$ is the standard deviation of the Gaussian panel at that smoothing scale.
• Figure 2: Top: Cross-correlation coefficient between bispectra in Zeldovich boxes and those measured from the same Gaussian initial conditions. The shape of this cross-correlation coefficient is in close agreement with the prediction from \ref{['eqn:rhozag']}. Bottom: Ratio of the empirically measured cross-correlation coefficients and the predictions from \ref{['eqn:rhozagapprox']}, in blue and \ref{['eqn:rhozag']}, in orange. The grey bands denote $\pm 5\%$ residuals. The left and right panels show the same data, but as a function of $k_{\rm sph} = \sqrt{(k_1^2 + k_2^2 + k_3^2)/3}$ on left and the triangle index representation on the right. The dominant behavior comes from the derived $k_{\rm sph}$ dependence. The gray band in the residuals denotes $\pm 5\%$ deviations from unity.
• Figure 3: Bispectrum cross-correlation coefficient between the Zeldovich density field and Eulerian realizations, including the 'Zeldovich tree level' bispectrum of \ref{['eqn:bktreeza']}. The cross-correlation coefficient is enhanced by a term $\propto k_{\rm sph}^2,$ as argued analytically. We also see a significant tightening of the cross-correlation coefficients for different triangle configurations
• Figure 4: Top panel: Relative deviation of tree level bispectrum averaged over $N=1000$ lattice-based realizations of the $\delta^{(2)}$ field compared to analytic predictions. The blue points show the tree level bispectrum when evaluated at the central $(k_1, k_2, k_3)$ of the bin, the orange curve shows the prediction when using the 'effective triangle' in \ref{['eqn:keff']} and the green points show the result from averaging over all discrete triangles. The error induced due to not bin-averaging can reach 30% for triangles in this scheme where $\Delta k = 2k_f$. Bottom panel: The same panel, but zoomed to show the size of spread when the discrete bin average is made. The gray bands indicate 1% scatter around zero. Error bars correspond to error bars on the mean of $N=1000$ realizations.
• Figure 5: Left: Ratio of bispectra in Eulerian perturbation theory, to order $\mathcal{O}(\epsilon^2)$ divided by the $N$-body bispectrum, as a function of $k_{\rm sph}^2 = (k_1^2 + k_2^2 + k_3^2)/3$. Each bispectrum is averaged over $N=1000$ realizations. The blue points indicate the Gaussian bispectrum, orange points include the tree-level contribution, and the green points include terms up to the bispectrum diagrams $B^{221}$ and $B^{311}$, which have zero mean. Right: Cross-correlation coefficients between the full $N$-body bispectrum and different realizations of the bispectrum measured in lattice Eulerian perturbation theory. Despite their zero-mean, the Gaussian term, as well as the fifth-order-in-$\delta$ terms, both contribute to increase the cross-correlation coefficient with the $N$-body bispectrum, and thus the total variance reduction achievable with EPT. The dashed vertical line in the right panel corresponds to the $k_{\rm sph}$ mode at which one e-folding of decorrelation is expected, according to \ref{['eqn:rhozagapprox']}.