Infrared Nonlocal Cosmological Models

Infrared nonlocal cosmology turns past curvature into an effective source for the present metric. A local auxiliary-field system is computationally useful only when its data reproduce the selected inverse operator. This page derives that restriction for the representative $R\Box^{-2}R$ model and supplies a reproducible benchmark: local evolution is compared with a direct retarded convolution on a matter-era test history.

Required background. PDEs and Green Functions for Gravitation supplies retarded kernels and direct convolution; Gravitational Actions and Boundary Terms supplies the constrained localization and its boundary data.

Helpful background. Einstein Equations supplies FLRW curvature and the Friedmann equations; Formulations of General Relativity supplies background and perturbative initial-data counting; Gravity Beyond General Relativity supplies the branch and validity checks needed before interpreting acceleration.

The RR background system

The following is a formal effective functional: vary first and impose the retarded prescription on the resulting equations. Ordinary variation of a single-field action with a retarded kernel would symmetrize that kernel, as shown on the causality page. Use the reduced Planck mass $M_{\rm Pl}^{-2}=8\pi G$ and write

$$\Gamma_{\mathrm{RR}} =\frac{M_{\mathrm{Pl}}^2}{2}\int d^4x\sqrt{-g} \left[ R-\frac{m^2}{6}R\Box_{\mathrm{ret}}^{-2}R \right],$$

with

$$U=-\Box_{\mathrm{ret}}^{-1}R, \qquad S=-\Box_{\mathrm{ret}}^{-1}U.$$

The localized metric equation is

$$G_{ab}=\frac{m^2}{6}K_{ab}+8\pi G T_{ab},$$

where

$$\begin{aligned} K_{ab}={}&2SG_{ab}-2\nabla_a\nabla_bS+2g_{ab}\Box S\\ &+g_{ab}\nabla_cS\nabla^cU-\frac12g_{ab}U^2\\ &-\left(\nabla_aS\nabla_bU+\nabla_bS\nabla_aU\right). \end{aligned}$$

The auxiliary equations are

$$\Box U=-R, \qquad \Box S=-U.$$

The fields have dimensions $[U]=1$ and $[S]={\rm length}^2$. A useful independent consistency check is the divergence:

$$\nabla^aK_{ab} =-(R+\Box U)\nabla_bS-(U+\Box S)\nabla_bU=0.$$

The cancellation uses the contracted Bianchi identity and the commutator of covariant derivatives acting on $\nabla_bS$. It fails if an auxiliary sign is changed while the metric equation is left unchanged.

For homogeneous fields and signature $(-,+,+,+)$,

$$\ddot U+3H\dot U=R, \qquad \ddot S+3H\dot S=U.$$

The mass $m$ controls when $K_{ab}$ becomes important. Define $8\pi G T^{\mathrm{NL}}_{ab}=m^2K_{ab}/6$. For a spatially flat FLRW solution, acceleration requires the total combination $\rho_m+3p_m+\rho_{\mathrm{NL}}+3p_{\mathrm{NL}}\lt0$. A negative nonlocal contribution must be large enough to overcome the matter contribution; its sign alone is insufficient. That inequality cannot be inferred from the sign of the action term alone. The full background equations and their self-accelerating solutions are derived in Maggiore and Mancarella 2014, §§ II and VI, Eqs. (2)–(7) and (30)–(35).

For example, the $00$ equation gives the actual background constraint:

$$3H^2=8\pi G\rho_m+\frac{m^2}{6} \left(6SH^2+6H\dot S-\dot S\dot U+\frac12U^2\right).$$

This explains why prescribing $a(t)$ tests the auxiliary solver but does not solve the self-consistent cosmology: the same $U$ and $S$ must also satisfy this Friedmann constraint together with matter conservation.

Retarded data are branch data

Let $t_*$ lie in a radiation era and define both inverse operators with zero data:

$$U(t_*)=\dot U(t_*)=S(t_*)=\dot S(t_*)=0.$$

Because $R=0$ in an exactly flat, radiation-dominated FLRW phase, the fields remain zero there. Once the matter trace turns on, the unique retarded solutions are

$$U(t)=\int_{t_*}^{t}\frac{dt'}{a^3(t')} \int_{t_*}^{t'}dt''\,a^3(t'')R(t''),$$ $$S(t)=\int_{t_*}^{t}\frac{dt'}{a^3(t')} \int_{t_*}^{t'}dt''\,a^3(t'')U(t'').$$

These nested integrals are the direct-convolution benchmark. Adding arbitrary solutions of $\ddot U+3H\dot U=0$ or $\ddot S+3H\dot S=0$ changes the inverse prescription and can change the cosmological branch. The background evolution and radiation-to-matter onset are discussed in Belgacem et al. 2018, §§ 3.1–3.2, equations (3.1)–(3.22).

Controlled matter-era benchmark

To test the implementation without mixing in a nonlinear Friedmann solve, fix dimensionless time $1\le t\le4$ and the exact matter-era history

as follows. If $\tau$ is physical time and $\tau_0$ the start of the benchmark, then $t=\tau/\tau_0$, the displayed $H$ and $R$ denote $\tau_0H_{\rm phys}$ and $\tau_0^2R_{\rm phys}$, and the numerical $S$ denotes $S_{\rm phys}/\tau_0^2$. Normalize $a(1)=1$ and impose all four auxiliary data to vanish at $t=1$. The prescribed history is

$$a(t)=t^{2/3}, \qquad H=\frac{2}{3t}, \qquad R=6(\dot H+2H^2)=\frac{4}{3t^2}.$$

The first auxiliary field has the analytic solution

$$U_{\mathrm{exact}}(t) =\frac43\left(\ln t+\frac1t-1\right).$$

The second auxiliary also admits an exact independent answer:

$$S_{\rm exact}(t) =\frac{2t^2}{9}\ln t-\frac{11t^2}{27} +\frac{2t}{3}+\frac{2}{27t}-\frac13.$$

Both $S_{\rm exact}(1)$ and $S'_{\rm exact}(1)$ vanish. Substitution into $S''+2S'/t=U$ verifies the signs without reusing the numerical integrator.

We integrate $(U,\dot U,S,\dot S)$ with fixed-step fourth-order Runge–Kutta. An independent nested SciPy quadrature evaluates both retarded convolutions at $t=4$. The equation residual is computed from centered finite differences rather than from the right-hand side used by the integrator.

The retained reference environment is CPython 3.9.6, NumPy 2.0.2, and SciPy 1.13.1, using IEEE-754 binary64 arithmetic. Save the following executable source exactly as rr_background.py. This program produces the seven-column refinement table printed below. The final command exits nonzero if the pinned versions differ, endpoint errors fail to converge, the finest errors exceed tolerance, the residual is too large, or retarded initial data drift.

import argparse
import math
import sys

import numpy as np
import scipy
from scipy.integrate import quad

PINNED = ("3.9.6", "2.0.2", "1.13.1")


def curvature(t):
    return 4.0 / (3.0 * t * t)


def u_exact(t):
    return (4.0 / 3.0) * (math.log(t) + 1.0 / t - 1.0)


def rhs(t, y):
    u, v, s, w = y
    return np.array([v, curvature(t) - 2.0 * v / t,
                     w, u - 2.0 * w / t])


def rk4(n, t0=1.0, t1=4.0):
    h = (t1 - t0) / n
    ts = np.linspace(t0, t1, n + 1)
    ys = np.zeros((n + 1, 4))
    for i in range(n):
        t = ts[i]
        y = ys[i]
        k1 = rhs(t, y)
        k2 = rhs(t + h / 2.0, y + h * k1 / 2.0)
        k3 = rhs(t + h / 2.0, y + h * k2 / 2.0)
        k4 = rhs(t + h, y + h * k3)
        ys[i + 1] = y + h * (k1 + 2.0 * k2
                              + 2.0 * k3 + k4) / 6.0
    udd = (ys[2:, 0] - 2.0 * ys[1:-1, 0]
           + ys[:-2, 0]) / h**2
    ud = (ys[2:, 0] - ys[:-2, 0]) / (2.0 * h)
    residual = np.max(np.abs(
        udd + 2.0 * ud / ts[1:-1] - curvature(ts[1:-1])
    ))
    return ts, ys, residual


def direct_convolution(t1=4.0):
    def inner_u(tp):
        return quad(lambda s: s * s * curvature(s), 1.0, tp,
                    epsabs=1e-13, epsrel=1e-13)[0]

    u = quad(lambda tp: inner_u(tp) / (tp * tp), 1.0, t1,
             epsabs=1e-13, epsrel=1e-13)[0]

    def inner_s(tp):
        return quad(lambda s: s * s * u_exact(s), 1.0, tp,
                    epsabs=1e-13, epsrel=1e-13)[0]

    s = quad(lambda tp: inner_s(tp) / (tp * tp), 1.0, t1,
             epsabs=1e-13, epsrel=1e-13)[0]
    return u, s


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--check", action="store_true")
    args = parser.parse_args()
    actual = (sys.version.split()[0], np.__version__,
              scipy.__version__)
    if actual != PINNED:
        raise SystemExit(
            f"environment mismatch: expected {PINNED}, got {actual}"
        )
    u_conv, s_conv = direct_convolution()
    rows = []
    for n in (32, 64, 128, 256):
        ts, ys, residual = rk4(n)
        u_error = abs(ys[-1, 0] - u_conv)
        s_error = abs(ys[-1, 2] - s_conv)
        boundary = np.max(np.abs(ys[0]))
        rows.append((n, ys[-1, 0], u_error, ys[-1, 2],
                     s_error, residual, boundary))
        print(n, *(f"{x:.12e}" for x in rows[-1][1:]))
    finest = rows[-1]
    convergent = (
        finest[2] < 1e-8
        and finest[4] < 1e-8
        and finest[5] < 2e-4
        and finest[6] == 0.0
        and rows[-1][2] < rows[-2][2] < rows[-3][2]
        and rows[-1][4] < rows[-2][4] < rows[-3][4]
    )
    if args.check and not convergent:
        raise SystemExit("convergence or benchmark check failed")


if __name__ == "__main__":
    main()

Run the saved source in a clean shell:

python3 --version
python3 -c "import numpy, scipy; print(numpy.__version__, scipy.__version__)"
python3 rr_background.py --check

Convergence and independent comparison

The columns $\Delta U$ and $\Delta S$ compare the local endpoint values with the direct nested convolution. The residual is $\max\lvert\ddot U+2\dot U/t-R\rvert$ over interior nodes.

Steps $N$	$U(4)$	$U$ error	$S(4)$	$S$ error	Maximum residual	Initial mismatch
32	0.848393712265	1.231e-6	0.762382852323	2.902e-6	1.309e-3	0
64	0.848392554736	7.324e-8	0.762380125670	1.750e-7	4.021e-4	0
128	0.848392485955	4.462e-9	0.762379961372	1.072e-8	1.110e-4	0
256	0.848392481769	2.753e-10	0.762379951312	6.632e-10	2.911e-5	0

Halving the step reduces the endpoint errors by approximately sixteen, as expected for RK4. The independently formed finite-difference residual is only second order and decreases by approximately four. The different convergence orders are a useful guard against a self-comparison that reuses the integrator formula. The analytic value

$$U(4)=\frac43\left(\ln4-\frac34\right) =0.848392481493\ldots$$

agrees with the direct convolution and the refined local solution.

Perturbations inherit the prescription

Write $U=\bar U+\delta U$ and perturb the defining equation. To first order,

$$\bar\Box\,\delta U =-\delta R-(\delta\Box)\bar U,$$ $$\bar\Box\,\delta S =-\delta U-(\delta\Box)\bar S.$$

The perturbations of the auxiliaries are sourced by metric perturbations and use the perturbed retarded kernel. They are not arbitrary extra isocurvature data. Setting freely chosen $\delta U$ and $\delta S$ on the initial surface would again enlarge the model. A full prediction must then evolve the coupled metric, matter, and auxiliary perturbations and check kinetic signs, large-scale growth, and the range in which the effective description applies.

Failure diagnostics

Reject a run if refinement changes the solution branch, the retarded initial data drift, direct convolution disagrees, or the discrete equation residual stagnates. In a self-consistent background solve, also monitor the Friedmann constraint, stiffness near rapid transitions, cancellation between $K_{ab}$ terms, and sensitivity to the initial surface. Starting too late can omit curvature memory; starting in a phase where the assumed $R=0$ relation is violated can inject a boundary artifact.

The benchmark validates the auxiliary implementation and prescription on a controlled history. It does not validate a particular cosmological fit, establish a fundamental origin for $m$, or prove stability of the coupled nonlinear RR solution. In particular, changing the initial surface or retaining a homogeneous auxiliary mode can mimic a change in late-time effective energy density; such a run must be labeled as a different branch definition rather than folded into the numerical uncertainty of this one.

Exercises

Derive the nested retarded integrals for $U$ and $S$ on FLRW. For $a=t^{2/3}$ and zero data at $t=1$, obtain $U_{\mathrm{exact}}$ and explain which two independent comparisons in the retained calculation test the local solver.

Solution

Multiplying the first equation by $a^3$ gives

$$\frac{d}{dt}(a^3\dot U)=a^3R.$$

Two integrations with $U(1)=\dot U(1)=0$ produce the first nested integral. Repeating the operation on $d(a^3\dot S)/dt=a^3U$ gives the second.

For $a^3=t^2$ and $R=4/(3t^2)$,

$$\frac{d}{dt}(t^2\dot U)=\frac43.$$

Thus

$$\dot U=\frac43\left(\frac1t-\frac1{t^2}\right), \qquad U=\frac43\left(\ln t+\frac1t-1\right).$$

The first independent comparison is the endpoint agreement of RK4 with SciPy’s nested quadrature, which evaluates the Green-function convolution rather than the local ODE. The analytic $U(t)$ is a second benchmark. The centered finite-difference residual adds a discretization check whose second-order convergence differs from the fourth-order endpoint convergence. Together they test solution value, equation satisfaction, and the retarded boundary data.

References

• Belgacem, Enis, Yves Dirian, Stefano Foffa, and Michele Maggiore (2018), “Nonlocal Gravity: Conceptual Aspects and Cosmological Predictions,” Journal of Cosmology and Astroparticle Physics 2018(03), 002. DOI; Open PDF.
• Maggiore, Michele, and Michele Mancarella (2014), “Nonlocal Gravity and Dark Energy,” Physical Review D 90, 023005. DOI; Open PDF.