Unified Lagrangian Perturbation Theory (ULPT)

• ULPT is a Lagrangian perturbation framework that separates intrinsic nonlinear density evolution from convective displacement mapping, clearly isolating infrared (IR) effects.
• The approach enables nonperturbative resummation of long-wavelength displacements, yielding IR-safe expressions for power spectra, bispectra, and BAO damping with efficient numerical implementations.
• ULPT unifies treatments of biased tracers, redshift-space distortions, and reconstruction, extending its applicability to modified gravity and advanced large-scale structure analyses.

Unified Lagrangian Perturbation Theory (ULPT) is a Lagrangian formulation of cosmological perturbation theory in which nonlinear density fluctuations are reorganized into an intrinsic source sector and a displacement-mapping sector. In the formulation developed for large-scale structure correlators, the intrinsic source is the Jacobian deviation of the Lagrangian-to-Eulerian map, while the displacement-mapping factor captures the convective remapping induced by bulk flows; this separation makes infrared (IR) physics explicit, enables nonperturbative resummation of long-wavelength displacements, and yields IR-safe expressions for power spectra, bispectra, and related observables. Across the recent literature, the term also denotes a broader unification in which one universal spatial operator basis is used across cosmological models, with model dependence isolated in time-dependent coefficients. In that combined sense, ULPT is both a field-level reorganization of LPT and a unifying perturbative language for real space, redshift space, reconstruction, bias, modified gravity, and higher-order correlators (Sugiyama, 29 Aug 2025, Sugiyama, 24 Aug 2025, Marinucci et al., 2024, Sugiyama, 24 Feb 2026).

1. Conceptual scope and nomenclature

In the recent large-scale-structure literature, ULPT denotes two closely related unifications. The first is structural: the matter or tracer density is written as a product of a source field and a displacement-mapping factor, so that intrinsic nonlinear growth and convective transport are treated as distinct physical ingredients. The second is operator-theoretic: the displacement field is expanded in a universal basis of spatial operators fixed by symmetry and conservation laws, while cosmology, gravity model, or multi-species effects enter only through time-dependent coefficients determined by ordinary differential equations.

The structural formulation emphasizes that the density field should be defined in Eulerian space and only then rewritten in Lagrangian variables so that the Jacobian structure remains explicit. The model-independent formulation emphasizes that rotational and translational invariance, extended Galilean invariance, the equivalence principle, locality in the single-stream regime, mass conservation, and momentum conservation are sufficient to determine the allowed operator content of the displacement. These are not competing definitions. Rather, they organize different aspects of the same Lagrangian program: one isolates the physical origin of IR effects in correlators, while the other isolates the universal spatial content of perturbation theory across theories beyond $\Lambda$CDM (Sugiyama, 24 Aug 2025, Marinucci et al., 2024).

2. Field-level formulation and IR-safe decomposition

The basic map is

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

with Jacobian

$$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$

Mass conservation in the single-stream regime implies

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

ULPT introduces the Jacobian deviation

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$

and rewrites the Fourier-space density as

$$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$

This representation makes the two sectors explicit. The factor $\delta_J$ encodes intrinsic nonlinear growth. The exponential $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$ maps that source from Lagrangian to Eulerian space and contains the effect of long-wavelength bulk flows. The power spectrum can then be organized as

$$P(k)=e^{-\overline{\Sigma}(k)}\int d^3r\,e^{-i\mathbf{k}\cdot\mathbf{r}}\,e^{\Sigma(k,r,\mu)}\,\xi_J(r),$$

with $\xi_J(r)$ the source correlation function, $$\mathbf{x}=\mathbf{q}+\boldsymbol{\Psi}(\mathbf{q}),$$0, and

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

A useful decomposition is

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

where $$\mathbf{x}=\mathbf{q}+\boldsymbol{\Psi}(\mathbf{q}),$$3 is the convolution-free part and $$\mathbf{x}=\mathbf{q}+\boldsymbol{\Psi}(\mathbf{q}),$$4 is the convolution-containing part.

The IR-safe character of ULPT follows from keeping the displacement-mapping factor in its full exponential form rather than expanding it as a series. In the Gaussian limit,

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

with

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

More generally, the cumulant-resummed exponential absorbs the long-wavelength displacement sector nonperturbatively. Standard perturbation theory fails at small scales because it expands this exponential and thereby inherits spurious IR sensitivity and poor convergence. ULPT removes that failure mode by construction (Sugiyama, 29 Aug 2025, Sugiyama, 24 Aug 2025).

3. One-loop power spectrum, BAO physics, and numerical realization

At one loop in real space, ULPT organizes the nonlinear matter spectrum as

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

where the source sector is

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

and the displacement-mapping correction is generated by the linear part of the displacement cumulant acting on the linear source correlation. In this formulation,

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

so the source term is explicitly regularized against the IR-sensitive contribution that would otherwise appear in the standard $$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$0-type term. The final one-loop expression is

$$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$1

The same split clarifies the baryon acoustic oscillation mechanism. For the wiggle component,

$$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$2

with

$$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$3

or equivalently

$$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$4

In ULPT, BAO damping is the exponential suppression produced by the displacement-mapping factor, while the source sector, especially the $$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$5-type contribution proportional to $$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$6, produces mild peak sharpening in configuration space. The observed BAO pattern is therefore the combined effect of displacement-induced smearing and intrinsic nonlinear source evolution.

A fast numerical implementation follows directly from this factorization. The published algorithm uses FFTLog for the Hankel transforms, FAST-PT for the one-loop SPT building blocks, and a truncated expansion of the angular dependence in the convolution-containing term. Truncating that expansion at $$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$7–$$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$8 is sufficient for sub-percent convergence up to $$J(\mathbf{q})=\det\!\left(\delta_{ij}+\Psi_{i,j}(\mathbf{q})\right).$$9. Typical runtimes are $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$0 seconds on a standard laptop for 450 $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$1-points spanning $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$2 to $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$3. Across 100 sampled cosmologies, the one-loop ULPT matter spectrum agrees with the Dark Emulator and Euclid Emulator 2 at the $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$4–$$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$5 level up to $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$6 for $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$7, without nuisance parameters; at $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$8, the agreement is $$1+\delta(\mathbf{x}(\mathbf{q}))=\frac{1}{J(\mathbf{q})}.$$9 up to the same scale. In configuration space, the two-point correlation function remains accurate down to $\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$0. The corresponding open-source implementation is released as ulptkit (Sugiyama, 29 Aug 2025).

4. Bias, redshift-space distortions, and reconstruction

A central extension of ULPT is that the same source–mapping split applies to biased tracers, redshift space, and post-reconstruction fields. In redshift space, the line-of-sight map is

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$1

with

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$2

The redshift-space tracer density is then written as

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$3

At linear order, this reproduces the Kaiser limit,

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$4

Reconstruction is treated as a pure displacement update. If $\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$5 is the reconstruction shift, then

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$6

and analogously in redshift space

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$7

The source sector remains unchanged: reconstruction modifies only the displacement-mapping factor. This is the basis for the single-Gaussian BAO damping pre- and post-reconstruction, and for the residual Gaussian suppression in pre/post cross spectra,

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$8

For galaxy bias, ULPT adopts the same Galileon operator basis that appears in the intrinsic dark-matter source. The key second- and third-order invariants are

$\delta_J(\mathbf{q})\equiv 1-J(\mathbf{q}),$9

$$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$0

and

$$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$1

Within ULPT,

$$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$2

Because these Galileon operators have vanishing ensemble and volume averages, the bias expansion is renormalization-free at the field level. At second order, ULPT predicts

$$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$3

In the one-loop biased-tracer implementation, three parameters suffice for correlation functions in configuration space and four parameters are needed for power spectra in Fourier space. Joint fits to halo auto- and cross-spectra from the Dark Emulator, covering nine redshift-mass combinations with 100 cosmologies each, show that a single set of bias parameters reproduces both spectra within $$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$4 up to $$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$5 for $$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$6–2, and to $$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$7 for $$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$8; the same parameters match two-point correlation functions down to $$\widetilde{\delta}(\mathbf{k})=\int d^3q\,e^{-i\mathbf{k}\cdot\mathbf{q}}\,e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}(\mathbf{q})}\,\delta_J(\mathbf{q}).$$9 (Sugiyama, 24 Aug 2025, Sugiyama, 4 Sep 2025).

5. Bispectrum, $\delta_J$0-point structure, and operator-basis unification

ULPT extends directly to the dark-matter bispectrum by promoting the source–mapping decomposition from two-point to three-point statistics. For bispectrum components labeled by a pair of independent wavevectors, the general ULPT form is

$\delta_J$1

At one loop, the exponent is kept nonperturbatively at linear order while the source bispectrum is computed to tree plus one loop. In the auto-bispectrum, the coincident-point identity $\delta_J$2 implies exact nonperturbative IR cancellation. In mixed pre/post reconstruction bispectra, that identity is violated, and a universal residual Gaussian factor appears instead: for example,

$\delta_J$3

ULPT also supplies an analytic description of BAO damping in the bispectrum. The mixed wiggle–smooth pieces acquire single-leg Gaussian factors such as $\delta_J$4 or $\delta_J$5, while the wiggle–wiggle sector carries a triangle-dependent structure,

$\delta_J$6

where $\delta_J$7 contains the additional BAO-scale modulation associated with two correlated wiggle-bearing separations. In $\delta_J$8CDM that modulation is small, so familiar Gaussian-damped templates are recovered as accurate approximations.

A broader sense of ULPT appears in the bootstrap formulation of the displacement field. There the displacement is written as

$\delta_J$9

with scalar operators $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$0 and transverse vector operators $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$1 generated recursively from lower-order Hessians and $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$2-tensor contractions. The basis is universal for theories sharing the symmetries of $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$3CDM, including Horndeski-type modified gravity and multiple non-relativistic species. The explicit enumeration reaches sixth order, with five independent scalars and three vectors at $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$4, fifteen scalars and eleven vectors at $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$5, and fifty-two scalars and thirty-five vectors at $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$6. This universal operator basis is conceptually aligned with the correlator-level ULPT decomposition: in both cases, the unification lies in separating universal spatial structure from theory-dependent or observable-dependent evolution (Sugiyama, 24 Feb 2026, Chen et al., 2024, Marinucci et al., 2024).

6. Antecedents, generalizations, and limitations

ULPT emerged from earlier developments in Lagrangian resummation and higher-order LPT. Convolution Lagrangian perturbation theory exponentiated the full $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$7 dependence of the Gaussian displacement kernel, recovered the Zel’dovich approximation as the lowest order for matter, and improved on earlier resummation schemes for real-space and redshift-space correlation functions. Independent work on recursion relations showed how LPT displacement kernels can be generated iteratively from the standard SPT kernels, while maintaining the exact Jacobian relation between displacement and density. Fourth-order LPT established that even with a purely longitudinal first-order perturbation, transverse contributions appear at third order and feed into six longitudinal-mode equations and four transverse-mode equations at fourth order. These higher-order results are relevant to ULPT because one-loop bispectrum calculations are closed by 4LPT, while still higher-loop resummations require correspondingly higher-order displacement information.

ULPT has also been generalized away from the simplest $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$8CDM dust setting. In modified gravity with scale-dependent gravitational strength, LPT acquires an effective $e^{-i\mathbf{k}\cdot\boldsymbol{\Psi}}$9 and frame-lagging terms generated when Eulerian Poisson and Klein–Gordon operators are rewritten in Lagrangian coordinates. Those frame-lagging terms are essential: without them, the theory does not recover the $$P(k)=e^{-\overline{\Sigma}(k)}\int d^3r\,e^{-i\mathbf{k}\cdot\mathbf{r}}\,e^{\Sigma(k,r,\mu)}\,\xi_J(r),$$0CDM limit at large scales and IR cancellations fail. A different generalization recasts LPT in relativistic cosmology using Cartan coframes, where the gravitoelectric subsystem reproduces the Newtonian hierarchy under the Minkowski restriction and provides an all-order relativistic extension of the Lagrangian hierarchy.

The current domain of validity remains the quasi-linear, single-stream regime. For the one-loop matter spectrum, the published ULPT implementation begins to deviate by a few percent at $$P(k)=e^{-\overline{\Sigma}(k)}\int d^3r\,e^{-i\mathbf{k}\cdot\mathbf{r}}\,e^{\Sigma(k,r,\mu)}\,\xi_J(r),$$1 and/or $$P(k)=e^{-\overline{\Sigma}(k)}\int d^3r\,e^{-i\mathbf{k}\cdot\mathbf{r}}\,e^{\Sigma(k,r,\mu)}\,\xi_J(r),$$2, where two-loop corrections become important. The same limitation appears in biased-tracer applications, where the usable $$P(k)=e^{-\overline{\Sigma}(k)}\int d^3r\,e^{-i\mathbf{k}\cdot\mathbf{r}}\,e^{\Sigma(k,r,\mu)}\,\xi_J(r),$$3-range narrows as the linear bias increases. In redshift space, higher multipoles remain more sensitive to small-scale velocity effects and model details than monopole-level statistics. More generally, shell crossing, baryonic physics, nonlocal-in-time effects, and beyond-quasi-static relativistic or modified-gravity corrections lie outside the one-loop formulations presently validated. Even so, the common conclusion of these developments is that ULPT supplies a single perturbative language in which bulk-flow resummation, Jacobian nonlinearities, bias operators, redshift-space mappings, reconstruction, and theory extensions can be organized within the same Lagrangian framework (Carlson et al., 2012, Tatekawa, 2012, Aviles et al., 2017, Alles et al., 2015).