Cosmological Fluctuations in Delta Gravity

Abstract

About 70% of the Universe is Dark Energy, but the physics community still does not know what it is. Delta gravity (DG) is an alternative theory of gravitation that could solve this cosmological problem. Previously, we studied the Universe’s accelerated expansion, where DG was able to explain the SNe-Ia data successfully. In this work, we computed the cosmological fluctuations in DG that give rise to the CMB through a hydrodynamic approximation. We calculated the gauge transformations for the metric and the perfect fluid to present the equations of the evolution of cosmological fluctuations. This provided the necessary equations to solve the scalar TT power spectrum in a semi-analytical way. These equations are useful for comparing the DG theory with astronomical observations and thus being able to constrain the DG cosmology.

1. Introduction

In the past decade, there has been a surge of interest in cosmology due to the increasingly precise observational constraints that can help elucidate the physics underlying the Universe. Despite the mounting evidence for cosmological phenomena such as the acceleration of the Universe attributed to Dark Energy (DE) and the presence of non-visible matter known as Dark Matter (DM) [1,2,3], the physics community has yet to provide a comprehensive explanation for their nature.

The standard cosmological model, known as $\Lambda$CDM, describes the composition of the Universe, where 69% of the energy density corresponds to DE, 26% corresponds to DM, and the remaining 5% is composed of ordinary matter and light [1]. This model has successfully accounted for various observations, including those of Type Ia Supernovae (SNe-Ia) and the cosmic microwave background (CMB), and it has been validated through cosmological simulations that depict the formation of large-scale structures [4,5].

However, the $\Lambda$CDM model exhibits inconsistencies between its description of the early and late Universe [6]. These inconsistencies manifest in different cosmological parameters, such as the Hubble constant [7,8], the curvature [9,10], and the $S_8$ tension [11].

The Planck team measured the local expansion rate through the cosmic microwave background (CMB) radiation and obtained a value of $H_0=67.37\pm 0.54$ km/s/Mpc, which aligns with a flat $\Lambda$CDM model [1]. However, the SH0ES collaboration independently measured a higher value of $H_0=73.52\pm 1.62$ km/s/Mpc for the local Universe [7], creating a discrepancy with the Planck value exceeding $3.5\sigma$. Importantly, this tension between the early and late Universe persists even without considering the Planck CMB data or the SH0ES distance ladder [6].

Furthermore, the H0LiCOW collaboration derived a direct measurement of $H_0=72.5_{-2.3}^{+2.1}$ km/s/Mpc based on lensing time delays, which exhibits a moderate tension with the Planck value [12]. Additionally, a constraint obtained from the Big Bang nucleosynthesis (BBN) combined with baryon acoustic oscillation (BAO) data yielded $H_0=66.98\pm 1.18$ km/s/Mpc, which is inconsistent with the SH0ES measurement [6].

Other studies have attempted to explain the discrepancy by suggesting that the Hubble constant determined from nearby SNe-Ia may differ from that measured from the CMB due to cosmic variance, with a potential difference of $\pm 0.8$ percent at $1\sigma$ statistical significance. However, this variation does not account for the observed discrepancy between SNe-Ia and CMB measurements [13]. In an extreme case, observers situated in the centers of vast cosmic voids might measure a Hubble constant biased high by 5 percent from SNe-Ia.

Since the initial publication highlighting the $H_0$ tension [2], numerous inquiries have emerged regarding the source of this discrepancy. One suggestion is that errors in the calibration of Cepheids, which contribute to systematic errors, could be responsible. However, this potential error has been thoroughly discussed and dismissed by Riess et al. [14].

Other publications have explored possible solutions to the observed acceleration of the Universe, including anisotropies at local scales. Wang et al. [15], using SNe-Ia data, found evidence of anisotropies associated with the direction and amplitude of the bulk flow. Nonetheless, the impact of dipolar distribution of dark energy cannot be ruled out at high redshifts. Similarly, another publication [16] suggests that the anisotropies in cosmic acceleration may be linked to the nature of Dark Energy, implying that the perceived cosmic acceleration deduced from supernovae could be an artifact of our non-Copernican perspective rather than evidence of a dominant “dark energy” component in the Universe. Sun et al. [17] conclude that even in the presence of anisotropy, Dark Energy cannot be entirely ruled out. While such proposals could potentially explain variations in local measurements, including different values for the local Hubble constant, they could contradict the analyses conducted by Planck using the $\Lambda$CDM model. The Dark Energy component is crucial for the evolution of CMB photons from the last scattering surface until the present, and altering the sum over $\Omega$ for each component in the Universe would lead to significant changes. Various other suggestions concerning discrepancies have emerged not only related to SNe-Ia measurements but also within the Planck data itself. The presence of anisotropies in these measurements has been a subject of debate due to high uncertainties and inconsistent results. Hypotheses proposing the possibility of a Universe with less Dark Energy [18] have also been put forward.

Another potential source of error in local measurements could be the inhomogeneity in local density [19,20]. However, in this scenario, the presence of local structures does not seem to impede the possibility of measuring the Hubble constant with a precision of 1%, and there is no evidence of a Hubble constant change corresponding to an inhomogeneity.

Today, there are different methods to obtain the Hubble constant, including the use of SNe-II, ref. [21]. In this research, SNe-II were employed as standard candles to obtain an independent measurement of the Hubble constant. The resulting value was $H_0=75.8_{-4.9}^{+5.2}$ km /s/Mpc. The local $H_0$ value is higher than the value derived from the early Universe, with a confidence level of 95%. The researchers concluded that there is no evidence that SNe-Ia are the source of the $H_0$ tension. In another publication analyzing SNe-Ia as standard candles in the near-infrared, it was concluded that $H_0$ = 72.8 ± 1.6 (statistical) ±2.7 (systematic) km/s/Mpc. This study also suggested that the tension in the competing $H_0$ distance ladders is likely not a result of supernova systematics.

Other proposals have tried to reconcile Planck and SNe-Ia data, including modifications to the physics of the DE. In other words, introducing an equation of state of the interacting dark energy component, where w is allowed to vary freely, could solve the $H_0$ tension [22]. Additionally, a decaying dark matter model has been proposed to alleviate the $H_0$ and ${\sigma }_8$ anomalies [23]; in their work, they reduce the tension for both measurements when only consider Planck CMB data and the local SH0ES prior on $H_0$. However, when BAOs and the JLA supernova dataset are included, their model is weakened.

Other disagreements are related to inconsistencies with curvature (and other parameters needed to describe the CMB) [10], or they are related to the tension between measurements of the amplitude of the power spectrum of density perturbations (inferred using CMB) and directly measured by large-scale structure (LSS) on smaller scales [11]. Extensions of $\Lambda$CDM models have been considered [24] in an attempt to solve the tension of $H_0$. However, they concluded that none of these extended models can convincingly resolve the $H_0$ tension. For a full scope of the Hubble tension, please see [25]. Through the time, the tension between Planck and SNe-Ia persists [1,14], where the $H_0$ is the most significant tension. Furthermore, the Universe is composed principally by DE, but we still do not know what it is.

Over the past decades, various proposals have been made to explain the observed acceleration of the Universe. These proposals involve the inclusion of additional fields in approaches such as Quintessence, Chameleon, Vector Dark Energy or Massive Gravity, the addition of higher-order terms in the Einstein–Hilbert action, such as $f\left(R\right)$ theories and Gauss–Bonnet terms, and the introduction of extra dimensions for a modification of gravity on large scales [26]. Other interesting possibilities include the search for non-trivial ultraviolet fixed points in gravity (asymptotic safety, [27]) and the notion of induced gravity [28,29,30,31]. The first possibility uses exact renormalization-group techniques [32,33] together with a lattice and numerical techniques, such as Lorentzian triangulation analysis [34]. Induced gravity proposes that gravitation is a residual force produced by other interactions.

Delta gravity (DG) is an extension of General Relativity (GR) where new fields are added to the Lagrangian through a new symmetry [35,36,37,38]. The main properties of this model at the classical level follow: (a) It agrees with GR outside the sources and with adequate boundary conditions. In particular, the causal structure of delta gravity in a vacuum is the same as in General Relativity, satisfying all standard tests automatically. (b) When studying the evolution of the Universe, it predicts acceleration without a cosmological constant or additional scalar fields. The Universe ends in a BigRip, which is similar to the scenario considered in [39]. (c) The scale factor agrees with the standard cosmology at early times and show acceleration only at late times. Therefore, we expect that density perturbations should not have large corrections at the moment of last scattering (denoted by $t_{ls}$).

It was noticed in [40] that the Hamiltonian of delta models is not bounded from below. Phantom cosmological models [39,41] also have this property. The present model could provide an arena to study the quantum properties of a phantom field, since the model has a finite quantum effective action. In this respect, the advantage of the present model is that being a gauge model, it could give us the possibility of solving the problem of lack of unitarity using standard techniques of gauge theories such as the BRST method [36]. However, we are not concerned about this feature in this work, because we are considering DG as a phenomenological model that interpolates the observations of the early with the late Universe.

This theory predicts an accelerating Universe without a cosmological constant $\Lambda$ and a Hubble parameter $H_0=74.47\pm 1.63$ Km/s/Mpc [42] when fitting SN-Ia Data, which is in agreement with SH0ES.

On the other hand, temperature correlations provide us with information about the constituents of the Universe, including baryonic and dark matter. Typically, these calculations are performed using software such as CMBFast [43,44] or CAMB1 [45]. These codes employ Boltzmann equations for the fluids and their interactions, yielding well-established results that are consistent with Planck measurements [1].

Nevertheless, one can obtain a good approximation of this complex problem [46,47]. In this work, we use an analytical method that consists of two steps instead of studying the evolution of the scalar perturbations using Boltzmann equations. First, we use a hydrodynamic approximation, which assumes photons and baryonic plasma as a fluid in thermal equilibrium at recombination time when there is a high rate of collisions between free electrons and photons. Second, we study the propagation of photons [35] by radial geodesics from the moment when the Universe switches from opaque to transparent at time $t_{ls}$ until now.

In this research, we develop the theory of scalar perturbations at first order. We discuss the gauge transformations in an extended Friedmann–Lemaître–Robertson–Walker (FRLW) Universe. Then, we show how to obtain an expression for temperature fluctuations, and we demonstrate that they are gauge invariant, which is crucial from a theoretical point of view. With this result, we derive a formula for the scalar contribution to temperature multipole coefficients. This formula is useful to test the theory, and it could indicate the physical consequence of the “delta matter” introduced in this theory. This work has been incorporated as a part of the Ph.D. thesis [48], where more details can be found.

The CMB provides cosmological constraints crucial for testing a model. Many cosmological parameters can be obtained directly from the CMB power spectrum, such as $h^2{\Omega }_b, h^2{\Omega }_c, 100\theta ,\tau , A_s$ and $n_s$ [1], while others can be derived from constraining CMB observation with SNe-Ia or BAOs. By studying the CMB anisotropies, we can address two aspects: the compatibility between the CMB power spectrum and DG fluctuations and the compatibility between CMB and SNe-Ia in the DG theory. In [49], we already fitted Planck satellite’s data with a DG model using Markov Chain Monte Carlo analysis. We also studied the compatibility between SNe Ia and CMB observation in this framework. We obtained the scalar CMB TT power spectrum and the fitted parameters needed to explain both SNe-Ia data and CMB measurements. The results are in reasonable agreement with both observations considering the analytical approximation. We also discussed whether the Hubble constant and the accelerating Universe are in concordance with the observational evidence in the DG context. With this in mind, the aim of this work is to present the full theoretical scheme for the scalar perturbation theory.

The paper is organized as follows: In Section 2, we introduce the definition of DG and its equations of motion. We then review some implications of the first law of thermodynamics, which will allow us to interpret the physical quantities of DG. Before finishing this section, we state the ansatz that the moment of equality between matter and radiation was equal in both DG and GR, and we discuss its implications. In Section 3, we study the gauge transformation for small perturbations of both geometrical and matter fields. We choose a gauge and present the gauge-invariant equations of motion for small perturbations. In Section 4, we study the evolution of cosmological perturbations, solving the equations partially when the Universe is dominated by radiation and when it is dominated by matter. In Section 5, we derive the formula for temperature fluctuation. Here, we find that this fluctuation can be expressed in three independent and gauge-invariant terms. In Section 6, we obtain a formula for temperature multipole coefficients for scalar modes. We also present preliminary numerical results for the power spectrum of the CMB. Finally, we provide conclusions and remarks.

For notation, we will use the Riemann tensor:

$$\begin{array}{c}\hfill R_{\beta \mu \nu }^{\alpha }={\partial }_{\mu }{\Gamma }_{\nu \beta }^{\alpha }-{\partial }_{\nu }{\Gamma }_{\mu \beta }^{\alpha }+{\Gamma }_{\mu \gamma }^{\alpha }{\Gamma }_{\nu \beta }^{\gamma }-{\Gamma }_{\nu \gamma }^{\alpha }{\Gamma }_{\mu \beta }^{\gamma },\end{array}$$

(1)

where the Ricci Tensor is given by $R_{\mu \nu }=R_{\mu \alpha \nu }^{\alpha }$, the Ricci scalar $R=g^{\mu \nu }{R}_{\mu \nu }$ and:

$$\begin{array}{c}\hfill {\Gamma }_{\mu \nu }^{\alpha }=\frac{1}{2}{g}^{\alpha \beta }({\partial }_{\nu }{g}_{\beta \mu }+{\partial }_{\mu }{g}_{\nu \beta }-{\partial }_{\beta }{g}_{\mu \nu })\end{array}$$

(2)

is the usual Christoffel symbol. Finally, the covariant derivative is given by:

$$\begin{array}{c}\hfill D_{\nu }{A}_{\mu }\equiv A_{\mu ;\nu }=A_{\mu ,\nu }-{\Gamma }_{\mu \nu }^{\alpha }{A}_{\alpha }.\end{array}$$

(3)

So, it is defined with the usual metric $g_{\mu \nu }$.

2. Definition of Delta Gravity

In this section, we will present the action as well all the symmetries of the model and derive the equations of motion. These approaches are based on the application of a variation called $\tilde{\delta }$, and it has the usual properties of a variation such as:

$$\begin{array}{ccc}\hfill \tilde{\delta }\left(AB\right)& =& \left(\tilde{\delta }A\right)B+A\left(\tilde{\delta }B\right),\hfill \\ \hfill \tilde{\delta }\delta A& =& \delta \tilde{\delta }A,\hfill \\ \hfill \tilde{\delta }\left({\Phi }_{\mu }\right)& =& {(\tilde{\delta }\Phi )}_{\mu },\hfill \end{array}$$

(4)

where $\delta$ is another variation. The main point of this variation is that when it is applied on a field (function, tensor, etc), it produces new elements that we define as $\tilde{\delta }$ fields, which we treat as an entirely new independent object from the original, $\tilde{\Phi }=\tilde{\delta }(\Phi )$. We use the convention that a tilde tensor is equal to the $\tilde{\delta }$ transformation of the original tensor when all its indexes are covariant. (For more detail about $\tilde{\delta }$, please see Appendix A.1.)

Now, we will present the $\tilde{\delta }$ prescription for a general action. The extension of the new symmetry is given by:

$$S_0=\int d^{n}x{\mathcal{L}}_0(\varphi ,{\partial }_i\varphi )\to S=\int d^{n}x({\mathcal{L}}_0(\varphi ,{\partial }_i\varphi )+\tilde{\delta }{\mathcal{L}}_0(\varphi ,{\partial }_i\varphi )),$$

(5)

where $S_0$ is the original action and S is the extended action in Delta Gauge Theories. When we apply this formalism to the Einstein–Hilbert action of GR, we obtain [35]

$$\begin{array}{c}\hfill S=\int d^{4}x\sqrt{-g}\left(\frac{R}{2\varkappa }+L_M-\frac{1}{2\varkappa }\left(G^{\alpha \beta }-\varkappa T^{\alpha \beta }\right){\tilde{g}}_{\alpha \beta }+{\tilde{L}}_M\right),\end{array}$$

(6)

where $\varkappa =\frac{8\pi G}{c^4}$ (hereafter, we set $c=1$), ${\tilde{g}}_{\mu \nu }=\tilde{\delta }{g}_{\mu \nu }$, $L_M$ is the matter Lagrangian and:

$$\begin{array}{c}\hfill T^{\mu \nu }=\frac{2}{\sqrt{-g}}\frac{\delta }{\delta g_{\mu \nu }}\left[\sqrt{-g}{L}_M\right],\end{array}$$

(7)

$$\begin{array}{c}\hfill {\tilde{L}}_M={\tilde{\varphi }}_I\frac{\delta L_M}{\delta {\varphi }_I}+\left({\partial }_{\mu }{\tilde{\varphi }}_I\right)\frac{\delta L_M}{\delta \left({\partial }_{\mu }{\varphi }_I\right)},\end{array}$$

(8)

where $\tilde{\varphi }=\tilde{\delta }\varphi$ are the $\tilde{\delta }$ matter fields or “delta matter” fields. The equations of motion are given by the variation of $g_{\mu \nu }$ and ${\tilde{g}}_{\mu \nu }$. It is easy to see that we obtain the usual Einstein’s equations varying the action (6) with respect to ${\tilde{g}}_{\mu \nu }$. On the other hand, variations with respect to $g_{\mu \nu }$ give the equations for ${\tilde{g}}_{\mu \nu }$:

$$\begin{array}{ccc}\hfill F^{\left(\mu \nu \right)\left(\alpha \beta \right)\rho \lambda }{D}_{\rho }{D}_{\lambda }{\tilde{g}}_{\alpha \beta }& +& \frac{1}{2}{R}^{\alpha \beta }{\tilde{g}}_{\alpha \beta }{g}^{\mu \nu }+\frac{1}{2}R{\tilde{g}}^{\mu \nu }-R^{\mu \alpha }{\tilde{g}}_{\alpha }^{\nu }-R^{\nu \alpha }{\tilde{g}}_{\alpha }^{\mu }+\frac{1}{2}{\tilde{g}}_{\alpha }^{\alpha }{G}^{\mu \nu }\hfill \\ & =& \frac{\varkappa }{\sqrt{-g}}\frac{\delta }{\delta g_{\mu \nu }}\left[\sqrt{-g}\left(T^{\alpha \beta }{\tilde{g}}_{\alpha \beta }+2{\tilde{L}}_M\right)\right],\hfill \end{array}$$

(9)

with:

$$\begin{array}{ccc}\hfill F^{\left(\mu \nu \right)\left(\alpha \beta \right)\rho \lambda }& =& P^{\left(\right(\rho \mu \left)\right(\alpha \beta \left)\right)}{g}^{\nu \lambda }+P^{\left(\right(\rho \nu \left)\right(\alpha \beta \left)\right)}{g}^{\mu \lambda }-P^{\left(\right(\mu \nu \left)\right(\alpha \beta \left)\right)}{g}^{\rho \lambda }-P^{\left(\right(\rho \lambda \left)\right(\alpha \beta \left)\right)}{g}^{\mu \nu },\hfill \\ \hfill P^{\left(\right(\alpha \beta \left)\right(\mu \nu \left)\right)}& =& \frac{1}{4}\left(g^{\alpha \mu }{g}^{\beta \nu }+g^{\alpha \nu }{g}^{\beta \mu }-g^{\alpha \beta }{g}^{\mu \nu }\right),\hfill \end{array}$$

(10)

where $\left(\mu \nu \right)$ denotes the totally symmetric combination of $\mu$ and $\nu$. It is possible to simplify (9) (see [35]) to obtain the following system of equations:

$$\begin{array}{ccc}\hfill G^{\mu \nu }& =& \varkappa T^{\mu \nu },\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill F^{\left(\mu \nu \right)\left(\alpha \beta \right)\rho \lambda }{D}_{\rho }{D}_{\lambda }{\tilde{g}}_{\alpha \beta }+\frac{1}{2}{g}^{\mu \nu }{R}^{\alpha \beta }{\tilde{g}}_{\alpha \beta }-\frac{1}{2}{\tilde{g}}^{\mu \nu }R& =& \varkappa {\tilde{T}}^{\mu \nu },\hfill \end{array}$$

(12)

where ${\tilde{T}}^{\mu \nu }=\tilde{\delta }{T}^{\mu \nu }$. The energy momentum conservation now is given by

$$\begin{array}{ccc}{D}_{\nu }{T}^{\mu \nu }\hfill & =& 0,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill D_{\nu }{\tilde{T}}^{\mu \nu }& =& \frac{1}{2}{T}^{\alpha \beta }{D}^{\mu }{\tilde{g}}_{\alpha \beta }-\frac{1}{2}{T}^{\mu \beta }{D}_{\beta }{\tilde{g}}_{\alpha }^{\alpha }+D_{\beta }\left({\tilde{g}}_{\alpha }^{\beta }{T}^{\alpha \mu }\right).\hfill \end{array}$$

(14)

Then, we are going to work with Equations (11)–(14). However, as the perturbation theory in the standard sector is well known (see [47]), we will focus on the DG sector.

One important result of DG is that photons follow geodesic trajectories given by the effective metric ${\mathbf{g}}_{\mu \nu }=g_{\mu \nu }+{\tilde{g}}_{\mu \nu }$ [35], and for an FRLW Universe, these metrics take the form (with constant curvature parameter $k=0$)

$${\overline{g}}_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-d{t}^2+a^2\left(t\right)(d{x}^2+d{y}^2+d{z}^2),$$

(15)

and

$${\tilde{\overline{g}}}_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-3F\left(t\right)d{t}^2+F\left(t\right)a^2\left(t\right)(d{x}^2+d{y}^2+d{z}^2),$$

(16)

where $F\left(t\right)$ is a time-dependent function which is determined by the solution of the unperturbed equations system, and $a\left(t\right)$ is the standard scale factor, which in Section 4 we will show is no longer the physical scale factor of the Universe. To obtain the form of ${\overline{g}}_{\mu \nu }$ and ${\tilde{\overline{g}}}_{\mu \nu }$, first we impose isotropy and homogeneity, and then, we apply the harmonic gauge $g^{\mu \nu }{\Gamma }_{\mu \nu }^{\alpha }=0$ and its tilde version (for details, see [37]). One of the implications of this effective metric is that geometry is now described by a new tridimensional metric given by [35]2 (latin indexes run from 1 to 3)

$$\begin{array}{ccc}\hfill d{l}^2& =& {\gamma }_{ij}d{x}^{i}d{x}^j,\hfill \\ \hfill {\gamma }_{ij}& =& \frac{g_{00}}{{\mathbf{g}}_{00}}\left({\mathbf{g}}_{ij}-\frac{{\mathbf{g}}_{i0}{\mathbf{g}}_{j0}}{{\mathbf{g}}_{00}}\right),\hfill \end{array}$$

(17)

while the proper time is defined by $g_{\mu \nu }$. In this case, t is the cosmic time.

2.1. DG and Thermodynamics

Now, we will study some implications of thermodynamics in cosmology for DG. Equation (17) defines the modified scale factor of this theory:

$$a_{DG}\left(t\right)=a\left(t\right)\sqrt{\frac{1+F\left(t\right)}{1+3F\left(t\right)}}.$$

(18)

Then, the volume of a cosmological sphere is now

$$V=\frac{4}{3}\pi r^3{a}_{DG}^3.$$

Any physical fluid has a density given by

$${\rho }_{DG}=\frac{U}{V},$$

(19)

where U is the internal energy and V is the volume. From the first law of thermodynamics, we have

$$\frac{dU}{dt}=T\frac{d\mathcal{S}}{dt}-P_{DG}\frac{dV}{dt}.$$

(20)

We will assume that the Universe evolved adiabatically; this means $\dot{\mathcal{S}}=0$. Then, we obtain the well-known relation for the energy conservation

$${\dot{\rho }}_{DG}=-3H_{DG}\left({\rho }_{DG}+P_{DG}\right),$$

(21)

with $H_{DG}={\dot{a}}_{DG}/a_{DG}$. In order to know the evolution of $\rho$, we need an equation of state $P\left(\rho \right)$. In [38], they showed that $H_{DG}\left(t\right)$ replaces the first Friedmann equation. Now, we know that the second Friedmann equation is the thermodynamics statement that the Universe evolves adiabatically, so the physical densities must satisfy Equation (21). If we assume $P=\omega \rho$, we found

$${\rho }_{DG}{a}_{DG}^{3(1+\omega )}={\rho }_{DG0}{a}_{DG0}^{3(1+\omega )},$$

(22)

where ${\rho }_0$ is the density at present. A crucial point in this theory is that the GR field Equations (11) and (13) are valid; then, we also have a similar relation for the densities of GR but with the standard scale factor $a\left(t\right)$, explicitly

$${\rho }_{GR}{a}^{3(1+\omega )}={\rho }_{GR0}{a}_0^{3(1+\omega )}.$$

(23)

Then, we can relate both densities by the ratio between them

$$\frac{{\rho }_{DG}}{{\rho }_{GR}}{\left(\sqrt{\frac{1+F\left(t\right)}{1+3F\left(t\right)}}\right)}^{3(1+\omega )}=constant\left(\omega \right).$$

(24)

This ratio will be vitally important when we study the perturbations of the system. Because we will study the evolution of fractional perturbations at the last-scattering time defined as

$${\delta }_{GR\alpha }=\frac{\delta {\rho }_{GR\alpha }}{{\overline{\rho }}_{GR\alpha }+{\overline{p}}_{GR\alpha }},$$

(25)

where $\alpha$ runs between $\gamma$, $\nu$, B and D (photons, neutrinos, baryons and dark matter, respectively). If we consider the results from [38], at the moment of last-scattering ($T\sim 3000K$), we obtain

$$\sqrt{\frac{1+F\left(t_{ls}\right)}{1+3F\left(t_{ls}\right)}}\sim 1.$$

(26)

This mean that at that moment, the physical density was proportional to the densities of GR, and without a loss of generality, we can take

$${\delta }_{DG\alpha }\left(t_{ls}\right)={\delta }_{GR\alpha }\left(t_{ls}\right)\equiv {\delta }_{\alpha }\left(t_{ls}\right),$$

(27)

as it will be introduced in Section 4. In fact, Equation (26) is valid for a wide range of times, from the beginning of the Universe ($z\to \infty$) until $z\sim 10$, so this approximation is valid in the study of primordial perturbations in DG when using the equations of GR. On the other hand, the number density (number of photons over the volume) at equilibrium with matter at temperature T is

$$n_T\left(\nu \right)d\nu =\frac{8\pi {\nu }^{2}d\nu }{e^{\frac{h\nu }{k_{B}T}}-1};$$

(28)

After decoupling, photons travel freely from the surface of last scattering to us. So, the number of photons is conserved

$$dN=n_{T_{ls}}\left({\nu }_{ls}\right)d{\nu }_{ls}d{V}_{ls}=n_T\left(\nu \right)d\nu dV,$$

(29)

as frequencies are redshifted by $\nu ={\nu }_{ls}{a}_{DG}\left(t_{ls}\right)/a_{DG}$, and the volume $V=V_{ls}{a}_{DG}^3/a_{DG}^3\left(t_{ls}\right)$. We find that in order to keep the form of a black body distribution, temperature in the number density should evolve as $T=T_{ls}{a}_{DG}\left(t_{ls}\right)/a_{DG}$.

2.2. Equality Time $t_{EQ}$

After concluding this section, there is an ansatz that we need to propose in order to be completely consistent when solving the cosmological perturbation theory in the next section. This is about when the radiation was equal to the non-relativistic matter. We state that the moment when radiation and matter were equal at some $t_{EQ}$ is the same in GR as in DG. The implication of this statement is the following: let us consider the ratio of the matter and radiation densities of GR (23)

$$\begin{array}{c}\hfill \frac{{\rho }_{GRM}}{{\rho }_{GRR}}=\frac{Y}{C},\end{array}$$

(30)

We remind that $C={\Omega }_R/{\Omega }_M$. Then, the moment of equality in GR corresponds to $Y_{EQ}=C$. On the other hand, if we consider the same ratio but now between the physical densities using (22), we obtain

$$\begin{array}{c}\hfill \frac{{\rho }_{DG;M}}{{\rho }_{DGR}}=\frac{Y_{DG}}{C_{DG}},\end{array}$$

(31)

where $C_{DG}={\Omega }_{DGR}/{\Omega }_{DGM}$. Then, in the equality, we need to impose $Y_{DG}\left(Y_{EQ}\right)=C_{DG}$, explicitly

$$C_{DG}=C\frac{\sqrt{\frac{1+F\left(C\right)}{1+3F\left(C\right)}}}{\sqrt{\frac{1+F\left(1\right)}{1+3F\left(1\right)}}},$$

(32)

if we take the value from [42] (they used $L_2$ instead of L, but these are the same quantity also), $C\sim {10}^{-4}$ and $L\sim 0.45$ implies $F\left(C\right)\sim {10}^{-3}<<1$ and $F\left(1\right)\sim -L/3$, then

$$C_{DG}=C\sqrt{\frac{1-L}{1-L/3}}.$$

(33)

This means that the total density of matter and radiation today depend explicitly on the geometry measured with L [42].

3. Perturbation Theory

Now, we perturbed the metric as the following

$$\begin{array}{c}\hfill g_{\mu \nu }={\overline{g}}_{\mu \nu }+h_{\mu \nu },\end{array}$$

(34)

$$\begin{array}{c}\hfill {\tilde{g}}_{\mu \nu }={\tilde{\overline{g}}}_{\mu \nu }+{\tilde{h}}_{\mu \nu }.\end{array}$$

(35)

Then, we follow the standard method, known as Scalar–Vector–Tensor decomposition [50]. This decomposition depends on four scalar functions ($A,B,E$ and H), two vector functions ($C_i$ and $G_i$), and one tensor function ($D_{ij}$) (with their respective delta part). This process allows us to study those sectors independently. Therefore, the perturbations are

$$h_{00}=-E{h}_{i0}=a\left[\frac{\partial H}{\partial x^i}+G_i\right]h_{ij}=a^2\left[A{\delta }_{ij}+\frac{{\partial }^{2}B}{\partial x^i\partial x^j}+\frac{\partial C_i}{\partial x^j}+\frac{\partial C_j}{\partial x^i}+D_{ij}\right],$$

(36)

where

$$\frac{\partial C_i}{\partial x^i}=\frac{\partial G_i}{\partial x^i}=0\frac{\partial D_{ij}}{\partial x^j}=0D_{ii}=0.$$

(37)

This decomposition must be equivalent for ${\tilde{h}}_{\mu \nu }$ (by group theory):

$${\tilde{h}}_{00}=-\tilde{E}{\tilde{h}}_{i0}=a\left[\frac{\partial \tilde{H}}{\partial x^i}+{\tilde{G}}_i\right]{\tilde{h}}_{ij}=a^2\left[\tilde{A}{\delta }_{ij}+\frac{{\partial }^2\tilde{B}}{\partial x^i\partial x^j}+\frac{\partial {\tilde{C}}_i}{\partial x^j}+\frac{\partial {\tilde{C}}_j}{\partial x^i}+{\tilde{D}}_{ij}\right],$$

(38)

with

$$\frac{\partial {\tilde{C}}_i}{\partial x^i}=\frac{\partial {\tilde{G}}_i}{\partial x^i}=0\frac{\partial {\tilde{D}}_{ij}}{\partial x^j}=0{\tilde{D}}_{ii}=0.$$

(39)

If we replace perturbations in (11)–(14), we obtain the equations for the perturbations. However, there are degrees of freedom that we have to take into account to have physical solutions. In the next subsection, we show how to choose a gauge to delete the nonphysical solutions.

3.1. Choosing a Gauge

Under a space–time coordinate transformation, the metric perturbations transform as 3

$$\Delta h_{\mu \nu }\left(x\right)=-{\overline{g}}_{\lambda \nu }\left(x\right)\frac{\partial {ϵ}^{\lambda }}{\partial x^{\mu }}-{\overline{g}}_{\mu \lambda }\left(x\right)\frac{\partial {ϵ}^{\lambda }}{\partial x^{\nu }}-\frac{\partial {\overline{g}}_{\mu \nu }}{\partial x^{\lambda }}{ϵ}^{\lambda }.$$

(40)

In more detail,

$$\begin{array}{ccc}\hfill \Delta h_{ij}& =& -\frac{\partial {ϵ}_i}{\partial x^j}-\frac{\partial {ϵ}_j}{\partial x^i}+2a\dot{a}{\delta }_{ij}{ϵ}_0,\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill \Delta h_{i0}& =& -\frac{\partial {ϵ}_i}{\partial t}-\frac{\partial {ϵ}_0}{\partial x^i}+2\frac{\dot{a}}{a}{ϵ}_i,\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill \Delta h_{00}& =& -2\frac{\partial {ϵ}_0}{\partial t}.\hfill \end{array}$$

(43)

For delta perturbations, we obtain

$$\Delta {\tilde{h}}_{\mu \nu }=-{\tilde{\overline{g}}}_{\mu \lambda }\frac{\partial {ϵ}^{\lambda }}{\partial x^{\nu }}-{\tilde{\overline{g}}}_{\lambda \nu }\frac{\partial {ϵ}^{\lambda }}{\partial x^{\mu }}-\frac{\partial {\tilde{\overline{g}}}_{\mu \nu }}{\partial x^{\lambda }}{ϵ}^{\lambda }-{\overline{g}}_{\mu \lambda }\frac{\partial {\widetilde{ϵ}}^{\lambda }}{\partial x^{\nu }}-{\overline{g}}_{\lambda \nu }\frac{\partial {\widetilde{ϵ}}^{\lambda }}{\partial x^{\mu }}-\frac{\partial {\overline{g}}_{\mu \nu }}{\partial x^{\lambda }}{\widetilde{ϵ}}^{\lambda }.$$

(44)

In more detail,

$$\begin{array}{ccc}\hfill \Delta {\tilde{h}}_{ij}& =& -F\frac{\partial {ϵ}_i}{\partial x^j}-F\frac{\partial {ϵ}_j}{\partial x^i}-\frac{\partial {\widetilde{ϵ}}_j}{\partial x^i}-\frac{\partial {\widetilde{ϵ}}_i}{\partial x^j}+\left[{ϵ}_0\left(2Fa\dot{a}+\dot{F}{a}^2\right)+2{\widetilde{ϵ}}_{0}a\dot{a}\right]{\delta }_{ij},\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill \Delta {\tilde{h}}_{i0}& =& -F\frac{\partial {ϵ}_i}{\partial t}-3F\frac{\partial {ϵ}_0}{\partial x^i}-\frac{\partial {\widetilde{ϵ}}_i}{\partial t}-\frac{\partial {\widetilde{ϵ}}_0}{\partial x^i}+2F\frac{\dot{a}}{a}{ϵ}_i+2\frac{\dot{a}}{a}{\widetilde{ϵ}}_i,\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill \Delta {\tilde{h}}_{00}& =& -3{ϵ}_0\dot{F}-6F\frac{\partial {ϵ}_0}{\partial t}-2\frac{\partial {\widetilde{ϵ}}_0}{dt},\hfill \end{array}$$

(47)

where $ϵ$ and $\widetilde{ϵ}=\tilde{\delta }ϵ$ define the coordinates transformation. In addition, we raised and lowered the index using ${\overline{g}}_{\mu \nu }$, so ${ϵ}^0=-{ϵ}_0$, ${\widetilde{ϵ}}^0=-{\widetilde{ϵ}}_0$, ${ϵ}^i=a^{-2}{ϵ}_i$ and ${\widetilde{ϵ}}^j=a^{-2}{\widetilde{ϵ}}_j$. Following the standard procedure, we decompose the spatial part of ${ϵ}^{\mu }$ and ${\widetilde{ϵ}}^{\mu }$ into the gradient of a spatial scalar plus a divergenceless vector:

$$\begin{array}{ccc}\hfill {ϵ}_i& =& {\partial }_i{ϵ}^S+{ϵ}_i^V,{\partial }_i{ϵ}^V=0,\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill {\widetilde{ϵ}}_i& =& {\partial }_i{\widetilde{ϵ}}^S+{\widetilde{ϵ}}_i^V,{\partial }_i{\widetilde{ϵ}}^V=0.\hfill \end{array}$$

(49)

Thus, we can compare Equations (36) and (38) with (41)–(43) and (45)–(47) to obtain the gauge transformations of the metric components:

$$\begin{array}{ccc}\hfill \Delta A& =& \frac{2\dot{a}}{a}{ϵ}_0,\Delta B=-\frac{2}{a^2}{ϵ}^S,\hfill \\ \hfill \Delta C_i& =& -\frac{1}{a^2}{ϵ}_i^V,\Delta D_{ij}=0,\Delta E=2{\dot{ϵ}}_0,\hfill \\ \hfill \Delta H& =& \frac{1}{a}\left(-{ϵ}_0-{\dot{ϵ}}^S+\frac{2\dot{a}}{a}{ϵ}^S\right),\Delta G_i=\frac{1}{a}\left(-{\dot{ϵ}}_i^V+\frac{2\dot{a}}{a}{ϵ}_i^V\right),\hfill \end{array}$$

(50)

and

$$\begin{array}{ccc}\hfill \Delta \tilde{A}& =& \left(\frac{2\dot{a}F}{a}+\dot{F}\right){ϵ}_0+2\frac{\dot{a}}{a}{\widetilde{ϵ}}_0,\Delta \tilde{B}=-\frac{2}{a^2}\left(F{ϵ}^S+{\widetilde{ϵ}}^S\right),\hfill \\ \hfill \Delta {\tilde{C}}_i& =& -\frac{1}{a^2}\left(F{ϵ}_i^V+{\widetilde{ϵ}}_i^V\right),\Delta {\tilde{D}}_{ij}=0,\Delta \tilde{E}=6F{\dot{ϵ}}_0+3\dot{F}{ϵ}_0+2{\dot{\widetilde{ϵ}}}_0,\hfill \\ \hfill \Delta \tilde{H}& =& \frac{1}{a}\left(-3F{ϵ}_0-{\widetilde{ϵ}}_0-F{\dot{ϵ}}^S-{\dot{\widetilde{ϵ}}}^S+\frac{2F\dot{a}}{a}{ϵ}^S+\frac{2\dot{a}}{a}{\widetilde{ϵ}}^S\right),\hfill \\ \hfill \Delta {\tilde{G}}_i& =& \frac{1}{a}\left(-F{\dot{ϵ}}_i^V-{\dot{\widetilde{ϵ}}}_i^V+\frac{2F\dot{a}}{a}{ϵ}_i^V+\frac{2\dot{a}}{a}{\widetilde{ϵ}}_i^V\right).\hfill \end{array}$$

(51)

There are different scenarios in which we can continue with the calculations when we impose conditions on the parameters ${ϵ}_{\mu }$ and ${\widetilde{ϵ}}_{\mu }$. However, before discussing this, we will study the gauge transformation of energy-momentum tensors $T_{\mu \nu }$ and ${\tilde{T}}_{\mu \nu }$.

3.2. $T_{\mu \nu }$ and ${\tilde{T}}_{\mu \nu }$

Now, we will decompose the energy-momentum tensors $T_{\mu \nu }$ and ${\tilde{T}}_{\mu \nu }$ in the same way. For a perfect fluid, we would have (for more details, see [37])

$$T_{\mu \nu }=p{g}_{\mu \nu }+(\rho +p)u_{\mu }{u}_{\nu },$$

(52)

while for ${\tilde{T}}_{\mu \nu }$ (here $u_{\mu }^T=e_{\mu \alpha }{\tilde{u}}^{\alpha }$, where $e_{\mu \alpha }$ is the Vierbein $g_{\mu \nu }=e_{\mu \alpha }{e}_{\nu \beta }{\eta }^{\alpha \beta }$, with $\eta$ being the Minkowski metric); for a full derivation, please see [35,37]):

$${\tilde{T}}_{\mu \nu }=\tilde{p}{g}_{\mu \nu }+p{\tilde{g}}_{\mu \nu }+(\tilde{\rho }+\tilde{p})u_{\mu }{u}_{\nu }+(\rho +p)\left(\frac{1}{2}({\tilde{g}}_{\mu \alpha }{u}_{\nu }{u}^{\alpha }+{\tilde{g}}_{\nu \alpha }{u}_{\mu }{u}^{\alpha })+u_{\mu }^T{u}_{\nu }+u_{\mu }{u}_{\nu }^T\right),$$

(53)

where

$$g^{\mu \nu }{u}_{\mu }{u}_{\nu }=-1,$$

(54)

$$g^{\mu \nu }{u}_{\mu }{u}_{\nu }^T=0.$$

(55)

The tensors $g_{\mu \nu }$ and ${\tilde{g}}_{\mu \nu }$ are defined in (15) and (16), respectively; besides, we consider

$$\begin{array}{ccc}\hfill p& =& \overline{p}+\delta p,\hfill \\ \hfill \rho & =& \overline{\rho }+\delta \rho ,\hfill \\ \hfill u_{\mu }& =& {\overline{u}}_{\mu }+\delta u_{\mu },\hfill \\ \hfill \tilde{p}& =& \tilde{\overline{p}}+\delta \tilde{p},\hfill \\ \hfill \tilde{\rho }& =& \tilde{\overline{\rho }}+\delta \tilde{\rho },\hfill \\ \hfill u_{\mu }^T& =& {\overline{u}}_{\mu }^T+\delta u_{\mu }^T.\hfill \end{array}$$

(56)

Usually, the equation of state is given by $p\left(\rho \right)$, so we could reduce this system. For now, we will work in the generic case. When we work in the frame ${\overline{u}}_{\mu }=(-1,0,0,0)$, we have ${\overline{u}}_{\mu }^T=0$, and the normalization conditions (54) and (55) give

$$\begin{array}{ccc}\hfill \delta u^0& =& \delta u_0=\frac{h_{00}}{2}\hfill \\ \hfill \delta u_0^T& =& \delta u_T^0=0\hfill \end{array}$$

(57)

while $\delta u_i$ and $\delta u_i^T$ are independent dynamical variables (note that $\delta u^{\mu }\equiv \delta \left(g^{\mu \nu }{u}_{\nu }\right)$ is not given by ${\overline{g}}^{\mu \nu }\delta u_{\nu }$. The same is true for $\delta u_T^{\mu }$). Then, the first-order perturbations for both energy-momentum tensors (a perfect fluid) are

$$\delta T_{\mu \nu }=\overline{p}{h}_{\mu \nu }+\delta p{\overline{g}}_{\mu \nu }+(\overline{p}+\overline{\rho })({\overline{u}}_{\mu }\delta u_{\nu }+\delta u_{\mu }{u}_{\nu })+(\delta p+\delta \rho ){\overline{u}}_{\mu }{\overline{u}}_{\nu },$$

(58)

Therefore,

$$\delta T_{ij}=\overline{p}{h}_{ij}+a^2{\delta }_{ij}\delta p,,\delta T_{i0}=\overline{p}{h}_{i0}-(\overline{p}+\overline{\rho })\delta u_i,,\delta T_{00}=-\overline{\rho }{h}_{00}+\delta \rho .$$

(59)

While

$$\begin{array}{ccc}\hfill \delta {\tilde{T}}_{\mu \nu }& =& \tilde{\overline{p}}{h}_{\mu \nu }+\delta \tilde{p}{\overline{g}}_{\mu \nu }+\overline{p}{\tilde{h}}_{\mu \nu }+\delta p{\tilde{\overline{g}}}_{\mu \nu }+(\tilde{\overline{\rho }}+\tilde{\overline{p}})({\overline{u}}_{\mu }\delta u_{\nu }+\delta u_{\mu }{\overline{u}}_{\nu })\hfill \\ & +& (\delta \tilde{\rho }+\delta \tilde{p}){\overline{u}}_{\mu }{\overline{u}}_{\nu }+(\overline{\rho }+\overline{p})\left\{\frac{1}{2}\left[{\tilde{\overline{g}}}_{\mu \alpha }({\overline{u}}_{\nu }\delta u^{\alpha }+\delta u_{\nu }{\overline{u}}^{\alpha })+{\tilde{h}}_{\mu \alpha }{\overline{u}}_{\nu }{\overline{u}}^{\alpha }\right.\right.\hfill \\ & +& \left.\left.{\tilde{\overline{g}}}_{\nu \alpha }({\overline{u}}_{\mu }\delta u^{\alpha }+\delta u_{\mu }{\overline{u}}^{\alpha })+{\tilde{h}}_{\nu \alpha }{\overline{u}}_{\mu }{\overline{u}}^{\alpha }\right]+{\overline{u}}_{\mu }^T\delta u_{\nu }+\delta u_{\mu }^T{\overline{u}}_{\nu }+{\overline{u}}_{\mu }\delta u_{\nu }^T+\delta u_{\mu }{\overline{u}}_{\nu }^T\right\}\hfill \\ & +& (\delta \rho +\delta p)\left\{\frac{1}{2}\left[{\tilde{\overline{g}}}_{\mu \alpha }{\overline{u}}_{\nu }{\overline{u}}^{\alpha }+{\tilde{\overline{g}}}_{\nu \alpha }{\overline{u}}_{\mu }{\overline{u}}^{\alpha }\right]+{\overline{u}}_{\mu }^T{\overline{u}}_{\nu }+{\overline{u}}_{\mu }{\overline{u}}_{\nu }^T\right\},\hfill \end{array}$$

(60)

and

$$\begin{array}{ccc}\hfill \delta {\tilde{T}}_{00}& =& -\tilde{\overline{\rho }}{h}_{00}-\overline{\rho }{\tilde{h}}_{00}+3F\delta \rho +\delta \tilde{\rho },\hfill \\ \hfill \delta {\tilde{T}}_{i0}& =& \tilde{\overline{p}}{h}_{i0}+\overline{p}{\tilde{h}}_{i0}-(\tilde{\overline{\rho }}+\tilde{\overline{p}})\delta u_i+(\overline{\rho }+\overline{p})\left\{\frac{1}{2}[F{h}_{i0}-{\tilde{h}}_{i0}-4F\delta u_i]-\delta u_i^T\right\},\hfill \\ \hfill \delta {\tilde{T}}_{ij}& =& \tilde{\overline{p}}{h}_{ij}+\delta \tilde{p}{a}^2{\delta }_{ij}+\overline{p}{\tilde{h}}_{ij}+\delta pF{a}^2{\delta }_{ij},\hfill \end{array}$$

(61)

where we used $\delta u^{\alpha }\equiv \delta \left(g^{\alpha \beta }{u}_{\beta }\right)={\overline{g}}^{\alpha \beta }\delta u_{\beta }-h^{\alpha \beta }{\overline{u}}_{\beta }$.

Generally, we decompose $\delta u_i\left(\delta u_i^T\right)$ into the gradient of a scalar velocity potential $\delta u\left(\delta \tilde{u}\right)$ and a divergenceless vector $\delta u_i^V\left(\delta {\tilde{u}}_i^V\right)$, and the dissipative corrections to the inertia tensor are added as follows

$$\begin{array}{ccc}\hfill \delta T_{ij}& =& \overline{p}{h}_{ij}+a^2\left[{\delta }_{ij}\delta p+{\partial }_i{\partial }_j{\pi }^S+{\partial }_i{\pi }_j^V+{\partial }_j{\pi }_i^V+{\pi }_{ij}^T\right],\hfill \end{array}$$

(62)

$$\begin{array}{ccc}\hfill \delta T_{i0}& =& \overline{p}{h}_{i0}-(\overline{p}+\overline{\rho })\left({\partial }_i\delta u+\delta u_i^V\right),\hfill \end{array}$$

(63)

$$\begin{array}{ccc}\hfill \delta T_{00}& =& -\overline{\rho }{h}_{00}+\delta \rho ,\hfill \end{array}$$

(64)

and

$$\begin{array}{ccc}\hfill \delta {\tilde{T}}_{ij}& =& \tilde{\overline{p}}{h}_{ij}+a^2\left[{\delta }_{ij}\delta \tilde{p}+{\partial }_i{\partial }_j{\tilde{\pi }}^S+{\partial }_i{\tilde{\pi }}_j^V+{\partial }_j{\tilde{\pi }}_i^V+{\tilde{\pi }}_{ij}^T\right]+\overline{p}{\tilde{h}}_{ij}\hfill \\ & +& F{a}^2\left[{\delta }_{ij}\delta p+{\partial }_i{\partial }_j{\pi }^S+{\partial }_i{\pi }_j^V+{\partial }_j{\pi }_i^V+{\pi }_{ij}^T\right],\hfill \\ \hfill \delta {\tilde{T}}_{i0}& =& \tilde{\overline{p}}{h}_{i0}+\overline{p}{\tilde{h}}_{i0}-(\tilde{\overline{\rho }}+\tilde{\overline{p}})({\partial }_i\partial u+\partial u_i^V)\hfill \end{array}$$

(65)

$$\begin{array}{ccc}& +& (\overline{\rho }+\overline{p})\left\{\frac{1}{2}[F{h}_{i0}-{\tilde{h}}_{i0}-4F({\partial }_i\delta u+\delta u_i^V)]-{\partial }_i\delta \tilde{u}+\delta {\tilde{u}}_i^V\right\},\hfill \end{array}$$

(66)

$$\begin{array}{ccc}\hfill \delta {\tilde{T}}_{00}& =& -\tilde{\overline{\rho }}{h}_{00}-\overline{\rho }{\tilde{h}}_{00}+3F\delta \rho +\delta \tilde{\rho },\hfill \end{array}$$

(67)

where ${\pi }_i^V\left({\tilde{\pi }}_i^V\right)$, ${\pi }_{ij}^T\left({\tilde{\pi }}_{ij}^T\right)$ and $\delta u_i^V\left(\delta {\tilde{u}}_i^V\right)$ satisfy similar conditions to (37) and (39). These conditions are (expressed before as $C_i\left({\tilde{C}}_i\right)$, $D_{ij}\left({\tilde{D}}_{ij}\right)$$G_i\left({\tilde{G}}_i\right)$):

$${\partial }_i{\pi }_i^V={\partial }_i{\tilde{\pi }}_i^V={\partial }_i\delta u_i^V={\partial }_i\delta {\tilde{u}}_i^V=0{\partial }_i{\pi }_{ij}^T={\partial }_i{\tilde{\pi }}_{ij}^T=0,{\pi }_{ii}^T={\tilde{\pi }}_{ii}^T=0.$$

(68)

3.3. Gauge Transformations for the Energy-Momentum Tensors

The gauge transformation for $T_{\mu \nu }$ is given by

$$\Delta \delta T_{\mu \nu }\left(x\right)=-{\overline{T}}_{\lambda \nu }\left(x\right)\frac{\partial {ϵ}^{\lambda }}{\partial x^{\mu }}-{\overline{T}}_{\mu \lambda }\left(x\right)\frac{\partial {ϵ}^{\lambda }}{\partial x^{\nu }}-\frac{\partial {\overline{T}}_{\mu \nu }}{\partial x^{\lambda }}{ϵ}^{\lambda },$$

(69)

where the components are

$$\begin{array}{ccc}\hfill \Delta \delta T_{ij}& =& -\overline{p}\left(\frac{\partial {ϵ}_i}{\partial x^j}+\frac{\partial {ϵ}_j}{\partial x^i}\right)+\frac{\partial }{\partial t}\left(a^2\overline{p}\right){\delta }_{ij}{ϵ}_0,\hfill \end{array}$$

(70)

$$\begin{array}{ccc}\hfill \Delta \delta T_{i0}& =& -\overline{p}\frac{\partial {ϵ}_i}{\partial t}+\overline{\rho }\frac{\partial {ϵ}_0}{\partial x^i}+2\overline{p}\frac{\dot{a}}{a}{ϵ}_i,\hfill \end{array}$$

(71)

$$\begin{array}{ccc}\hfill \Delta \delta T_{00}& =& 2\overline{\rho }\frac{\partial {ϵ}_0}{\partial t}+\dot{\overline{\rho }}{ϵ}_0.\hfill \end{array}$$

(72)

While the gauge transformation of $\delta {\tilde{T}}_{\mu \nu }$ is given by

$$\Delta \delta {\tilde{T}}_{\mu \nu }=-{\tilde{\overline{T}}}_{\mu \lambda }\frac{\partial {ϵ}^{\lambda }}{\partial x^{\nu }}-{\tilde{\overline{T}}}_{\lambda \nu }\frac{\partial {ϵ}^{\lambda }}{\partial x^{\mu }}-\frac{\partial {\tilde{\overline{T}}}_{\mu \nu }}{\partial x^{\lambda }}{ϵ}^{\lambda }-{\overline{T}}_{\mu \lambda }\frac{\partial {\widetilde{ϵ}}^{\lambda }}{\partial x^{\nu }}-{\overline{T}}_{\lambda \nu }\frac{\partial {\widetilde{ϵ}}^{\lambda }}{\partial x^{\mu }}-\frac{\partial {\overline{T}}_{\mu \nu }}{\partial x^{\lambda }}{\widetilde{ϵ}}^{\lambda },$$

(73)

where the components are

$$\begin{array}{ccc}\hfill \Delta \delta {\tilde{T}}_{ij}& =& -(\tilde{\overline{p}}+\overline{p}F)\frac{\partial {ϵ}_i}{\partial x^j}-(\tilde{\overline{p}}+\overline{p}F)\frac{\partial {ϵ}_j}{\partial x^i}-\overline{p}\frac{\partial {\widetilde{ϵ}}_j}{\partial x^i}-\overline{p}\frac{\partial {\widetilde{ϵ}}_i}{\partial x^j}+\left[{ϵ}_0\frac{\partial }{\partial t}\left[a^2\left(\tilde{\overline{p}}+\overline{p}F\right)\right]+\frac{\partial }{\partial t}\left(a^2\overline{p}\right){\widetilde{ϵ}}_0\right]{\delta }_{ij}\hfill \end{array}$$

(74)

$$\begin{array}{ccc}\hfill \Delta \delta {\tilde{T}}_{i0}& =& -(\tilde{\overline{p}}+\overline{p}F)\frac{\partial {ϵ}_i}{\partial t}+(\tilde{\overline{\rho }}+3F\overline{\rho })\frac{\partial {ϵ}_0}{\partial x^i}-\overline{p}\frac{\partial {\widetilde{ϵ}}_i}{\partial t}+\overline{\rho }\frac{\partial {\widetilde{ϵ}}_0}{\partial x^i}+2(\tilde{\overline{p}}+\overline{p}F)\frac{\dot{a}}{a}{ϵ}_i+2\overline{p}\frac{\dot{a}}{a}{\widetilde{ϵ}}_i\hfill \end{array}$$

(75)

$$\begin{array}{ccc}\hfill \Delta \delta {\tilde{T}}_{00}& =& {ϵ}_0\frac{\partial }{\partial t}\left(\tilde{\overline{\rho }}+3F\overline{\rho }\right)+2(\tilde{\overline{\rho }}+3F\overline{\rho })\frac{\partial {ϵ}_0}{\partial t}+\dot{\overline{\rho }}{\widetilde{ϵ}}_0+2\overline{\rho }\frac{\partial {\widetilde{ϵ}}_0}{dt}.\hfill \end{array}$$

(76)

${ϵ}_i$ and ${\widetilde{ϵ}}_i$ were decomposed in (48) to write these gauge transformations in terms of the scalar, vector and tensor components. The transformations (41)–(43) and (45)–(47) with (70) and (71) and (74)–(76) give the gauge transformation for the pressure, energy density and velocity potential:

$$\Delta \delta p=\dot{\overline{p}}{ϵ}_0,\Delta \delta \rho =\dot{\overline{\rho }}{ϵ}_0,\Delta \delta u=-{ϵ}_0.$$

(77)

The other ingredients of the energy-momentum tensor are gauge invariants:

$$\Delta {\pi }^S=\Delta {\pi }_i^V=\Delta {\pi }_{ij}^T=\Delta \delta u_i^V=0.$$

(78)

Nevertheless, the other transformations are

$$\begin{array}{ccc}\hfill \Delta \delta \tilde{\rho }& =& \frac{\partial }{\partial t}(\tilde{\overline{\rho }}+3F\overline{\rho }){ϵ}_0+2(\tilde{\overline{\rho }}+3F\overline{\rho })\dot{{ϵ}_0}+\dot{\overline{\rho }}{\widetilde{ϵ}}_0+2\overline{\rho }{\dot{\widetilde{ϵ}}}_0-\tilde{\overline{\rho }}\Delta E\hfill \\ & -& 3F\overline{\rho }\Delta \tilde{E}-3F\Delta \delta \rho ,\hfill \end{array}$$

(79a)

$$\begin{array}{ccc}\hfill \Delta \delta \tilde{p}& =& \frac{1}{a^2}\frac{\partial }{\partial t}\left[a^2\left(\tilde{\overline{p}}+\overline{p}F\right)\right]{ϵ}_0+\frac{1}{a^2}\frac{\partial }{\partial t}\left(a^2\overline{\rho }\right){\widetilde{ϵ}}_0-\tilde{\overline{p}}\Delta A-\overline{p}F\Delta \tilde{A}-F\Delta \delta p,\hfill \end{array}$$

(79b)

$$\begin{array}{ccc}\hfill \Delta \delta \tilde{u}& =& \frac{1}{(\overline{\rho }+\overline{p})}\left\{(\tilde{\overline{p}}+\overline{p}F){\dot{ϵ}}^S-(\tilde{\overline{\rho }}+3F\overline{\rho }){ϵ}_0+\overline{p}{\dot{\widetilde{ϵ}}}^S-\overline{\rho }{\widetilde{ϵ}}_0-2(\tilde{\overline{p}}+\overline{p}F)\frac{\dot{a}}{a}{ϵ}^S\right.\hfill \\ & -& \left.2\overline{p}\frac{\dot{a}}{a}{\widetilde{ϵ}}^S+\tilde{\overline{p}}a\Delta H+\overline{p}a\Delta \tilde{H}-(\overline{\rho }+\overline{p})\left[\frac{1}{2}(1-F)a\Delta \tilde{H}+2F\Delta \delta u\right]\right\},\hfill \end{array}$$

(79c)

$$\begin{array}{ccc}\hfill \Delta \delta {\tilde{u}}_i^V& =& \frac{1}{(\overline{\rho }+\overline{p})}\left\{(\tilde{\overline{p}}+\overline{p}F){\dot{ϵ}}_i^V+\overline{p}{\dot{\widetilde{ϵ}}}_i^V-2(\tilde{\overline{p}}+\overline{p}F)\frac{\dot{a}}{a}{ϵ}_i^V-2\overline{p}\frac{\dot{a}}{a}{\widetilde{ϵ}}_i^V+\tilde{\overline{p}}a\Delta G_i\right.\hfill \\ & +& \left.\overline{p}a\Delta {\tilde{G}}_i-\frac{1}{2}(\overline{\rho }+\overline{p})(1-F)a\Delta {\tilde{G}}_i\right\},\hfill \end{array}$$

(79d)

$$\begin{array}{ccc}\hfill \Delta {\delta \tilde{\pi }}^S& =& -\frac{2}{a^2}(\tilde{\overline{p}}+\overline{p}F){ϵ}^S-2\frac{\overline{p}}{a^2}{\widetilde{ϵ}}^S-\tilde{\overline{p}}\Delta B-\overline{p}F\Delta \tilde{B},\hfill \end{array}$$

(79e)

$$\begin{array}{ccc}\hfill \Delta {\delta \tilde{\pi }}_i^V& =& -\frac{1}{a^2}(\tilde{\overline{p}}+\overline{p}F){ϵ}_i^V-\frac{\overline{p}}{a^2}{\widetilde{ϵ}}_i^V-\tilde{\overline{p}}\Delta C_i-\overline{p}F\Delta {\tilde{C}}_i,\hfill \end{array}$$

(79f)

$$\begin{array}{ccc}\hfill \Delta {\delta \tilde{\pi }}_{ij}& =& 0.\hfill \end{array}$$

(79g)

The results given in (50), (51) and (77) are used to obtain

$$\begin{array}{ccc}\hfill \Delta \delta \tilde{\rho }& =& \dot{\tilde{\overline{\rho }}}{ϵ}_0+(\dot{\overline{\rho }}-3F\overline{\rho }){\widetilde{ϵ}}_0,\hfill \end{array}$$

(80a)

$$\begin{array}{ccc}\hfill \Delta \delta \tilde{p}& =& \dot{\tilde{\overline{p}}}{ϵ}_0+\dot{\overline{p}}{\widetilde{ϵ}}_0,\hfill \end{array}$$

(80b)

$$\begin{array}{ccc}\hfill \Delta \delta \tilde{u}& =& \left[(1-3F)\frac{F}{2}-\frac{(\tilde{\overline{p}}+\tilde{\overline{\rho }})}{(\overline{p}+\overline{\rho })}\right]{ϵ}_0-\frac{1}{2}(1+F){\widetilde{ϵ}}_0-(1-F)\frac{\dot{a}}{a}\left(F{ϵ}^S+{\widetilde{ϵ}}^S\right)\hfill \\ & +& \frac{1}{2}(1-F)\left(F{\dot{ϵ}}^S+{\dot{\widetilde{ϵ}}}^S\right),\hfill \end{array}$$

(80c)

$$\begin{array}{ccc}\hfill \Delta {\delta \tilde{u}}_i^V& =& \frac{1}{2}(1-F)[F{\dot{ϵ}}_i^V+{\dot{\widetilde{ϵ}}}_i^V-2\frac{\dot{a}}{a}F{ϵ}_i^V-2\frac{\dot{a}}{a}{\widetilde{ϵ}}_i^V],\hfill \end{array}$$

(80d)

$$\begin{array}{ccc}\hfill \Delta {\delta \tilde{\pi }}^S& =& 0,\hfill \end{array}$$

(80e)

$$\begin{array}{ccc}\hfill \Delta {\delta \tilde{\pi }}_i^V& =& 0,\hfill \end{array}$$

(80f)

$$\begin{array}{ccc}\hfill \Delta {\delta \tilde{\pi }}_{ij}& =& 0.\hfill \end{array}$$

(80g)

As we said before, there are different choices for the $ϵ$ and $\widetilde{ϵ}$ parameters to fix all the gauge freedoms. The most common and well-known gauges are the Newtonian gauge and synchronous gauge. The former fix ${ϵ}^S$ such that $B=0$, and we choose ${ϵ}_0$ such that $H=0$ (in Equation (50)). In DG, this choice is extended, imposing similar conditions in (51) for ${\widetilde{ϵ}}^S$ and ${\widetilde{ϵ}}_0$, such that $\tilde{B}=\tilde{H}=0$. There is no remaining freedom to make a gauge transformation in this scenario. Nevertheless, in this work, we will use the synchronous gauge, where we will choose ${ϵ}_0$ such that $E=0$ and ${ϵ}^S$ such that $H=0$ (similar conditions for ${\widetilde{ϵ}}_0$ and ${\widetilde{ϵ}}^S$). In the next section, we present the perturbed equations of motion in this frame, and we discuss the suitability of this choice for our purposes.

3.4. Fields Equations and Energy Momentum Conservations in Synchronous Gauge

Under this gauge fixing, the perturbed Einstein Equation (11) reads (at first order)4:

$$-4\pi G(\delta \rho +3\delta p+{\nabla }^2{\pi }^S)=\frac{1}{2}\left(3\ddot{A}+{\nabla }^2\ddot{B}\right)+\frac{\dot{a}}{2a}\left(3\dot{A}+{\nabla }^2\dot{B}\right).$$

(81)

where ${\nabla }^2\equiv {\partial }_x^2+{\partial }_y^2+{\partial }_z^2$. The energy-momentum conservation gives

$$\begin{array}{ccc}\hfill \delta p+{\nabla }^2{\pi }^S+{\partial }_0\left[(\overline{\rho }+\overline{p})\delta u\right]+\frac{3\dot{a}}{a}(\overline{\rho }+\overline{p})\delta u& =& 0,\hfill \end{array}$$

(82)

$$\begin{array}{ccc}\hfill \delta \dot{\rho }+\frac{3\dot{a}}{a}(\delta \rho +\delta p)+{\nabla }^2\left[a^{-2}(\overline{\rho }+\overline{p})\delta u+\frac{\dot{a}}{a}{\pi }^S\right]+\frac{1}{2}(\overline{\rho }+\overline{p}){\partial }_0\left[3A+{\nabla }^{2}B\right]& =& 0.\hfill \end{array}$$

(83)

We define

$$\Psi \equiv \frac{1}{2}\left[3A+{\nabla }^{2}B\right],$$

(84)

then,

$$-4\pi G{a}^2(\delta \rho +3\delta p+{\nabla }^2{\pi }^S)=\frac{\partial }{\partial t}\left(a^2\dot{\Psi }\right),$$

(85)

$$\delta \dot{\rho }+\frac{3\dot{a}}{a}(\delta \rho +\delta p)+{\nabla }^2\left[a^{-2}(\overline{\rho }+\overline{p})\delta u+\frac{\dot{a}}{a}{\pi }^S\right]+\frac{1}{2}(\overline{\rho }+\overline{p})\dot{\Psi }=0.$$

(86)

The unperturbed Einstein equations correspond to the Friedmann equations. In the delta sector, computations give the non-perturbed equations:

$$3\dot{F}\frac{\dot{a}}{a}=\varkappa (3F\overline{\rho }+\tilde{\overline{\rho }})$$

(87)

and

$$12F\frac{\ddot{a}}{a}+6F{\left(\frac{\dot{a}}{a}\right)}^2+3\dot{F}\frac{\dot{a}}{a}-3\ddot{F}=\varkappa (\tilde{\overline{\rho }}+3\tilde{\overline{p}}+3F\overline{\rho }+3F\overline{p}).$$

(88)

The perturbed contribution (at first order) is

$$\begin{array}{c}\hfill \left[2\dot{F}\frac{\dot{a}}{a}+\ddot{F}\right]\left[3A+{\nabla }^{2}B\right]+\left[6F\frac{\dot{a}}{a}+\frac{5}{2}\dot{F}\right]\left[3\dot{A}+{\nabla }^2\dot{B}\right]-\left[2\frac{\dot{a}}{a}\right]\left[3\dot{\tilde{A}}+{\nabla }^2\dot{\tilde{B}}\right]\\ \hfill +3F\left[3\ddot{A}+{\nabla }^2\ddot{B}\right]-\left[3\ddot{\tilde{A}}+{\nabla }^2\ddot{\tilde{B}}\right]=\varkappa \left(3\delta \tilde{p}+\delta \tilde{\rho }+F\delta \rho +3F\delta p+{\nabla }^2\tilde{\pi }+F{\nabla }^2\pi \right)\end{array}$$

(89)

In addition, the 00 component of delta energy-momentum conservation in (14) gives

$$\begin{array}{c}\hfill \delta \dot{\tilde{\rho }}+\frac{3\dot{a}}{a}(\delta \tilde{\rho }+\delta \tilde{p})+\frac{3\dot{F}}{2}(\delta \rho +\delta p)+{\nabla }^2\left[\frac{(\tilde{\overline{\rho }}+\tilde{\overline{p}})}{a^2}\delta u+\frac{(\overline{\rho }+\overline{p})F}{a^2}\delta u+\frac{(\overline{\rho }+\overline{p})}{a^2}\delta \tilde{u}\right]\\ \hfill +\frac{(\tilde{\overline{\rho }}+\tilde{\overline{p}})}{2}{\partial }_0[3A+{\nabla }^{2}B]+\frac{(\overline{\rho }+\overline{p})}{2}{\partial }_0[3\tilde{A}+{\nabla }^2\tilde{B}]-\frac{(\overline{\rho }+\overline{p})}{2}{\partial }_0\left(F[3A+{\nabla }^{2}B]\right)=0,\end{array}$$

(90)

while the $i0$ component gives

$$\begin{array}{c}\hfill \delta \tilde{p}+{\partial }_0\left[(\tilde{\overline{\rho }}+\tilde{\overline{p}})\delta u\right]+{\partial }_0\left[(\overline{\rho }+\overline{p})\delta \tilde{u}\right]-{\partial }_0\left[(\overline{\rho }+\overline{p})F\delta u\right]+3(\overline{\rho }+\overline{p})\dot{F}\delta u\\ \hfill +\frac{3\dot{a}}{a}(\overline{\rho }+\overline{p})\delta \tilde{u}+\frac{3\dot{a}}{a}(\tilde{\overline{\rho }}+\tilde{\overline{p}})\delta u-\frac{3\dot{a}}{a}F(\overline{\rho }+\overline{p})\delta u=0.\end{array}$$

(91)

Analogous to the standard sector, we define

$$\tilde{\Psi }\equiv \frac{1}{2}\left[3\tilde{A}+{\nabla }^2\tilde{B}\right],$$

(92)

then the gravitational equation becomes

$$\begin{array}{c}\hfill \left[2\dot{F}\frac{\dot{a}}{a}+\ddot{F}\right]a^2\Psi +\left[6F\frac{\dot{a}}{a}+\frac{5}{2}\dot{F}\right]a^2\dot{\Psi }+3F{a}^2\ddot{\Psi }-\frac{d}{dt}\left(a^2\dot{\tilde{\Psi }}\right)=\frac{\varkappa }{2}\left(3\delta \tilde{p}+\delta \tilde{\rho }+F\delta \rho \right.\\ \hfill +\left.3F\delta p+{\nabla }^2\tilde{\pi }+F{\nabla }^2\pi \right).\end{array}$$

(93)

Now, the delta energy conservation is given by

$$\begin{array}{c}\hfill \delta \dot{\tilde{\rho }}+\frac{3\dot{a}}{a}(\delta \tilde{\rho }+\delta \tilde{p})+\frac{3\dot{F}}{2}(\delta \rho +\delta p)+{\nabla }^2\left[\frac{(\tilde{\overline{\rho }}+\tilde{\overline{p}})}{a^2}\delta u+\frac{(\overline{\rho }+\overline{p})F}{a^2}\delta u+\frac{(\overline{\rho }+\overline{p})}{a^2}\delta \tilde{u}\right]\\ \hfill +(\tilde{\overline{\rho }}+\tilde{\overline{p}})\dot{\Psi }+(\overline{\rho }+\overline{p})\dot{\tilde{\Psi }}-(\overline{\rho }+\overline{p}){\partial }_0(F\Psi )=0.\end{array}$$

(94)

The study of the non-perturbed sector was already treated in Alfaro et al. and applied to the supernovae observations [37,42]. We will consider these results when necessary. For now, we only need the expression for the time-dependent function $F\left(t\right)$, which is

$$F\left(Y\right)=-\frac{LY}{3}\sqrt{Y+C},$$

(95)

where $Y\equiv Y\left(t\right)=a\left(t\right)/a_0$ is the quotient between the scale factor at a time t over the scale factor in the actuality (which for our purposes, we will consider equal to one). L$(\sim 0.45)$ and C$(\sim {10}^{-4})$ are the new parameters of DG that are already determined by supernova data [37,42]. We have to remark that our definition of $\Psi$ is not the usual, since the standard definition is with the time derivative of fields A and B, respectively. In the delta sector, the combinations of these fields appear explicitly without a time derivative, so if the reader wants to compare results with other works, he or she should take into consideration this definition to analyze the gauge. In the next section, we will discuss the evolution of the cosmological fluctuations, which will help us compute the scalar contribution to the CMB.

4. Evolution of Cosmological Fluctuations

Until now, we have developed the perturbation theory in DG; now, we are interested in studying the evolution of the cosmological fluctuations to have a physical interpretation of the delta matter fields, which this theory naturally introduces. Even in the standard cosmology, the system of equations that describes these perturbations are complicated to allow analytic solutions, and there are comprehensive computer programs for this task, such as CMBfast [43,44] and CAMB [45]. However, such computer programs cannot give a clear understanding of the physical phenomena involved. Nevertheless, some good approximations allow computing the spectrum of the CMB fluctuations with a rather good agreement with these computer programs [46,47]. In particular, we are going to extend the Weinberg approach for this task. This method consists of two main aspects: first, the hydrodynamic limit, which assumes that near recombination time photons were in local thermal equilibrium with the baryonic plasma; then, photons could be treated hydro-dynamically, such as plasma and cold dark matter. Second, a sharp transition from thermal equilibrium to complete transparency at the moment $t_{ls}$ of the last scattering.

Since we will reproduce this approach, we consider the Universe’s standard components, which means photons, neutrinos, baryons, and cold dark matter. Then, the task is to understand the role of their own delta counterpart. We will also neglect both anisotropic inertia tensors and took the usual state equation for pressures and energy densities and perturbations. As we will treat photons and delta photons hydrodynamically, we will use $\delta u_{\gamma }=\delta u_B$ and $\delta {\tilde{u}}_{\gamma }=\delta {\tilde{u}}_B$. Finally, as the synchronous scheme does not completely fix the gauge freedom, one can use the remaining freedom to put $\delta u_D=0$, which means that cold dark matter evolves at rest with respect to the Universe expansion. In our theory, the extended synchronous scheme also has extra freedom, which we will use to choose $\delta {\tilde{u}}_D=0$ as its standard part. Now, we will present the equations for both sectors. However, we will provide more detail in the delta sector because Weinberg [47] already calculates the solution of Einstein’s equations. Einstein’s equations and its energy-momentum conservation in Fourier space are5

$$\begin{array}{ccc}\hfill \frac{d}{dt}\left(a^2{\dot{\Psi }}_q\right)& =& -4\pi G{a}^2\left(\delta {\rho }_{Dq}+\delta {\rho }_{Bq}+2\delta {\rho }_{\gamma q}+2\delta {\rho }_{\nu q}\right),\hfill \end{array}$$

(96)

$$\begin{array}{ccc}\hfill \delta {\dot{\rho }}_{\gamma q}+4H\delta {\rho }_{\gamma q}-(4q/3a){\overline{\rho }}_{\gamma }\delta u_{\gamma q}& =& -(4/3){\overline{\rho }}_{\gamma }{\dot{\Psi }}_q,\hfill \end{array}$$

(97)

$$\begin{array}{ccc}\hfill \delta {\dot{\rho }}_{Dq}+3H\delta {\rho }_{Dq}& =& -{\overline{\rho }}_D{\dot{\Psi }}_q,\hfill \end{array}$$

(98)

$$\begin{array}{ccc}\hfill \delta {\dot{\rho }}_{Bq}+3H\delta {\rho }_{Bq}-(q/a){\overline{\rho }}_B\delta u_{\gamma q}& =& -{\overline{\rho }}_B{\dot{\Psi }}_q,\hfill \end{array}$$

(99)

$$\begin{array}{ccc}\hfill \delta {\dot{\rho }}_{\nu q}+4H\delta {\rho }_{\nu q}-(4q/3a){\overline{\rho }}_{\nu }\delta u_{\nu q}& =& -(4/3){\overline{\rho }}_{\nu }{\dot{\Psi }}_q,\hfill \end{array}$$

(100)

where $H\equiv \dot{a}/a$. It is useful to rewrite these equations in term of the dimensionless fractional perturbation

$${\delta }_{\alpha q}=\frac{\delta {\rho }_{\alpha q}}{{\overline{\rho }}_{\alpha }+{\overline{p}}_{\alpha }},$$

(101)

where $\alpha$ can be $\gamma$, $\nu$, B and D (photons, neutrinos, baryons and dark matter, respectively). $a^4{\overline{\rho }}_{\gamma }$, $a^4{\overline{\rho }}_{\nu }$, $a^3{\overline{\rho }}_D$, $a^3{\overline{\rho }}_B$ are time-independent quantities; then, (96)–(100) are

$$\begin{array}{ccc}\hfill \frac{d}{dt}\left(a^2{\dot{\Psi }}_q\right)& =& -4\pi G{a}^2\left({\overline{\rho }}_D{\delta }_{Dq}+{\overline{\rho }}_B{\delta }_{Bq}+\frac{8}{3}{\overline{\rho }}_{\gamma }{\delta }_{\gamma q}+\frac{8}{3}{\overline{\rho }}_{\nu }\delta {}_{\nu q}\right),\hfill \end{array}$$

(102a)

$$\begin{array}{ccc}\hfill {\dot{\delta }}_{\gamma q}-(q^2/a^2)\delta u_{\gamma q}& =& -{\dot{\Psi }}_q,\hfill \end{array}$$

(102b)

$$\begin{array}{ccc}\hfill {\dot{\delta }}_{Dq}& =& -{\Psi }_q,\hfill \end{array}$$

(102c)

$$\begin{array}{ccc}\hfill {\dot{\delta }}_{Bq}-(q^2/a^2)\delta u_{\gamma q}& =& -{\dot{\Psi }}_q,\hfill \end{array}$$

(102d)

$$\begin{array}{ccc}\hfill {\dot{\delta }}_{\nu q}-(q^2/a^2)\delta u_{\nu q}& =& -{\dot{\Psi }}_q,\hfill \end{array}$$

(102e)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\left(\frac{\left(1+R\right)\delta u_{\gamma q}}{a}\right)& =& -\frac{1}{3a}{\delta }_{\gamma q},\hfill \end{array}$$

(102f)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\left(\frac{\delta u_{\nu q}}{a}\right)& =& -\frac{1}{3a}{\delta }_{\nu q},\hfill \end{array}$$

(102g)

where $R=3{\overline{\rho }}_B/4{\overline{\rho }}_{\gamma }$. By the other side, in the delta sector, we will use a dimensionless fractional perturbation. However, this perturbation is defined as the delta transformation of (101)6,

$${\tilde{\delta }}_{\alpha q}\equiv \tilde{\delta }{\delta }_{\alpha q}=\frac{\delta {\tilde{\rho }}_{\alpha q}}{{\overline{\rho }}_{\alpha }+{\overline{p}}_{\alpha }}-\frac{{\tilde{\overline{\rho }}}_{\alpha }+{\tilde{\overline{p}}}_{\alpha }}{{\overline{\rho }}_{\alpha }+{\overline{p}}_{\alpha }}{\delta }_{\alpha q}.$$

(103)

In [37], they found

$$\frac{{\tilde{\overline{\rho }}}_R}{{\overline{\rho }}_R}=-2F\left(a\right)and\frac{{\tilde{\overline{\rho }}}_M}{{\overline{\rho }}_M}=-\frac{3}{2}F\left(a\right).$$

(104)

We will assume that this quotient holds for every component. In addition, using the result that $a^4{\tilde{\overline{\rho }}}_{\gamma }/F$, $a^4{\tilde{\overline{\rho }}}_{\nu }/F$, $a^3{\tilde{\overline{\rho }}}_D/F$, $a^3{\tilde{\overline{\rho }}}_B/F$ are time independent, the equations for the delta sector are

$$\begin{array}{c}\hfill \left[2\dot{F}\frac{\dot{a}}{a}+\ddot{F}\right]a^2{\Psi }_q+\left[6F\frac{\dot{a}}{a}+\frac{5}{2}\dot{F}\right]a^2{\dot{\Psi }}_q+3F{a}^2{\ddot{\Psi }}_q-\frac{d}{dt}\left(a^2{\dot{\tilde{\Psi }}}_q\right)\\ \hfill =\frac{\varkappa }{2}{a}^2\left[{\overline{\rho }}_D{\tilde{\delta }}_{Dq}+{\overline{\rho }}_B{\tilde{\delta }}_{Bq}+\frac{8}{3}{\overline{\rho }}_{\gamma }{\tilde{\delta }}_{\gamma q}+\frac{8}{3}{\overline{\rho }}_{\nu }{\tilde{\delta }}_{\nu q}-\frac{F}{2}\left({\overline{\rho }}_D{\delta }_{Dq}+{\overline{\rho }}_B{\delta }_{Bq}\right)\right.& -& \left.\frac{8}{3}F\left({\overline{\rho }}_{\gamma }{\delta }_{\gamma q}+{\overline{\rho }}_{\nu }{\delta }_{\nu q}\right)\right],\hfill \end{array}$$

(105a)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\delta }}}_{\gamma q}-\frac{q^2}{a^2}\left(\delta {\tilde{u}}_{\gamma q}+F\delta u_{\gamma q}\right)+{\dot{\tilde{\Psi }}}_q-{\partial }_0\left(F{\Psi }_q\right)& =& 0,\hfill \end{array}$$

(105b)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\delta }}}_{Dq}+{\dot{\tilde{\Psi }}}_q-{\partial }_0\left(F{\Psi }_q\right)& =& 0,\hfill \end{array}$$

(105c)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\delta }}}_{Bq}-\frac{q^2}{a^2}\left(\delta {\tilde{u}}_{\gamma q}+F\delta u_{\gamma q}\right)+{\dot{\tilde{\Psi }}}_q-{\partial }_0\left(F{\Psi }_q\right)& =& 0,\hfill \end{array}$$

(105d)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\delta }}}_{\nu q}-\frac{q^2}{a^2}\left(\delta {\tilde{u}}_{\nu q}+F\delta u_{\nu q}\right)+{\dot{\tilde{\Psi }}}_q-{\partial }_0\left(F{\Psi }_q\right)& =& 0,\hfill \end{array}$$

(105e)

$$\begin{array}{ccc}\frac{{\tilde{\delta }}_{\gamma q}}{3a}+\frac{d}{dt}\left(\frac{(1+R)\delta {\tilde{u}}_{\gamma q}}{a}\right)+2F\frac{d}{dt}\left(\frac{(R-\tilde{R})\delta u_{\gamma q}}{a}\right)\hfill & & \\ \hfill -F\frac{d}{dt}\left(\frac{(1+R)\delta u_{\gamma q}}{a}\right)-2\dot{F}(\tilde{R}-R)\frac{\delta u_{\gamma q}}{a}& =& 0,\hfill \end{array}$$

(105f)

$$\begin{array}{ccc}\hfill \frac{{\tilde{\delta }}_{\nu q}}{3a}+\frac{d}{dt}\left(\frac{\delta {\tilde{u}}_{\nu q}}{a}\right)-F\frac{d}{dt}\left(\frac{\delta u_{\nu q}}{a}\right)& =& 0,\hfill \end{array}$$

(105g)

with $\tilde{R}=3{\tilde{\overline{\rho }}}_B/4{\tilde{\overline{\rho }}}_{\gamma }$. Due to the definition of tilde fractional perturbation (103), solutions for (105a)–(105g) can be obtained easily, putting all solutions of GR equal to zero; then, the system is exactly equal to the system of Equations (102a)–(102g) and the solutions of tilde perturbations in the homogeneous system are exactly equal to the GR solutions. Then, we only need to “turn on” the GR source and find the complete solutions just like a forced-system. We will impose initial conditions to find solutions valid up to recombination time. At sufficiently early times, the Universe was dominated by radiation, and as Friedmann equations are valid in our theory (in particular the first equation), we can use a good approximation given by $a\propto \sqrt{t}$ and $8\pi G{\overline{\rho }}_R/3=1/4t^2$, while R and $\tilde{R}\ll 1$. Here

$${\overline{\rho }}_M\equiv {\overline{\rho }}_D+{\overline{\rho }}_B,{\overline{\rho }}_R\equiv {\overline{\rho }}_{\gamma }+{\overline{\rho }}_{\nu }.$$

(106)

We are interested in adiabatic solutions in the sense that all the ${\delta }_{\alpha q}$ and ${\tilde{\delta }}_{\alpha q}$ become equal at very early times. So, we make the ansatz:

$$\begin{array}{c}\hfill {\delta }_{\gamma q}={\delta }_{\nu q}={\delta }_{Bq}={\delta }_{Dq}={\delta }_q,\delta u_{\gamma q}=\delta u_{\nu q}=\delta u_q,\end{array}$$

(107)

$$\begin{array}{c}\hfill {\tilde{\delta }}_{\gamma q}={\tilde{\delta }}_{\nu q}={\tilde{\delta }}_{Bq}={\tilde{\delta }}_{Dq}={\tilde{\delta }}_q,\delta {\tilde{u}}_{\gamma q}=\delta {\tilde{u}}_{\nu q}=\delta {\tilde{u}}_q.\end{array}$$

(108)

Finally, we drop the term $q^2/a^2$ because we are considering very early times. Then, Equations (102a)–(102g) become

$$\begin{array}{ccc}\hfill \frac{d}{dt}\left(t{\Psi }_q\right)& =& -\frac{1}{t}{\delta }_q,\hfill \end{array}$$

(109)

$$\begin{array}{ccc}\hfill {\dot{\delta }}_q& =& -{\Psi }_q,\hfill \end{array}$$

(110)

and

$$\frac{d}{dt}\left(\frac{\delta u_q}{\sqrt{t}}\right)=-\frac{1}{t}{\delta }_q.$$

(111)

While Equations (105a)–(105g) become

$$\begin{array}{c}\hfill \left[2\dot{F}\frac{\dot{a}}{a}+\ddot{F}\right]a^2{\Psi }_q+\left[6F\frac{\dot{a}}{a}+\frac{5}{2}\dot{F}\right]a^2{\dot{\Psi }}_q+3F{a}^2{\ddot{\Psi }}_q\end{array}$$

$$\begin{array}{ccc}\hfill -\frac{d}{dt}\left(a^2{\dot{\tilde{\Psi }}}_q\right)& =& \frac{a^2}{t^2}({\tilde{\delta }}_q-F{\delta }_q),\hfill \end{array}$$

(112)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\delta }}}_q+{\dot{\tilde{\Psi }}}_q-{\partial }_0\left(F{\Psi }_q\right)& =& 0,\hfill \end{array}$$

(113)

$$\begin{array}{ccc}\hfill \frac{{\tilde{\delta }}_q}{3a}+\frac{d}{dt}\left(\frac{\delta {\tilde{u}}_q}{a}\right)-F\frac{d}{dt}\left(\frac{\delta u_q}{a}\right)& =& 0.\hfill \end{array}$$

(114)

An inspection of Equation (95) shows that at this era, for $a\ll C$, we have $F\propto -La\sqrt{C}/3$. In addition, in DG, time can be integrated from the first Friedmann equation with only radiation and matter, and one obtains:

$$t\left(Y\right)=\frac{2\sqrt{1+C}}{3H_0}\left(\sqrt{Y+C}(Y-2C)+2C^{\frac{3}{2}}\right),$$

(115)

We recall that $Y=a/a_0=a$ assuming $a_0=1$, $H_0={\dot{a}}_0/a_0$ is the usual Hubble parameter which we recall is no longer the physical Hubble parameter. Thus, radiation era time and $a\left(t\right)$ were related by $a\left(t\right)={(3H_0\sqrt{C}/\sqrt{1+C})}^{1/2}{t}^{1/2}$. This complete system consisting of Equations (109)–(111) and Equations (112)–(114) has an analytical solution:

$${\delta }_{\gamma q}={\delta }_{Bq}={\delta }_{Dq}={\delta }_{\nu q}=\frac{q^2{t}^2{\mathcal{R}}_q}{a^2},$$

(116)

$${\dot{\Psi }}_q=-\frac{t{q}^2{\mathcal{R}}_q}{a^2},$$

(117)

$$\delta u_{\gamma q}=\delta u_{\nu q}=-\frac{2t^3{q}^2{\mathcal{R}}_q}{9a^2},$$

(118)

where7

$$q^2{\mathcal{R}}_q\equiv -a^{2}H{\Psi }_q+4\pi G{a}^2\delta {\rho }_q+q^{2}H\delta u_q,$$

(119)

is a gauge-invariant quantity, which take a time-independent value for $q/a\ll H$. Here, $H=\dot{a}/a$ is the GR definition of the Hubble parameter, which we recall is no longer the physical one. On the other hand, we obtain

$$\begin{array}{ccc}\hfill {\tilde{\delta }}_q& =& -\frac{L\sqrt{C}{q}^2{\mathcal{R}}_q{t}^2}{3a},\hfill \end{array}$$

(120)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\Psi }}}_q& =& \frac{L\sqrt{C}{q}^2{\mathcal{R}}_{q}t}{a},\hfill \end{array}$$

(121)

$$\begin{array}{ccc}\hfill \delta {\tilde{u}}_q& =& \frac{L\sqrt{C}{q}^2{\mathcal{R}}_q{t}^3}{a}.\hfill \end{array}$$

(122)

We will talk about these initial conditions later. Note that Equations (102b)–(102d) give

$$\frac{d}{dt}\left({\delta }_{\gamma }-{\delta }_B\right)=0.$$

(123)

This implies that if we start from adiabatic solutions, ${\delta }_{\gamma }={\delta }_B$ is true for all the Universe evolution (the same happens for its delta version from Equations (105b)–(105d)).

Matter Era

In this era, we use $a\propto t^{2/3}$, then (still using $R=\tilde{R}=0$), we have

$$\frac{d}{dt}\left(a^2{\Psi }_q\right)=-4\pi G{\overline{\rho }}_D{a}^2{\delta }_{Dq},$$

(124a)

$${\dot{\delta }}_{Dq}=-{\Psi }_q,$$

(124b)

$$\frac{d}{dt}\left(\frac{(1+R)\delta u_{\gamma q}}{a}\right)=-\frac{1}{3a}{\delta }_{\gamma q},$$

(124c)

$$\frac{d}{dt}\left(\frac{\delta u_{\nu q}}{a}\right)=-\frac{1}{3a}{\delta }_{\nu q}.$$

(124d)

For the delta sector,

$$\begin{array}{c}\hfill \left[2\dot{F}\frac{\dot{a}}{a}+\ddot{F}\right]a^2{\Psi }_q+\left[6F\frac{\dot{a}}{a}+\frac{5}{2}\dot{F}\right]a^2{\dot{\Psi }}_q+3F{a}^2{\ddot{\Psi }}_q,\end{array}$$

$$\begin{array}{c}\hfill -\frac{d}{dt}\left(a^2{\dot{\tilde{\Psi }}}_q\right)=\frac{2a^2}{3t^2}\left({\tilde{\delta }}_{Dq}-F\frac{{\delta }_{Dq}}{2}\right),\end{array}$$

(124e)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\delta }}}_{\gamma q}-\frac{q^2}{a^2}\left(\delta {\tilde{u}}_{\gamma q}+F\delta u_{\gamma q}\right)+{\dot{\tilde{\Psi }}}_q-{\partial }_0\left(F{\Psi }_q\right)& =& 0,\hfill \end{array}$$

(124f)

$$\begin{array}{ccc}\hfill {\dot{\tilde{\delta }}}_{Dq}+{\dot{\tilde{\Psi }}}_q-{\partial }_0\left(F{\Psi }_q\right)& =& 0,\hfill \end{array}$$

(124g)

$$\begin{array}{ccc}\hfill \frac{{\tilde{\delta }}_{\gamma q}}{3a}+\frac{d}{dt}\left(\frac{\delta {\tilde{u}}_{\gamma q}}{a}\right)-F\frac{d}{dt}\left(\frac{\delta u_{\gamma q}}{a}\right)& =& 0.\hfill \end{array}$$

(124h)

where (in this era),

$$a\left(t\right)={\left(\frac{3H_0}{2\sqrt{1+C}}\right)}^{2/3}{t}^{2/3},$$

(125)

$$\begin{array}{ccc}\hfill F\left(t\right)& \propto & -\frac{L}{3}a{\left(t\right)}^{3/2}.\hfill \end{array}$$

(126)

It is remarkable that in the GR sector, there are exact solutions given by

$$\begin{array}{ccc}\hfill {\delta }_{Dq}& =& \frac{9q^2{t}^2{\mathcal{R}}_q\mathcal{T}\left(\kappa \right)}{10a^2},\hfill \end{array}$$

(127)

$$\begin{array}{ccc}\hfill {\dot{\Psi }}_q& =& -\frac{3q^{2}t{\mathcal{R}}_q\mathcal{T}\left(\kappa \right)}{5a^2},\hfill \end{array}$$

(128)

$$\begin{array}{c}\hfill {\delta }_{\gamma q}={\delta }_{\nu q}=\frac{3{\mathcal{R}}_q}{5}\left[\mathcal{T}\left(\kappa \right)-\mathcal{S}\left(\kappa \right)cos\left(q{\int }_0^t\frac{dt}{\sqrt{3}a}+\Delta \left(\kappa \right)\right)\right],\end{array}$$

(129)

$$\begin{array}{c}\hfill \delta u_{\gamma q}=\delta u_{\nu q}=\frac{3t{\mathcal{R}}_q}{5}\left[-\mathcal{T}\left(\kappa \right)+\mathcal{S}\left(\kappa \right)\frac{a}{\sqrt{3}qt}sin\left(q{\int }_0^t\frac{dt}{\sqrt{3}a}+\Delta \left(\kappa \right)\right)\right].\end{array}$$

(130)

where $\mathcal{T}\left(\kappa \right)$, $\mathcal{S}\left(\kappa \right)$ and $\Delta \left(\kappa \right)$ are time-independent dimensionless functions of the dimensionless re-scaled wave number

$$\kappa \equiv \frac{q\sqrt{2}}{a_{EQ}{H}_{EQ}}$$

(131)

$a_{EQ}$ and $H_{EQ}$ are, respectively, the Robertson–Walker scale factor and the expansion rate at matter radiation equally. These are referred to as transfer functions. (These functions can only depend on $\kappa$ because they need to be independent of spatial coordinates’ normalization and are dimensionless. A comprehensive analysis of the behavior of these functions can be found in [47].) Conversely, delta perturbations do not have an exact solution, and numerical calculations are necessary to determine them. However, in this work, we will not present numerical solutions and instead focus on estimating the initial conditions of the perturbations at the end of this section.

In order to obtain all transfer functions, we have to compare solutions with the full equation system (with ${\rho }_B={\tilde{\rho }}_B=0$). To do this task, let us make the change of variable $y\equiv a/a_{EQ}=a/C$; this means

$$\frac{d}{dt}=\frac{H_{EQ}}{\sqrt{2}}\frac{\sqrt{1+y}}{y}\frac{d}{dy}.$$

(132)

In addition, we will use the following parametrization for all perturbations

$$\begin{array}{c}\hfill {\delta }_{Dq}={\kappa }^2{\mathcal{R}}_q^{0}d\left(y\right)/4,{\delta }_{\gamma q}={\delta }_{\nu q}={\kappa }^2{\mathcal{R}}_q^{0}r\left(y\right)/4,\\ \hfill {\dot{\Psi }}_q=({\kappa }^2{H}_{EQ}/4\sqrt{2}){\mathcal{R}}_q^{0}f\left(y\right),\delta u_{\gamma q}=\delta u_{\nu q}=({\kappa }^2\sqrt{2}/4H_{EQ}){\mathcal{R}}_q^{0}g\left(y\right),\end{array}$$

and

$$\begin{array}{c}\hfill {\tilde{\delta }}_{Dq}={\kappa }^2{\mathcal{R}}_q^0\tilde{d}\left(y\right)/4,{\tilde{\delta }}_{\gamma q}={\tilde{\delta }}_{\nu q}={\kappa }^2{\mathcal{R}}_q^0\tilde{r}\left(y\right)/4,\\ \hfill {\dot{\tilde{\Psi }}}_q=({\kappa }^2{H}_{EQ}/4\sqrt{2}){\mathcal{R}}_q^0\tilde{f}\left(y\right),\delta {\tilde{u}}_{\gamma q}=\delta {\tilde{u}}_{\nu q}=({\kappa }^2\sqrt{2}/4H_{EQ}){\mathcal{R}}_q^0\tilde{g}\left(y\right).\end{array}$$

Then, Equations (124a)–(124d) and Equations (124e)–(124h) become

$$\begin{array}{ccc}\hfill \sqrt{1+y}\frac{d}{dy}\left(y^{2}f\left(y\right)\right)& =& -\frac{3}{2}d\left(y\right)-\frac{4r\left(y\right)}{y},\hfill \end{array}$$

(133a)

$$\begin{array}{ccc}\hfill \sqrt{1+y}\frac{d}{dy}r\left(y\right)-\frac{{\kappa }^{2}g\left(y\right)}{y}& =& -yf\left(y\right),\hfill \end{array}$$

(133b)

$$\begin{array}{ccc}\hfill \sqrt{1+y}\frac{d}{dy}d\left(y\right)& =& -yf\left(y\right),\hfill \end{array}$$

(133c)

$$\begin{array}{ccc}\hfill \sqrt{1+y}\frac{d}{dy}\left(\frac{g\left(y\right)}{y}\right)& =& -\frac{r\left(y\right)}{3},\hfill \end{array}$$

(133d)

and

$$\begin{array}{c}\hfill -\left[(1+2y)y{F}^{\prime }\left(y\right)+y(1+y)F^{″}\left(y\right)\right]d\left(y\right)+\left[6F\left(y\right)+\frac{5}{2}y{F}^{\prime }\left(y\right)\right]y\sqrt{1+y}f\left(y\right)\\ \hfill +3F\left(y\right)y^2\sqrt{1+y}{f}^{\prime }\left(y\right)-\sqrt{1+y}\frac{d}{dy}\left(y^2\tilde{f}\left(y\right)\right)=\frac{3\tilde{d}\left(y\right)}{2}+\frac{4\tilde{r}\left(y\right)}{y}\\ \hfill -\frac{3F\left(y\right)d\left(y\right)}{4}-\frac{4F\left(y\right)r\left(y\right)}{y}\end{array},$$

(133e)

$$\begin{array}{ccc}\hfill \sqrt{1+y}\frac{d}{dy}\tilde{d}\left(y\right)=-y\tilde{f}\left(y\right)-\sqrt{1+y}\frac{d}{dy}d\left(y\right)& & ,\hfill \end{array}$$

(133f)

$$\begin{array}{ccc}\hfill \sqrt{1+y}\frac{d}{dy}\tilde{r}\left(y\right)=\frac{{\kappa }^2}{y}[\tilde{g}\left(y\right)+F\left(y\right)g\left(y\right)]-y\tilde{f}\left(y\right)-\sqrt{1+y}\frac{d}{dy}d\left(y\right)& & ,\hfill \end{array}$$

(133g)

$$\begin{array}{ccc}\hfill \sqrt{1+y}\frac{d}{dy}\left(\frac{\tilde{g}\left(y\right)}{y}\right)=-\frac{\tilde{r}\left(y\right)}{3}+\sqrt{1+y}F\left(y\right)\frac{d}{dy}\left(\frac{g\left(y\right)}{y}\right)& & .\hfill \end{array}$$

(133h)

In this notation, the initial conditions are

$$\begin{array}{ccc}\hfill d\left(y\right)& =& r\left(y\right)\to y^2,\hfill \\ \hfill f\left(y\right)& \to & -2,\hfill \\ \hfill g\left(y\right)& \to & -\frac{y^4}{9}.\hfill \end{array}$$

For the delta sector,

$$\begin{array}{ccc}\hfill \tilde{d}\left(y\right)& =& \tilde{r}\left(y\right)\to -\frac{L{C}^{3/2}}{3}{y}^3,\hfill \\ \hfill \tilde{f}\left(y\right)& \to & \sqrt{2}L{C}^{3/2}y,\hfill \\ \hfill \tilde{g}\left(y\right)& \to & \frac{L{C}^{3/2}}{2}{y}^5.\hfill \end{array}$$

From supernovae fit, we know that $C\sim {10}^{-4}$ and $L\sim 0.45$ [37,42]; thus, we can estimate that fluctuations of “delta matter” at the beginning of the Universe were much smaller than fluctuations of standard matter. For example, at $y\sim {10}^{-3}$, the ratio between components of the Universe is $|{\tilde{\delta }}_{\alpha }/{\delta }_{\alpha }|\sim {10}^{-10}$.

We do not show numerical solutions here because the aim of this work is to trace a guide for future work, in particular, in the numeric computation of multipole coefficients for temperature fluctuations in the CMB. However, we will derive the equations to calculate that computation.

5. Derivation of Temperature Fluctuations

It is possible to find expressions analogous to temperature fluctuations usually obtained by Boltzmann equations by studying photons propagation in FRLW-perturbed coordinates, with the condition ${\overline{g}}_{i0}=0$8. For DG, the metric which photons follow is given by

$$\begin{array}{ccc}\hfill {\mathbf{g}}_{00}& =& -((1+3F\left(t\right))+E(\mathbf{x},t)+\tilde{E}(\mathbf{x},t)),{\mathbf{g}}_{i0}=0,\hfill \\ \hfill {\mathbf{g}}_{ij}& =& a^2\left(t\right)(1+F\left(t\right)){\delta }_{ij}+h_{ij}(\mathbf{x},t)+{\tilde{h}}_{ij}(\mathbf{x},t),\hfill \end{array}$$

(134)

A ray of light propagating to the origin of the FRLW coordinate system, from a direction $\widehat{n}$, will have a comoving radial coordinate r related with t by

$$\begin{array}{c}\hfill 0={\overline{\mathbf{g}}}_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-(1+3F\left(t\right)+E(r\widehat{n},t)+\tilde{E}(r\widehat{n},t))d{t}^2\\ \hfill +(a^2\left(t\right)(1+F\left(t\right))+h_{rr}(r\widehat{n},t)+{\tilde{h}}_{rr}(r\widehat{n},t))d{r}^2,\end{array}$$

(135)

in other words,

$$\begin{array}{ccc}\hfill \frac{dr}{dt}& =& -{\left(\frac{(1+3F\left(t\right))+E+\tilde{E}}{a^2\left(t\right)(1+F\left(t\right))+h_{rr}+{\tilde{h}}_{rr}}\right)}^{1/2}\hfill \\ & \simeq & -\frac{1}{a_{DG}\left(t\right)}+\frac{(h_{rr}+{\tilde{h}}_{rr})}{2(1+3F\left(t\right))a_{DG}^3\left(t\right)}-\frac{E+\tilde{E}}{2(1+3F\left(t\right))a_{DG}\left(t\right)},\hfill \end{array}$$

(136)

where $a_{DG}\left(t\right)$ is the modified scale factor given by

$$a_{DG}\left(t\right)=a\left(t\right)\sqrt{\frac{1+F\left(t\right)}{1+3F\left(t\right)}}.$$

(137)

Now, we will use the approximation of a sharp transition between the opaque and transparent Universe at a moment $t_{ls}$ of last scattering at red shift $z\simeq 1090$. With this approximation, the relevant term at the first order in Equation (136) is

$$r\left(t\right)=\left[s\left(t\right)+{\int }_{t_{ls}}^t\frac{d{t}^{\prime }}{a_{DG}\left(t^{\prime }\right)}N\left(s\left(t^{\prime }\right)\widehat{n},t^{\prime }\right)\right],$$

(138)

where

$$N(\mathbf{x},t)\equiv \frac{1}{2(1+3F)}\left[\frac{h_{rr}(\mathbf{x},t)+{\tilde{h}}_{rr}(\mathbf{x},t)}{a_{DG}^2}-E(\mathbf{x},t)-\tilde{E}(\mathbf{x},t)\right],$$

(139)

and $s\left(t\right)$ is the zero-order solution for the radial coordinate. $s\left(t\right)=r_{ls}$ when $t=t_{ls}$:

$$s\left(t\right)=r_{ls}-{\int }_{t_{ls}}^t\frac{d{t}^{\prime }}{a_{DG}\left(t^{\prime }\right)}={\int }_t^{t_0}\frac{d{t}^{\prime }}{a_{DG}\left(t^{\prime }\right)}.$$

(140)

If a ray of light arrives to $r=0$ at a time $t_0$, then Equation (138) gives

$$0=s\left(t_0\right)+{\int }_{t_{ls}}^t\frac{d{t}^{\prime }}{a_{DG}\left(t^{\prime }\right)}N\left(s\left(t^{\prime }\right)\widehat{n},t^{\prime }\right)=r_{ls}+{\int }_{t_{ls}}^{t_0}\frac{dt}{a_{DG}\left(t\right)}\left(N\left(s\left(t\right)\widehat{n},t\right)-1\right).$$

(141)

A time interval $\delta t_{ls}$, between the departure of successive rays of light at time $t_{ls}$ of last scattering, produces an interval of time $\delta t_0$ between the arrival of the rays of light at $t_0$, which are given by the variation of Equation (141):

$$\begin{array}{c}\hfill 0=\frac{\delta t_{ls}}{a_{DG}\left(t_{ls}\right)}\left[1-N(r_{ls}\widehat{n},t_{ls})+{\int }_{t_{ls}}^{t_0}\frac{dt}{a_{DG}\left(t\right)}{\left(\frac{\partial N\left(r\left(t\right)\widehat{n},t\right)}{\partial r}\right)}_{r=s\left(t\right)}\right]\\ \hfill +\delta t_{ls}(\partial u_{\gamma }^r(r_{ls}\widehat{n},t_{ls})+\partial {\tilde{u}}_{\gamma }^r(r_{ls}\widehat{n},t_{ls}))+\frac{\delta t_0}{a_{DG}\left(t_0\right)}\left[-1+N(0,t_0)\right].\end{array}$$

(142)

The velocity terms of the photon–gas or photon–electron–nucleon arise because of the variation with respect to the time of the radial coordinate $r_{ls}$ described by Equation (141). The exchange rate of $N(s\left(t\right)\widehat{n},t)$ is

$$\frac{d}{dt}N\left(s\left(t\right)\widehat{n},t\right)={\left(\frac{\partial }{\partial t}N(r\widehat{n},t)\right)}_{r=s\left(t\right)}-\frac{1}{a_{DG}\left(t\right)}{\left(\frac{\partial N(r\widehat{n},t)}{\partial r}\right)}_{r=s\left(t\right)},$$

then,

$$\begin{array}{c}\hfill 0=\frac{\delta t_{ls}}{a_{DG}\left(t_{ls}\right)}\left[1-N(0,t_{ls})+{\int }_{t_{ls}}^{t_0}dt{\left(\frac{\partial N\left(r\widehat{n},t\right)}{\partial t}\right)}_{r=s\left(t\right)}\right]\\ \hfill +\delta t_{ls}(\partial u_{\gamma }^r(r_{ls}\widehat{n},t_{ls})+\partial {\tilde{u}}_{\gamma }^r(r_{ls}\widehat{n},t_{ls}))+\frac{\delta t_0}{a_{DG}\left(t_0\right)}\left[-1+N(0,t_0)\right].\end{array}$$

(143)

This result gives the ratio between the time intervals between ray of lights that are emitted and received. However, we are interested in this ratio but for the proper time, which in DG is defined with the original metric $g_{\mu \nu }$:

$$\delta {\tau }_L=\sqrt{1+E(r_{ls},t_{ls})}\delta t_{ls},\delta {\tau }_0=\sqrt{1+E(0,t_0)}\delta t_0,$$

(144)

At first order, it gives the ratio between a received frequency and an emitted one:

$$\begin{array}{c}\hfill \frac{{\nu }_0}{{\nu }_L}=\frac{\delta {\tau }_L}{\delta {\tau }_0}=\frac{a_{DG}\left(t_{ls}\right)}{a_{DG}\left(t_0\right)}\left[1+\frac{1}{2}\left(E(r_{ls}\widehat{n},t)-E(0,t_0)\right)\right.\\ \hfill -\left.{\int }_{t_{ls}}^{t_0}{\left(\frac{\partial }{\partial t}N(r\widehat{n},t)\right)}_{r=s\left(t\right)}dt-a_{DG}\left(t\right)(\delta u_{\gamma }^r(r_{ls}\widehat{n},t)+\delta {\tilde{u}}_{\gamma }^r(r_{ls}\widehat{n},t))\right].\end{array}$$

(145)

In [42], we defined the physical scale factor as $Y_{DG}\left(t\right)\equiv a_{DG}\left(t\right)/a_{DG}\left(t_0\right)$. Thus, we recover the standard expression for the redshift. The observed temperature at the present time $t_0$ from direction $\widehat{n}$ is

$$T\left(\widehat{n}\right)=\left(\frac{{\nu }_0}{{\nu }_L}\right)(\overline{T}\left(t_{ls}\right)+\delta T(r_{ls}\widehat{n},t_{ls})),$$

(146)

In the absence of perturbations, the observed temperature in all directions should be

$$T_0=\left(\frac{a_{DG}\left(t_{ls}\right)}{a_{DG}\left(t_0\right)}\right)\overline{T}\left(t_{ls}\right),$$

(147)

Therefore, the ratio between the observed temperature shift that comes from direction $\widehat{n}$ and the unperturbed value is

$$\begin{array}{ccc}\hfill \frac{\Delta T\left(\widehat{n}\right)}{T_0}& \equiv & \frac{T\left(\widehat{n}\right)-T_0}{T_0}=\frac{{\nu }_0{a}_{DG}\left(t_0\right)}{{\nu }_L{a}_{DG}\left(t_{ls}\right)}-1+\frac{\delta T(r_{ls}\widehat{n},t_{ls})}{\overline{T}\left(t_{ls}\right)}\hfill \\ & =& \frac{1}{2}\left(E(r_{ls}\widehat{n},t)-E(0,t_0)\right)-{\int }_{t_{ls}}^{t_0}dt{\left(\frac{\partial }{\partial t}N(r\widehat{n},t)\right)}_{r=s\left(t\right)}\hfill \\ & -& a_{DG}\left(t\right)(\delta u_{\gamma }^r(r_{ls}\widehat{n},t)+\delta {\tilde{u}}_{\gamma }^r(r_{ls}\widehat{n},t))+\frac{\delta T(r_{ls}\widehat{n},t_{ls})}{\overline{T}\left(t_{ls}\right)}.\hfill \end{array}$$

(148)

For scalar perturbations in any gauge with ${\mathbf{h}}_{i0}=0$, the metric perturbations are

$$\begin{array}{c}\hfill h_{00}=-E,h_{ij}=(1+F)a^2\left[A{\delta }_{ij}+\frac{{\partial }^{2}B}{\partial x^i\partial x^j}\right],\\ \hfill {\tilde{h}}_{00}=-\tilde{E},{\tilde{h}}_{ij}=(1+F)a^2\left[\tilde{A}{\delta }_{ij}+\frac{{\partial }^2\tilde{B}}{\partial x^i\partial x^j}\right].\end{array}$$

(149)

In addition, for scalar perturbations, the radial velocity of the photon fluid and the delta versions are given in terms of the velocity potentials $\delta u_{\gamma }$ and $\delta {\tilde{u}}_{\gamma }$, respectively,

$$\begin{array}{c}\hfill \delta u_{\gamma }^r={(\overline{g}+\overline{\tilde{g}})}^{r\mu }\frac{\partial \delta u_{\gamma }}{\partial x^{\mu }}=\frac{1}{(1+F\left(t\right))a^2}\frac{\partial \delta u_{\gamma }}{\partial r},\\ \hfill \delta {\tilde{u}}_{\gamma }^r={(\overline{g}+\overline{\tilde{g}})}^{r\mu }\frac{\partial \delta {\tilde{u}}_{\gamma }}{\partial x^{\mu }}=\frac{1}{(1+F\left(t\right))a^2}\frac{\partial \delta {\tilde{u}}_{\gamma }}{\partial r}.\end{array}$$

(150)

Then, Equation (148) gives the scalar contribution to temperature fluctuations

$$\begin{array}{ccc}\hfill {\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}^S& =& \frac{1}{2}\left(E(r_{ls}\widehat{n},t)-E(0,t_0)\right)-{\int }_{t_{ls}}^{t_0}dt{\left(\frac{\partial }{\partial t}N(r\widehat{n},t)\right)}_{r=s\left(t\right)}\hfill \\ & -& \frac{1}{(1+3F\left(t\right))a_{DG}}\left(\frac{\partial \delta u_{\gamma }(r_{ls}\widehat{n},t)}{\partial r}+\frac{\partial \delta {\tilde{u}}_{\gamma }(r_{ls}\widehat{n},t)}{\partial t}\right)\hfill \\ & +& \frac{\delta T(r_{ls}\widehat{n},t_{ls})}{\overline{T}\left(t_{ls}\right)},\hfill \end{array}$$

(151)

where

$$N=\frac{1}{2}\left[A+\frac{{\partial }^{2}B}{\partial r^2}+\left(\tilde{A}+\frac{{\partial }^2\tilde{B}}{\partial r^2}\right)-\frac{E}{1+3F}-\frac{\tilde{E}}{1+3F}\right].$$

(152)

In the next step, we will study the gauge transformations of these fluctuations. The following identity for the fields B and $\tilde{B}$ will be useful:

$${\left(\frac{{\partial }^2\dot{B}}{\partial r^2}\right)}_{r=s\left(t\right)}=-{\left(\frac{d}{dt}\left[a_{DG}\frac{\partial \dot{B}}{\partial r}+a_{DG}{\dot{a}}_{DG}\dot{B}+a_{DG}^2\ddot{B}\right]+\frac{\partial }{\partial t}\left[a_{DG}{\dot{a}}_{DG}\dot{B}+a_{DG}^2\ddot{B}\right]\right)}_{r=s\left(t\right)}.$$

(153)

Then, the temperature fluctuations are described by

$${\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}^S={\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{early}^S+{\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{late}^S+{\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{ISW}^S$$

(154)

where

$$\begin{array}{ccc}\hfill {\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{early}^S& =& -\frac{1}{2}{a}_{DG}\left(t_{ls}\right){\dot{a}}_{DG}\left(t_{ls}\right)\dot{B}(r_{ls}\widehat{n},t_{ls})-\frac{1}{2}{a}_{DG}^2\left(t_{ls}\right)\ddot{B}(r_{ls}\widehat{n},t_{ls})+\frac{1}{2}E(r_{ls}\widehat{n},t_{ls})\hfill \\ \hfill +\frac{\delta T\left(r_{ls}\widehat{n}\right)}{\overline{T}\left(t_{ls}\right)}& -& a_{DG}\left(t_{ls}\right)\left[\frac{\partial }{\partial r}{\left(\frac{1}{2}\dot{B}(r\widehat{n},t_{ls})+\frac{1}{(1+3F\left(t_{ls}\right))a_{DG}^2\left(t_{ls}\right)}\delta u_{\gamma }(r\widehat{n},t_{ls})\right)}_{r=r_{ls}}\right]\hfill \\ & -& \left\{\left(\frac{1}{2}{a}_{DG}\left(t_{ls}\right){\dot{a}}_{DG}\left(t_{ls}\right)\dot{\tilde{B}}(r_{ls}\widehat{n},t_{ls})+\frac{1}{2}{a}_{DG}^2\left(t_{ls}\right)\ddot{\tilde{B}}(r_{ls}\widehat{n},t_{ls})\right)\right.\hfill \\ \hfill +a_{DG}\left(t_{ls}\right)& \times & \left.\left[\frac{\partial }{\partial r}{\left(\frac{1}{2}\dot{\tilde{B}}(r\widehat{n},t_{ls})+\frac{1}{(1+3F\left(t_{ls}\right))a_{DG}^2\left(t_{ls}\right)}\delta {\tilde{u}}_{\gamma }(r\widehat{n},t_{ls})\right)}_{r=r_{ls}}\right]\right\},\hfill \end{array}$$

(155)

$$\begin{array}{ccc}\hfill {\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{late}^S& =& \frac{1}{2}{a}_{DG}\left(t_0\right){\dot{a}}_{DG}\left(t_0\right)\dot{B}(0,t_0)+\frac{1}{2}{a}_{DG}^2\left(t_0\right)\ddot{B}(0,t_0)-\frac{1}{2}E(0,t_0)\hfill \\ & +& a_{DG}\left(t_0\right)\left[\frac{\partial }{\partial r}{\left(\frac{1}{2}\dot{B}(r\widehat{n},t_0)+\frac{1}{(1+3F\left(t_0\right))a_{DG}^2\left(t_0\right)}\delta u_{\gamma }(r\widehat{n},t_0)\right)}_{r=0}\right]\hfill \\ & +& \left\{\left(\frac{1}{2}{a}_{DG}\left(t_0\right){\dot{a}}_{DG}\left(t_0\right)\dot{\tilde{B}}(0,t_0)+\frac{1}{2}{a}_{DG}^2\left(t_0\right)\ddot{\tilde{B}}(0,t_0)\right)\right.\hfill \\ \hfill +a_{DG}\left(t_0\right)& \times & \left.\left[\frac{\partial }{\partial r}{\left(\frac{1}{2}\dot{\tilde{B}}(r\widehat{n},t_0)+\frac{1}{(1+3F\left(t_0\right))a_{DG}^2\left(t_0\right)}\delta {\tilde{u}}_{\gamma }(r\widehat{n},t_0)\right)}_{r=r_{ls}}\right]\right\},\hfill \end{array}$$

(156)

$$\begin{array}{ccc}\hfill {\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{ISW}^S& =& -\frac{1}{2}{\int }_{t_{ls}}^{t_0}dt\left\{\frac{\partial }{\partial t}\left[a_{DG}^2\left(t\right)\ddot{B}(r\widehat{n},t)+a_{DG}\left(t\right){\dot{a}}_{DG}\left(t\right)\dot{B}(r\widehat{n},t)+A(r\widehat{n},t)\right.\right.\hfill \\ & -& \frac{E(r\widehat{n},t)}{1+3F\left(t\right)}+a_{DG}^2\left(t\right)\ddot{\tilde{B}}(r\widehat{n},t)+a_{DG}\left(t\right){\dot{a}}_{DG}\left(t\right)\dot{\tilde{B}}(r\widehat{n},t)\hfill \\ & +& \left.\left.\tilde{A}(r\widehat{n},t)-\frac{\tilde{E}(r\widehat{n},t)}{1+3F\left(t\right)}\right]\right\},\hfill \end{array}$$

(157)

The “late” term is the sum of independent direction terms and a term proportional to $\widehat{n}$, which was added to represent the local anisotropies of the gravitational field and the local fluid. In GR, these terms only contribute to the multipole expansion for $l=0$ and $l=1$. Thus, we will ignore their contribution to our derivation of the temperature fluctuations multipoles coefficients.

5.1. Gauge Transformations

We are going to study the gauge transformations for photons propagating in the metric ${\mathbf{g}}_{\mu \nu }$ for a parameter ${\mathit{ϵ}}_{\mu }$. Then, the transformations are

$$\begin{array}{ccc}\hfill \Delta A& =& \frac{2\dot{a}}{(1+F)a}\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F},\Delta B=-\frac{2}{1+F}\frac{{\mathit{ϵ}}^{\mathbf{S}}}{(1+F)a^2},\hfill \\ \hfill \Delta C_i& =& -\frac{1}{1+F}\frac{{\mathit{ϵ}}_{\mathbf{i}}^{\mathbf{V}}}{(1+F)a^2},\Delta D_{ij}=0,\Delta E=2\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F}\right),\hfill \\ \hfill \Delta H& =& -\frac{1}{\sqrt{1+F}a}\left[a^2\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}^{\mathbf{S}}}{(1+F)a^2}\right)+\frac{{\mathit{ϵ}}_{\mathbf{0}}}{(1+3F)}\right],\Delta G_i=-\frac{a}{\sqrt{1+F}}\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}_{\mathbf{i}}^{\mathbf{V}}}{(1+F)a^2}\right).\hfill \end{array}$$

(158)

and

$$\begin{array}{ccc}\hfill \Delta \tilde{A}& =& \frac{1}{(1+F)a^2}\left[\frac{\partial }{\partial t}\left(F{a}^2\right)\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F}\right],\Delta \tilde{B}=-\frac{1}{(1+F)a^2}\left[\frac{2F}{1+F}{\mathit{ϵ}}^{\mathbf{S}}\right],\hfill \\ \hfill \Delta {\tilde{C}}_i& =& -\frac{F}{1+F}\frac{{\mathit{ϵ}}_{\mathbf{i}}^{\mathbf{V}}}{(1+F)a^2},\Delta {\tilde{D}}_{ij}=0,\Delta \tilde{E}=6F\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F}\right)+\frac{3\dot{F}}{1+3F}{\mathit{ϵ}}_{\mathbf{0}}\hfill \\ \hfill \Delta \tilde{H}& =& -\frac{1}{\sqrt{1+F}a}\left[F{a}^2\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}^{\mathbf{S}}}{(1+F)a^2}\right)+\frac{3F{\mathit{ϵ}}_{\mathbf{0}}}{(1+3F)}\right],\hfill \\ \hfill \Delta {\tilde{G}}_i& =& -\frac{1}{\sqrt{1+F}a}\left[F{a}^2\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}_{\mathbf{i}}^{\mathbf{V}}}{(1+F)a^2}\right)\right].\hfill \end{array}$$

(159)

Now, considering the sum of the perturbations, we obtain

$$\begin{array}{ccc}\hfill \Delta A+\Delta \tilde{A}& =& \frac{1}{(1+F)a^2}\frac{\partial }{\partial t}\left[(1+F)a^2\right]\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F},\hfill \end{array}$$

(160a)

$$\begin{array}{ccc}\hfill \Delta B+\Delta \tilde{B}& =& -\frac{2{\mathit{ϵ}}^{\mathbf{S}}}{(1+F)a^2},\hfill \end{array}$$

(160b)

$$\begin{array}{ccc}\hfill \Delta E+\Delta \tilde{E}& =& 2(1+3F)\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F}\right)+\frac{3\dot{F}}{1+3F}{\mathit{ϵ}}_{\mathbf{0}},\hfill \end{array}$$

(160c)

$$\begin{array}{ccc}\hfill \Delta H+\Delta \tilde{H}& =& -\frac{1}{\sqrt{1+F}a}\left[(1+F)a^2\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}^{\mathbf{S}}}{(1+F)a^2}\right)+{\mathit{ϵ}}_{\mathbf{0}}\right],\hfill \end{array}$$

(160d)

$$\begin{array}{ccc}\hfill \Delta C_i+\Delta {\tilde{C}}_i& =& -\frac{{\mathit{ϵ}}_{\mathbf{i}}^{\mathbf{V}}}{(1+F)a^2},\hfill \end{array}$$

(160e)

$$\begin{array}{ccc}\hfill \Delta G_i+\Delta {\tilde{G}}_i& =& -\frac{1}{\sqrt{1+F}a}\left[(1+F)a^2\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}_{\mathbf{i}}^{\mathbf{V}}}{(1+F)a^2}\right)\right].\hfill \end{array}$$

(160f)

Now, we will study the gauge transformations that preserve the condition ${\mathbf{g}}_{i0}=g_{io}+{\tilde{g}}_{i0}=0$. This means that $\Delta H+\Delta \tilde{H}=0$. This gives a solution for ${\mathit{ϵ}}_{\mathbf{0}}$ given by

$$\begin{array}{c}\hfill {\mathit{ϵ}}_{\mathbf{0}}=-(1+F)a^2\frac{\partial }{\partial t}\left(\frac{{\mathit{ϵ}}^{\mathbf{S}}}{(1+F)a^2}\right).\end{array}$$

(161)

When we study how the “ISW” term transforms under this type of transformation, we found that $\Delta ISW=0$. While for the “early” term, we should note that temperature perturbations transform as

$$\Delta \delta T(r_{ls}\widehat{n},t)=\dot{\overline{T}}\left(t\right)\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F},$$

(162)

With this expression and $\overline{T}{a}_{DG}=cte$, we finally obtain

$$\frac{\Delta \delta T(r_{ls}\widehat{n},t)}{\overline{T}\left(t_{ls}\right)}=-\frac{{\dot{a}}_{DG}}{a_{DG}}\frac{{\mathit{ϵ}}_{\mathbf{0}}}{1+3F}.$$

(163)

This results implies that the “early” term is invariant under this gauge transformation. Note that this gauge transformation is equivalent to the previously discussed in Section 1, because we can always take $\mathit{ϵ}$ as a combination of $ϵ$ and $\widetilde{ϵ}$. Then, we remark that temperature fluctuations are gauge invariant under scalar transformations that leave ${\mathbf{g}}_{i0}=0$.

5.2. Single Modes

We will assume that since the last scattering until now, all the scalar contributions are dominated by a unique mode such that any perturbation $X(\mathbf{x},t)$ could be written as

$$X(\mathbf{x},t)=\int d^{3}q\alpha \left(\mathbf{q}\right)e^{i\mathbf{q}·\mathbf{x}}{X}_q\left(t\right),$$

(164)

where $\alpha \left(\mathbf{q}\right)$ is a stochastic variable, which is normalized such that

$$\langle \alpha \left(\mathbf{q}\right){\alpha }^{\ast }\left({\mathbf{q}}^{\prime }\right)\rangle ={\delta }^3(\mathbf{q}-{\mathbf{q}}^{\prime }).$$

(165)

Then, Equations (155) and (157) become

$$\begin{array}{ccc}\hfill {\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{early}^S& =& \int d^{3}q\alpha \left(\mathbf{q}\right)e^{i\mathbf{q}·\widehat{n}r\left(t_{ls}\right)}\left(\mathcal{F}\left(q\right)+\tilde{\mathcal{F}}\left(q\right)+i\widehat{q}·\widehat{n}(\mathcal{G}\left(q\right)+\tilde{\mathcal{G}}\left(q\right))\right),\hfill \end{array}$$

(166)

$$\begin{array}{ccc}\hfill {\left(\frac{\Delta T\left(\widehat{n}\right)}{T_0}\right)}_{ISW}^S& =& -\frac{1}{2}{\int }_{t_0}^{t_1}dt\int d^{3}q\alpha \left(\mathbf{q}\right)e^{i\mathbf{q}·\widehat{n}s\left(t\right)}\frac{d}{dt}\left[a_{DG}^2\left(t\right){\ddot{B}}_q\left(t\right)+a_{DG}\left(t\right){\dot{a}}_{DG}\left(t\right){\dot{B}}_q\left(t\right)\right.\hfill \\ & +& A_q\left(t\right)-\frac{E_q\left(t\right)}{1+3F\left(t\right)}+\left(a_{DG}^2\left(t\right){\ddot{\tilde{B}}}_q\left(t\right)+a_{DG}\left(t\right){\dot{a}}_{DG}\left(t\right){\dot{\tilde{B}}}_q\left(t\right)+{\tilde{A}}_q\left(t\right)\right.\hfill \\ & -& \left.\left.\frac{{\tilde{E}}_q\left(t\right)}{1+3F\left(t\right)}\right)\right],\hfill \end{array}$$

(167)

where

$$\begin{array}{ccc}\hfill \mathcal{F}\left(q\right)& =& -\frac{1}{2}{a}_{DG}^2\left(t\right){\ddot{B}}_q\left(t_{ls}\right)-\frac{1}{2}{a}_{DG}\left(t\right){\dot{a}}_{DG}\left(t_{ls}\right){\dot{B}}_q\left(t_{ls}\right)+\frac{1}{2}{E}_q\left(t_{ls}\right)+\frac{\delta T_q\left(t_{ls}\right)}{\overline{T}\left(t_{ls}\right)},\hfill \end{array}$$

(168)

$$\begin{array}{ccc}\hfill \tilde{\mathcal{F}}\left(q\right)& =& -\frac{1}{2}{a}_{DG}^2\left(t\right){\ddot{\tilde{B}}}_q\left(t_{ls}\right)-\frac{1}{2}{a}_{DG}\left(t_{ls}\right){\dot{a}}_{DG}\left(t_{ls}\right){\dot{\tilde{B}}}_q\left(t_{ls}\right),\hfill \end{array}$$

(169)

$$\begin{array}{ccc}\hfill \mathcal{G}\left(q\right)& =& -q\left(\frac{1}{2}{a}_{DG}\left(t_{ls}\right){\dot{B}}_q\left(t_{ls}\right)+\frac{1}{(1+3F\left(t_{ls}\right))a_{DG}\left(t_{ls}\right)}\delta u_{\gamma }\left(t_{ls}\right)\right),\hfill \end{array}$$

(170)

$$\begin{array}{ccc}\hfill \tilde{\mathcal{G}}\left(q\right)& =& -q\left(\frac{1}{2}{a}_{DG}\left(t_{ls}\right){\dot{\tilde{B}}}_q\left(t_{ls}\right)+\frac{1}{(1+3F\left(t_{ls}\right))a_{DG}\left(t_{ls}\right)}\delta {\tilde{u}}_{\gamma }\left(t_{ls}\right)\right).\hfill \end{array}$$

(171)

These functions are called form factors. We emphasize that combination given by $\mathcal{F}\left(q\right)+\tilde{\mathcal{F}}\left(q\right)$ and $\mathcal{G}\left(q\right)+\tilde{\mathcal{G}}\left(q\right)$, and the expressions inside the integral are gauge invariants under gauge transformations that preserve ${\mathbf{g}}_{i0}$ equal to zero.

6. Coefficients of Multipole Temperature Expansion: Scalar Modes

As an application of the previous results, we will study the contribution of the scalar modes for temperature–temperature correlation, which is given by:

$$C_{TT,l}=\frac{1}{4\pi }\int d^2\widehat{n}\int d^2{\widehat{n}}^{\prime }{P}_l(\widehat{n}·{\widehat{n}}^{\prime })\langle \Delta T\left(\widehat{n}\right)\Delta T\left({\widehat{n}}^{\prime }\right)\rangle ,$$

(172)

where $\Delta T\left(\widehat{n}\right)$ is the stochastic variable which gives the deviation of the average of observed temperature in direction $\widehat{n}$, and $\langle \dots \rangle$ denotes the average over the position of the observer. However, the observed quantity is

$$C_{TT,l}^{obs}=\frac{1}{4\pi }\int d^2\widehat{n}\int d^2{\widehat{n}}^{\prime }{P}_l(\widehat{n}·{\widehat{n}}^{\prime })\Delta T\left(\widehat{n}\right)\Delta T\left(\widehat{n^{\prime }}\right),$$

(173)

Nevertheless, the mean square fractional difference between this equation and Equation (172) is $2/(2l+1)$, and therefore, it may be neglected for $l\gg 1$. In order to calculate these coefficients, we use the following expansion in spherical harmonics

$$e^{i\widehat{q}·\widehat{n}\rho }=4\pi \sum _{l=0}^{\infty }\sum _{m=-l}^{m=l}{i}^l{j}_l\left(\rho \right)Y_l^m\left(\widehat{n}\right)Y_l^{m\ast }\left(\widehat{q}\right),$$

(174)

where $j_l\left(\rho \right)$ represents the spherical Bessel’s functions. Using this expression in Equation (166), and replacing the factor $i\widehat{q}·\widehat{n}$ for time derivatives of Bessel’s functions, the scalar contribution of the observed T–T fluctuations in direction $\widehat{n}$ is

$${(\Delta T\left(\widehat{n}\right))}^S=\sum _{lm}{a}_{T,lm}^S{Y}_l^m\left(\widehat{n}\right),$$

(175)

where

$$\begin{array}{c}\hfill a_{T,lm}^S=4\pi i^l{T}_0\int d^{3}q\alpha \left(\mathbf{q}\right)Y_l^{l\ast }\left(\widehat{q}\right)\left[j_l\left(q{r}_{ls}\right)(\mathcal{F}\left(q\right)+\tilde{\mathcal{F}}\left(q\right))+j_l^{\prime }\left(q{r}_{ls}\right)(\mathcal{G}\left(q\right)+\tilde{\mathcal{G}}\left(q\right))\right],\end{array}$$

(176)

and $\alpha \left(\mathbf{q}\right)$ is a stochastic parameter for the dominant scalar mode. It is normalized such that

$$\langle \alpha \left(\mathbf{q}\right){\alpha }^{\ast }\left({\mathbf{q}}^{\prime }\right)\rangle ={\delta }^3(\mathbf{q}-{\mathbf{q}}^{\prime }).$$

(177)

Inserting this expression in Equation (172), we obtain

$$C_{TT,l}^S=16{\pi }^2{T}_0^2{\int }_0^{\infty }{q}^{2}dq{\left[j_l\left(q{r}_{ls}\right)(\mathcal{F}\left(q\right)+\tilde{\mathcal{F}}\left(q\right))+j_l^{\prime }\left(q{r}_{ls}\right)(\mathcal{G}\left(q\right)+\tilde{\mathcal{G}}\left(q\right))\right]}^2.$$

(178)

Now, we will consider the case $l\gg 1$. In this limit, we can use the following approximation for Bessel’s functions9:

$$\begin{array}{c}\hfill j_l\left(\rho \right)\to \left\{\begin{array}{cc}cos\left(b\right)cos\left[\nu (tanb-b)-\pi /4\right]/\left(\nu \sqrt{sinb}\right)& \rho >\nu ,\\ 0& \rho <\nu ,\end{array}\right.\end{array}$$

(179)

where $\nu \equiv l+1/2$, and $cosb\equiv \nu /\rho$, with $0\le b\le \pi /2$. In addition, for $\rho >\nu \gg 1$, the phase $\nu (tanb-b)$ is a function of $\rho$ that grows very fast; then, the derivatives of Bessel’s functions only act in its phase:

$$\begin{array}{c}\hfill j_l^{\prime }\left(\rho \right)\to \left\{\begin{array}{cc}-cos\left(b\right)\sqrt{sinb}sin\left[\nu (tanb-b)-\pi /4\right]/\nu & \rho >\nu ,\\ 0& \rho <\nu .\end{array}\right.\end{array}$$

(180)

Using these approximations in Equation (178) and changing the variable from q to $b={cos}^{-1}(\nu /q{r}_{ls})$, we obtain

$$\begin{array}{ccc}\hfill C_{TT,l}^S& =& \frac{16{\pi }^2{T}_0^2\nu }{r_{ls}^3}{\int }_0^{\pi /2}\frac{db}{{cos}^{2}b}\hfill \\ & & \times \left[\left(\mathcal{F}\left(\frac{\nu }{r_{ls}cosb}\right)+\tilde{\mathcal{F}}\left(\frac{\nu }{r_{ls}cosb}\right)\right)cos[\nu (tanb-b)-\pi /4]\right.\hfill \\ & & -{\left.sinb\left(\mathcal{G}\left(\frac{\nu }{r_{ls}cosb}\right)+\tilde{\mathcal{G}}\left(\frac{\nu }{r_{ls}cosb}\right)\right)sin[\nu (tanb-b)-\pi /4]\right]}^2.\hfill \end{array}$$

(181)

When $\nu \gg 1$, the functions $cos\left[\nu \right(tanb-b)-\pi /4]$ and $sin\left[\nu \right(tanb-b)-\pi /4]$ oscillate very rapidly; then, the squared average of its values are $1/2$, while the averaged cross-terms are zero. Using $l\approx \nu$, and changing the integration variable from b to $\beta =1/cosb$, Equation (181) becomes

$$\begin{array}{ccc}\hfill l(l+1)C_{TT,l}^S& =& \frac{8{\pi }^2{T}_0^2{l}^3}{r_{ls}^3}{\int }_1^{\infty }\frac{\beta d\beta }{\sqrt{{\beta }^2-1}}\hfill \\ & & \times \left[{\left(\mathcal{F}\left(\frac{l\beta }{r_{ls}}\right)+\tilde{\mathcal{F}}\left(\frac{l\beta }{r_{ls}}\right)\right)}^2+\frac{{\beta }^2-1}{{\beta }^2}{\left(\mathcal{G}\left(\frac{l\beta }{r_{ls}}\right)+\tilde{\mathcal{G}}\left(\frac{l\beta }{r_{ls}}\right)\right)}^2\right].\hfill \end{array}$$

(182)

Note that $d_A=r_{ls}{a}_{DG}\left(t_{ls}\right)$ is the angular diameter distance of the last scattering surface. To calculate the CMB power spectrum, we need to know the value of ${\dot{\tilde{B}}}_q$. We use the off-diagonal equation from the delta sector to obtain it. This gives:

$${\dot{\tilde{A}}}_q={\dot{A}}_{q}F+A_q\dot{F}-2a^2(\rho +p)\delta u_q-a^2(\tilde{\rho }+\tilde{p})\delta u_q-(\rho +p)\delta {\tilde{u}}_q,$$

(183)

so if we use this equation with the definition of $\tilde{\Psi }$

$${\dot{\tilde{\Psi }}}_q=\frac{1}{2}(3{\dot{\tilde{A}}}_q-q^2{\dot{\tilde{B}}}_q),$$

(184)

it allows us to find $\dot{\tilde{B}}$. Now, we will use the approximation that perturbations of a gravitation field are dominated by perturbations of dark matter density. In this regime ${\dot{A}}_q\left(t_{ls}\right)=0$ and in the synchronous gauge, the velocity perturbations for dark matter are zero; then,

$${\dot{\tilde{A}}}_q\left(t_{ls}\right)=A_q\left(t_{ls}\right)\dot{F}\left(t_{ls}\right),$$

(185)

and

$${\dot{\tilde{B}}}_q\left(t_{ls}\right)=\frac{3}{q^2}{A}_q\left(t_{ls}\right)\dot{F}\left(t_{ls}\right)-\frac{2{\dot{\tilde{\Psi }}}_q\left(t_{ls}\right)}{q^2}\Rightarrow {\ddot{\tilde{B}}}_q\left(t_{ls}\right)=\frac{3}{q^2}{A}_q\left(t_{ls}\right)\ddot{F}\left(t_{ls}\right)-\frac{2{\ddot{\tilde{\Psi }}}_q\left(t_{ls}\right)}{q^2},$$

(186)

where

$$\begin{array}{ccc}\hfill q^2{A}_q& =& 8\pi G{a}^2\delta {\rho }_{Dq}-2H{a}^2{\dot{\Psi }}_q\hfill \\ & =& 3H^2{a}^2{\delta }_{Dq}-2H{a}^2{\dot{\Psi }}_q.\hfill \end{array}$$

(187)

In GR, ${\dot{B}}_q=-2{\dot{\Psi }}_q/q^2$, and ${\dot{\Psi }}_q\propto t^{-1/3}$ implies ${\ddot{B}}_q=2{\dot{\Psi }}_q/3t{q}^2$. Therefore, the usual form factors are:

$$\begin{array}{ccc}\hfill \mathcal{F}\left(q\right)& =& \frac{1}{3}{\delta }_{\gamma q}\left(t_{ls}\right)+\frac{{\dot{\Psi }}_q\left(t_{ls}\right)}{q^2}\left(a_{DG}\left(t_{ls}\right){\dot{a}}_{DG}\left(t_{ls}\right)-\frac{2}{3}\frac{a_{DG}^2\left(t_{ls}\right)}{t_{ls}}\right),\hfill \end{array}$$

(188)

$$\begin{array}{ccc}\hfill \mathcal{G}\left(q\right)& =& -q\frac{\delta u_{\gamma q}\left(t_{ls}\right)}{(1+3F\left(t_{ls}\right))a_{DG}\left(t_{ls}\right)}+\frac{a_{DG}\left(t_{ls}\right){\dot{\Psi }}_q\left(t_{ls}\right)}{q}.\hfill \end{array}$$

(189)

where we have used $\delta T_q/\overline{T}=\delta {\rho }_{\gamma q}/4{\overline{\rho }}_{\gamma }={\delta }_{\gamma q}/3$. Nevertheless, for the “delta” contribution, ${\dot{\tilde{\Psi }}}_q$ and ${\ddot{\tilde{\Psi }}}_q$ satisfy the same relation as the standard case. Due to our decomposition, the tilde expresions are

$$\begin{array}{ccc}\hfill \tilde{\mathcal{F}}\left(q\right)& =& -\frac{3}{2}\frac{A_q\left(t_{ls}\right)}{q^2}(a_{DG}^2\left(t_{ls}\right)\ddot{F}\left(t_{ls}\right)+a_{DG}\left(t_{ls}\right){\dot{a}}_{DG}\left(t_{ls}\right)\dot{F}\left(t_{ls}\right))\hfill \\ & +& \frac{{\dot{\tilde{\Psi }}}_q\left(t_{ls}\right)}{q^2}\left(a_{DG}\left(t_{ls}\right){\dot{a}}_{DG}\left(t_{ls}\right)-\frac{2}{3}\frac{a_{DG}^2\left(t_{ls}\right)}{t_{ls}}\right),\hfill \end{array}$$

(190)

$$\begin{array}{ccc}\hfill \tilde{\mathcal{G}}\left(q\right)& =& -q\frac{\delta {\tilde{u}}_{\gamma q}\left(t_{ls}\right)}{(1+3F\left(t_{ls}\right))a_{DG}\left(t_{ls}\right)}+\frac{a_{DG}\left(t_{ls}\right){\dot{\tilde{\Psi }}}_q\left(t_{ls}\right)}{q}.\hfill \end{array}$$

(191)

Unfortunately, due to all the approximations we have used, we need to add some corrections to the solutions of the GR sector. After that, we will be able to find the numerical solutions for DG perturbations. The first consideration is that in the set of equations presented in the matter era, we have used $R=3{\overline{\rho }}_B/4{\rho }_{\gamma }=0$, which is not valid in this era. Corrections to the solutions can be calculated using a WKB approximation for perturbations10 [47]. The second consideration that we must include in the solution of photons perturbations is the so-called Silk damping11 [52,53], which takes into account the viscosity and heat conduction of the relativistic medium. Moreover, the transition from opaque to a transparent Universe at the last scattering moment was not instantaneous, but it could be considered a Gaussian. This effect is known as Landau damping12. We must recall that the physical geometry now is described by $Y_{DG}\left(t\right)=a_{DG}\left(t\right)/a_{DG}(t=0)$, so the expression for both Silk and Landau effects has to be expressed in this geometry. With these considerations, the solutions of perturbations are given by:

$$\begin{array}{ccc}\hfill {\dot{\Psi }}_q\left(t_{ls}\right)& =& -\frac{3q^2{t}_{ls}{\mathcal{R}}_q^o\mathcal{T}\left(\kappa \right)}{5a^2\left(t_{ls}\right)},\hfill \end{array}$$

(192)

$$\begin{array}{ccc}\hfill {\delta }_{\gamma q}\left(t_{ls}\right)& =& \frac{3{\mathcal{R}}_q^o}{5}\left[\mathcal{T}\left(\kappa \right)(1+3R_{ls})-{(1+R_{ls})}^{-1/4}{e}^{-q^2{d}_D^2/a_{ls}^2}\right.\hfill \\ & \times & \left.\mathcal{S}\left(\kappa \right)cos\left({\int }_0^{t_{ls}}\frac{qdt}{\sqrt{3(1+R(t\left)\right)}a\left(t\right)}+\Delta \left(\kappa \right)\right)\right],\hfill \end{array}$$

(193)

$$\begin{array}{ccc}\hfill \delta u_{\gamma q}\left(t_{ls}\right)& =& \frac{3{\mathcal{R}}_q^o}{5}\left[-t_{ls}\mathcal{T}\left(\kappa \right)+\frac{a\left(t_{ls}\right)}{\sqrt{3}q{(1+R_{ls})}^{3/4}}{e}^{-q^2{d}_D^2/a_{ls}^2}\right.\hfill \\ & \times & \left.\mathcal{S}\left(\kappa \right)sin\left({\int }_0^{t_{ls}}\frac{qdt}{\sqrt{3(1+R(t\left)\right)}a\left(t\right)}+\Delta \left(\kappa \right)\right)\right],\hfill \end{array}$$

(194)

Here, we used an approximation given by $a_{DG}\left(t_{ls}\right)\approx a\left(t_{ls}\right)\propto t^{2/3}$, and the error of this approximation is of the order ${10}^{-4}\%$.

$$\begin{array}{ccc}\hfill {\dot{\Psi }}_q\left(t_{ls}\right)& =& -\frac{3q^2{t}_{ls}{\mathcal{R}}_q^o\mathcal{T}\left(\kappa \right)}{5a_{DG}^2\left(t_{ls}\right)},\hfill \end{array}$$

(195)

$$\begin{array}{ccc}\hfill {\delta }_{\gamma q}\left(t_{ls}\right)& =& \frac{3{\mathcal{R}}_q^o}{5}\left[\mathcal{T}\left(\kappa \right)(1+3R_{ls})-{(1+R_{ls})}^{-1/4}{e}^{-q^2{d}_D^2/a_{DG}^2\left(t_{ls}\right)}\right.\hfill \\ & \times & \left.\mathcal{S}\left(\kappa \right)cos\left(q{\int }_0^{t_{ls}}\frac{dt}{\sqrt{3(1+R(t\left)\right)}{a}_{DG}\left(t\right)}+\Delta \left(\kappa \right)\right)\right],\hfill \end{array}$$

(196)

$$\begin{array}{ccc}\hfill \delta u_{\gamma q}\left(t_{ls}\right)& =& \frac{3{\mathcal{R}}_q^o}{5}\left[-t_{ls}\mathcal{T}\left(\kappa \right)+\frac{a_{DG}\left(t_{ls}\right)}{\sqrt{3}q{(1+R_{ls})}^{3/4}}{e}^{-q^2{d}_D^2/a_{DG}^2\left(t_{ls}\right)}\right.\hfill \\ & \times & \left.\mathcal{S}\left(\kappa \right)sin\left(q{\int }_0^{t_{ls}}\frac{dt}{\sqrt{3(1+R(t\left)\right)}{a}_{DG}\left(t\right)}+\Delta \left(\kappa \right)\right)\right],\hfill \end{array}$$

(197)

where

$$\begin{array}{c}\hfill d_D^2=d_{Silk}^2+d_{Landau}^2,\end{array}$$

(198)

$$\begin{array}{c}\hfill d_{Silk}^2=Y_{DG}^2\left(t_{ls}\right){\int }_0^{t_{ls}}\frac{t_{\gamma }}{6Y_{DG}^2(1+R)}\left\{\frac{16}{15}+\frac{R^2}{(1+R)}\right\}dt,\end{array}$$

(199)

$$\begin{array}{c}\hfill d_{Landau}^2=\frac{{\sigma }_t^2}{6(1+R_{ls})},\end{array}$$

(200)

where $t_{\gamma }$ is the mean free time for photons and $R=3{\overline{\rho }}_B/4{\overline{\rho }}_{\gamma }=3h^2{\Omega }_B{Y}_{DG}/4h^2{\Omega }_{\gamma }$. In order to evaluate the Silk damping, we have

$$t_{\gamma }=\frac{1}{n_e{\sigma }_T},$$

(201)

where $n_e$ is the number density of electrons and ${\sigma }_T$ is the Thomson cross-section.

On the other hand

$$\begin{array}{ccc}\hfill q{\int }_0^{r_{ls}}{c}_{s}dr& =& q{\int }_0^{t_{ls}}\frac{dt}{\sqrt{3(1+R(t\left)\right)}{a}_{DG}\left(t\right)}\equiv q{r}_{ls}^{SH}\hfill \\ & =& \frac{q}{a_{DG}\left(t_{ls}\right)}·\left(a_{DG}\left(t_{ls}\right)r_{ls}^{SH}\right)=\frac{q}{a_{DG}\left(t_{ls}\right)}·d_H\left(t_{ls}\right)\hfill \end{array}$$

(202)

where $c_s$ is the speed of sound, $r_{ls}^{SH}$ is the sound horizon radial coordinate and $d_H$ is the horizon distance.

With all this approximation, the transfers functions were simplified to the following expressions:

$$\begin{array}{ccc}\hfill \mathcal{F}\left(q\right)& =& \frac{1}{3}{\delta }_{\gamma q}\left(t_{ls}\right)+\frac{a_{DG}^2\left(t_{ls}\right){\dot{\Psi }}_q\left(t_{ls}\right)}{3q^2{t}_{ls}},\hfill \end{array}$$

(203)

$$\begin{array}{ccc}\hfill \mathcal{G}\left(q\right)& =& -q\frac{\delta u_{\gamma q}\left(t_{ls}\right)}{(1+3F\left(t_{ls}\right))a_{DG}\left(t_{ls}\right)}+\frac{a_{DG}\left(t_{ls}\right){\dot{\Psi }}_q\left(t_{ls}\right)}{q},\hfill \end{array}$$

(204)

where $A_q\left(t_{ls}\right)={\mathcal{R}}_q^o\mathcal{T}\left(\kappa \right)$. Then, we replaced the GR solutions, and we obtain

$$\begin{array}{ccc}\hfill \mathcal{F}\left(q\right)& =& \frac{{\mathcal{R}}_q^o}{5}\left[3\mathcal{T}(q{d}_T/a_{DG}\left(t_{ls}\right))R_{ls}-{(1+R_{ls})}^{-1/4}{e}^{-q^2{d}_D^2/a_{DG}^2\left(t_{ls}\right)}\right.\hfill \\ & \times & \left.\mathcal{S}(q{d}_T/a_{DG}\left(t_{ls}\right))cos\left(q{d}_H/a_{DG}\left(t_{ls}\right)+\Delta (q{d}_T/a_{DG}\left(t_{ls}\right))\right)\right],\hfill \end{array}$$

(205)

$$\begin{array}{ccc}\hfill \mathcal{G}\left(q\right)& =& \frac{\sqrt{3}{\mathcal{R}}_q^o}{5{(1+R_{ls})}^{3/4}}{e}^{-q^2{d}_D^2/a_{DG}^2\left(t_{ls}\right)}\hfill \\ & \times & \mathcal{S}(q{d}_T/a_{DG}\left(t_{ls}\right))sin\left(q{d}_H/a_{DG}\left(t_{ls}\right)+\Delta (q{d}_T/a_{DG}\left(t_{ls}\right))\right),\hfill \end{array}$$

(206)

where $\kappa =q{d}_T/a_{ls}$ (defined in Equation (131)) and

$$d_T\left(t_{ls}\right)\equiv \frac{\sqrt{2}{a}_{DG}\left(t_{ls}\right)}{a_{EQ}{H}_{EQ}}=\frac{a_{DG}\left(t_{ls}\right)\sqrt{{\Omega }_R}}{H_0{\Omega }_M}=\frac{a_{DG}\left(t_{ls}\right)}{100h}\sqrt{C(C+1)}.$$

(207)

The final consideration that we must include is that due to the reionization of hydrogen at $z_{reion}=10$ by ultraviolet light coming from the first generation of massive stars, photons of the CMB have a probability of being scattered $1-exp(-{\tau }_{reion})$. CMB has two contributions. The non-scattered photons provide the first contribution, where we have to correct by a factor given by $exp(-{\tau }_{reion})$. The scattered photons provide the second contribution, but the reionization occurs at $z\ll z_L$ affecting only low ls. We are not interested in this effect, and therefore, we will not include it. Measurements show that in GR, $exp(-2{\tau }_{reion})\approx 0.8$.

On the other hand, we will use a standard parametrization of ${\mathcal{R}}_q^0$ given by

$$|{\mathcal{R}}_q^0{|}^2=N^2{q}^{-3}{\left(\frac{q/R_0}{{\kappa }_{\mathcal{R}}}\right)}^{n_s-1},$$

(208)

where $n_s$ could vary with the wave number. It is usual to take ${\kappa }_{\mathcal{R}}=0.05$ Mpc${}^{-1}$.

Note that $d_A\left(t_{ls}\right)=r_{ls}{a}_{DG}\left(t_{ls}\right)$ is the angular diameter distance of the last scattering surface.

$$\begin{array}{ccc}\hfill d_A\left(t_{ls}\right)& =& a_{DG}\left(t_{ls}\right){\int }_{t_{ls}}^{t_0}\frac{d{t}^{\prime }}{a_{DG}\left(t^{\prime }\right)}=\frac{a_{DG}\left(t_0\right)}{1+z_{ls}}{\int }_{t_{ls}}^{t_0}\frac{d{t}^{\prime }}{a_{DG}\left(t^{\prime }\right)}=\frac{1}{1+z_{ls}}{\int }_{t_{ls}}^{t_0}\frac{d{t}^{\prime }}{Y_{DG}\left(t^{\prime }\right)}\hfill \\ & =& \frac{1}{1+z_{ls}}{\int }_{Y_{ls}}^1\frac{d{Y}^{\prime }}{Y_{DG}\left(Y^{\prime }\right)}\frac{dt}{d{Y}^{\prime }}=\frac{d_L\left(t_{ls}\right)}{{(1+z_{ls})}^2}.\hfill \end{array}$$

(209)

This is consistent with the luminosity distance definition [38]. Then, when we set $q=\beta l/r_{ls}$, we obtain

$$\begin{array}{ccc}\hfill |{\mathcal{R}}_{\beta l/r_{ls}}^0{|}^2& =& N^2{\left(\frac{\beta l}{r_{ls}}\right)}^{-3}{\left(\frac{\beta l}{{\kappa }_{\mathcal{R}}{r}_{ls}}\right)}^{n_s-1}=N^2{\left(\frac{\beta l}{r_{ls}}\right)}^{-3}{\left(\frac{\beta l{a}_{DG}\left(t_{ls}\right)}{{\kappa }_{\mathcal{R}}{r}_{ls}{a}_{DG}\left(t_{ls}\right)}\right)}^{n_s-1}\hfill \\ & =& N^2{\left(\frac{\beta l}{r_{ls}}\right)}^{-3}{\left(\frac{\beta l{a}_{DG}\left(t_{ls}\right)}{{\kappa }_{\mathcal{R}}{d}_A\left(t_{ls}\right)}\right)}^{n_s-1}\equiv N^2{\left(\frac{\beta l}{r_{ls}}\right)}^{-3}{\left(\frac{\beta l}{l_R}\right)}^{n_s-1}.\hfill \end{array}$$

Using similar computations for the other distances, the final form of the form factors is given by

$$\begin{array}{ccc}\hfill \mathcal{F}\left(q\right)& =& \frac{{\mathcal{R}}_q^o}{5}\left[3\mathcal{T}(\beta l/l_T)R_{ls}-{(1+R_{ls})}^{-1/4}{e}^{-{\beta }^2{l}^2/l_D^2}\right.\hfill \\ & \times & \left.\mathcal{S}(\beta l/l_T)cos\left(\beta l/l_H+\Delta (\beta l/l_T)\right)\right],\hfill \end{array}$$

(210)

$$\begin{array}{ccc}\hfill \mathcal{G}\left(q\right)& =& \frac{\sqrt{3}{\mathcal{R}}_q^o}{5{(1+R_{ls})}^{3/4}}{e}^{-{\beta }^2{l}^2/l_D^2}\mathcal{S}(\beta l/l_T)sin\left(\beta l/l_H+\Delta (\beta l/l_T)\right),\hfill \end{array}$$

(211)

where

$$l_R=\frac{{\kappa }_{\mathcal{R}}{d}_A\left(t_{ls}\right)}{a_{DG}\left(t_{ls}\right)},l_H=\frac{d_A\left(t_{ls}\right)}{d_H\left(t_{ls}\right)},l_T=\frac{d_A\left(t_{ls}\right)}{d_T\left(t_{ls}\right)},l_D=\frac{d_A\left(t_{ls}\right)}{d_D\left(t_{ls}\right)}.$$

(212)

To summarize, for reasonably large values of l (say $l>20$), CMB multipoles are given by

$$\begin{array}{ccc}\hfill \frac{l(l+1)C_{TT,l}^S}{2\pi }& =& \frac{4\pi T_0^2{l}^{3}exp(-2{\tau }_{reion})}{r_{ls}^3}{\int }_1^{\infty }\frac{\beta d\beta }{\sqrt{{\beta }^2-1}}\hfill \\ & \times & \left[{\left(F\left(\frac{l\beta }{r_{ls}}\right)+\tilde{F}\left(\frac{l\beta }{r_{ls}}\right)\right)}^2+\frac{{\beta }^2-1}{{\beta }^2}{\left(G\left(\frac{l\beta }{r_{ls}}\right)+\tilde{G}\left(\frac{l\beta }{r_{ls}}\right)\right)}^2\right].\hfill \end{array}$$

(213)

The structure of Equation (213) is remarkable, where the delta sector contributes additively inside the integral. If we set all the delta sector equal to zero, we recover the result directly for scalar temperature–temperature multipole coefficients in GR given by Weinberg. Numerical solutions are needed to compute the solution for the perturbations. In a preliminary numerical solution, we obtained the temperature power spectrum of the CMB following Weinberg’s approach [47]. We obtain:

$$\begin{array}{cc}\hfill C=5.25\times {10}^{-4},& {\Omega }_M=0.13038,{\Omega }_B=0.02228\end{array}$$

where ${\Omega }_M$ is the density of non-relativistic matter, and ${\Omega }_B$ is the baryon density. Figure 1 displays the temperature power spectrum for a specific set of cosmological parameters. These cosmological parameters were obtained through exploratory analysis without any statistical examination. In other words, preliminary values were tested to determine if there was any possibility of obtaining reasonable results for the temperature power spectrum using the delta gravity equations. Consequently, the plot does not incorporate error bars.

In the work conducted in ApJ [49], the shape of the temperature power spectrum is determined by five free parameters. These parameters were explored using a modified adaptive Metropolis MCMC algorithm. The statistical study, which aimed at determining the optimal parameter space to match the observed data of the temperature power spectrum, is complicated. The complexity arises from the involved equations of delta gravity, which encompass integrals that are computationally intensive, particularly when the MCMC algorithm performs numerous calculation cycles. Moreover, additional equations accounting for the Landau damping effect and other physical considerations during the last scattering epoch, specific to the delta gravity model, need to be included.

To overcome the computational challenges, an adaptive step method was employed for each parameter independently. Additionally, pre-generated interpolation tables for each integral involved in the calculation were utilized to reduce the computational time during each execution of the MCMC cycle. The endeavor of obtaining optimal parameters for the power spectrum involves an extensive study that significantly differs from this work in terms of physical considerations, statistical approach and numerical methods required to achieve that goal.

7. Conclusions

We discussed the implications of the first law of thermodynamics using the modified geometry of this model. We distinguished the physical densities from the GR densities in terms of which scale factor they dilute. However, knowing the solutions of the GR sector is enough for us to know about the behavior of the physical densities. In addition, if we consider that the number of photons is conserved after the moment of decoupling, the black body distribution should keep the form, and that means that temperature is redshifted with the modified scale factor $Y_{DG}$. Finally, we stated the ansatz that the moment of equality between radiation and matter was the same in GR and in DG, and we showed its implications in some parameters of the theory. We had developed the theory of perturbations for delta gravity and its gauge transformations. Following Weinberg [47], we used the synchronous gauge which leaves a residual gauge transformation, which can be used to set $\delta u_D=0$ (and also $\delta {\tilde{u}}_D=0$). Then, we computed the equations for cosmological perturbations using the hydrodynamic approximation, which we solved for the radiation era, while for a matter-dominated Universe, we presented the equations with the respective initial conditions. As in GR, we found an expression for temperature fluctuations in DG, studying the photon propagation in an effective metric from the moment of the last scattering until now. We found that those temperature fluctuations can be split into three independent terms: an early term which only depends on the moment of the last scattering $t_{ls}$, an ISW term that includes the evolution of gravitational fields from the last scattering to the present, and a late term which depends on the actual value for those fields. We compute the gauge transformations which leaves ${\mathbf{g}}_{i0}=0$, and we found that those three terms are separately gauge invariants. Then, we derived the TT multipole coefficients for scalar modes for large l ($l>200$), where we found that DG affects additively, which could have an observational effect that could be compared with Plank results and give a physical meaning for the so-called “delta matter”. As we mentioned in Section 2.1, differences between physicals and GR fluctuations will occur near the present (low z). Therefore, we should study the full scalar contribution of the multipole coefficients for low ls, where we expect relevant differences between both models. This task requires a full numerical derivation of perturbations, which is beyond the aim of this work. With the full scalar expression for the CMB power spectrum coefficients, we can find the shape of the spectrum. The full analysis and results can be found in [49]. The principal purpose of this work was to present the derivation of this expression and obtain physical insight into the cosmological fluctuations in delta gravity.