CosmoFit.theory

Build a cosmological model from an action rather than from an already-solved expansion history.

Every model in cosmology.models encodes the result of a derivation somebody did by hand: an E(z) typed in from a paper. cosmology.custom.define_model() lowers the barrier to adding one, but not the work – it still wants E(z).

This subpackage takes the input instead of the output. Give it a gravitational action on an FLRW metric, and it reduces the action to a point-like Lagrangian, varies the lapse to get the Friedmann constraint, solves that constraint for E(z), and hands back an ordinary Cosmology subclass that every dataset, likelihood, sampler and plot in this library already knows how to use.

>>> from CosmoFit import Fitter
>>> from CosmoFit.theory import Action
>>>
>>> model = Action(
...     "T + alpha * (-T)**b",
...     geometry="teleparallel",
...     params={
...         "alpha": {"default": 1.0, "bounds": (0.0, 20.0)},
...         "b": {"default": 0.1, "bounds": (-3.0, 0.99)},
...     },
...     closure="alpha",
... ).build("PowerLawFT")
>>>
>>> fit = Fitter(
...     model=model,
...     datasets=["cc", "desi"],
...     free_params=["H0", "Omega_m", "b"],
... )
>>> fit.best_fit()

An action can also carry dynamical scalar fields, in which case the expansion history is integrated rather than solved pointwise:

>>> quintessence = Action(
...     "R",
...     fields={"phi": "X - V0*exp(-lam*phi)"},
...     params={
...         "V0": {"default": 2.1, "bounds": (0.05, 50.0)},
...         "lam": {"default": 0.5, "bounds": (0.0, 1.7)},
...     },
...     closure="V0",
... ).build("ExponentialQuintessence")

A general f(R) is fourth-order and needs its own reduction, which happens automatically:

>>> starobinsky = Action(
...     "R - 2*Lam + alpha_fr*R**2",
...     params={
...         "Lam": {"default": 2.1, "bounds": (0.0, 6.0)},
...         "alpha_fr": {"default": 1e-3, "bounds": (1e-6, 1.0)},
...     },
... ).build("Starobinsky")

The gravitational sector may couple to a field, which is scalar-tensor gravity rather than a field on top of General Relativity – and changes how structure grows as well as how the universe expands:

>>> scalar_tensor = Action(
...     "(1 + xi*phi**2)*R",
...     fields={"phi": "X - V0"},
...     params={
...         "xi": {"default": 0.02, "bounds": (-0.5, 0.5)},
...         "V0": {"default": 2.1, "bounds": (0.05, 20.0)},
...     },
...     closure="V0",
...     growth="quasi_static",     # G_eff/G_N is not 1 here
... ).build("ScalarTensor")

See minisuperspace for the reduction itself, curvature for the Lagrange-multiplier route a general f(R) takes and why it integrates backwards, solve for how a transcendental constraint is solved on the physical branch, and fields for why a field’s initial conditions are set early rather than today.

The action

User-facing action specification, and the compiler that turns one into a fittable Cosmology.

minisuperspace does the physics – it reduces an action on FLRW to a Friedmann constraint. This module wraps that in something a user can actually write down, and closes the loop by solving the constraint for E(z) and handing back a model class that every existing dataset, likelihood, fitter and plot in CosmoFit already knows how to use.

The distinction from cosmology.custom.define_model() is where the derivation happens. define_model takes E(z), which means the user has already done the variational calculus by hand. Action takes the action, and does that calculus.

Example

Lambda-CDM, written as an action rather than as an answer:

>>> from CosmoFit.theory import Action
>>> LCDM = Action("R - 2*Lam", closure="Lam").build("LCDM_from_action")
>>> c = LCDM(H0=67.4, Omega_m=0.315)
>>> c.E(0.5)
1.2795...

The exponential f(Q) model of Anagnostopoulos, Basilakos & Saridakis, whose Friedmann equation this rederives:

>>> FQ = Action(
...     "Q * exp(lam * Q0 / Q)",
...     geometry="symmetric",
...     params={"lam": {"default": 0.3}},
...     closure="lam",
...     growth="quasi_static",
... ).build("FQ_from_action")
class CosmoFit.theory.action.Fluid(w, parameter)[source]

Bases: object

A perfect fluid with constant equation of state.

Parameters:

  • w (float) – Equation of state p / rho. 0 for pressureless matter, 1/3 for radiation.

  • parameter (str) – Name of the cosmological parameter carrying this fluid’s density parameter today. Its value enters the action as rho_0 = 3 Omega (see the units note in minisuperspace).

class CosmoFit.theory.action.Action(gravity, *, geometry=None, fluids=('matter',), params=None, closure=None, growth='gr', fields=None, z_init=3000.0)[source]

Bases: object

A gravitational action on an FLRW background, and the machinery to turn it into a fittable model.

Parameters:

  • gravity (str) –

    The gravitational Lagrangian f, as an expression in the geometry scalar (R, T or Q), the cosmological parameters, and any parameters declared in params. R0/T0/Q0 are available as the scalar’s value today for a matter-free de Sitter normalization (T0 = -6, Q0 = 6 in H0 = 1 units), which is how the f(Q) literature writes its models.

    Written in the sign convention of minisuperspace: an undeformed f (just "R", "T" or "Q") is exactly General Relativity.

  • geometry ({"metric", "teleparallel", "symmetric"}, optional) – Which geometric formulation the action is built on – curvature, torsion, or non-metricity. Inferred from which scalar symbol appears in gravity when only one does.

  • fluids (sequence, optional) – Matter content, as names from STANDARD_FLUIDS or Fluid instances. Defaults to pressureless matter alone.

  • params (dict, optional) – Model parameters beyond the standard set, in the same {name: {"default": ..., "bounds": ..., "label": ...}} form cosmology.custom.define_model() takes.

  • closure (str, optional) – Name of the one parameter fixed by requiring E(0) = 1 rather than fit. Every consistent model needs this condition satisfied somehow: in Lambda-CDM it is what makes Omega_de0 = 1 - Omega_m - Omega_k, and in the exponential f(Q) model it is what makes lam a function of Omega_m rather than a free parameter. If the action satisfies it identically, leave this unset; build() checks and complains either way.

  • growth ({"gr", "quasi_static"}, optional) –

    How the linear growth of structure responds to the modification, i.e. what mu = G_eff/G_N the generated model reports.

    "gr" (default) leaves mu = 1: correct whenever the action modifies the background only in the sense that the extra gravitational degrees of freedom do not propagate into sub-horizon clustering – and the safe default, because it is what every dark-energy model assumes.

    "quasi_static" asks for the sub-horizon quasi-static result, which depends on what the action modified:

    • a deformed geometry scalar (f(T), f(Q)) gives mu = 1/f', what this library’s hand-written FQExponential and FRTLinear already use;

    • a field coupled to the curvature gives the scalar-tensor result of Boisseau, Esposito-Farese, Polarski & Starobinsky (2000), mu = (2F + 4 F_phi^2) / (F (2F + 3 F_phi^2)) with F = df/dR, evaluated on the field’s own solution so that it moves as the field rolls.

    Either way it is an additional physical assumption on top of the action – a statement about perturbations, which a background action does not by itself determine – so it must be asked for explicitly. Asking for it where nothing is modified is an error rather than a no-op.

  • fields (dict, optional) –

    Scalar fields, mapping name -> Lagrangian density written in terms of the field and its kinetic scalar X ("X - V0*exp(-lam*phi)" for exponential quintessence, or any other L(X, phi) for k-essence). More than one is fine; each gets its own equation of motion.

    A field’s name is also in scope in gravity, which is how scalar-tensor gravity is written – "(1 + xi*phi**2)*R" couples the field to curvature rather than adding it on top of General Relativity, and brings the 3 H dF/dt term into the Friedmann equation.

    A field’s expansion history is integrated rather than solved pointwise – see fields – and it adds two parameters, <name>_i and d<name>_i, its value and dphi/dN at z_init. Because those are given early rather than today, E(0) = 1 becomes a shooting condition, so closure is required.

  • z_init (float, optional) – Redshift at which a field’s initial conditions are set, and the earliest redshift the resulting model can be evaluated at. Default 3000, which covers recombination. Ignored by an action with no fields.

property field_lagrangians: dict

Each field’s Lagrangian density, as {name: (expression, X_symbol)}.

lagrangian()[source]

The reduced point-like Lagrangian of this action, together with the minisuperspace it lives on.

Return type:

tuple[Minisuperspace, Expr]

constraint()[source]

The Friedmann constraint of this action, as an expression in E2 (the squared dimensionless Hubble rate) and z, vanishing on-shell.

Returns:

(expr, E2, z) ((sympy expression, Symbol, Symbol))

field_equations()[source]

Equation of motion for each scalar field in this action (each vanishing on-shell), in the gauge N = 1.

Provided for inspection – build() integrates this system rather than returning it, so this is the way to see what was actually derived.

Return type:

dict

closure_equation()[source]

The condition E(0) = 1, as an expression in the model parameters that must vanish.

Every model has to satisfy this – H(z=0) is H0 by definition of H0. What differs is whether the action satisfies it identically (nothing to do), or fixes one parameter in terms of the others.

build(name, label=None)[source]

Compile this action into a Cosmology subclass.

The returned class is an ordinary CosmoFit model: pass it to Fitter and every dataset, likelihood, sampler and plot works on it unchanged.

An action carrying dynamical scalar fields is integrated rather than solved pointwise – see fields – and gains two parameters per field, <name>0 and d<name>0, its value and dphi/dN at a = 1.

Parameters:

  • name (str) – Model name, used for the class and for MODEL_NAME.

  • label (str, optional) – LaTeX name for figures, as in cosmology.custom.define_model().

Return type:

type

field_system()[source]

The compiled equations of motion of this action, for the case where it carries dynamical fields.

Returns (system, args) – see theory.fields.build_system(). Cached: assembling it means solving a linear system symbolically.

curvature_system()[source]

The compiled equations of motion of a general f(R) action. Cached: assembling it means a symbolic solve.

property coupling: Expr

F, the coefficient of R/2 in this action – which is d f / d R. Constant for a minimally coupled field, and a function of the field for scalar-tensor gravity.

property is_non_minimal: bool

Whether the gravitational sector couples to a field – scalar-tensor gravity, where the field sets the strength of gravity rather than sitting on top of it.

property is_fourth_order: bool

Whether this action is a general f(R) – nonlinear in the Ricci scalar, and so fourth-order in the metric.

Those need the Lagrange-multiplier reduction of curvature; everything else goes through the ordinary one.

property namespace: dict

Every symbol a user expression of this action may mention.

solver_arguments(symbols)[source]

The parameter symbols a compiled expression needs, in a fixed order – and a check that every one of them is actually a parameter the generated model will carry.

Return type:

tuple

The reduction

Minisuperspace reduction: from an action to the Friedmann equation.

Every model in cosmology.models was written the other way round – somebody derived E(z) by hand from a Lagrangian and typed the result in. This module does the derivation itself, so a model can be specified the way it is specified in a paper (an action and a metric) rather than the way it is specified in code (an already-solved expansion history).

The method is the standard minisuperspace one. Write FLRW with an explicit lapse N(t),

ds^2 = -N(t)^2 dt^2 + a(t)^2 dSigma_k^2,

substitute it into the action to get a point-like Lagrangian L(N, Ndot, a, adot, addot, fields, fielddots), and vary. The lapse is what makes this work: it is a non-dynamical gauge degree of freedom, so varying it produces a constraint rather than an evolution equation – and that constraint, evaluated at the gauge choice N = 1, is exactly the Friedmann equation. Dropping the lapse from the metric (writing N = 1 from the start) loses it, which is why it has to be carried through the whole reduction and only set to 1 at the very end.

Because L contains addot for curvature-based gravity, the reduction integrates by parts first, discarding the resulting total time derivative – legitimate exactly when L is linear in addot, which holds for General Relativity, a non-minimally coupled F(phi) R, and every f(T) / f(Q) (whose scalars contain no addot at all). It fails for a general f(R), which is genuinely fourth-order; reduce_order() raises rather than silently discarding a term that is not a total derivative, and curvature handles that case by making R an independent variable first – after which the reduction here applies again unchanged.

Units

Everything here is dimensionless: kappa = 8 pi G = 1, H0 = 1, time in units of 1/H0. In these units the GR Friedmann equation is 3 H^2 = rho, so a fluid with today’s density parameter Omega_i has rho_i0 = 3 Omega_i, and spatial curvature enters as k = -Omega_k.

Sign conventions

The three geometry scalars carry different signs in the literature, and picking one at random silently flips the sense of every modification. They are fixed here by one requirement: the undeformed action f = R, f = T or f = Q must reproduce General Relativity exactly. That gives

R = 6 [ addot/(a N^2) + adot^2/(a^2 N^2)

  • adot Ndot/(a N^3) + k/a^2 ], L_g = +N a^3 f(R)/2

T = -6 adot^2/(a^2 N^2) (= -6 H^2), L_g = +N a^3 f(T)/2

Q = +6 adot^2/(a^2 N^2) (= +6 H^2), L_g = -N a^3 f(Q)/2

and test_theory.py asserts the GR limit for all three rather than trusting the table.

CosmoFit.theory.minisuperspace.GEOMETRIES = ('metric', 'teleparallel', 'symmetric')

The three geometry sectors this module can reduce.

CosmoFit.theory.minisuperspace.GEOMETRY_SCALAR = {'metric': 'R', 'symmetric': 'Q', 'teleparallel': 'T'}

Symbol each sector’s gravitational Lagrangian is a function of.

class CosmoFit.theory.minisuperspace.Minisuperspace(fields=(), curvature=False)[source]

Bases: object

The symbols and functions of one FLRW minisuperspace.

Holds the time coordinate, the lapse N(t), the scale factor a(t), the curvature constant k and any scalar fields, so the pieces below all talk about the same objects.

Parameters:

  • fields (sequence of str, optional) – Names of scalar fields living on this background. Each becomes a sympy.Function of t.

  • curvature (bool, optional) – Carry the Ricci scalar as an independent dynamical variable R(t) rather than as shorthand for a combination of a and its derivatives. That is what makes a general f(R) reducible – see curvature.

property coordinates

a, every field, and the curvature variable where there is one.

Type:

The dynamical coordinates

d(expr, n=1)[source]

n-th time derivative of expr.

Parameters:

n (int)

geometry_scalar(geometry)[source]

The curvature / torsion / non-metricity scalar of this minisuperspace, in the sign convention fixed in the module docstring.

Returns:

(expr, sign) ((sympy expression, int)) – expr is the scalar; sign multiplies the gravitational Lagrangian N a^3 f(scalar) / 2 so that an undeformed f gives General Relativity.

Parameters:

geometry (str)

CosmoFit.theory.minisuperspace.fluid_lagrangian(ms, densities)[source]

Minisuperspace Lagrangian of a set of perfect fluids.

A fluid with constant equation of state w has rho(a) = rho_0 a^{-3(1+w)}, and contributes -sqrt(-g) rho = -N a^3 rho(a) = -N rho_0 a^{-3w}. Pressureless matter (w = 0) therefore contributes a term with no scale-factor dependence at all, and radiation (w = 1/3) one going as 1/a – the familiar forms.

Parameters:

  • densities (dict) – Maps w (as a sympy-compatible number) to today’s density rho_0 (in units where 3 H0^2 = 1, i.e. rho_0 = 3 Omega).

  • ms (Minisuperspace)

Return type:

Expr

CosmoFit.theory.minisuperspace.gravity_lagrangian(ms, geometry, f, scalar_symbol)[source]

Minisuperspace Lagrangian of the gravitational sector, sign * N a^3 f(scalar) / 2 with the scalar substituted for its FLRW expression (kappa = 1).

f is given as an expression in scalar_symbol (plus any model parameters and fields); the substitution happens here so that everything downstream differentiates through ordinary sympy chain rules, with no abstract Function left to trip over.

Parameters:
Return type:

Expr

CosmoFit.theory.minisuperspace.field_lagrangian(ms, lagrangians)[source]

Minisuperspace Lagrangian of the scalar-field sector.

Each field’s Lagrangian density is given as an expression in the field itself and its kinetic scalar X; on FLRW,

X = -g^{mu nu} d_mu phi d_nu phi / 2 = phidot^2 / (2 N^2),

so a canonical field is X - V(phi) and a k-essence one is any other function of (X, phi). The contribution to the point-like Lagrangian is sqrt(-g) L = N a^3 L(X, phi).

Parameters:

  • lagrangians (dict) – Maps field name -> (lagrangian_expr, X_symbol).

  • ms (Minisuperspace)

Return type:

Expr

CosmoFit.theory.minisuperspace.reduce_order(L, ms)[source]

Remove second time derivatives from L by integrating by parts and discarding the total derivative.

For a Lagrangian linear in qddot, writing L = A qddot + B and dropping d(A qdot)/dt leaves L - A qddot - Adot qdot, which has the same equations of motion. This is what turns the Einstein-Hilbert term (which contains addot) into the familiar -3 a adot^2 / N.

Raises:

ValueError – If L is nonlinear in some qddot – the case of a general f(R), which is genuinely fourth-order. Then the discarded piece is not a total derivative and dropping it would quietly change the theory, so this refuses instead.

Parameters:
Return type:

Expr

CosmoFit.theory.minisuperspace.euler_lagrange(L, q, t)[source]

Euler-Lagrange expression dL/dq - d/dt (dL/dqdot).

Parameters:

L (Expr)

Return type:

Expr

CosmoFit.theory.minisuperspace.friedmann_constraint(L, ms)[source]

The Friedmann equation: vary L with respect to the lapse, then fix the gauge N = 1.

Returned as an expression that vanishes on-shell (the equation is constraint == 0). The overall normalization is whatever the variation produces – only the zero set matters – so a comparison against a textbook form should allow an arbitrary non-zero factor.

Parameters:
Return type:

Expr

CosmoFit.theory.minisuperspace.field_equations(L, ms)[source]

Equation of motion for each scalar field, in the gauge N = 1. Each vanishes on-shell.

Parameters:
Return type:

dict

Solving the constraint

Solving a Friedmann constraint for E(z).

action produces an expression C(E2, z; params) that vanishes on-shell. Turning that into the E(z) every distance, likelihood and plot in CosmoFit calls means solving it, and the constraint of an interesting model is rarely a polynomial – the exponential f(Q) model’s is transcendental, and its hand-written counterpart in this library needs a Lambert W to invert. So there are two paths here:

  • a closed form, when sympy can solve the constraint for E2 and returns exactly one branch. Vectorized and exact.

  • a continuation solve otherwise. Newton’s method on the whole redshift array at once, walked out from z = 0 – where the closure condition guarantees E2 = 1 – along a ladder of intermediate redshifts s*z, each step seeded by the last.

The continuation is what makes branch selection well-posed. A transcendental constraint generally has several roots, and a root-finder handed the whole redshift range at once can land on a different branch at different z and return a discontinuous, physically meaningless E(z) without failing. Starting from the known root at z = 0 and never letting the solution jump follows one branch by construction – the same reason FQExponential has to pick a specific Lambert W branch rather than any root of its equation.

The derivative dE/dz is never finite-differenced: implicit differentiation of the constraint gives

dE2/dz = -(dC/dz) / (dC/dE2)

exactly, from expressions sympy already has.

CosmoFit.theory.solve.compile_constraint(C, E2, z, args)[source]

Lambdify a constraint and its two partial derivatives.

Only the numerator of the constraint is used. A constraint assembled from a Lagrangian routinely carries a denominator in a (equivalently 1 + z) which never vanishes on the physical domain, and clearing it keeps Newton’s method away from spurious poles.

Parameters:

  • args (sequence of Symbol) – Parameter symbols, in the order the returned callables expect them after (E2, z).

  • C (Expr)

  • E2 (Symbol)

  • z (Symbol)

Returns:

(value_and_slope, df_dz) (tuple of callables) – value_and_slope(E2, z, *params) returns the pair [C, dC/dE2]; df_dz(E2, z, *params) returns dC/dz.

CosmoFit.theory.solve.closed_form(C, E2, z, args)[source]

Try to solve the constraint for E2 symbolically.

Returns a lambdified E2(z, *params) if sympy finds exactly one solution, else None – more than one branch is treated as “no closed form” rather than resolved by guesswork, and left to the continuation solve, which picks the branch physically.

Parameters:

  • C (Expr)

  • E2 (Symbol)

  • z (Symbol)

CosmoFit.theory.solve.solve_E2(z, functions, values, *, seed=None)[source]

Solve C(E2, z) = 0 for E2 on a whole redshift array.

Walks out from z = 0, where E2 = 1 by the closure condition, along _CONTINUATION_STEPS intermediate redshifts, running Newton’s method at each. Every point of the array is advanced together, so this costs a fixed number of vectorized passes rather than a root-find per redshift.

Parameters:

  • z (ndarray) – Redshifts. May be in any order, and may include negative values (the future), which the continuation handles the same way.

  • functions (tuple) – (value_and_slope, df_dz) from compile_constraint().

  • values (tuple) – Parameter values, matching the order given to compile_constraint().

  • seed (ndarray, optional) – Starting guess. When given, the continuation ladder is skipped and Newton runs directly, on the understanding that the seed is already near the right branch – the residual is checked either way, and a failed direct solve falls back to the full ladder. Defaults to None, i.e. walking out from E2 = 1 at z = 0.

Returns:

ndarray – E2 at each redshift.

Raises:

RuntimeError – If Newton fails to converge, rather than returning a number that merely looks like a solution.

CosmoFit.theory.solve.dE2_dz(z, E2, functions, values)[source]

dE2/dz by implicit differentiation of the constraint – exact, given E2 already solved, and free of the truncation error a finite difference would carry into background.q().

Dynamical fields

Integrating an action that carries dynamical scalar fields.

solve handles actions whose Friedmann constraint determines E(z) on its own – f(T), f(Q), General Relativity with a cosmological constant. Add a scalar field and that stops being true: the constraint now involves the field and its velocity, which have their own equation of motion, and the expansion history has to be integrated rather than solved pointwise.

The system

In e-folds N = ln a, with u = dphi/dN, the state is (H, phi, u) and the equations come out of what minisuperspace already derived, with no further physics typed in:

  • the Friedmann constraint C(a, H, phi, u) = 0

  • one Euler-Lagrange equation per field

The trick is that the constraint is a first integral of the equations of motion – so instead of separately deriving the acceleration equation (varying a), this differentiates the constraint along the solution,

dC/dN = a dC/da + (dC/dH) H’ + sum_i [ (dC/dphi_i) u_i

  • (dC/du_i) u_i’ ] = 0,

and solves that together with the field equations for the unknowns H' and u_i'. The two routes are equivalent by the Bianchi identity, and this one reuses expressions already in hand. It also leaves the constraint as an independent check on the integration: it is imposed only in the initial conditions, so how far it drifts along the solution measures the error.

Initial conditions, and why they are set early

Each field contributes two parameters – phi_i and dphi_i, its value and dphi/dN – imposed at an early time z = z_init, and the system is integrated forwards from there.

Setting them at a = 1 instead would be far more convenient: H = 1 holds there by the definition of H0, so the closure condition E(0) = 1 would be algebraic and no shooting would be needed at all. It is also wrong. Integrating backwards from a field at rest today, the Hubble friction -3u that damps the field forwards in time becomes anti-friction, the generic past solution is kinetic-dominated, and rho_phi grows as a^-6: for exponential quintessence with the potential normalized to today’s dark-energy density, that gives E(2.5) = 9.3 where Lambda-CDM gives 3.7. Nothing about the calculation fails – it is the correct past of those initial conditions, and those initial conditions are not a universe.

So the field’s state is specified where a quintessence model actually specifies it, early and typically frozen (dphi_i = 0), and E(0) = 1 becomes a shooting condition on whichever parameter Action(closure=...) names. Forwards is also the numerically stable direction, since the friction term that amplified error backwards now damps it.

The cost is that closure needs a handful of integrations rather than one evaluation, and that E(z) is only defined out to z_init – beyond which there is no history, and asking for one is an error rather than an extrapolation.

class CosmoFit.theory.fields.FieldSystem(constraint, equations, state, derivatives, args, densities=None)[source]

Bases: object

The equations of motion of an action with dynamical fields, compiled to a right-hand side solve_ivp can integrate.

Parameters:

  • constraint (sympy expression) – The Friedmann constraint, already written in terms of the state symbols.

  • equations (list of sympy expressions) – One Euler-Lagrange equation per field, in the same terms.

  • state (tuple) – (a, H, phis, us) – the symbols the above are written in.

  • derivatives (tuple) – (dH, dus) – the symbols standing for dH/dN and each du_i/dN.

  • args (tuple of Symbol) – Parameter symbols, in the order the compiled callables take them.

symbolic_constraint

the closure condition E(0) = 1 is this same constraint at a = 1, H = 1, which is how a field action reuses Action(closure=…) unchanged.

Type:

Kept symbolically as well

H_at(a, fields, velocities, values)[source]

H at scale factor a for a given field state, from the Friedmann constraint.

Return type:

float

CosmoFit.theory.fields.build_system(action)[source]

Assemble the FieldSystem of an Action.

Everything here is substitution: the reduction and the variation already happened in minisuperspace. What changes is the independent variable, from coordinate time to e-folds, via

adot = a H, addot = a (H^2 + H H’), phidot = H u, phiddot = H u H’ + H^2 u’

with ' meaning d/dN.

Return type:

tuple[FieldSystem, list]

class CosmoFit.theory.fields.History(N_lo, N_hi, splines, drift)[source]

Bases: object

A solved expansion history: H and dH/dN as functions of N = ln a, over the range that has been integrated so far.

Held as a cubic Hermite spline with derivatives taken from the equations of motion rather than estimated from the samples – the same construction the rest of the library uses for its distance integrals.

splines

H first, then each field, then each field velocity. The fields are kept because the dark-energy density and pressure are read off them rather than off E(z).

Type:

One spline per state variable

drift

Largest relative drift of the Friedmann constraint over the solution. The constraint is imposed only in the initial conditions, so this is an independent measure of the integration error.

state(N)[source]

Every state variable at N, in the solver’s order.

CosmoFit.theory.fields.initial_state(system, values, fields, velocities, a_i)[source]

The state vector [H, *phis, *us] at a = a_i, with H taken from the Friedmann constraint.

CosmoFit.theory.fields.expansion_today(system, values, fields, velocities, a_i)[source]

H at a = 1 from integrating forwards out of a_i.

This is the residual the closure condition drives to 1, so it runs without t_eval – only the endpoint is wanted, and asking for samples along the way would pay for output nobody reads.

Return type:

float

CosmoFit.theory.fields.integrate(system, values, fields, velocities, a_i, N_hi)[source]

Integrate the field system forwards from a_i and build the interpolating history.

Parameters:

  • system (FieldSystem)

  • values (tuple) – Parameter values, in the order system expects.

  • fields (sequence of float) – The field state at a_i.

  • velocities (sequence of float) – The field state at a_i.

  • a_i (float) – Initial scale factor, 1/(1 + z_init).

  • N_hi (float) – Upper end of N = ln a to cover. Positive when the future is wanted (w(z) and the deceleration parameter are plotted there).

Return type:

History

Fourth-order actions

Fourth-order gravity: a general f(R).

minisuperspace reduces an action by integrating the addot in the Einstein-Hilbert term away by parts, which is legitimate only while the Lagrangian is linear in it. A general f(R) is not, and the reduction refuses – for good reason: the term it would discard is not a total derivative, and dropping it returns a different theory with nothing to show for it.

The Lagrange-multiplier route

The standard way round is to stop treating R as shorthand for a combination of a and its derivatives, and make it an independent variable held to that combination by a multiplier:

S = (1/2) integral dt { N a^3 f(R) - lambda (R - R_geom) }.

Varying lambda returns R = R_geom; varying R gives lambda = N a^3 f'(R). Substituting that back,

L = (1/2) N a^3 [ f(R) - f’(R) R + f’(R) R_geom ],

which is linear in addot – because R_geom is, and it now appears only multiplied by f'(R). So the ordinary reduction applies again, at the cost of one extra dynamical variable. That extra variable is the theory’s fourth order, made visible.

The system

What comes out is smaller than “fourth order” suggests. In e-folds, with the state (H, R):

  • varying R gives back R = 6(2H^2 + H dH/dN), i.e. an explicit dH/dN;

  • the Friedmann constraint is linear in dR/dN, so it gives an explicit dR/dN.

Two first-order equations. The price is where dR/dN divides by f''(R): as the theory approaches General Relativity that vanishes, the equation stiffens, and the integration gets slower without ever becoming wrong. That is measured in tests/test_theory_curvature.py rather than asserted.

Direction of integration

Backwards, from today. That is the opposite of what fields does, and the difference is not an inconsistency – it is what the two systems actually do.

A scalar field integrated backwards runs away: the Hubble friction that damps it forwards becomes anti-friction, and the generic past solution is kinetic-dominated. The f(R) scalaron does not. A relative perturbation of 1e-8 in today’s R changes E(z) by less than that out to z = 1100, and loosening the integrator’s tolerance a hundredfold moves it by 1e-8. Both are measured. So there is no shooting here and no early initial condition: H = 1 at a = 1 holds by construction, and the model’s extra freedom is carried by R today.

CosmoFit.theory.curvature.is_higher_order(f, scalar)[source]

Whether f is nonlinear in the geometry scalar, and so needs this module rather than the ordinary reduction.

Parameters:

  • f (Expr)

  • scalar (Symbol)

Return type:

bool

CosmoFit.theory.curvature.multiplier_lagrangian(ms, f, scalar)[source]

The gravitational Lagrangian of f(R) with R promoted to an independent variable, as set out in the module docstring.

Parameters:

  • ms (Minisuperspace) – Must carry a curvature variable (ms.R).

  • f (sympy expression) – f in terms of scalar.

  • scalar (Symbol) – The symbol f is written in.

Return type:

Expr

class CosmoFit.theory.curvature.CurvatureSystem(constraint, curvature_equation, acceleration, state, derivatives, args)[source]

Bases: object

An f(R) background, compiled to a right-hand side in (H, R) over N = ln a.

Built from the two expressions the reduction produces: the equation that comes from varying R (explicit in dH/dN) and the Friedmann constraint (linear in dR/dN).

CosmoFit.theory.curvature.build_system(action)[source]

Assemble the CurvatureSystem of an f(R) Action.

CosmoFit.theory.curvature.integrate(system, values, R_today, N_lo, N_hi)[source]

Integrate the f(R) background outwards from a = 1, where H = 1 by the definition of H0.

Parameters:

  • system (CurvatureSystem)

  • values (tuple) – Parameter values, in the order system expects.

  • R_today (float) – The Ricci scalar today, in units of H0^2. This is the model’s extra initial condition – the freedom a fourth-order theory has and General Relativity does not.

  • N_lo (float) – Range of N = ln a to cover.

  • N_hi (float) – Range of N = ln a to cover.

Return type:

History