%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  New template code for TAMU Theses and Dissertations starting Fall 2016.  
%
%
%  Author: Sean Zachary Roberson
%  Version 3.17.09
%  Last Updated: 9/21/2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                          CHAPTER II
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\chapter{\uppercase{Conformally and Disformally coupled Scalar-Tensor Theory}\myfnt{Reprinted with permission from "Dark matter relics and the expansion rate in scalar-tensor theories'' by Bhaskar Dutta, Esteban Jimenez, Ivonne Zavala, 2017. PhysRevD.96.103506, Copyright 2017 by IOP Publishing Ltd and Sissa Medialab.}} \label{Ch:Conformal}


We are interested in scalar-tensor theories coupled to matter both conformally and disformally \cite{Bekenstein}. 
Our motivation comes from theories with extra dimensions and in particular string theory compactifications,
where several additional scalar fields appear, from closed and open string theory sectors of the theory \cite{KWZ}. 
Our approach in this work nonetheless will be phenomenological and therefore our equations will be simplified.  
However, we present the more general set-up, which can accommodate a realization from concrete string theory compactifications in  Appendix \ref{App1}. 

As was mentioned in chapter \ref{Ch:Intro}, ST theories are often formulated in one of two frames of reference, namely, the Jordan or Einstein frames.  The respective advantages of these two frames are that the scalar couplings enter through either the gravitational sector (Jordan frame) or the matter sector (Einstein frame), leaving the other sector unaffected. In the Jordan frame matter fields, $\Psi$, are coupled directly to the metric, $\tilde g_{\mu\nu}$, which means that the matter sector of the action can be written as $S_{Matter}=S_{Matter}(\tilde g_{\mu\nu},\Psi)$. Thus, this frame is more convenient for particle physics considerations because the usual observables, e.g. a mass, have their standard interpretation. However, the scalar field couples to the gravitational sector producing rather cumbersome gravitational field equations.

 On the other hand, in the Einstein frame, the matter piece of the action becomes  $S_{Matter}=S_{Matter}(g_{\mu\nu},\phi,\partial_\mu\phi,\Psi)$. This implies that physical quantities associated with particles (i.e. mass) measured in this frame have a spacetime dependency. However, the gravitational field equations take their standard form, where the Einstein tensor is proportional to the total energy-momentum tensor.  

The most common strategy followed in the literature ~\cite{Catena,MW} is to determine the cosmic evolution in the Einstein frame, where the cosmological equations take a more straightforward form, and then transform the results over to the Jordan frame. As was already hinted out before, the effect of modified gravity will enter the computation of particle physics processes through the expansion rate, $\tilde H$, in the Jordan frame.  

The action we  want to consider is given by: 
\bea\label{S1}
 S_{EH}=\frac{1}{2\kappa^2} \!\!\int{\!d^4x\sqrt{-g}\,R}
- \!\!\int{\!d^4x\sqrt{- g} \left[\frac{1}{2} (\partial\phi)^2+V(\phi)\right]} 
+S_{Matter} \,. 
\eea
%
where $S_{Matter}=- \!\!\int{\!d^4x\sqrt{-\tilde g} \,{\cal L}_{M}(\tilde g_{\mu\nu},\Psi) }$, with ${\cal L}_{M}$ the lagrangian density of the matter sector.

Here the disformally coupled metric is given by
%
\be \label{gtilde}
\tilde g_{\mu\nu} = C(\phi) g_{\mu\nu} + D(\phi) \partial_\mu \phi \partial_\nu \phi\,,
\ee
%
and the  inverse by:
\be\label{gtildein}
\tilde g^{\mu\nu} = \frac{1}{C}\left[ g^{\mu\nu} - \frac{D\,\partial^\mu \phi \partial^\nu \phi}{C+D(\partial\phi)^2}\right]\,.
\ee
%

Moreover, $\kappa^2=M_{P}^{-2}=8\pi G$, but keep in mind that $G$ is not in general equal to Newton's constant as measured by e.g.~local experiments. Further, $C(\phi), D(\phi)$ are functions of $\phi$, which can be identified as a conformal and disformal couplings of the scalar to the metric, respectively (note that the conformal coupling is dimensionless, whereas the disformal one has units of $mass^{-4}$).
%

The equations of motion obtained from  (\ref{S1}) are:
%
\be\label{EM}
R_{\mu\nu} -\frac{1}{2}g_{\mu\nu} R = \kappa^2\left(T^\phi_{\mu\nu} + T_{\mu\nu}\right)\,,
\ee
where, in the  frame  relative to $g_{\mu\nu}$, the energy-momentum tensors are defined as
% 
\bea\label{EM1}
T^\phi_{\mu\nu} = -\frac{2}{\sqrt{-g}} \frac{\delta S_{\phi}}{\delta g^{\mu\nu}}  \,,\qquad \quad %\\
T_{\mu\nu}= -\frac{2}{\sqrt{-g}} \frac{\delta \left(-\sqrt{-\tilde g} \,{\cal L}_{M}\right) }{\delta g^{\mu\nu}}\,,
\eea
% 
and we model the energy-momentum tensor for matter and both dark components as perfect fluids, that is:
\be\label{EM2}
T_{\mu\nu}^i = P_i g_{\mu\nu} +(\rho_i + P_i) u_\mu u_\nu\,,
\ee
where $\rho_i$, $P_i$ are the energy density and pressure for each  fluid $i$ with equation of state $P_i/\rho_i = \omega_i$. For the scalar field, the energy-momentum tensor takes the form:
%
\be\label{EMphi}
T_{\mu\nu}^{\phi} = - g_{\mu\nu} \left[ \frac{1}{2} (\partial\phi)^2+ V \right] 
+  \partial_\mu\phi \, \partial_\nu \phi\,,
\ee
%
and one can define the energy density and pressure of the scalar field as:
\bea\label{rhoP}
\rho_\phi =- \frac{1}{2}(\partial\phi)^2 +  V  \,, \qquad P_\phi = - \frac{1}{2}(\partial\phi)^2 - V  \,.
\eea
%
 
Finally the equation of motion for the scalar field dark energy becomes:
%
\bea\label{Eqphi}
&&\hskip-1cm - \nabla_\mu \nabla^\mu \phi \! + V' 
- \frac{T^{\mu\nu}}{2}\!\left[\frac{C'}{C} g_{\mu\nu}  +\frac{D'}{C}\partial_\mu\phi\partial_\nu\phi\right]
+\nabla_\mu \left[\frac{D}{C}T^{\mu\nu} \partial_\nu\phi \right] =0 \,.
\eea
%

%
 Due to the nontrivial coupling, the individual conservation equations for the two fluids are modified. However, the conservation equation for the full system is preserved and given in the usual way by
 % 
\be
\nabla_\mu \left(T^{\mu\nu}_{\phi} + T^{\mu\nu}\right) =0\,. 
\ee
%

Thus using (\ref{EMphi})  and the equation of motion for the scalar field we can write 
%
\be\label{conserva}
\nabla_\mu T^{\mu\nu}_\phi = Q \,\partial^\nu \phi = - \nabla_\mu T^{\mu\nu}\,,
\ee
where
\be
Q\equiv \nabla_\mu \left[\frac{D}{C} \,T^{\mu\lambda} \,\partial_\lambda \phi\right] - \frac{T^{\mu\nu} }{2} \left[\frac{C'}{C} g_{\mu\nu} +
\frac{D'}{C} \,\partial_\mu\phi \,\partial_\nu\phi\right]\,.
\ee

\bigskip 

%
In the Jordan, as defined above,  matter is conserved,
%
\be\label{conservaD}
\tilde \nabla_\mu \tilde T^{\mu\nu} =0  \,,
\ee
%
where $\tilde \nabla_\mu$ is the covariant derivative computed with respect to the disformal metric \eqref{gtilde} with  the Christoffel symbols   given by
%
\bea\label{Chris}
 \tilde \Gamma^\mu_{\alpha\beta} =
\Gamma^\mu_{\alpha\beta} + \frac{C'}{C}\, \delta^\mu_{(\alpha}\partial_{\beta)}\phi 
- \gamma^2 \, \frac{C'}{2C} \,\partial^\mu\phi \, g_{\alpha\beta}  
+\frac{D}{C}\,\gamma^2\,\partial^\mu\phi \left[\!\nabla_\alpha\nabla_\beta \phi + \!
 \left(\frac{D'}{2D} -\!\frac{C'}{C}\right) \!\partial_\alpha\phi \partial_\beta \phi \right]\,, \nonumber\\
\eea
%
and we have introduced the ``Lorentz factor" $\gamma$ defined as 
\be\label{Lorentz}
\gamma = \frac{1}{\sqrt{1+\frac{D}{C} (\partial\phi)^2}}\,.
\ee
%

In this frame,  the energy-momentum tensor is defined as 
%
\be\label{tmunuD}
\tilde T^{\mu\nu}  = \frac{2}{\sqrt{-\tilde g}} \frac{\delta S_{M}}{\delta \tilde g_{\mu\nu} } \,
\ee
and  the disformal energy-momentum tensor can be written as:
\be\label{pfluiD}
\tilde T^{\mu\nu}  = (\tilde \rho +\tilde P) \tilde u^{\mu}  \tilde u^{\mu} + \tilde P\, \tilde g^{\mu\nu}\,,
\ee
%
where   $\tilde u^{\mu} = C^{-1/2} \gamma \,u^{\mu} $. 
Using \eqref{tmunuD}, we  obtain a relation between the energy momentum tensor in both frames as:
\be\label{tmunuDD}
\tilde T^{\mu\nu} = C^{-3} \gamma \, T^{\mu\nu}\,.
\ee
%

Further using  \eqref{pfluiD} we  arrive at a relation among the energy densities and pressures in both frames, given by
\be\label{eqrhos}
\tilde \rho = C^{-2} \gamma^{-1} \rho \,, \qquad \tilde P = C^{-2}\gamma\, P,
\ee
and therefore the equations of state in both frames are related by $\tilde \omega = \omega\, \gamma^2$. Note that in the pure conformal case, $D=0$, $\gamma=1$ and therefore $\tilde \omega = \omega $.
%%%%%%%%%%%%%%%%

\section{Cosmological equations}

Consider an homogeneous and isotropic  FRW metric $g_{\mu\nu}$, 
\be
ds^2 = - dt^2 + a(t)^2 dx_i dx^i  \,,
\ee
where $a(t)$ is the scale factor.  In this background, the Einstein and  Klein-Gordon equations  become, respectively
\bea
&& H^2 =\frac{\kappa^2}{3} \left[\rho_\phi +\rho\right]\,, \label{friedmann1}\\
&& \dot H + H^2 = -\frac{\kappa^2}{6}\left[ \rho_\phi+ 3P_\phi +\rho +3 P \right]\,,\label{friedmann2}\\
&& \ddot \phi +3H\dot\phi +  V_{,\phi}+Q_0 =0 \label{kgSimple}\,,
\eea
%
where, $H= \frac{\dot a}{a}$, dots are derivatives with respect to $t$ and we have denoted $V_{,\phi} \equiv \frac{dV}{d\phi}$. Also the Lorentz factor becomes    $$\gamma= (1-D \,\dot\phi^2/C)^{-1/2}.$$
The continuity equations for the scalar field and matter are  given by
\bea 
&&\dot\rho_\phi + 3H(\rho_{\phi}+P_{\phi}) = -Q_0\dot\phi\,, \label{cont}\\
&&\dot\rho + 3H(\rho+P) = Q_0\,\dot\phi\,,\label{cont1}
\eea
where $Q_0$ is given by
%
%
\bea
Q_0 = \rho \left[ \frac{D}{C} \,\ddot \phi + \frac{D}{C} \,\dot \phi \left(\!3H + \frac{\dot \rho}{\rho} \right) \!+ \!\left(\!\frac{D_{,\phi}}{2C}-\frac{D}{C}\frac{C_{,\phi}}{C}\!\right) \dot\phi^2 +\frac{C_{,\phi}}{2\,C} (1-3\,\omega)
\right]. \nonumber \\
\eea
%
Using \eqref{cont1} we can rewrite this in a more compact  and useful form as 
\be\label{Q0}
Q_0 = \rho\left( \frac{\dot \gamma}{\dot \phi\, \gamma} + \frac{C_{,\phi}}{2C}  (1-3\,\omega \,\gamma^2) -3H\omega\,\frac{(\gamma^{2}-1) }{\dot \phi}\right)\,.
\ee

\ni Plugging this into  the  
(non-)conservation equation for dark matter \eqref{cont1},  gives:
\be\label{conservaDM}
\dot \rho + 3H (\rho + P\,\gamma^{2}) = \rho \left[\frac{\dot \gamma}{\gamma} + \frac{C_{,\phi} }{2C} \,\dot\phi\, (1-3\,\omega \gamma^2)\right]\,.
\ee
%
Using the relations for the physical proper time and the scale factors in the two frames, given by
\be\label{tildea}
\tilde a = C^{1/2} a   \,, \qquad \quad d\tilde \tau = C^{1/2} \gamma^{-1} d\tau\,,
\ee
we can define the disformal-frame Hubble parameter $\tilde H \equiv \frac{d \ln{\tilde a}}{d\tilde \tau}$, as
%
\be\label{tildeH}
\tilde H = \frac{\gamma}{C^{1/2}}\left[ H + \frac{C_{,\phi}}{2C}\dot \phi \right]\,,
\ee
so that \eqref{conservaD} takes  the standard form in terms of $\tilde H$:
\be\label{tilderho2}
\frac{d{\tilde \rho}}{d\tilde\tau}  +  3 \tilde H( \tilde \rho +\tilde P)  =0 \,.
\ee
%
Equations \eqref{tildeH} and \eqref{tilderho2} give the background evolution equations for the modified expansion rate and  matter's density evolution. 

\section{Master equations} \label{ME}

In order to solve the cosmological equations, it is convenient to replace time derivatives with derivatives with respect to the number of e-folds $N$, defined as $ N = \ln a/a_0$  and define  $\lambda  = \frac{V}{ \rho}(=\frac{\tilde V}{\tilde \rho})$.
With these definitions, we can rewrite the Friedmann equation \eqref{friedmann1} and $Q_0$ as:
\bea 
H^2 &=& \frac{\kappa^2 \rho}{3} \frac{(1+\lambda)}{\left(1- \frac{\kappa^2\phi'^2}{6 }\right)} \label{fried2}\,, \\
\!\!\frac{Q_0}{\rho}& =& \frac{\gamma^2 H^2}{2}\Bigg[ \frac{2D}{C}  \phi'' -\frac{2D}{C}  \phi' \!\left(
3\,\omega + \frac{\kappa^2\phi'^2}{2} + \frac{3(1+\omega) B}{2(1+\lambda)} \right) 
+ \left(\frac{D}{C}\right)_{\!\!,\phi} \!\!\phi'^2  + \frac{C_{\!,\phi}}{H^2 C}(\gamma^{-2}-3\omega) \Bigg], \nonumber \\
\eea
where  here we  denote $'= d/dN$. Note also that \eqref{fried2} implies that $\kappa\,\phi' \leq\pm \sqrt{6}$. 


Using these equations and further defining a dimensionless scalar field $\varphi = \kappa\,\phi$,  we can rewrite \eqref{friedmann2} and \eqref{kgSimple} as:
%
\bea
&& H' = - H \left[ \frac{3B}{2(1+\lambda)}   (1+ \omega) + \frac{\varphi'^2}{2}\right], \label{Hprime}\\ \nonumber \\
&&  \varphi'' \left[1\!+\! \frac{3H^2\gamma^2 B}{\kappa^2(1+\lambda)} \frac{D}{C} \right]  
+ 3\,\varphi'  \left[1- \omega \frac{3H^2\gamma^2 B}{\kappa^2(1+\lambda)} \frac{D}{C} \right] +
\frac{H'}{H} \varphi' \left[1+ \frac{3H^2\gamma^2 B}{\kappa^2(1+\lambda)} \frac{D}{C} \right]  
 \nonumber \\
&& \hskip1.5cm  + \frac{3B}{1+\lambda} \alpha(\varphi)(1-3 \,\omega \gamma^2) 
  + \frac{3B\lambda}{(1+\lambda)} \frac{V_{,\varphi}}{V}  + \frac{3 H^2 \gamma^2 B}{\kappa^2(1+\lambda)} \frac{D}{C}\left[ (\delta(\varphi) - \alpha(\varphi)) \,\varphi'^2 \right]   =0\,, 
  \nonumber \\
  \label{phiHeq}
\eea
where we  defined: 
\bea
&&B \equiv 1-\frac{\varphi'^2}{6},  \\
&& \gamma^{-2} = 1- \frac{H^2}{\kappa^2}\frac{D}{C} \varphi'^2\,,  \label{gammaH}\\
&& \alpha(\varphi) = \frac{d \ln C^{1/2}}{d\varphi}, \label{alphaeq}\\
&& \delta(\varphi) = \frac{d \ln D^{1/2}}{d\varphi} \label{deltaeq}\,.
\eea

One can solve the system of coupled equations above for $H$ and $\varphi$ as functions of $N$. However, in some cases it is simpler to use \eqref{Hprime} into \eqref{phiHeq}  and solve the following disformal master equation: 
%
\bea\label{mastereqsimple}
 \frac{2(1+\lambda)}{3 B} \,\varphi'' &+& \left(2\lambda +1 -\omega \right)\varphi'  + 2 \lambda\, \frac{d\ln  V}{d\varphi}   +2 (1 -3\,\omega\,\gamma^{2})\,\alpha(\varphi)
  \nonumber \\
&+ & \,  \frac{2\gamma^2(1+\lambda)}{3B}\frac{D \rho}{C} \left(\varphi'' -3\,\varphi'
\left[ \omega + \frac{\varphi'^2}{6} +  \frac{(1+\omega) B}{2(1+\lambda)}\right] 
+ \frac{C}{2D}\left(\frac{D}{C}\right)_{\!,\varphi} \varphi'^2 \right) \!=0 \,,  \nonumber \\
\eea
with $\gamma$ given by:
\be
\gamma^{-2} = 1- \frac{(1+\lambda)}{3B}\frac{D \rho}{C} \varphi'^2\,.
\ee

\ni From \eqref{mastereqsimple} we see that the conformal case is recovered for  $D=0$ when the second line vanishes. Moreover, the disformal piece appears always together with derivatives of the scalar field, as expected and also nontrivially coupled to the energy density. This complicates considerably the analysis of the disformal case, as we will see below.  

%
\section{Modified expansion rate}\label{Sec:Mod2}

The effect of the expansion rate during the early time evolution due to the  presence of a scalar field can be extracted from the Hubble parameter evolution in the  Jordan frame defined as:
$$
\tilde H = d (\log \tilde a )/d\tilde \tau, 
$$ 
which can be written using \eqref{tildea} as:
\be\label{Htilde}
\tilde H = \frac{H\gamma}{C^{1/2}} \left(1+ \alpha(\varphi) \varphi' \right)\,, 
\ee
 where remember that  $\gamma$ depends on $H$  (or $\rho$) as seen from \eqref{gammaH}, while in the pure conformal case  $D=0$ and $\gamma = 1$. Note that in principle, the factor $(1+ \alpha(\varphi) \varphi' )$ can be positive or negative, indicating an expansion or contraction modified rate. We stick to positive definite values for this factor and therefore only modified expansion rates, though, in principle, one could have a brief contraction period during the early universe evolution, before the onset of BBN\footnote{See  \cite{Contraction} for a review on scenarios with a possible contraction phase in the early universe.}. 
 % 
 Moreover, notice that while $\tilde H$ can grow during the cosmological evolution, the null energy condition (NEC) is not violated. This is because the Einstein frame expansion rate $H$ is dictated by the energy density $\rho$ and pressure $p$, which obey the NEC and therefore $\dot H<0$ during the whole evolution, as it should (see for example \cite{Creminelli}).

 
We further want to relate the modified expansion rate to the expected expansion rate in general relativity (GR), that is: 
\be\label{HGR}
H_{GR}^2 = \frac{\kappa_{GR}^2}{3} \, \tilde \rho\,.
\ee
%
%
We can do this be using the Friedmann equation \eqref{fried2} and the relation between the energy densities \eqref{eqrhos} to write
\be\label{Heq2}
\gamma^{-1} H^2 =  \frac{\kappa^2}{\kappa_{GR}^2} \frac{ C^2\,(1+\lambda)}{B} \, H_{GR}^2  \,.
 \ee
Using the definition of $\gamma$ (see \eqref{gammaH}) into this equation, one finds a cubic equation for $H^2$ in terms of all the other parameters. The real positive solution to that equation can then be replaced into  \eqref{Htilde} to find the modified expansion rate $\tilde H$, which will thus be a complicated function of $H_{GR}$ as we now see. 
%
The cubic equation for $H$ takes the form:
%
\be\label{eqforH2}
d_1 H^6-H^4 + d_2^2=0\,,
\ee
where 
\bea\label{Eqad}
d_1=\frac{D}{C}\frac{\varphi'^2}{\kappa^2} \,, \hspace{1cm} d_2=
\frac{\kappa^2}{\kappa_{GR}^2}\frac{C^2 (1+\lambda)H_{GR}^2}{B}\,.%
\eea
 The solutions  to \eqref{eqforH2} can be written as
%
\be\label{Hdis}
H^2=\frac{1}{3 d_1}\left(1+\left(\frac{2}{\Delta}\right)^{1/3}+\left(\frac{\Delta}{2}\right)^{1/3}\right)\,,
\ee

\ni with $\Delta=2-27d_1^2d_2^2+d_1d_2\sqrt{27(27d_1^2d_2^2-4)}$. The other two solutions can be  obtained by replacing 
$$ \left(\frac{2}{\Delta}\right)^{1/3} \to e^{2\pi i/3} \left(\frac{2}{\Delta}\right)^{1/3}
\qquad {\rm and } \qquad \left(\frac{2}{\Delta}\right)^{1/3} \to e^{4\pi i/3} \left(\frac{2}{\Delta}\right)^{1/3}\,.$$
%
We are interested in real positive solutions for $H^2$. One possibility to get this is to have the imaginary part of $(\Delta/2)^{1/3}$ vanish by requiring  that $\Delta>0$, which is impossible. 
Therefore, the way to obtain real solutions for $H$ is to have the imaginary parts of $(\Delta/2)^{1/3}$ and  $(\Delta/2)^{-1/3}$ cancel each other, leaving a real positive solution. 

For  this, we need that  $27d_1^2d_2^2 \leq 4$, which implies the following  relation between the conformal and disformal functions: 

% 
\be\label{ConditionD}
\frac{3\sqrt{3} \,D C\, \varphi'^2 (1+\lambda)}{B} \frac{H^2_{GR}}{\kappa_{GR}^2}\leq 2\,.
\ee
%
Under this condition, we can rewrite $\Delta$ as:
$$\Delta = 2-27d_1^2d_2^2+ i d_1d_2\sqrt{27(4- 27d_1^2d_2^2)}\,, $$
which allows us to define a complex number $Z\equiv \Delta/2$ and it is easy to check that $\bar Z = 2/\Delta$ and thus
 $|Z|^2 =1$. Denoting further $Z_i$  with $i=1,2,3$ denoting the three solutions to $H$ as explained above, the solutions for $H$, \eqref{Hdis} takes the simple form:
 \be\label{Hdis2}
H^2_i=\frac{1}{3 d_1}\left[1+Z_i^{1/3} + \bar Z_i^{1/3}\right] \,,
\ee
and remember that we are interested only in the real positive solution. We can now plug in \eqref{Hdis2}, as well as the expression for $\gamma$ in terms of $H$ into the Jordan frame expansion rate \eqref{Htilde}, can be written as:
%
\be\label{Htilde2}
\tilde H^2 = \frac{\kappa^2}{\kappa^2_{GR}} \frac{  \gamma^3\,  C(1+\lambda) (1+\alpha(\varphi)\varphi')^2 }{B} \,
\,H_{GR}^2\,, 
\ee
where there is a non-trivial dependence of $H_{GR}$ encoded in 
\be\label{gammai}
\gamma_i = \frac{1}{3d_1d_2}\left[1+Z_i^{1/3} + \bar Z_i^{1/3}\right] \,.
\ee
%
In the conformal case, $D=0$, $\gamma =1$ and therefore \eqref{Htilde2} is simply
%
\be\label{ratioH}
\tilde H^2 =  \frac{\kappa^2}{\kappa_{GR}^2}\frac{C (1+\lambda) (1+\alpha(\varphi) \varphi')^2}{B}\, H^2_{GR} \,.
\ee
From  this relation we define a speed-up parameter $\xi$, which will be useful below to measure the departures from the $GR$ expansion rate result:
\be\label{xieq}
\xi \equiv \frac{\tilde H}{H_{GR}}\,. 
\ee

We now revisit the conformal case, discussed originally in \cite{Catena}. We first solve (numerically) the master equation for the scalar field \eqref{mastereqsimple} in order to compute the modified expansion rate $\tilde H$ and compare it with the standard expansion rate, $H_{GR}$. 

Before solving the master equation \eqref{mastereqsimple}, we would like to write it in terms of Jordan  frame quantities $\tilde \omega = \omega \gamma^2$, $\tilde \rho = C^{-2}\gamma^{-1} \rho$. Moreover, the number of e-folds $N$ can be expressed in terms of Jordan frame quantities as follows. In this frame, the entropy is conserved and is given by $\tilde S=\tilde a\, \tilde s$, where $\tilde s= \frac{2\pi}{45} g_s(\tilde T) \tilde T^3$. So, the conservation of entropy and \eqref{tildea} show that $N$ is a function of temperature and the scalar field as:
\be\label{N}
N\equiv\ln\frac{a}{a_0}=\ln\left[\frac{\tilde T_0}{\tilde T}\left(\frac{g_s(\tilde T_0)}{g_s(\tilde T)}\right)^{1/3}\right]+\ln\left[\frac{C_0}{C}\right]^{1/2}\,.
\ee
Therefore, we can introduce the  parameter, $\tilde N$, defined as
\be\label{Ntilde}
\tilde N \equiv \ln\left[\frac{\tilde T_0}{\tilde T}\left(\frac{g_s(\tilde T_0)}{g_s(\tilde T)}\right)^{1/3}\right]\,,
\ee
and transform to  derivatives w.r.t.~$\tilde N$  (assuming well behaved functions): 
\bea\label{NtotildeN}
%&&
\varphi' =  \frac{1}{\left(1-\alpha(\varphi) \frac{d\varphi}{d\tilde N}\right)} \frac{d\varphi}{d\tilde N} \,,\qquad %\\
%&& 
\varphi'' = \frac{1}{\left(1-\alpha(\varphi) \frac{d\varphi}{d\tilde N}\right)^3} \left( \frac{d^2\varphi}{d\tilde N^2} + \frac{d\alpha}{d\varphi}  \left(\frac{d\varphi}{d\tilde N}\right)^3\right)\,. \\ \nonumber 
\eea
%
In a slight abuse of notation and to keep expressions  neat, in what follows we  denote derivatives w.r.t.~$\tilde N$ with a prime $'$. 


\section{Pure conformal case}\label{Sec:PureC}

We start with the pure conformal case. That is, we take  $D(\phi) =0$ in \eqref{mastereqsimple} and therefore $\gamma=1$ (and  $\tilde \omega = \omega$). Moreover, during the radiation and matter dominated eras, of interest for us,  the potential energy of the scalar field is subdominant and thus, we take  $\lambda \sim 0$. 
Therefore the  master equation \eqref{mastereqsimple} simplifies to:
\bea\label{masterConf}
 \frac{2}{3 (1-\varphi'^2/6)} \,\varphi'' &+& \left(1 -\tilde \omega \right)\varphi'  + 2 (1 -3\,\tilde \omega)\,\alpha(\varphi) =0, 
\eea
%
which in terms of derivatives wrt $\tilde N$ takes the form:
%
\bea\label{masterConf2}
\frac{1 }{3B\left[1-\alpha(\varphi) \varphi' \right]^3} \left(\varphi'' + \frac{d\alpha}{d\varphi}  \left(\varphi'\right)^3\right)
+ \frac{(1-\tilde \omega)  }{\left[1-\alpha(\varphi)\varphi'\right]} \,\varphi'
+(1-3\,\tilde \omega) \, \alpha(\varphi)=0 \,, \nonumber \\
\eea
%
%
where  $B= 1-\frac{\left(\varphi'\right)^2 }{6\left(1-\alpha(\varphi) \varphi'\right)^2} $. 
Using the relation between $\tilde H$ and $H_{GR}$ defined in \eqref{ratioH}, we can write the speed-up parameter as 
%
\be
\xi = \frac{\tilde H}{H_{GR}} = \frac{C^{1/2}(\varphi)}{C^{1/2}(\varphi_0)} \frac{1}{\left(1-\alpha(\varphi) \varphi'\right)\sqrt{B}} \frac{1}{\sqrt{1+ \alpha^2(\varphi_0)}}\,,
\ee
%
where we have used the relation between the bare gravitational constant and that measured by local experiments for  conformally coupled theories \cite{Nordtvedt:1970uv}:  
\be\label{ratioG}
\kappa_{GR}^2 = \kappa^2 C(\varphi_0) [1+\alpha^2(\varphi_0)]\,,
\ee
where $\varphi_0$ is the value of the scalar field at present time. 


\subsection{Expansion rate modification}\label{ERM}

%

The scalar equation in the conformal case \eqref{masterConf}, as function of $N$ (for  $\lambda =0$)  contains a term which can be interpreted as an effective potential, dictated by $V_{eff} = \ln C^{1/2}$. For a strictly radiation dominated era, $\tilde \omega =1/3$, the effective potential term vanishes and we are left with an equation that can be solved analytically \cite{Coc}, giving  $\varphi' \propto e^{-N}$. That is, any initial velocity will rapidly go to zero (remember that from the Friedmann equation \eqref{fried2}, $\varphi'$ is constrained to be $|\varphi'| \lesssim  \sqrt{6}$).  
%
Therefore we explore the effects of having a non-zero initial velocity in our analysis below. 
Since the scalar field is expressed in Planck units, we focus on order one or smaller field variations $\Delta\varphi$. One can check, using the analytic solution to \eqref{masterConf} deep in the radiation era, that for initial velocities $\varphi'_0 \ll \pm \sqrt{6}$, the total field displacement is of order $\Delta \varphi \sim \varphi'_0$ \cite{Coc}.  However, given that we don't know much about the theory before BBN, we  explore different initial values for $(\varphi_0, \varphi'_0)$ and study their 
%
consequences. In particular we explore initial values $\varphi_0$ and $\varphi'_0 \in (-1.0, 1.0)$.
 
We now concentrate on an explicit conformal factor. We  use  the same conformal factor as that studied in \cite{Catena}, which is given by:
\be\label{Cfactor}
C(\varphi) = (1+ b \,e^{-\beta\, \varphi})^2\,,
\ee
with the values $b=0.1$, $\beta=8$, which have been shown to satisfy the constraints imposed by tests of gravity, for the parameters $\alpha, \beta$, $\xi$. As we will see, the requirement of reaching the  GR expansion rate value by the time of the onset of BBN drives these parameters to very small values, which are thus consistent with the constraints from gravity for their values today. 


As we discussed above, in the equation of motion for $\varphi$, with $\tilde\omega\neq 1/3 $, the conformal factor acts as an effective potential on which the scalar field moves, damped by the Hubble friction (see \eqref{kgSimple}). 
Since any initial velocity $\varphi'$ goes rapidly to zero deep in the radiation era, in the subsequent evolution of $\varphi$,  the term $(\varphi' )^2$  in the master equation will be negligible. In this regime, the equation is that of a particle moving in an effective potential with a damping term. Therefore, one can understand the evolution of $\varphi$  from the form of the effective potential ($V_{eff}=\ln C^{1/2}$) and the initial conditions used. For the conformal function we are considering \eqref{Cfactor}, one sees that there is a set of initial conditions that will give rise to an interesting behavior in $\varphi$, and therefore an interesting modified expansion rate in the Jordan frame $\tilde H$ \eqref{ratioH}, as we now explain.


In general, both the initial position and velocity of the scalar field can take any value, positive and negative. In the runaway effective potential dictated by $\ln C^{1/2}$ for the conformal factor \eqref{Cfactor} we consider here, we have the following possibilities. i) The scalar field starting somewhere up in the runaway effective potential with zero initial velocity. In this case, the scalar field will roll-down the potential, eventually stopping due to Hubble friction, at some constant value of $\varphi$, which depends on its initial value. So long as $(1+\alpha(\varphi) \,\varphi')$  stays positive (see \eqref{Htilde}),  $C$  will evolve rapidly towards $1$.
 More generally, the initial velocity can be different from zero. If the initial velocity is positive, the behavior will be similar to the previous case. The field will roll-down the effective potential towards its final terminal value. 
 

ii) A more interesting possibility arises when one allows for negative initial velocities. In this case, the field will start rolling-up the effective potential towards smaller values of the field, eventually turning back down and moving towards its terminal value. 
%
It is easy to see that an interesting effect happens when the field starts at an initial positive value. Given a sufficient initial negative velocity the field will move towards negative values until its velocity becomes zero and then positive again, as it rolls back down the effective potential. This change in sing for the scalar evolution will produce a pick in the conformal function that will give rise to a non-trivial modification of Jordan's frame expansion rate $\tilde H$, as we are looking for. 
As mentioned before, we are interested in (sub-)Planckian initial values $(\varphi_0, \varphi'_0)$, such that $\tilde H>0$. With these requirements, one can see that given an initial negative velocity, there is a suitable initial value of the scalar field such that the behavior just described holds and the expansion rate $\tilde H$ has an interesting evolution before the onset of BBN. At late times the conformal function goes to one and the GR expansion is recovered. 
 We show this behaviour explicitly in Figures \ref{FieldFig} and \ref{CFig} where we plot the numerical solution for the evolution of $\varphi$ and  $C(\varphi)$  as functions of the temperature. In these plots we find $\varphi=\varphi(\tilde T)$ by first solving \eqref{masterConf2} numerically with initial conditions $ (\varphi_0,\varphi'_0) = (0.2,  -0.99)$  and then use \eqref{Ntilde} to express $\varphi(\tilde N)$ as function of $\tilde T$.
 %
 As we can see,  the conformal factor  starts growing  towards a maximum value as $\varphi$ moves to negative values, to rapidly drop down towards its GR value at $C\to 1$ as $\varphi$ moves down the effective potential towards positive values. This non-trivial effect will give rise to the possibility of re-annihilation, as we discuss below.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
\centerline{
\includegraphics[width=.65\textwidth]{Plots_phi_initial_02V2.pdf}}
%
\caption{Typical evolution of the scalar field  as temperature decreases. The initial values are $(\varphi,d\,\varphi/d\,\tilde N) = (0.2,-0.994)$.}
%\end{center}
\label{FieldFig}
\end{figure}


\begin{figure}[h!]
\centerline{
\includegraphics[width=.65\textwidth]{CevolutionV2Log.pdf}}
%\vspace*{-1cm}
\caption{Behaviour of the conformal factor, $C(\varphi)$ as a function of the temperature for the same initial values as in Fig.~\ref{FieldFig}.}
\label{CFig}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Based on the discussion above, we have solved the master equation \eqref{masterConf2}, to find the scalar field as a function of $\tilde N$ for various initial conditions,
where we see the interesting behavior explained above. The resulting modified expansion rate and its comparison with the standard case are shown in Figure \ref{HubblesFig} for the same initial conditions as in Figures \ref{FieldFig} and \ref{CFig}. In our numerical exploration, we choose initial conditions for which the notch in the expansion rate
(see Fig.~\ref{HubblesFig}) occurs closer to the BBN time. This has interesting consequences for the dark matter annihilation, as we discuss in chapter \ref{Ch:DM}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
%\vspace*{-5mm}%\hspace{-3mm}
\centerline{
\includegraphics[width=.8\textwidth]{HubbleRatesV2.pdf}}
%\vspace*{-1cm}
\caption{Comparing the Hubble expansion rate $\tilde H$ in the Jordan Frame with the standard Hubble expansion rate $H_{GR}$. 
The presence of the scalar field enhances and decreases the expansion rate during the radiation dominated era. This plot corresponds to
initial conditions given by $(\varphi_0,\varphi'_0) = (0.2,-0.994)$.}
\label{HubblesFig}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{The equation of state parameter, $\tilde \omega$}\label{Sec:EqState}

%
When solving the master equation \eqref{masterConf2}, we have taken into account an important effect that occurs during the radiation dominated era. Deep in this epoch, the equation of state is given by $\tilde\omega=1/3$. When a particle species in the cosmic soup becomes non-relativistic, $\tilde\omega$ differs slightly from $1/3$. 
When the temperature of the universe drops below the rest mass of each of the particle types, there are non-zero contributions to $1-3\tilde \omega$. This activates the 
effective potential, which can be seen in the last term of \eqref{masterConf2}, and displaces, or ``kicks`` the field along $V_{eff}$.

To examine this effect in more detail, we start by writing $1-3\,\tilde \omega$ during the early stages of the universe as in \cite{Catena} and \cite{MW}
%
\be\label{eqomega}
1-3\,\tilde\omega = \frac{\tilde \rho - 3\, \tilde p}{\tilde \rho} =\sum_{A} \frac{\tilde \rho_A - 3\tilde p_A}{\tilde \rho} +\frac{\tilde \rho_m}{\tilde \rho}\,,
\ee
where the sum runs over all particles in thermal equilibrium during the radiation dominated era and $\tilde \rho_m$ is the contribution
from the non-relativistic decoupled and pressureless matter. The summation over all the particle is responsible for the kicking effect discussed above. Then, a kick 
function is defined as
\be\label{kick}
\Sigma(\tilde T) \equiv \sum_{A} \frac{\tilde \rho_A - 3\tilde p_A}{\tilde \rho}\,,
\ee
where the energy density $\tilde\rho_A$ and pressure $\tilde p_A$ of each type $A$ of particle are given by
\be\label{erhoA}
\tilde \rho_A(\tilde T) =\frac{g_A}{2\pi^2}\int^\infty_{m_A}\frac{\left(E^2-m_A^2\right)^{1/2}}{\exp(E/\tilde T)\pm 1}E^2dE\,,
\ee
\be\label{epA}
\tilde p_A(\tilde T) =\frac{g_A}{6\pi^2}\int^\infty_{m_A}\frac{\left(E^2-m_A^2\right)^{3/2}}{\exp(E/\tilde T)\pm 1}E^2dE\,,
\ee
with $g_A$ being the number of internal degrees of freedom of species of type $A$ and the plus (minus) sign in the integral corresponds to fermions (bosons).

After using $\tilde \rho \simeq \pi^2 g_{eff}(\tilde T) \tilde T^4/30$, \eqref{erhoA} and \eqref{epA} $\Sigma$ becomes
\be\label{kick1}
\Sigma(\tilde T)  = \sum_{A} \frac{15}{\pi^4} \frac{g_A}{g_{eff}(\tilde T)} \,y_A^2 \int_{y_A}^{\infty} dx \frac{\sqrt{x^2-y^2_A}}{e^x\pm 1}\,,
\ee 
where $g_{eff}(\tilde T)$ is the total number of relativistic degrees of freedom\footnote{To calculate $g_{eff}$ we follow the numerical procedure described in Appendix A of \cite{Erickcek:2013dea}.} and $y_A=m_A/\tilde T$. 

To compute \eqref{kick1}, we consider the Standard Model particle spectrum. In particular, we take into account the top quark, the Higgs boson, Z boson, W bosons, 
bottom quark, tau lepton, charm quark, charged pions, neutral pion, muon lepton and the electron (See table \ref{spectrum} for details on the masses and internal degrees of freedom per species). As we show in Figure \ref{sigmafig}, $\Sigma$ is mostly zero, except when the kicks happen. Then, during the radiation dominated era, we compute the equation of state parameter from  \eqref{eqomega}   as  $\tilde\omega=(1-\Sigma(\tilde T))/3$\footnote{$\tilde \rho$ is dominated by radiation, thus $\tilde \rho_m/\tilde \rho$ is negligible.}.


\begin{table}
\begin{center}
 \begin{tabular}{|lll|}
 \hline
  Particle &Mass (GeV)&$g_A$\\ \hline
    $t$ &173.2 &12\\ 
    $H$ &125.1 &1\\
    $Z$ &91.19 &3\\
    $W^\pm$ &80.39 &6\\
    $b$ &4.18 &12\\
    $\tau$ &1.78 &4\\
    $c$ &1.27 &12\\
    $\pi^0$ &0.140 &12\\
    $\pi^\pm$ &0.135 &2\\
    $\mu$ &0.106 &4\\
    $e$ &0.000511 &4\\ \hline
 \end{tabular}
\end{center}
\caption{Spectrum of particles used to calculate the kick function \eqref{kick1}. For each particle, we show its mass and number of internal degrees of freedom, $g_A$.}
\label{spectrum}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h!]
\centerline{
\includegraphics[width=.8\textwidth]{plotsigma.pdf}}
\caption{Thermal evolution of the kick function during the radiation dominated era. Outside the interval of temperatures shown, $\Sigma$ vanishes.}
\label{sigmafig}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In Figure \ref{plotomegaradiation} we show the evolution of $\tilde \omega$ between 10 TeV and 10 eV. This figure shows four troughs, which are the ``kicks'' mentioned above. 
Each kick corresponds to the transition of one or more particles to the non-relativistic regime. For example, the trough at around 0.5 MeV is due to the electron, while the one
at around 100 GeV is due to the heavy particles ($t$, $H$, $Z$ and $W$).

Towards the end of the radiation era, approaching the transition to the matter dominated era, \eqref{eqomega} takes the approximate form:
\be\label{lateomega}
 1-3\,\tilde\omega \simeq \frac{\tilde \rho_{m}}{\tilde \rho_m + \tilde \rho_r} \simeq \frac{1}{1+ \tilde T/\tilde T_{eq}} \,,
\ee
%
where $\tilde T_{eq} \sim {\cal O}(10^{-9})$GeV is the temperature at matter-radiation equality, that is, $\tilde \rho_m(\tilde T_{eq}) = \tilde \rho_r(\tilde T_{eq})$. 
We can now combine \eqref{kick} and  \eqref{lateomega} to compute the thermal evolution of \eqref{eqomega} in the radiation dominated and matter dominated
eras and use it in the master equation.
%

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h!]
\centerline{
\includegraphics[width=.65\textwidth]{plotomegaradiationV2.pdf}}
\caption{Evolution of $\tilde\omega$ in \eqref{eqomega} as function of temperature during the radiation dominated era.}
\label{plotomegaradiation}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
\subsection{Parameter Constraints}\label{Sec:ParCons2}


In scalar-tensor theories of gravity, there are some constraints on the parameters that need to be taken into account. 
Deviations from GR can be parametrised in terms of the post-Newtonian parameters $\gamma_{PN}$ and $\beta_{PN}$, which are given in terms of  $\alpha(\varphi_0)$ defined in \eqref{alphaeq} and  its derivative $\alpha'_0= d\alpha/ d\varphi |_{\varphi_0}$  as  \cite{Esposito,Will}:
\be \label{PostNewt}
\gamma_{PN} -1 = -\frac{2\alpha_0^2}{1+\alpha_0^2} \,, \qquad \beta_{PN} -1 = \frac12\frac{\alpha'_0\alpha_0^2}{(1+\alpha_0^2)^2} \,.
\ee
%
Solar system  tests of gravity, including the perihelion shift of Mercury, Lunar Laser Ranging experiments, and the measurements of the Shapiro time delay   by the Cassini spacecraft \cite{Lunar,Shapiro,Cassini} indicate that $\alpha_0$ should be very small, with values $\alpha_0^2 \lesssim 10^{-5}$, while binary pulsar observations impose that $\alpha'_0\gtrsim -4.5$.
%
The last constraint applies to the speed-up factor $\xi$, which has to be of order 1 before the onset of BBN. In our  examples we have 
$\alpha_0^2 \simeq 2\times10^{-5}$, $\alpha'_0 > 0$ and $\xi\approx1.05$.

\section{Conformal case plus a simple disformal coupling }\label{Sec:CoD}

%
We now discuss briefly the effect of the disformal factor in the metric \eqref{gtilde} to the expansion rate of the universe, $\tilde H$, and compare it to the conformal modification to $\tilde H$\footnote{We leave a detailed exploration for a future publication.}. Hence, we explore $D(\phi)\neq0$ for the same
conformal factor studied before, that is, $C(\varphi) = (1+ b \,e^{-\beta\, \varphi})^2$ for $b=0.1$, $\beta =8$. To investigate these modifications, we first need to look at the scalar field evolution with temperature.

%
In the pure conformal case studied above, we found the thermal evolution of the scalar field by solving the master equation \eqref{masterConf2} numerically, which is  \eqref{mastereqsimple} for 
$D(\phi)=0$.  However, to study the effects of the disformal factor on the scalar field, 
it is more convenient to solve the system of two coupled equations \eqref{Hprime} and
\eqref{phiHeq}. Using these equations we find solutions for the dimensionless scalar field $\varphi$, and for the expansion rate in the Einstein frame $H$.

Notice that solving the system of coupled equations or solving the master equation to find the thermal evolution of the scalar field are equivalent methods (as we have explicitly checked), because \eqref{mastereqsimple} it is nothing but a 
combination \eqref{Hprime} and \eqref{phiHeq}. However, while in the pure conformal case the master equation can be made independent of $H$ (or $\rho$), this is not the case for the more general disformal case, as we can see in \eqref{mastereqsimple}.


In the same way, as for the conformal case, we are interested mainly in the radiation and matter eras and therefore we can neglect the potential energy of the scalar field. Thus, we consider $V\sim 0$ and $\lambda=0$. Also, while solving the coupled equations we have to express $\omega$ in the Jordan frame by using $\tilde \omega = \omega \gamma^2$ and transform all derivatives w.r.t.
$N$ to derivatives w.r.t. $\tilde N$ by using \eqref{NtotildeN}.


With this information, we solve the system of coupled equations numerically to find the dimensionless scalar field $\varphi$ and the Hubble parameter $H$, as functions of the number of e-folds  $\tilde N$ (and the temperature). 
We choose the same initial conditions for the scalar field and its derivative as in the conformal case and to obtain the initial condition for $H$, we use \eqref{Hdis2}.

Once we have the solutions for $\varphi$ and $H$ as functions of temperature,  we can go back to \eqref{Htilde} and \eqref{gammaH}  to obtain the expansion rate for the disformal model. As an example, in Figure \ref{Hdisformal} we show the effects of a disformal factor given by  $D(\varphi) = D_0\,\varphi^2$ with $D_0 =-4.9\times 10^{-14}$.
In this plot, we illustrate the effect of the disformal contribution on the expansion rate  ($\tilde H_{Disformal}$) and compare it to the modified expansion rate for the conformal case ($\tilde H_{Conformal}$) and the standard case ($H_{GR}$).
 We use the same initial conditions as in Figures \ref{FieldFig} and \ref{CFig} for the scalar field and its derivative.

 Also, it is important to mention that for the case shown the parameter constraints described in section \ref{Sec:ParCons2} are satisfied. In particular we find  $\alpha_0^2 \simeq 2\times10^{-5}$, $\alpha'_0 > 0$ and $\xi\approx1.02$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h!]
\centerline{
\includegraphics[width=.65\textwidth]{Hdisformal.pdf}}
\caption{Comparing the modified expansion rate of the universe in the disformal and conformal scenarios for the same initial conditions as in Fig.~\ref{HubblesFig}.}
\label{Hdisformal}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
From our example, with $C$ and $D$ as indicated above, we can clearly see the differences from the disformally modified expansion rate $\tilde H_{Disformal}$  compared to the conformally modified and standard case, $H_{GR}$. The evolution of 
$\tilde H_{Disformal}$ is similar to that of  $\tilde H_{Conformal}$, having an enhancement and a decrement compared to the standard expansion rate $H_{GR}$.
Moreover, the main differences with respect to the conformal modification are the position of the notch and its shape. The notch is moved to higher temperatures and it becomes a little bit sharper.

%
These differences between the expansion rates can be understood from \eqref{phiHeq}. First, in this equation, we see that the factor
$\frac{3H^2\gamma^2BD}{\kappa^2(1+\lambda)C}$ in the coefficients of $\varphi''$, $\varphi'$ and $\varphi'^2$ vanishes when $D=0$. For the disformal example shown in Figure \ref{Hdisformal},  this factor is a very small correction to the equation, which is reflected in the slight shape modification
of $\tilde H_{Disformal}$ compared to $\tilde H_{Conformal}$. Second, the term proportional to $\delta(\varphi)$ plays a more important role, being responsible for the shifting of the notch. 

