# Cosmological Constraints on 4-Dimensional Einstein-Gauss-Bonnet Gravity

C. M. A. Zanoletti  B. R. Hull  C. D. Leonard  and R. B. Mann

###### Abstract

4-Dimensional Einstein-Gauss-Bonnet (4DEGB) gravity has garnered significant
attention in the last few years as a phenomenological competitor to general
relativity. We consider the theoretical and observational implications of this
theory in both the early and late universe, (re-)deriving background and
perturbation equations and constraining its characteristic parameters with
data from cosmological probes. Our investigation surpasses the scope of
previous studies by incorporating non-flat spatial sections. We explore
consequences of 4DEGB on the sound and particle horizons in the very early
universe, and demonstrate that 4DEGB can provide an independent solution to
the horizon problem for some values of its characteristic parameter $\alpha$.
Finally, we constrain an unexplored regime of this theory in the limit of
small coupling $\alpha$ (empirically supported in the post-Big Bang
Nucleosynthesis era by prior constraints). This version of 4DEGB includes a
geometric term that resembles dark radiation at the background level, but
whose influence on the perturbed equations is qualitatively distinct from that
of standard forms of dark radiation. In this limit, only one
beyond-$\Lambda$CDM degree of freedom persists, which we denote as
$\tilde{\alpha}_{C}$. Our analysis yields the estimate
$\tilde{\alpha}_{C}=(-9\pm 6)\times 10^{-6}$ thereby providing a new
constraint of a previously untested sector of 4DEGB.

Keywords. 4-Dimensional Einstein-Gauss-Bonnet gravity, modified gravity,
scalar-tensor theory, parameter estimation, horizon problem.

##  1 Introduction

Cosmology aims to unravel the profound mysteries surrounding dark matter, dark
energy and the early stages of cosmic evolution. Although general relativity
(GR) has for over 100 years enjoyed immense empirical success [1],
cosmological tensions (such as those in the $H_{0}$ and $S_{8}$ parameters
[2]) hint towards the possible need for new physics to explain modern
cosmological data in this era of precision cosmology. The prediction of
singularities is also a troublesome consequence (and some would argue an
impossibility) of GR.

A multitude of proposals have been put forth to address these unresolved
issues in cosmology [3]. One method that has garnered the attention of
theorists, which will be the subject of this paper, is modified gravity.
Stress testing GR by constraining alternative gravity theories with modern
observations is an important step forward in confirming or rejecting its
validity. Modified gravity has the ability to not only affect cosmology but
gravitational phenomena as a whole, giving it an important role in the quest
for a quantum theory of gravity. While this approach can be coupled with
modifications to canonical particle theory, in this study we assume the
Standard Model (SM) of particle physics.

Extensions of general relativity or alternative theories of gravity have been
around since Einstein’s discovery in 1915. This is a vast subject, and as such
the reader is referred to [4] for a concise and thorough review. One important
result for the study of gravitational theory is Lovelock’s theorem [5, 6],
which states that Einstein’s field equations are the unique 4-dimensional
2${}^{\textrm{nd}}$ order equations of motion for gravity that come from an
action principle that solely relies on the metric. To deviate from general
relativity one must therefore make modifications to the Einstein-Hilbert
action (1.1)

$S=\frac{1}{16\pi
G}\int\sqrt{-g}R\differential^{4}x+\int\mathcal{L}_{m}\left(g_{\mu\nu},\psi\right)\differential^{4}x$
(1.1)

in a non-trivial way. A common procedure is allowing for external fields
coupled to the curvature, as well as the addition of non-linear curvature
terms in (1.1). Of interest here is the case of Gauss-Bonnet gravity, where
one allows the unique set of quadratic invariants of the Riemann tensor [7, 5]

$\mathcal{G}=R^{2}-4R_{\nu}^{\mu}R_{\mu}^{\nu}+R_{\rho\lambda}^{\mu\nu}R_{\mu\nu}^{\rho\lambda}$
(1.2)

to be included in the action (1.1). Although the effects of Gauss-Bonnet
gravity have been extensively studied in higher dimensional space-times, this
invariant in 4 dimensions is a total derivative term, and as such does not
contribute to spacetime dynamics. While it was initially reported in [8] that
a pure Einstein-Gauss-Bonnet theory of gravity in 4 dimensions had been
discovered (which we shall refer to as the Glavin-Lin (GL) formulation), it
was established that there is no such theory [9, 10, 11]; see [12] for a
discussion on the formalism and shortcomings of [8].

These developments prompted others to consider how one would go about
constructing a 4-dimensional theory of gravity that contains the Gauss-Bonnet
term whilst retaining second order equations of motion and locality, and at
the same time avoiding so-called ghosts (theoretical instabilities leading to
unphysical consequences). It was shown that this was possible by allowing the
action to contain an additional scalar field, resulting in [13, 14]

$S=\int d^{4}x\sqrt{-g}\frac{1}{16\pi
G}\left[R-2\Lambda+\alpha\left(\phi\mathcal{G}+4G^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi-4(\partial\phi)^{2}\square\phi+2\left((\partial\phi)^{2}\right)^{2}\right)\right]+S_{m}$
(1.3)

which we shall refer to as 4-dimensional Einstein-Gauss-Bonnet (4DEGB)
gravity. This is no longer a metric theory, but instead is within the class of
modified gravity theories referred to as scalar-tensor theories and more
specifically the Horndeski class [15]. A generalization of this theory
containing additional terms dependent on $\phi$ can be obtained from
dimensional reduction of Gauss-Bonnet gravity [16].

Here we investigate some of the implications of 4DEGB gravity for cosmological
observations, contributing to the already extensive efforts [17] at
constraining this theory. Specifically, we consider a largely untested sector
of 4DEGB parameterized by $C$, a free parameter that comes from solutions to
the scalar field equations of motion derived from the action (1.3), in
addition to the coupling constant $\alpha$. Since the GL formulation is
obtained by considering the $D\to 4$ limit of various $D$-dimensional
solutions, it has no counterpart for the parameter $C$.

Several papers [12, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] have investigated
the cosmological, astrophysical and solar system implications of the GL
formulation (and also, to a lesser extent, of 4DEGB theory). Tests from the
Weak-Field limit provide the tightest constraints

$\alpha<10^{10}\;\textrm{m}^{2}\implies\tilde{\alpha}\lesssim 10^{-43}$ (1.4)

from LAGEOS satellite observations [12], where we define the dimensionless
parameter $\tilde{\alpha}$ as $\tilde{\alpha}=\alpha H_{0}^{2}$. Requiring
that atomic nuclei not be shielded by a horizon imposes the strongest bound on
negative $\alpha$ [26], yielding

$-10^{-83}\lesssim\tilde{\alpha}\lesssim 10^{-43}.$ (1.5)

The terms in the field equations coming from 4DEGB gravity will be of the same
order as those from GR for $\tilde{\alpha}\sim\mathcal{O}(1)$ or
$\alpha\sim\mathcal{O}(10^{52}\textrm{ m}^{2})$. This means that in 4DEGB the
$\alpha$ parameter alone provides negligible corrections to Einstein gravity
in all regimes except in the very early universe, extremely close to a black
hole [17], or inside highly compact objects [28].

The equations describing the Weak-Field limit are recovered almost exclusively
by applying the Post-Newtonian expansion, which produces equations with no
counterpart for the parameter $C$. Consequently, the constraints (1.4) on
$\alpha$ from the Weak-Field limit are valid for both the GL formulation and
4DEGB. This means that the parameter $\alpha$ is highly constrained in 4DEGB
[12], while few conclusive constraints on the parameter $C$ exist (see Section
4.1 for a short overview).

The remainder of our paper is structured as follows: in Section 2 we describe
the theoretical background and derive the field and perturbation equations for
4DEGB. In Section 3 we introduce the remaining effects of the parameter
$\alpha$, which will still impact universe dynamics at very early times even
when priors (1.5) are enforced. We show that the sound horizon in 4DEGB
diverges in the limit of early time, requiring the explicit introduction of an
early-time cut-off. We also demonstrate how 4DEGB can provide an alternative
to inflation in solving the horizon problem for some values of the parameter
$\alpha$.

In Section 4.1 we introduce a useful limit that allows us, in some regimes, to
parameterize 4DEGB through a single additional degree of freedom compared to
GR (the parameter $\tilde{\alpha}_{C}\equiv\frac{\alpha C^{2}}{H_{0}^{2}}$).
We show that, in the small-$\alpha$ limit at times substantially past
nucleosynthesis, 4DEGB is described by a parameter that at the background
level introduces a “dark radiation” term and, at the level of the perturbation
equations, modifies structure growth. In Section 4.2 we make use of this
formulation to place novel constraints on this heretofore minimally tested
sector of 4DEGB using cosmic microwave background (CMB) data, finding
$\tilde{\alpha}_{C}\equiv\frac{\alpha C^{2}}{H_{0}^{2}}=(-9\pm 6)\times
10^{-6}$. Section 4.3 provides a qualitative analysis of the perturbative
behaviour of the theory by analysing trends in the 4DEGB linear matter power
spectrum, produced by modifying a GR Boltzmann solver. We conclude in Section
5.

##  2 Theory & Solutions

###  2.1 Theoretical Background

The 4DEGB theory (1.3) is obtained using an approach similar to that of the
$D\to 2$ limit of general relativity [29]. In this method [13, 14] one begins
with the action

$S_{D}^{\mathrm{GB}}=\frac{1}{16\pi G_{D}}\left[\int
d^{D}x\sqrt{-g}(R-2\Lambda)+\frac{\alpha}{D-4}\left(\int
d^{D}x\sqrt{-\tilde{g}}\tilde{\mathcal{G}}-\int
d^{D}x\sqrt{-g}\mathcal{G}\right)\right]+S_{m}$ (2.1)

with no prior assumptions about the curvature of the $D$-dimensional space,
where $\tilde{g}_{\mu\nu}=\exp(\phi)g_{\mu\nu}$ and $\tilde{\mathcal{G}}$ is
the Gauss-Bonnet term (1.2) for the metric $\tilde{g}_{\mu\nu}$; $S_{m}$ is
the matter action. Writing $\tilde{g}_{\mu\nu}$ in terms of ${g}_{\mu\nu}$ and
taking the limit $D\to 4$ yields the action (1.3).

An alternate method for the construction of a 4DEGB theory is the Kaluza-Klein
like approach [16, 30], which begins with the action

$S_{D}=\frac{1}{16\pi G_{D}}\left[\int
d^{D}x\sqrt{-g}(R+\alpha\mathcal{G})\right]+S_{m}$ (2.2)

assuming a spacetime

$ds_{D}^{2}=ds_{4}^{2}+e^{2\phi}d\Sigma_{D-4}^{2}$ (2.3)

where the $(D-4)$-dimensional space with line element $d\Sigma_{D-4}^{2}$ is
maximally symmetric and is parameterized by the curvature $\lambda$, and
$\phi$ is a scalar field. Dimensional reduction to $D=4$ yields

$\displaystyle S=\int d^{4}x\sqrt{-g}$ $\displaystyle\frac{1}{16\pi
G}\left[R-2\Lambda+\alpha\left(\phi\mathcal{G}+4G^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi-4(\partial\phi)^{2}\square\phi+2\left((\partial\phi)^{2}\right)^{2}\right)\right.$
(2.4) $\displaystyle\left.-2\lambda
Re^{-2\phi}-12\lambda(\partial\phi)^{2}e^{-2\phi}-6\lambda^{2}e^{-4\phi}\right]+S_{m}$

with manifest dependence on the curvature parameter $\lambda$.

The variation of (2.4) with respect to the metric yields the gravitational
field equations [13, 16]

$\displaystyle\mathcal{E}_{\mu\nu}\equiv\Lambda
g_{\mu\nu}+G_{\mu\nu}+\alpha\left[\phi
H_{\mu\nu}-2R\left[\left(\nabla_{\mu}\phi\right)\left(\nabla_{\nu}\phi\right)+\nabla_{\nu}\nabla_{\mu}\phi\right]+8R_{(\mu}^{\sigma}\nabla_{\nu)}\nabla_{\sigma}\phi+8R_{(\mu}^{\sigma}\left(\nabla_{\nu)}\phi\right)\left(\nabla_{\sigma}\phi\right)\right.$
(2.5)
$\displaystyle\quad-2G_{\mu\nu}\left[(\nabla\phi)^{2}+2\square\phi\right]-4\left[\left(\nabla_{\mu}\phi\right)\left(\nabla_{\nu}\phi\right)+\nabla_{\nu}\nabla_{\mu}\phi\right]\square\phi-\left[g_{\mu\nu}(\nabla\phi)^{2}-4\left(\nabla_{\mu}\phi\right)\left(\nabla_{\nu}\phi\right)\right]\left(\nabla^{2}\right)^{2}$
$\displaystyle+8\left(\nabla_{(\mu}\phi\right)\left(\nabla_{\nu)}\nabla_{\sigma}\phi\right)\nabla^{\sigma}\phi-4g_{\mu\nu}R^{\sigma\rho}\left[\nabla_{\sigma}\nabla_{\rho}\phi+\left(\nabla_{\sigma}\phi\right)\left(\nabla_{\rho}\phi\right)\right]+2g_{\mu\nu}(\square\phi)^{2}-2g_{\mu\nu}\left(\nabla_{\sigma}\nabla_{\rho}\phi\right)\left(\nabla^{\sigma}\nabla^{\rho}\phi\right)$
$\displaystyle\quad-4g_{\mu\nu}\left(\nabla^{\sigma}\phi\right)\left(\nabla^{\rho}\phi\right)\left(\nabla_{\sigma}\nabla_{\rho}\phi\right)+4\left(\nabla_{\sigma}\nabla_{\nu}\phi\right)\left(\nabla^{\sigma}\nabla_{\mu}\phi\right)+4R_{\mu\nu\sigma\rho}\left[\left(\nabla^{\sigma}\phi\right)\left(\nabla^{\rho}\phi\right)+\nabla^{\rho}\nabla^{\sigma}\phi\right]$
$\displaystyle\left.\quad+3\lambda^{2}e^{-4\phi}g_{\mu\nu}-2\lambda
e^{-2\phi}\left(G_{\mu\nu}+2\left(\nabla_{\mu}\phi\right)\left(\nabla_{\nu}\phi\right)+2\nabla_{\nu}\nabla_{\mu}\phi-2g_{\mu\nu}\square\phi+g_{\mu\nu}(\nabla\phi)^{2}\right)\right]=8\pi
GT_{\mu\nu}$
$\begin{array}[]{r}H_{\mu\nu}=2\left[RR_{\mu\nu}-2R_{\mu\alpha\nu\beta}R^{\alpha\beta}+R_{\mu\alpha\beta\sigma}R_{\nu}^{\alpha\beta\sigma}-2R_{\mu\alpha}R_{\nu}^{\alpha}\right.\\\
\left.-\frac{1}{4}g_{\mu\nu}\left(R_{\alpha\beta\rho\sigma}R^{\alpha\beta\rho\sigma}-4R_{\alpha\beta}R^{\alpha\beta}+R^{2}\right)\right]\end{array}$

while the variation with respect to the scalar field yields

$\displaystyle\mathcal{E}_{\phi}\equiv$
$\displaystyle-\mathcal{G}+8G^{\mu\nu}\nabla_{\nu}\nabla_{\mu}\phi+8R^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-8(\square\phi)^{2}+8(\nabla\phi)^{2}\square\phi+16\nabla^{a}\phi\nabla^{\nu}\phi\nabla_{\nu}\nabla_{\mu}\phi$
(2.6)
$\displaystyle+8\nabla_{\nu}\nabla_{\mu}\phi\nabla^{\nu}\nabla^{\mu}\phi-24\lambda^{2}e^{-4\phi}-4\lambda
Re^{-2\phi}+24\lambda e^{-2\phi}\left[(\nabla\phi)^{2}-\square\phi\right]=0.$

We are reserving Greek indices for spacetime indices whereas Latin indices
will be used for spatial ones. If $\lambda=0$ then these field equations
reduce to those obtained from the action (1.3). One elegant property of this
theory is that the trace of (2.5) yields

$T=4\Lambda-R-\frac{\alpha}{2}\mathcal{G}$ (2.7)

upon using (2.6), where $T$ represents the trace of the stress energy tensor.

In this paper we shall endeavor to place cosmological constraints on solutions
to action (1.3), for which $\lambda=0$.

###  2.2 FLRW Solutions

In this section we will discuss the familiar FLRW metric in the context of
4DEGB gravity. We will first discuss the background spacetime and the
associated field equations with a perfect fluid stress energy tensor. We will
then move onto the perturbed spacetime, looking at scalar perturbations only.
We will treat the perturbations of the stress energy tensor to be those of an
anisotropic fluid, again with only scalar perturbations.

The cosmological background and scalar perturbation equations in 4DEGB gravity
have been previously studied [12, 24]; in this section, we re-visit their
construction afresh. We present our derivation in detail, which leads in the
case of scalar perturbations to a more complete set of equations than
previously presented in the literature. Furthermore, we expand upon the
existing literature by explicitly including the possibility of a non-zero
spatial curvature in the FLRW metric, with the resulting equations presented
in Appendix A.

####  2.2.1 Background spacetime

We begin with the familiar FLRW line element in co-moving coordinates with
arbitrary spatial curvature

$\differential s^{2}=-\differential t^{2}+a(t)^{2}\left(\frac{\differential
r^{2}}{1-kr^{2}}+r^{2}\left(\differential\theta^{2}+\sin(\theta)^{2}\differential\phi^{2}\right)\right)$
(2.8)

where we parameterize to $k=+1,-1,0$ respectively corresponding to positive,
negative, or flat spatial curvature, with $a(t)$ the scale factor. Isotropy
and homogeneity assumptions imply $\phi=\bar{\phi}(t)$ from which we obtain

$(a\ddot{\bar{\phi}}+\dot{a}\dot{\bar{\phi}}-\ddot{a})\left(a^{2}\dot{\bar{\phi}}^{2}-2a\dot{a}\dot{\bar{\phi}}+\dot{a}^{2}+k\right)=0$
(2.9)

from (2.6) using (2.8). We have suppressed the time dependence on
$\bar{\phi}(t)$ and $a(t)$ and the dot represents a derivative with respect to
coordinate time $t$, with overbars denoting background quantities. We will
also make use of the standard Hubble parameter

$H\equiv\frac{\dot{a}}{a}.$

The solutions to (2.9) are

$\displaystyle\dot{\bar{\phi}}$ $\displaystyle=H\pm\frac{\sqrt{-k}}{a}$
(2.10a) $\displaystyle\dot{\bar{\phi}}$ $\displaystyle=H+\frac{A}{a}$ (2.10b)

where the former comes from the second term in (2.9); clearly $k\leq 0$ for a
real solution. The latter solution (2.10b) comes from the first term in (2.9),
with $A$ an arbitrary constant. This is the most general solution, and we will
henceforth take it to be the solution for the background scalar field.

For the field equations (2.5) we will use the perfect fluid stress energy with
components

$T^{\mu}_{\phantom{\mu}\nu}=\text{diag}(-\bar{\rho},\bar{p},\bar{p},\bar{p})$

where both $\bar{\rho}$ and $\bar{p}$ are understood to be functions of
coordinate time $t$. Using this, the $tt$ component of (2.5) yields

$\frac{\left(-\dot{\bar{\phi}}^{4}a^{3}+4\dot{\bar{\phi}}^{3}\dot{a}a^{2}-2\dot{\bar{\phi}}^{2}\left(3\dot{a}^{2}+k\right)a+4\dot{a}\dot{\bar{\phi}}\left(\dot{a}^{2}+k\right)\right)\alpha}{a^{3}}+\frac{\dot{a}^{2}+k}{a^{2}}=\frac{8\pi
G\bar{\rho}}{3}$ (2.11)

which through the use of (2.10b) and the Hubble parameter becomes

$\alpha\left(H^{2}-\frac{A^{2}}{a^{2}}\right)\left(H^{2}+\frac{A^{2}+2k}{a^{2}}\right)+H^{2}+\frac{k}{a^{2}}=\frac{8\pi
G\bar{\rho}}{3}.$ (2.12)

Without loss of generality we can set $A^{2}=-k+C$, obtaining

$\left(H^{2}+\frac{k}{a^{2}}\right)^{2}\alpha+H^{2}+\frac{k}{a^{2}}=\frac{8\pi
G\bar{\rho}}{3}+\frac{\alpha C^{2}}{a^{4}}$ (2.13)

where the arbitrary constant $C$ now acts as a geometrical “dark radiation”
parameter [12]. If $C=0$ we recover the generalization of the Friedmann
equations obtained in in [8].

The other generalized Friedmann equation can be obtained from any one of the
diagonal terms of (2.5), yielding

$\frac{\left(4\left(H^{2}a^{4}+ka^{2}\right)\dot{H}+3H^{4}a^{4}+2kH^{2}a^{2}+A^{2}\left(A^{2}+2k\right)\right)\alpha}{a^{2}}+3H^{2}a^{2}+2\dot{H}a^{2}+k=-8\pi
Ga^{2}\bar{p}$ (2.14)

or, setting $A^{2}=-k+C$,

$\frac{\left(H^{2}a^{2}+k\right)\left(3H^{2}a^{2}+4\dot{H}a^{2}-k\right)\alpha}{a^{2}}+3H^{2}a^{2}+2\dot{H}a^{2}+k=-8\pi
Ga^{2}\bar{p}-\frac{\alpha C^{2}}{a^{2}}.$ (2.15)

The conservation of stress energy yields

$\displaystyle\nabla_{\mu}T^{\mu\nu}=0\quad\rightarrow$
$\displaystyle\dot{\bar{\rho}}a+3\dot{a}\bar{\rho}+3\dot{a}\bar{p}=0$ (2.16)
$\displaystyle\dot{\bar{\rho}}=\frac{-3(\bar{\rho}+\bar{p})\dot{a}}{a}$

and the system (2.13), (2.15), and (2.16) comprises the background equations
for 4DEGB. Note that (2.15) is a consequence of (2.13) and (2.16), as in
general relativity.

Before we consider perturbations about the background solution, it will be
useful to consider these equations in conformal time, where the line element
(2.8) is now written as

$\mathrm{d}s^{2}=a(\eta)^{2}\left(-\mathrm{d}\eta^{2}+\frac{\mathrm{d}r^{2}}{1-kr^{2}}+r^{2}\left(\mathrm{~{}d}\theta^{2}+\sin(\theta)^{2}\mathrm{~{}d}\phi^{2}\right)\right)$
(2.17)

with

$\differential\eta=\frac{\mathrm{d}t}{a(t)}\;.$ (2.18)

Letting prime denote a derivative with respect to $\eta$, the Hubble parameter
becomes

$\mathcal{H}(\eta)\equiv\frac{a^{\prime}(\eta)}{a(\eta)}\quad\quad\Rightarrow\mathcal{H}(\eta)=a(t)H(t)$
(2.19)

with $\bar{\rho}$ and $\bar{p}$ now functions of $\eta$, along with the scalar
field $\bar{\phi}$, with solution

$\bar{\phi}^{\prime}=\mathcal{H}+A\Rightarrow\bar{\phi}=\ln(a)+A\eta+B\;.$
(2.20)

In conformal time the generalized Friedmann equations are now written as

$\displaystyle\frac{\left(\mathcal{H}^{2}+k\right)^{2}\alpha}{a^{2}}+\mathcal{H}^{2}+k=\frac{8\pi
Ga^{2}\bar{\rho}}{3}+\frac{\alpha C^{2}}{a^{2}}$ (2.21)
$\displaystyle\frac{\left(\mathcal{H}^{2}+k\right)\left(4\mathcal{H}^{\prime}-\mathcal{H}^{2}-k\right)\alpha}{a^{2}}+\mathcal{H}^{2}+2\mathcal{H}^{\prime}+k=-8\pi
Ga^{2}\bar{p}-\frac{\alpha C^{2}}{a^{2}}$
$\displaystyle\bar{\rho}^{\prime}=-\frac{3a^{\prime}(\bar{p}+\bar{\rho})}{a}.$

The use of a conformal or co-moving metric is a matter of taste. We shall
employ the former (unless otherwise noted) since $\bar{\phi}^{\prime}$ depends
on the scale factor only when in terms of the Hubble parameter, as we see in
(2.20).

####  2.2.2 Perturbed Spacetime

We will work in the context of the Conformal Newtonian gauge, and only
consider scalar perturbations. All leading order equations are linear order in
the perturbative quantities, which depend on all spacetime coordinates. The
perturbed version of the line element (2.17) is

$\mathrm{d}s^{2}=a^{2}\left(-(1+2\Phi)\mathrm{d}\eta^{2}+(1-2\Psi)\left(\frac{\mathrm{d}r^{2}}{1-kr^{2}}+r^{2}\left(\mathrm{~{}d}\theta^{2}+\sin(\theta)^{2}\mathrm{~{}d}\phi^{2}\right)\right)\right)$
(2.22)

where $\Phi$ and $\Psi$ are the metric scalar perturbations. We follow the
convention of [17] in the labelling of $\Phi$ and $\Psi$. Perturbations of the
stress energy [12, 31] are given by $\delta\rho$ and $\delta p$ for the
density and pressure respectively, as well as the velocity perturbation $v$
and the anisotropic fluid perturbation $\pi_{ij}$ through the following [12,
31]

$\delta T_{0}^{0}=-\delta\rho,\quad\delta
T_{i}^{0}=(1+w)\bar{\rho}\partial_{i}v,\quad\delta
T_{j}^{i}=\delta_{j}^{i}\delta p+\pi_{j}^{i}$ (2.23)

where $\pi$ is defined via

$\pi_{ij}=\left(\nabla_{i}\partial_{j}-\frac{1}{3}\gamma_{ij}\nabla^{2}\right)\Pi$
(2.24)

and where $\Pi$ is a scalar function representing the anisotropy, $\nabla_{i}$
represents a covariant derivative***In the absence of spatial curvature one
may replace $\nabla_{i}$ with $\partial_{i}$ in $\pi_{ij}$, consistent with
[12, 31]. , and $\nabla^{2}$ is the Laplacian with respect to the spatial
metric $\gamma_{ij}$.

The relative energy density perturbation and speed of sound are respectively
defined as

$\bm{\delta}=\frac{\delta\rho}{\bar{\rho}}\qquad c_{s}^{2}=\frac{\delta
p}{\delta\rho}$ (2.25)

and

$\phi\rightarrow\bar{\phi}+\delta\phi$ (2.26)

defines the perturbation of the scalar field.

We now consider the perturbed equations of motion, starting with the scalar
field equation (2.6). Its leading order correction is

$\alpha\left(A^{2}+k\right)\left(\nabla^{2}(\delta\phi-\Phi)+3A(\Psi^{\prime}+\Phi^{\prime})-3(\Psi^{\prime\prime}+\delta\phi^{\prime\prime})\right)=0$
(2.27)

upon using (2.20) and (2.22). Note that it has no dependence on the scale
factor $a$, which would not be the case - as stated above - if we were working
in co-moving coordinates. Setting $A^{2}=-k+C$, this becomes

$\alpha
C\left(\nabla^{2}(\delta\phi-\Phi)+3\sqrt{-k+C}\left(\Psi^{\prime}+\Phi^{\prime}\right)-3\left(\Psi^{\prime\prime}+\delta\phi^{\prime\prime}\right)\right)=0.$
(2.28)

Similar to the discussion above, when $C=0$ the equation is trivially
satisfied and no information about the perturbation $\delta\phi$ is given.

Henceforth we focus on the situation in which $C\neq 0$ and $k=0$ (with the
$k\neq 0$ case presented in Appendix A). Turning now to the gravitational
field equations (2.5), the $\mathcal{E}_{\eta\eta}=8\pi GT_{\eta\eta}$
component is

$2\alpha A^{2}\left(3\Phi
A^{2}-6A(\Psi^{\prime}+\delta\phi^{\prime})-2\nabla^{2}(\Psi+\delta\phi)\right)+\mathcal{A}(2\nabla^{2}\Psi-6\mathcal{H}\Psi^{\prime})-6\alpha\Phi\mathcal{H}^{4}=8\pi
Ga^{2}(a^{2}\delta\rho+2a^{2}\Phi\rho)$ (2.29)

at leading order, where $\delta
T_{\eta\eta}=a^{2}\delta\rho+2a^{2}\Phi\bar{\rho}$ and we follow [12] in
defining

$\mathcal{A}\equiv 2\alpha\mathcal{H}^{2}+a^{2}\qquad\mathcal{B}\equiv
2\alpha\mathcal{H}^{2}-a^{2}$ (2.30)

for convenience. From the zeroth order equation we identify

$\bar{\rho}=\frac{-3A^{4}\alpha+3\alpha\mathcal{H}^{4}+3a^{2}\mathcal{H}^{2}}{8a^{4}\pi
G}$

and inserting this into (2.29) we find

$4\alpha
A^{2}\left(3A^{2}\Phi-3A(\Psi^{\prime}+\delta\phi^{\prime})-\nabla^{2}(\Psi+\delta\phi)\right)+2\mathcal{A}\left(\nabla^{2}\Psi-3\mathcal{H}\Psi^{\prime}-3\Phi\mathcal{H}^{2}\right)=8\pi
Ga^{4}\delta\rho.$ (2.31)

For any time-spatial component $\mathcal{E}_{\eta i}=8\pi GT_{\eta i}$ we
obtain

$\left(4\alpha
A^{2}(A(\Phi-\delta\rho)-\Psi^{\prime}-\delta\rho^{\prime})+2\mathcal{A}\Psi^{\prime}+2\mathcal{H}\mathbf{\mathcal{A}}\Phi\right),_{i}=-8\pi
Ga^{4}(1+w)\rho\partial_{i}v$

to leading order, which readily integrates to

$4\alpha
A^{2}(A(\Phi-\delta\rho)-\Psi^{\prime}-\delta\rho^{\prime})+2\mathcal{A}\Psi^{\prime}+2\mathcal{H}\mathbf{\mathcal{A}}\Phi=-8\pi
Ga^{4}(1+w)\rho v.$ (2.32)

For the spatial-spatial components of $\mathcal{E}_{ij}=8\pi GT_{ij}$ where
$i\neq j$ we have

$\left(2\alpha(A^{2}(\Psi+\Phi)+2\mathcal{H}^{\prime}\Psi)-\mathcal{B}\Psi-\mathcal{A}\Phi\right),_{ij}=8\pi
Ga^{4}\Pi,_{ij}$

to leading order, which again easily integrates into

$2\alpha\left(A^{2}(\Psi+\Phi)+2\mathcal{H}^{\prime}\Psi\right)-\mathcal{B}\Psi-\mathcal{A}\Phi=8\pi
Ga^{4}\Pi$ (2.33)

where we have set any integration functions of $\eta$ to zero. This can be
done without any loss of generality as the constraints must hold for all time;
as such we can tune them to any value we wish. We note that for time-spatial
and mixed spatial components there are no zeroth order counterparts.

Lastly we consider the diagonal components, $\mathcal{E}_{ii}=8\pi GT_{ii}$.
Since we lose homogeneity and isotropy at leading order, the equations are no
longer equivalent for $i=x,y,z$, unlike the zeroth order case. Due to this
feature we take the trace, $\mathcal{E}^{i}_{i}=8\pi GT^{i}_{i}$ which to
leading order can be written as

$\delta G_{xx}+\delta G_{yy}+\delta G_{zz}=8\pi G(3a^{2}(\delta
p-2\Psi\bar{p})).$ (2.34)

Now isolating

$\bar{p}=\frac{-A^{4}\alpha+\alpha
H^{4}-4\dot{H}H^{2}\alpha-H^{2}a^{2}-2\dot{H}a^{2}}{8a^{4}\pi G}$ (2.35)

from any zeroth order diagonal equation and inserting into $\delta
G_{xx}+\delta G_{yy}+\delta G_{zz}=8\pi G(3a^{2}(\delta p-2\Psi\bar{p}))$ we
get the final expression

$\displaystyle
2\alpha\left(3A^{2}\Phi-3A(\delta\phi^{\prime}-\Phi^{\prime})-\nabla^{2}(\Phi+\Psi)-3(\delta\phi^{\prime\prime}+\Psi^{\prime\prime})\right)+\mathcal{A}(3\mathcal{H}\Phi^{\prime}+\nabla^{2}\Phi+3\Psi^{\prime\prime})+\mathcal{B}\nabla^{2}\Psi$
(2.36)
$\displaystyle+6\mathcal{H}^{\prime}\mathcal{C}\Phi+6\mathcal{H}\mathcal{D}\Psi^{\prime}-3\mathcal{H}^{2}\mathcal{B}\Phi-4\mathcal{H}^{\prime}\alpha\nabla^{2}\Psi=12\pi
Ga^{4}\delta p$

where $\mathcal{C}=4\alpha\mathcal{H}^{2}+a^{2}$ and
$\mathcal{D}=2\mathcal{H}^{\prime}\alpha+a^{2}$. We note here a subtle point
regarding the sum of the diagonal components, which yields three equations
that are not identical at leading order. However making use of (2.33) renders
these equations equivalent and hence the sum may be taken without any loss of
generality. This subtle feature sometimes gets taken for granted, even when
discussing General Relativity.

At this point we have 7 unknown scalar functions –
$\Phi,\Psi,\delta\phi,v,\delta\rho,\delta p$ and $\Pi$ – and 5 equations:
(2.28) from the scalar field equation, and (2.31),(2.32),(2.33), and (2.36)
from the gravitational field equations (2.5). The remaining two equations come
from the conservation of stress-energy $\nabla_{\mu}T^{\mu\nu}=0$. At zeroth
order we retrieve only one equation, the cosmological fluid equation

$\bar{\rho}^{\prime}=-3\mathcal{H}(\bar{p}+\bar{\rho})$ (2.37)

obtained above. The first order $\nabla_{\mu}T^{\mu\eta}$ equation is

$\bar{\rho}\left[(1+w)\nabla^{2}v-3(\Psi^{\prime}+2\mathcal{H}\Phi)\right]-3\Psi^{\prime}\bar{p}+\delta\rho^{\prime}-2(3\mathcal{H}\bar{p}+\bar{\rho}^{\prime})\Phi+3\mathcal{H}(\delta
p+\delta\rho)=0$ (2.38)

and after inserting the zeroth order solution (2.37) this becomes

$\delta\rho^{\prime}=(1+w)(3\Psi^{\prime}-\nabla^{2}v)\bar{\rho}-3\mathcal{H}(\delta\rho+\delta
p).$ (2.39)

This can be rewritten as

$\bm{\delta}^{\prime}=3\mathcal{H}\bm{\delta}\left(w-c_{s}{}^{2}\right)+(1+w)\left(3\Psi^{\prime}-\nabla^{2}v\right)$
(2.40)

using
$\bm{\delta}\bar{\rho}=\delta\rho\Rightarrow\delta\rho^{\prime}=\bm{\delta}^{\prime}\bar{\rho}+\bm{\delta}\bar{\rho}^{\prime}$.
Lastly, any $\nabla_{\mu}T^{\mu i}=0$ component - they are all equivalent - in
conjunction with (2.37) yields

$\partial_{i}(2\nabla^{2}\Pi+3\delta
p+3(1+w)\bar{\rho}[(1-3w)v+v^{\prime}+\Phi])=0$ (2.41)

whose divergence is

$2\nabla^{4}\Pi+3\nabla^{2}\delta
p+3(1+w)\bar{\rho}[(1-3w)\nabla^{2}v+\nabla^{2}v^{\prime}+\nabla^{2}\Phi]=0,$
(2.42)

a form that will prove useful later on.

Upon setting $A^{2}=C$, the full set of equations for $k=0$ is

$\mathclap{\begin{aligned} &\alpha
C\left(\nabla^{2}(\delta\phi-\Phi)+3\sqrt{C}\left(\Psi^{\prime}+\Phi^{\prime}\right)-3\left(\Psi^{\prime\prime}+\delta\phi^{\prime\prime}\right)\right)=0\\\
&4\alpha
C\left(3C\Phi-3\sqrt{C}\left(\Psi^{\prime}+\delta\phi^{\prime}\right)-\nabla^{2}(\Psi+\delta\phi)\right)+2\mathcal{A}\left(\nabla^{2}\Psi-3\mathcal{H}\Psi^{\prime}-3\Phi\mathcal{H}^{2}\right)=8\pi
Ga^{4}\delta\rho\\\ &4\alpha
C\left(\sqrt{C}(\Phi-\delta\rho)-\Psi^{\prime}-\delta\rho^{\prime}\right)+2\mathcal{A}\Psi^{\prime}+2\mathcal{H}\mathcal{A}\Phi=-8\pi
Ga^{4}(1+w)\bar{\rho}v\\\
&2\alpha\left(C(\Psi+\Phi)+2\mathcal{H}^{\prime}\Psi\right)-\mathcal{B}\Psi-\mathcal{A}\Phi=8\pi
Ga^{4}\Pi\\\ &\begin{aligned}
&2\alpha\left(3C\Phi-3\sqrt{C}\left(\delta\phi^{\prime}-\Phi^{\prime}\right)-\nabla^{2}(\Phi+\Psi)-3\left(\delta\phi^{\prime\prime}+\Psi^{\prime\prime}\right)\right)+\mathcal{A}\left(3\mathcal{H}\Phi^{\prime}+\nabla^{2}\Phi+3\Psi^{\prime\prime}\right)\\\
&+\mathcal{B}\nabla^{2}\Psi+6\mathcal{H}^{\prime}\mathcal{C}\Phi+6\mathcal{H}\mathcal{D}\Psi^{\prime}-3\mathcal{H}^{2}\mathcal{B}\Phi-4\mathcal{H}^{\prime}\alpha\nabla^{2}\Psi=12\pi
Ga^{4}\delta p\end{aligned}\\\
&\bm{\delta}^{\prime}=3\mathcal{H}\bm{\delta}\left(w-\frac{\delta
p}{\delta\rho}\right)+(1+w)\left(3\Psi^{\prime}-\nabla^{2}v\right)\\\
&2\nabla^{4}\Pi+3\nabla^{2}\delta
p+3(1+w)\bar{\rho}\left[(1-3w)\nabla^{2}v+\nabla^{2}v^{\prime}+\nabla^{2}\Phi\right]=0.\end{aligned}}$
(2.43)

The full set of equations with $k\neq 0$ are given in Appendix A.

In what follows, we will consider these equations from two perspectives.
First, we will look at the full solutions at the background level outlined by
equations (2.21), which describe phenomena in the very early Universe. We
shall then investigate an unexplored region of parameter space in which
$\alpha$ is small but $C$ is unconstrained. Mathematically this is most easily
studied by taking $\alpha\rightarrow 0$ but keeping $\alpha C^{2}$ finite.
This will have interesting implications for cosmological observables at CMB
redshifts and later times.

##  3 The very early universe

When $a$ is small or $H$ is large the $\alpha H^{4}$ term in (2.13) will be
non-negligible (even for small $\alpha$). Therefore, in the very early
universe the parameter $\alpha$ will have an effect on background dynamics. We
will primarily explore the effects of $\alpha$ on the sound and particle
horizons in the early universe and on the horizon problem.

###  3.1 Diverging horizons close to the Big Bang

One striking feature of 4DEGB is its effect on the comoving sound horizon
$r_{s}(z)$, where $a=\frac{1}{1+z}$. This is defined as the comoving distance
a sound wave could travel from the beginning of the universe in some time
$t(z)$. The sound horizon at recombination is a fundamental scale shaped by
early universe physics that leaves an imprint on the clustering of matter in
the universe [32]. The equation describing $r_{s}(z)$ is

$r_{s}(z)=\int^{\infty}_{z}\frac{c_{s}(z^{\prime})}{H(z^{\prime})}\text{d}z^{\prime}=c\int^{a}_{0}\frac{1}{H(a^{\prime})a^{\prime
2}\sqrt{3(1+R)}}\text{d}a^{\prime}$ (3.1)

where $c_{s}(z)=\frac{c}{\sqrt{3(1+R)}}$ is the sound speed and
$R=\frac{3\rho_{b}}{4\rho_{r}}=3a/4$ (when assuming collisonless cold dark
matter) where $\rho_{b}$ and $\rho_{r}$ are the energy densities of baryons
and radiation respectively [33].

Let us see what effect $\tilde{\alpha}$ has on this equation. We notice that
(2.13) gives us an expression for $H$:

$H(a)=\frac{H_{0}}{\sqrt{2\tilde{\alpha}}}\sqrt{\frac{2\alpha
k}{a^{2}}-1+\sqrt{1+4\tilde{\alpha}\left(\frac{\Omega_{m0}}{a^{3}}+\frac{\Omega_{r0}}{a^{4}}+\Omega_{\Lambda
0}+\frac{\tilde{\alpha}\tilde{C}^{2}}{a^{4}}\right)}}.$ (3.2)

As a first step, this then allows us to find an analytic solution to (3.1)
close to the Big Bang (in the limit $a\rightarrow 0$). Within GR we have the
standard result

$\begin{array}[]{lr}H\approx H_{0}\sqrt{\Omega_{r0}}a^{-2}\\\ R\rightarrow
0\end{array}\Big{\\}}\implies
r_{s}\approx\frac{c}{H_{0}}\int^{a}_{0}\frac{1}{\sqrt{3\Omega_{r0}}}\text{d}a^{\prime}\approx\frac{1}{H_{0}}\frac{1}{\sqrt{3\Omega_{r0}}}a\rightarrow
0$

whereas for a 4DEGB universe

$\begin{array}[]{lr}H\approx\mathcal{K}H_{0}a^{-1}\\\ R\rightarrow
0\end{array}\Big{\\}}\implies
r_{s}\approx\frac{1}{\sqrt{3}H_{0}\mathcal{K}}\int^{a}_{0}\frac{1}{a^{\prime}}\text{d}a^{\prime}\approx\frac{1}{\sqrt{3}H_{0}\mathcal{K}}\left[\text{ln}a\right]^{a}_{0}$

where
$\mathcal{K}=\sqrt[4]{\frac{\Omega_{r0}+\tilde{\alpha}\tilde{C}^{2}}{\tilde{\alpha}}}$
is a constant. This simplified version of Equation (3.1) diverges close to the
Big Bang for a 4DEGB universe. In general, in a 4DEGB universe $r_{s}$
diverges logarithmically as $a\to 0$ for any $\alpha\neq 0$ if we use Equation
(3.1) to define the sound horizon. This seems to suggest 4DEGB predicts an
infinite comoving sound horizon as long as $\alpha\neq 0$.

This effect is perhaps not as severe as it may seem at first glance. First,
there will in general be a cutoff at small $a$ at the Planck time
$t_{\textrm{Planck}}=5.39\times 10^{-44}\textrm{ s}$, since quantum
gravitational effects will become important. More generally, recall that what
is physically relevant is an effective sound horizon, defined as the distance
traveled by a sound wave that has started propagating after the end of
inflation. Waves that left their sources earlier than the end of inflation
underwent incredibly large Doppler shifts as the universe expanded, so it is
reasonable to consider them to be completely unobservable in practice. If this
is the case, the lower bound of the integral in (3.1) is some scale factor at
the end of inflation rather than at the beginning of the universe. This then
gives us a new equation for our sound horizon:

$r_{s}(a)=\int^{a}_{a_{\text{min}}}\frac{c}{H(a^{\prime})a^{\prime
2}\sqrt{3(1+R)}}\text{d}a^{\prime}$ (3.3)

where $a_{\text{min}}$ is close to the end of inflation,
$a_{\text{min}}\approx 10^{-28}$. We now encounter no divergences when
calculating $r_{s}$. In GR+$\Lambda$CDM, the contribution to the integral from
early times is negligible and the exact value of the lower bound isn’t
important. In 4DEGB, the contributions are large at early times and an exact
value of this lower bound is needed to find the sound horizon.

We can solve this apparent issue by assuming that Big Bang Nucleosynthesis
(BBN) has to occur at the same time (and be described by similar physics) in
4DEGB as in GR because of primordial abundance constraints.

Figure 1: The dimensionless sound horizon as a function of redshift $z$ in GL
4DEGB compared to $\Lambda$CDM, with its residual. We assume fiducial
cosmological parameters from [34] and $a_{\text{min}}\approx 10^{-28}$. We
have set $C=0$ for convenience, but non-zero $C$ will not affect the sound
horizon significantly. We display the sound horizon for values of
$\alpha\lesssim 10^{17}\implies\tilde{\alpha}\lesssim 10^{-36}$ [17].

This requires $\alpha\lesssim 10^{17}\implies\tilde{\alpha}\lesssim 10^{-36}$
[17]. We can then show that for this small value of $\tilde{\alpha}$, if the
end of inflation occurs around $a=a_{\text{min}}\approx 10^{-28}$, the exact
lower bound of (3.3) isn’t important. This small (but non-zero) $\alpha$ now
causes a minimal modification to the sound horizon, which will be virtually
the same in 4DEGB and GR as illustrated in Figure 1. Because the divergence of
the integral in (3.3) is logarithmic, in general the lower cutoff does not
meaningfully change the result. Even taking the cutoff $a_{\textrm{min}}$ to
be at the Planck time (before which any classical theory of gravity is
believed to be invalid), the difference between the $\Lambda$CDM and 4DEGB
solutions would still be negligible for $\tilde{\alpha}\lesssim 10^{-36}$.

In spite of this argument, having to impose a cutoff scale in the equation for
the sound horizon is a theoretically unsatisfying characteristic of this
theory. We have no firm knowledge of the correct $a_{\text{min}}$, and the
fact that we have to pick one is a complication of the theory.

###  3.2 The Horizon Problem

Another interesting feature of 4DEGB in the early Universe is that in the
limit of early time, the first equation in (2.21) reduces to:

$\alpha\mathcal{H}^{4}=\Omega_{r0}H_{0}^{2}+\alpha C^{2}$ (3.4)

and thus it follows that

$a(\eta)=a_{\text{min}}\,\text{exp}\left[H_{0}\sqrt[4]{\frac{\Omega_{r0}+\tilde{\alpha}\tilde{C}^{2}}{\tilde{\alpha}}}(\eta-\eta_{\text{min}})\right]=a_{\text{min}}\,\text{exp}\left[\mathcal{K}H_{0}(\eta-\eta_{\text{min}})\right]$
(3.5)

where $\eta_{\text{min}}$ is the conformal time at which the cutoff scale
factor $a_{\text{min}}$ occurs. In terms of comoving time $\eta$, (3.5) means
the universe isn’t accelerating or decelerating ($\ddot{a}=0$), just expanding
at a constant velocity. This implies that the comoving Hubble radius
$(aH)^{-1}$ is a constant, as shown in Figure 2.

Figure 2: The comoving Hubble radius as a function of scale factor. Its
logarithmic integral is the comoving particle horizon $\eta$. In general
$\mathcal{H}^{-1}=(aH)^{-1}\propto a^{n}$, where for inflationary models
$n=-1$, for radiation domination $n=1$ and in 4DEGB $n=0$. $a_{e}$ is the
transition scale factor between inflation domination and radiation domination.
Here $\lambda$ is an example of a primordial perturbation exiting the horizon
during inflation and re-entering during radiation domination. We can see the
same effect doesn’t occur in 4DEGB. We assume fiducial cosmological parameters
from [34] with $C=0$ and with $\tilde{\alpha}=10^{-112}$ and
$\tilde{\alpha}=10^{-44}$ for the blue and blue-dashed lines respectively. The
cutoff for the integration in (3.6) occurs at scale factors below the lowest
limit of the x-axis for all investigated values of $\alpha$.

Let us now compare the comoving Hubble radius $(aH)^{-1}$ in 4DEGB to that in
$\Lambda$CDM, and see the effect this has on the horizon problem. $\Lambda$CDM
without inflation predicts that patches on the CMB are causally disconnected
for $\theta\gtrsim 1$°, contrary to the observed uniformity of the CMB
temperature map (which is 1 part in $10^{5}$). This is called the horizon
problem. If the universe had an early epoch during which the comoving Hubble
radius decreased or stayed constant, then there would be more conformal time
between the beginning of the universe and the CMB, so the past light cones of
two points on the CMB would have had more time to be in causal contact. In the
standard cosmological model, this is achieved by introducing an inflationary
era at very early times.

Alternatively, as mentioned above, 4DEGB also has this effect at early times
(again, see Figure 2). As we see from Equation (3.5), 4DEGB is characterized
by an early epoch where the Hubble radius $(aH)^{-1}=\dot{a}^{-1}$ is a
constant. We then have a transition from an $\alpha$-dominated to a radiation-
dominated universe. In this scenario, points in the CMB have overlapping past
light cones and therefore originated from a causally connected region of space
(see Figure 3).

Figure 3: The 4DEGB solution compared to the inflationary solution to the
horizon problem, visualized by looking at distances of interest as a function
of conformal time ($\eta$). The plot assumes the exaggerated value of
$\tilde{\alpha}=10^{-7}$ to show the qualitative behaviour of the theory (and
to aid the visual distinction between CMB and reheating space-like surfaces).
We assume fiducial cosmological parameters from [34]. The $\Lambda$CDM+GR
solution is plotted so that it matches 4DEGB at late conformal times. The
reheating space-like surface is defined here as the conformal time at which
the scale factor for the $\Lambda$CDM+GR solution goes to zero.

Another way of thinking about this is by looking at the comoving particle
horizon $\eta(z)$, defined as:

$\eta(z)=\int^{z_{\text{max}}}_{z}\frac{c}{H(z^{\prime})}\text{d}z^{\prime}=\int^{a}_{a_{\text{min}}}\frac{c}{H(a^{\prime})a^{\prime
2}}\text{d}a^{\prime}=\int^{a}_{a_{\text{min}}}\frac{c}{H(a^{\prime})a^{\prime}}\text{dln}a^{\prime}.$
(3.6)

This is the logarithmic integral of the inverse Hubble radius; it diverges
logarithmically as $a_{\text{min}}\rightarrow 0$. A diverging comoving
particle horizon (a diverging integral for (3.6)) means any two particles will
have been in causal contact in some point in the past. However, as established
above in Section 3.1, $a_{\rm min}$ must be finite and therefore Equation 3.6
does not diverge. As a result, 4DEGB won’t necessarily provide a solution to
the horizon problem. Whether it does depends on how large the value of
$\tilde{\alpha}$ is. If we enforce $a_{\text{min}}=a(t_{\text{planck}})$,
4DEGB provides an independent solution to the horizon problem when
$\tilde{\alpha}\gtrsim 5\times 10^{-11}$. Thus, while 4DEGB solves the horizon
problem in principle, in practice this is a much larger $\tilde{\alpha}$ than
current constraints (1.5) allow.

4DEGB doesn’t provide an early period of acceleration and $\Omega_{k}(z)$ is a
constant in the early 4DEGB universe. Therefore, while 4DEGB can in principle
be sufficient to solve the horizon problem, it can never solve the flatness
problem. 4DEGB without inflation is also not sufficient to have primordial
fluctuations exit and then re-enter the horizon as in $\Lambda$CDM+GR with
inflation (again, see Figure 2).

##  4 Constraints from cosmological probes

###  4.1 A new formalism for cosmological parameter constraints

Although the cosmological impact of the parameter $\alpha$ alone is relatively
weak other than in the very early Universe (see Section 3), several studies
have nevertheless placed constraints on $\alpha$ using late-time cosmological
probes [25, 35]. While these do provide valuable independent constraints on
$\alpha$, they are uninformative compared to bounds due to other effects. More
importantly, we have seen in Section 3.1 that to use observational information
that relates to the sound horizon we must impose a strong prior from BBN of
$|\alpha|\lesssim 10^{17}\text{m}^{2}$ [17]. This is the only meaningful way
of using this information without assuming an inflationary model and doing a
more detailed analysis of the sound horizon in 4DEGB. If we take this as a
prior for a cosmological analysis using e.g. CMB and BAO data, the resulting
constraints on $\alpha$ will then be completely prior-dominated.

Given the stringent constraints on $\alpha$ from Equation (1.5), the GL
formulation is very highly constrained. However, the full 4DEGB theory, with
an additional free parameter $C$, has not been confronted with data (for one
exception, see [23], which discusses how one might constrain $C$ from
theoretical considerations and existing observational bounds, but does not
perform a full cosmological parameter inference analysis). Motivated by the
lack of conclusive constraints on the parameter $C$ in the literature, and by
the stringent observational constraints on $\alpha$ as described above, we now
introduce a version of the 4DEGB equations in the limit of small $\alpha$.

We develop the equations of 4DEGB in the limit where terms in $\alpha$,
$\alpha A$, $\alpha A^{2}$ or $\alpha A^{3}$ are negligible compared to terms
in $\alpha A^{4}=\alpha C^{2}$. In this small-$\alpha$ limit and with $k=0$,
Equation (3.2) takes the form:

$H(z)=H_{0}\sqrt{\Omega_{m0}(1+z)^{3}+(\Omega_{r0}+\tilde{\alpha}_{C})(1+z)^{4}+\Omega_{\Lambda
0}}$ (4.1)

where we define a new dimensionless parameter

$\tilde{\alpha}_{C}\equiv\tilde{\alpha}\tilde{C}^{2}$ (4.2)

where $\tilde{C}\equiv CH_{0}^{-2}$, and $\alpha_{C}\equiv\alpha C^{2}$ such
that $\tilde{\alpha}_{C}=\alpha_{C}H_{0}^{-2}$.

In the small-$\alpha$ limit, this will be the only parameter to enter both
background and perturbation equations. We assume $\tilde{\alpha}_{C}$ is
finite and seek to empirically constrain it. At times substantially past
nucleosynthesis, the background is then equivalent to that of a $\Lambda$CDM
universe with a “dark radiation” term that depends on $\tilde{\alpha}_{C}$.

We can now look at the set of equations (2.43) in the same limit. These reduce
to:

$\displaystyle
k^{2}\delta\varphi+3\delta\varphi^{\prime\prime}=3(\Phi^{\prime}+\Psi^{\prime})$
(4.3) $\displaystyle
12\alpha_{C}\Phi-12\alpha_{C}\delta\varphi^{\prime}-2a^{2}(k^{2}\Psi+3\mathcal{H}\Psi^{\prime}+3\mathcal{H}^{2}\Phi)=8\pi
Ga^{4}\delta\rho$ $\displaystyle\Psi^{\prime}+\mathcal{H}\Phi=-4\pi
Ga^{2}(1+\omega)\rho v$ $\displaystyle\Phi=\Psi+8\pi Ga^{2}\Pi$ $\displaystyle
3\Psi^{\prime\prime}+3\mathcal{H}(\Phi^{\prime}+2\Psi^{\prime})+k^{2}(\Psi-\Phi)+(6\mathcal{H}^{\prime}-3\mathcal{H}^{2})\Phi=12\pi
Ga^{2}\delta\rho$
$\displaystyle\bm{\delta}^{\prime}=3\mathcal{H}\bm{\delta}\left(w-\frac{\delta
p}{\delta\rho}\right)+(1+w)\left(3\Psi^{\prime}-\nabla^{2}v\right)$
$\displaystyle 2\nabla^{4}\Pi+3\nabla^{2}\delta
p+3(1+w)\bar{\rho}\left[(1-3w)\nabla^{2}v+\nabla^{2}v^{\prime}+\nabla^{2}\Phi\right]=0$

where $\delta\varphi=\frac{\delta\phi}{\sqrt{C}}$ is the same order of
magnitude as $\Phi$ with respect to the constant $C$. For a derivation of the
equations above, we refer the reader to Appendix B.

This limit of 4DEGB at the background level behaves like a $\Lambda$CDM
universe with an additional contribution to the dark relativistic degrees of
freedom. Interestingly, the parameter $\tilde{\alpha}_{C}$ couples this
background dark radiation term (which has a purely geometric origin) to the
modified equations for structure growth. This offers the prospect of
differentiating 4DEGB from modifications to the standard model of particle
physics that cause dark radiation, since the equations for growth will be
different.

We can also look at equations (146c) and (146d) from [17] in the same limit.
These reduce to:

$\begin{split}&a^{2}\partial^{2}\dot{F_{i}}+16\pi
G(1+\omega)a^{4}\bar{\rho}v_{i}=0\\\ &a^{2}\ddot{\gamma_{ij}}\hskip
6.0pt+2\mathcal{H}a^{2}\dot{\gamma_{ij}}\hskip
6.0pt-a^{2}\partial^{2}\gamma_{ij}=8\pi Ga^{4}\Pi_{ij}.\end{split}$ (4.4)

These tell us that the tensor and vector perturbations reduce to those from
$\Lambda$CDM in the small-$\alpha$ limit, even though the scalar perturbations
are modified. Then any data that depends only on the vector and tensor
perturbations (such as data from gravitational waves) won’t enforce additional
constraints on our parameter $\tilde{\alpha}_{C}$.

In the following subsections we use cosmological probes of background
equations to constrain our parameter $\tilde{\alpha}_{C}$. In Section 4.3 we
also introduce some qualitative effects of the modifications to the scalar
perturbations on the power spectrum.

###  4.2 Background constraints on $\tilde{\alpha}_{C}$ as a non-interacting
dark radiation with a geometrical origin

Tests of dark radiation have been carried out extensively in the literature
(see, e.g., [36, 37, 38, 39, 34]). We can take advantage of existing
constraints and map them onto the question of constraining the
$\tilde{\alpha}_{C}$ parameter of 4DEGB in the small-$\alpha$ limit. We first
place background constraints by using datasets that don’t require the
assumption of a GR+$\Lambda$CDM structure growth history.

We will find it productive here and in the following section to examine how
analyses that consider additional beyond-standard model neutrinos can equally
constrain a cosmological model with standard model neutrinos and 4DEGB
gravity, parameterized by $\tilde{\alpha}_{C}$ in the small-$\alpha$ limit.
Flexibility in the number of such early-Universe relativistic degrees of
freedom is often parameterized by:

$\Omega_{\nu\text{0}}=N_{\text{eff}}\frac{7}{8}\left(\frac{4}{11}\right)^{7/8}\Omega_{\gamma
0}$ (4.5)

where $N_{\text{eff}}$ is the effective number of relativistic neutrino
species (equal to 3.044 for standard model neutrinos only) and $\Omega_{\gamma
0}$ is the dimensionless photon density parameter. Constraints on $N_{\rm
eff}$ from this $\Lambda$CDM+$N_{\text{eff}}$ model exist in the literature
from a number of probes of the cosmological background dynamics (see, e.g.,
[36, 37]), and we will take advantage of this here.

At the level of the background only, we can treat $\tilde{\alpha}_{C}$ as if
it were a contribution from an additional massless particle that neither
interacts nor decays. Doing so, we can write down an equivalence between
$\tilde{\alpha}_{C}$ in 4DEGB and $N_{\rm eff}$ in GR:

$\tilde{\alpha}_{C}=(N_{\text{eff}}-3.044)\frac{7}{8}\left(\frac{4}{11}\right)^{7/8}\Omega_{\gamma
0}=\Delta
N_{\text{eff}}\frac{7}{8}\left(\frac{4}{11}\right)^{7/8}\Omega_{\gamma 0}$
(4.6)

It is important to emphasize that $\tilde{\alpha}_{C}$ has a geometrical
origin — we are not proposing to test a model with a varying number of
neutrinos. We are testing a 4DEGB universe with standard model neutrinos (with
$N_{\text{eff, SM}}=3.044$). This means we will allow, for example, negative
values of $\tilde{\alpha}_{C}$ which correspond to $N_{\text{eff}}<3.044$, but
the number of physical neutrino species $N_{\text{eff, SM}}$ will always be
the standard 3.044. We emphasize also that this equivalence should only be
seen as a valid description of background cosmological observables; those that
are sensitive to the growth of structure require consideration of
modifications to the perturbation equations.

With this in mind, we can initially consider constraints from
${}^{4}\mathrm{He}$ abundances (which give a conservative bound of
$|N_{\text{eff}}|\leq 3.5$ [36]). This maps to an upper bound on
$\tilde{\alpha}_{C}$ of:

$|\tilde{\alpha}_{C}|\lesssim 10^{-5}.$

####  4.2.1 Constraints on $\tilde{\alpha}_{C}$ from the CMB

Next, we will look at data from the Atacama Cosmology Telescope (ACT),
specifically that taken with its second-generation polarization-sensitive
receiver ACTPol [40]. This dataset provides information about the CMB
anisotropies, which in general depend on the perturbative behavior of the
theory as well as the background behaviour. In practice, the effect of small
modifications to the perturbation equations (which will change the Sachs-Wolfe
and Integrated Sachs-Wolfe effects) is apparent only at angular multipoles
$\ell\lesssim 200$ [41]. The ACTPol analysis that we build on here takes a
minimum multipole of 600 in the CMB temperature auto-spectrum, and 350 in the
E-mode polarization auto-spectrum as well as the cross-spectrum between E-mode
polarization and temperature, which means we can treat it effectively as a
background-only probe. We use posterior chains of cosmological parameters
obtained with the actpollite likelihood for a $\Lambda$CDM+$N_{\text{eff}}$
model, which results in $N_{\text{eff}}=2.42\pm 0.43$ (using the maximum
marginalized likelihood at a 68% confidence level) [42].

Figure 4:  ACTPol posterior results for a $\Lambda$CDM universe with varying
$N_{\rm eff}$. We map this analysis onto our scenario of interest, with
$N_{\rm eff}$ playing an equivalent role to $\tilde{\alpha}_{C}$ in 4DEGB for
background-only probes. $\tilde{\alpha}_{C}$ is consistent with zero within
$1.34\sigma$. We are displaying the $1\sigma$, $2\sigma$ and $3\sigma$
confidence regions in the contour plot, and the $1\sigma$ and $2\sigma$
regions in the 1D posterior probability distribution for $\tilde{\alpha}_{C}$.

This translates to a constraint for 4DEGB (with standard neutrinos) of:

$\tilde{\alpha}_{C}=(-3.5\pm 2.4)\times 10^{-6}h^{2}=(-9\pm 6)\times 10^{-6}$
(4.7)

at 68% confidence level, which is consistent with zero at the 1.34$\sigma$
level. $h$ is the dimensionless Hubble constant defined by $H_{0}=100h\textrm{
km}\hskip 1.79997pt\textrm{s}^{-1}\textrm{Mpc}^{-1}$ and we have used the the
maximum posterior value for $h$ from the ACTPol analysis in question in our
calculations.

As we can see in Figure 4, the analysis of [42] prefers a lower
$N_{\text{eff}}$ than the standard model (or in 4DEGB a negative
$\tilde{\alpha}_{C}$). An interesting feature of 4DEGB is that it provides a
modified gravity framework that mimics a continuous variation of
$N_{\text{eff}}$ from its $\Lambda$CDM value of 3.044, including values
$N_{\text{eff}}<3.044$. Few models are well-motivated theoretically for
$N_{\text{eff}}<3.044$ [34], so 4DEGB provides an interesting new precedent,
which should be borne in mind should future constraints also favour a negative
value for $\Delta N_{\rm eff}$. While we cannot draw conclusions from all-sky
CMB Planck data without accounting for a modification to the scalar
perturbation equations, we note here for completeness that this preference for
a lower value of $N_{\text{eff}}$ does not seem to persist in Planck data [34]
for a $\Lambda\text{CDM+}\text{GR}$ model.

###  4.3 Effect of modified perturbations in 4DEGB on the matter power
spectrum

Having considered the background-only effect of $\tilde{\alpha}_{C}$, we now
examine the growth of perturbations in this small $\alpha$ limit of 4DEGB. We
incorporate the 4DEGB background and perturbation equations within a bespoke
cosmological Boltzmann code to study the qualitative behaviour of
small-$\alpha$ 4DEGB in a perturbed universe. We use a modified version of
equations (4.1) and (4.3) (see Appendix C) to recover the matter power
spectrum $P(k)$. The Boltzmann code we use is constructed following the
methodology of [43]. In taking this approach, we prioritize ease of
implementation rather than the computational efficiency which would be gained
by modifying the highly optimised public packages CAMB [44] or CLASS [45]. We
validate our Boltzmann code against these packages in the GR case and find
excellent agreement.

Figure 5: Present-day ($z=0$) linear matter power spectrum for small-$\alpha$
4DEGB. The middle (solid black) curve is the power spectrum for a $\Lambda$CDM
model. The values of $\tilde{\alpha}_{C}$ displayed here are chosen to be of
the same order of magnitude as background constraints of Section 4.2. The
dotted (green) and dash-double dotted (gray) curves show the
$\Lambda$CDM$+N_{\textrm{eff}}$ model with extra neutrino densities equivalent
to the plotted values of $\tilde{\alpha}_{C}$ ($\Delta\Omega_{\nu
0}=\tilde{\alpha}_{C}=8\times 10^{-6}$ and $-8\times 10^{-6}$ respectively),
which correspond to $N_{\textrm{eff}}=5.46$ and $0.62$ respectively. When
varying $N_{\rm eff}$, we fix
$\\{Y_{p},\Omega_{M}^{0},\Omega_{B}^{0},\Omega_{R}^{0},h,A_{s},n_{s}\\}$ to be
the fiducial cosmological parameters from [34]. We approximate all neutrinos
as massless.

The resulting power spectra can be seen in Figure 5. We see that the
difference in $P(k)$ between 4DEGB and $\Lambda$CDM with the equivalent
$N_{\textrm{eff}}$ value is non-negligible for $\tilde{\alpha}_{C}\approx\pm
8\times 10^{-6}$. The difference in the paired spectra is larger at large
scales (low $k$); to understand this, we turn to the second equation in (4.3),
which suggests that the $\tilde{\alpha}_{C}$ term is more dominant when $k$ is
smaller. Beyond this large-scale effect, there is also a scale-dependent
deviation between the two theories as we move to smaller scales (larger $k$).
Combined, these effects suggest that there is potential for the degeneracy
between this sector of 4DEGB and $\Lambda$CDM$+N_{\textrm{eff}}$ to be broken
by observational probes sensitive to changes in the power spectrum (such as
weak gravitational lensing, galaxy clustering, CMB spectra at lower multipoles
which probe perturbation effects, or CMB lensing). Although the fractional
deviations between the two models in Figure 5 are small, upcoming surveys such
as Euclid [46] and the Rubin Observatory Legacy Survey of Space and Time
(LSST, [47]) may be capable of observationally distinguishing these percent-
level effects. Such an analysis is beyond the scope of this paper but is of
great interest for future work.

##  5 Conclusion

4DEGB gravity has broad theoretical interest as a competitor to general
relativity, primarily because it is a 4-dimensional higher-curvature theory of
gravity that enjoys all the basic features of higher-dimensional Lovelock
gravity. It has also provoked substantial observational interest as a theory
with few extra parameters compared to $\Lambda$CDM$+$GR, which has resulted in
the Glavin-Lin version being extensively investigated, and its single
parameter $\alpha$ being well-constrained. Here we have for the first time
studied the implications of scalar-tensor 4DEGB, providing the first
observational constraints on a geometric term parametrized by a constant $C$.

Considering a background spacetime that is homogeneous and isotropic and a
perturbed spacetime that only depends on scalar functions, we find that the
stress-energy perturbations are that of a perfect fluid with anisotropic
stress. This in turn dictates that there are 7 unknown functions correcting
the cosmic background. We have provided the first derivation of the complete
set of equations necessary for computing these 7 unknown functions for both
flat and curved spatial hyper surfaces. The conservation of stress-energy
provides two equations of motion (that do not differ from GR), while the
gravitational and scalar field equations provide the other five. The complete
set forms a nonlinear coupled PDE system that does not have a closed form
solution. These equations become the first order perturbation equations of
general relativity in the appropriate limit.

We have explored some consequences of $\alpha$ in the very early universe,
where the small, allowable values of this parameter can still have important
cosmological effects. We find that the sound horizon in 4DEGB depends on the
behaviour of the universe at early times, meaning that the value of the
parameter $\alpha$ and the exact redshift at which sound propagation begins to
matter both affect the measured solution for the sound horizon. We show that
this issue can be alleviated somewhat, though this requires imposing a minimum
scale factor after which 4DEGB becomes relevant. Furthermore, we have shown
that 4DEGB can resolve the horizon problem: the particle horizon diverges at
early times, meaning that all particles will have been in causal contact at
some very early time. This solution fails for very low values of $\alpha$ if
we enforce a cutoff to this classical theory at the Planck time. This is
because the aforementioned divergence is logarithmic and therefore is very
slow to diverge even for a very small minimum time. This theory also fails to
act as a replacement for an inflation-type model in the early Universe because
it does not resolve the flatness problem.

Unlike the Glavin-Lin formulation, full 4DEGB contains a dark radiation-like
term at the background level, whose associated parameter $C$ has remained
largely unconstrained in the literature before now. Introducing a new
parameter $\tilde{\alpha}_{C}=\frac{\alpha_{C}}{H_{0}^{2}}=\frac{\alpha
C^{2}}{H_{0}^{2}}$ we have shown that, under the assumption that $\alpha$ is
small (justified by prior constraints on this parameter [17]), the equations
describing the theory at times substantially past inflation depend only on
$\tilde{\alpha}_{C}$. The late-time behaviour of the theory can then be
constrained as a combination of background “dark radiation” effects and a
modification to structure growth. We have shown that the 4DEGB parameter is
constrained using CMB data from ACTPol to

$\tilde{\alpha}_{C}=\tilde{\alpha}\tilde{C}^{2}=(-9\pm 6)\times
10^{-6}\implies\alpha_{C}=(-10\pm 7)\times 10^{-58}h^{2}\text{
m}^{-2}\approx(-4\pm 3)\times 10^{-58}\text{ m}^{-2}$ (5.1)

where we have used the conversion factor $H_{0}^{2}=(1.17\times
10^{-52})h^{2}\text{ m}^{-2}$. This implies $\alpha_{C}$ is consistent with
zero within $1.34\sigma$. Interestingly, ACTPol data seems to prefer a lower
value of $N_{\textrm{eff}}$ than that suggested from the standard model (or in
4DEGB a negative $\tilde{\alpha}_{C}$). 4DEGB enjoys the interesting feature
of being a geometrical alternative to continuously varying $N_{\text{eff}}$,
including values of $N_{\text{eff}}<3.044$ that are much less theoretically
motivated in other settings. This may prove interesting in the future, should
further cosmological observables favour such values of $N_{\rm eff}$ in a
$\Lambda$CDM+$N_{\text{eff}}$ model.

While we completed a qualitative analysis of the power spectrum in this paper,
modifying an efficient public Boltzmann solver to compute this quantity would
enable further constraints on $\alpha_{C}$ using observables sensitive to
cosmic structure growth. This would enable us to break the degeneracy between
this branch of 4DEGB and GR+$\Lambda$CDM with a varying $N_{\rm eff}$. Figure
5 shows that the power spectrum (and therefore $\sigma_{8}$ and $S_{8}$) will
be modified in 4DEGB compared to $\Lambda$CDM (or
$\Lambda$CDM+$N_{\textrm{eff}}$). It would therefore be interesting to
investigate in future work what impact this sector of 4DEGB might have on the
$S_{8}$ tension.

Finally, we note that there are other potential avenues of exploration of the
cosmological consequences of 4DEGB in the very early universe. We have not
here discussed the implications of a 4DEGB universe in terms of the behaviour
and evolution of super-horizon perturbations, or the implications regarding
early-universe singularities or lack thereof in these theories. A universe
dominated by 4DEGB at early times might have an effect on the scale dependence
or non-gaussianity of the primordial power spectrum which would be interesting
to investigate. As yet no such calculation has been attempted. The work
presented here on the sound horizon and its divergent behaviour as
$a\rightarrow 0$ suggests that these may be fruitful avenues for further work.

## Acknowledgments

This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada and in part by the Philip Robinson Cosmology PhD
studentship. We thank Ian Moss and Ghazal Geshnizjani for helpful discussions,
and Pedro G. Ferreira for guidance in writing the original GR version of the
Boltzmann code modified for this work. We also thank Chris Harrison for kindly
allowing the use of his research computing server and Erminia Calabrese for
providing helpful assistance with accessing ACTPol data products. The Python
libraries SciPy [48], NumPy [49], PyCCL†††https://github.com/LSSTDESC/CCL [50]
and corner.py [51] were significant in enabling this work to be done.

Author contributions: CMAZ performed the analysis of Sections 3 and 4
(including writing the corresponding code) and wrote a significant portion of
the text of this paper. BRH conducted the theoretical calculations of Section
2 and wrote a significant portion of the text of this paper. CDL provided
technical and conceptual guidance and expertise most pertinent to Sections 3
and 4, checked some calculations, wrote the original GR version of the
Boltzmann code used in Section 4, wrote some text for this paper, and revised
the paper text. RBM provided technical and conceptual guidance and expertise
most pertinent to Section 2, checked some calculations, wrote some text for
this paper, and revised the paper text.

For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## References

  * Will [2014] Clifford M. Will.  The confrontation between general relativity and experiment.  _Living Reviews in Relativity_ , 17(1), jun 2014.  doi: 10.12942/lrr-2014-4.  URL https://doi.org/10.12942%2Flrr-2014-4.
  * et al. [2021] Eleonora Di Valentino et al.  Cosmology intertwined III: f$\sigma$8 and S8.  _Astroparticle Physics_ , 131:102604, sep 2021.  doi: 10.1016/j.astropartphys.2021.102604.
  * Hu and Wang [2023] Jian-Ping Hu and Fa-Yin Wang.  Hubble tension: The evidence of new physics, 2023.
  * Clifton et al. [2012] Timothy Clifton, Pedro G. Ferreira, Antonio Padilla, and Constantinos Skordis.  Modified gravity and cosmology.  _Physics Reports_ , 513(1-3):1–189, mar 2012\.  doi: 10.1016/j.physrep.2012.01.001.  URL https://doi.org/10.1016%2Fj.physrep.2012.01.001.
  * Lovelock [1971] D. Lovelock.  The Einstein tensor and its generalizations.  _J. Math. Phys._ , 12:498–501, 1971.  doi: 10.1063/1.1665613.
  * Lovelock [1972] D. Lovelock.  The four-dimensionality of space and the einstein tensor.  _J. Math. Phys._ , 13:874–876, 1972.  doi: 10.1063/1.1666069.
  * Lanczos [1938] Cornelius Lanczos.  A Remarkable property of the Riemann-Christoffel tensor in four dimensions.  _Annals Math._ , 39:842–850, 1938.  doi: 10.2307/1968467.
  * Glavan and Lin [2020] Draž en Glavan and Chunshan Lin.  Einstein-gauss-bonnet gravity in four-dimensional spacetime.  _Physical Review Letters_ , 124(8), feb 2020.  doi: 10.1103/physrevlett.124.081301.  URL https://doi.org/10.1103%2Fphysrevlett.124.081301.
  * Gürses et al. [2020] Metin Gürses, Tahsin Çağrı Şişman, and Bayram Tekin.  Is there a novel Einstein–Gauss–Bonnet theory in four dimensions?  _Eur. Phys. J. C_ , 80(7):647, 2020.  doi: 10.1140/epjc/s10052-020-8200-7.
  * Ai [2020] Wen-Yuan Ai.  A note on the novel 4D Einstein–Gauss–Bonnet gravity.  _Commun. Theor. Phys._ , 72(9):095402, 2020.  doi: 10.1088/1572-9494/aba242.
  * Shu [2020] Fu-Wen Shu.  Vacua in novel 4D Einstein-Gauss-Bonnet Gravity: pathology and instability?  _Phys. Lett. B_ , 811:135907, 2020.  doi: 10.1016/j.physletb.2020.135907.
  * Clifton et al. [2020] Timothy Clifton, Pedro Carrilho, Pedro G. S. Fernandes, and David J. Mulryne.  Observational constraints on the regularized 4d einstein-gauss-bonnet theory of gravity.  _Physical Review D_ , 102(8), oct 2020.  doi: 10.1103/physrevd.102.084005.  URL https://doi.org/10.1103%2Fphysrevd.102.084005.
  * Hennigar et al. [2020] Robie A. Hennigar, David KubizÅák, Robert B. Mann, and Christopher Pollack.  On taking the D → 4 limit of Gauss-Bonnet gravity: theory and solutions.  _Journal of High Energy Physics_ , 2020(7):27, July 2020.  doi: 10.1007/JHEP07(2020)027.
  * Fernandes et al. [2020] Pedro G. S. Fernandes, Pedro Carrilho, Timothy Clifton, and David J. Mulryne.  Derivation of Regularized Field Equations for the Einstein-Gauss-Bonnet Theory in Four Dimensions.  _Phys. Rev. D_ , 102(2):024025, 2020.  doi: 10.1103/PhysRevD.102.024025.
  * Horndeski [1974] Gregory Walter Horndeski.  Second-order scalar-tensor field equations in a four-dimensional space.  _Int. J. Theor. Phys._ , 10:363–384, 1974.  doi: 10.1007/BF01807638.
  * Lü and Pang [2020] H. Lü and Yi Pang.  Horndeski gravity as $d\rightarrow 4$ limit of gauss-bonnet.  _Physics Letters B_ , 809:135717, oct 2020.  doi: 10.1016/j.physletb.2020.135717.  URL https://doi.org/10.1016%2Fj.physletb.2020.135717.
  * Fernandes et al. [2022] Pedro G S Fernandes, Pedro Carrilho, Timothy Clifton, and David J Mulryne.  The 4d einstein–gauss–bonnet theory of gravity: a review.  _Classical and Quantum Gravity_ , 39(6):063001, feb 2022.  doi: 10.1088/1361-6382/ac500a.  URL https://doi.org/10.1088%2F1361-638%2Fac500a.
  * Banerjee et al. [2021a] Ayan Banerjee, Takol Tangphati, and Phongpichit Channuie.  Strange Quark Stars in 4D Einstein-Gauss-Bonnet Gravity.  _Astrophys. J._ , 909(1):14, 2021a.  doi: 10.3847/1538-4357/abd094.
  * Banerjee et al. [2021b] Ayan Banerjee, Takol Tangphati, Daris Samart, and Phongpichit Channuie.  Quark Stars in 4D Einstein–Gauss–Bonnet Gravity with an Interacting Quark Equation of State.  _Astrophys. J._ , 906(2):114, 2021b.  doi: 10.3847/1538-4357/abc87f.
  * Tangphati et al. [2021a] Takol Tangphati, Anirudh Pradhan, Abdelghani Errehymy, and Ayan Banerjee.  Quark stars in the Einstein–Gauss–Bonnet theory: A new branch of stellar configurations.  _Annals Phys._ , 430:168498, 2021a.  doi: 10.1016/j.aop.2021.168498.
  * Tangphati et al. [2021b] Takol Tangphati, Anirudh Pradhan, Abdelghani Errehymy, and Ayan Banerjee.  Anisotropic quark stars in Einstein-Gauss-Bonnet theory.  _Phys. Lett. B_ , 819:136423, 2021b.  doi: 10.1016/j.physletb.2021.136423.
  * Pretel et al. [2022] Juan M. Z. Pretel, Ayan Banerjee, and Anirudh Pradhan.  Electrically charged quark stars in 4D Einstein–Gauss–Bonnet gravity.  _Eur. Phys. J. C_ , 82(2):180, 2022.  doi: 10.1140/epjc/s10052-022-10123-4.
  * Feng et al. [2021] Jia-Xi Feng, Bao-Min Gu, and Fu-Wen Shu.  Theoretical and observational constraints on regularized 4d einstein-gauss-bonnet gravity.  _Physical Review D_ , 103(6), mar 2021.  doi: 10.1103/physrevd.103.064002.  URL https://doi.org/10.1103%2Fphysrevd.103.064002.
  * Haghani [2020] Zahra Haghani.  Growth of matter density perturbations in 4d einstein–gauss–bonnet gravity.  _Physics of the Dark Universe_ , 30:100720, dec 2020.  doi: 10.1016/j.dark.2020.100720.  URL https://doi.org/10.1016%2Fj.dark.2020.100720.
  * Wang and Mota [2021] Deng Wang and David Mota.  4D Gauss–Bonnet gravity: Cosmological constraints, $H_{0}$ tension and large scale structure.  _Physics of the Dark Universe_ , 32:100813, may 2021.  doi: 10.1016/j.dark.2021.100813.
  * Charmousis et al. [2022] C. Charmousis, A. Lehé bel, E. Smyrniotis, and N. Stergioulas.  Astrophysical constraints on compact objects in 4d einstein-gauss-bonnet gravity.  _Journal of Cosmology and Astroparticle Physics_ , 2022(02):033, feb 2022.  doi: 10.1088/1475-7516/2022/02/033.  URL https://doi.org/10.1088%2F1475-7516%2F2022%2F02%2F033.
  * Aoki et al. [2021] Katsuki Aoki, Mohammad Ali Gorji, and Shinji Mukohyama.  Cosmology and gravitational waves in consistent $d\to 4$ einstein-gauss-bonnet gravity, 2021.
  * [28] Michael Gammon, Sarah Rourke, and Robert B. Mann.  Quark Stars in 4DEGB.  _(work in preparation)_.
  * Mann and Ross [1993] R. B. Mann and S. F. Ross.  The $D\to 2$ limit of general relativity.  _Classical and Quantum Gravity_ , 10(7):1405–1408, July 1993.  doi: 10.1088/0264-9381/10/7/015.
  * Kobayashi [2020] Tsutomu Kobayashi.  Effective scalar-tensor description of regularized Lovelock gravity in four dimensions.  _JCAP_ , 07:013, 2020.  doi: 10.1088/1475-7516/2020/07/013.
  * Malik and Wands [2009] Karim A. Malik and David Wands.  Cosmological perturbations.  _Physics Reports_ , 475(1-4):1–51, may 2009.  doi: 10.1016/j.physrep.2009.03.001.  URL https://doi.org/10.1016%2Fj.physrep.2009.03.001.
  * eBOSS Collaboration [2021] eBOSS Collaboration.  Completed SDSS-IV extended baryon oscillation spectroscopic survey: Cosmological implications from two decades of spectroscopic surveys at the apache point observatory.  _Physical Review D_ , 103(8), apr 2021.  doi: 10.1103/physrevd.103.083533.  URL https://doi.org/10.1103%2Fphysrevd.103.083533.
  * Zhai and Wang [2019] Zhongxu Zhai and Yun Wang.  Robust and model-independent cosmological constraints from distance measurements.  _Journal of Cosmology and Astroparticle Physics_ , 2019(07):005–005, jul 2019.  doi: 10.1088/1475-7516/2019/07/005.  URL https://doi.org/10.1088%2F1475-7516%2F2019%2F07%2F005.
  * Collaboration [2020] Planck Collaboration.  Planck 2018 results. $\text{VI}.$ cosmological parameters.  _Astronomy & Astrophysics_, 641:A6, sep 2020.  doi: 10.1051/0004-6361/201833910.  URL https://doi.org/10.1051%2F0004-6361%2F201833910.
  * García-Aspeitia and Hernández-Almada [2021] Miguel A. García-Aspeitia and A. Hernández-Almada.  Einstein–gauss–bonnet gravity: Is it compatible with modern cosmology?  _Physics of the Dark Universe_ , 32:100799, 2021.  ISSN 2212-6864.  doi: https://doi.org/10.1016/j.dark.2021.100799.  URL https://www.sciencedirect.com/science/article/pii/S2212686421000303.
  * Steigman [2007] Gary Steigman.  Primordial Nucleosynthesis in the Precision Cosmology Era.  _Annual Review of Nuclear and Particle Science_ , 57(1):463–491, nov 2007.  doi: 10.1146/annurev.nucl.56.080805.140437.  URL https://doi.org/10.1146%2Fannurev.nucl.56.080805.140437.
  * Schöneberg et al. [2019] Nils Schöneberg, Julien Lesgourgues, and Deanna C. Hooper.  The BAO$+$BBN take on the hubble tension.  _Journal of Cosmology and Astroparticle Physics_ , 2019(10):029–029, oct 2019.  doi: 10.1088/1475-7516/2019/10/029.  URL https://doi.org/10.1088%2F1475-7516%2F2019%2F10%2F029.
  * Riemer-Sørensen et al. [2013] S. Riemer-Sørensen, D. Parkinson, and T. M. Davis.  What is half a neutrino? reviewing cosmological constraints on neutrinos and dark radiation.  _Publications of the Astronomical Society of Australia_ , 30, 2013\.  doi: 10.1017/pas.2013.005.  URL https://doi.org/10.1017%2Fpas.2013.005.
  * Collaboration [2023] DES Collaboration.  Dark energy survey year 3 results: Constraints on extensions to $\mathrm{\Lambda}\mathrm{CDM}$ with weak lensing and galaxy clustering.  _Phys. Rev. D_ , 107:083504, Apr 2023.  doi: 10.1103/PhysRevD.107.083504.  URL https://link.aps.org/doi/10.1103/PhysRevD.107.083504.
  * Thornton et al. [2016] RJ Thornton, PAR Ade, S Aiola, FE Angile, M Amiri, JA Beall, DT Becker, HM Cho, SK Choi, P Corlies, et al.  The atacama cosmology telescope: the polarization-sensitive actpol instrument.  _The Astrophysical Journal Supplement Series_ , 227(2):21, dec 2016.  doi: 10.3847/1538-4365/227/2/21.  URL https://dx.doi.org/10.3847/1538-4365/227/2/21.
  * Sugiyama [2014] Naoshi Sugiyama.  Introduction to temperature anisotropies of Cosmic Microwave Background radiation.  _Progress of Theoretical and Experimental Physics_ , 2014(6):06B101, 06 2014.  ISSN 2050-3911.  doi: 10.1093/ptep/ptu073.  URL https://doi.org/10.1093/ptep/ptu073.
  * Aiola et al. [2020] Simone Aiola, Erminia Calabrese, Loïc Maurin, Sigurd Naess, Benjamin L. Schmitt, Maximilian H. Abitbol, Graeme E. Addison, Peter A. R. Ade, David Alonso, Mandana Amiri, et al.  The atacama cosmology telescope: DR4 maps and cosmological parameters.  _Journal of Cosmology and Astroparticle Physics_ , 2020(12):047–047, dec 2020.  doi: 10.1088/1475-7516/2020/12/047.  URL https://doi.org/10.1088%2F1475-7516%2F2020%2F12%2F047.
  * Callin [2006] Petter Callin.  How to calculate the CMB spectrum.  _arXiv e-prints_ , art. astro-ph/0606683, June 2006.  doi: 10.48550/arXiv.astro-ph/0606683.  URL https://ui.adsabs.harvard.edu/abs/2006astro.ph..6683C.
  * Lewis et al. [2000] Antony Lewis, Anthony Challinor, and Anthony Lasenby.  Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models.  _The Astrophysical Journal_ , 538(2):473, aug 2000\.  doi: 10.1086/309179.  URL https://dx.doi.org/10.1086/309179.
  * Blas et al. [2011] Diego Blas, Julien Lesgourgues, and Thomas Tram.  The cosmic linear anisotropy solving system (CLASS). part II: Approximation schemes.  _Journal of Cosmology and Astroparticle Physics_ , 2011(07):034–034, jul 2011.  doi: 10.1088/1475-7516/2011/07/034.  URL https://doi.org/10.1088%2F1475-7516%2F2011%2F07%2F034.
  * Collaboration [2022] Euclid Collaboration.  Euclid preparation. I. The Euclid Wide Survey.  _Astronomy & Astrophysics_, 662:A112, June 2022.  doi: 10.1051/0004-6361/202141938.  URL https://ui.adsabs.harvard.edu/abs/2022A&A...662A.112E.
  * Ivezić et al. [2019] Željko Ivezić et al.  LSST: From Science Drivers to Reference Design and Anticipated Data Products.  _The Astrophysical Journal_ , 873(2):111, 2019\.  doi: 10.3847/1538-4357/ab042c.  URL https://dx.doi.org/10.3847/1538-4357/ab042c.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.  URL https://rdcu.be/b08Wh.
  * Harris et al. [2020] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, et al.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.  doi: 10.1038/s41586-020-2649-2.  URL https://doi.org/10.1038/s41586-020-2649-2.
  * Chisari et al. [2019] Nora Elisa Chisari, David Alonso, Elisabeth Krause, C. Danielle Leonard, Philip Bull, Jérémy Neveu, Antonio Villarreal, Sukhdeep Singh, Thomas McClintock, John Ellison, Zilong Du, Joe Zuntz, Alexander Mead, Shahab Joudaki, Christiane S. Lorenz, Tilman Tröster, Javier Sanchez, Francois Lanusse, Mustapha Ishak, Renée Hlozek, Jonathan Blazek, Jean-Eric Campagne, et al.  Core cosmology library: Precision cosmological predictions for LSST.  _The Astrophysical Journal Supplement Series_ , 242(1):2, may 2019.  doi: 10.3847/1538-4365/ab1658.  URL https://doi.org/10.3847%2F1538-4365%2Fab1658.
  * Foreman-Mackey [2016] Daniel Foreman-Mackey.  corner.py: Scatterplot matrices in python.  _The Journal of Open Source Software_ , 1(2):24, jun 2016.  doi: 10.21105/joss.00024.  URL https://doi.org/10.21105/joss.00024.
  * Dodelson and Schmidt [2020] S. Dodelson and F. Schmidt.  _Modern Cosmology_.  Elsevier Science, 2020.  ISBN 9780128159484.  URL https://books.google.co.uk/books?id=GGjfywEACAAJ.

##  Appendix A Equations of Motion for $k\neq 0$ spatial sections

Here we will present the gravitational and scalar field equations for non-flat
($k\neq 0$) spatial curvature where we compute in conformal time with the line
element outlined by (2.17). While some of the presentations will be redundant,
we nonetheless will proceed with clarity in our calculations. The background
calculations are carried out in Section 2.2.1 so we will not repeat them here,
and instead just note that the solution to (2.6) is given as

$\bar{\phi}^{\prime}=\mathcal{H}+A\Rightarrow\bar{\phi}=\ln(a)+A\eta+B$ (A.1)

which used in conjunction with (2.5) yields the finalized gravitational field
equations

$\displaystyle\frac{\left(\mathcal{H}^{2}+k\right)^{2}\alpha}{a^{2}}+\mathcal{H}^{2}+k=\frac{8\pi
Ga^{2}\bar{\rho}}{3}+\frac{\alpha C^{2}}{a^{2}}$ (A.2)
$\displaystyle\frac{\left(\mathcal{H}^{2}+k\right)\left(4\mathcal{H}^{\prime}-\mathcal{H}^{2}-k\right)\alpha}{a^{2}}+\mathcal{H}^{2}+2\mathcal{H}^{\prime}+k=-8\pi
Ga^{2}\bar{p}-\frac{\alpha C^{2}}{a^{2}}$

where we have made the substitution $A^{2}=-k+C$. The conservation of stress-
energy $\nabla_{\mu}T^{\mu\nu}=0$ yields only one equation for the background,
that being

$\bar{\rho}^{\prime}=-3\mathcal{H}(\bar{p}+\bar{\rho})$ (A.3)

For the first order equations we will present first with the arbitrary
constant$A$ and then make the substitution $A^{2}=-k+C$, placing a box around
the finalized equation. Beginning with the scalar field equation,which is
found in Section 2.2.2

$\alpha\left(A^{2}+k\right)\left(\nabla^{2}(\delta\phi-\Phi)+3A\left(\Psi^{\prime}+\Phi^{\prime}\right)-3\left(\Psi^{\prime\prime}+\delta\phi^{\prime\prime}\right)\right)=0$
(A.4) $\boxed{\alpha
C\left(\nabla^{2}(\delta\phi-\Phi)+3\sqrt{-k+C}\left(\Psi^{\prime}+\Phi^{\prime}\right)-3\left(\Psi^{\prime\prime}+\delta\phi^{\prime\prime}\right)\right)=0.}$
(A.5)

For the gravitational field equations, we will consider the $\eta\eta$
component first. The stress energy tensor remains the same $\delta
T_{\eta\eta}=a^{2}\delta\rho+2a^{2}\Phi\bar{\rho}$. Computing
$\delta\mathcal{E}_{\eta\eta}=8\pi G(a^{2}\delta\rho+2a^{2}\Phi\bar{\rho})$
and using the first equation of (A.2) we arrive at the simplified expression

$\displaystyle 4\alpha
A^{2}\left[3A^{2}\Phi-3A(\delta\phi^{\prime}+\Psi^{\prime})-\nabla^{2}(\delta\phi+\Psi)-3k(\Psi-\Phi)\right]-4\alpha
k\left[3A(\delta\phi^{\prime}+\Psi^{\prime})+\nabla^{2}\delta\phi\right]+2\mathcal{A}\nabla^{2}\Psi$
(A.6) $\displaystyle-6\mathcal{H}(\mathcal{A}+2\alpha
k)\Psi^{\prime}-12\alpha\mathcal{H}^{2}\left((k+\mathcal{H}^{2})\Phi-k\Psi\right)+6a^{2}(k\Psi-\mathcal{H}^{2}\Phi)=8\pi
Ga^{4}\delta\rho.$ $\mathclap{\boxed{\begin{aligned}
&2(\mathcal{A}+2k\alpha)(\nabla^{2}\Psi-3\mathcal{H}\Psi^{\prime}+3(k\Psi-\Phi\mathcal{H}^{2}))+4\alpha
C(3\Phi
C-\nabla^{2}(\Psi+\delta\phi)-3\sqrt{-k+C}(\delta\phi^{\prime}+\Psi^{\prime}))=8\pi
a^{4}G\delta\rho.\end{aligned}}}$

Any one of the time-spatial equations $\delta\mathcal{E}_{0i}=8\pi G\delta
T_{0i}$ has no corresponding zeroth order equation so calculating this
directly yields

$\mathclap{4\alpha
A^{2}\left[A(\Phi_{,i}-\delta\phi_{,i})-(\delta\phi^{\prime}_{,i}+\psi^{\prime}_{,i})\right]+4\alpha
kA(\Phi_{,i}-\delta\phi_{,i})+2\mathcal{H}(\mathcal{A}+2\alpha
k)\Phi_{,i}+2\mathcal{A}\Psi^{\prime}_{,i}-4\alpha
k\delta\phi^{\prime}_{,i}=-8\pi Ga^{4}(1+w)\bar{\rho}v_{,i}}$
$\boxed{2\mathcal{H}(\mathcal{A}+2\alpha k)(\Phi+\Psi^{\prime})+4\alpha
C(\sqrt{-k+C}(\Phi-\delta\phi)-\delta\phi^{\prime}-\Psi^{\prime})=-8\pi
G(1+w)a^{4}\bar{\rho}v.}$ (A.7)

For a mixed spatial-spatial component $\delta\mathcal{E}_{ij}=8\pi G\delta
T_{ij}$ for $i\neq j$ we have

$\nabla_{i}\partial_{j}(2\alpha
A^{2}(\Phi+\Psi)+4\alpha\mathcal{H}^{\prime}\Psi-\mathcal{A}\Phi-\mathcal{B}\Psi)=8\pi
Ga^{4}\nabla_{i}(\partial_{j}\Pi)$ (A.8)
$\boxed{2\alpha(-k+C)(\Phi+\Psi)+4\alpha\mathcal{H}^{\prime}\Psi-\mathcal{A}\Phi-\mathcal{B}\Psi=8\pi
Ga^{4}\Pi}.$ (A.9)

For the diagonal spatial components we employ the same procedure as with the
flat case. We consider all the components individually and compute the trace,
$\delta\mathcal{E}_{rr}+\delta\mathcal{E}_{\theta\theta}+\delta\mathcal{E}_{\phi\phi}=8\pi
G(\delta T_{rr}+\delta T_{\theta\theta}+\delta T_{\phi\phi})$; using the 2nd
equation of (A.2), this gives

$\displaystyle 4\alpha
A^{2}\left(3A^{2}\Phi+3A(\Phi^{\prime}-\delta\phi^{\prime})-\nabla^{2}(\Phi+\Psi)-3(\Psi^{\prime\prime}+\delta\phi^{\prime\prime})+3k(\Phi-\Psi)\right)+12k\alpha
A(\Phi^{\prime}-\delta\phi^{\prime})$ (A.10)
$\displaystyle+4\mathcal{H}^{\prime}\left(6\mathcal{H}\alpha\Psi^{\prime}+3(4\mathcal{H}^{2}+a^{2}+2k\alpha)\Phi-6k\alpha\Psi-2\nabla^{2}\Psi\alpha\right)+2\mathcal{A}\nabla^{2}\Phi+6\mathcal{H}(\mathcal{A}+2k\alpha)\Phi^{\prime}+2\mathcal{B}\nabla^{2}\Psi$
$\displaystyle-12\alpha
k\delta\phi^{\prime\prime}+6\mathcal{A}\Psi^{\prime\prime}+12\mathcal{H}a^{2}\Psi^{\prime}-12\mathcal{H}^{2}(\Phi\mathcal{H}^{2}+k(\Phi-\Psi))+6a^{2}(\mathcal{H}^{2}\Phi-k\Psi)=24\pi
Ga^{4}\delta p$ $\mathclap{\boxed{\begin{aligned}
&(\mathcal{B}+2k\alpha)\nabla^{2}\Psi+(\mathcal{A}+2k\alpha)(\nabla^{2}\Phi+3\Psi^{\prime\prime})+2\alpha
C\left(-\nabla^{2}\Phi-3(\delta\phi^{\prime\prime}+\Psi^{\prime\prime})-3k(\Phi+\Psi)-3\sqrt{-k+C}\delta\phi^{\prime}\right)\\\
&2\mathcal{H}^{\prime}\left(-2\alpha\nabla^{2}\Psi+12\alpha\mathcal{H}\Phi+3a^{2}\Phi+6\alpha
k\Phi-6\alpha k\Psi+6\alpha\mathcal{H}\Psi^{\prime}\right)+6\alpha
C^{2}\Phi+6\mathcal{H}a^{2}\Psi^{\prime}+3a^{2}(\mathcal{H}^{2}-k\Psi)\\\
&+3\left(2\alpha\mathcal{H}^{3}+2\alpha(-k+C)^{3/2}+\mathcal{H}a^{2}+2\mathcal{H}\alpha
k+2k\alpha\sqrt{-k+C}\right)\Phi^{\prime}-6\alpha(\mathcal{H}^{2}+k)(\mathcal{H}^{2}\Phi-\Psi
k)=12\pi Ga^{4}\delta p.\end{aligned}}}$

We now need to consider the conservation of stress energy,
$\nabla_{\mu}T^{\mu\nu}=0$. Considering the time component first
$\nabla_{\mu}T^{\mu\eta}=0$, we obtain the same zeroth and first order
equations as in the flat case; hence

$\boxed{\delta\rho^{\prime}=(1+w)\left(3\Psi^{\prime}-\nabla^{2}v\right)\bar{\rho}-3\mathcal{H}(\delta\rho+\delta
p)}\;.$ (A.11)

For the spatial components there are some slight modifications. For any
spatial coordinate $i$, the equation $\nabla_{\mu}T^{\mu i}=0$ becomes

$\partial_{i}\left(\nabla^{2}\Pi+(3w+3)\bar{\rho}v^{\prime}-3(1+w)(3w-1)\mathcal{H}\bar{\rho}v+3(1+w)\bar{\rho}\Phi+3\delta
p+6k\Pi\right)=0.$ (A.12)

Taking the spatial divergence of this yields

$\nabla^{i}\left(\partial_{i}\left(\nabla^{2}\Pi+(3w+3)\bar{\rho}v^{\prime}-3(1+w)(3w-1)\mathcal{H}\bar{\rho}v+3(1+w)\bar{\rho}\Phi+3\delta
p+6k\Pi\right)\right)=0$ (A.13)
$\nabla^{2}\left(\nabla^{2}\Pi+(3w+3)\bar{\rho}v^{\prime}-3(1+w)(3w-1)\mathcal{H}\bar{\rho}v+3(1+w)\bar{\rho}\Phi+3\delta
p+6k\Pi\right)=0$ (A.14)
$\boxed{\nabla^{4}\Pi+3(1+w)\bar{\rho}\nabla^{2}v^{\prime}-3(1+w)(3w-1)\mathcal{H}\bar{\rho}\nabla^{2}v+3(1+w)\bar{\rho}\nabla^{2}\Phi+3\nabla^{2}\delta
p+6k\nabla^{2}\Pi=0}$ (A.15)

##  Appendix B Scalar Perturbations in the small-$\alpha$ limit

We start with equations (2.43). We derive an approximation for these scalar
perturbations in the following limits:

  1. 1.

$\alpha$ is small, meaning $\alpha/H_{0}^{2}<<1$. This is a reasonable
approximation given [12] gives us priors
$-10^{-83}\lesssim\tilde{\alpha}\lesssim 10^{-43}$. We can then write
$\mathcal{A}=\mathcal{C}=\mathcal{D}=a^{2}$ and $\mathcal{B}=-a^{2}$.

  2. 2.

$C$ is large. Specifically, we assume large enough that $\alpha_{C}=\alpha
C^{2}$ is finite and affects late universe dynamics, but small enough that
$\alpha C^{1/2}$, $\alpha C$ and $\alpha C^{3/2}$ terms can be neglected.

  3. 3.

$\alpha$ is small enough that the background equation reduces to
$H^{2}=\frac{8\pi G\rho}{3}+\frac{\alpha_{C}}{a^{4}}$.

Using these approximations we can rewrite equation (2.43a) as:

$k^{2}\delta\phi+3\delta\phi^{\prime\prime}=3\sqrt{C}(\Phi^{\prime}+\Psi^{\prime})$
(B.1)

We can do this because we know $\nabla^{2}\Phi<<\sqrt{C}\Phi^{\prime}$ and
$\Psi^{\prime\prime}<<\sqrt{C}\Psi^{\prime}$. This suggests that
$\delta\varphi=\frac{\delta\phi}{\sqrt{C}}$ is the same order of magnitude as
the metric perturbations $\Phi$ and $\Psi$, so that we can re-write our
equation as:

$k^{2}\delta\varphi+3\delta\varphi^{\prime\prime}=3(\Phi^{\prime}+\Psi^{\prime})$
(B.2)

Then using the same approximations we can rewrite equation (2.43b) as:

$12\alpha_{C}\Phi-12\alpha_{C}\delta\varphi^{\prime}-2a^{2}(k^{2}\Psi+3\mathcal{H}\Psi^{\prime}+3\mathcal{H}^{2}\Phi)=8\pi
Ga^{4}\delta\rho$ (B.3)

The remaining equations reduce to their GR equivalent:

$\displaystyle\Psi^{\prime}+\mathcal{H}\Phi=-4\pi Ga^{2}(1+\omega)\rho v$
(B.4) $\displaystyle\Phi=\Psi+8\pi Ga^{2}\Pi$ $\displaystyle
3\Psi^{\prime\prime}+3\mathcal{H}(\Phi^{\prime}+2\Psi^{\prime})+k^{2}(\Psi-\Phi)+(6\mathcal{H}^{\prime}-3\mathcal{H}^{2})\Phi=12\pi
Ga^{2}\delta\rho$
$\displaystyle\bm{\delta}^{\prime}=3\mathcal{H}\bm{\delta}\left(w-\frac{\delta
p}{\delta\rho}\right)+(1+w)\left(3\Psi^{\prime}-\nabla^{2}v\right)$
$\displaystyle 2\nabla^{4}\Pi+3\nabla^{2}\delta
p+3(1+w)\bar{\rho}\left[(1-3w)\nabla^{2}v+\nabla^{2}v^{\prime}+\nabla^{2}\Phi\right]=0$

##  Appendix C Coupled Differential Equations for the Boltzmann Solver

The following derivations follow the procedures outlined in [43, 52].

The Boltzmann equations for photons, photon polarization, neutrinos, dark
matter and baryons are independent of the cosmological model (they describe
the equations in any given expanding spacetime with perturbations). They will
therefore be the same in $\Lambda$CDM and 4DEGB. These can be written as (see
equation (22) in [43]):

$\begin{split}&\Theta^{\prime}+ik\mu\Theta=\Psi^{\prime}-ik\Phi-\tau^{\prime}\left[\Theta_{0}-\Theta+ik\mu
v_{b}-\frac{1}{2}\mathcal{P}_{2}\Pi\right]\\\
&\Theta_{P}^{\prime}+ik\mu\Theta_{P}=-\tau^{\prime}\left[\Theta_{P}+\frac{1}{2}(1-\mathcal{P}_{2})\Pi\right]\\\
&\delta^{\prime}-k^{2}v=3\Psi^{\prime}\\\ &v^{\prime}+\mathcal{H}v=-\Phi\\\
&\delta_{b}^{\prime}-k^{2}v_{b}=3\Psi^{\prime}\\\
&v_{b}^{\prime}+\mathcal{H}v_{b}=-\Phi+\frac{\tau^{\prime}R}{k}(v_{b}+3\Theta_{1})\\\
&\mathcal{N}^{\prime}+ik\mu\mathcal{N}=\Psi^{\prime}-ik\mu\Phi\end{split}$
(C.1)

where all quantities are defined as in [43], except for
$\Phi_{\textrm{\cite[cite]{[\@@bibref{Number}{callin2006calculate}{}{}]}}}=-\Psi$,
$\Psi_{\textrm{\cite[cite]{[\@@bibref{Number}{callin2006calculate}{}{}]}}}=\Phi$
and
$v_{\textrm{\cite[cite]{[\@@bibref{Number}{callin2006calculate}{}{}]}}}=kv$.
Care should also be taken when comparing derivatives (the prime symbol in [43]
indicates a derivative with respect to $x=\textrm{ln}a$). To explicitly solve
the coupled ODEs we need to use the multipole expansion of the first two
equations in (C.1).

Next, by looking at (B.3) and the second equation in (B.4) (when including
photon and neutrino perturbations) we find that the relevant Einstein
equations take the form:

$\begin{split}&k^{2}\Psi+3\mathcal{H}(\Psi^{\prime}+\mathcal{H}\Phi)=\frac{6\alpha_{C}(\Phi-\delta\varphi^{\prime})}{a^{2}}+4\pi
Ga^{2}\left[\rho\delta+\rho_{b}\delta_{b}+4\rho_{\gamma}\Theta_{0}+4\rho_{\nu}\mathcal{N}_{0}\right]\\\
&\Phi=\Psi-\frac{12H_{0}^{2}}{k^{2}a^{2}}\left[\Omega_{\gamma
0}\Theta_{2}+\Omega_{\nu 0}\mathcal{N}_{2}\right]\end{split}$ (C.2)

We also need an additional differential equation for $\delta\varphi$ which is
given by (B.2):

$k^{2}\delta\varphi+3\varphi^{\prime\prime}=3(\Phi^{\prime}+\Psi^{\prime})$
(C.3)

###  C.1 Initial conditions for the coupled ODEs

In the early universe for small $k$, the first equation in (C.2) reduces to:

$\begin{split}&\frac{\Psi^{\prime}}{\eta}+\frac{\Phi}{\eta^{2}}-\frac{2\tilde{\alpha}_{C}}{(\Omega_{r0}+\tilde{\alpha}_{C})\eta^{2}}(\Phi-\delta\varphi^{\prime})=-2\frac{8\pi
Ga^{2}\rho}{3}(\frac{\rho_{\gamma 0}}{\rho}\Theta_{0}+\frac{\rho_{\nu
0}}{\rho}\mathcal{N}_{0})\implies\\\
&\frac{\Psi^{\prime}}{\eta}+\frac{\Phi}{\eta^{2}}-\frac{2\tilde{\alpha}_{C}}{(\Omega_{r0}+\tilde{\alpha}_{C})\eta^{2}}(\Phi-\delta\varphi^{\prime})=-\frac{2}{\eta^{2}}(\frac{\Omega_{r0}}{\Omega_{r0}+\tilde{\alpha}_{C}})((1-f_{\nu})\Theta_{0}+f_{\nu}\mathcal{N}_{0})\end{split}$
(C.4)

where we have used the background equation $\frac{8\pi
Ga^{2}\rho}{3}=\mathcal{H}^{2}-\frac{H_{0}^{2}\tilde{\alpha}_{C}}{a^{2}}$, so
that $\mathcal{H}=\frac{1}{\eta}$ and
$a=H_{0}\sqrt{\Omega_{r0}+\tilde{\alpha}_{C}}\eta$, and where
$f_{\nu}=\frac{\rho_{\nu}}{\rho_{r}+\rho_{\nu}}=\frac{1}{\frac{8}{7N_{\textrm{eff}}}(\frac{11}{4})^{4/3}+1}$
when all neutrinos are massless. Now we use the fact that for small $k$,
$\Psi^{\prime}=\Theta_{0}^{\prime}=\mathcal{N}_{0}^{\prime}$ together with the
derivative of (C.4)$\times\eta^{2}$ to get:

$\eta\Psi^{\prime\prime}+\Psi^{\prime}+(1-\frac{2\tilde{\alpha}_{C}}{\tilde{\alpha}_{C}+\Omega_{r0}})\Phi^{\prime}+\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}+\tilde{\alpha}_{C}}\delta\varphi^{\prime\prime}=-\frac{2\Omega_{r0}}{\Omega_{r0}+\tilde{\alpha}_{C}}\Psi^{\prime}$
(C.5)

Next we employ the ansatz $\Psi=(1+\epsilon)\Phi$ in the early universe, where
$\epsilon$ is a constant to be determined. Then, using (C.2) for small $k$
(i.e. $\delta\varphi^{\prime\prime}=(2+\epsilon)\Phi^{\prime}$), we get:

$\begin{split}&(1+\epsilon)\eta\Phi^{\prime\prime}+\left(2+\epsilon-\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}+\tilde{\alpha}_{C}}+\frac{2\Omega_{r0}}{\Omega_{r0}+\tilde{\alpha}_{C}}(1+\epsilon)+\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}+\tilde{\alpha}_{C}}(2+\epsilon)\right)\Phi^{\prime}=0\implies\\\
&(1+\epsilon)\eta\Phi^{\prime\prime}+(4+3\epsilon)\Phi^{\prime}=0\end{split}$
(C.6)

This equation is solved by $\Phi=\eta^{P}$ for
$P=0,1-\frac{4+3\epsilon}{1+\epsilon}$, which are respectively a growing and a
decaying mode for $\epsilon>-1$.This means the early-universe growing mode has
solution $\Phi=\Phi_{i}$ (defined at some early time $\eta_{i}$) which is a
constant. This implies $\delta\varphi^{\prime}=\delta\varphi^{\prime}_{i}$ is
also a constant in the early universe (at some early time $\eta_{i}$).

We can now find the initial conditions for other quantities as a function of
$\Phi_{i}$ and $\delta\varphi^{\prime}_{i}$ both. We start with the bottom
line of equation (C.4) with $\Psi^{\prime}=0$ and
$\Theta_{0,i}=\mathcal{N}_{0,i}$ to get:

$\begin{split}&(1-\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}+\tilde{\alpha}_{C}})\Phi_{i}+\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}+\tilde{\alpha}_{C}}\delta\varphi^{\prime}_{i}=-\frac{2\Omega_{r0}}{\Omega_{r,0}+\tilde{\alpha}_{C}}\Theta_{0,i}\implies\\\
&\boxed{\Theta_{0,i}=-\frac{1}{2}\Big{\\{}\frac{\Omega_{r0}-\tilde{\alpha}_{C}}{\Omega_{r0}}\Phi_{i}+\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}}\delta\varphi^{\prime}_{i}\Big{\\}}}\end{split}$
(C.7)

Next, we use the third equation in (C.1) as well as the multipole expansion of
the first equation in (C.1) (see the first equation in (22) from [43]) with
adiabatic initial conditions ($\frac{\delta}{\Theta_{0}}$ is constant) to get:

$\boxed{\delta_{i}=3\Theta_{0,i}}$ (C.8)

We also use the first equation in (22) from [43] with
$\Phi^{\prime}=0,3\Theta_{1}+v_{b}\approx 0$ and $\Theta_{2}\ll 1$ to get:

$\Theta_{1,i}=\frac{k\eta}{3}(\Phi_{i}+\Theta_{0,i})\implies\boxed{\Theta_{1,i}=\frac{k\Phi_{i}}{6\mathcal{H}}\Big{\\{}\frac{\Omega_{r0}+\tilde{\alpha}_{C}}{\Omega_{r0}}-\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}}\frac{\delta\varphi^{\prime}_{i}}{\Phi_{i}}\Big{\\}}}$
(C.9)

The third equation in (C.1) as well as the first equation in (22) from [43]
(with adiabatic initial conditions) also give:

$-k^{2}v_{i}=(3\Phi^{\prime}-\delta^{\prime})_{i}=(3\Theta_{0}^{\prime}+3k\Theta_{1}-\delta^{\prime})_{i}\approx
3k\Theta_{1,i}\implies\boxed{\Theta_{1,i}=-\frac{kv_{i}}{3}}$ (C.10)

To find the value of $\epsilon$, we first note that the seventh and eighth
equations in (22) from [43] in the early universe imply:

$\mathcal{N}_{2}^{\prime\prime}\approx\frac{2k}{5}\mathcal{N}_{1}^{\prime}=\frac{2k^{2}}{15}(\mathcal{N}_{0}-2\mathcal{N}_{2}+\Phi)$
(C.11)

and the second equation in (C.2) with $\Theta_{2}<<\mathcal{N}_{2}$ and
$a=H_{0}\sqrt{\Omega_{r0}+\tilde{\alpha}_{C}}\eta$ implies:

$\mathcal{N}_{2}=-\frac{\Omega_{r0}+\tilde{\alpha}_{C}}{\Omega_{r0}}\frac{(k\eta)^{2}(\Phi-\Psi)}{12f_{\nu}}$
(C.12)

Differentiating (C.12) and equating to (C.11) (with
$\mathcal{N}_{2}\ll|\Phi+\mathcal{N}_{0}|$) we get:

$\mathcal{N}_{0}=\frac{5}{4f_{\nu}}\left(\frac{\Omega_{r0}+\tilde{\alpha}_{C}}{\Omega_{r0}}\right)\Psi-\left(\frac{5}{4f_{\nu}}\left(\frac{\Omega_{r0}+\tilde{\alpha}_{C}}{\Omega_{r0}}\right)+1\right)\Phi$
(C.13)

which, using
$\mathcal{N}_{0,i}=\Theta_{0,i}=-\frac{1}{2}\Big{\\{}\frac{\Omega_{r0}-\tilde{\alpha}_{C}}{\Omega_{r0}}\Phi_{i}+\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}}\delta\varphi^{\prime}_{i}\Big{\\}}$
implies:

$\begin{split}&\Psi_{i}=\frac{4f_{\nu}}{5}\left(\frac{\Omega_{r0}}{\Omega_{r0}+\tilde{\alpha}_{C}}\right)\left\\{\left(1-\frac{\Omega_{r0}-\tilde{\alpha}_{C}}{2\Omega_{r0}}+\frac{5}{4f_{\nu}}\left(\frac{\Omega_{r0}+\tilde{\alpha}_{C}}{\Omega_{r0}}\right)\right)\Phi_{i}-\frac{\tilde{\alpha}_{C}}{\Omega_{r0}}\delta\varphi^{\prime}_{i}\right\\}\implies\\\
&\boxed{\Psi_{i}=\left(\frac{2f_{\nu}}{5}+1\right)\Phi_{i}-\frac{4f_{\nu}\tilde{\alpha}_{C}}{5(\Omega_{r0}+\tilde{\alpha}_{C})}\delta\varphi^{\prime}_{i}}\end{split}$
(C.14)

This satisfies the requirements of the ansatz for sufficiently small
$\tilde{\alpha}_{C}$, $f_{\nu}$ and $\delta\varphi^{\prime}_{i}$. It also
allows us to find:

$\boxed{\mathcal{N}_{2,i}=\frac{k^{2}a^{2}}{12H_{0}^{2}\Omega_{\nu
0}\left(\frac{5}{2f_{\nu}}+1\right)}\left(\Psi_{i}+\frac{2\tilde{\alpha}_{C}}{\Omega_{r0}+\tilde{\alpha}_{C}}\delta\varphi^{\prime}_{i}\right)},\boxed{\mathcal{N}_{\ell}=\frac{k}{(2\ell+1)\mathcal{H}}\mathcal{N}_{\ell-1},\ell>2}$
(C.15)

The derivation for the moments of the photon polarization perturbations can be
found in [43]. Turning now to the scalar field, we notice that in the early
universe (C.3) is the equation for a simple harmonic oscillator with frequency
$\omega=\frac{k}{\sqrt{3}}$. This will have a general homogeneous solution:

$\delta\varphi=\delta\varphi_{i}\textrm{cos}\left(\frac{k}{\sqrt{3}}\eta\right)+\frac{\sqrt{3}\delta\varphi^{\prime}_{i}}{k}\textrm{sin}\left(\frac{k}{\sqrt{3}}\eta\right)$
(C.16)

where $\delta\varphi_{i}$ is the initial value of the scalar field at some
early time $\eta_{i}$ and $\delta\varphi_{i}^{\prime}$ is the initial value of
its derivative. When $\Phi^{\prime}\neq 0$, there will also be a particular
integral solution to equation (C.3). In the early universe, if
$\Phi^{\prime}\lesssim\mathcal{O}(\eta)$, this particular solution will be
$\delta\varphi_{i}\lesssim\mathcal{O}(\eta^{3})\implies\left(\frac{\textrm{d}\delta\varphi}{\textrm{d}(\textrm{ln}a)}\right)_{i}=\mathcal{O}(\eta_{i}^{3})\approx
0$.

In the interest of simplicity we look at the special initial condition in
which the scalar field is initially at rest at a minimum of its effective
potential, with $\delta\varphi_{i}=\delta\varphi_{i}^{\prime}=0$. However, we
acknowledge that in general this theory is not necessarily constructed to
behave like GR in the early universe. Therefore an interesting avenue for
future research could be the exploration of the effect of different initial
conditions on perturbations in 4DEGB.

We can then write down the relevant initial conditions as:

$\begin{split}&\Theta_{0,i}=-\frac{1}{2}\left(\frac{\Omega_{r0}-\tilde{\alpha}_{C}}{\Omega_{r0}}\Phi_{i}\right)\\\
&\delta_{m,i}=3\Theta_{0,i}\\\
&\Theta_{1,i}=\frac{k\Phi_{i}}{6\mathcal{H}(\eta_{i})}\left(\frac{\Omega_{r0}+\tilde{\alpha}_{C}}{\Omega_{r0}}\Phi_{i}\right)\\\
&v_{m,i}=-\frac{3}{k}\Theta_{1,i}\\\ &\mathcal{N}_{0,i}=\Theta_{0,i}\\\
&\left(\frac{\textrm{d}\delta\varphi}{\textrm{d}(\textrm{ln}a)}\right)_{i}=\mathcal{O}(\eta_{i}^{3})\approx
0\\\
&\delta\varphi_{i}=\frac{1}{3}\left(\frac{\textrm{d}\delta\varphi}{\textrm{d}(\textrm{ln}a)}\right)_{i}\approx
0\end{split}$ (C.17)

and all other initial conditions have the same form as their GR equivalent as
described in [43].