PySCo-EFT: Python Toolkit for EFT Analyses

• PySCo-EFT is a Python suite of computational tools that integrates EFT amplitude bounds with cosmological Dark Energy reconstruction and N-body simulation techniques.
• It employs dispersive sum-rule methods combined with data-driven inversion of H(z) measurements to enforce causality and reconstruct key EFT functions.
• The toolkit features a high-performance PM N-body solver using advanced multigrid and finite-difference methods to accurately simulate nonlinear large-scale structure dynamics.

PySCo-EFT is a suite of Python-based computational tools and N-body simulation codes designed for Effective Field Theory (EFT) analyses in cosmology and quantum field theory. It encompasses both pipelines for reconstructing the EFT of Dark Energy (EFTofDE) from observational data and for establishing rigorous amplitude bounds in scalar EFTs, as well as high-performance numerical solvers for nonlinear large-scale structure simulations. The PySCo-EFT project integrates key developments from dispersive amplitude bounds, cosmological data inversion methods, and particle-mesh numerics, validated by direct comparison with theory and benchmark codes.

1. Conceptual Foundations and Scope

PySCo-EFT ("Python Scalar Code for Effective Field Theory") provides a modular framework for both theoretical and simulation-driven EFT studies. The package focuses on two principal domains:

• Dispersive Amplitude Bounds: Implements the sum-rule and positivity-based approach to bounding EFT coefficients of scalar field theories under causality and unitarity, as developed in "Extremal Effective Field Theories" (Caron-Huot et al., 2020).
• Cosmological EFT of Dark Energy: Supports the data-driven reconstruction and numerical simulation of the EFTofDE, including dark energy background reconstruction from $H(z)$ measurements (Okamatsu et al., 2 Dec 2025) and large-scale structure evolution with nonlinear screening (Ganjoo et al., 16 Apr 2026).

The architecture is characterized by modular Python/Numba codebases, statistical interfaces for covariance propagation, and, in its simulation variant, a parallelized particle-mesh (PM) N-body solver for scalar-tensor theories with nonlinear interactions and screening.

2. Mathematical Formalism of the Upstream Modules

2.1 Dispersive Sum-Rule and SDP Pipeline

The amplitude-bounds module of PySCo-EFT sets up the forward-scattering amplitude for massless scalar $2\to2$ processes using a twice-subtracted dispersion relation: $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$ Positivity ($\rho(s) \ge 0$) gives rise to a Hankel-type matrix constraint: $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$ alongside crossing-symmetry (“null”) constraints that are imposed as linear equalities over the $a_{2k}$. The resulting problem is a semidefinite program (SDP) that can be solved for optimal upper and lower bounds on any low-energy Wilson coefficient. The practical code structure employs CVXPY/CVXOPT for SDP formulation and solution, and SciPy’s ConvexHull for outlining the extremal, allowed region (Caron-Huot et al., 2020). Extremal points correspond to analytic simple spectra (such as $\delta$-peaks and stu-pole models).

2.2 EFTofDE Background and Phenomenological Functions

Within the EFTofDE formalism, the background dynamics in unitary-gauge scalar-tensor theories are controlled by time-dependent functions $M^2(t)$ (Planck mass squared), $\Lambda(t)$ (effective cosmological constant), and $c(t)$ (scalar-kinetic coefficient). These tadpole-like EFT functions are related to the Friedmann equations as: $2\to2$0

$2\to2$1

A common reparametrization uses the $2\to2$2-basis:

• $2\to2$3 (Planck mass run-rate)
• $2\to2$4 (kineticity)
• $2\to2$5, $2\to2$6 (braiding, tensor speed excess)

Reconstruction of $2\to2$7, $2\to2$8, $2\to2$9 from $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$0 and $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$1—with observed $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$2 obtained via cosmic chronometers and Gaussian process regression—uses inversion formulae (e.g., eqs. (3.9)–(3.12) in (Okamatsu et al., 2 Dec 2025)) with all time derivatives mapped to $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$3-derivatives. These routines are implemented in Python (NumPy/SciPy), following a strict module separation between data acquisition, regression, differentiation, and nonlinear inversion.

3. PySCo-EFT N-Body Simulation Code for EFTofDE

3.1 Numerical Core and Data Structures

The PySCo-EFT PM code extends the original PySCo architecture with modules for EFTofDE perturbations and background evolution (Ganjoo et al., 16 Apr 2026). Key features include:

• Particle positions and velocities stored in Numba-typed float64 arrays for JIT compilation and parallelization.
• Density assignment to a uniform grid using cloud-in-cell (CIC) interpolation.
• Effective grid-level fields for gravitational potential and the scalar degree of freedom.
• Parallelization is realized with thread-level splitting of cell/particle loops, no MPI.
• All EFTofDE coefficient functions are precomputed on a grid in $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$4 and interpolated per-step.

3.2 Physical Equations and Algorithmic Methods

The EFTofDE PM code solves:

• The modified Poisson equation for the gravitational potential $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$5,

$$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$6

where $$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$7.

• A nonlinear scalar field equation incorporating Vainshtein screening,

$$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$8

These equations are solved using a full-approximation-storage (FAS) multigrid method with local quadratic updates and under-relaxed Jacobi smoothing. The nonlinear Vainshtein bracket is discretized using finite-difference stencils for all second derivatives to maintain accuracy in high-density regimes.

4. Validation, Benchmarking, and Numerics

Validation includes:

• Agreement with linear theory on large scales ($$A(s) = \sum_{k=0}^\infty a_{2k} \Lambda^{2k} s^k,\qquad a_{2k} = \frac{1}{\pi} \int_{M^2}^\infty \frac{ds'}{s'^{k+1}} \Im A(s').$$9) with power spectrum ratios $\rho(s) \ge 0$0 agreeing to <0.4% between code and analytic growth prediction.
• Cross-validation with ECOSMOG-EFT (RAMSES-AMR) finds agreement in $\rho(s) \ge 0$1 to <1% up to $\rho(s) \ge 0$2 for both linearised and fully nonlinear runs with $\rho(s) \ge 0$3 particles in $\rho(s) \ge 0$4Mpc boxes.

Extensive resolution and parameter convergence tests quantify the effects of mass resolution, finite-volume, refinement thresholds, smoothing parameters, and starting redshift on $\rho(s) \ge 0$5. For example, at $\rho(s) \ge 0$6, the aggregate numerical error is kept below 2% for all tested settings.

5. EFTofDE Parameter Dependence and Screening

The effect of model parameters is systematically explored:

• Varying $\rho(s) \ge 0$7 determines the large-scale power spectrum boost and the importance of the Vainshtein mechanism. For $\rho(s) \ge 0$8, the larger the negative braiding, the stronger the linear boost.
• With $\rho(s) \ge 0$9, the effective $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$0 is enhanced at late times, increasing small-scale clustering; for $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$1, small-scale growth is suppressed.
• Screening is quantified by $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$2. For $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$3, linearized and full nonlinear equations agree to <1% up to $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$4. For $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$5, fully nonlinear solution is essential; the difference reaches up to 35% in $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$6 at $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$7 for the largest $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$8 tested (Ganjoo et al., 16 Apr 2026).

6. Data-Driven EFTofDE Reconstruction Pipeline

The background module provides a pipeline:

1. Loading $$M_{rs} = \sum_{k=0}^{K-1} a_{2k} \mu_k^{(r,s)},\qquad \mu_k^{(r,s)} = \int_1^\infty dx\, x^{-(k+r+s+1)} > 0,\qquad M \succeq 0,$$9 data from cosmic chronometers.
2. Gaussian process regression for interpolation and error propagation of $a_{2k}$0 and $a_{2k}$1.
3. Inversion of $a_{2k}$2 to obtain $a_{2k}$3 (either constant, or using a parameterized $a_{2k}$4), then $a_{2k}$5 and $a_{2k}$6 using analytic formulae, propagating uncertainties throughout.
4. Validation is performed by comparing residuals $a_{2k}$7 and normalized $a_{2k}$8 to $a_{2k}$9CDM predictions and computing $\delta$0-type tests.
5. For quintessence, reconstruction proceeds via the mapping $\delta$1 with diagnostics on the scalar field’s reality and potential slope (Okamatsu et al., 2 Dec 2025).

7. Algorithmic Summary and Usage

Representative module functions, all verifiable in the cited materials, include:

• Dispersive amplitude bounds: PySCoEFT class with solve_bound() to maximize/minimize $\delta$2 given positivity/null constraints (Caron-Huot et al., 2020).
• EFTofDE background: reconstruct_EFT_functions(z,H,dHdz,...) for inversion, compute_dM2dz(z, M2, alpha_M) for derivatives.
• Particle mesh code: command-line specification of cosmological and EFTofDE parameters, grid, particle number, redshift, and Jacobi smoothing; $\delta$3 cpu-hr per 256$\delta$4 run with percent-level accuracy at $\delta$5 (Ganjoo et al., 16 Apr 2026).

The pipeline is encapsulated in self-contained Python modules—enabling researchers to perform causality/analyticity-constrained amplitude bounding, model-independent background reconstruction from cosmological data, and efficient PM-based N-body simulation of EFTofDE models, including nonlinear scalar screening and observable effects on the matter power spectrum.