One-Loop Nonlinear Matter Power Spectrum from Unified Lagrangian Perturbation Theory: Fast Computation and Comparison with Emulators

TL;DR

This paper introduces a fast, accurate, and IR-safe one-loop nonlinear matter power spectrum computation method using Unified Lagrangian Perturbation Theory, validated against emulators and suitable for large-scale structure analysis.

Contribution

The authors develop a novel ULPT-based framework that improves convergence and computational efficiency for nonlinear power spectrum modeling, with open-source implementation.

Findings

ULPT achieves 2-3% accuracy up to k ≈ 0.4 h/Mpc for z ≥ 0.5

The method is computationally efficient, taking about 2 seconds per evaluation

ULPT accurately reproduces BAO features and is validated against emulators.

Abstract

We present a fast and accurate formulation for computing the nonlinear matter power spectrum at one-loop order based on Unified Lagrangian Perturbation Theory (ULPT). ULPT decomposes the density field into the Jacobian deviation, capturing intrinsic nonlinear growth, and the displacement-mapping factor, accounting for large-scale distortions due to bulk flows. This structural separation leads to a natural division of the power spectrum into a source term and a displacement-mapping factor, ensuring infrared (IR) safety by construction. We implement an efficient numerical algorithm using FFTLog and FAST-PT, achieving approximately 2-second evaluations on a standard laptop. The results are validated against simulation-based emulators, including the Dark Emulator and Euclid Emulator 2. Across 100 sampled cosmologies, ULPT agrees with emulator predictions at the 2--3\% level up to \( k \simeq 0.4\,h\,\mathrm{Mpc}^{-1} \) for \( z \geq 0.5 \), without any nuisance parameters. Similar agreement is found in configuration space, where the two-point correlation function remains accurate down to \( r \simeq 10\,h^{-1}\mathrm{Mpc} \). Compared to standard perturbation theory, which fails at small scales due to series expansion of the displacement factor, ULPT maintains convergence by preserving its full exponential form. We also clarify the mechanism of BAO damping: exponential suppression by displacement and peak sharpening by nonlinear growth. The combination accurately reproduces BAO features seen in simulations. ULPT thus offers a robust, IR-safe, and computationally efficient framework for modeling large-scale structure in galaxy surveys. The numerical implementation developed in this work is publicly released as the open-source Python package \texttt{ulptkit} (https://github.com/naonori/ulptkit).

Tables

Table 1. Table 1: Coefficients c α ​ β ​ ℓ c_{\alpha\beta\ell} and c J , α ​ β ​ ℓ c_{{\rm J},\alpha\beta\ell} used in the Legendre expansion of mode-coupling kernels relevant to the 22-type term in the one-loop power spectrum. The coefficients c α ​ β ​ ℓ c_{\alpha\beta\ell} appear in the standard SPT expression for the one-loop power spectrum, while c J , α ​ β ​ ℓ c_{{\rm J},\alpha\beta\ell} correspond to the contributions retained in the ULPT source term. In particular, the IR-sensitive component ( α , β , ℓ ) = ( 2 , − 2 , 0 ) (\alpha,\beta,\ell)=(2,-2,0) is excluded in the ULPT formulation, leading to automatic IR regularization.

$\alpha$	$\beta$	$\ell$	$c_{\alpha \beta \ell }$	$c_{J,\alpha \beta \ell }$
0	0	0	$\frac{1219}{1470}$	$\frac{242}{735}$
0	0	2	$\frac{671}{1029}$	$\frac{671}{1029}$
0	0	4	$\frac{32}{1715}$	$\frac{32}{1715}$
1	-1	1	$\frac{62}{35}$	$\frac{27}{35}$
1	-1	3	$\frac{8}{35}$	$\frac{8}{35}$
2	-2	0	$\frac{1}{6}$	$0$
2	-2	2	$\frac{1}{3}$	$0$

Table 2. Table 2: Parameter ranges for Dark Emulator and Euclid Emulator 2

Parameter	Dark Emulator	Euclid Emulator 2
${\omega }_b={\Omega }_b{h}^2$	$[0.02114,0.02336]$	$[0.02133,0.02373]$
${\omega }_c={\Omega }_c{h}^2$	$[0.10782,0.13178]$	$[0.11293,0.15537]$
${\Omega }_{\mathrm{de}}$ or ${\Omega }_m$	$[0.54752,0.82128]$ (DE)	$[0.23116,0.38666]$ (${\Omega }_m$)
$\ln \left({10}^{10}{A}_s\right)$	$[2.4752,3.7128]$	$[2.398,3.178]$
$n_s$	$[0.9163,1.0127]$	$[0.928,1.000]$
$w_0$	$[-1.2,-0.8]$	$[-1.3,-0.5]$
$w_a$	—	$[-1.73,0.5]$
$\sum m_{\nu }$ [eV]	fixed (0.06 eV)	$[0.0,0.6]$
$h$	—	$[0.645,0.795]$

Equations

$$\rho (\mathbf{x})=\bar{\rho }\left[1+\delta (\mathbf{x})\right].$$

$$\mathbf{x}=\mathbf{q}+\mathbf{\Psi }(\mathbf{q}),$$

$$d^{3}x=J(\mathbf{q})d^{3}q,J(\mathbf{q})\equiv \det \left(\frac{\partial \mathbf{x}}{\partial \mathbf{q}}\right).$$

$$\rho (\mathbf{x})d^{3}x=\bar{\rho }{d}^{3}q.$$

$$\delta (\mathbf{q}+\mathbf{\Psi }(\mathbf{q}))=\frac{1}{J(\mathbf{q})}-1.$$

$$\delta (\mathbf{x})=\int d^3{x}^{\prime }\delta ({\mathbf{x}}^{\prime }){\delta }_{\mathrm{D}}(\mathbf{x}-{\mathbf{x}}^{\prime }),$$

$$\delta (\mathbf{x})=\int d^{3}q{\delta }_{\mathrm{J}}(\mathbf{q}){\delta }_{\mathrm{D}}(\mathbf{x}-\mathbf{q}-\mathbf{\Psi }(\mathbf{q})),$$

$${\delta }_{\mathrm{J}}(\mathbf{q})$$

$$=1-J(\mathbf{q}).$$

$$\tilde{\delta }(\mathbf{k})=\int d^{3}q{e}^{-i\mathbf{k}\cdot \mathbf{q}}{e}^{-i\mathbf{k}\cdot \mathbf{\Psi }(\mathbf{q})}{\delta }_{\mathrm{J}}(\mathbf{q}),$$

$$P(\mathbf{k})=\int d^{3}r{e}^{-i\mathbf{k}\cdot \mathbf{r}}⟨e^{-i\mathbf{k}\cdot \left[\mathbf{\Psi }(\mathbf{q})-\mathbf{\Psi }({\mathbf{q}}^{\prime })\right]}{\delta }_{\mathrm{J}}(\mathbf{q}){\delta }_{\mathrm{J}}({\mathbf{q}}^{\prime })⟩,$$

$$X\equiv -i\mathbf{k}\cdot \left[\mathbf{\Psi }(\mathbf{q})-\mathbf{\Psi }({\mathbf{q}}^{\prime })\right],$$

$$P(\mathbf{k})$$

$$\times \left[{⟨e^X{\delta }_{\mathrm{J}}(\mathbf{q}){\delta }_{\mathrm{J}}({\mathbf{q}}^{\prime })⟩}_{\mathrm{c}}+{⟨e^X{\delta }_{\mathrm{J}}(\mathbf{q})⟩}_{\mathrm{c}}{⟨e^X{\delta }_{\mathrm{J}}({\mathbf{q}}^{\prime })⟩}_{\mathrm{c}}\right],$$

$$⟨e^X⟩=\exp \left[-\overline{\Sigma }(\mathbf{k})+\Sigma (\mathbf{k},\mathbf{r})\right],$$

$$-\overline{\Sigma }(\mathbf{k})+\Sigma (\mathbf{k},\mathbf{r})=\sum _{m=2}^{\infty }\frac{1}{m!}{⟨{\left[-i\mathbf{k}\cdot (\mathbf{\Psi }(\mathbf{q})-\mathbf{\Psi }({\mathbf{q}}^{\prime }))\right]}^m⟩}_{\mathrm{c}}.$$

$$\overline{\Sigma }(\mathbf{k})\equiv \Sigma (\mathbf{k},\mathbf{r}=0).$$

$${\xi }_{\mathrm{J}}(\mathbf{r})\equiv {⟨e^X{\delta }_{\mathrm{J}}(\mathbf{q}){\delta }_{\mathrm{J}}({\mathbf{q}}^{\prime })⟩}_{\mathrm{c}}+{⟨e^X{\delta }_{\mathrm{J}}(\mathbf{q})⟩}_{\mathrm{c}}{⟨e^X{\delta }_{\mathrm{J}}({\mathbf{q}}^{\prime })⟩}_{\mathrm{c}}.$$

$${\xi }_{\mathrm{J}}(\mathbf{r})$$

$$+⟨(e^X-1){\delta }_{\mathrm{J}}(\mathbf{q}){⟩}_{\mathrm{c}}⟨(e^X-1){\delta }_{\mathrm{J}}({\mathbf{q}}^{\prime }){⟩}_{\mathrm{c}},$$

$$P(\mathbf{k})=e^{-\overline{\Sigma }(\mathbf{k})}\int d^{3}r{e}^{-i\mathbf{k}\cdot \mathbf{r}}{e}^{\Sigma (\mathbf{k},\mathbf{r})}{\xi }_{\mathrm{J}}(\mathbf{r}).$$

$$P_{\mathrm{J}}(\mathbf{k})\equiv \int d^{3}r{e}^{-i\mathbf{k}\cdot \mathbf{r}}{\xi }_{\mathrm{J}}(\mathbf{r}).$$

$$P(\mathbf{k})=P_{\mathrm{J}}(\mathbf{k})+P_{\mathrm{DM}}(\mathbf{k}),$$

$$P_{\mathrm{DM}}(\mathbf{k})=\int d^{3}r{e}^{-i\mathbf{k}\cdot \mathbf{r}}\left[e^{-\overline{\Sigma }(\mathbf{k})+\Sigma (\mathbf{k},\mathbf{r})}-1\right]{\xi }_{\mathrm{J}}(\mathbf{r}).$$

$$P(\mathbf{k})=P_{\mathrm{CF}}(\mathbf{k})+P_{\mathrm{CC}}(\mathbf{k}),$$

$$P_{\mathrm{CF}}(\mathbf{k})$$

$$P_{\mathrm{CC}}(\mathbf{k})$$

$${\tilde{f}}_{\ell }(k)=(-i{)}^{\ell }(4\pi ){\int }_0^{\infty }dr{r}^2{j}_{\ell }(kr)f(r),$$

$$f(r)=i^{\ell }\int \frac{dk{k}^2}{2{\pi }^2}{j}_{\ell }(kr){\tilde{f}}_{\ell }(k).$$

$$\Sigma (\mathbf{k},\mathbf{r})$$

$${\xi }_{\mathrm{J}}(r)$$

$$P_{\mathrm{J}}(k)=P_{\mathrm{lin}}(k)+P_{\mathrm{J}}^{(1\text{-loop})}(k),$$

$${\Sigma }^{(\mathrm{lin})}(\mathbf{k},\mathbf{r})=\int \frac{d^3{k}^{\prime }}{(2\pi {)}^3}{e}^{i{\mathbf{k}}^{\prime }\cdot \mathbf{r}}{\left(\frac{\mathbf{k}\cdot {\mathbf{k}}^{\prime }}{k^{\prime 2}}\right)}^2{P}_{\mathrm{lin}}(k^{\prime }).$$

Full text

I Introduction

The large-scale structure (LSS) of the Universe provides a powerful observational window into the physics of cosmic expansion and structure formation. Among various statistical observables, the matter power spectrum plays a central role by connecting the primordial fluctuations of the early Universe to the present-day distribution of matter.

To model the nonlinear evolution of the matter power spectrum, standard perturbation theory (SPT) has been widely employed (for a review, see Ref. [1]). However, SPT suffers from poor convergence at small scales due to long-wavelength displacements and nonlinear mode-coupling effects.

Lagrangian perturbation theory (LPT) offers a physically motivated alternative by treating the displacement vector as the primary dynamical variable. This perspective allows partial resummation of infrared (IR) contributions and has led to successful extensions such as the Zel’dovich approximation (ZA) [2], convolution LPT (CLPT) [3] (see also [4]), and IR-resummed approaches [5, 6, 7, 8, 9, 10, 11, 12], all of which have demonstrated accurate predictions for the nonlinear damping of the baryon acoustic oscillation (BAO) feature [13, 14].

Despite recent developments, many existing LPT-based frameworks still face several limitations. While the conventional LPT provides a unified treatment of density fluctuations in real and redshift space, it does not offer a consistent formulation that incorporates post-reconstruction space on an equal footing. In particular, it has been pointed out that IR safety is violated in reconstructed fields [15].

Moreover, when incorporating galaxy bias in LPT [16], one must account not only for the displacement field but also for biased density fluctuations defined in Lagrangian space. This additional ingredient partially breaks the structural unity of the LPT formulation.

Spectroscopic galaxy surveys measure the large-scale structure of the Universe in redshift space, where the observed galaxy density field is inevitably modulated by line-of-sight peculiar velocities. These distortions, known as redshift-space distortions (RSD) [17], introduce anisotropic features that must be accurately modeled to extract unbiased cosmological information from clustering measurements.

In parallel, density-field reconstruction techniques [18], originally developed to sharpen the BAO feature by undoing large-scale bulk flows, have recently been shown to improve constraints on cosmological parameters derived from the broadband shape of the power spectrum [19, 20], beyond the BAO feature. Together, these observational advances underscore the need for a unified theoretical framework that consistently describes density fluctuations across real space, redshift space, and post-reconstruction space.

To address these challenges, we have recently developed Unified Lagrangian Perturbation Theory (ULPT) [21], a new perturbative framework that enables a consistent and unified description of the density field in all three coordinate systems: real space, redshift space, and post-reconstruction space. Importantly, the ULPT formulation maintains IR safety across all three representations, ensuring that large-scale displacements are consistently resummed.

Unlike conventional approaches such as ZA, ULPT reformulates LPT by explicitly separating the density fluctuation into two physically distinct components: the Jacobian deviation, which encapsulates intrinsic linear and nonlinear growth (including bias effects), and the displacement-mapping effect, which accounts for large-scale convective distortions due to coordinate shifts.

This decomposition brings several theoretical advantages:

• •

exact IR cancellation is realized nonperturbatively;

• •

the framework naturally recovers analytically derived IR-resummed models for BAO damping;

• •

the exponential damping observed in cross spectra between pre- and post-reconstruction fields emerges directly from the displacement-mapping term.

Furthermore, the ULPT formulation is inherently extensible. The effects of Lagrangian galaxy bias can be consistently encoded through additive operators in the Jacobian deviation. The structure of RSD is incorporated by modifying both the displacement vector and the Jacobian to account for velocity components along the line of sight. Finally, reconstruction enters the formalism as a shift applied solely to the displacement vector, leaving the Jacobian deviation unaffected. As such, ULPT provides a unified theoretical structure that allows bias, RSD, and reconstruction to be incorporated consistently within the same formalism.

While ULPT offers a theoretically robust foundation, its practical performance in modeling broadband power spectra, especially beyond the BAO scale, remains to be tested. This is of particular importance because broadband features also encode cosmological information beyond that captured by the BAO peak.

In this paper, we assess the predictive accuracy of ULPT by computing the one-loop matter power spectrum in real space and comparing it with both SPT and high-precision numerical predictions from the Dark Emulator [22]. To further validate the robustness of our results, we perform an independent cross-check using the Euclid Emulator 2 [23], which has been carefully calibrated over a broad range of cosmological models.

In addition, we develop a fast and numerically stable algorithm for evaluating ULPT power spectra, based on efficient computation of the displacement-mapping integrals via Hankel transforms. This implementation provides a computationally efficient foundation for applying ULPT to future large-scale structure analyses.

Throughout this work, we adopt the fiducial cosmology implemented in the Dark Emulator suite, which is consistent with the Planck 2015 best-fit $\Lambda$CDM model [24]. The cosmological parameters are specified as follows: physical baryon density $\omega_{b}\equiv\Omega_{b}h^{2}=0.02225$, physical cold dark matter density $\omega_{c}\equiv\Omega_{c}h^{2}=0.1198$, dark energy density $\Omega_{\mathrm{de}}=0.6844$, scalar spectral index $n_{s}=0.9645$, amplitude of primordial curvature perturbations $\ln(10^{10}A_{s})=3.094$, dark energy equation-of-state parameter $w_{0}=-1$, and total neutrino mass $\sum m_{\nu}=0.06\,\mathrm{eV}$. The Hubble parameter is then determined to be $h=0.6727$ from the flatness condition.

This paper is organized as follows. In Sec. II, we review the theoretical foundation of ULPT, emphasizing its decomposition into Jacobian deviation and displacement-mapping effects. Section III presents a fast and numerically stable algorithm for evaluating the convolution integrals arising in the ULPT formulation. In Sec. IV, we compute the displacement correlation functions, which serve as key inputs to these integrals. Section V describes the one-loop computation of the source power spectrum using a modified FAST-PT algorithm. In Sec. VI, we assess the convergence of the displacement-mapping expansion and determine an optimal truncation scheme. Section VII analyzes the relative roles of the source term and the displacement-mapping factor in shaping the full one-loop power spectrum. In Sec. VIII, we summarize the simulation-based emulators used for validation, along with their parameter coverage. Section IX focuses on the nonlinear damping of the BAO feature. The main results are presented in Sec. X, where ULPT predictions are compared with emulator outputs across various cosmological models. Section XI outlines future directions, including extensions to biased tracers, RSD, and reconstruction, as well as two-loop improvements. We conclude in Sec. XII with a brief summary.

II Unified Lagrangian Formalism

In this section, we briefly review the core structure of ULPT [21]. Throughout this paper, we omit the explicit time (or redshift) dependence of physical quantities for notational simplicity, unless stated otherwise.

II.1 ULPT formulation of the Density Contrast

Let $\rho(\bm{x})$ denote the dark matter density field at Eulerian position $\bm{x}$. The density contrast $\delta(\bm{x})$, defined with respect to the mean density $\bar{\rho}$, is given by

[TABLE]

The mapping between Lagrangian coordinates $\bm{q}$ and Eulerian coordinates $\bm{x}$ is described by

[TABLE]

where $\bm{\Psi}(\bm{q})$ is the Lagrangian displacement vector. The Jacobian relation between the volume elements is given by

[TABLE]

Mass conservation implies

[TABLE]

From this, the density contrast evaluated at the Eulerian position $\bm{x}=\bm{q}+\bm{\Psi}(\bm{q})$ can be written as

[TABLE]

To express the Eulerian density contrast $\delta(\bm{x})$ in terms of Lagrangian variables, we start from the identity

[TABLE]

where $\delta_{\rm D}$ denotes the three-dimensional Dirac delta function. We then apply the coordinate transformation $\bm{x}^{\prime}=\bm{q}+\bm{\Psi}(\bm{q})$, yielding

[TABLE]

where we have defined the Jacobian deviation as

[TABLE]

Taking the Fourier transform of Eq. (7), we obtain the expression for the Fourier-space density contrast:

[TABLE]

where we denote Fourier-transformed quantities with a tilde.

Equation (9) provides the starting point for Unified Lagrangian Perturbation Theory (ULPT) in real space. In this formulation, the density contrast is explicitly decomposed into two physically distinct components:

• •

Jacobian deviation: $\delta_{\rm J}(\bm{q})$ characterizes the intrinsic linear and nonlinear growth of structure.

• •

Displacement-mapping effect: the exponential factor $e^{-i\bm{k}\cdot\bm{\Psi}(\bm{q})}$ describes the nonlinear coordinate remapping of the Jacobian field induced by displacements fields.

This structural decomposition, which separates the density contrast into the Jacobian deviation and the displacement-mapping effect, can be extended to a variety of observational contexts beyond the case of dark matter in real space. Specifically:

• •

For biased tracers, an additional contribution representing the biased fluctuation is linearly added to the Jacobian deviation, while the displacement-mapping effect retains the original structure.

• •

For RSD, both the Jacobian deviation and the displacement vector are modified to incorporate the effects of line-of-sight peculiar velocities. These modifications preserve the overall structure of the ULPT formulation while encoding the anisotropic nature of RSD.

• •

For density-field reconstruction, the displacement vector is modified by adding a reconstruction-induced shift field, whereas the Jacobian deviation is left unaltered. This structure ensures that the linear density field is preserved under reconstruction, as theoretically expected.

In this paper, we focus on the case of dark matter in real space for clarity.

II.2 ULPT Formulation of the Power Spectrum

The power spectrum in ULPT is computed as the Fourier transform of the two-point ensemble average of the density contrast. Starting from the expression

[TABLE]

where $\bm{r}=\bm{q}-\bm{q}^{\prime}$ is the separation vector between two Lagrangian positions, and $\langle\cdots\rangle$ denotes the ensemble average.

Defining

[TABLE]

the power spectrum can be rewritten using cumulant expansion as

[TABLE]

where $\langle\cdots\rangle_{\rm c}$ denotes the connected part of the ensemble average.

The first factor, $\left\langle e^{X}\right\rangle$, depends only on the displacement field and is statistically uncorrelated with the Jacobian deviation field. We refer to this as the displacement-mapping factor. Its cumulant expansion takes the form

[TABLE]

where the exponent is defined as

[TABLE]

When considering the difference of a single displacement field, the displacement variance $\overline{\Sigma}(\bm{k})$ is defined as the zero-separation limit:

[TABLE]

We define the remaining part of the integrand as the source correlation function:

[TABLE]

Expanding $e^{X}$ perturbatively, we obtain

[TABLE]

where we used $\langle\delta_{\rm J}\rangle=0$. This expression shows that the leading (i.e., linear-order) contribution to the source correlation function $\xi_{\rm J}$ arises from the two-point correlation of the Jacobian deviation. Higher-order contributions to $\xi_{\rm J}$ originate from two distinct sources: the nonlinear structure inherent in the Jacobian deviation’s two-point correlation itself, and cross-correlations between the Jacobian deviation and the displacement field.

Combining the displacement-mapping factor and the source correlation function, the full power spectrum in ULPT is compactly expressed as

[TABLE]

Although already introduced in the course of the derivation, the two physically distinct contributions to the ULPT power spectrum are restated here to highlight their conceptual significance:

• •

Source correlation function $\xi_{\rm J}(r)$: the two-point correlation of the Jacobian deviation field, incorporating both linear and nonlinear intrinsic effects, as well as cross-correlations with the displacement field.

• •

Displacement-mapping factor $e^{-\overline{\Sigma}(\bm{k})+\Sigma(\bm{k},\bm{r})}$: an exponential modulation capturing the nonlinear influence of the displacement field on the clustering pattern. This factor isolates the convective coordinate remapping effect that is statistically uncorrelated with the intrinsic source field.

This structural separation plays a central role throughout the ULPT framework and enables a consistent treatment of IR effects.

We also define the source power spectrum as the Fourier transform of the source correlation function:

[TABLE]

With this definition, the total power spectrum can be systematically decomposed into the source power spectrum and an additional modulation term arising from the displacement-mapping factor acting on the source correlation function:

[TABLE]

where the subscript “DM” stands for the “displacement-mapping” contribution. The explicit expression for $P_{\rm DM}(\bm{k})$ reads

[TABLE]

Alternatively, one may classify the contributions to the power spectrum based on whether they arise from convolution integrals associated with the displacement-mapping factor. In this case, the power spectrum is decomposed as

[TABLE]

where “CF” denotes the “convolution-free“ part and ”CC“ denotes the ”convolution-containing“ part. These are given explicitly by

[TABLE]

Since we are working in real space and assuming statistical isotropy, both $P(k)$ and $\xi_{\rm J}(r)$ depend only on the magnitudes of their respective arguments, $k=|\bm{k}|$ and $r=|\bm{r}|$. Furthermore, the displacement variance $\Sigma(\bm{k},\bm{r})$ depends only on three scalar quantities: $k$, $r$, and the scalar product $\hat{\bm{k}}\cdot\hat{\bm{r}}$. This property significantly simplifies the numerical evaluation of the convolution integrals, particularly when using spherical coordinates or Hankel transforms, as will be explicitly demonstrated in Section III.

III Fast Algorithms for ULPT Convolution Integrals

We now present an efficient approximation scheme for evaluating the convolution integrals arising from the displacement-mapping factor in the ULPT formulation. This scheme closely follows the approach developed in Ref. [4] for computing the ZA power spectrum as well as LPT power spectra including higher-order corrections, and recasts the relevant integrals in terms of Hankel transforms involving spherical Bessel functions.

Let $f(r)$ be an arbitrary spherically symmetric function. The Hankel transform of order $\ell$, commonly used in cosmology, is defined as

[TABLE]

where $j_{\ell}(x)$ denotes the spherical Bessel function of order $\ell$. The corresponding inverse transform is given by

[TABLE]

These transforms can be efficiently evaluated using FFTLog-based algorithms [25]. In our implementation, we make use of the publicly available mcfit package [26].

The goal of this paper is to compute the ULPT power spectrum including one-loop corrections. Our approximation exactly reproduces the one-loop result of SPT, while resumming higher-order terms that originate from the convolution structure of the displacement-mapping factor.

To this end, we expand the exponent of the displacement-mapping factor at linear order and compute the source correlation function up to one-loop order:

[TABLE]

At linear order, the source correlation function $\xi_{\rm J}^{(\mathrm{lin})}(r)$ coincides with the standard matter correlation function, and we denote it simply by $\xi_{\mathrm{lin}}(r)$ in the following.

The corresponding source power spectrum reads

[TABLE]

where $P_{\rm lin}(k)$ is the linear dark matter power spectrum.

The linear-order displacement variance is given by

[TABLE]

Performing the angular integrals analytically, it becomes

[TABLE]

with

[TABLE]

Here, $\mathcal{L}_{\ell}$ is the Legendre polynomial of order $\ell$. The variance $\overline{\Sigma}(k)$ is defined via the zero-separation limit:

[TABLE]

leading to

[TABLE]

Substituting these results into the general expression for the ULPT power spectrum, we obtain:

[TABLE]

where $\mu=\hat{\bm{k}}\cdot\hat{\bm{r}}$.

We now decompose the total power spectrum into the CF and CC contributions, as defined in Eq. (23):

[TABLE]

To further simplify the CC term, we expand the angular dependence in the exponent:

[TABLE]

which leads to the decomposition

[TABLE]

with the zeroth and higher-order terms given by

[TABLE]

Performing the angular integrals analytically, we recast these expressions as Hankel transforms:

[TABLE]

The angular integrals $I^{[n]}(k,r)$ are defined by

[TABLE]

and for $n=1$ to $4$ they evaluate to

[TABLE]

These results demonstrate that the convolution contributions can be systematically evaluated using Hankel transforms.

IV Correlation Functions of the Displacement Field

In this section, we present the numerical evaluation of the displacement correlation functions $\sigma^{2}_{\ell}(r)$ for $\ell=0$ and $\ell=2$, as defined in Eq. (31). These quantities characterize the scale-dependent contributions of the displacement field and play a central role in evaluating the convolution integral induced by the displacement-mapping factor in ULPT.

The results shown here correspond to redshift $z=0$, where nonlinear effects are most pronounced. The correlation functions are efficiently computed using the FFTLog algorithm [25, 26]. However, at small separations $r$, the finite $k_{\rm max}$ cutoff in the input linear power spectrum induces spurious oscillations in the numerical results. To suppress such artifacts, we adopt a conservative value of $k_{\rm max}=100\,h\,\mathrm{Mpc}^{-1}$, which stabilizes the evaluation down to sub-megaparsec scales.

Figure 1 shows that numerical instabilities become significant for $r\lesssim 0.75\,h^{-1}\mathrm{Mpc}$, where both $\sigma_{0}^{2}(r)$ and $\sigma_{2}^{2}(r)$ exhibit artificial oscillations. Nevertheless, the limiting values at $r=0$ are known analytically: $\sigma_{0}^{2}(0)=\bar{\sigma}^{2}$, as given in Eq. (32), and $\sigma_{2}^{2}(0)=0$. We therefore apply a linear interpolation between $r=0$ and $r=0.75\,h^{-1}\mathrm{Mpc}$ to smoothly connect the stable analytic values at the origin to the numerically computed values at the cutoff scale. This procedure effectively eliminates the unphysical FFT-induced fluctuations in the small-$r$ regime and ensures numerical stability in subsequent integrations.

Figure 2 shows the displacement correlation functions $\sigma_{\ell}^{2}(r)$ over a wide range of comoving separations, spanning from $r=0.1\,h^{-1}\mathrm{Mpc}$ to $10^{4}\,h^{-1}\mathrm{Mpc}$. For $r<0.75\,h^{-1}\mathrm{Mpc}$, the curves are linearly interpolated between the numerical values at $r=0.75\,h^{-1}\mathrm{Mpc}$ and the analytic limits at $r=0$, following the same procedure used in Fig. 1.

As required by definition, the monopole component $\sigma_{0}^{2}(r)$ asymptotically approaches $\bar{\sigma}^{2}$ at small separations, but rapidly decreases between $r\sim 10\,h^{-1}\mathrm{Mpc}$ and $r\sim 100\,h^{-1}\mathrm{Mpc}$, eventually approaching zero toward large scales. In contrast, the quadrupole component $\sigma_{2}^{2}(r)$ exhibits a broad negative trough around $r\sim 50\,h^{-1}\mathrm{Mpc}$, and vanishes both in the small- and large-scale limits.

When comparing the relative amplitudes of the monopole and quadrupole components, it is important to note that the monopole contribution always appears in the combination $\sigma_{0}^{2}(r)-\bar{\sigma}^{2}$ in the convolution integrals. Therefore, the relevant quantity to assess the hierarchy is the ratio $|\sigma_{2}^{2}(r)/[\sigma_{0}^{2}(r)-\bar{\sigma}^{2}]|$. This ratio remains below $\sim 0.2$ for $r\gtrsim 10\,h^{-1}\mathrm{Mpc}$, indicating that the quadrupole contribution is significantly suppressed relative to the effective monopole term over the scales of interest. This provides numerical support for the perturbative expansion scheme adopted in Sec. III, in which the $\sigma_{2}^{2}$-dependent terms are systematically expanded.

V Fast Evaluation of the One-Loop Source Power Spectrum

In this section, we develop a fast and numerically stable algorithm for computing the one-loop source power spectrum $P_{\rm J}$ within the ULPT framework. We begin with a brief overview of the publicly available FAST-PT algorithm [27, 28] in Sec. V.1, and in Sec. V.2 we describe how the same algorithm can be applied to ULPT by replacing the standard SPT kernel with its ULPT counterpart.

V.1 Review of the FAST-PT Algorithm

FAST-PT [27] is a Python-based tool designed to efficiently compute convolution integrals which arise in cosmological perturbation theory. The method is particularly well-suited for applications involving the decomposition of angular dependencies in mode-coupling kernels into Legendre polynomials, enabling analytic evaluation of angular integrals. The remaining radial integrals, expressed as Hankel transforms involving spherical Bessel functions, can be rapidly evaluated using FFTLog techniques. Consequently, the overall computational cost is reduced to $\mathcal{O}(N\log N)$.

For the one-loop matter power spectrum in SPT, the total correction is decomposed as

[TABLE]

where $P^{(22)}$ represents a double convolution of the linear power spectrum and $P^{(13)}$ captures the coupling between linear and third-order terms. In FAST-PT, these terms are regularized to avoid IR divergences:

[TABLE]

where $\bar{\sigma}^{2}$ is the displacement variance given by Eq. (32). This procedure preserves the total one-loop correction while canceling unphysical IR contributions in each term.

For the $P^{(13)}_{\text{reg}}$ term, the regularized expression takes the form

[TABLE]

where the kernel $Z_{\mathrm{reg}}(r)$ is given by

[TABLE]

Although this integral naively appears to scale as $\mathcal{O}(N^{2})$, it can be rewritten as a discrete convolution using logarithmic variable transformations and then evaluated in $\mathcal{O}(N\log N)$ operations using FFT techniques.

The $P^{(22)}$ term poses additional challenges due to its dependence on two distinct wavevectors. FAST-PT circumvents this by transforming the convolution to configuration space using Legendre expansions. In configuration space, $P^{(22)}$ is evaluated as

[TABLE]

where

[TABLE]

and the second-order kernel is

[TABLE]

with $\mu_{12}=\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}$ denoting the cosine of the angle between the wavevectors.

The kernel squared can be expanded in Legendre polynomials:

[TABLE]

where the coefficients $c_{\alpha\beta\ell}$ are listed in Table 1. This expansion allows the angular integrations to be performed analytically. The resulting correlation function is then written as

[TABLE]

with each term defined by

[TABLE]

These integrals are efficiently evaluated using Hankel transforms, leading to substantial computational speedup compared to brute-force multidimensional integration.

Among these terms, the component with $(\alpha,\beta,\ell)=(2,-2,0)$ is known to be IR-sensitive. To address this, FAST-PT applies a regularization procedure by analytically subtracting the divergent component at the integrand level. The regularized form of $J_{2,-2,0}(r)$ is given by

[TABLE]

where the subtraction of unity ensures convergence in the IR limit, since $[j_{0}(k_{2}r)-1]=(k_{2}r)^{2}/6+\cdots$ for small $k_{2}$.

The regularization terms such as $k^{2}\bar{\sigma}^{2}P_{\mathrm{lin}}(k)$ correspond to the effect of long-wavelength bulk flows. These displacements induce apparent power at large scales, which must cancel in the final, physically meaningful power spectrum due to statistical homogeneity and Galilean invariance. While the total one-loop power spectrum does exhibit this cancellation, individual components such as $P_{13}$ and $P_{22}$ do not, and must therefore be regularized explicitly.

From a numerical standpoint, even though the linear matter power spectrum in a $\Lambda$CDM cosmology behaves as $P_{\mathrm{lin}}(k)\propto k^{n_{s}}$ with $n_{s}\simeq 0.96$ in the $k\to 0$ limit, and thus poses no true divergence, the FFT-based approach used in FAST-PT introduces a logarithmic discretization. This results in a log-periodic effective power spectrum with artificial power at extremely large or small scales. Consequently, these spurious contributions must be analytically subtracted before integration to avoid numerical instabilities.

V.2 FAST-PT Implementation with ULPT-Specific Kernels

In our implementation, we retain the core algorithmic structure of FAST-PT and adapt it to the ULPT framework by substituting the standard $F_{2}$ kernel with the ULPT-specific version. No modification to the underlying numerical techniques is required.

In the ULPT formulation, the SPT one-loop power spectrum can be decomposed, following Eq. (20), into two distinct contributions:

[TABLE]

where the displacement-mapping (DM) contribution is given by

[TABLE]

Accordingly, it is not necessary to directly evaluate Eq. (17), which defines the source correlation function. Instead, the one-loop source power spectrum can be obtained by subtracting the displacement-mapping contribution from the known one-loop result in SPT:111At two-loop order, the source power spectrum is given by

$P_{\rm J}^{(2\text{-loop})}(k)=P^{(2\text{-loop})}(k)-P_{\rm DM}^{(2\text{-loop})}(k),$

(56)

where the two-loop DM contribution can be expressed in terms of the linear and one-loop results as

$\displaystyle P_{\rm DM}^{(2\text{-loop})}(k)$ $\displaystyle=\int d^{3}r\,e^{-i\bm{k}\cdot\bm{r}}\Bigg\{\left[-\overline{\Sigma}^{\rm(lin)}(\bm{k})+\Sigma^{\rm(lin)}(\bm{k},\bm{r})\right]\xi_{\rm J}^{(1\text{-loop})}(r)$

$\displaystyle\quad+\frac{1}{2}\left[-\overline{\Sigma}^{\rm(lin)}(\bm{k})+\Sigma^{\rm(lin)}(\bm{k},\bm{r})\right]^{2}\xi_{\rm lin}(r)$

$\displaystyle\quad+\left[-\overline{\Sigma}^{(1\text{-loop})}(\bm{k})+\Sigma^{(1\text{-loop})}(\bm{k},\bm{r})\right]\xi_{\rm lin}(r)\Bigg\}.$

[TABLE]

As in the SPT framework in Eq. (42), Eq. (55) can be decomposed into two components: the so-called 13-type and 22-type terms, corresponding to the contributions without and with mode-coupling integrals, respectively. In the ULPT formulation, the term involving $-\overline{\Sigma}^{\rm(lin)}(\bm{k})$ corresponds to the 13-type contribution, while the term involving $\Sigma^{\rm(lin)}(\bm{k},\bm{r})$ corresponds to the 22-type contribution.

Using Eqs. (29) and (32), the individual DM terms can be expressed as:

[TABLE]

where

[TABLE]

Accordingly, the 13-type and 22-type components of the source power spectrum can be defined by subtracting the DM contributions from the SPT terms:

[TABLE]

where

[TABLE]

It is important to note that the source 13-type term, $P_{\rm J}^{(13)}$, coincides with the regularized version of the SPT 13-type contribution as defined in FAST-PT:

[TABLE]

Therefore, for this term, we can directly adopt the result provided by FAST-PT.

The 22-type component of the DM contribution, $P_{\rm DM}^{(22)}$, is computed in configuration space following the methods outlined in Sec. V.1. Its inverse Fourier transform, given in Eq. (59), is decomposed in terms of Legendre polynomials as follows:

[TABLE]

where the functions $J_{\alpha\beta\ell}(r)$ are defined in Eq. (52). Substituting this result into Eq. (61), the inverse Fourier transform of the source 22-type component becomes

[TABLE]

where the coefficients $c_{{\rm J},\alpha\beta\ell}$ are listed in Table 1.

Notably, in $P_{\rm J}^{(22)}$, the IR-sensitive term $J_{2,-2,0}(r)$ appears with a zero coefficient. As a result, the source power spectrum $P_{\rm J}^{(22)}$ is automatically regularized.

This natural regularization arises from the fact that IR-sensitive contributions, such as those due to long-wavelength bulk flows, are fully captured by the displacement-mapping factor. These effects are thus separated from the source power spectrum and handled analytically.

In particular, numerical instabilities in the IR regime are effectively absorbed into the displacement-mapping exponent, $\sigma_{0}^{2}(r)$ in Eq. (31). However, because $\sigma_{0}^{2}(r)$ always appears in the combination $\sigma_{0}^{2}(r)-\bar{\sigma}^{2}$, the integrands involve the form $[j_{0}(kr)-1]$, ensuring convergence in the IR limit.

As shown in Sec. IV, in this work, the function $\sigma_{0}^{2}(r)$ is computed directly using FFTLog, implemented via the mcfit package [26]. The zero-separation variance $\bar{\sigma}^{2}$ is evaluated independently. As a result, we do not apply the analytic regularization procedure described above, which may lead to numerical instability at small values of $r$. To mitigate this issue, we extend the linear input power spectrum up to a maximum wavenumber of $k_{\rm max}=100\,h\,{\rm Mpc}^{-1}$. Furthermore, for $r<0.75\,h^{-1}{\rm Mpc}$, we linearly interpolate $\sigma_{0}^{2}(r)$ between its value at $r=0.75\,h^{-1}{\rm Mpc}$ and the limiting value $\bar{\sigma}^{2}$ at $r=0$. Similar numerical instabilities are also observed in $\sigma_{2}^{2}(r)$, and we apply the same interpolation procedure to this quantity as well. Since our final results target scales up to $k\simeq 0.4\,h\,{\rm Mpc}^{-1}$, we expect that contributions from $r<0.75\,h^{-1}{\rm Mpc}$ do not significantly affect the accuracy.

VI Convergence of the Expansion Scheme for Displacement-Mapping Integrals

In this section, we validate the convergence of the expansion method for evaluating the convolution integrals induced by the displacement-mapping factor, as proposed in Sec. III. Specifically, we perform a detailed analysis at redshift $z=0$, where nonlinear effects are maximized, and identify the number of expansion terms in $P_{\rm CC}$ required to achieve sufficient accuracy. Since the nonlinear contribution becomes smaller at higher redshifts, convergence confirmed at $z=0$ guarantees even better accuracy at $z>0$, thereby justifying the use of the same approximation scheme at higher redshifts.

As shown in Eq. (22), the total power spectrum can be decomposed into two distinct components based on whether they involve convolution integrals: the convolution-free term $P_{\rm CF}$, and the convolution-containing term $P_{\rm CC}$. The upper panel of Fig. 3 displays $P_{\rm CF}/P_{\rm nw}$ and the cumulative sums $\sum_{n=0}^{N_{\rm C}}P^{[n]}_{\rm CC}/P_{\rm nw}$ for $N_{\rm C}=0,1,2,3,4$, where $P_{\rm nw}$ denotes the no-wiggle linear power spectrum without BAO [29, 30, 31].

In accordance with Eq. (23), $P_{\rm CF}$ asymptotically approaches the linear power spectrum on large scales and exhibits exponential suppression on small scales due to the displacement variance. Conversely, $P_{\rm CC}$ is subdominant on large scales but dominates the nonlinear shape of the power spectrum on small scales. Together, these two components constitute the full one-loop ULPT power spectrum.

Regarding the convergence of the $P_{\rm CC}$ expansion, we find that truncating at $N_{\rm C}=1$ already captures the majority of the nonlinear correction. As $N_{\rm C}$ increases to 2, 3, and 4, the contribution from higher-order terms decreases monotonically, indicating good convergence of the series expansion.

To assess this quantitatively, the lower panel of Fig. 3 shows the relative error with respect to the reference result obtained by summing up to $N_{\rm C}=4$, which can be regarded as sufficiently converged to represent the total ULPT power spectrum:

[TABLE]

where we define $P_{\rm ULPT}\equiv P_{\rm CF}+\sum_{n=0}^{N_{\rm C}=4}P_{\rm CC}^{[n]}$. The figure demonstrates that for $N_{\rm C}=1$, the deviation reaches about $0.4\%$ at $k=0.2\,h\,{\rm Mpc}^{-1}$, and exceeds $2\%$ at $k=0.4\,h\,{\rm Mpc}^{-1}$. Increasing the expansion to $N_{\rm C}=2$, the error reduces to $0.1\%$ at $k=0.2\,h\,{\rm Mpc}^{-1}$ and to approximately $1.2\%$ at $k=0.4\,h\,{\rm Mpc}^{-1}$. At $N_{\rm C}=3$, the approximation achieves sub-percent accuracy over the entire range, with a maximum error of only $0.025\%$ at $k=0.4\,h\,{\rm Mpc}^{-1}$.

Although not shown in the figure, we confirm that the convergence improves further at higher redshifts. For example, at $z=0.5$, the error with $N_{\rm C}=2$ is below $0.03\%$ at $k=0.2\,h\,{\rm Mpc}^{-1}$ and around $0.4\%$ at $k=0.4\,h\,{\rm Mpc}^{-1}$. At $z=1.0$, the respective errors decrease further to $0.01\%$ and $0.13\%$, confirming the robustness of the expansion scheme.

Since current galaxy surveys such as the Dark Energy Spectroscopic Instrument (DESI) [32] typically utilize power spectrum measurements up to $k\simeq 0.2\,h\,{\rm Mpc}^{-1}$, we conclude that truncating the expansion at $N_{\rm C}=2$ suffices for most practical applications in cosmological data analysis. Nevertheless, to adopt a conservative stance throughout this paper, we present numerical results based on the more stringent choice of $N_{\rm C}=3$.

The full pipeline begins by obtaining the linear power spectrum from the Dark Emulator, followed by evaluating the one-loop correction using FAST-PT as described in Sec. V, and performing the convolution integrals of the displacement-mapping factor in Sec. III to compute the final one-loop ULPT power spectrum. When executed on a standard MacBook Pro, the runtime for each stage of the computation is as follows: retrieving the linear power spectrum takes approximately 0.73 seconds, while the FAST-PT evaluation requires about 0.016 seconds. The convolution integrals for the displacement-mapping factor take 1.0 seconds for $N_{\rm C}=2$, 1.75 seconds for $N_{\rm C}=3$, and 2.7 seconds for $N_{\rm C}=4$. Consequently, the total runtime amounts to approximately 1.8 seconds, 2.5 seconds, and 3.4 seconds for $N_{\rm C}=2$, 3, and 4, respectively. These timings are measured using 450 logarithmically spaced $k$-points in the range $10^{-4}\,h\,\mathrm{Mpc}^{-1}\leq k\leq 100\,h\,\mathrm{Mpc}^{-1}$, which provides sufficient resolution and coverage for computing the two-point correlation function via mcfit. Since a truncation at $N_{\rm C}=2$ or $3$ is sufficient to achieve sub-percent accuracy in the convolution expansion, we conclude that the full ULPT power spectrum can be computed in approximately 2 seconds with practically negligible loss of precision. This level of computational efficiency highlights the practical applicability of the method for rapid and repeated model evaluations.

VII Source Power Spectrum vs. Displacement-Mapping Factor

In this section, we investigate the impact of the displacement-mapping factor on the one-loop power spectrum. Specifically, we compare the one-loop source power spectrum $P_{\rm J}$ with the full one-loop ULPT power spectrum $P_{\rm ULPT}$, which includes convolution integrals arising from the displacement-mapping factor.

For reference, we also compute the one-loop SPT power spectrum $P_{\rm SPT}$. While the computation follows standard procedures implemented in FAST-PT, its structural relationship to the ULPT formulation is clarified in Eq. (54): both share the same one-loop source power spectrum, but differ in how the displacement-mapping factor is treated. In SPT, this factor is expanded to $\mathcal{O}(P_{\rm lin})$, whereas ULPT retains its full exponential form. Comparing $P_{\rm SPT}$ and $P_{\rm ULPT}$ thus allows us to isolate and quantify the consequences of expanding versus preserving the displacement-mapping factor.

The upper panels of Fig. 4 show the ratios $P_{\rm J}/P_{\rm nw}$, $P_{\rm ULPT}/P_{\rm nw}$, and $P_{\rm SPT}/P_{\rm nw}$ at redshifts $z=0$, $0.5$, and $1.0$, where $P_{\rm nw}$ denotes the no-wiggle linear power spectrum. All spectra include one-loop corrections.

We observe that $P_{\rm ULPT}$ and $P_{\rm SPT}$ begin to diverge around $k=0.1\,h\,{\rm Mpc}^{-1}$ at redshift $z=0.0$ (left panel of Fig. 4), where nonlinear effects are most pronounced. As $k$ increases, the SPT prediction grows rapidly and eventually overshoots, while the ULPT result shows a more moderate nonlinear growth. This divergence arises entirely from the treatment of the displacement-mapping factor: both approaches share the same one-loop source term $P_{\rm J}$, but SPT expands the exponential, whereas ULPT retains its full nonperturbative form. The excessive growth of $P_{\rm SPT}$ at small scales reflects the breakdown of the expansion and underscores the importance of preserving the exponential structure.

Comparing $P_{\rm ULPT}$ and $P_{\rm J}$, we find that the overall shape of the ULPT power spectrum is primarily governed by the source term. The displacement-mapping factor modulates its amplitude across a broad range of scales, either suppressing or enhancing it depending on redshift and wavenumber.

To quantify this modulation, the lower panels of Fig. 4 plot the ratio $P_{\rm ULPT}/P_{\rm J}$ at redshifts $z=0$, $0.5$, and $1.0$. The results show that the displacement-mapping factor modifies the source power spectrum by a few percent, with the sign and magnitude of the modulation depending on both scale and redshift. At $z=0$, the effect is purely suppressive, reaching a level of 3–5% up to $k\simeq 0.4\,h\,{\rm Mpc}^{-1}$. At higher redshifts ($z=0.5$ and $z=1.0$), the modulation exhibits a qualitatively different trend: while a suppression of 2–3% persists at large scales ($k\lesssim 0.2\,h\,{\rm Mpc}^{-1}$), the displacement-mapping factor enhances the source power spectrum by up to 5% at smaller scales ($k\gtrsim 0.2\,h\,{\rm Mpc}^{-1}$). This enhancement becomes more pronounced at higher redshift, with the most significant amplification observed at $z=1.0$. These results indicate a redshift-dependent transition in the displacement-mapping effect, from overall suppression at low redshift to scale-dependent enhancement at high redshift.

Importantly, these ratios exhibit oscillatory features at the BAO scale. This is because $P_{\rm ULPT}$ and $P_{\rm J}$ encode different BAO behavior. In particular, the displacement-mapping factor damps the BAO signal in addition to modulating the broadband shape, resulting in residual oscillations in their ratio.

A more detailed analysis of nonlinear BAO damping will be presented in Sec. IX.

VIII Simulation-Based Emulators: Dark Emulator & Euclid Emulator 2

VIII.1 Overview and Motivation

To evaluate the validity of ULPT predictions across a wide range of cosmological models, we employ simulation-based emulators that provide fast and accurate estimates of the nonlinear matter power spectrum. In this work, we primarily use the Dark Emulator [22], which was originally developed to compute halo statistics such as halo–halo and halo–matter power spectra with high efficiency and accuracy. Although Dark Emulator also provides predictions for the nonlinear matter power spectrum, this component has received less validation in the literature because the emulator was designed mainly for halo-related observables.

In the present study, we use only the nonlinear matter power spectrum from the Dark Emulator, without utilizing its halo modules. This choice is made to ensure consistency in future developments, where ULPT may be extended to compute one-loop corrections for halo statistics. Using the same emulator for both matter and halo predictions will allow for seamless validation within a unified framework.

However, given the relatively limited validation of the matter power spectrum component in the Dark Emulator, we perform an independent cross-check using a second emulator, Euclid Emulator 2 [23], which has been carefully validated to achieve sub-percent accuracy across a wide parameter range, including models with massive neutrinos and time-dependent dark energy.

VIII.2 Simulation Specifications

The Dark Emulator is based on the Dark Quest $N$-body simulation suite, which covers 101 flat $w$CDM cosmologies. Each model is simulated with a high-resolution run using $2048^{3}$ particles in a $2\,h^{-1}\mathrm{Gpc}$ box. To suppress statistical noise, the fiducial Planck cosmology is simulated with 14 independent realizations, enabling a robust estimate of cosmic variance. Accordingly, most of our comparisons between ULPT and emulator predictions are first performed using this fiducial model.

In contrast, Euclid Emulator 2 is trained on 127 cosmological models, each simulated with a paired-and-fixed (P+F) initial condition scheme. For each cosmology, two simulations are performed using $3000^{3}$ particles in a $1\,h^{-1}\mathrm{Gpc}$ box, resulting in a total of 254 simulations. The P+F method is known to significantly suppress sample variance, particularly on large scales. This design suggests that Euclid Emulator 2 exhibits inherently lower sample variance than Dark Emulator across most relevant cosmological models.

VIII.3 Cross-Validation Between Emulators

To assess the reliability of the Dark Emulator in cosmologies beyond the fiducial model, we perform a direct cross-comparison with Euclid Emulator 2. Using the common parameter region specified in Table 2, we randomly sample 100 cosmological parameter sets that lie within the overlapping validity range of both emulators. For each model, we compute the nonlinear matter power spectrum from both emulators and examine their fractional difference.

Figure 5 displays the relative deviation between the two emulator predictions for these 100 cosmological models at redshift $z=0$, where nonlinear effects are strongest. While some individual models show deviations exceeding 5% at certain scales, the mean deviation remains within $\pm 1\%$, indicating that the two emulators are statistically consistent across the sampled parameter space.

The substantial scatter observed for individual cosmologies is attributed primarily to the larger statistical variance inherent in the Dark Emulator, as anticipated from the simulation specifications in Sec. VIII.2. In particular, for scales $k\lesssim 0.05\,h\,\mathrm{Mpc}^{-1}$, even the mean deviation between the two emulators can reach up to 2%, again reflecting the sample variance in Dark Emulator.

As we demonstrate in Sec. X.2, the ULPT predictions on large scales agree well with the Euclid Emulator 2 results, confirming that the perturbative framework remains reliable in this regime.

IX Nonlinear Damping of BAO in the Correlation Function

In this section, we investigate how ULPT describes the nonlinear damping of the BAO feature in the two-point correlation function. The behavior at smaller scales is deferred to Sec. X.1.

IX.1 Displacement Mapping and Nonlinear BAO

To elucidate the role of the displacement-mapping factor in modeling BAO suppression, we compare several correlation functions at one-loop order. Specifically, Fig. 6 presents four results: (i) the source correlation function $\xi_{\rm J}$, which excludes the displacement-mapping effect; (ii) the one-loop SPT correlation function $\xi_{\rm SPT}$, where the displacement-mapping factor is expanded to leading order; (iii) the ULPT correlation function $\xi_{\rm ULPT}$, which retains the full exponential structure; and (iv) the nonlinear correlation function computed from the Dark Emulator. All computations are performed using the fiducial cosmological parameters introduced in Sec. I, for which the Dark Emulator exhibits minimal variance.

This comparison highlights the critical role of the displacement-mapping factor. The source-only result $\xi_{\rm J}$ exhibits no suppression of the BAO feature. The SPT result $\xi_{\rm SPT}$, obtained by expanding the displacement-mapping factor, distorts the shape of the BAO peak. In contrast, the ULPT prediction $\xi_{\rm ULPT}$ accurately reproduces the nonlinear damping effect and shows excellent agreement with the emulator result. These findings emphasize the importance of preserving the full exponential structure of the displacement-mapping factor to correctly capture the nonlinear behavior of BAO.

IX.2 Structure of One-Loop Corrections in ULPT

To further clarify the internal structure of ULPT, the top panel of Fig. 7 compares the linear correlation functions from SPT and ULPT with the full one-loop ULPT prediction. The linear ULPT correlation function is computed by evaluating the source correlation function at linear order while fully incorporating the displacement-mapping factor. Although this approximation already captures the general shape of BAO suppression, it slightly underestimates the peak amplitude compared to the one-loop ULPT and emulator results. This discrepancy indicates that one-loop corrections to the source term $\xi_{\rm J}$ mildly enhance the BAO peak.

The bottom panel of Fig. 7 compares the linear SPT result with the one-loop source correlation function. Since both exclude the displacement-mapping factor, their difference reflects purely the nonlinear evolution contained in $\xi_{\rm J}$. The one-loop source correlation function exhibits a small but noticeable increase in the BAO peak relative to the linear SPT prediction.

As shown in Eq. (60), the one-loop source power spectrum correction $P_{\rm J}^{(\text{1-loop})}$ consists of a 13-type term and a 22-type term. The 13-type term is proportional to the linear power spectrum and retains the BAO signal, while the 22-type term involves mode-coupling integrals that tend to smear out oscillatory features. Therefore, the enhancement of the BAO peak in $\xi_{\rm J}$ is primarily attributed to the 13-type term.

To gain further insight into the 22-type contribution, we now examine how mode-coupling integrals impact the nonlinear damping of the BAO feature. The BAO signal represents an isotropic oscillation, and its structure can be smoothed through angular averaging inherent in the mode-coupling process. This smoothing arises from destructive interference between oscillatory modes with different phases. However, the efficiency of this phase cancellation depends nontrivially on the detailed form of the linear power spectrum and the specific structure of the coupling kernel. For example, if the SPT kernel $F_{2}$ is replaced by an arbitrary kernel, it is not guaranteed that the resulting integrals will fully suppress the BAO signal. In practice, for $\Lambda$CDM-like cosmologies, this smoothing is known to occur and provides a subdominant but non-negligible contribution to the total nonlinear power spectrum [33].

In particular, when the coupling kernel contains scale-dependent terms such as $(k_{1}/k_{2})^{2}$ or $k_{1}/k_{2}$, as listed in Table 1, the integrals become sensitive to specific phases in the infrared limit $k_{2}\to 0$. This sensitivity weakens the cancellation of BAO oscillations and results in a residual imprint in the final spectrum. Accordingly, such terms in the SPT 22-type contribution $P^{(22)}$ have been used in analytic models to characterize nonlinear BAO damping. The widely used IR-resummed model expresses the power spectrum as [7, 8]

[TABLE]

where $P_{\rm w}$ and $P_{\rm nw}$ are the oscillatory (wiggle) and smooth (no-wiggle) components of the linear power spectrum, respectively. The damping scale $\bar{\sigma}_{\rm BAO}^{2}$ is given by

[TABLE]

with $r_{\rm BAO}\sim 105\,h^{-1}\mathrm{Mpc}$ being the characteristic BAO scale. The appearance of spherical Bessel functions $j_{0}$ and $j_{2}$ in Eq. (67) reflects the $(k_{1}/k_{2})^{2}$ dependence of the relevant terms in $P^{(22)}$, particularly those corresponding to $(\alpha,\beta,\ell)=(2,-2,0)$ and $(2,-2,2)$ in Table 1. If all BAO-related contributions in the mode-coupling integrals are neglected, the resulting IR-resummed model reduces to a simplified form that does not involve spherical Bessel functions such as $j_{0}$ or $j_{2}$ [5].

In the ULPT framework, the 22-type source contribution $P_{\rm J}^{(22)}$ is automatically regularized as shown in Sec. V.2, and terms involving $(k_{1}/k_{2})^{2}$ do not appear, as shown in Table 1. These components are effectively absorbed into the displacement-mapping factor, which encapsulates the dominant nonlinear damping of the BAO. The remaining terms in $P_{\rm J}^{(22)}$ include only those with weaker scale dependence such as $k_{1}/k_{2}$ and scale-independent contributions. We thus conclude that the residual BAO signal in $P_{\rm J}^{(22)}$ is clearly subdominant compared to that in the 13-type source contribution $P_{\rm J}^{(13)}$, which is directly proportional to the linear power spectrum and retains the full oscillatory structure.

To verify this interpretation, we compute a modified version of the source power spectrum in which the linear power spectrum in the 13-type term is replaced by its no-wiggle counterpart:

[TABLE]

The corresponding correlation function, denoted the 1-loop $\xi_{\rm J,nw}$, is plotted in the lower panel of Fig. 7. This result closely matches the linear SPT prediction and shows no enhancement of the BAO peak, confirming that the 13-type contribution in the source term is essential for capturing the nonlinear sharpening of the BAO feature.

In summary, a complete description of BAO in the context of ULPT requires two complementary effects:

• •

exponential suppression of the BAO amplitude by the displacement-mapping factor;

• •

mild enhancement of the BAO peak due to the 13-type contribution in the 1-loop source power spectrum correction, $P_{\rm J}^{(13)}$.

IX.3 IR-Resummed Model Derived from ULPT

One of the notable features of ULPT is its ability to naturally reproduce the structure of IR-resummed models, which are based on the decomposition of the linear power spectrum into wiggly and no-wiggly components (see, for example, Refs. [7, 8] for standard treatments). In Ref. [21], the linear IR-resummed model, Eq. (66), was derived directly from the ULPT formulation as the simplest case. In the present work, we generalize this derivation to construct a more comprehensive formulation of IR-resummed models within the ULPT framework. As a concrete example, we derive an IR-resummed expression incorporating one-loop corrections and compare it with the full one-loop prediction obtained directly from ULPT.

Within the ULPT framework, the IR-resummed model is constructed so that the displacement-mapping factor is responsible solely for describing the nonlinear damping of the BAO feature, while the overall shape of the power spectrum is governed by SPT. To proceed, we decompose the source power spectrum into a BAO-carrying wiggle part, $P_{\mathrm{J,w}}$, and a smooth no-wiggle part, $P_{\mathrm{J,nw}}$, as

[TABLE]

The corresponding source correlation function is similarly decomposed as

[TABLE]

Substituting this decomposition into the general expression for the ULPT power spectrum [Eq. (18)] yields

[TABLE]

where $\mu=\hat{r}\cdot\hat{k}$.

To simplify the convolution integrals arising from the displacement-mapping factor, we introduce an approximation. The wiggle part of the source correlation function $\xi_{\mathrm{J,w}}(r)$ is sharply peaked around the BAO scale and vanishes elsewhere, resembling a delta function. Therefore, we approximate the radial dependence of the displacement-mapping factor by fixing it at $r=r_{\mathrm{BAO}}\sim 105\,h^{-1}\mathrm{Mpc}$, and adopt the conventional choice $\mu=1$ for the angular dependence. This yields

[TABLE]

where we used the relation $\int d^{3}r\,e^{-i\bm{k}\cdot\bm{r}}\,\xi_{\mathrm{J,w}}(r)=P_{\mathrm{J,w}}(k)$.

Next, we consider the no-wiggle component. In the spirit of IR-resummed models, we assume that the displacement-mapping factor primarily affects the BAO feature, while the broadband shape of the power spectrum is governed by SPT. Accordingly, the no-wiggle contribution can be expressed as

[TABLE]

which effectively corresponds to a finite-order SPT prediction evaluated using the no-wiggle source correlation function, and does not contain any BAO oscillations.

Combining the wiggle and no-wiggle results given in Eqs. (72) and (73), the general form of the IR-resummed model in the ULPT framework is given by

[TABLE]

where $P_{\mathrm{J,w}}(k)$ captures the BAO oscillatory component modulated by the displacement-mapping factor, and $P_{\mathrm{SPT,nw}}(k)$ corresponds to the broadband shape computed from SPT.

We now derive the explicit form of the IR-resummed model at one-loop order. We begin by evaluating the displacement-mapping exponent in the linear approximation. Using Eq. (30), we obtain

[TABLE]

where

[TABLE]

is the standard BAO damping scale consistent with Eq. (67).

Next, we consider the decomposition of the source power spectrum into wiggle and no-wiggle components. The linear matter power spectrum is split as

[TABLE]

The one-loop 13-type source correction $P_{\mathrm{J}}^{(13)}(k)$ is proportional to $P_{\mathrm{lin}}(k)$, and can thus be decomposed using Eq. (77) as

[TABLE]

As discussed in Sec. IX.2, the BAO contribution from the 22-type source term $P_{\rm J}^{(22)}$ is subdominant and can be neglected in the wiggle part. Therefore, the one-loop wiggle contribution to the source power spectrum is given by

[TABLE]

where we have used the identity $P_{\mathrm{J}}^{(13)}=P^{(13)}+k^{2}\bar{\sigma}^{2}P_{\mathrm{lin}}$ from Eqs. (58) and (61).

On the other hand, the no-wiggle part of the IR-resummed model is obtained by evaluating the SPT prediction up to one-loop order, replacing all occurrences of the linear power spectrum $P_{\mathrm{lin}}$ with the no-wiggle spectrum $P_{\mathrm{nw}}$. The resulting expression reads

[TABLE]

where

[TABLE]

Here, $P^{(13)}$ and $P^{(22)}$ represent the one-loop contributions in SPT, and the notation $[P_{\mathrm{lin}}\to P_{\mathrm{nw}}]$ indicates that these terms are evaluated by replacing the linear power spectrum with its no-wiggle counterpart.

Substituting Eqs. (75), (79), and (80) into the general IR-resummed expression given in Eq. (74), we obtain the final one-loop IR-resummed model:

[TABLE]

Figure 8 presents a comparison between the one-loop ULPT prediction and the one-loop IR-resummed model in configuration space. The two results exhibit excellent agreement, confirming that the IR-resummed model accurately reproduces the nonlinear damping of the BAO feature.

X Comparison with Emulators

We now assess the predictive accuracy and practical applicability of ULPT by comparing its one-loop predictions against simulation-based emulators. In Sec. X.1, we examine the fiducial cosmology, for which the Dark Emulator provides high-precision predictions owing to multiple independent realizations. This setting enables a stringent benchmark test of the ULPT framework in both Fourier and configuration space. In Sec. X.2, we extend the comparison across a broader cosmological parameter space, using 100 randomly sampled models within the overlapping validity range of the Dark Emulator and Euclid Emulator 2. This analysis allows us to evaluate the robustness of ULPT under realistic parameter variations relevant for cosmological inference.

X.1 Fiducial Cosmology

We begin by examining the fiducial cosmology adopted in the Dark Emulator suite, corresponding to the Planck 2015 best-fit $\Lambda$CDM model. This model has been simulated using 14 independent $N$-body realizations, substantially reducing sample variance and yielding reliable emulator predictions for benchmarking theoretical models.

In this subsection, we compare the ULPT prediction for the nonlinear matter power spectrum with the Dark Emulator output, and also present the corresponding two-point correlation function obtained via inverse Hankel transformation.

X.1.1 $z=0$

Figure 10 presents a three-way comparison between the nonlinear matter power spectrum predicted by the Dark Emulator ($P_{\rm Dark}$), the one-loop SPT prediction ($P_{\rm SPT}$), and the one-loop ULPT prediction ($P_{\rm ULPT}$) at redshift $z=0$. In the upper panel, each spectrum is divided by the corresponding no-wiggle power spectrum ($P_{\rm nw}$), while the lower panel shows the fractional differences of the perturbative models relative to $P_{\rm Dark}$.

On large scales ($k\lesssim 0.1\,h\,{\rm Mpc}^{-1}$), SPT and ULPT both reproduce the emulator within approximately three percent, indicating that linear growth and next leading-order corrections are accurately captured where bulk displacements remain small.

At small scales ($0.1\lesssim k\lesssim 0.4\,h\,{\rm Mpc}^{-1}$), the SPT prediction deviates rapidly from the emulator, while ULPT remains accurate to within five percent. This behavior supports the analysis in Sec. VII, where it was shown that the primary source of inaccuracy in SPT arises from the expansion of the displacement-mapping factor, which is avoided in ULPT by retaining its exponential form.

Configuration-space results are shown in Fig. 10, where the two-point correlation functions obtained via inverse Hankel transforms of $P_{\rm ULPT}$ and $P_{\rm Dark}$ are compared. ULPT agrees with the emulator to within five percent down to $r\simeq 10\,h^{-1}{\rm Mpc}$ and reproduces both the position and amplitude of the BAO peak at $r\approx 105\,h^{-1}{\rm Mpc}$ with sub-percent precision. This level of agreement in real space mirrors the consistency observed in Fourier space, thereby confirming the internal coherence of the ULPT framework.

X.1.2 $z=0.5$ and $z=1.0$

Because ongoing galaxy surveys such as DESI and the Prime Focus Spectrograph (PFS) [34] target the redshift range $z\gtrsim 0.5$, Figures 12 and 12 extend the comparison to $z=0.5$ and $z=1.0$. For the power spectrum (Fig. 12), ULPT agrees with the emulator to within three percent at $z=0.5$ and within two percent at $z=1.0$ up to $k\simeq 0.4\,h\,{\rm Mpc}^{-1}$. This progressive improvement with redshift is a natural consequence of the reduced nonlinear gravitational effects at earlier cosmic times.

The correlation functions shown in Fig. 12 exhibit a similarly close agreement: ULPT matches the emulator to better than three percent at $z=0.5$ and two percent at $z=1.0$ for separations $r\gtrsim 10\,h^{-1}\mathrm{Mpc}$. A mild ($\lesssim 0.5\%$) excess appears near $r\approx 50\,h^{-1}\mathrm{Mpc}$, but the overall agreement remains consistent with that seen in the power spectrum comparison.

X.2 Various Cosmologies

We next investigate the performance of ULPT across a diverse set of cosmological models. Specifically, we compare the ULPT predictions for the nonlinear matter power spectrum against those of the Dark Emulator and Euclid Emulator 2, using a set of $100$ cosmologies randomly sampled from the parameter region where both emulators are valid. This test is particularly relevant for applications such as Markov Chain Monte Carlo (MCMC) analyses, where repeated evaluations of the power spectrum across a wide parameter space are required.

Figures 14 and 14 present the fractional deviation, $(P_{\rm ULPT}/P_{\rm emu})-1$, for each of the sampled cosmologies at redshifts $z=0.5$ and $z=1.0$.

Dark Emulator (Fig. 14).

At both redshifts, the magenta curves show the mean difference, while the grey lines represent individual models. The mean agreement remains within $\sim 3\%$ at $z=0.5$ and $2\%$ at $z=1.0$ up to $k\simeq 0.4\,h\,\mathrm{Mpc}^{-1}$, matching the precision already achieved for the fiducial model in Fig. 12. However, the individual Dark Emulator realizations display considerable scatter, with deviations exceeding $5\%$ in some cases, due to sample variance (see also Fig. 5). A systematic $\sim 2\%$ excess of the Dark Emulator prediction relative to ULPT is also visible at very large scales ($k\lesssim 0.05\,h\,\mathrm{Mpc}^{-1}$), and is attributed to the same sample variance effect.

Given that ULPT provides a variance-free theoretical prediction with 2–3% accuracy across the full range of scales, it serves as a more stable and reliable baseline than the individual emulator outputs, especially in parameter or scale regimes where sample variance is non-negligible.

Euclid Emulator 2 (Fig. 14).

Replacing the Dark Emulator with Euclid Emulator 2 removes the large-scale bias: no systematic offset is observed for $k\lesssim 0.05\,h\,\mathrm{Mpc}^{-1}$, and the mean curve closely follows the ULPT prediction. This excellent agreement on large scales is expected, as perturbation theory is most reliable in this regime. It confirms the accuracy of Euclid Emulator 2 and, reciprocally, supports the validity of the ULPT framework. The ULPT predictions shown here are computed using the same pipeline as in Fig. 14, with the linear matter power spectrum consistently taken from the Dark Emulator.

Owing to its paired-and-fixed simulation strategy, Euclid Emulator 2 exhibits a typical variance of $\lesssim 1\%$ across models. After averaging, its predictions remain consistent with those of the Dark Emulator, and the agreement with ULPT again stays within $3\%$ at $z=0.5$ and $2\%$ at $z=1.0$ up to $k\simeq 0.4\,h\,\mathrm{Mpc}^{-1}$.

Implications for parameter inference.

Because modern MCMC analyses typically require $\mathcal{O}(10^{5})$ evaluations of the matter power spectrum, a fast, variance-free, and percent-level accurate theory is highly desirable. The consistency demonstrated across $100$ cosmological models shows that ULPT meets these requirements while remaining computationally competitive with emulator-based approaches.

These findings, summarized in Figs. 14 and 14, demonstrate that ULPT enables accurate and efficient cosmological parameter estimation across a broad parameter space.

XI Future Prospects

The results presented in this work establish one-loop ULPT as a fast and percent-level accurate description of the nonlinear matter power spectrum across the six-dimensional $w$CDM parameter space supported by current $N$-body emulators. Several extensions naturally follow.

First, since ULPT contains no phenomenological fitting parameters, it can be applied beyond the emulator training volume, offering reliable predictions in cosmologies not covered by simulations. The formalism also admits straightforward extensions to redshift-space distortions and density-field reconstruction, as both effects preserve the structural decomposition into Jacobian and displacement-mapping components.

Incorporating galaxy bias is another promising direction. The original ULPT framework [21] allows for a bias expansion that avoids renormalization by including additive Galileon operators in the Jacobian deviation. Extending the current implementation to biased tracers will enable predictions for halo clustering statistics that can be directly compared to emulator-based models.

Finally, extending ULPT to two loops is expected to further improve accuracy toward the sub-percent level. Since the residual $2$–$3\%$ discrepancies observed at one loop match the typical size of two-loop corrections, such an extension will provide a crucial test of the framework’s robustness and its ability to surpass the limitations of standard perturbation theory.

XII Conclusion

In this work, we have presented a fast and accurate formulation for computing the nonlinear matter power spectrum at one-loop order within the framework of Unified Lagrangian Perturbation Theory (ULPT). The key feature of ULPT is the structural decomposition of the density field into two physically distinct components: the Jacobian deviation, which describes intrinsic nonlinear growth, and the displacement-mapping effect, which captures large-scale convective distortions.

In this framework, the power spectrum naturally separates into two contributions: a source term that encodes correlations between the displacement field and the Jacobian deviation, and a displacement-mapping factor that encapsulates the uncorrelated impact of large-scale flows. This structural separation not only facilitates physical interpretation but also ensures infrared (IR) safety by construction, playing a central role in the resummation of long-wavelength modes.

We have demonstrated that the breakdown of standard perturbation theory (SPT) at small scales can be traced to its expansion of the displacement-mapping factor, which inherently possesses an exponential form. In contrast, ULPT retains this structure without expansion, leading to improved convergence and a more physically faithful description of nonlinear evolution. Notably, ULPT and SPT share the same one-loop source term $P_{\rm J}^{(\text{1-loop})}$, but differ in their treatment of the displacement-mapping contribution: SPT expands it, while ULPT resums it. This distinction accounts for the superior performance of ULPT at nonlinear scales.

We have implemented ULPT numerically using a combination of FFTLog and FAST-PT techniques, achieving rapid evaluations suitable for cosmological applications; the full one-loop power spectrum can be computed in approximately 2 seconds on a standard laptop. The predictive accuracy of ULPT has been validated against both the Dark Emulator and Euclid Emulator 2 across 100 cosmological models. In particular, for redshifts $z\geq 0.5$, ULPT predictions match simulation-based emulators at the 2–3% level up to $k\simeq 0.4\,h\,\mathrm{Mpc}^{-1}$, without invoking any nuisance parameters. Comparable accuracy is also observed in configuration space, where the two-point correlation function remains consistent with emulator results down to $r\simeq 10\,h^{-1}\mathrm{Mpc}$.

Our analysis also elucidates the mechanism responsible for nonlinear baryon acoustic oscillation (BAO) damping. In ULPT, the displacement-mapping factor induces an exponential suppression of the BAO amplitude by smearing out oscillatory modes, while the nonlinear source term, particularly the 13-type component, mildly sharpens the peak through constructive growth. This complementary interplay reproduces the characteristic BAO features observed in simulations. In contrast, SPT distorts the BAO signal because its expansion of the displacement-mapping factor fails to capture the correct nonlinear structure. This again highlights the importance of preserving the full exponential form, not only for the broadband shape of the power spectrum, but also for accurately describing the nonlinear evolution of the BAO peak.

Looking ahead, ULPT provides a solid theoretical foundation for several future developments. These include consistent extensions to redshift-space distortions, density-field reconstruction, and galaxy bias, all of which can be seamlessly incorporated into the existing formalism. Since ULPT contains no phenomenological fitting parameters, it can also be applied beyond the emulator training volume, offering reliable predictions in cosmologies not covered by simulations. Furthermore, extending ULPT to two-loop order is a natural next step toward achieving sub-percent theoretical precision and may offer new insights into the limitations of standard perturbative approaches.

Acknowledgements.

N.S. acknowledges financial support from JSPS KAKENHI Grant No. 25K07343, administratively hosted by the National Astronomical Observatory of Japan. The author is grateful to the developers of the publicly available Python packages used in this work, including mcfit [26], which was used for Hankel transforms, and FAST-PT [27, 28], which enabled efficient evaluation of one-loop corrections. Nonlinear reference spectra were obtained from the Dark Emulator [22] and the Euclid Emulator 2 [23]. N.S. also acknowledges the use of ChatGPT (OpenAI) for assistance in language refinement, literature exploration, and Python scripting during the preparation of this manuscript. The author is thankful to the referee for their constructive comments, which have significantly improved the quality and clarity of the paper.