\documentclass[a4paper,11pt]{article}
\usepackage[margin=1.27cm]{geometry}
\usepackage[dvipsnames]{xcolor}
\usepackage{subfiles,amsmath,amssymb,mathtools,bm,multirow,siunitx,graphicx,subcaption,tabularx,longtable,booktabs,hyperref,diagbox,colortbl,makecell,float,bbm,slashed,macros,soul,cancel}
\usepackage[most]{tcolorbox}
\newlength\tindent
\setlength{\tindent}{\parindent}
\setlength{\parindent}{0pt}
\renewcommand{\indent}{\hspace*{\tindent}}
\urlstyle{same}
\hypersetup{colorlinks=true, linkcolor=blue, filecolor=magenta, urlcolor=red, citecolor=cyan}

\numberwithin{figure}{section}
\numberwithin{equation}{section}
\numberwithin{table}{section}
\usepackage{titling}

\graphicspath{{./Figures/}}

\title{\textbf{Modified gravitational waves}}
\author{}
\date{}

\begin{document}

\maketitle

\begin{abstract}



\end{abstract}

\section{Introduction}

Lovelock's theorem states that
the only second-order, 
local gravitational field equations
that can be derived from
an action containing only the metric tensor
in four dimensions
are the Einstein field equations (EFEs)
~\cite{Lovelock71,Lovelock72}.
Therefore, breaking any of the conditions stated in the theorem modifies gravity,
i.e.,
adding additional fields other than the four-dimensional metric,
having theories in a dimension other than four,
including higher-than-second-order derivatives,
or introducing non-locality.
See Ref.~\cite{Clifton+11} for a comprehensive review on modified gravity in general.

Among these approaches,
Horndeski gravity is a generic class of four-dimensional models that includes a scalar field 
and contains only second-order derivatives in the field equations~\cite{Horndeski74}.
It has been shown that in Horndeski gravity,
the linear perturbations on an FLRW background
can be characterized by four free parameters
separable from the full theory~\cite{Bellini+14}.
These free parameters are
\begin{itemize}
\item Kineticity $\alpK$, 
indicating the kinetic energy of scalar perturbations.
\item Braiding $\alpB$,
indicating the mixing of the scalar and tensor kineticities.
\item Planck-mass run rate $\alpM$,
denoting the time-dependence of the effective Planck mass.
\item Tensor speed excess $\alpT$,
measuring the deviation of the speed of gravity from the speed of light.
\end{itemize}
In terms of GWs,
only $\alpM$ and $\alpT$ enter.
Moreover, although the graviton is massless in Horndeski gravity,
a massive graviton can be included in certain beyond-Horndeski gravity models.
\red{Alberto: We should include some motivation to consider 
such types of theories, for example to explain dark energy, and 
also current constraints on these modifications. In particular, 
to describe the constraints in the prerecombination epoch, using CMB observations, for example, and as a way to alleviate
the Hubble tension.}

\section{GW equation in Horndeski gravity}
\label{GW_eq}
%YH: should we say Horndeski here? the modified equation can be from other theories too, right?
%AR: here we are 'only' considering Horndeski theory since we are focusing on the effect of the Horndeski parameters alpT and alpM
%YH: okay. then should we still have \mg term in the equation?
%YH: it was here as I thought it might be good to compare the modifications from \mg and \alpM -- both change slopes in IR
%YH: but we could also get rid of it to focus on \alpM and \alpT only

The GW equation describes the tensor-mode perturbations $\hij$ over the isotropic and
homogeneous expanding Universe, described by the Friedmann-Lema\^itre-Robertson-Walker
(FLRW) metric, such that the line element is
\begin{equation}
\dd s^2 = \dd t^2 + a ^2 \bigl[\delta_{ij} + \hij (\xx, t)\bigr] \dd \xx^2.
\end{equation}
In general relativity, the GW equation in Fourier space reads
\begin{equation}
\ddot{\tilde h}_\ij (\kk, t) + 3H\dot {\tilde h}_\ij (\kk, t) + k^2 \tilde h_\ij (\kk, t) =
2\kap \tilde T_\ij^{\rm TT} (\kk, t),
\end{equation}
where $H\equiv\dot a/a$ is the Hubble parameter,
$\kap = 8\pi\GN$ is the gravitational coupling constant, and a dot indicates
time derivatives.
The source term corresponds to the traceless-transverse (TT) projection of the stress-energy tensor
$\tilde T^{\rm TT}_{ij} = \Lambda_{ijlm} \tilde T_{lm}$.
In Horndeski gravity~\cite{Horndeski74},
which is the most general Lorentz-invariant scalar-tensor gravity that involves up to second derivatives in the Lagrangian,
the modified parameters $\alpM$ and $\cT^2 = 1 + \alpT$ appear in the GW equation.
Additionally, including the graviton mass $\mg$, the GW equation reads~\cite{Nunes+18}
\begin{equation}
\ddot{\tilde h}_\ij (\kk, t) + (3 + \alpM)H\dot{\tilde h}_\ij (\kk, t) +
(\cT^2 k^2 + \mg^2)\tilde h_\ij (\kk, t) = 2\kap \tilde T_\ij^{\rm TT} (\kk, t),
\label{eqn:GW_mod}
\end{equation}
where 
$\cT$ denotes the GW propagation speed,
and $\alpM$ is the run rate of the effective Planck mass $M_\eff^2$, parameterized as
\begin{equation}
\alpM = \frac{\dd\ln{M_\eff^2}}{\dd\ln{a}}. 
\end{equation}
We introduce conformal time $\eta$, such that $\dd t = a \dd \eta$ and, hence,
we transform the time derivatives to derivatives with respect to conformal time (indicated by a
prime)
\begin{equation}
a \dot h_\ij = h_\ij',\quad
a^2 \ddot h_\ij = h_\ij'' -
\HH h_\ij',\quad
\dot a = \HH,\quad
\ddot a = \frac{a''}{a^2} - \frac{\HH^2}{a},
\end{equation}
where $\hhh\equiv a'/a$ is the conformal Hubble rate.
We can rewrite \Eq{eqn:GW_mod} as
\begin{equation}
\tilde h_\ij'' (\kk, \eta) + (2 + \alpM)\hhh \tilde h_\ij' (\kk, \eta) +
a^2 \Big(\cT^2 k^2 + \mg^2\Big)\tilde h_\ij (\kk, t) = 2\kap a^2 \tilde T_\ij^{\rm TT} (\kk, t).
\label{eqn:GW_conf}
\end{equation}

Next, adopting comoving strains, energy-momentum tensor,
wave numbers, and graviton mass, such that $\bar h_\ij = ah_\ij$,
$\bar T_\ij = a^4 T_\ij$, $\bar k = a k$, and 
$\bar m_{\rm g} = a m_{\rm g}$ then
\begin{equation}
%T_\ij = \frac{\bar T_\ij}{a^4},\quad
%h_\ij = \frac{\bar h_\ij}{a},\quad
a h_\ij' = \bar h_\ij' - \HH \bar h_\ij,\quad
a h_\ij'' = \bar h_\ij'' - 2 \HH \bar h_\ij' +
\biggl(2 \HH^2 - \frac{a''}{a} \biggr) \bar h_{ij},
\end{equation}
and \Eq{eqn:GW_conf} becomes
\begin{equation}
\bar h_\ij'' + 
\alpM\hhh\bar h_\ij' + 
\biggl(\cT^2\bar k^2 + \bar m_{\rm g}^2 - \alpM\hhh^2 - \frac{a''}{a}\biggr)\bar h_\ij
= 2\kap\frac{\bar T_\ij}{a}.
\label{eqn:GW_conf_comove}
\end{equation}
%
% \begin{equation}
%     \frac{a''}{a^3} a^2 \bar h_{ij} \propto T_{ij} a^3
%     \Rightarrow \frac{a''}{a^3} \bar h_{ij}
%     \propto a
% \end{equation}
Note that \Eq{eqn:GW_mod} does not depend
explicitly on $a$.
On the other hand, due to the inclusion of conformal
time $\eta$, scaled strains, and comoving variables,
the terms in
\Eq{eqn:GW_conf_comove} scale with the value of the
scale factor as $a^3$, which is only physically defined as
a ratio of its value at different epochs.
However, the final solution $h_{ij}$ is independent
of the reference value at which we set $a = 1$
by convention, although
the scaled strains $\bar h_{ij}$ are inversely proportional
to the reference value of $a$.

Finally,
we normalize the comoving stress tensor against the critical energy density $\eee_\crit^*$, and
the time derivatives, Hubble rate, wave numbers,
and graviton mass against the Hubble rate
at the initial time of our simulations, which belongs to the
radiation era,
i.e.,
\begin{equation}
\bar{\bar T}_\ij = \bar T_\ij/\eee_\crit^* = \kappa
\bar T_\ij/3 H_*^2 ,\quad
%\p_{\bar\eta} = \p_{(\hhh_*\eta)} = \hhh_*^{-1}\p_{\eta}
%\Map
\p_{\eta} = \HH_*\p_{\bar\eta}%
%= (a_* H_*) \, \p_{\bar\eta}
,\quad
%\bar\hhh = \hhh/\hhh_*,\quad
\bar\hhh = \hhh/\HH_*,\quad
%\bar k = k/\hhh_*,\quad
\bar{\bar k} = \bar k/\HH_*,\quad
%\bar\mg = \mg/\hhh_*,
\bar{\bar{m}}_{\rm g} = \bar{m}_{\rm g}/\HH_*.
\end{equation}
Dropping the overbars,
then \Eq{eqn:GW_conf_comove} yields the normalized
GW equation,
\begin{align}
h_\ij'' + 
\alpM\hhh h_\ij' + 
\biggl(\cT^2 k^2 + \mg^2 - \alpM\hhh^2 - \frac{a''}{a}\biggr) h_\ij
& = \frac{6}{a} a_*^2 \, T_\ij,
\label{modGR}
\end{align}
where the prime denotes time derivative with respect to the
normalized conformal time $\bar \eta = \HH_* \eta$.
This equation is proportional to $a$ via $h_\ij$ on the
left hand side and via $a_*^2/a$ on the right hand side.
Hence, it depends
on the normalization of the scale factor $a$.
Finally, we set $a_* = 1$ for simplicity.
Hence, after solving \Eq{modGR}, we recover the physical
strains at a scale $a$
by multiplying $a_*/a$ to the scaled strains $h_\ij$.

%AR: please feel free to delete the commented text if you agree with the corresponding modifications in this section.
%YH: the GW equation or Friedmann section? both have commented text

\section{Friedmann equations}

In previous works using the {\sc Pencil Code}, we have 
restricted the simulations to specific eras (e.g., 
radiation-dominated or matter-dominated) using
an equation of state $p = w \rho$ with constant $w$
and ignoring the evolution of degrees of freedom (DOFs) within
the simulation.
In GR, the previous formulation is justified since the evolution
of $h_\ij$ when the source is inactive ($T_\ij \sim 0$) can
be approximated to dilute as $h_\ij \propto a^{-1}$, if one
neglects evolving relativistic DOFs and
transitions between radiation, matter, and dark energy 
dominations \cite{Saikawa:2018rcs}.
In the present work, we focus in including such effects
and solving the GW equation
in modified gravity, which present richer dynamics even when
the source is inactive, from the time of their generation
up to present time.

\red{Alberto: made some changes and moved this section above,
after GW equation. Added some introductory text above.}

This leads to discontinuities in the time evolution of
$a$ as a function of $\eta$ at the times where the value
of $w$ (defined as $p = w \rho$) changes discontinuously.

\red{Alberto: text continues here:}

To find the time evolution of the scale factor $a$ 
we solve Friedmann equations for a generic equation
of state $p = w \rho$,
\begin{equation}
    %\frac{\ddot{a}}{a} = - \frac{4\pi G}{3} \bigl(\rho + 3 p) =
    %- \frac{4\pi G}{3} \rho \bigl(1 + 3w\bigr),
    \frac{\ddot{a}}{a} = - \frac{\kappa}{6} \, \bigl(\rho + 3 p) =
    - \frac{\kappa}{6} \, \rho \bigl(1 + 3w\bigr),
    %\label{friedmann1} %\\ % \Rightarrow \rho = \rho_0 \, a^{-3}
    %e^{-3 \int a^{-1} w(a) \dd a},
    \quad 
    \biggl(\frac{\dot a}{a} \biggr)^2 =
    %\frac{8\pi G}{3} \rho.
    \frac{\kappa}{3} \, \rho.
    \label{friedmann}
    % \Rightarrow a' = 
    % C a^{1\over2} e^{-{3\over 2}
    % \int a^{-1} w(a) \dd a},
\end{equation}
%where $C$ is a constant that depends on $\rho_0$.
The total energy density can be expressed in terms
of $a$ as
\begin{equation}
    %\rho (a) = \frac{3 H_0^2}{8\pi G} \, \Omega (a),
     \rho (a) = \rho_{{\rm c}, 0} \, \Omega (a) = 
     \frac{3 H_0^2}{\kappa} \, \Omega (a),
\end{equation}
where $\rho_{{\rm c}, 0}$ and $H_0$ are the critical energy density
and the Hubble rate, respectively, at present time.
This 
allows us to write the Friedmann equations as
\begin{equation}
    \frac{\ddot a}{a} = - \frac{1}{2} H_0^2 \, \Om (a)
    \bigl[1 + 3 w (a)\bigr],
    \quad
    \biggl(\frac{\dot a}{a} \biggr)^2 =
    H_0^2 \, \Om(a).
    \label{friedmann2}
\end{equation}
The ratio of total energy to current critical energy
densities is
\begin{equation}
\Omega (a) = \frac{\rho (a)}{\rho_{{\rm c}, 0}} =
    \Omega_\rad + \Omega_\mat + \Omega_\Lam = 
     \biggl({a\over a_0}\biggr)^{-4} \,
    \frac{g_*}{g_*^0} \biggl(\frac{g_{\rm S}}{g_{\rm S}^0}
    \biggr)^{-{4\over3}} \, \Omega_{\rad,0} + \biggl({a\over a_0}\biggr)^{-3} \, \Omega_{\mat,0} + \Omega_{\Lambda,0},
    \label{rho_a}
\end{equation}
with $\Om_{\Lam,0} \simeq 0.6841$,
%$\Om_{\rad,0} = 9.28656\times10^{-5}$~\cite{Planck18}.
$\Om_{\rad,0} \simeq 4.16 \times10^{-5} \, h^{-2}$, and
 $\Om_{\mat,0} = 1 - \Om_{\rad, 0} - \Om_{\Lam, 0}
 \simeq 0.316$
%$\Om_{\rad,0} = 5.5088\times10^{-5}$
(for $H_0 = 100 \, h$ km/s/Mpc)~\cite{Planck18},
%
% The radiation energy density is omitting the variation of the
% relativistic and adiabatic degrees of freedom in the early universe,
% which leads to
% \begin{equation}
%     \Om_\rad (a, g_*, g_{\rm S}) = a^{-4} \, \Om_{\rad, 0} \,
%      \frac{g_*}{g_*^0} \biggl(\frac{g_{\rm S}}{g_{\rm S}^0}
%     \biggr)^{-{4\over3}},
% \end{equation}
%AR: added now to equation above
%being $g_*^0 \simeq 3.36$ and $g_{\rm S}^0 \simeq 3.91$ the
and $g_*^0 \simeq 3.36$ and $g_{\rm S}^0 \simeq 3.91$ are 
the reference
%relativistic and adiabatic degrees of freedom at present time
relativistic and adiabatic DOFs at present time
[see section about radiation-dominated (RD) era
below for further details, in particular \Eq{Omrad0}].

% \begin{figure}
%     \centering
%     \caption{Ratio of total to
%     present time
%     critical energy density $\Omega$ as a function of
%     the scale factor $a$.
%     It takes into account the dependence of $\Om(a)$ on the
%     relativistic and adiabatic degrees of freedom during the
%     radiation-dominated era and it uses $h = 0.6732$.}
%     \label{rho_vs_a}
% \end{figure}
\begin{figure}
    \centering
    \includegraphics[width=.49\textwidth]{rho_vs_a}
    \includegraphics[width=.49\textwidth]{g_vs_a}
    %\includegraphics[width=.49\textwidth]{rho_vs_a_rat}
    %\includegraphics[width=.6\textwidth]{rho_vs_a_comp}
    \caption{{\em Left panel: }ratio of total to
    present time
    critical energy density $\Omega$ as a function of
    the scale factor $a$.
    It takes into account the dependence of $\Om(a)$ on the
    %relativistic and adiabatic degrees of freedom during the
    relativistic and adiabatic DOFs during the
    radiation-dominated era (right panel)
    and it uses $h = 0.6732$.}
    % {\em Right panel:} evolution with the scale factor $a$
    % of the relativistic $g_*$ and adiabatic
    % $g_{\rm S}$ degrees of freedom in the early universe
    % (left panel) and ratio that needs to be multiplied to 
    % $a^{-4} \, \Om_{\rad, 0}$ to take into account their effect
    % on the radiation energy density (right panel),
    % ignoring the neutrinos decoupling at later times.}}
    \label{Om_vs_a}
\end{figure}

\boxbegin
The radiation energy today consists of CMB photons and relic neutrinos.
So the total radiation energy density is $\Om_{\rad,0} = \Om_{\gamma,0} + \Om_{\nu,0}$, with
\begin{equation}
\Om_{\gamma,0} = \frac{\rho_{\gamma,0}}{\rho_\crit},\quad
\Om_{\nu,0}    = \frac{\rho_{\nu,0}}{\rho_\crit}.
\end{equation}
The photon energy $\rho_{\gamma,0}$ can be calculated via
\begin{equation}
\rho_{\gamma,0} = \frac{\aB T_\CMB^4}{c^2},
\end{equation}
where $\aB = \pi^2\kB^4/(15\hbar^3c^3)\approx7.5657\times10^{-16}\rm Jm^{-3}K^{-4}$ is the radiation constant,
calculated from Boltzmann constant $\kB$
and reduced Planck constant $\hbar$.
Given $T_\CMB\approx2.726\rm K$, $\rho_{\gamma,0}$ can be calculated to further give $\Om_{\gamma,0}\approx4.93812\times10^{-5}$.
For neutrinos, their energy $\rho_\nu$ can be related to $\rho_\gamma$ as~\cite{Planck18}
\begin{equation}
\rho_\nu = N_\eff\Big(\frac{7}{8}\Big)\Big(\frac{4}{11}\Big)^{4/3}\rho_\gamma,
\end{equation}
%where $N_\eff\sim3.046$ is today's relativistic degrees of freedom (DOF).
where $N_\eff\sim3.046$ is today's relativistic DOFs.
So $\Om_{\nu,0}\approx3.41604\times10^{-5}$.
From here we have $\Om_{\rad,0} = \Om_{\gamma,0} + \Om_{\nu,0}\approx8.35417\times10^{-5}$.
\boxend
\FFig{Om_vs_a} shows the function $\Omega (a)$, computed
from \Eq{rho_a}, using $h = 0.6372$ \cite{Planck18}
(left panel) and the evolution of the degrees
of freedom during the RD era (right panel).
We find that it evolves approximately
proportional to $a^{-4}$ during 
the RD era, to $a^{-3}$ during
matter-dominated (MD) era, and it is
constant during dark energy
domination ($\Lambda$D).
%During RD, the evolution of degrees of freedom induces some
During RD, the evolution of DOFs induces some
modifications with respect to the $a^{-4}$ evolution.
In the transition between different eras, $\Om (a)$ evolves
smoothly among the asymptotic behaviors.
%Note that we take here the convention $a_0 = 1$.

The characteristic equation of state $w (a)$ corresponding to
the energy density of \Eq{rho_a} can be computed combining
Friedmann equations [see \Eq{friedmann2}],
% \begin{equation}
%     \frac{\ddot a}{a} = - \frac{1}{2} \Omega(a) H_0^2
%     \bigl(1 + 3w\bigr), \quad 
%     \biggl(\frac{\dot a}{a} \biggr)^2 = H_0^2 \, \Omega (a),
%     \label{friedmann2}
% \end{equation}
by taking the time derivative of the second equation and
substituting the first equation.
This yields
\begin{gather}
    % 2 \, \frac{\dot a}{a} \biggl(\frac{\ddot a}{a} - \frac{\dot a^2}{a^2} \biggr) = H_0^2 \, \dot \Omega (a)
    % \Rightarrow
    \frac{\dot \Om (a)}{\Om (a)} = 
    -3 (1 + w) \frac{\dot a}{a},
    \label{Om_a_dot}
\end{gather}
which allows us to compute $w (a)$ using \Eq{rho_a},
\begin{equation}
    w (a)= \frac{{1 \over 3} \Om_{\rad} (a) 
    - \Om_{\Lam, 0}}{\Om (a)},
    %= \frac{{1 \over 3} a^{-4} \Om_{\rad}^0}{\Om_{\rad}^0 
    %a^{-4} + \Om_{\mat}^0 a^{-3}}
\end{equation}
% \begin{equation}
%     w (a) - 1/3 = \frac{{1 \over 3} \Om_{\rad} (a) 
%     - \Om_{\Lam, 0}}{\Om (a)} = \frac{{1 \over 3} a^{-4} \Om_{\rad}^0}{\Om_{\rad}^0 a^{-4} + \Om_{\mat}^0 a^{-3}}
%     - \frac{1}{3} = \frac{-{1\over3} \Om_{\mat}^0 a^{-3}}{\Om_{\rad}^0 a^{-4} + \Om_{\mat}^0 a^{-3}} \approx -
%     {1\over3} \frac{\Om_{\mat, 0}}{\Om_{\rad, 0}} a,
% \end{equation}
and is shown in \Fig{w_vs_a}.
We find the ultrarelativistic equation of
state during RD era $w = 1/3$, $w = 0$ for MD, and
$w = -1$ for $\Lambda$D, as expected.
In the intermediate times, the function $w (a)$ transitions
smoothly.
Using $w (a)$, we can compute the scale factors
at which equipartition occurs, $a_\EQ$, and when dark energy becomes
dominant, $a_\Lam$, by taking $w (a_\EQ) = 1/6$ and $w (a_\Lam) = -0.5$, which yields (for $h = 0.6372$)
\begin{equation}
    \frac{a_\EQ}{a_0} = 2.9 \times 10^{-4}, \quad
    \frac{a_\Lam}{a_0} = 0.772.
\end{equation}


\begin{figure}
    \centering
    \includegraphics[width=.49\textwidth]{w_vs_a}
    \includegraphics[width=.5\textwidth]{w_vs_a_rad_1o3}
    \caption{{\em Left panel:} equation of state $w$ such that $p = w
    \rho$ at different values of the scale factor
    $a$ computed from Friedmann equations.
    The solid and dashed lines use $h = 0.6732$ and $h = 1$,
    respectively.
    {\em Right panel: }Deviation from $w = 1/3$ in the equation
    of state during RD era.
    The solid line shows the exact result and the dashed
    line the approximation given in 
    \Eq{Delta_w}, ignoring the variation of degrees of
    freedom.}
    \label{w_vs_a}
\end{figure}

We now compute the time evolution of $a$, using Friedmann
equations [see \Eq{friedmann2}]: 
\begin{equation}
    \dot{a} = H_0 \, a \sqrt{\Omega (a)} \Rightarrow
    \dd \bigl(H_0 t\bigr) = \frac{\dd a}{a \sqrt{\Omega (a)}}
    \Rightarrow H_0 \bigl( t - t_{\rm ini} \bigr) =
    \int_{a_{\rm ini}}^a
    \frac{\dd a}{a \sqrt{\Omega (a)}}.
    \label{dota}
\end{equation}
This allows us to compute $t (a)$ and then invert the
relation to obtain $a (t)$.
Similarly, in conformal time we have $\dot a = a'/a$, so we solve the following integral and invert to obtain $a (\eta)$.
\begin{equation}
    H_0 \bigl( \eta - \eta_{\rm ini} \bigr) \, a_0 = a_0 \int_{a_{\rm ini}}^a
    \frac{\dd a}{a^2 \sqrt{\Omega (a)}},
    \label{dota}
\end{equation}
where we include the reference value $a_0$ to take into account
general convention in which, in general, $a_0 \neq 1$.
In \Eq{dota}, both terms are independent on the specific 
normalization of the scale factor.

\FFig{a_vs_t} shows the time evolutions $a (t)$ and $a (\eta)$.
The initial times $t_{\rm ini}$ and $\eta_{\rm ini}$ are set to 
zero at $a_{\rm ini} = 10^{-20} a_0$ for the numerical integration,
much earlier than the EWPT ($a \sim 10^{-16} a_0$).
Note that $t$ and $\eta$ here correspond to
physical times, and the latter is expressed as $a_0 \, \eta$,
which allows us to scale it when we use a convention
with $a_0 \neq 1$.
From the solution to Friedmann equations, we can compute the
time of equipartition, dark energy domination, and present time
(using $h = 0.6372$ and $a_0 = 1$)
\begin{align}
    t_\EQ \simeq 3.46 \times 10^{-6}/H_0, &\, \quad
    t_\Lam \simeq 0.71/H_0, \quad
    t_0 \simeq 0.95/H_0, \\
    \eta_\EQ \simeq 0.025/H_0, &\, \quad
    \eta_\Lam \simeq 2.90/H_0, \quad
    \eta_0 \simeq 3.18/H_0.
\end{align}

\begin{figure}
    \centering
    \includegraphics[width=.49\textwidth]{a_vs_eta}
    \includegraphics[width=.49\textwidth]{a_vs_t}
    \caption{Time evolution of the scale factor with
    conformal
    (left panel) and cosmic (right panel) times.
    The solid and dashed lines use $h = 0.6732$ and $h = 1$,
    respectively.}
    \label{a_vs_t}
\end{figure}

The derivatives of the scale factor with respect to 
conformal time that appear in the GW \Eq{modGR}
can be directly computed from Friedmann
equations [see \Eq{friedmann2}],
\begin{equation}
    a_0^{-1} \, \HH = H_0\, \frac{a}{a_0} \sqrt{\Omega (a)},
    \quad 
    a_0^{-2} \, \frac{a''}{a} = \frac{1}{2}
    \biggl(\frac{a}{a_0}\biggr)^2 \Omega (a)
    \bigl[1 - 3 w (a) \bigr] H_0^2.
    \label{eqn:a_derivatives_omega}
\end{equation}
In these expressions, we have included powers of $a_0$ to allow
for general convention $a_0 \neq 1$.
The functions $a'/a$ and $a''/a$ are shown in \Fig{ap_app_eta}
and can be directly rescaled to normalizations with $a_0 \neq 1$.

\begin{figure}
    \centering
    \includegraphics[width=.45\textwidth]{ap_vs_eta}
    \includegraphics[width=.45\textwidth]{app_vs_eta}
    \caption{Variation of $\HH = a'/a$ (left panel) and $a''/a$ 
    (right panel) with the scale factor $a$.
    % Both functions are computed using $h = 0.6372$ and taking
    % into account the modification of degrees of freedom in the
    % early universe.}
    Both functions are computed using $h = 0.6372$.}
    \label{ap_app_eta}
\end{figure}
%
% Additionally, \Eq{Om_a_dot} allows us to compute $\Om(a)$
% and $a(t)$ analytically
% at the different epochs by taking $w = {1\over3}$, $0$, $-1$
% for the RD, MD and $\Lam$D eras, respectively. We study each
% of these epochs in the following subsections.

\subsection{Radiation-dominated (RD) era}

In this section, we study the solution to Friedmann equations
in the early Universe, when it is dominated by
radiation.
In the current work, we will consider scenarios in which the
GW background is produced at some time during the RD era, so
we present here some approximations that hold during this 
epoch, and compare the results to those obtained by solving
numerically Friedmann equations.
% Hence, the normalization presented in \Eq{modGR} is justified
% for this type of simulations, as we show below.

During the RD era, $w \approx {1\over 3}$ [see \Eq{Om_a_dot} and
\Fig{w_vs_a}].
Hence, using \Eq{Om_a_dot},
we find
\begin{equation}
    \frac{\dd \Om}{\Om} = - 4 \frac{\dd a}{a} \Rightarrow
    \Om_\rad (a) = \Om_{\rad} (a_*) 
    \biggl(\frac{a}{a_{*}}\biggr)^{-4} \approx
    \biggl(\frac{H_{*}}{H_0}\biggr)^2 \biggl(\frac{a}{a_{*}}\biggr)^{-4},
    \label{Om_rad}
\end{equation}
where $\Om_{\rad} (a_*)$ is the energy density at the reference scale factor
$a_{*}$, which belongs to the RD epoch, i.e., $a_*/a_0 \lesssim 10^{-5}$
[see \Fig{Om_vs_a} or \Fig{w_vs_a}].
Since radiation dominates, then we can assume that $\Om_{\rad} (a_*) \approx \Om (a_*)$.
Now we use \Eqs{friedmann2}{Om_rad} to find
\begin{equation}
    \frac{\dot a}{a} \approx H_{*}
    \biggl(\frac{a}{a_{*}}\biggr)^{-2}
    \Rightarrow a \dd a = H_{*} a_{*}^2 \Rightarrow
    \frac{a}{a_{*}} =
    \sqrt{\biggl(\frac{a_\ini}{a_*}\biggr)^2 +
    2 H_{*} \bigl(t - t_{\ini} \bigr)}.
\end{equation}
We can choose $t_\ini = a_\ini = 0$ as the integration
constants to get
\begin{equation}
    \frac{a}{a_*} = \sqrt{2 H_* t}.
\end{equation}
This allows us to compute $t_* = (2 H_*)^{-1}$.
%
% We can choose $t_{\rad, *}^{-1} = 2 H_{\rad, *}$
% as the integration constant
% to simplify the previous expression to
% \begin{equation}
%     \frac{a_\rad (t)}{a_{\rad, *}} = \sqrt{2 H_{\rad, *} t} = 
%     \sqrt{t/t_{\rad, *}} = \sqrt{\bar t},
% \end{equation}
% where we have defined the normalized
% $\bar t = t/t_{\rad, *}$.
%
In terms of the conformal time, we find
\begin{equation}
    \dd a = H_{*} \, a_{*}^2 \dd \eta \Rightarrow 
    \frac{a}{a_{*}} = 
    \frac{a_\ini}{a_*} + \HH_{*} \bigl(\eta - \eta_{\ini}
    \bigr),
\end{equation}
where we have used the conformal Hubble rate $\HH_{*} = a_{*} H_{*}$.
Again, choosing the integration constant
$t_\ini = a_\ini = 0$, we find
\begin{equation}
    \frac{a}{a_*} = \HH_* \eta = \bar \eta.
\end{equation}
% such that $\eta_{\rad, *} = \HH_{\rad, *}^{-1}$
% we find
% \begin{equation}
%     \frac{a_\rad (\eta)}{a_{\rad, *}} = \eta/\eta_{\rad, *} = \bar \eta,
% \end{equation}
Hence, during RD era, the conformal time
is $\eta = \HH^{-1}$ and, in particular, $\eta_* = \HH_*^{-1}$.
This justifies the use of the normalization introduced
in \Sec{modGR} [see \Eq{GW_eq}].
%
% where 
% we have defined the normalized $\bar \eta = \eta/\eta_{\rad, *}$.
% Then, the energy density evolves as
% \begin{equation}
%     \Om_\rad = \Om_{\rad, *} \bar t^{-2} = 
%      \Om_{\rad, *} \bar \eta^{-4}.
%     \label{Omrad_bareta}
% \end{equation}

Note that the presented calculations in this section ignore
the variations in the number of relativistic
%degrees of freedom during the RD epoch and their effect in
DOFs during the RD epoch and their effect in
Friedmann equations.
In general, the radiation energy density can be computed as
\begin{equation}
    \Om_\rad (T, g_*) =
    \frac{\pi^2}{30} \frac{g_* T^4}{\rho_{\rm c, 0}}
    \simeq 4.064 \times 10^{55} \, h^{-2} \frac{g_*}{100}
    \biggl(\frac{T}{100 \GeV} \biggr)^4,
    \label{rad_Om}
\end{equation}
where $g_*$ is the number of degrees of the relativistic
species and
we have taken the reference values $T \sim 100 \GeV$ and
$g_* \sim 100$, corresponding to the EWPT.
The value $\Om_{\rad, 0}$ at present time can
be computed using \Eq{rad_Om},
\begin{equation}
    \Om_{\rad, 0} = \frac{\pi^2}{30}
    \frac{g_*^0 T_0^4}
    {\rho_{\rm c}} = \frac{4 \pi^3 G}{45}
    \frac{g_*^0 T_0^4}{H_0^2} \simeq 2.47
    \times 10^{-5} \, h^{-2} \, \frac{g_*^0}{2},
\end{equation}
where present time values are $T_0 = 2.75$ K
and $g_*^0 = 2$.
This value only takes into account photons, and neglects
the contribution from neutrinos $\Om_\nu$, which
contribute to the radiation energy
at times before neutrinos become massive.
Since during RD era, neutrinos are relativistic, we 
want to incorporate them by redefining the relativistic
%degrees of freedom (using $N_{\rm eff} = 3$),
DOFs (using $N_{\rm eff} = 3$),
\begin{equation}
    g_*^0 = 2 \Biggl[1 + \, N_{\rm eff} \, \frac{7}{8} \biggl(\frac{4}{11}
    \biggr)^{4\over3} \Biggr] \simeq 3.36.
\end{equation}
This yields $\Om_{\rad, 0} \simeq 4.16 \times 10^{-5} \, h^{-2}$, used in \Eq{rho_a} as the reference value.
Note that this leads to an error in the actual radiation
energy density after neutrinos become massive.
However, at this and later times,
the contribution of radiation
to the total energy density is negligible and hence,
it does not affect to our calculations.

The radiation energy density [see \Eq{rad_Om}] can
be expressed in terms of $a$ if we take into account
the adiabatic expansion of the Universe, such that
\begin{equation}
    \frac{a}{a_0} = \frac{T_0}{T} \biggl(\frac{g_{\rm S}^0}
    {g_{\rm S} (T)} \biggr)^{1\over 3} \Rightarrow
    \biggl(\frac{T}{100 \GeV}\biggr)^4 \simeq 4.038 \times 
    10^{-61} \biggl(\frac{a}{a_0}\biggr)^{-4} \biggl(
    \frac{100}{g_{\rm S}}\biggr)^{4\over3},
    \label{T_vs_a}
\end{equation}
where the present time value of the adiabatic
%degrees of freedom is $g_{\rm S}^0 = 3.91$.
DOFs is $g_{\rm S}^0 = 3.91$.
%and we
%use $g_{\rm S} \approx g_*$, which holds at large $T$.
Substituting this into \Eq{rad_Om}, we find
\begin{equation}
    \Om_\rad (a, g_*, g_{\rm S})
    \simeq 1.64 \times 10^{-5} \, h^{-2}
    \biggl(\frac{a}{a_0}\biggr)^{-4}
    \biggl(\frac{g_*}{100}\biggr)
    \biggl(\frac{g_{\rm S}}{100}\biggr)^{-{4\over 3}}.
    \label{Omrad2}
\end{equation}
This expression allows us to use $g_* = g_*^0 \simeq
3.36$, $g_{\rm S} = g_{\rm S}^0 \simeq 3.91$, and $a = a_0$ to find
\begin{equation}
    \Om_{\rad} (a) = \biggl(\frac{a}{a_0}
    \biggr)^{-4} \frac{g_*}{g_*^0} 
    \biggl(\frac{g_{\rm S}}{g_{\rm S}^0} \biggr)^{-{4\over3}}
    \, \Om_{\rad, 0},
    \label{Omrad0}
\end{equation}
which has been used in \Eq{rho_a}.
% If we compare \Eq{Omrad2} to \Eq{Om_rad}, 
% \begin{equation}
%     \Om_\rad (\bar \eta, g_*, g_{\rm S}) =
%     \Om_{\rad} (a_*) \, \bar \eta^{-4}
%     \biggl(\frac{g_*}{g_* (a_*)} \biggr)
%     \biggl(\frac{g_{\rm S}}{g_{\rm S} (a_*)}
%     \biggr)^{-{4\over 3}},
%     \label{Omrad_eta_corr}
% \end{equation}
% where
% \begin{equation}
%     \Om_{\rad, *} = \biggl(\frac{H_*}{H_0} \biggr)^2
%     \simeq 1.64 \times 10^{-5} \, h_0^{-2} \biggl(
%     \frac{a_*}{a_0}\biggr)^{-4}.
% \end{equation}
Since $g_*$ and $g_{\rm S}$ evolve very slowly (compared to $a^{-4}$),
ignoring this dependence
leads to small errors.
In \Fig{errors_rad}, we compare the exact results of $\Om(a)$, $a(\eta)$,
and $a(t)$ obtained by solving
numerically Friedmann equations with those obtained
analytically, by ignoring the variations
%of the degrees of freedom.
of the DOFs.
% and just correct \Eq{Omrad_bareta} by including the factors $g_*/g_*^0$
% and $(g_{\rm S}/g_{\rm S}^0)^{-{4\over3}}$
% [see \Eq{Omrad_eta_corr}].
% The variation of $g_*$ and $g_{\rm S}$ with $a$, the correction
% added to $\Om_{\rad}$ and the comparison of the latter with the 
% one ignoring the evolution of $g$ are shown in \Fig{Om_vs_a}.
In particular, the right panel of \Fig{errors_rad} shows the
error in assuming $\eta = \HH^{-1}$ during the RD era due to
%the time variations of the degrees of freedom.
the time variations of the DOFs.
However, this error is never larger than 0.2 during the RD 
epoch.

During the RD era, the second conformal time derivative of $a$ 
is identically zero if we take $w = {1 \over 3}$
[see \Eq{eqn:a_derivatives_omega}].
The deviations from this value
(due to the matter content) are
\begin{equation}
    \Delta w = 1 - 3 w(a) =  3 \, \frac{\Om_\mat}{\Om_\rad
    + \Om_\mat} \simeq 3\, \frac{\Om_\mat}
    {\Om_\rad} = 3 \, \frac{a}{a_0}
    \frac{g_*^0}{g_*} \biggl(\frac{g_{\rm S}^0}{g_{\rm S}}
    \biggr)^{-{4\over3}} \, \frac{\Om_{\mat, 0}}{\Om_{\rad, 0}}.
    \label{Delta_w}
\end{equation}
These deviations are small during RD era [see \Fig{w_vs_a}]
and are usually neglected.
The resulting second conformal time derivative can be obtained
plugging \Eq{Delta_w} into \Eq{eqn:a_derivatives_omega},
to find
\begin{equation}
    a_0^{-2}  \, \frac{a''}{a} \approx \frac{3}{2} \biggl(
    \frac{a}{a_0} \biggr)^2 \Om_\mat \, H_0^2 
    = \frac{3}{2} \biggl(\frac{a}{a_0}\biggr)^{-1}
    \Om_{\mat, 0} \, H_0^2.
    \label{app_dev}
\end{equation}
Hence, $a''/a$ decays as $a^{-1}$ with the scale factor or
as $\eta^{-1}$ with conformal time, as it is shown in \Fig{ap_app_eta}.


\begin{figure}
    \centering
    \includegraphics[width=.49\textwidth]{errors_rad}
    \includegraphics[width=.49\textwidth]{Heta_vs_a}
    \caption{{\em Left panel:} ratio of the conformal and cosmic
    times obtained solving Friedmann equations numerically
    and analytically by neglecting the variation of degrees
    of freedom (no g).
    The modification on the energy density is $\Om_\rad = 
    \Om_\rad^{\rm no \ g} (g_*/g_*^0) (g_{\rm S}/g_{\rm S}^0
    )^{-{4\over3}}$.
    {\em Right panel: } Deviations of the numerical results
    of $\HH \eta$ from 1 due to the inclusion of time
    %dependent degrees of freedom.}
    dependent DOFs.}
    \label{errors_rad}
\end{figure}

\subsection{Initial conditions for the GW equation}

We have presented the GW equation in \Sec{GW_eq},
and the normalized GW \Eq{modGR}.
To solve the GW equation, we first choose the scale factor
$a_*/a_0$ at which we want to start our simulation.
Then, we find $\HH_*$ from the numerical solution to Friedmann
equations (see \Fig{ap_app_eta}).
We set as the initial time (normalized conformal time)
of our simulation $\bar \eta_* = 1 + \varepsilon$.
We have shown that during the RD era,
$\varepsilon$ is a small
correction [see right panel of \Fig{errors_rad}] and, hence,
$\bar \eta_* \approx 1$.
We then evolve the equation in $\bar \eta$ and at each time
step we compute the normalized $\bar \HH = \HH/\HH_*$ and
$a''/a$ (normalized by $\HH_*^{-2}$).
Finally, we set the final time of our simulations our present time
$\bar \eta_0 = \eta_0 \, \HH_*$.
Hence, the GW equation is evolved along the different eras using
the smoothed values obtained by solving Friedmann equations.

In previous simulations, we have ignored the term $a''/a$
during the RD era since it is very small compared to the other
terms of the GW equation.
This can be seen by estimating the value of the term
$h_\ij''$ in \Eq{eqn:GW_conf_comove} as $h_\ij''/h_\ij \sim 
\eta^{-2}$, which is
\begin{equation}
    \eta^{-2} \approx \HH^2 = a^2 H^2 \approx
    a^2 H_0^2 \, \Om_\rad \sim \biggl(\frac{a}{a_0}\biggr)^{-2}
    H_0^2 \, \Om_{\rad, 0}
\end{equation}
The term $a''/a$ can be estimated using \Eq{app_dev}.
Hence, the ratio of this term over $\eta^{-2}$ is of order
\begin{equation}
    \frac{a''}{a} \eta^2 \sim \frac{3}{2}
    \frac{a}{a_0} \frac{\Om_{\mat, 0}}{\Om_{\rad, 0}}
    \simeq 5 \times 10^3 \frac{a}{a_0}.
    \label{app_relevance}
\end{equation}
Therefore, the term $a''/a$
is smaller than 1 and decreases linearly when $a/a_0 \lesssim
2 \times 10^{-4}$.
Hence, it is in general justified to omit this value during the
RD era.
In any casse, in the current work, we keep its value,
obtained from the numerical solution to the Friedmann
equations [see right panel of \Fig{ap_app_eta}] to have a
smooth function along our simulation and avoid discontinuities.
In addition, we run the simulations up to present time and, hence,
we will evolve $\bar \eta$ into the MD and $\Lambda$D eras.

Let us consider the case in which the initial time of the
simulations corresponds to the electroweak phase transition
at a temperature scale $T = 100 \GeV$.
At this $T$, $g_{\rm S} \simeq 103.8$
and $a_*/a_0 \simeq 7.87 \times 10^{-16}$ [from \Eq{T_vs_a} and
right panel of \Fig{Om_vs_a}].
From the numerical solutions to Friedmann equations we find
$\HH_* \simeq 1.66 \times 10^{-5} \,$Hz and $\varepsilon
\simeq 3.4 \times 10^{-3}$, and the present time in units of
$\HH_*$ is $\bar \eta_0 \simeq  2.4 \times 10^{13}$.
The normalized time $\bar \eta$ as a function of
$a/a_*$, and the evolution of $\bar \HH$ and $\HH_*^{-2}
(a''/a)$ with $\bar \eta$ are shown in \Fig{EWPT}, where present
time is indicated as $a_0$.
In the plot of $a''/a$ it is shown $\bar \eta^{-2}$ for comparison,
which corresponds to the amplitude of the second time derivatives.
It can be seen that the former term is much smaller than $\bar 
\eta^{-2}$ until we reach equipartition, as estimated in
\Eq{app_relevance}.

\begin{figure}
    \centering
    \includegraphics[width=.49\textwidth]{a_vs_eta_EWPT}\\
    \includegraphics[width=.49\textwidth]{Hbar_vs_eta_EWPT}
    \includegraphics[width=.49\textwidth]{app_vs_eta_EWPT}
    \caption{Normalization taking the electroweak phase
    transition scale as the initial time $\eta_*$.
    {\em Upper panel:} evolution of the scale
    factor $a/a_*$ in terms of the normalized $\bar \eta = \HH_*
    \eta$. {\em Lower panels:} evolution of the normalized
    $\bar \HH = \HH/\HH_*$ (left panel) and $\HH_*^{-2} (a''/a)$
    (right panel)
    as a function of $\bar \eta$.}
    \label{EWPT}
\end{figure}

\subsubsection*{Matter-dominated (MD) era ($w = 0$)}

In this case, we have
\begin{equation}
    \frac{\dd \Om}{\Om} = -3 \frac{\dd a}{a} \Rightarrow 
    \Om_{\mat} (a) = \Om_{\mat, *}
    \biggl(\frac{a}{a_{\mat, *}}\biggr)^{-3} = 
    \biggl(\frac{H_{\mat, *}}{H_0}\biggr)^2
    \biggl(\frac{a}{a_{\mat, *}}\biggr)^{-3},
\end{equation}
where the reference scale factor $a_{\mat, *}$
and energy density $\Om_{\mat, *}$ belong to the
matter-dominated (MD) era, c.f., $a/a_0 \in (\sim 10^{-3}$--$0.3)$, and, again, we can assume that 
$\Om_{\mat, *} = \rho_{\rm crit}^*/\rho_{\rm crit}^0$.
The evolution of the scale factor with time is
\begin{equation}
    \frac{\dot a}{a} = H_{\mat, *} \biggl(
    \frac{a}{a_{\mat, *}}\biggr)^{-{3\over2}} \Rightarrow 
    \frac{a}{a_{\mat, *}} = \biggl[1 + \twothird
    H_{\mat, *} \bigl(t - t_{\mat, *} 
    \bigr)\biggr]^{2\over3}.
\end{equation}
We choose the integration constant such that
$t_{\mat, *}^{-1} = \threehalf H_{\mat, *}$ to write
\begin{equation}
    \frac{a_\mat (t)}{a_{\mat, *}} = \biggl(\frac{t}{t_{\rm mat, *}}\biggr)^{2\over 3}.
\end{equation}
In terms of conformal time, we have
\begin{equation}
    \frac{\dd a}{a^{1\over2}} = \HH_{\mat, *}
    a_{\mat, *}^{{1\over2}} \dd \eta \Rightarrow
    \frac{a}{a_{\mat, *}}
    =\biggl[ 1 + \half \HH_{\mat, *}
    \bigl(\eta - \eta_{\mat, *}\bigr) \biggr]^2.
\end{equation}
The choice of the integration constant such that $\eta_{\mat, *}^{-1} =
\half \HH_{\mat, *}$, yields
\begin{equation}
    \frac{a (\eta)}{a_{\mat, *}} =
    \biggl(\frac{\eta}{\eta_{\mat, *}}
    \biggr)^2.
\end{equation}

\subsubsection*{Dark energy-dominated ($\Lambda$D)
era ($w = -1$)}

In the dark energy-dominated era, which is only approached
at $a_*/a_0 \gtrsim 10$, we have
\begin{equation}
    \frac{\dd \Om}{\Om} = 0 \Rightarrow \Om_{\Lam} (a)
    = \Om_{\Lam, *} = \biggl(\frac{H_{\Lam, *}}{H_0}
    \biggr)^2,
\end{equation}
where $\Om_{\Lam, *}$ and $H_{\Lam, *}$ are constant
values with time.
The scale factor evolves as
\begin{equation}
    \frac{\dot a}{a} = H_{\Lam, *} \Rightarrow 
    \frac{a(t)}{a_{\Lam, *}} = e^{H_{\Lam, *} (
    t - t_{\Lam, *})} = e^{\bar t - 1},
\end{equation}
where we have used $t_{\Lam, *} = H_{\Lam, *}^{-1}$
and defined $\bar t = t/t_{\Lam, *}$.
In terms of the conformal time,
\begin{equation}
    \frac{a'}{a^2} = H_{\Lam, *} \Rightarrow
    \frac{a (\eta)}{a_{\Lam, *}} = \frac{1}{1 -
    \HH_{\Lam, *} \bigl(\eta - \eta_{\Lam, *}\bigr)}
     = \frac{1}{2 - \bar \eta},
\end{equation}
where we hace chosen $\eta_{\Lam, *} = \HH_{\Lam, *}^{-1}$
and defined $\bar \eta = \eta/\eta_{\Lam, *}$.

% We can then take that the equation of state as a function of $a$
% is
% \begin{equation}
%     w = \frac{{1\over3} \Omega_\rad(a) - \Omega_\Lambda (a)}
%     {\rho (a)},
% \end{equation}
% which takes into account that radiation has an equation of
% state $p = {1\over3} \rho$, matter $p = 0$, and dark energy
% $p = - \rho$.
% \blue{
% How did you arrive at Eq.~(5.12)?
% And do we need Eq.~(5.12)?
% I thought we would still put in $w = 1/3, 0, -1$
% by hand at the right time,
% but solve for $a(\eta)$ from Friedmann?
% Because Eq.~(5.12) would still require us knowing $a(\eta)$ to begin with, right?
%}
%\begin{equation}
%a' = \sqrt{\frac{8\pi\GN}{3}\rho}a^2,\quad
%a'' = \frac{(a')^2}{a} - \frac{4\pi\GN}{3}a^3\rho(1 + 3w),
%\end{equation}
%where $\rho$ is Eq.~\ref{eqn:rho_a}.
%Initial conditions: $\eta_* = 1, a_* = a_*' = 1$.
%\begin{figure}[H]
%\centering
%\includegraphics[scale = 0.55]{Figures/rho_w_a.png}\\
%\includegraphics[scale = 0.55]{Figures/a_t_eta.png}
%\caption{Dots in the upper panels indicate today.}
%\end{figure}

\section{Evolution matrix}

Similar to what has been done in Ref.~\cite{RoperPol+18} and Ref.~\cite{BS21},
we implement a matrix solver for the modified GW equation.
We start by rewriting Eqn.~\ref{eqn:GW_mod_rad} as
\begin{equation}
h_\ij'' + \sig h_\ij' + \om^2h_\ij = \ggg T_\ij,
\label{eqn:GW_mod_rewrite}
\end{equation}
where $\ggg\equiv6/\eta$, and
\begin{equation}
\sig\equiv\alpM\hhh\quad,\quad
\om^2\equiv(1 + \alpT)k^2 - \alpM\hhh^2 - \frac{a''}{a}.
\end{equation}
Then from the characteristic equation
\begin{equation}
\lam^2 + \sig\lam + \om^2 = 0,
\end{equation}
the eigenvalues $\lam$ can be obtained as
\begin{equation}
\lam_{1,2} = -\frac{1}{2}\Big(\sig\mp\sqrt{\sig^2 - 4\om^2}\Big).
\end{equation}
Then we take the ansatz
\begin{align}
h_\ij(\eta + \del\eta) & = h_\ij(\eta) + \ccc_\ij e^{\lam_1\del\eta} + \ddd_\ij e^{\lam_2\del\eta},\\
h_\ij'(\eta + \del\eta) & = h_\ij'(\eta) + \ccc_\ij\lam_1 e^{\lam_1\del\eta} + \ddd_\ij\lam_2 e^{\lam_2\del\eta},\\
h_\ij''(\eta + \del\eta) & = h_\ij''(\eta) + \ccc_\ij\lam_1^2 e^{\lam_1\del\eta} + \ddd_\ij\lam_2^2 e^{\lam_2\del\eta},
\end{align}
where $\ccc_\ij$ and $\ddd_\ij$ are amplitude coefficients.
For $\del\eta\ll1$, 
assume $T_\ij(\eta)\approx T_\ij(\eta + \del\eta)$,
%and $h_\ij'(\eta)\approx h_\ij''(\eta)\approx0$
%i.e., the strain derivatives vanish at the start of the time step.
%AR: changed
and $\lambda_{1,2}' \approx \lambda_{1,2}''
\approx 0$.
Then Eqn.~\ref{eqn:GW_mod_rewrite} becomes
\begin{equation}
\om^2h_\ij(\eta)\approx\ggg T_\ij(\eta)\approx \om^2h_\ij(\eta + \del\eta)\approx\ggg T_\ij(\eta + \del\eta),
\end{equation}
which gives
\begin{equation}
h_\ij(\eta)\approx h_\ij(\eta + \del\eta)\approx \om^{-2}\ggg T_\ij(\eta + \del\eta).
\end{equation}
This can be inserted into the ansatz and yield
\begin{align}
h_\ij(\eta + \del\eta) - \om^{-2}\ggg T_\ij(\eta + \del\eta) &
= \ccc_\ij e^{\lam_1\del\eta} + \ddd_\ij e^{\lam_2\del\eta} \\
h_\ij'(\eta + \del\eta) &
= \ccc_\ij\lam_1 e^{\lam_1\del\eta} + \ddd_\ij\lam_2 e^{\lam_2\del\eta} \\
h_\ij''(\eta + \del\eta) &
= \ccc_\ij\lam_1^2 e^{\lam_1\del\eta} + \ddd_\ij\lam_2^2 e^{\lam_2\del\eta}.
\end{align}
Equivalently in matrix form for $h_\ij$ and $h_\ij'$,
the ansatz is
\begin{equation}
\begin{pmatrix}
\om h - \om^{-1}\ggg T \\
h'
\end{pmatrix}_\ij^{\eta + \del\eta}
=
\begin{pmatrix}
\om e^{\lam_1\del\eta} & \om e^{\lam_2\del\eta} \\
\lam_1e^{\lam_1\del\eta} & \lam_2e^{\lam_2\del\eta}
\end{pmatrix}
\begin{pmatrix}
\ccc \\
\ddd
\end{pmatrix}_\ij.
\end{equation}
So the amplitude coefficients can be obtained as
\begin{align}
\begin{pmatrix}
\ccc \\
\ddd
\end{pmatrix}_\ij
& =
\begin{pmatrix}
\om e^{\lam_1\del\eta} & \om e^{\lam_2\del\eta} \\
\lam_1e^{\lam_1\del\eta} & \lam_2e^{\lam_2\del\eta}
\end{pmatrix}^{-1}
\begin{pmatrix}
\om h - \om^{-1}\ggg T \\
h'
\end{pmatrix}_\ij^{\eta + \del\eta}\\
& = 
\frac{1}{\om(\lam_2 - \lam_1)e^{(\lam_1 + \lam_2)\del\eta}}
\begin{pmatrix}
\lam_2e^{\lam_2\del\eta} & -\om e^{\lam_2\del\eta} \\
-\lam_1e^{\lam_1\del\eta} & \om e^{\lam_1\del\eta}
\end{pmatrix}
\begin{pmatrix}
\om h - \om^{-1}\ggg T \\
h'
\end{pmatrix}_\ij^{\eta + \del\eta}\\
\lim_{\del\eta\map0}
\begin{pmatrix}
\ccc \\
\ddd
\end{pmatrix}_\ij
& =
\frac{1}{\om(\lam_2 - \lam_1)}
\begin{pmatrix}
\lam_2 & -\om \\
-\lam_1 & \om
\end{pmatrix}
\begin{pmatrix}
\om h - \om^{-1}\ggg T \\
h'
\end{pmatrix}_\ij^{\eta},
\end{align}
where, in the last step,
we used the fact that $\ccc_\ij,\ddd_\ij$ are time-independent.
Therefore, the time evolution of the relevant quantities 
can be obtained via a matrix multiplication as
\begin{equation}
\begin{pmatrix}
\om h - \om^{-1}\ggg T \\
h'
\end{pmatrix}_\ij^{\eta + \del\eta}
=
\mmm
\begin{pmatrix}
\om h - \om^{-1}\ggg T \\
h'
\end{pmatrix}_\ij^{\eta},
\label{eqn:GW_evol_mtx_GAFD}
\end{equation}
where
\begin{equation}
\mmm =
\frac{1}{\lam_1 - \lam_2}
\begin{pmatrix}
\lam_1e^{\lam_2\del\eta} - \lam_2e^{\lam_1\del\eta} 
& \om(e^{\lam_1\del\eta} - e^{\lam_2\del\eta}) \\
\om^{-1}\lam_1\lam_2(e^{\lam_2\del\eta} - e^{\lam_1\del\eta})
& \lam_1e^{\lam_1\del\eta} - \lam_2e^{\lam_2\del\eta}
\end{pmatrix}.
\label{eqn:evol_mtx_full}
\end{equation}
In the case of $\sig = 0, \om^2 = k^2$,
the eigenvalues become $\lam_{1,2} = \pm ik$ 
and the evolution matrix reduces to
\begin{equation}
\mmm = 
\begin{pmatrix}
\cos{k\del\eta} & \sin{k\del\eta}\\
-\sin{k\del\eta} & \cos{k\del\eta}
\end{pmatrix}.
\label{eqn:evol_mtx_simp}
\end{equation}
Note that so far these are written in line with the convention of Ref.~\cite{RoperPol+18}.
However, in the {\sc Pencil Code}~\cite{pencil},
the relevant evolution equation (\Eq{eqn:GW_evol_mtx_GAFD}) is implemented as
\begin{equation}
\begin{pmatrix}
h - \om^{-2}\ggg T \\
\om^{-1}h'
\end{pmatrix}_\ij^{\eta + \del\eta}
=
\mmm
\begin{pmatrix}
h - \om^{-2}\ggg T \\
\om^{-1}h'
\end{pmatrix}_\ij^{\eta},
\label{eqn:GW_evol_mtx_pencil}
\end{equation}
where $\mmm$ is the same as in \Eq{eqn:evol_mtx_full}.

\section{WKB formulation}

Modifications of the GW propagation,
i.e., when $T_\ij = 0$,
can be studied using the WKB formulation~\cite{Nishizawa17,Arai+17}
\begin{equation}
h_\ij = \aaa e^{i\bbb}h_\ij^\GR,
\end{equation}
where $h_\ij^\GR$ is the standard strain in GR.
The total strain $h$ reads
\begin{equation}
h = e^{-\ddd}e^{-ik\Del T}h^\GR,
\end{equation}
where $h_\ij = he_\ij$ and
\begin{align}
\ddd &
= \frac{1}{2}\int_{\eta_*}^{\eta_0}\alpM\hhh d\eta'
= \frac{1}{2}\int_0^{z_*}\frac{\alpM}{1 + z'}dz',\\
\Del T &
= \int_{\eta_*}^{\eta_0}\Big((1 - \cT) - \frac{a^2\mg^2}{2k^2}\Big)d\eta'
= \int_0^{z_*}\frac{1}{\hhh}\Big(\frac{1-\cT}{1 + z'} - \frac{\mg^2}{2k^2(1 + z')^3}\Big)dz'.
\end{align}
\boxbegin
The following change of variables is used in writing the two integrals above in terms of $z$:
\begin{equation}
\int\hhh d\eta'
= \int\frac{1}{a}\frac{da}{d\eta'}d\eta'
= \int(1+z)da
= \int(1+z)d\Big(\frac{1}{1+z}\Big)
= \int(1+z)\frac{-dz}{(1+z)^2}
= \int \frac{-dz}{1 + z}.
\end{equation}
\boxend
Here $\eta_*$ is the time of GW production,
and $z$ is redshift such that $a = 1/(1 + z)$.
To determine the integrands and integration bounds,
we look at the relevant quantities
during RD, MD, and $\Lam$D.

\red{
WKB approximation for the energy spectrum can be obtained by considering
\begin{equation}
\EGW (k) \sim h' {h'}^*.
\end{equation}
The derivatives are
\begin{equation}
 h' = - \ddd' e^{-\ddd}e^{-ik \Delta T} h^\GR - ik \Delta T'
 e^{-\ddd} e^{-i k \Delta T} h^\GR + e^{-\ddd} e^{-ik\Delta T} 
 {h'}^\GR,
\end{equation}
where
\begin{equation}
    \ddd' = \frac{1}{2} \alpM \HH, \quad \Delta T' = 
    \bigl(1 - \cT\bigr) - \frac{a^2 \mg^2}{2 k^2}.
\end{equation}
Hence,
\begin{align}
    h' = &\, - \frac{1}{2} \alpM \HH e^{-\ddd}e^{-ik \Delta T} h^\GR - ik \biggl(\bigl(1 - \cT\bigr) - \frac{a^2 \mg^2}{2 k^2}\biggr)
    e^{-\ddd} e^{-i k \Delta T} h^\GR + e^{-\ddd} e^{-ik\Delta T} 
    {h'}^\GR \\
    {h'}^* = &\, - \frac{1}{2} \alpM \HH e^{-\ddd}e^{ik \Delta T} {h^\GR}^* + ik \biggl(\bigl(1 - \cT\bigr) - \frac{a^2 \mg^2}{2 k^2}\biggr)
    e^{-\ddd} e^{i k \Delta T} {h^\GR}^* + e^{-\ddd} e^{ik\Delta T} 
    {{h'}^\GR}^*,
\end{align}
and we can write
\begin{align}
    \EGW (k) \sim  &\, \frac{1}{4} \alpM^2 \HH^2 e^{-2\ddd} h^\GR
    {h^\GR}^* - \frac{1}{2} \alpM \HH e^{-2\ddd} \bigl(
    h^\GR {h'^\GR}^* + {h'^\GR}  {h^\GR}^* \bigr)
    + k^2 \biggl(\bigl(1 - \cT\bigr) - \frac{a^2 \mg^2}{2k^2}
    \biggr)^2 e^{-2\ddd} h^\GR {h^\GR}^* \nonumber \\ &\, + ik \biggl(\bigl(1 - \cT\bigr)
    - \frac{a^2 \mg^2}{2k^2} \biggr) e^{-2\ddd} \bigl({h^\GR}^* h'^\GR  -  h^\GR {h'^\GR}^*\bigr) + e^{-2\ddd} h'^\GR {h'^\GR}^*.
\end{align}
We now assume that ${h'}^\GR \sim k h^\GR$ to get
\begin{align}
    \EGW (k) \sim  &\, \biggl[ \frac{1}{4 k^2} \alpM^2 \HH^2 e^{-2\ddd} - \frac{1}{k} \alpM \HH e^{-2\ddd}
    + \biggl(\bigl(1 - \cT\bigr) - \frac{a^2 \mg^2}{2k^2}
    \biggr)^2 e^{-2\ddd} + e^{-2\ddd} \biggr] h'^\GR {h'^\GR}^*.
\end{align}
Therefore, the GW spectrum in modified GR due to the parameters
$\alpM$ and $\cT$ is
\begin{equation}
    \EGW (k) = \biggl[ \frac{1}{4 k^2} \alpM^2 \HH^2 -
    \frac{1}{k} \alpM \HH + \biggl(\bigl(1 - \cT\bigr) -
    \frac{a^2 \mg^2}{2k^2}
    \biggr)^2  + 1 \biggr] \EGW^\GR (k)
    e^{-2\ddd}.
\label{eqn:EGW_mod}
\end{equation}
Omitting $\cT = \mg = 0$, we can observe that the GW spectrum can
be exactly zero if the following condition holds
\begin{equation}
    \frac{\alpM^2}{k^2} - \frac{4}{\HH} \frac{\alpM}{k} +
    \frac{4}{\HH^2} = 0 \Rightarrow \frac{\alpM}{k} =
    \frac{2}{\HH} \pm \sqrt{\frac{4}{\HH^2} - \frac{4}{\HH^2}} =
    \frac{2}{\HH} \Rightarrow k_0 = \frac{\alpM}{2\HH}.
\end{equation}
}
In the case $\mg = 0$,
\begin{equation}
\frac{\alpM}{k} = \frac{2\pm2\sqrt{-(1 - \cT)^2}}{\hhh}
\Map
k_0 = \frac{\alpM\hhh}{2}.
\end{equation}

\boxbegin
Summary of WKB formulation:
\begin{itemize}
\item IR modification of $k^{-2}$ upon the GR solution is to be expected.
\item Such IR modification depends on $\alpM$, but not $\alpT$.
\item For $\alpM > 0$, there exists critical wavenumber at $k_0 = \alpM\hhh/2$ where $\EGW$ vanishes.
\blue{
\item The same $k_0$ also marks where $k^{-2}$ modification kicks in.
\item As time increases, $\hhh$ decreases, so $k_0$ also decreases.
\item Eventually, when time becomes too large, the $k^{-2}$ slope modification becomes too small to be shown on whichever $k$ range we start with.}
\end{itemize}
\begin{figure}[H]
\centering
\includegraphics[scale = 0.55]{Figures/WKB_time_evol.png}
\caption{Time evolution of a mock energy spectrum using WKB approximation,
where the gray curve indicates the initial spectrum in both panels.
Left panel: $\alpT = -0.99$ (solid) and $\alpT = 99$ (dashed).
Right panel: $\alpM = -0.1$ (solid) and $\alpM = 1$ (dashed).}
\label{fig:WKB_time_evol}
\end{figure}
\boxend

\section{Stability conditions}

For a function $\pi$ of the form
\begin{equation}
A(\eta)\ddot\pi + B(\eta)\dot\pi + \Big(C(\eta) + k^2 D(\eta)\Big)\pi + E(\eta,k) = 0,
\end{equation}
we require~\cite{Hu+13,Creminelli+08}
\begin{itemize}
\item $A > 0$ to avoid ghost instability. We have $A = 1$ so this is automatically satisfied.
\item $\cs^2\equiv D/A\geq0$ to avoid gradient instability. In addition, we require $\cs^2\leq1$ for Lorentz invariance.
\blue{
Here $\cs^2 = \cT^2/a^2$. So we need
$-1\leq\alpT\leq a^{-2} - 1$.
This naturally suggests a parameterization choice for $\alpT$,
i.e., the maximum modification is allowed by setting 
$\alpT(a) = a^{-2} - 1$.
}
\item $m_\pi^2\equiv C/A\geq0$ to avoid tachyonic instability~\cite{Gumrukcuoglu+13}.
\blue{
Recall that the stability conditions should be applied to the original equation \Eq{eqn:GW_mod} and not the scaled \Eq{eqn:GW_conf_comove}.
So the non-negativity of graviton mass,
$m_\pi^2 = \mg^2\geq0$,
satisfies tachyonic stability.
}
\end{itemize}

\section{Parameterizations of $\alpM$ and $\alpT$}

For $\alpT$, unless $\alpT$ takes on unrealistic values (see Fig.~\ref{fig:WKB_time_evol}),
it does not change the spectral shapes or amplitudes.
Therefore, we will not further consider time- or frequency-dependent~\cite{deRham+18} parameterizations of $\alpT$.
The present day value of the tensor speed excess $\alpTz$ is tightly constrained.
So any modifications in this parameter should only take place in the early universe.
This behavior could be parameterized by
\begin{equation}
\alpT(\eta) = \alpTz\tanh\Big(1 - \Big(\frac{a(\eta)}{a_0}\Big)^\tau\Big),
\end{equation}
where $\tau > 0$ is a constant.
This ensure that at late times, $\alpT$ vanishes to recover GR.
And $\alpT$ approaches some nonzero constant at early times.

Instead, we focus on $\alpM$ and its various parameterization choices.

The modified parameters can be written down in accordance with the specific modified gravity models~\cite{Odintsov+22}.
But here we parameterize phenomenologically.

The simplest consideration would be $\alpM = \alpMz = $ constant.

For the Planck-mass run rate,
the modifications could be either early- or late-time.
If late-time modifications are considered,
then the parameterization could be
\begin{equation}
\alpM(\eta) = 
\begin{cases}
\cM\Om_\Lam(\eta)
& \text{(choice 1)},\\
\alpMz a^n(\eta)
& \text{(choice 2)}.
\end{cases}
\end{equation}
Here in choice 1, $\cM$ is a constant
and $\Om_\Lam$ is the dark energy density.
This ensure that the modifications are negligible in the early universe,
but become more important in the late universe.
For choice 2,
\begin{equation}
\begin{dcases}
\alpMz < 0 \quad\text{for}\quad n > \frac{5}{2},\\
\alpMz > 0 \quad\text{for}\quad 0 < n < 1 + \frac{3}{2}\Om_{\mat,0}.
\end{dcases}
\end{equation}

In choice 2, we have (for $n \neq 0, -1$)
\begin{equation}
    {\cal D} = \frac{1}{2} \alpM^{0} \int_{z_0}^{z_*}
    \frac{1}{(1 + z)^{1 + n}} \dd z = - \frac{1}{2n}
    \alpM^{0} \Bigl[ \bigl(1 + z_*\bigr)^{-n - 2}
    - \bigl(1 + z_0\bigr)^{-n - 2} \Bigr].
\end{equation}
Hence, this reads
\begin{equation}
    {\cal D} (z_*) = \frac{1}{2n} \alpM^0 \Biggl[1 -
    \biggl(\frac{a_0}{a_*}\biggr)^{-n - 2} \, \Biggr]
    \approx \frac{1}{2n} \alpM^0,
\end{equation}
where the approximation holds whenever $a_0/a_* \gg 1$ and $n > -1$.

\subsubsection{Parameterize $\alpM(\eta) = \alpMz a^n(\eta)$}

\begin{align}
\int_{\eta_*}^{\eta_0}(\alpM\hhh)'d\eta 
& = \int_{\eta_*}^{\eta_\EQ}(\alpM'\hhh + \alpM\hhh')d\eta
+ \int_{\eta_\EQ}^{\eta_\Lam}(\alpM'\hhh + \alpM\hhh')d\eta
+ \int_{\eta_\Lam}^{\eta_0}(\alpM'\hhh + \alpM\hhh')d\eta\\
& = \alpMz\Bigg[
\ccc_\RD^n\eta^{n-1}\Big|_{\eta_*}^{\eta_\EQ}
+ 2\ccc_\MD^n\eta^{2n-1}\Big|_{\eta_\EQ}^{\eta_\Lam}
- \ccc_\LD^n\eta^{-n-1}\Big|_{\eta_\Lam}^{\eta_0}
\Bigg]
\end{align}

\subsubsection{Parameterize $\alpM(\eta) = \alpMz\Om_\Lam(\eta)$}

\begin{align}
\int_{\eta_*}^{\eta_0}(\alpM\hhh)'d\eta 
& = \int_{\eta_*}^{\eta_\EQ}(\alpM'\hhh + \alpM\hhh')d\eta
+ \int_{\eta_\EQ}^{\eta_\Lam}(\alpM'\hhh + \alpM\hhh')d\eta
+ \int_{\eta_\Lam}^{\eta_0}(\alpM'\hhh + \alpM\hhh')d\eta\\
& = \alpMz\rho_\Lam\Bigg[
\int_{\eta_*}^{\eta_\EQ}
\Big(\frac{2\rho_{\mat,0}a_0^3}{\ccc_\RD^3\eta^5} + \frac{3\rho_{\rad,0}a_0^4}{\ccc_\RD^4\eta^6} - \frac{\rho_\Lam}{\eta^2}\Big)
\Big/
\Big(\rho_\Lam + \rho_{\mat,0}\Big(\frac{a_0}{\ccc_\RD\eta}\Big)^3 + \rho_{\rad,0}\Big(\frac{a_0}{\ccc_\RD\eta}\Big)^4\Big)^2d\eta\\\nonumber
& + \int_{\eta_\EQ}^{\eta_\Lam}
\Big(\frac{10\rho_{\mat,0}a_0^3}{\ccc_\MD^3\eta^8} + \frac{14\rho_{\rad,0}a_0^4}{\ccc_\MD^4\eta^{10}} - \frac{2\rho_\Lam}{\eta^2}\Big)
\Big/
\Big(\rho_\Lam + \rho_{\mat,0}\Big(\frac{a_0}{\ccc_\RD\eta}\Big)^3 + \rho_{\rad,0}\Big(\frac{a_0}{\ccc_\RD\eta}\Big)^4\Big)^2d\eta\\\nonumber
& + \int_{\eta_\Lam}^{\eta_0}
\Big(\frac{4\rho_{\mat,0}a_0^3}{\ccc_\LD^3\eta^{-1}} + \frac{5\rho_{\rad,0}a_0^4}{\ccc_\LD^4\eta^{-2}} + \frac{\rho_\Lam}{\eta^2}\Big)
\Big/
\Big(\rho_\Lam + \rho_{\mat,0}\Big(\frac{a_0}{\ccc_\LD\eta^{-1}}\Big)^3 + \rho_{\rad,0}\Big(\frac{a_0}{\ccc_\LD\eta^{-1}}\Big)^4\Big)^2d\eta
\Bigg]
\end{align}

\subsection{Frequency-dependent $\alpT$}

\begin{equation}
\cT^2(k) = (1 + \alpT) = \Big(1 + \sum_n\beta_n\Big(\frac{k}{k_*}\Big)^n\Big)^2
\end{equation}

\begin{equation}
\alpT(k) = \Big(1 + \sum_n\beta_n\Big(\frac{k}{k_*}\Big)^n\Big)^2 - 1
\end{equation}

\section{Run parameters}

\section{Constraints}

\begin{table}[H]
\centering
\begin{tabular}{l|ll}
         & Lower bound        & Upper bound \\\hline
$\alpMz$ & $-0.35$            & $0.35$ \\
$\alpTz$ & $-6\times10^{-15}$ & $1.4\times10^{-15}$
\end{tabular}
\caption{Current constraints on the modified parameters.}
\label{tab:constraints}
\end{table}

\section{Simulations: $1024^3$ no source}

\begin{figure}[H]
\centering
\includegraphics[scale = 0.55]{Figures/EGW_1024.png}
\caption{$\EGW(k)$ for different values of $\alpM$ and $\alpT$.
Black curves indicate GR solution. 
Dashed curves indicate WKB approximation for the corresponding $\alpM$ or $\alpT$.}
\label{fig:EGW1024}
\end{figure}

\blue{From Fig.~\ref{fig:EGW1024},
we see that the IR regime for $\alpM\neq0$ is the most interesting.
Because here the deviation from GR is the most significant,
and such deviations are not captured by WKB.
On the other hand,
$\alpT\neq0$ cases do not introduce much changes to the GR spectrum.
So in general we probably shouldn't focus too much on $\alpT$.
Maybe no need to consider time or frequency dependent $\alpT$ at all?}

\section{Simulations: $2048^1$ no source}

\begin{figure}[H]
\centering
\includegraphics[scale = 0.47]{Figures/sp_alpT_const.png}
\caption{$\alpT\in\{-0.9, -0.1, 0.1, 0.9\}$.}
\label{fig:sp_alpT_const}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale = 0.47]{Figures/sp_alpM_const.png}
\caption{$\alpM\in\{-1, -0.5, -0.1, 0.1\}$.}
\label{fig:sp_alpM_const}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale = 0.47]{Figures/sp_alpM_scale_factor.png}
\caption{$\alpM = \alpMz(a/a_*)^1$ with $\alpMz\in\{-0.1, -0.01, 0.01, 0.1\}$.}
\label{fig:sp_alpM_scale_factor}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale = 0.47]{Figures/sp_alpM_dark_energy.png}
\caption{$\alpM = \alpMz\Om_\Lam$ with $\alpMz\in\{-1,-0.1,0.1,1\}$.}
\label{fig:sp_alpM_dark_energy}
\end{figure}

\section*{Appendix (bin..)}

\section{Analytical solution 1}

Let's start considering the sourceless solution with the modifications introduced by $\alpM$, i.e., setting $\alpT = \mg = 0$.
Then, \Eq{modGR} yields
\begin{equation}
    \hij'' + \alpM \HH \hij' + \biggl(k^2 - \alpM \HH^2 - \frac{a''}{a}
    \biggr) \hij = 0.
\end{equation}
We use the ansatz $\hij = A \hij^\GR e^{-B}$, such that the time derivatives are
\begin{align}
    \hij' = &\, A {\hij'}^\GR e^{-B} - A B' \hij^\GR e^{-B} + A'
    \hij^\GR e^{-B}, \\
    \hij'' = &\, A {\hij''}^\GR e^{-B} - 2 A B' {\hij'}^\GR e^{-B}
    + A \bigl({B'}^2 - B''\bigr) \hij^\GR e^{-B} + 2 A'{\hij'}^\GR
    e^{-B} - 2 A'B' \hij^\GR e^{-B} \nonumber \\
    + &\, A'' \hij^\GR e^{-B},
\end{align}
where $\hij^\GR$ is the solution to the GW equation with $\alpM = 0$.
Hence, the GW equation can be written as
\begin{align}
    A {\hij''}^\GR e^{-B} &\, -  2 A B' {\hij'}^\GR e^{-B} + A \bigl({B'}^2 - B''\bigr) \hij^\GR e^{-B} + 2 A'{\hij'}^\GR
    e^{-B} - 2 A'B' \hij^\GR e^{-B} + A'' \hij^\GR e^{-B} \nonumber \\
    &\, + \alpM \HH \Bigl(A {\hij'}^\GR e^{-B} - A B' \hij^\GR e^{-B}
    + A'\hij^\GR e^{-B} \Bigr) + \biggl(k^2 - \alpM \HH^2 - \frac{a''}{a}
    \biggr) A \hij^\GR e^{-B} = 0, \\
     - 2 A B' {\hij'}^\GR &\, + A \bigl({B'}^2 - B''\bigr) \hij^\GR
    + 2 A'{\hij'}^\GR - 2 A'B' \hij^\GR + A'' \hij^\GR
    \nonumber \\
    &\, + \alpM \HH \Bigl(A {\hij'}^\GR - A B' \hij^\GR + A' \hij^\GR \Bigr) - \alpM \HH^2 A
    \hij^\GR = 0, \\
    {\hij'}^\GR \Bigl[-2 A B' + &\,  2 A' + \alpM \HH A \Bigr] + \hij^\GR
    \Bigl[A {B'}^2 - A B'' - 2 A'B' + A'' - \alpM \HH A B' + \alpM \HH A'
    - \alpM \HH^2 A \Bigr] = 0.
\end{align}
This yields two equations (for $A = 1$)
\begin{equation}
    B' = \frac{1}{2} \alpM \HH, \quad
    B'' - {B'}^2 + B' \alpM \HH + \alpM \HH^2 = 0,
    \label{Aeq1}
\end{equation}
which does not have a compatible solution, and for generic $A$,
\begin{gather}
    -2 B' + 2 \frac{A'}{A} + \alpM \HH = 0, \\
    B'^2 - B'' - 2 \frac{A'}{A} B' + \frac{A''}{A} - \alpM \HH
    B' + \alpM \HH \frac{A'}{A} - \alpM \HH^2 = 0.
\end{gather}
From the first equation, we get $B' = A'/A + \half \alpM \HH$, which,
substituting in the second equation, leads to
{\scriptsize
\begin{equation}
    \biggl(\frac{A'}{A}\biggr)^2 + \frac{1}{4} \alpM^2 \HH^2 
    + \alpM \HH \frac{A'}{A} - \frac{A''}{A} - \biggl(\frac{{A'}}{A}\biggr)^2
    - \frac{1}{2} \alpM' \HH - \frac{1}{2} \alpM \HH' -
    2 \biggl(\frac{A'}{A} \biggr)^2 - \alpM \HH \frac{A'}{A} +
    \frac{A''}{A} - \alpM \HH \frac{A'}{A} - \frac{1}{2} \alpM^2
    \HH^2 + \alpM\HH \frac{A'}{A} - \alpM \HH^2 = 0.
\end{equation}}
\begin{gather}
    \frac{1}{4} \alpM^2 \HH^2 + 2 \frac{{A'}^2}{A^2} + \frac{1}{2}
    \alpM' \HH + \frac{1}{2} \alpM \HH' + \alpM \HH^2 = 0, \nonumber \\
    \frac{A'}{A} = \frac{i}{2} \sqrt{\frac{1}{2} \alpM^2 \HH^2 +
    \alpM \HH' + \alpM' \HH + 2 \alpM\HH^2}.
\end{gather}
Since $\HH' + \HH^2 = a''/a$, we can rewrite the solution as
\begin{equation}
    \frac{A'}{A} = \frac{i}{2} \sqrt{\frac{1}{2} \alpM^2 \HH^2 +
    \alpM\HH^2 + 
    \alpM \frac{a''}{a} + \alpM' \HH}.
\end{equation}
We can now assume that $\alpM$ is constant and that we are in the RD epoch,
such that $a'' = 0$ to get rid of the last two terms:
\begin{equation}
    \frac{A'}{A} = \frac{i}{2} \HH \sqrt{\frac{1}{2} \alpM^2 + \alpM}
    \Rightarrow A = e^{\frac{i}{2}\sqrt{\half \alpM^2 + \alpM}
    \int^\tau \HH \dd \tau'}.
\end{equation}
Otherwise, the general solution is
\begin{equation}
    A = e^{\frac{i}{2}\int^\tau
    \sqrt{\half \alpM^2 \HH^2 + \alpM \HH^2 +  \alpM \frac{a''}{a} +
    \alpM' \HH} \dd \tau'}.
\end{equation}
Going back to $B$, we substitute now the solution of $A$ to get
\begin{equation}
    B = \frac{i}{2} \int^\tau \sqrt{\frac{1}{2} \alpM^2 \HH^2 + \alpM \HH^2 +
    \alpM\frac{a''}{a} + \alpM' \HH} \dd \tau' + \frac{1}{2}
    \int^\tau \alpM \HH \dd \tau'.
\end{equation}
\red{
The general solution is then
\begin{equation}
    \hij = A \hij^\GR e^{-B} = e^{\int A'/A \dd \tau} e^{-\int A'/A \dd \tau - \half\int \alpM\HH \dd \tau} \hij^\GR = e^{-\half \int \alpM \HH \dd \tau} \hij^\GR = e^{-B} \hij^\GR,
\end{equation}
which holds for any $A$ such that $A'/A$ is not a function of A, as it
is in this case!!
Hence, the solution should be equivalent to setting $A = 1$, but this leads to an absurd [see \Eq{Aeq1}]!!
}

\blue{
\Eq{Aeq1} is the result of checking against the comoving GW equation (with scaled strain),
whereas we should check against the unscaled equation.
See the next section for this,
which also doesn't seem to be right.
This is probably because we took the ansatz different from the one in the literature.
Correcting this leads to the same results as in the papers -- see the next-next section below.
However, this doesn't explain issue \#3.}

\section{Analytical solution 2}

Use ansatz (tensor indices suppressed)
\begin{align}
h & = Ah^\GR e^{-B},\\
h' & = \Big(Ah'^\GR - AB'h^\GR + A'h^\GR\Big)e^{-B},\\
h'' & = \Big(Ah''^\GR - 2 AB'h'^\GR
+ A(B'^2 - B'')h^\GR + 2 A'h'^\GR - 2 A'B'h^\GR + A''h^\GR\Big)e^{-B}.
\end{align}
Substitute back into \Eq{eqn:GW_conf} (not \Eq{eqn:GW_conf_comove}):
\begin{align}
0 & 
= Ah''^\GR - 2 AB'h'^\GR
+ A(B'^2 - B'')h^\GR + 2 A'h'^\GR - 2 A'B'h^\GR + A''h^\GR\\\nonumber
& + (2 + \alpM)\hhh\Big(Ah'^\GR - AB'h^\GR + A'h^\GR\Big) + k^2Ah^\GR\\
& = h'^\GR(-2AB' + 2A')
+ h^\GR\Big[A(B'^2 - B'') - 2A'B' + A'' + (2 + \alpM)\hhh(-AB' + A')\Big].
\end{align}
This leads to
\begin{align}
0 & = -AB' + A',\\
0 & = A(B'^2 - B'') - 2A'B' + A'' + (2 + \alpM)\hhh(-AB' + A').
\end{align}

\section{Analytical solution 3}

Ansatz in Ref.~\cite{Nishizawa17}:
\begin{align}
h & = Ae^{iB},\\
h' & = (A' + iAB')e^{iB},\\
h'' & = (A'' + i2A'B' + iAB'' - AB'^2)e^{iB}.
\end{align}
Substitute back into \Eq{eqn:GW_conf} (not \Eq{eqn:GW_conf_comove}):
\begin{align}
0 & = A'' + i2A'B' + iAB'' - AB'^2 + (2 + \alpM)\hhh(A' + iAB') + k^2A\\
& = \frac{A''}{A} + i2B'\frac{A'}{A} + iB'' - B'^2 + (2 + \alpM)\hhh\Big(\frac{A'}{A} + iB'\Big) + k^2,
\end{align}
giving
\begin{equation}
0 = \frac{A''}{A} - B'^2 + (2 + \alpM)\hhh\frac{A'}{A} + k^2
\quad,\quad
0 = 2\frac{A'}{A} + \frac{B''}{B'} + (2 + \alpM)\hhh.
\end{equation}
We can get the expression for $B'$,
\begin{align}
    B' = &\, \sqrt{\frac{A''}{A} + (2 + \alpM) \HH \frac{A'}{A} + k^2}, \\
    \frac{B''}{B'} =  \frac{1}{2} \frac{{({B'}^2})'}{{B'}^2} = 
    \frac{1}{2} &\, \Biggl(\frac{A'''}{A} - \frac{A'' A'}{A^2} + \bigl(2 + \alpM
    \bigr) \HH \biggl[\frac{A''}{A} - \frac{{A'}^2}{A^2} \biggr] + 
    \alpM' \HH \frac{A'}{A} + \bigl(2 + \alpM\bigr) \HH' \frac{A'}{A} \Biggr)
    \nonumber \\
    \times &\,
    \Biggl(\frac{A''}{A} + \bigl(2 + \alpM \bigr) \HH \frac{A'}{A} + k^2
    \Biggr)^{-1}.
\end{align}
The equation becomes
\begin{gather}
    \Biggl(2 \frac{A'}{A} + \bigl(2 + \alpM\bigr)\HH \Biggr)
    \Biggl(\frac{A''}{A} + \bigl(2 + \alpM \bigr) \HH \frac{A'}{A} + k^2
    \Biggr) +  \nonumber \\ +
    \frac{1}{2} \Biggl(\frac{A'''}{A} - \frac{A'' A'}{A^2}
    + \bigl(2 + \alpM
    \bigr) \HH \biggl[\frac{A''}{A} - \frac{{A'}^2}{A^2} \biggr] + 
    \alpM' \HH \frac{A'}{A} + \bigl(2 + \alpM\bigr) \HH' \frac{A'}{A} \Biggr)
    \\
    \frac{A'''}{A} + 3 \frac{A'' A'}{A^2} + 3 \bigl(2 + \alpM\bigr) \HH \frac{A''}{A} +
    3 \bigl(2 + \alpM\bigr) \HH  \frac{{A'}^2}{A^2} +
    \frac{A'}{A} \Bigl[4 k^2 + 2 \bigl(2 + \alpM\bigr)^2 \HH^2
    + \alpM' \HH + \bigl(2 + \alpM\bigr)
    \HH' \Bigr]
     \nonumber \\
    + 2 \bigl(2 + \alpM\bigr) \HH k^2= 0.
\end{gather}

\boxbegin
\blue{
\st{
Issue 1:
Since we adopt $a\map\bar a = a/a_*$ 
and $\eta\map\bar\eta = \eta/\eta_*$,
shouldn't we also consider
$z\map\bar z = z/z_*$ such that $z_0 = 0$ and $z_* = 1$?
Because if we have $\bar z = a_0/a - 1$ above, which is the same as $\bar a = a_0/(1 + z)$,
then we're essentially having a different scaling/redefinition for the scale factor as $a\map\bar a = a a_0\neq a/a_*$, right?
}
}
\boxend

\boxbegin
\blue{
Issue 2: For $\alpM =$ constant and $\alpT = 0$, 
the modifications are $h = e^{-\ddd}h^\GR$.
Given Eqn.~\ref{eqn:D_in_a}, we have
\begin{equation}
h = e^{-\frac{\alpM}{2}\ln\frac{a_0}{a_*}}h^\GR
= \Big(\frac{a_0}{a_*}\Big)^{-\alpM/2}h^\GR
= a^{-\alpM/2}h^\GR,
\end{equation}
where in the last step
we have used $a_* = 1$ 
and set $a_0 = a$ as a scale factor within RD.
This means
\begin{align}
h' & = a^{-\alpM/2}\Big(-\frac{\alpM}{2a}h^\GR + h'^\GR\Big),
\label{eqn:hprime_hGR}\\
h'' & = a^{-\alpM/2}\Big(\frac{\alpM^2 + 2\alpM}{4a^2}h^\GR - \frac{\alpM}{a}h'^\GR + h''^\GR\Big).
\label{eqn:hpprime_hGR}
\end{align}
Recall GW propagation during RD:
\begin{equation}
h''^\GR + k^2 h^\GR = 0
\quad\text{and}\quad
h'' + \frac{\alpM}{a}h' + (k^2 - \frac{\alpM}{a^2})h = 0,
\end{equation}
the latter of which could be rewritten by using Eqns.~\ref{eqn:hprime_hGR} and~\ref{eqn:hpprime_hGR}
\begin{align}
0 & = \Big(\frac{\alpM^2 + 2\alpM}{4a^2}h^\GR - \frac{\alpM}{a}h'^\GR + h''^\GR\Big)
+ \frac{\alpM}{a}\Big(-\frac{\alpM}{2a}h^\GR + h'^\GR\Big)
+ \Big(k^2 - \frac{\alpM}{a^2}\Big)h^\GR\\
0 & = h''^\GR + \Big(\frac{\alpM^2 + 2\alpM}{4a^2} - \frac{\alpM^2}{2a^2} + k^2 - \frac{\alpM}{a^2}\Big)h^\GR\\
0 & = h''^\GR + \Big(k^2 - \frac{\alpM^2 + 2\alpM}{4a^2}\Big)h^\GR \\
0 & = \Big(\frac{\alpM^2 + 2\alpM}{4a^2}\Big)h^\GR.
\end{align}
Does this imply that either $\alpM = 0$ or $h^\GR = 0$?
This doesn't seem reasonable.
Where did I do wrong?
}
\boxend

\boxbegin
\blue{
Issue 3: Testing a 0d sourceless case at arbitrary $k$, with an initial values of $h_\ini = 1, h'_\ini = 0$.
Result is below.
\begin{figure}[H]
\centering
\includegraphics[scale=0.7]{Figures/0d_ts_alpM.png}
\caption{
\blue{
Time evolution of GW strain $h(k = 0.5)$ in GR (red),
in modified gravity with $\alpM = -0.5$ solved in Python (blue solid)
and IDL (green solid).
The blue dashed curve is $e^{-\ddd}h^\GR$,
i.e., the red curve multiplied by a factor of $e^{-\ddd}$.
Here $\ddd$ is obtained with \Eq{eqn:D_in_a}.
However, the blue dashed curve does not overlap
the blue or green solid curves
although they should.
}}
\end{figure}
}
\boxend

\section{One-dimensional test cases for constant $\alpM$}

Here we show 1d test cases for a few constant values of $\alpM$,
where the initial strain spectrum is flat ($h_0(k) = 1$ for all $k$)
and the source is turned off ($T = 0$).

\begin{figure}[H]
\centering
\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{Figures/1d_ts_alpM_pos.png}
\caption{Time evolution of $h(k)$ where $k = 0.2$ and $5.0$ for $\alpM\in\{0.1,0.3,0.5,0.9\}$.}
\end{subfigure}
\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{Figures/1d_ts_hspec_alpM_pos.png}
\caption{Time evolution of total strain $h(\eta)$ (left), and strain spectrum $h(k)$ at $\eta = 7.0$ (right).}
\end{subfigure}
\end{figure}

\begin{figure}[H]
\centering
\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{Figures/1d_ts_alpM_neg.png}
\caption{Time evolution of $h(k)$ where $k = 0.2$ and $5.0$ for $\alpM\in\{-0.1,-0.3,-0.5,-0.9\}$.}
\end{subfigure}
\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{Figures/1d_ts_hspec_alpM_neg.png}
\caption{Time evolution of total strain $h(\eta)$ (left), and strain spectrum $h(k)$ at $\eta = 7.0$ (right).}
\end{subfigure}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale = 0.8]{Figures/a_eta.png}
\caption{Scale factor as a function of conformal time.}
\label{fig:a_eta}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale = 0.55]{Figures/rho_om_eta.png}
\caption{Evolution of the fractional energy contents.}
\label{fig:rho_om_eta}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale=0.6]{idl/p.pdf}
\end{figure}

\section{Field redefinitions in different epochs}

Table~\ref{tab:fac_evol_era} shows the scaling relations of relevant quantities during the three cosmic epochs.
But to quantify the exact expressions,
we start by writing
\begin{equation}
a_\RD(\eta) = \ccc_\RD\eta,\quad
a_\MD(\eta) = \ccc_\MD\eta^2,\quad
a_\LD(\eta) = \ccc_\LD\eta^{-1}.
\end{equation}
\begin{table}[H]
\centering
\begin{tabular}{l|l|l|l|l|l}
Epoch   &$a(t)$             & $a(\eta)$
& $\hhh\equiv a'/a$         & $a''/a$ 
& $z(t)$ \\\hline
RD      &$\propto t^{1/2}$  & $\propto\eta^1$    
& $=\eta^{-1}$              & $=0$ 
& $\propto t^{-1/2}$ \\
MD      &$\propto t^{2/3}$  & $\propto\eta^2$    
& $=2\eta^{-1}$             & $=2\eta^{-2}$ 
& $\propto t^{-2/3}$ \\
$\Lam$D &$\propto e^{H_0t}$ & $\propto\eta^{-1}$
& $=-\eta^{-1}$             & $=2\eta^{-2}$
& $\propto e^{-H_0t}$
\end{tabular}
\caption{Relevant quantities in different cosmic epochs.}
\label{tab:fac_evol_era}
\end{table}

\subsection{During radiation domination (RD)}

Let $*$ be the subscript indicating the onset of RD,
then we have $\eta_* = \hhh_*^{-1}$ and $a_* = a_\RD(\eta_*) = \ccc_\RD\hhh_*^{-1}$.
This fixes the coefficient
\begin{equation}
\ccc_\RD = a_*\hhh_*,
\end{equation}
which in turn fixes the expression for $a_\RD$ to be
\begin{equation}
a_\RD(\bar\eta) = a_*\bar\eta\quad\forall\quad\bar\eta_*\leq\bar\eta\leq\bar\eta_\EQ,
\end{equation}
where $\bar\eta\equiv\hhh_*\eta$ is the normalized conformal time,
and $\bar\eta_\EQ$ marks the radiation-matter equality.

Furthermore, we set $a_* = 1$ as the onset of RD.
So $a' = 1, a'' = 0$ and $\hhh_* = 1$,
and \Eq{modGR2} reduces to
\begin{equation}
h_\ij'' + 
\frac{\alpM}{\bar\eta}h_\ij' + 
\Big((1 + \alpT) k^2 - \frac{\alpM}{\bar\eta^2}\Big)h_\ij
= \frac{6}{\bar\eta}T_\ij.
\label{eqn:GW_mod_rad}
\end{equation}

\subsection{During matter domination (MD)}

At radiation-matter equality, we have
$a_\EQ = a_\MD = \ccc_\MD\eta_\EQ^2 = a_\RD = \ccc_\RD\eta_\EQ = a_*\hhh_*\eta_\EQ$.
This gives the coefficient
\begin{equation}
\ccc_\MD = \frac{a_*\hhh_*}{\eta_\EQ}
\end{equation}
and the scale factor
\begin{equation}
a_\MD(\bar\eta) = \frac{a_*\hhh_*}{\eta_\EQ}\eta^2 = \frac{a_*}{\bar\eta_\EQ}\bar\eta^2
\quad\forall\quad\bar\eta_\EQ\leq\bar\eta\leq\bar\eta_\Lam,
\end{equation}
where $\bar\eta_\Lam$ is the onset of $\LD$ and
\begin{equation}
\bar\eta_\EQ = \hhh_*\eta_\EQ = \frac{a_\EQ}{a_*}.
\end{equation}
Numerically, take the redshift at radiation-matter equality as $z_\EQ = 3600$,
then the corresponding scale factor and conformal time can be obtained as
$a_\EQ/a_0 = 1/(1 + z_\EQ) = 0.0002777$ and $\bar\eta_\EQ = 0.0002777a_0$.
In MD, the GW equation is
\begin{equation}
h_\ij'' + 
\frac{2\alpM}{\bar\eta}h_\ij' + 
\Big((1 + \alpT) k^2 - \frac{4\alpM}{\bar\eta^2} - \frac{2}{\bar\eta^2}\Big)h_\ij
= \frac{6\bar\eta_\EQ}{\bar\eta^2}T_\ij.
\label{eqn:GW_mod_mat}
\end{equation}

\subsection{During dark energy domination ($\LD$)}

At the onset of $\LD$, we have $a_\Lam = a_\MD = \ccc_\MD\eta_\Lam^2 = a_\LD = \ccc_\LD\eta_\Lam^{-1}$. 
So the coefficient is
\begin{equation}
\ccc_\LD = \frac{a_*\hhh_*\eta_\Lam^3}{\eta_\EQ}
\end{equation}
and the expressions for $a_\LD$ becomes
\begin{align}
a_\LD = \frac{a_*\bar\eta_\Lam^3}{\bar\eta_\EQ}\bar\eta^{-1}
\quad\forall\quad\bar\eta_0\leq\bar\eta\leq\bar\eta_\Lam,
\end{align}
where
\begin{equation}
\bar\eta_\Lam = \sqrt{\frac{a_\Lam}{a_*}\bar\eta_\EQ}.
\end{equation}
Here take the redshift $z_\Lam = 0.4$,
then we have $a_\Lam/a_0 = 1/(1 + z_\Lam) = 0.7143$
and $\bar\eta_\Lam = 0.01408a_0$.
In $\LD$, the GW equation is
\begin{equation}
h_\ij'' - 
\frac{\alpM}{\bar\eta}h_\ij' + 
\Big((1 + \alpT) k^2 - \frac{\alpM}{\bar\eta^2} - \frac{2}{\bar\eta^2}\Big)h_\ij
= \frac{6\bar\eta_\EQ\bar\eta}{\bar\eta_\Lam^3}T_\ij.
\label{eqn:GW_mod_lam}
\end{equation}

\subsection{Numerical values of $\bar\eta_\EQ, \bar\eta_\Lam$ and $\bar\eta_0$}

So the critical points in time determining the three epochs depend on $a_0$,
which can be obtained via the conservation of $g(T)a^3T^3$ assuming adiabatic expansion during RD.
So $g_0a_0^3T_0^3 = g_*a_*^3T_*^3$, giving
\begin{equation}
\frac{a_0}{a_*} = \frac{T_*}{T_0}\Big(\frac{g_*}{g_0}\Big)^{1/3}.
\end{equation}
Today we have $T_0\sim2.33\times10^{-13}\GeV$ and $g_0\sim3.3$.
So assuming EWPT at $a_*$, then $T_*\sim100\GeV$ and $g_*\sim106$,
giving $a_0\sim1.3643\times10^{15}$.
If we assume QCDPT at $a_*$, then $T_*\sim0.15\GeV$ and $g_*\sim15$,
which give $a_0\sim1.0664\times10^{12}$.
Using these values of $a_0$, the critical points in time $\bar\eta_\EQ, \bar\eta_\Lam$ and $\bar\eta_0$ can be obtained in the table below.
\begin{table}[H]
\centering
\begin{tabular}{l|l|l|l|l}
& $\bar\eta_*$ 
& $\bar\eta_\EQ$    
& $\bar\eta_\Lam$  
& $\bar\eta_0$ \\\hline
EWPT  
& 1
& $3.789\times10^{11}$ 
& $1.9215\times10^{13}$ 
& $1.3725\times10^{13}$ \\
QCDPT 
& 1
& $2.961\times10^{8}$  
& $1.5019\times10^{10}$ 
& $1.0728\times10^{10}$
\end{tabular}
\end{table}

\section{Evolution of different energy contents}

We would like to obtain the evolution of dark energy (possibly also matter) through the epochs,
as they can be useful for parameterizing $\alpM$ later.
The energy density of the contents can be written as
\begin{equation}
\rho_\Lam = \frac{\Lam}{8\pi\GN},\quad
\rho_\mat(a) = \rho_{\mat,0}\Big(\frac{a_0}{a}\Big)^3,\quad
\rho_\rad(a) = \rho_{\rad,0}\Big(\frac{a_0}{a}\Big)^4.
\end{equation}
In conformal time, these are
\begin{equation}
\rho_\Lam = \frac{\Lam}{8\pi\GN},\quad
\rho_\mat(\bar\eta) =
\begin{dcases}
\rho_{\mat,0}\Big(\frac{a_0}{\ccc_\RD\bar\eta}\Big)^3\\
\rho_{\mat,0}\Big(\frac{a_0}{\ccc_\MD\bar\eta^2}\Big)^3\\
\rho_{\mat,0}\Big(\frac{a_0}{\ccc_\LD\bar\eta^{-1}}\Big)^3
\end{dcases},\quad
\rho_\rad(\bar\eta) = 
\begin{dcases}
\rho_{\rad,0}\Big(\frac{a_0}{\ccc_\RD\bar\eta}\Big)^4\\
\rho_{\rad,0}\Big(\frac{a_0}{\ccc_\MD\bar\eta^2}\Big)^4\\
\rho_{\rad,0}\Big(\frac{a_0}{\ccc_\LD\bar\eta^{-1}}\Big)^4
\end{dcases}.
\end{equation}
This is related to the dimensionless (fractional) density as
\begin{equation}
\Om_i = \frac{\rho_i}{\rho_\crit}\quad(i\in\{\Lam,\mat,\rad\}).
\end{equation}
Take $\Om_{\Lam,0} = 0.6889, \Om_{\mat,0} = 0.3111, \Om_{\rad,0} = 9.28656\times10^{-5}$~\cite{Planck18} as the current day values.
The critical density is calculated as $\rho_\crit = 3H^2/(8\pi\GN)$.
So by taking $H_0 = 2.1928\times10^{-18}\rm s^{-1}$,
we have $\rho_{\crit,0} = 9.4135\times10^{-27}\rm kg/m^3$.
\begin{equation}
\rho_{\Lam,0} = \rho_\Lam = 6.4849\times10^{-27}\rm kg/m^3,\quad
\rho_{\mat,0} = 2.9285\times10^{-27}\rm kg/m^3,\quad
\rho_{\rad,0} = 8.7419\times10^{-31}\rm kg/m^3.
\end{equation}
Using these, we can obtain the evolution history of $\Om_i$ as
\begin{equation}
\Om_\Lam(\bar\eta) = \frac{\rho_\Lam}{\rho_\crit(\bar\eta)},\quad
\Om_\mat(\bar\eta) = \frac{\rho_\mat(\bar\eta)}{\rho_\crit(\bar\eta)},\quad
\Om_\rad(\bar\eta) = \frac{\rho_\rad(\bar\eta)}{\rho_\crit(\bar\eta)},
\end{equation}
where
\begin{equation}
\rho_\crit(a) = \rho_\Lam + \rho_\mat + \rho_\rad
= \rho_\Lam + \rho_{\mat,0}\Big(\frac{a_0}{a}\Big)^3 + \rho_{\rad,0}\Big(\frac{a_0}{a}\Big)^4,
\end{equation}
or in conformal time,
\begin{equation}
\rho_\crit(\bar\eta) =
\begin{dcases}
\rho_\Lam + \rho_{\mat,0}\Big(\frac{a_0}{\ccc_\RD\bar\eta}\Big)^3 + \rho_{\rad,0}\Big(\frac{a_0}{\ccc_\RD\bar\eta}\Big)^4 
&\forall\quad\bar\eta_*\leq\bar\eta\leq\bar\eta_\EQ\\
\rho_\Lam + \rho_{\mat,0}\Big(\frac{a_0}{\ccc_\MD\bar\eta^2}\Big)^3 + \rho_{\rad,0}\Big(\frac{a_0}{\ccc_\MD\bar\eta^2}\Big)^4 
&\forall\quad\bar\eta_\EQ\leq\bar\eta\leq\bar\eta_\Lam\\
\rho_\Lam + \rho_{\mat,0}\Big(\frac{a_0}{\ccc_\LD\bar\eta^{-1}}\Big)^3 + \rho_{\rad,0}\Big(\frac{a_0}{\ccc_\LD\bar\eta^{-1}}\Big)^4 
&\forall\quad\bar\eta_0\leq\bar\eta\leq\bar\eta_\Lam
\end{dcases}
\end{equation}

\subsection{In terms of $\eta$}

Therefore, the factor $\ddd$ and $\Del T$ should be calculated as a sum of three pieces
\begin{equation}
\ddd(\eta) = \frac{1}{2}\Big(
\int_\RD\frac{\alpM}{\eta'}d\eta'
+ \int_\MD\frac{2\alpM}{\eta'}d\eta'
- \int_\LD\frac{\alpM}{\eta'}d\eta'\Big),
\end{equation}
where the integration limits are the conformal time in the corresponding epoch.
Given that $d\eta = a^{-1}dt$,
these can be calculated explicitly as
\begin{equation}
\eta(t) - 1 =
\begin{dcases}
\int_{t_*}^{t}\ccc_\RD^{-1} t'^{-1/2}dt'
= 2\ccc_\RD^{-1}
\big(t^{1/2} - t_*^{1/2}\big)
& \quad\forall\quad\RD:\quad t_*\leq t\leq t_\EQ,\\
\eta(t_\EQ) + 
\int_{t_\EQ}^{t}\ccc_\MD^{-1} t'^{-2/3}dt'
= \eta(t_\EQ) + 3\ccc_\MD^{-1}
\big(t^{1/3} - t_\EQ^{1/3}\big)
& \quad\forall\quad\MD:\quad t_\EQ\leq t\leq t_\Lam,\\
\eta(t_\Lam) + 
\int_{t_\Lam}^{t}\ccc_\LD^{-1}e^{-H_0t'}dt'
= \eta(t_\Lam) + \ccc_\LD^{-1}H_0^{-1}
\big(e^{-H_0t_\Lam} - e^{-H_0t}\big)
& \quad\forall\quad\LD:\quad t_\Lam\leq t\leq t_0,
\end{dcases}
\end{equation}
where $t_*$ is the time of GW generation approximately at $10^{-10}$ second for EWPT and $10^{-5}$ second for QCDPT,
$t_\EQ\approx4.7\times10^4$ years (or $z_\EQ\approx3600$)
is the time of matter-radiaiton equality,
$t_\Lam\approx9.8\times10^9$ years (or $z_\Lam\approx0.4$)
is when the universe started accelerated expansion.
The $\eta(t) - 1$ on the LHS is to enforce $\eta_*(t_*) = 1$ (otherwise it would be 0).
The coefficient $\ccc_\LD$ is such that $a_0 = \ccc_\LD e^{H_0 t_0} = 1$.
And the coefficients $\ccc_\MD$ and $\ccc_\RD$ are obtained by matching the boundary conditions at $t_\Lam$ and $t_\EQ$,
i.e.,
$a(t_\Lam) = \ccc_\LD e^{H_0t_\Lam} = \ccc_\MD t_\Lam^{2/3}$
and $a(t_\EQ) = \ccc_\MD t_\EQ^{2/3} = \ccc_\RD t_\EQ^{1/2}$.
Explicitly, these constants take the values
\begin{equation}
\ccc_\LD = e^{-H_0t_0},\quad
\ccc_\MD = e^{-H_0(t_0 -t_\Lambda)}
t_\Lambda^{-2/3},\quad
\ccc_\RD = e^{-H_0(t_0 -t_\Lambda)}
t_\Lambda^{-2/3}t_\EQ^{1/6}.
\end{equation}
\boxbegin
\textbf{Exact expressions of scale factor as functions of conformal time.}
\begin{equation}
a(\eta) =
\begin{dcases}
\ccc_\RD\Big(\frac{1}{2}\ccc_\RD\eta + t_*^{1/2}\Big)
&\forall\quad\RD:\quad\eta_*\leq\eta\leq\eta_\EQ,\\
\ccc_\MD\Big(\frac{1}{3}\ccc_\MD(\eta - \eta_\EQ) + t_\EQ^{1/3}\Big)^2 
&\forall\quad\MD:\quad\eta_\EQ\leq\eta\leq\eta_\Lam,\\
\ccc_\LD\Big(e^{-H_0 t_\Lam} - \ccc_\LD H_0(\eta - \eta_\Lam)\Big)^{-1}
&\forall\quad\LD:\quad\eta_\Lam\leq\eta\leq\eta_0,
\end{dcases}
\end{equation}
\boxend
After obtaining $a\in[a_*, a_0]$ and $\eta\in[\eta_*, \eta_0]$
in the convention $a_0 = 1$ and $\eta_*\neq1$,
we normalize the scale factor against its initial value $a_*$
such that $\bar a\in[\bar a_*, \bar a_0]$,
where $\bar a_* = a_*/a_* = 1$ and $\bar a_0 = a_0/a_* = a_*^{-1}$.
Similarly we normalize the conformal time such that
$\bar\eta_* = \eta_*/\eta_* = 1$ and $\bar\eta_0 = \eta_0/\eta_*$.
Then we drop the overbars.

\boxbegin
At $\eta_\EQ$,
$a'_\RD = \ccc_\RD\neq a'_\MD = 2\ccc_\MD\eta_\EQ$
and $a''_\RD = 0\neq a''_\MD = 2\ccc_\MD$.

At $\eta_\Lam$,
$a'_\MD = 2\ccc_\MD\eta_\Lam\neq a'_\LD = -\ccc_\LD\eta_\Lam^{-2}$
and $a''_\MD = 2\ccc_\MD = a''_\LD = 2\ccc_\LD\eta_\Lam^{-3}$.

So in the previous section, only $a$ and $a''$ at $\eta_\Lam$ are continuous,
whereas other quantities such as
$a'$ at both $\eta_\EQ$ and $\eta_\Lam$,
and $a''$ at $\eta_\EQ$ are discontinuous.
\boxend

\subsection{Constant $\alpM$ (in terms of $z$)}

\begin{equation}
\ddd(z_*)
= \frac{\alpM}{2}\Big(\ln(1 + z_*) - \ln(1 + z_0)\Big),
\end{equation}
Since $z = a_0/a - 1$,
then in our convention, $z_* = a_0/a_* - 1$
and $z_0 = 0$.
Hence, we can write
\begin{equation}
%    \ddd (z_*) = \frac{\alpM}{2} \ln \frac{a_0}{a_*} = 
    \ddd (a_*) = \frac{\alpM}{2} \ln \frac{a_0}{a_*} = 
    34.76511 \, \frac{\alpM}{2} \biggl(\ln \frac{T_*}{100 \GeV}
    + \frac{1}{3} \ln \frac{g_*}{100} \biggr)
\label{eqn:D_in_a}
\end{equation}

\subsection{Constant $\alpT$}

In this case, the factor $\Delta T$ is
\begin{align}
    \Delta T = &\, \bigl(1 - \sqrt{1 + \alpT} \bigr)
    \int_{z_0}^{z_*} \frac{1}{\HH} \frac{1}{1 + z'} \dd z' \nonumber \\
    = &\, 
    \bigl(1 - \sqrt{1 + \alpT} \bigr)
    \Biggl[- a_0^{-1} \, \ccc_\LD \int_{z_0}^{z_\Lam}
    \dd z' + a_0^{1/2} \, \ccc_\MD \int_{z_\Lam}^{z_\EQ} \frac{1}{(1 + z')^{3/2}} \dd z' + a_0 \, \ccc_\RD
    \int^{z_*}_{z_\EQ} \frac{1}{(1 + z')^2} \dd z'
    \biggr] \nonumber \\
    = &\, \bigl(1 - \sqrt{1 + \alpT} \bigr)
    \biggl(- a_0^{-1} \, \ccc_\LD \, z_\Lam - 2 \, a_0^{1/2} \,
    \ccc_\MD
    \Bigl[ \bigl(1 + z_\EQ\bigr)^{-1/2} - \bigl(1 + z_\Lam\bigr)^{-1/2}\Bigr] - a_0\, \ccc_\RD \biggl[
    \frac{1}{1 + z_*} - \frac{1}{1 + z_\EQ}\biggr]
    \biggr) \nonumber \\
    = &\, \bigl(1 - \sqrt{1 + \alpT} \bigr)
    \biggl(- a_0^{-1} \, \ccc_\LD z_\Lam - 2 \, a_0^{1/2} \, \ccc_\MD
    \Biggl[\sqrt{\frac{a_\EQ}{a_0}} - 
    \sqrt{\frac{a_\Lambda}{a_0}} \, \Biggr] - a_0\, \ccc_\RD
    \biggl[\frac{a_* - a_\EQ}{a_0} \biggr] \biggr),
\end{align}
where we have used $z_0 = 0$ and the constants $\ccc$ are defined such that
$\eta_\MD = \ccc_\MD a^{1/2}$, $\eta_\RD = \ccc_\RD a$,
$\eta_\LD = - \ccc_\LD a^{-1}$.
Their values are $\ccc_\RD = (a_* \HH_*)^{-1}$, $\ccc_\MD =
a_\EQ^{1/2} \, (a_* \HH_*)^{-1}$, and $\ccc_\LD = - \sqrt{a_\Lambda^3 a_\EQ} \, (a_* \HH_*)^{-1}$.
Hence,
\begin{equation}
    \Delta T = \bigl(1 - \sqrt{1 + \alpT} \bigr) \, \HH_*^{-1}
    \Biggl( \frac{\sqrt{a_\Lam a_\EQ}}
    {a_*} \, \frac{a_0 - a_\Lam}{a_0} - 2 \,
    \frac{a_\EQ}{a_*} + 2 \frac{a_\Lam}{a_0} -
    \frac{a_* - a_\EQ}{a_*}  \Biggr).
\end{equation}
We can set the normalization $a_0 = 1$, such that
$a_\EQ = 0.0002777$ and $a_\Lam = 0.7142857$ and leave $a_*/a_0$
as a parameter that will allow us to consider different times
of GW production within the RD epoch.
\begin{align}
    \Delta T = &\, (1 - \sqrt{1 + \alpT}) \, \HH_*^{-1}
    \Biggl(0.004024 \, \frac{a_0}{a_*} - 1 - 0.0002777 \,
    \frac{a_0}{a_*} + 1.42857 \Biggr) \nonumber \\
    = &\, (1 - \sqrt{1 + \alpT}) \, \HH_*^{-1}
    \Biggl(0.0037462\, \frac{a_0}{a_*} + 0.42857 \Biggr).
\end{align}
For an energy scale within the RD epoch (e.g., at the EWPT), we can compute $\HH_*$
and $a_*$,
\begin{equation}
    \HH_*^{-1} = 6.0697 \cdot 10^4 \s \, \biggl(\frac{100 \GeV}{T_*}
    \biggr) \biggl(\frac{g_*}{100}\biggr)^{1/6}, \quad
    \frac{a_0}{a_*} = 1.254 \cdot 10^{15} \biggl(\frac{T_*}{100 \GeV}
    \biggr)\biggl(\frac{g_*}{100}\biggr)^{1/3},
\end{equation}
where we have assumed that the entropic and relativistic dof are
the same at the time of generation.
Finally, we get
\begin{equation}
    \Delta T = (1 - \sqrt{1 + \alpT}) \biggl(\frac{g_*}{100}\biggr)^{1/2} 
    \Biggl[
    2.851384 \cdot 10^{17} 
    + 2.6012913 \cdot 10^4 \biggl(\frac{T_*}{100 \GeV}
    \biggr)\biggl(\frac{g_*}{100}\biggr)^{-1/6} \, \Biggr]
\end{equation}

\boxbegin
Given $\eta_\MD = \ccc_\MD a^{1/2}, \eta_\RD = \ccc_\RD a, \eta_\LD = -\ccc_\LD a^{-1}$ above,
write
\begin{equation}
a(\eta) = 
\begin{dcases}
\ccc_\RD^{-1}\eta & 
\forall\quad\RD:\quad\eta_*\leq\eta\leq\eta_\EQ\\
\ccc_\MD^{-2}\eta^2 & 
\forall\quad\MD:\quad\eta_\EQ\leq\eta\leq\eta_\Lam\\
-\ccc_\LD\eta^{-1} &
\forall\quad\LD:\quad\eta_0\leq\eta\leq\eta_\Lam
\end{dcases}
\end{equation}
\red{
Check continuity of $a$ at $\eta_\EQ$ and $\eta_\Lam$.
\begin{align}
a(\eta_\EQ) &
= \ccc_\RD^{-1}\eta_\EQ 
= \hhh_*a_*\eta_\EQ, \quad 
a(\eta_\EQ)
= \ccc_\MD^{-2}\eta_\EQ^2 
= a_\EQ^{-1}(\hhh_*a_*)^2\eta_\EQ^2 \Rightarrow
a_\EQ = \HH_* a_* \eta_\EQ ,\\
a(\eta_\Lam) &
= \ccc_\MD^{-2}\eta_\Lam^2 
= a_\EQ^{-1}(\hhh_* a_*)^2\eta_\Lam^2, \quad 
a(\eta_\Lam)
= -\ccc_\LD\eta_\Lam^{-1} 
= \sqrt{a_\Lam^3 a_\EQ}(a_*\hhh_*)^{-1}\eta_\Lam^{-1} \Rightarrow
a_\Lam = a_\EQ^{-1} (a_* \HH_*)^2 \eta_\Lam^2.
\end{align}
}
\boxend

\begin{align}
\Del T(z_*) & = \int_{z_0}^{z_*}\frac{1}{\hhh}\Big(\frac{-\alpT}{1+z} - \frac{\mg^2}{2k^2(1+z)^3}\Big)dz\\
& = \int_{z_0}^{z_*}\frac{1}{aH}\Big(\frac{-\alpT}{1+z} - \frac{\mg^2}{2k^2(1+z)^3}\Big)dz\\
& = \int_{z_0}^{z_*}(1+z)H_0^{-1}\Big(\Om_{\Lam,0} + \Om_{\mat,0}(1+z)^3 + \Om_{\rad,0}(1+z)^4\Big)^{-1/2}\Big(\frac{-\alpT}{1+z} - \frac{\mg^2}{2k^2(1+z)^3}\Big)dz,
\end{align}
where we have used
\begin{equation}
\frac{H(z)}{H_0} = \sqrt{\Om_{\Lam,0} + \Om_{\mat,0}\Big(\frac{a_0}{a}\Big)^3 + \Om_{\rad,0}\Big(\frac{a_0}{a}\Big)^4}
= \sqrt{\Om_{\Lam,0} + \Om_{\mat,0}(1+z)^3 + \Om_{\rad,0}(1+z)^4}.
\end{equation}
Here $\Om_{\Lam,0} = 0.6889, \Om_{\mat,0} = 0.3111, \Om_{\rad,0} = 9.28656\times10^{-5}$.

\bibliographystyle{naturemag}
\bibliography{references}

\end{document}