\documentclass[a4paper,11pt]{article}
\pdfoutput=1
\usepackage{jcappub,macros,graphicx,float,mathtools,caption,subcaption} \graphicspath{{./Figures/}}
\usepackage[T1]{fontenc}

\title{\boldmath Modified propagation of gravitational waves from the early radiation era}

\author[a,b]{Yutong He,}
\author[c,d,e]{Alberto Roper Pol}
\author[a,b,e,f]{and Axel Brandenburg}

\affiliation[a]{Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfv\'ens v\"ag 12, 10691 Stockholm, Sweden}
\affiliation[b]{The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden}
\affiliation[c]{D\'epartement de Physique Th\'eorique, Universit\'e de Gen\`eve, CH-1211 Gen\`eve, Switzerland}
\affiliation[d]{Universit\'e Paris Cit\'e, CNRS, Astroparticule et Cosmologie, Paris, F-75013, France}
\affiliation[e]{School of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili Ave, Tbilisi, GE-0194, Georgia}
\affiliation[f]{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213, USA
\\ \today}

\emailAdd{yutong.he@su.se}
\emailAdd{alberto.roperpol@unige.ch}
\emailAdd{brandenb@nordita.org}

\abstract{We study the propagation of cosmological gravitational wave
background (GWB) produced by hydromagnetic sources in the early radiation era until the present day.
Compared to standard general relativity (GR),
we study the effects that parameters of modified theories of gravity, such as 
the run rate of the effective Planck mass $\alpM$
and the tensor speed excess $\alpT$,
have on the present-day GW spectrum using numerical simulations beyond the WKB approximation.
While $\alpT$ makes relatively insignificant changes to the GR solution,
$\alpM$ can potentially introduce an enhancement to the spectral slopes of the energy spectrum in the low-frequency regime
depending on the chosen time evolution parameterization of $\alpM$.
This effect is additional to the damping or growth
occurring equally at all scales that can be predicted by the WKB approximation.
We discuss the observational implications in light of the recent observations by
pulsar timing array collaborations
and future detectors such as SKA, LISA, and DECIGO.}

\begin{document}
%\preprint{NORDITA-2022-080}
\maketitle
\flushbottom

\section{Introduction}
\label{sec:introduction}

The present day Hubble constant, $H_0$, is measured to be around $H_0\sim74\,{\rm km/s/Mpc}$ by astrophysical tests using type Ia supernovae~\cite{Riess:2019cxk,Riess:2020fzl,Riess:2021jrx}, lensed quasars~\cite{Wong:2019kwg}, and megamaser-hosting galaxies~\cite{Pesce:2020xfe}.
On the other hand, cosmological tests of $H_0$ constrain its value to $H_0\sim67\,{\rm km/s/Mpc}$ from
cosmic microwave background (CMB)~\cite{Planck:2018vyg} and baryon acoustic oscillations (BAO)~\cite{BOSS:2016wmc} experiments, assuming the standard $\Lam$CDM ($\Lam$ cold dark matter) model of cosmology.
These measurements suggest the presence of a $4$--$5\,\sig$ deviation
discrepancy between the early- and late-universe measurements of
$H_0$~\cite{Verde:2019ivm}, known as the Hubble tension.
Besides the well-known $H_0$ tension, a number of other observational discrepancies within the $\Lam$CDM model have been reported,
such as the $S_8$ tension, where $S_8\propto\sig_8\,\Om_{\mat,0}^{1/2}$ characterizes the structure growth rate $\sig_8$ and today's matter density $\Om_{\mat,0}$~\cite{Joudaki+17}.
A recent summary of over a dozen cosmological tensions and their varied
significance can be found in ref.~\cite{Perivolaropoulos+21}.

In view of these tensions, $\Lam$CDM, although being an extremely successful model, does not provide a
complete picture and it might be necessary to go beyond it.
Since general relativity (GR) is responsible for the gravity sector of $\Lam$CDM,
modifying $\Lam$CDM often requires modifying GR.
Within the broad subject of modified gravity, a theorem due to Lovelock proves the uniqueness of GR under four 
assumptions, and therefore provides four general approaches of modifying GR by relaxing each of them.
The 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, modified gravity could be obtained by adding additional fields
other than the metric, constructing 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.

Among the plethora of modified gravity mechanisms, we focus on scalar-tensor theories and the phenomenological effects in terms of the modified propagation of gravitational waves (GWs) therein.
The choice of scalar-tensor theories is due to their versatility.
In other words, apart from extending GR by the inclusion of additional fields, they can also be viewed as the four-dimensional effective description of certain higher-dimensional models after a suitable dimensional reduction, or higher-derivative models after reformulation.
Therefore, many modified gravity models following three of the four approaches motivated by circumventing Lovelock's theorem can be comprised into scalar-tensor theories.
In addition, generic four-dimensional scalar-tensor models with second-order derivatives can be conveniently summarized in the Horndeski class~\cite{Horndeski74,Deffayet:2011gz,Kobayashi:2011nu}.
As a result, Horndeski gravity includes a vast number of different models such as 
Brans-Dicke~\cite{BD61}, quintessence~\cite{Ratra:1987rm}, $f(R)$~\cite{DeFelice+10}, $f(G)$~\cite{Carroll:2004de}, $k$-essence~\cite{Armendariz-Picon+00}, kinetic gravity braiding~\cite{Deffayet:2010qz}, galileon~\cite{Nicolis:2008in}, etc.
Moreover, certain subclass within the Horndeski theories has been shown to ease the Hubble tension~\cite{Petronikolou:2021shp}.

At the level of linear perturbation, Horndeski gravity introduces two parameters to the standard GW
equation---%
the tensor speed excess $\alpT$ and the running of the effective Planck mass $\alpM$.
Although the multimessenger observation of the binary neutron star merger GW170817 and its gamma-ray burst GRB 170817A~\cite{LIGOScientific:2017zic} has put a tight constraint on the speed of GWs, i.e., $\alpT\lesssim\ooo(10^{-15})$, and therefore the theory space of Horndeski gravity~\cite{Ezquiaga:2017ekz},
a large number of models of Horndeski gravity remains viable.
%AR: Here we could add a comment about the fact that this test has been done at the ground-based detector
%AR: frequencies but it might be frequency-dependent, allowing larger values in the LISA or PTA bands.
%AR: I can add something after reviewing the rest of the paper!
%AB: yes, Yutong, you could add such a comment
%YH: perhaps it's better to leave as it is here in the intro? and later during the parameterization section we add a brief statement (which we have) that frequency-dependent alpT isn't ruled out
%AB: ok
%AR: I think a brief statement here would still be good even if we mention it in section 4
In fact, this shows the constraining power of GWs as a phenomenological probe.
Comparing to $\alpT$, the constraints on $\alpM$ are much less stringent.
Therefore, in the present work, we focus on exploring $\alpM$ and its effects on the energy spectrum of GWs.
This subclass of Horndeski theories is also known as reduced Horndeski theories.
Specifically, we initialize a GW spectrum based on previous studies of GWs sourced by primordial magnetohydrodynamic (MHD) turbulent
fields in the early radiation era---%
at the electroweak or QCD phase transitions (EWPT or QCDPT)---%
then propagate it through cosmic history to obtain the present-day relic spectrum.
This choice of the GW spectral shape is for convenience but we show that the obtained
results can be applied to predict the expected GW spectrum from other sources in reduced
Horndeski theories.
We then analytically and numerically compare the relic spectra in GR and Horndeski gravity, and discuss their potential
observational implications.
The numerical solutions in this study are obtained using the {\sc Pencil Code}~\cite{pencil},
which has been a tool for simulations of GWs from primordial turbulent sources since the implementation of a GW solver~\cite{RoperPol+18}.

Finally, note that an additional propagating degree of freedom (DOF) appears in Horndeski gravity as a consequence of the scalar field.
However, we restrict the current study to tensor-mode perturbations,
since the amplitude of the scalar mode remains subdominant~\cite{Takeda+21},
and the discussions of extra DOFs as a smoking gun for modified gravity overall might deserve a separate study.

We introduce the propagation of tensor-mode perturbations described by the GW equation in
modified theories of gravity in \Sec{mod_GWs}.
In \Sec{sec:WKB}, we introduce the WKB approximation, which is commonly used to describe approximate solutions
to the GW equation in modified GR.
Then, in \Sec{sec:param}, we present some common temporal parameterizations of the Horndeski parameter $\alpM$
through the cosmological history of the Universe that have been used in the literature,
and we present the numerical simulations that we perform to solve the GW equation using the {\sc Pencil Code},
along with the numerical solutions to the Friedmann equations in \Sec{sec:numerical_solutions}.
We study, in particular, the additional effect on the GW spectrum that is not found using the WKB
approximation.
Finally, we discuss potential observational implications of reduced Horndeski theories compared to GR
for the different parameterizations of $\alpM$ in \Sec{sec:obs_implications} and conclude in \Sec{sec:conclusions}.

Throughout the paper, we set $c = 1$, use the metric signature $(-+++)$, and define the gravitational coupling constant
$\kappa = 8 \pi G_{\rm N}$.
We indicate with a prime
derivatives with respect to conformal time normalized by the
conformal Hubble rate at the time of GW generation $\eta \HH_*$ 
and with a dot derivatives with respect to cosmic time $t$.
Both times are related via the
scale factor as $a \dd \eta = \dd t$.

\section{Modified GWs on FLRW background}
\label{mod_GWs}

\subsection{GW equation in Horndeski theories}

The homogeneous and isotropic background is described by the Friedmann-Lema\^itre-Robertson-Walker (FLRW) metric.
Including tensor perturbations, the FLRW line element reads
\begin{equation}
\dd s^2 = a^2 \bigl(\dd \eta^2 + [\delta_{ij} + a^{-1}
\hij (\xx, \eta)]\dd x^i\dd x^j\bigr),
\label{ds2}
\end{equation}
where $\hij = a\hij^\phys$ are
the strains obtained by scaling the physical strains $\hij^\phys$ with the scale factor $a$.
In GR, the GW equation in Fourier space\footnote{We use the Fourier convention
\begin{equation}
<<<<<<< .mine
    {\tilde h}_{\ij}'' (\kk, t) +
||||||| .r48
    {\tilde h}_\ij'' (\kk, t) +
=======
    \tilde h (\kk) = \int h (\xx) \, e^{-i \kk \cdot \xx} \dd^3 \xx,
    \quad
    h (\xx) = \frac{1}{(2\pi)^3}\int \tilde h (\kk) \, e^{i\kk \cdot \xx} \dd^3 \kk.
\end{equation}} 
reads \cite{Caprini:2019pxz,RoperPol+18}
\begin{equation}
    {\tilde h}_\ij'' (\kk, \eta) +
>>>>>>> .r103
    \biggl(k^2 - \frac{a''}{a} \biggr) \tilde h_{ij} 
    (\kk, \eta) = \frac{6}{a} \tilde T_{ij}^{\rm TT} (\kk, \eta),
    \label{GW_GR}
\end{equation}
where
$\HH \equiv a'/a$ is the
conformal Hubble rate normalized by $\HH_*$,
and an asterisk refers to a reference time in the
early universe, e.g., the time at which the GWs were
generated.
$\tilde T_{ij}^{\rm TT} = \Lambda_{ijlm} \tilde T_{lm}$
is the traceless-transverse (TT) projection of the
normalized stress energy tensor, i.e., divided by the
radiation energy density $\EErad = 3H^2/\kappa$, being $H\equiv\dot a/a$ the Hubble rate at time $\eta$.
The projection operator is
$\Lambda_{ijlm} (\hat \kk) = P_{il} P_{jm} - \half P_{ij} P_{lm}$,
$P_{ij} (\hat \kk) = \delta_{ij} - \hat k_i \hat k_j$, and
$\hat \kk = \kk/k$.
The scale factor at $\eta_*$ is set to unity,
and the wave numbers $k$ are also normalized by
the conformal Hubble rate.

In Horndeski gravity~\cite{Horndeski74},
the action consists of the metric tensor and an additional scalar field $\phi$.
Its linear perturbations on the FLRW metric
can be characterized by four free parameters
separable from the full theory~\cite{Bellini+14}.
These free parameters are all arbitrary functions of the scalar 
field and its kineticity $\nabla_\mu\phi\,\nabla^\mu\phi$.
They respectively indicate the following:
\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}
Out of the four parameters, however,
only $\alpM$ and $\alpT$ enter and modify the GW equation.
As a result, the modifications generalize \Eq{GW_GR} to~\cite{Nunes+18}
\begin{equation}
\tilde h_{ij}'' (\kk, \eta) + \alpM \HH \tilde h_{ij}'(\kk, \eta)
+ \biggl(\cT^2 k^2 - \alpM \HH^2 -
\frac{a''}{a}\biggr) \tilde h_{ij}(\kk, \eta) = \frac{6}{a} \tilde
T_{ij}^{\rm TT} (\kk, \eta),
\label{eqn:GW_mod}
\end{equation}
where
$\cT$ denotes the GW propagation speed such that $1 + \alpT \equiv \cT^2$,
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}
The graviton is massless in Horndeski gravity.
However, a massive graviton could be included by adding a $\mg \neq 0$ in \Eq{eqn:GW_mod}
as done in ref.~\cite{He:2021bqm} but we omit it in the present work.

\subsection{Friedmann equations}

In previous works using the {\sc Pencil Code} (see \Sec{sec:numerical_solutions} for more details about the 
numerical simulations), the simulations are
restricted to specific eras, e.g.,
radiation (RD) or matter domination (MD), using
a constant equation of state (EOS),
defined to be $w\equiv p/\rho$,
being $p$ the pressure and $\rho$ the energy density,
such that $w = 1/3$ and $0$ during RD and MD, respectively.
We also ignored the thermal history of the universe during RD,
represented by the relativistic $g_*$ and adiabatic $\gS$ DOFs.
In GR, the previous formulation is justified since the evolution
<<<<<<< .mine
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
||||||| .r48
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
=======
of the physical strains $h_\ij^\phys$ when the source is inactive ($T_\ij \sim 0$) can
be approximated to dilute as $h_\ij^{\rm phys}
\propto a^{-1}$ if one
neglects the evolving relativistic DOFs and
>>>>>>> .r103
transitions between radiation, matter, and dark energy 
dominations ($\Lam$D) \cite{Saikawa:2018rcs}.
In the present work, we focus on including such effects
and solve the GW equation
in modified gravity, which presents richer dynamics even when
the source is inactive, from the time of GW generation
up to present time.

However, assuming a piece-wise EOS,
such that $w = 1/3$, $0$, and $-1$ during RD, MD, and 
$\Lam$D, respectively,
would create discontinuities in the time evolution of $a(\eta)$ and its derivatives needed in \Eq{eqn:GW_mod}.
Therefore, 
to find a smooth $a(\eta)$,
we directly solve the Friedmann equations,
\begin{equation}
    \frac{\ddot{a}}{a} = - \frac{\kappa}{6} \, \rho \bigl(1 + 3w\bigr),
    \quad 
    \biggl(\frac{\dot a}{a} \biggr)^2 =
    \frac{\kappa}{3} \, \rho,
    \label{friedmann}
\end{equation}
where the dots denote derivatives with respect to cosmic time $t$.
A more convenient form of \Eq{friedmann} can be expressed in terms of the energy density $\Om(a)$ and a smooth $w(a)$ 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}
where $\Om(a)$ is defined to be the ratio of total energy density to
the present-day critical energy density
$\rho_{\crit, 0}\equiv3H_0^2/\kap$, i.e.,
\begin{align}
\Om(a) = \frac{\rho (a)}{\rho_{\crit, 0}} = &\,
    \Omega_\rad (a) + \Omega_\mat (a) + \Omega_{\Lam,0} \nonumber \\
    = &\,
     \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{align}
Numerically we take the present-time values of $\Om_{\Lam,0} \simeq 0.684$,
$\Om_{\rad,0} \simeq 4.16 \times10^{-5} \, h^{-2}$, and
$\Om_{\mat,0} = 1 - \Om_{\rad, 0} - \Om_{\Lam, 0}
\simeq 0.316$, where $h$ takes into account the uncertainties on
the present-time Hubble rate $H_0 = 100 \, h$ km/s/Mpc.
We set $h \simeq 0.67$ for the numerical
studies, using the value observed from the CMB,
and $g_*^0 \simeq 3.36$ and $g_{\rm S}^0 \simeq 3.91$ are 
the reference
relativistic and adiabatic DOFs at the present time\footnote{Note that 
neutrinos' contribution to the radiation energy density is accounted for by taking $g_*^0\simeq 3.36$, instead of $g_*^0 = 2$ at the present day due to photons only.
This leads to an excess in the calculation of the
radiation energy density after neutrinos become massive.
However, this occurs when the radiation energy is subdominant and
hence, it does not affect our calculations.}.

The evolution of the relativistic and adiabatic DOFs 
as a function of the temperature during RD are taken
from ref.~\cite{Hindmarsh:2020hop} and expressed as a function
of $a$ by taking $a^3 T^3 \gS$ to be constant, following an
adiabatic expansion of the universe.

The characteristic EOS, corresponding to
the energy density of \Eq{rho_a}, can be computed combining
equations~(\ref{friedmann2})
by taking the time derivative of the second equation and
introducing the first equation.
This yields
\begin{gather}
        \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)= \bigg({1 \over 3} \Om_{\rad} (a) 
    - \Om_{\Lam, 0}\bigg) \, \Om^{-1} (a).
    \label{eqn:w_a}
\end{equation}
\EEq{Om_a_dot} justifies
the expected evolution of $\Om(a)$ used in \Eq{rho_a}: approximately 
proportional to $a^{-4}$ during 
the RD era, to $a^{-3}$ during
the MD era, and
constant during $\Lambda$D.
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.
\EEq{eqn:w_a} yields
$w = 1/3$ for RD, $w = 0$ for MD, and
$w = -1$ for $\Lambda$D, as expected.
In the intermediate times, the function $w (a)$ transitions
smoothly.
With \Eqs{rho_a}{eqn:w_a} explicitly in terms of $a$, 
we can compute the relevant quantities that appear in the modified GW \Eq{eqn:GW_mod},
given in \Eq{friedmann2}, and expressed in terms of the
normalized conformal time $\eta \HH_*$ as
\begin{equation}
   \HH = \frac{a H_0}{\HH_*} \sqrt{\Omega (a)},
    \quad 
    \frac{a''}{a} = \frac{1}{2} \HH^2
    \bigl[1 - 3 w (a) \bigr],
    \label{eqn:a_derivatives_omega}
\end{equation}
where the term $\HH_*$ appears due to our definition of the
normalized $\HH = a'/a$, and we set $a_* = 1$ for consistency
with our GW equation; see discussion below \Eq{GW_GR}.
Finally, since \Eq{eqn:GW_mod} is expressed in terms of conformal time $\eta$ but \Eq{eqn:a_derivatives_omega} is still written in terms of $a$, we would like to substitute the variables via $a(\eta)$, which can be obtained using \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 $\dd t = a \dd \eta$, so we solve
\begin{equation}
    H_0 \bigl( \eta - \eta_{\rm ini} \bigr) =
    \int_{a_{\rm ini}}^a
    \frac{\dd a}{a^2 \sqrt{\Omega (a)}}.
    \label{dota}
\end{equation}
For the numerical integration, we set $a_{\rm ini} = 10^{-20}$ at
$t_{\rm ini} = \eta_{\rm ini} = 0$, which yields accurate results
for all $a \gtrsim 10^{-19}$.

In the present work, we consider that the GWs are generated
during a phase transition
(in particular, at the EWPT or QCDPT) within the RD era.
Assuming adiabatic expansion of the universe, one can compute $a_*/a_0$ as a function of the temperature $T_*$ and the adiabatic DOFs,
\begin{equation}
    \frac{a_*}{a_0} \simeq 7.97 \cdot 10^{-16}
    \frac{T_*}{100 \, \GeV} \biggl(\frac{\gS}{100}\biggr)^{-{1\over3}}.
    \label{as_a0}
\end{equation}
Setting $a_* = 1$, such that $\HH_* = H_*$,
this gives a value of the Hubble rate (valid during the RD era),
\begin{equation}
    \HH_* = 2.066 \cdot 10^{10} \Hz \biggl(\frac{T_*}{100 \GeV}
    \biggr)^2 \biggl(\frac{g_*}{100} \biggr)^{1\over2},
    \label{HHs}
\end{equation}
with $\eta_* = \HH_*^{-1}$.
These results allow us to use the solutions from Friedmann equations
and adapt them to compute the variables that appear in \Eq{eqn:GW_mod}, normalized to the specific epoch of GW generation.

\section{WKB approximation}

\subsection{Solution of the GW equation}
\label{sec:WKB}

Modifications of the GW propagation,
in the absence of sources\footnote{%
In the absence of sources, the GW equation [see \Eq{eqn:GW_mod}] does
not depend on the wave vector $\kk$ but only on its modulus (the wave number) $k$.
Hence, the solution can simply be expressed as a function
of $k$.}
(i.e., $T_{ij} = 0$),
have been studied using the WKB
approximation \cite{Nishizawa17,Arai+17}.
The WKB solution can be obtained using the following ansatz
\begin{equation}
    \tilde h_{ij} (k, \eta) = h (k, \eta) \, e_{ij} =
    A e_{ij} \, e^{iB},
    \label{ansatz_WKB}
\end{equation}
<<<<<<< .mine
where $a^{-1} h_{\ij}^\GR$ are the standard strains in GR, and are solutions
to the following equation
||||||| .r48
where $a^{-1} h_\ij^\GR$ are the standard strains in GR, and are solutions
to the following equation
=======
where $A$ and $B$ are generic coefficients, and $e_{ij}$ is the polarization tensor.\footnote{
In the absence of sources, the GW propagation of any
polarization mode is the same, so we can just call $h$
the amplitude of each mode.
If the produced GW signal is polarized, then for each mode we
need to impose the corresponding initial conditions.}
Substituting this into \Eq{eqn:GW_mod}, one gets
\begin{align}
    &\, 2 \frac{A'}{A} + \frac{B''}{B'} + \alpM \HH = 0 \Rightarrow A (\eta) = e^{-\ddd}, \label{eq1_WKB}\\
    &\, {B'}^2 = \frac{A''}{A} + \alpM \HH \frac{A'}{A}
    + \cT^2 k^2 - \alpM \HH^2 - \frac{a''}{a} \approx
    \cT^2 k^2 \Rightarrow B = \pm k (\eta \HH_* - \Del T),
    \label{eq2_WKB}
\end{align}
where, following the WKB approximation, we neglect the
terms $A''/A$, $A'/A$, $\HH^2$, and $a''/a$, when compared
to $\cT^2 k^2$.
In addition, we also neglect the term $B''/B'
\sim \cT'/\cT \sim 0$.
We have defined the damping
factor $\ddd$ and the time delay due to the effective
GW speed $\Del T$
that appear in theories of modified gravity \cite{Nishizawa17}
>>>>>>> .r103
\begin{equation}
\ddd (\eta) = \frac{1}{2}\int_{1}^{\eta \HH_*}\alpM \hhh
\dd \eta', \quad
\Delta T (\eta) = \int_{1}^{\eta\HH_*} (1 - \cT) \dd \eta'.
\label{eqn:ddd_DelT}
\end{equation}
Hence, the solution is
\begin{equation}
    h (k, \eta) = e^{-\ddd} e^{\pm i k (\eta \HH_* - \Delta T)} =
    e^{-\ddd \mp i k \Delta T} h^\GR (k, \eta),
    \label{h_part_sol}
\end{equation}
where $\ddd = \Delta T = 0$ in GR.
For the initial conditions\footnote{In the current work, we focus on the propagation of a GW background after it has already
been generated and reached a stationary solution, and the source
is no longer active.
Since the propagation only depends on $k$, the initial conditions
can be computed from the 3D fields in Fourier space, after
shell integration over directions $\hat \kk$, for each polarization mode.}
$h (k, \eta_*) = h_* (k)$ and $h' (k, \eta_*) = h_*' (k)$, we get
\begin{align}
    h (k, \eta) = &\, e^{-\ddd} \biggl[h_* (k) \cos k (\eta \HH_* - 1 - \Del T) + \frac{h_*' (k) +
    \half \alpM^* h_* (k)}
    {k \cT} \sin k (\eta \HH_* - 1 - \Del T) \biggr],
    \nonumber \\
    = &\, e^{-\ddd} \biggl[h_* (k) \cos k \, \tilde c_{\rm T} (\eta \HH_* - 1 ) + \frac{h_*' (k) +
    \half \alpM^* h_* (k)}
    {k \cT} \sin k \, \tilde c_{\rm T} (\eta \HH_* - 1) \biggr],
    \label{h_modGR}
\end{align}
where $\alpM^* = \alpM (\eta_*)$ and $\tilde c_{\rm T}
(\eta) = \int \cT \dd \eta'/(\eta \HH_* - 1)$.
Note that the specific choice of initial conditions
does not allow to give a linear relation between $h$ and
$h^\GR$ as done in ref.~\cite{Nishizawa17} and, in general,
$h \neq e^{-\ddd - i k \Delta T}  h^\GR$, which is only true
when referring to the particular solution of the ODE, given
in \Eq{h_part_sol}.
This reduces to the solution in GR when $\cT = 1$ and $\alpM = 0$,
\begin{equation}
     h^{\rm GR} (k, \eta) = h_* (k) \cos k (\eta \HH_* - 1) + \frac{h_*' (k)}{k} \sin k (\eta \HH_* - 1).
    \label{GR_WKB}
\end{equation}
In reduced Horndeski theories with $\alpT = 0$,
the WKB solution yields
\begin{equation}
\tilde h_{ij} (k, \eta) = e^{-\ddd} \biggl[\tilde 
h_{ij}^\GR (k, \eta) + \half e_{ij} \frac{\alpM^* h_* (k)}{k} \sin k (\eta \HH_* - 1) \biggr].
\label{eqn:h_WKB}
\end{equation}
When $\alpM = 0$ and $\cT$ is a constant in time we find
\begin{equation}
    \tilde h_{ij} (k, \eta) = h_* (k) \cos \bigl[\cT k \bigl(
    \eta \HH_* - 1 \bigr)\bigr] + \frac{h_*'(k)}{k \cT}
    \sin \bigl[\cT k \bigl(\eta\HH_* - 1\bigr)\bigr].
    \label{sol_alpT}
\end{equation}

\subsection{Limitations of the WKB approximation}
\label{limitations_WKB}

To obtain the WKB solution we have neglected some terms in \Eq{eq2_WKB} that would otherwise appear in the solution.
These terms can be expressed in the following way
\begin{equation}
    \frac{A''}{A} + \alpM \HH \frac{A'}{A} - \alpM \HH^2
    - \frac{a''}{a}
    = - \half \alpM \Bigl(1 + \half \alpM \Bigr) \HH^2
    - \half \alpM' \HH - \frac{a''}{a} \Bigl(1 + \half \alpM\Bigr),
\end{equation}
where we have used $\HH' = a''/a - \HH^2$.
Hence, the WKB assumption could break down when at least one
of these terms is not negligible when compared to $k^2 \cT^2$.
We have one term that depends on $\alpM$, one term that
depends on the time evolution of $\alpM$, and one term
that depends on $a''/a$.
The latter appears also in GR when WKB is used to approximate the
solution for the scaled strains, which do not decay (note that
the decay of the physical strains is already absorbed by the
scale factor), as shown in \Eq{GR_WKB}.
This term is given in \Eq{eqn:a_derivatives_omega} and its upper
bound can be found using Friedmann equations (using $h=0.67$),
\begin{equation}
    \frac{a''}{a} \leq \frac{1}{2}
    \frac{\Om_{\mat, 0}}{\Om_{\rad, 0}}
    \frac{g_*}{g_*^0}
    \biggl(\frac{\gS}{\gS^0}\biggr)^{-{4\over3}}
    \HH \, \frac{a_*}{a_0} \lesssim 
    4.5 \times 10^3 \, \HH \, \frac{a_*}{a_0},
    \label{rel_app}
\end{equation}
where the first inequality is an approximated value valid during 
the RD era and it decays during MD.
Hence, the upper bound is valid at all times and its specific value depends
on $a_*/a_0$.
For example, for GWs generated at the EWPT or at the QCDPT, this leads to the
following critical wave numbers, at which the WKB approximation might
break down,
\begin{align}
    &\, \klimap^\EW \sim 1.87 \cdot 10^{-6} \, \HH \,
    \sqrt{\frac{\bigl|1 + \half \alpM\bigr|}{\cT}} \leq
    1.87 \cdot 10^{-6} \, \sqrt{\frac{\bigl|1 + \half \alpM\bigr|}{\cT}},
    \label{kcritEW}\\
    &\, \klimap^\QCD \sim 5.90 \cdot 10^{-5} \, \HH \, \sqrt{\frac{\bigl|1 + \half \alpM\bigr|}{\cT}} \leq
    5.90 \cdot 10^{-5} \, \sqrt{\frac{\bigl|1 + \half \alpM\bigr|}{\cT}}.
    \label{kcritQCD}
\end{align}
In both cases, this term is subdominant up to very large
superhorizon scales and it is bounded by the values in \Eqs{kcritEW}{kcritQCD} since $\HH = \eta_*/\eta$
during RD era as can be seen using \Eq{eqn:a_derivatives_omega}.
This limit can be modified by the inclusion of $\alpM$ and $\alpT$
but unless they take large values these modifications are negligible.

On the other hand, two additional limitations to
the WKB approximation appear due to the $\alpM$ parameter,
\begin{equation}
    \klimaM \sim \frac{\HH}{\cT}\sqrt{\biggl|\frac{\alpM}{2} \Bigl(1 + \half \alpM 
    \Bigr)\biggr|} \lesssim \frac{1}{\cT}\sqrt{\frac{|\alpM|}{2}}, \quad
    \klimaMp \sim \frac{\HH}{\cT}\sqrt{\frac{|\alpM'|}{2}} \leq
    \frac{1}{\cT}\sqrt{\frac{|\alpM'|}{2}},
    \label{WKB_limitations}
\end{equation}
where we can neglect the term $\half \alpM$ in front of 1 in the
first limit
for small values of $\alpM$.
Hence, the WKB limit can break down in Horndeski theories
around the horizon or at larger scales.
In general, we expect the limit from $\alpM$ to be more
restrictive than that from $\alpM'$ and to dominate
at the initial time when $\HH = 1$.
However, this depends on the parameterization of $\alpM$,
which can give different results for both limits; see \Sec{sec:param}.

\subsection{GW spectrum using the WKB approximation}
\label{GW_sp_WKB}

The spectrum of GW energy density can be expressed as
\begin{equation}
\OmGW (k, \eta) = \frac{1}{\rho_{\crit,0}}
\frac{\dd \ln \rho_\GW}
{\dd \ln k} = \frac{1}{6} \biggl(\frac{H_*}{H_0} \biggr)^2 \biggl(
\frac{a_*}{a_0} \biggr)^4 k S_{h'} (k, \eta),
\label{eqn:OmGW_Sh}
\end{equation}
where $2 S_{h'}$ is the spectrum\footnote{
Following ref.~\cite{RoperPol+18}, we define the spectrum
$S_{h'} (k, \eta)$ from the $+$ and $\times$ polarization
modes, giving an extra factor of 2 due to the property
\begin{equation}
    \tilde h'_{ij} (\kk, \eta) \, \tilde h^{'*}_{ij} (\kk, \eta)
    = 2 \bigl[h'_+ (\kk, \eta) \, h^{'*}_+ (\kk, \eta) + 
    h'_\times (\kk, \eta) h^{'*}_\times (\kk, \eta) \bigr].
\end{equation}}
of the scaled strain
derivatives,
\begin{equation}
    \bra{\tilde h_{ij}' (\kk, \eta) \, \tilde h_\ij^{'*} (\kk', \eta)} = (2 \pi)^6 
    \delta^3 (\kk - \kk') \frac{2 S_{h'} (k, \eta)}{4 \pi k ^2},
\label{eqn:Sh}
\end{equation}
and it can be computed using the WKB approximation.
Taking $\ddd' = \frac{1}{2}\alpM\hhh$,
the modified GW energy spectrum can be obtained using the WKB
solution for reduced Horndeski theories, given in \Eq{eqn:h_WKB},
\begin{align}
   \OmGW &\, (k, \eta) = e^{-2\ddd} \OmGW^\GR (k, \eta) \nonumber \\
   \times &\, \Biggl(1 + \half \alpT + \frac{{\alpM^*}^2 \alpM^2 \HH^2}{32 k^4 \cT^2} + \frac{\alpM^* \alpM^2 \HH^2}{8 k^3 \cT^2}
   + \frac{1}{8k^2} \Biggl[\alpM^2 \HH^2\biggl(1 + {1\over\cT^2}\biggr) + {\alpM^*}^2
   \Biggr] +
   \frac{\alpM^*}{2k} \Biggr).
   \label{OmGW_mod_general}
\end{align}
At late times, $\HH^2 = (\eta_*/\eta)^2$ decreases, so we can neglect the
terms that dilute with $\HH^2$,
\begin{equation}
    \OmGW (k, \eta \gg \eta_*) = \half e^{-2\ddd} \OmGW^\GR (k, \eta) \biggl[1 +
    \alpT +\Bigl(1 + \frac{\alpM^*}{2k} \Bigr)^2 \biggr].
    \label{eqn:EGW_mod}
\end{equation}
From \Eq{eqn:EGW_mod}, we can identify the IR and UV
limiting ranges of the spectrum, 
\begin{align}
\OmGW^\IR = \OmGW \Bigl(k \ll \half \alpM^*, \eta \gg \eta_* \Bigr)
= & \, \frac{{\alpM^*}^2}{8 k^2} \, e^{-2\ddd} \, \OmGW^\GR(k,\eta),\\
\OmGW^\UV = \OmGW \Bigl(k \gg \half \alpM^*, \eta \gg \eta_*\Bigr)
= & \, \Bigl(1 + \half \alpT\Bigr) \, e^{-2\ddd}  \, \OmGW^\GR(k,\eta).
\label{GW_UV_mod}
\end{align}
Note that the factor ${\alpM^*}^2$ is part of the IR limit,
indicating that the IR enhancement of the form $k^{-2}$ holds regardless of the sign of $\alpM$.
The critical $k_\crit = \half \alpM^*$ marks the relative strength of the
$\ooo(k^{-2})$ term in \Eq{eqn:EGW_mod},
which means it also indicates where the IR regime begins, i.e., $k^{-2}$ becomes dominant.

\EEq{OmGW_mod_general} shows that the GW spectrum in the IR regime can present up to $k^{-4}$
GR spectrum but this and other terms vanish as time
evolves since they are proportional to $\HH^2$.
Hence, at late times, we end up with the GR spectrum
amplified by $e^{-2\ddd}$ at all wave numbers
with a $k^{-2}$ enhancement in the IR regime,
proportional to the additional $\eigth {\alpM^*}^2$
factor, as shown in \Eq{eqn:EGW_mod}.
In the UV regime, we find the $e^{-2\ddd}$ enhancement and an
additional factor $1 + \half \alpT$.

We have found that, according to the WKB approximation,
the parameter $\alpM$ introduces changes in the spectral shape
at $k \leq k_\crit = \half \alpM^*$ that do not dilute as time
evolves and depend on the value of $\alpM$ only at the time of
GW generation.
On the other hand, note that, using \Eq{WKB_limitations}, the WKB approximation breaks down at
$k \sim \max(\klimaM, \klimaMp)$.
This means that when $\HH \sim 1$, the spectral changes occur at wave numbers around and below the
critical $k_\crit$,
where the WKB approximation might not be valid (for $\cT \sim 1$).

Finally, we note that $\alpT$ introduces no changes to the spectral shape.
Instead, its impact corresponds to an enhancement or depletion
of the GW spectrum by a factor
\begin{equation}
    \Delta \OmGW (k, \eta \gg \eta_*) =
    \frac{\OmGW - \OmGW^\GR e^{-2\ddd}}{\OmGW^\GR e^{-2\ddd}} =
    \half \biggl[\alpT (\eta) + \frac{\alpM^*}{k} + 
    \frac{{\alpM^*}^2}{4 k^2} \biggr],
    \label{DeltaGW}
\end{equation}
which includes the increase on amplitude from $\alpM$ via the
damping factor $\ddd$.
Note that we define the increase on $\OmGW(k)$ after
compensating $\OmGW^\GR$ by $e^{-2\ddd}$ to see the
additional effects of $\alpM \neq 0$ on the spectral
shape and the change in amplitude due to $\alpT \neq 0$.
The sign of $\alpT$ determines if the GW spectrum is 
amplified or decreased with respect to that obtained from
GR.

We investigate the resulting spectra in modified GR using
numerical simulations in \Sec{sec:numerical_solutions}
with the objective to test the validity of the WKB
approximation and its potential limitations and to confirm the resulting GW spectra
when Horndeski parameters are introduced due to modified
theories of gravity.

\section{Phenomenological parameterizations}
\label{sec:param}

In Sec.~\ref{GW_sp_WKB}, we noted that, according to the WKB 
approximation, $\alpT \neq 0$ induces modifications in the total GW energy density, but not
in its spectral shape.
Realistically, the tensor speed excess at the present day is observationally constrained to be $\alpTz\lesssim\ooo(10^{-15})$ 
by the binary neutron star merger GW170817 and its gamma-ray burst GRB 170817A~\cite{LIGOScientific:2017zic}.
This constraint can be circumvented if $\alpT$ is either frequency- or time-dependent, such that larger deviations could hide outside the LIGO-Virgo frequency 
band~\cite{deRham+18,Baker:2022rhh} or in the past.
In this study, we show that even using larger values of
$\alpT$, constant in time and frequency, the modifications to 
the GW spectrum are negligible, which is seen from the WKB
approximation and can be confirmed with numerical simulations.
%Hence, from now on, we will focus on two cases: {\em i)} 
%constant $\alpT$ with $\alpM = 0$, and {\em ii)}
%AB: I'd put a left bracket in both cases, but am not sure that \em is needed.
%AR: it looks a bit better to me with \em but it is fine by me either way. I don't know about the journal policy.
%YH: I don't have a preference
Hence, from now on, we will focus on two cases: {\em(i)} 
constant $\alpT$ with $\alpM = 0$, and {\em(ii)}
reduced Horndeski theories with $\alpT = 0$, allowing $\alpM$ to take
different values as a function of time.

In particular, we are interested in the perspectives of negative values of $\alpM$, 
as they would enhance the amplitude of the GW spectrum.
Various constraints on the present-day value $\alpMz$ exist in
the literature,
but are significantly less stringent than those on $\alpT$.
Analytically, $-1.0\lesssim\alpMz\lesssim-0.5$ is found (see figure~4 in~\cite{Lombriser+15}).
This is compatible with the results from Monte Carlo simulations,
%AB: here, Monte Carlo simulations are mentioned; should one explain this?
giving $|\alpMz|\lesssim\ooo(1)$ (see figure~2 in~\cite{Arai+17}).
On the other hand, the Planck Collaboration reports
a constraint on $\alpMz \gtrsim -0.1$
(see table 8 of ref.~\cite{Planck:2018vyg}).
Note that they parameterize $\alpM = \alpMz \, (a/a_0)^n$ with
$n \in (0.5, 1)$.
Taking these different constraints into account,
we explore a range of $\alpMz\in[-0.5, 0.3]$.

Besides the simplest consideration of a constant $\alpM$ in time, 
which we call choice 0,
its time-dependent forms can be written in accordance with
specific gravity models~\cite{Odintsov+22}.
Following refs.~\cite{Bellini+14,Nunes+18,Kennedy:2018gtx}, we choose phenomenological forms of $\alpM$
as simplified parameterizations motivated by effective descriptions
of modified 
gravity~\cite{Gleyzes:2017kpi}.
Specifically:
\begin{equation}
\alpM(\eta) = 
\begin{dcases}
\alpMz & (\text{choice 0}),\\
\alpMz\biggl[\frac{a(\eta)}{a_0}\biggr]^n & (\text{choice I}),\\
\alpMz \frac{1}{\Om(\eta)} & (\text{choice II}),\\
\alpMz \frac{1-\Om_\mat(\eta)/\Om(\eta)}{1-\Om_{\mat,0}} & (\text{choice III}),
\end{dcases}
\label{eqn:alpM_param}
\end{equation}
where choices II and III give a value of $\alpM$ proportional
to the percentage of dark energy density
and the combination of dark and radiation
energy densities, respectively, at each time $\eta$ compared to
%their relative amounts at present time.
%AB: added "the"
their relative amounts at the present time.
For choice I,
the following conditions on the parameter values are needed in order to ensure the theory's stability~\cite{Denissenya:2018mqs}:
\begin{equation}
\begin{dcases}
\alpMz > 0 : 0 < n < \threehalf \, \Om_{\mat,0} \simeq \half,\\
\alpMz < 0 : n > \threehalf.
\end{dcases}
\end{equation}

\begin{figure}[t]
\centering
\includegraphics[width=0.48\textwidth]{alphaM_params_EWPT.pdf}
\includegraphics[width=0.48\textwidth]{params_WKB_lims_EWPT.pdf}
\caption{Time evolution of $\alpM\hhh$ (left panel), which
contributes to the change of amplitude over time,
$h_{ij} \sim e^{-\ddd}$ with $\ddd = \half \int^\eta \alpM \HH \dd \eta'$,
and $\klimaM/\cT$ and $\klimaMp/\cT$ (right panel),
which are the
terms neglected compared to $k \cT$ under the WKB approximation;
see \Eq{WKB_limitations}.
All four parameterization choices (0 to III) are shown and,
for illustrative purposes, $n = 1$ is chosen.
The main figures show the full evolution whereas
the inset (in left panel) shows only times after the onset of matter domination.
we show in \Fig{fig:alpM_Hcal}, for the different time dependencies.
We have taken $\eta_*$ to correspond to the EWPT for the specific
values in the axes, which puts the present time at
$\eta_0/\eta_* \simeq 2.38 \times 10^{13}$.
}
\label{fig:alpM_Hcal}
\end{figure}

%To provide an intuition of the aforementioned $\alpM$ parameterizations,
%AB: of? Maybe on or "regarding"?
To provide an intuition on the aforementioned $\alpM$ parameterizations,
we show in \Fig{fig:alpM_Hcal} (left panel), for the different time dependencies
given in \Eq{eqn:alpM_param}, the evolution of $\alpM\hhh$, which
characterizes the growth or damping produced by $\alpM$, according to
the WKB approximation [see \Eq{eqn:ddd_DelT}].
We also show the time evolutions of $\sqrt{\alpM} \HH$ and
$\sqrt{\alpM'} \HH$ (right panel), which correspond to the wave numbers
at which the WKB approximation might break down,
according to \Eq{WKB_limitations}.
\FFig{fig:alpM_Hcal} shows the time evolution from the time
of generation (e.g., EWPT) up to present time.
We see that for choices 0 and III, $\alpM\hhh$ converges to the same values during RD
and $\Lam$D, and only becomes different during MD, as expected.
On the other hand, for choices I and II, $\alpM \HH$ is negligible
for all of RD and most of MD and rapidly increases later on, especially during $\Lam$D.
It converges for all parameterizations to the value $\alpMz$ at the present time.
We have found in \Sec{limitations_WKB} the values of $k$ at which
the assumptions made by the WKB approximations do not hold.
We see in \Fig{fig:alpM_Hcal} (right panel) that for choices 0 and III,
this limiting $k$ is around the horizon scale at the
time of generation, i.e., $\cT \klimaM \sim \sqrt{\alpM^*/2}$
and then it decreases into superhorizon scales since it
is proportional to $\HH = \eta_*/\eta$ (this is because
for these choices, $\alpM \simeq \alpM^*$ during the RD era).
The term $\alpM' \ll \alpM$ for choice III, as can be
seen in \Fig{fig:alpM_Hcal}.
The ratio between $\alpM$ and $\alpM'$ during the RD era
can be found to be
\begin{equation}
    \biggl(\frac{\alpM}{\alpM'}\biggr)^{1\over2}_{\rm III} \sim 10^2 \biggl(\frac{a_*}{a_0}\biggr)^{1\over2}
    \simeq 3 \times 10^{-6} \biggl(\frac{a_*}{8 \times 
    10^{-16}} \biggr)^{1\over2},
    \label{rel_alpMp}
\end{equation}
which only becomes of order $1$ towards the end of the RD era,
when $a_*/a_0 \sim 10^{-4}$, which corresponds to times
much later than the EWPT and the QCDPT.
For choices I and II, although the two terms are of the 
same order, ${\alpM}^{'*} \sim \alpM^*$, their value
at the time of GW generation is much smaller than
$\alpMz$ (for all $n > 0$ in the choice I).
For this reason, in these cases the WKB estimate is completely valid and we do not expect to observe
any relevant spectral change, since the IR enhancement
is determined by the value of $\alpM^*$, as seen in \Eq{DeltaGW}.

Note that, since many efforts of modifying gravity are aimed at
addressing the late-time acceleration of the universe, the
$\alpM$ parameterization choices are commonly constructed to be dominant in the late universe.
This is indeed the case for choices I and II here.
If one considers only the MD and $\Lam$D history of the universe,
choice III is also well-motivated to be relevant at late times.
However, here we explore the entire history of the Universe from well within RD,
which means that $\alpM\hhh$ becomes dominant for choice III both at early and late times.
We note that the physical motivations of the choices here can be potentially ambiguous.
Even though there exists a wide range of discussions on modified
gravity during inflation~\cite{Nunes+18,Yoshimura:2022kyd} and around recombination~\cite{Zumalacarregui:2020cjh,Braglia:2020auw},
where $\alpM$ is essentially a free function of the scalar field,
there seems to be a 
relative lack of numerical studies on the effects of $\alpM$
during RD
(note, however, a brief discussion of $\alpM = -1$ during RD in ref.~\cite{Linder:2021pek}).
Hence,
we emphasize that the aim of this work is to provide an
understanding of the phenomenological behavior of the GW spectrum due to $\alpM$ for signals produced during RD.

\section{Numerical solutions}
\label{sec:numerical_solutions}

To explore the limits and validity of the WKB approximation, we use
the {\sc Pencil Code} to numerically solve the GW equation under
modified gravity, given in \Eq{eqn:GW_mod}.
The {\sc Pencil Code} is a highly parallelized modular code that can be
used to solve various differential equations \cite{pencil}.
In the context of cosmological GWs, it has previously been used to study GWs generated by hydrodynamic and
MHD stresses in the early universe \cite{RoperPol:2019wvy}.
It uses a GW solver that advances the strains at each time step
sourced by the anisotropic stresses that are separately computed as the
solution to the MHD equations \cite{RoperPol+18}.
Previous numerical works solved the GW equation under GR, while in the
present work, we have extended the code to solve \Eq{eqn:GW_mod}; see
\App{sec:numerical_scheme} for further details on the numerical scheme.

\subsection{Initial condition and time stepping schemes}

In the current work, we focus primarily on the propagation rather than the production of GWs.
Therefore, we have adapted the {\sc Pencil Code} to evolve an initial GW spectrum
in the absence of sources
with a spectral shape and amplitude based on those obtained in previous studies (see, e.g., refs.~\cite{He:2021bqm,RoperPol:2022iel}).
Hence, \Eq{eqn:GW_mod} is solved in one-dimension,
i.e., in $k>0$, in order to improve the efficiency
of the code to study the propagation of a GW background
along the cosmological history of the universe.
We take the initial spectrum for the time derivative of the strains to be a smoothed double broken power law described by
\begin{equation}
S_{h'}(k, \eta_*) = S_{h'}^*
\frac{2^{1\over\alpha_2} \biggl[1 + \Bigl(\frac{\kf}{\kb}\Bigr)^{a - b}\biggr]
\Big(\frac{k}{k_*}\Big)^a}{\bigg[1 +
\Big(\frac{k}{\kb}\Big)^{(a - b)\alpha_1}\bigg]^{1\over \alpha_1}\bigg[1 +
\Big(\frac{k}{k_*}\Big)^{(b + c)\alpha_2}\bigg]^{1\over \alpha_2}
},
\label{sp_ini}
\end{equation}
where $S_{h'}^*$ is approximately the initial peak amplitude, 
$\eta_* = \HH_*^{-1}$ corresponds to the time of GW production
during the RD era,
and $\alpha_1 = \alpha_2 = 2$ are fixed smoothness parameters.
We choose the slope in the IR range to be $k^2$, set by $a = 2$
(as is expected for causal sources of GWs, as those produced
for example during a phase transition)
up to the break wave number $\kb = 1$, which corresponds to the horizon scale.
At intermediate wave numbers, $\kb \leq k \leq \kf$, the slope becomes
$b = 0$, as found, for example, for MHD turbulence in ref.~\cite{RoperPol:2019wvy}
and for sound waves in ref.~\cite{Hindmarsh:2019phv}, although in
%AB: could you please check? I think there advertise here a k^9 spectrum!
%AR: they have k^9 spectrum at small k and a short intermediate linear range ~ k
%AB: Maybe we should clarify this better that they talk about an intermediate range.
both cases, the break $\kb$ is not exactly at the horizon but depends
on the dynamics of the source; see, e.g., ref.~\cite{RoperPol:2022iel}.
The peak wave number is chosen to be $k_* = 10$,
followed by the slope in the UV range as $k^{-{11\over3}}$, set by $c = {11\over3}$.
This corresponds to the spectrum obtained for
Kolmogorov-like MHD turbulence \cite{RoperPol:2019wvy}.
Note, however, that this part of the spectrum corresponds to subhorizon
scales, which are described accurately by the WKB approximation
and, hence, the resulting spectral shape is not expected to be modified in this range.
On the other hand, around the horizon or at larger scales, depending
on the value of $\alpM^*$, the resulting GW spectral shape might be modified by the inclusion of
an additional IR branch $\sim k^{-2} S_{h'} (k)$, as predicted in \Sec{GW_sp_WKB} by the WKB approximation.
Finally, the resulting spectrum is expected to be enhanced at all wave numbers
by a factor $e^{-2\ddd} (1 + \half \alpT)$ owing to the presence
of non-zero Horndeski parameters $\alpM$ and $\alpT$
in modified gravity.
After compensating by this enhancement, which affects equally to all wave numbers, we model the spectral modifications using the function $\xi (k)$, such that $\OmGW (k) \propto \xi(k) \, \OmGW^\GR (k)$ (see \Tab{tab:run_params}).
\FFig{fig:sp_sup_WKB} shows the resulting GW spectrum estimated
using the WKB approximation for different values of $\alpM^*$.
Note, however, that the assumptions made by the WKB approximation break
down around the critical wave number at which we expect the relevant spectral modifications.

%YH: might be good to show the (very small) slope between k_break and k_peak caused by alpM (could be seen in the commented-out figure), in case we want to talk about it in the observation section and have something to refer to
%YH: if you also think so, could you add it to your figure?:)
\begin{figure}[t]
\centering
\includegraphics[width=.49\textwidth]{spectrum_WKB_late_times}
%\includegraphics[width=.6\textwidth]{sp_sup_WKB}
\caption{Expected WKB spectrum [see \Eq{eqn:EGW_mod}] of the strain derivatives $S_{h'} (k)$
at late times due to non-zero $\alpha_{\rm M}^*$ at the time of generation for an initial GW background characterized by a double broken power law, given in \Eq{sp_ini}.
The black curve indicates the GR solution and the colored curves 
correspond to different values of $\alpM^*$, with the solid and
dashed ones being positive and negative values, respectively.
Vertical dashed lines indicate 
the estimated value of $\klimaM$, below which the WKB approximation
might break.
The inset zooms
in on the behaviors around the horizon at $k = \kb = 1$.}
\label{fig:sp_sup_WKB}
\end{figure}

%AR: to check values to see if we agree
\begin{table}[t]
\centering
\begin{tabular}{c|c|r|r|r|c|l}
Parameters & Run & $\alpMz$ & $n$ & $\alpTz$ & $\Del\EEGW/\EEGW$ & $\xi(k)$\\\hline
% 46000cosmo_1D_off, 46000cosmo_1D_off_LD2 (GR solution)
Choice 0   & T0A & $0$     & --     & $-0.5$ & $-2.5\times10^{-1}$        & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpT_m05
($\alpT$)  & T0B & $0$     & --     & $-0.2$ & $-1.0\times10^{-1}$        & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpT_m02
           & T0C & $0$     & --     & $0.2$  & $\ \ \, 1.0\times10^{-1}$  & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpT_p02
           & T0D & $0$     & --     & $0.5$  & $\ \ \,2.5\times10^{-1}$   & $k^0\,\,\forall\,\,k$ \\\hline% 46000cosmo_1D_alpT_p05
Choice 0   & M0A & $-0.5$  & --     & 0      & $\ 3.3\times10^7$          &$k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m05, 46000cosmo_1D_alpM_m05_LD2, 46000cosmo_1D_alpM_m05_wav1e7, 46000cosmo_1D_alpM_m05_wav1e7_LD2
($\alpM$)  & M0B & $-0.3$  & --     & 0      & $\ 3.2\times10^4$          &$k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m03, 46000cosmo_1D_alpM_m03_LD2, 46000cosmo_1D_alpM_m03_wav1e7, 46000cosmo_1D_alpM_m03_wav1e7_LD2
           & M0C & $-0.1$  & --     & 0      & $\ 3.1\times10^1$          &$k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m01, 46000cosmo_1D_alpM_m01_LD2, 46000cosmo_1D_alpM_m01_wav1e7, 46000cosmo_1D_alpM_m01_wav1e7_LD2
           & M0D & $-0.01$ & --     & 0      & $\ \ \, 4.1\times10^{-1}$  &$k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m001, 46000cosmo_1D_alpM_m001_LD2, 46000cosmo_1D_alpM_m001_wav1e7, 46000cosmo_1D_alpM_m001_wav1e7_LD2
           & M0E & $0.1$   & --     & 0      & $-9.7\times10^{-1}$        & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_p01, 46000cosmo_1D_alpM_p01_LD2, 46000cosmo_1D_alpM_p01_wav1e7, 46000cosmo_1D_alpM_p01_wav1e7_LD2
           & M0F & $0.3$   & --     & 0      & $\!\!\! -1.0\times10^0$    & $k^0\,\,\forall\,\,k$ \\\hline% 46000cosmo_1D_alpM_p03, 46000cosmo_1D_alpM_p03_LD2, 46000cosmo_1D_alpM_p03_wav1e7, 46000cosmo_1D_alpM_p03_wav1e7_LD2
Choice I   & M1A & $-0.5$  & $2$    & 0      & $\ 2.5\times10^{-1}$       & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_m05_a2, 46000cosmo_1D_alpM_m05_a2_LD2
($\alpM$)  & M1B & $-0.3$  & $2$    & 0      & $\ 1.5\times10^{-1}$       & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_m03_a2, 46000cosmo_1D_alpM_m03_a2_LD2
           & M1C & $-0.1$  & $2$    & 0      & $\ 4.7\times10^{-2}$       & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_m01_a2, 46000cosmo_1D_alpM_m01_a2_LD2
           & M1D & $0.1$   & $0.4$  & 0      & $\,-2.2\times10^{-1}$      & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_p01_a04, 46000cosmo_1D_alpM_p01_a04_LD2
           & M1E & $0.3$   & $0.4$  & 0      & $\,-5.2\times10^{-1}$      & $k^0\,\,\forall\,\,k$ \\\hline% 46000cosmo_1D_alpM_p03_a04, 46000cosmo_1D_alpM_p03_a04_LD2
Choice II  & M2A & $-0.5$  & --     & 0      & $\ 2.9\times10^{-1}$      & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_m05_Lambda, 46000cosmo_1D_alpM_m05_Lambda_LD2
($\alpM$)  & M2B & $-0.3$  & --     & 0      & $\ 1.6\times10^{-1}$      & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_m03_Lambda, 46000cosmo_1D_alpM_m03_Lambda_LD2
           & M2C & $-0.1$  & --     & 0      & $\ \ \,5.3\times10^{-2}$  & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_m01_Lambda, 46000cosmo_1D_alpM_m01_Lambda_LD2
           & M2D & $0.1$   & --     & 0      & $-5.0\times10^{-2}$       & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_p01_Lambda, 46000cosmo_1D_alpM_p01_Lambda_LD2
           & M2E & $0.3$   & --     & 0      & $-1.4\times10^{-1}$       & $k^0\,\,\forall\,\,k$ \\\hline% 46000cosmo_1D_alpM_p05_Lambda, 46000cosmo_1D_alpM_p05_Lambda_LD2
Choice III & M3A & $-0.5$  & --     & 0      & $\ 3.5\times10^8$         & $k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m05_mat, 46000cosmo_1D_alpM_m05_mat_LD2, 46000cosmo_1D_alpM_m05_mat_wav1e7, 46000cosmo_1D_alpM_m05_mat_wav1e7_LD2
($\alpM$)  & M3B & $-0.3$  & --     & 0      & $\ 1.3\times10^5$         & $k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m03_mat, 46000cosmo_1D_alpM_m03_mat_LD2, 46000cosmo_1D_alpM_m03_mat_wav1e7, 46000cosmo_1D_alpM_m03_mat_wav1e7_LD2
           & M3C & $-0.1$  & --     & 0      & $\ 5.0\times10^1$         & $k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m01_mat, 46000cosmo_1D_alpM_m01_mat_LD2, 46000cosmo_1D_alpM_m01_mat_wav1e7, 46000cosmo_1D_alpM_m01_mat_wav1e7_LD2
           & M3D & $-0.01$ & --     & 0      & $\ \ \,4.8\times10^{-1}$  & $k_\IR^{-2}\,\,\&\,\,k_\UV^0$ \\% 46000cosmo_1D_alpM_m001_mat, 46000cosmo_1D_alpM_m001_mat_LD2, 46000cosmo_1D_alpM_m001_mat_wav1e7, 46000cosmo_1D_alpM_m001_mat_wav1e7_LD2
           & M3E & $0.1$   & --     & 0      & $-9.8\times10^{-1}$       & $k^0\,\,\forall\,\,k$ \\% 46000cosmo_1D_alpM_p01_mat, 46000cosmo_1D_alpM_p01_mat_LD2, 46000cosmo_1D_alpM_p01_mat_wav1e7, 46000cosmo_1D_alpM_p01_mat_wav1e7_LD2
           & M3F & $0.3$   & --     & 0      & $\!\!\!-1.0\times10^0$    & $k^0\,\,\forall\,\,k$   % 46000cosmo_1D_alpM_p03_mat, 46000cosmo_1D_alpM_p03_mat_LD2, 46000cosmo_1D_alpM_p03_mat_wav1e7, 46000cosmo_1D_alpM_p03_mat_wav1e7_LD2
\end{tabular}
%YH: to check spectra from choices 0 and III (down to k = 1e-7), if same, keep only choice 0
\caption{Parameters used for the numerical studies:
for all runs, $k_1 = 10^{-3}$ is the smallest wave number and
$N = 46\,000$ is the number of grid points in one dimension.
$\Del\EEGW/\EEGW\equiv(\EEGW - \EEGW^\GR)/\EEGW^\GR$ indicates the 
relative changes to the total energy density at the present time.
In the last column, we briefly state the changes in the final
energy spectrum represented by $\xi(k)$, such that
$\OmGW(k)\propto \xi(k)\OmGW^\GR(k)$.}
\label{tab:run_params}
\end{table}

Table~\ref{tab:run_params} summarizes the input parameters for the numerical studies.
The values of $\alpTz$ are chosen to be much larger than the current constraints
in order to later show the relative insignificance of
$\alpTz$ even with unrealistically large values.
The choices of $\alpMz$ are made in line with the limits discussed in \Sec{sec:param}.

For series T0 (T0A--T0D),
we evolve the solution entirely with increasing time steps,
such that $\eta_{\rm next} = \eta_{\rm current}(1 + \del n_\incr)$ 
with $\del n_\incr = 0.01$, leading to equidistant points in logarithmic time spacing.
For series M0 (M0A--M0F), M1 (M1A--M1D), M2 (M2A--M2D), and M3 (M3A--M3F), 
we keep the nonuniform time scheme during RD and MD but switch to linear
time steps during $\Lam$D such that $\eta_{\rm next} = \eta_{\rm current}
+ \del\eta$ with $\del\eta \, \HH_* = 5\times10^9$.
We chose such time evolutions to improve the accuracy of late-time results for reduced Horndeski models,
especially the modifications that they present in the IR limit (see last 
column in table~\ref{tab:run_params}).
We show in \App{sec:numerical_accuracy} that
decreasing the time step below $\delta \eta \, \HH_* = 5\times 10^9$
does not affect the IR range of the spectra, which
indicates that the observed modifications are not due to numerical accuracy.
Since series T0 does not exhibit $k$-dependent modifications, it
does not require the additional
computational effort.
The choice of time schemes and their numerical accuracy is further discussed
in \App{sec:numerical_accuracy}.

In the simulations, we consider the initial time to be the EWPT
with a temperature scale $T_* \sim 100\,$GeV and the number of relativistic and 
adiabatic DOFs are $\gS \approx g_* \sim 100$,
which yields the values of $a_*$ and $\HH_*$ in \Eqs{as_a0}{HHs}.
With these values, we can then solve \Eq{eqn:GW_mod} in units of
the normalized time $\eta \HH_*$, by mapping the parameterizations in \Eq{eqn:alpM_param} $\alpM(\eta) \map \alpM(\eta/\eta_*)$ and using the results from Friedmann equations for
$\HH$ and $a''/a$ (which are already normalized since $a'$ is
computed as the derivative with respect to $\eta \HH_*$),
given in \Eq{eqn:a_derivatives_omega}.
The effects on $S_{h'} (k, \eta)$ of the specific choice of the time at which the
GWs are generated only appear via the relative magnitude of the terms $a''/a$
and $\alpM'$ that involve time derivatives (normalized by $\HH_*$)
compared to $k$.
These terms have been parameterized in \Eqs{rel_app}{rel_alpMp},
respectively, and their magnitude has been discussed.
We have shown that the term $a''/a$ can only induce
modifications to the solution at scales several orders of magnitude
above the horizon scale, while $\alpM'$ is only of the order of
$\alpM$ for choice III (see \Sec{sec:param}) at very late times within
the RD era.
For other choices of $\alpM (\eta)$, either $\alpM' = 0$
(choice 0) or $\alpM$ itself is orders of magnitude below its present-time value $\alpMz$ (choices I and II).
We find that the WKB approximation
is expected to be valid around $k \gg \klimaM \sim \sqrt{|\alpM^*|/2}$ and, on the other hand,
it predicts an enhancement $\xi (k) = k^{-2}$ in the IR regime $k \lesssim k_\crit \leq \klimaM$, which is not on the range
of validity of the WKB estimate.
However, we characterize $\xi (k)$ using numerical simulations
and we find an IR enhancement $\xi (k) \propto k^{-2}$ that does not
in general follow the spectral shape predicted by WKB and that
can become shallower at smaller $k$.

The position of this spectral change is a fixed fraction of
the Hubble horizon at the
time of generation, determined by the value of $\alpM^*$, that
does not depend on the specific value of $\eta_*$.
Finally, the time of generation determines the range of frequencies
where we observe the signal as well as its amplitude via \Eq{eqn:OmGW_Sh}.

The peak amplitude $S_{h'}^*$ can be chosen to represent a specific
model.
For example, a value of $S_{h'}^* = 3.5 \times 10^{-10}$ is used in our runs (see table~\ref{tab:run_params}) and it would
produce a normalized initial total energy density\footnote{
The turbulent and GW energy densities are normalized to the 
radiation energy density, such that $\EEGW$ during RD can be
computed as
\begin{equation}
     \EEGW (\eta) = \frac{1}{6} \int_0^\infty S_{h'} (k, \eta) \dd k.
\end{equation}
}
$\EEGW^*\simeq 10^{-9}$,
which corresponds to a vortically turbulent source energy density of roughly
$\eee_\turb^* \simeq 3 q^{-1} \times 10^{-4}$.
This is related via
\begin{equation}
\EEGW^* = (q\eee_\turb^*/k_*)^2,
\end{equation}
where $q$ is an empirically determined coefficient for a 
specific type of turbulence source $\eee_\turb^*$, and is found
to be of the order of unity or larger (up to $\sim$\,5), depending on the production
mechanism of the source \cite{RoperPol:2019wvy,Brandenburg:2021bvg,RoperPol:2021xnd}.

\subsection{Time evolution}

For constant $\alpT \neq 0$ and $\alpM = 0$, the GW energy density
stays constant in time and its magnitude is
modified by the specific value of $\alpT$, as it can be predicted from the analytical solution to \Eq{eqn:GW_mod}
in the absence of sources during RD, given in \Eqs{sol_alpT}{eqn:EGW_mod} under the
WKB approximation.
We find excellent agreement between the WKB estimate and the
numerical solution of the GW spectra, enhanced by a factor $1 + \half \alpT$.
Note that this corresponds to a boost of energy for
$\alpT > 0$ and a depletion for $\alpT < 0$.
The relative changes upon the GR solutions are of order
$\alpT \sim \ooo(10^{-1})$
(see \Tab{tab:run_params}),
which are much larger than the constraint $\alpTz\lesssim\ooo(10^{-15})$.
Therefore, the changes due to realistic 
values of $\alpT$ on the total GW energy density are negligible.

\FFig{fig:ts_all} presents the time evolution of the total GW 
energy density $\EEGW(\eta)$ of the runs in series M0 (upper left), M1 (upper right), M2 (lower left), and M3 (lower right).
In these runs, $\alpM$ follows each of the parameterizations given
in \Eq{eqn:alpM_param} and $\alpT = 0$.
The numerical solutions (dashed lines) and the WKB 
approximations (solid lines) agree on the total GW energy
density time dependence, enhanced or depleted by a factor
$e^{-2\ddd}$.
%
\begin{figure}[t]
\centering
\includegraphics[width=0.49\textwidth]{ts_M0}
\includegraphics[width=0.484\textwidth]{ts_M1}
\includegraphics[width=0.515\textwidth]{ts_M2}
\includegraphics[width=0.475\textwidth]{ts_M3}
\caption{Time evolution of the total energy density $\EEGW$.
On the upper right panel, solid and dashed curves respectively
indicate $\alpM$ choice 0 and III, which overlap completely.
The two lower panels show the evolution during $\Lam$D, with the corresponding insets showing the full evolution.
In all panels, the WKB estimate is shown as square dots and the GR solution as the black solid line.
All runs are initialized at EWPT such that the present day is represented by $\eta_0/\eta_*\simeq 2.38\times10^{13}$.}
\label{fig:ts_all}
\end{figure}
%
\begin{figure}[t]
\centering
%\includegraphics[width=.496\textwidth]{time_series_alpM_choice0_error_WKB.pdf}
%\includegraphics[width=.496\textwidth]{time_series_alpM_choiceI_error_WKB.pdf}
\includegraphics[width=0.48\textwidth]{ts_diff_M0}
\includegraphics[width=0.46\textwidth]{ts_diff_M1}
\includegraphics[width=0.48\textwidth]{ts_diff_M2}
\includegraphics[width=0.46\textwidth]{ts_diff_M3}
\caption{Relative error on the time series
of the WKB approximation compared to the numerical simulations,
$\varepsilon^\WKB ({\cal E}_\GW)$.}
\label{fig:ts_error_WKB}
\end{figure}
Regardless of the specific parameterization,
we find an enhancement or a suppression of $\EEGW$ over time
for negative or positive values of $\alpMz$, respectively.

Choices 0 and III yield similar results, with a time evolution
$\EEGW (\eta) \sim \eta^{-\alpMz}$ and $\EEGW (\eta) \sim \eta^{-1.5 \, \alpMz}$ during RD era, and $\EEGW \sim \eta^{-2 \alpMz}$ and
$\EEGW (\eta) \sim \eta^0$ during MD, respectively.
This agrees with the expectation since choice III mostly differs from choice 0 during MD.
During RD, both parameterizations are almost the same and just differ
by a factor $1/(1 - \Om_{\mat, 0}) \simeq 1.46$, leading to the different
scalings.
In other words, the evolution of $\EEGW$,
proportional to $a^{-\alpM}$ when $\alpM$ is constant (e.g., for choice 0, and for choice III during RD era),
is determined by the cosmic expansion itself.
Note that choices 0 and III induce an enhancement or suppression in the total energy, $\Del\EEGW/\EEGW \equiv (\EEGW - \EEGW^\GR)/\EEGW^\GR$,
that heavily depends on the values of $\alpMz$,
i.e., $\Del\EEGW/\EEGW$ can range from $\ooo(10^{-1})$ at $\alpMz = -0.01$ 
to $\ooo(10^7)$ (choice 0) and $\ooo(10^8)$ (choice III) with $\alpMz = -0.5$.
In general, choice III leads to a larger enhancement due to larger
values of $\alpM$ during RD than choice 0.
The potential implications of such a large GW energy density
enhancement are discussed in \Sec{sec:obs_implications}.

%YH: rewrote the paragraph below to be in line with new results
For both choices I and II (see the upper right and lower left panels of \Fig{fig:ts_all}),
the modified GW solutions remain close to their GR counterparts for most of the time,
and rapidly depart from GR as $\eta$ enters $\Lam$D and approaches the present day.
%This is also expected, since for both scenarios, with $\alpM$ being
%respectively proportional to the scale factor and the dark energy density, 
%the difference remains negligibly small for most of the cosmic
%history but quickly becomes significant as the $\Lam$D starts.
%AR: slightly rephrased, I am not sure this is the reason, it is more like the values of \alpM are negligible till LD, right? what difference do you refer to?
%YH: yes, I meant \alpM ~0 till LD. now it reads more clearly:)
This is expected since $\alpM$ is
proportional to the scale factor and the dark energy density
for choices I and II, respectively, and hence, the values of $\alpM$ are negligibly small for most of the cosmic
history until $\Lam$D era.
%AR
For this reason, for the same values of $\alpMz$,
the final values of $\EEGW$ in both of these cases are significantly lower than those in choices 0 and III,
where the modifications are accumulated from RD onward.

In \Fig{fig:ts_all}, the differences between the WKB and the numerical solutions are indistinguishable.
We quantify in \Fig{fig:ts_error_WKB} the relative error between the two,
defined to be $\varepsilon^\WKB(\EEGW)\equiv[\EEGW^\WKB (\eta) - \EEGW (\eta)]/\EEGW(\eta)$.
In all panels, there exists a brief but relative large error region,
amounting to $\lesssim10\%$,
around the initial time $\eta_*$.
This is due to the sinusoidal oscillations of each $k$ mode that are present in GR and modified gravity alike [see equations \eqref{GR_WKB} and \eqref{eqn:h_WKB}].
After the initial oscillations settle down,
%$\varepsilon^\WKB(\EEGW)$ for choices 0 and III decrease over the majority of time in RD,
$\varepsilon^\WKB(\EEGW)$ for choices 0 and III decrease over the majority of time in RD.
%before increasing until the present day (for choice 0)
%or increasing to a plateau and reducing to zero again at the present day (for choice III).
%However, the late-time growth of errors in choice 0 is rather limited such that even with the largest value of $\alpMz = -0.5$,
%the present day relative error yields a tolerable $\varepsilon^\WKB(\EEGW)\sim3\%$.
%A similar amount of error around $3\%$ is also found as the maximum level of the late-time plateau in choice III.
Although they both increase somewhat later on,
the maximum error at the final time is only $\varepsilon^\WKB(\EEGW)\sim3\%$.
%AR: I think this is a bit too much technical details that don't seem very relevant, I think we could just say that at late times the maximum error is 3%.
%YH: modified above, how about now?
For choices I and II, due to the negligible impact of $\alpM$ during RD,
the GW solutions settle down to the same magnitude after the initial oscillations.
Therefore, the relative errors during RD remain zero.
During $\Lam$D, the errors grow as the effects of $\alpM$ become more significant.
But this is also limited, since at most $\varepsilon^\WKB(\EEGW)\lesssim1\%$ is found at the present day,
lower than in choices 0 and III.

\subsection{Energy spectrum}
\label{ssec:energy_spectrum}

\begin{figure}[t]
\centering
%AR: are choice O and III identical also in the numerical results??
%YH: added preliminary plots, please suppress and replace with better ones if you have those
\includegraphics[width=0.505\textwidth]{sp_sup_alpM0_neg}
\includegraphics[width=0.485\textwidth]{sp_sup_alpM0_pos}
%\includegraphics[width=.496\textwidth]{spectra_alpM_choice0.pdf}
%YH: suppressed spectrum_WKB_late_times.pdf, as you have shown it in the beginning of the section
%\includegraphics[width=.496\textwidth]{spectrum_WKB_late_times.pdf}
\caption{Saturated final energy spectra compensated by the growth on the
total GW energy density
$e^{2\ddd}S_{h'}(k)/S_{h'}^*$ for the runs with the choice 0 of 
$\alpM$ parameterization.
Left and right panels represent negative and positive values of $\alpMz$, respectively.
The GR solution is shown in black, while the dashed and dotted curves are the WKB approximation and numerical solutions, respectively.
Vertical dotted lines indicate the corresponding $\klimaM$.}
%AR: to prepare the figure for choice III and separate positive and negative values of \alpM
\label{fig:sp_all}
\end{figure}
In \Fig{fig:sp_all}, we show
the final energy spectrum normalized against the initial amplitude $S_{h'}^*$ and compensated by the factor $e^{-2\ddd}$
to study the changes on the spectral shape.
The saturated spectra, averaged over oscillations at late times, are shown for the choice 0 of $\alpM$
parameterization,
since the results from choice III exhibit similar spectral features,
%AR: to be confirmed if they are the same
and the rest of runs ($\alpT$ choice 0, $\alpM$ choices I and II) show the same spectral slopes as in GR.
As a reminder of the potential limitations of WKB approximation,
we mark with vertical lines the values of $\klimap^\EW$ (solid), $\klimap^\QCD$ (dashed), and $\klimaM$ (dash-dotted).
We note that $\klimaM$ occurs at wave numbers larger than the IR regime characterized by $k_\crit$.
This implies that the analytically calculated values of $k_\crit = |\alpM^*|/2$ can potentially be unreliable.
However, \Fig{fig:kcrit_alpM} demonstrates that $k_\crit$ is numerically found to coincide with the analytic estimations very well.
\begin{figure}[t]
\centering
\includegraphics[width=0.5\textwidth]{kcrit_alpM}
\caption{Relation between $k_\crit$ and
$|\alpM^*|$ found in the numerical simulations (red crosses), compared
to the one expected from the
WKB approximation (green squares).}
\label{fig:kcrit_alpM}
\end{figure}

%YH: added saturated spectra zoomed-in around the peak
\begin{figure}[t]
\centering
\includegraphics[width=0.505\textwidth]{sp_sup_alpM0A_peak}
\includegraphics[width=0.485\textwidth]{sp_sup_alpM0B_peak}
\includegraphics[width=0.505\textwidth]{sp_sup_alpM0C_peak}
\includegraphics[width=0.485\textwidth]{sp_sup_alpM0D_peak}
\includegraphics[width=0.505\textwidth]{sp_sup_alpM0E_peak}
\includegraphics[width=0.485\textwidth]{sp_sup_alpM0F_peak}
\caption{Saturated final energy spectra around between $\kb$ and $\kpeak$. In all panels, numerical results are represented by scattered dots, WKB approximation is shown in solid curves, and GR solution is in black.}
\label{fig:sp_sup_peak}
\end{figure}

Comparing the numerical results to the expected spectra obtained
using the WKB approximation [see \Fig{fig:sp_sup_WKB}],
we note two main differences.
%YH: "in the first place"?
%In first place, the IR enhancement found under the WKB approximation
In the first place, the IR enhancement found under the WKB approximation
only appears for negative values of $\alpMz$
whereas for positive values there are no significant modifications to the spectral shape.
In addition, the modifications to the spectral shape are slightly
different to those predicted by the WKB estimate.
At wave numbers close to where the IR regime starts, $k\lesssim k_\crit$,
the $k^{-2}$ enhancement is reproduced by the numerical results,
flattening that part of the original GR spectrum to $k^0$.
But as $k$ decreases, $k\ll k_\crit$,
the numerical spectra exhibit slopes steeper than $k^0$
that is not predicted by the WKB approximation, which maintains the $k^0$ slope throughout all scales $k\leq k_\crit$.
The departure in the IR regime occurs at different scales depending on the values of $\alpMz$.
For $\alpMz = -0.01$, this is the least obvious,
where the enhanced spectrum still keeps a roughly $k^0$ shape even at $k\ll k_\crit$.
For $\alpMz = -0.5$, however,
the difference becomes visible when $k$ is still relatively close to the start of the IR regime at $k\sim k_\crit$.

In \Fig{fig:sp_sup_peak}, we inspect the modified spectra between $\kb$ and $\kpeak$,
where the thick solid lines are from WKB approximation [recall the inset of \Fig{fig:sp_sup_WKB}].
In terms of the spectrum amplitudes, for all $\alpMz < -0.01$, we see that WKB underestimates in this subrange compared to the numerical results.
For $\alpMz\geq-0.01$, WKB approximation gradually approaches but does not seem to overshoot the numerical results.
In terms of the spectral slopes, both WKB and numerical solutions give the same results for all values of $\alpMz$ considered.
However, note that the changes in the spectral slopes are extremely small,
no more than $k^{0.1}$ for all $|\alpMz| < 0.5$.

\begin{figure}[t]
\centering
\includegraphics[width=0.505\textwidth]{sp_diff_alpM0_neg}
\includegraphics[width=0.485\textwidth]{sp_diff_alpM0_pos}
%\includegraphics[width=0.496\textwidth]{spectra_alpM_choice0_error_WKB}
\caption{
Relative error $\varepsilon^\WKB (S_{h'})$ in the saturated spectra at late times of the WKB approximation compared to the numerical simulations.
The vertical lines correspond to the wave number
$\klimaM$, below which the WKB estimate is expected to break down.}
\label{fig:sp_error_WKB}
\end{figure}

Finally, to quantify the departures of WKB approximation from numerical solutions, especially at $k\leq\klimaM$,
we show the relative errors in the spectra as $\varepsilon^\WKB (S_{h'})\equiv [S_{h'} (k)^\WKB - S_{h'} (k)]/S_{h'} (k)$ in
\Fig{fig:sp_error_WKB}, calculated at the present time.
We observe that for positive values of $\alpMz$, 
the errors grow monotonically as $k$ becomes smaller.
In fact, $\varepsilon^\WKB(S_{h'})\sim k^{-2}$ is found due to 
a consistent disagreement between spectral slopes $k^0$ predicted by WKB formalism and $\sim k^{-2}$ obtained numerically.
On the other hand, for negative values of $\alpMz$,
$\varepsilon^\WKB(S_{h'})$ also increase as $k$ becomes smaller.
But since $\alpMz < 0$ still produces certain degree of IR enhancement in numerical results,
the slope in the error spectrum does not become as steep and depends on the corresponding values of $\alpMz$.
In both panels of \Fig{fig:sp_error_WKB}, however,
the errors are limited at $\lesssim10\%$ above $\klimaM$.

\section{Observational implications}
\label{sec:obs_implications}

To infer the observational prospects of detecting a modified GWB, 
we convert the linear energy spectrum $S_{h'}(k)$ directly obtained from the code to the commonly used logarithmic energy spectrum $\OmGW(f)$ via \Eq{eqn:OmGW_Sh}.
The normalized wave numbers used in the code are converted to the
present-day physical frequencies via
\begin{equation}
f = \frac{k\hhh_*}{2\pi}\Big(\frac{a_*}{a_0}\Big).
\end{equation}
As a result, the GR solution in the form of the double broken power law used in the numerical simulations,
i.e., $S_{h'}(k)\propto k^2$ and $S_{h'}(k)\propto k^{-2}$ in the low- and high wave number regimes,
becomes $\OmGW(f)\propto f^3$ in the IR and $\OmGW(f)\propto f^{-1}$ in the UV regimes, respectively.

Choosing the cases with the most pronounced modifications,
we show in \Fig{fig:OmGW_obs} the results from series M0, 
%AB: figure ref not resolved
where the solid black curves show GR solutions and colored ones indicate $\alpMz = -0.5$ (blue), $-0.3$ (green), $-0.1$ (orange), and $-0.01$ (purple).
Assuming QCDPT as the initial time of GW production, the current-day $\OmGW(f)$ peaks around the nHz frequency band,
shown by the fainter set of colored curves.
If instead the initial time is around EWPT, then $\OmGW(f)$ today peaks closer to the mHz band,
illustrated by the curves in the same but darker colors.
For these, we extend further into the IR and UV regimes with dashed
lines with slopes $f^1$ and $f^{-8/3}$, respectively.
As a $f^2$ spectral enhancement upon the $f^3$ slope of the GR solution in the IR, 
all curves corresponding to $\alpMz < 0$ already exhibit $f^1$ towards their low-frequency tails.
On the other hand, we extend the UV tail with a slope of $f^{-8/3}$, in accordance to what is commonly found for turbulence-induced GWs.

%YH: added figure below (to be finalized after discussing the observation section here)
%YH: just some notes so not to forget later. given k_crit is generally below k_lim, maybe it's not a good idea to generate spectra only from WKB like below
%YH: would PTA be sensitive enough to see the tiny differences between k_break = 1 and k_peak = 10 due to different alpM?
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{OmGW_obs_WKB}
\caption{The present-day modified GW energy spectrum for parameter choice 0. 
Different colors of the shades represent different initial peak amplitude.
The shaded regions are bounded by $\alpMz = -0.01$ and $-0.1$.
Several detectors' power-law sensitivity (PLS) curves are shown for comparison: 
SKA, LISA, DECIGO, Big Bang Observatory (BBO), Einstein Telescope (ET), and Cosmic Explorer (CE) are in gray,
and NANOGrav and LVK are in light blue.
All PLS assumes a signal-to-noise ratio of 10 and an observation time of 4 years.}
\label{fig:OmGW_obs_WKB}
\end{figure}
%YH: commented out the old figures
%\begin{figure}[t]
%\centering
%\includegraphics[scale = 0.485]{Figures/OmGW_obs.png}\\
%\includegraphics[scale = 0.485]{Figures/OmGW_obs_2.png}
%\caption{The present-day modified GW energy spectrum for parameter choices 0 and III with 
%$\alpMz = -0.01$ (purple), $-0.1$ (orange), and $-0.3$ (green). 
%The GR solution is shown in black.
%Several detectors' power-law sensitivity (PLS) curves are shown for comparison: 
%SKA, LISA, DECIGO, Big Bang Observatory (BBO), Einstein Telescope (ET), and Cosmic Explorer (CE) are in gray,
%and NANOGrav and LVK are in light blue.
%All PLS assumes a signal-to-noise ratio of 10 and an observation time of 4 years.
%Upper and lower panels assume GW production from conservative and optimistic values of initial MHD source energy $\EEM^*\approx10^{-3}$ and $\EEM^*\approx10^{-1}$, respectively.}
%\label{fig:OmGW_obs}
%\end{figure}
For observational comparisons, 
we plot in gray the power law sensitivity (PLS) curves of proposed future detectors such as SKA~\cite{Moore:2014lga}, DECIGO, BBO~\cite{Schmitz:2020syl}, LISA~\cite{Caprini:2019pxz,RoperPol:2021xnd}, Einstein Telescope (ET), and Cosmic Explorer (CE)~\cite{Schmitz:2020syl}.
For current detectors, we show in light blue the $2\sig$-confidence contour assuming a 30-frequency power law using the NANOGrav 12.5-year dataset~\cite{NANOGrav:2020bcs} and the PLS of the LIGO-Virgo-KAGRA (LVK)~\cite{Schmitz:2020syl} network.

There is a number of implications to unfold here.
%First, assuming QCDPT origin and a conservative amount of initial GWs
%AB: inserted "a" for origin
First, assuming a QCDPT origin and a conservative amount of initial GWs
at $\EEGW^*\approx10^{-9}$ (upper panel of \Fig{fig:OmGW_obs}),
$\alpMz\gtrsim-0.5$ can already be excluded as it overshoots the NANOGrav contour.
Smaller values such as $\alpMz = -0.3$ and $-0.1$ amplify the otherwise
%unobserved GR solution to the NANOGrav and the expected SKA range.
%AB: range -> ranges
unobserved GR solution to the NANOGrav and the expected SKA ranges.
%Similarly for GWs from EWPT epoch, $\alpMz\lesssim-0.3$ could boost
%AB: inserted "the" for epoch
Similarly for GWs from the EWPT epoch, $\alpMz\lesssim-0.3$ could boost
%the GR solution otherwise under the LISA and DECIGO curves to be
the GR solution, otherwise under the LISA and DECIGO curves, to be
above them.
Assuming a more optimistic value of $\EEGW^*\approx10^{-5}$ (lower panel of \Fig{fig:OmGW_obs}),
which could be due to a larger amount of $\EEM^*\approx10^{-1}$,
the constraints are correspondingly stronger.
In this case, NANOGrav could rule out $\alpMz\gtrsim-0.3$ for GWs from
%QCDPT epoch, whereas for GWs from EWPT, a smaller $\alpMz\lesssim-0.1$
%AB: inserted "the" for epoch
the QCDPT epoch, whereas for GWs from EWPT, a smaller $\alpMz\lesssim-0.1$
would be sufficient to boost a GR solution to intersect with LISA and
DECIGO detection ranges.
Although these amplification effects could in principle provide upper and
lower bounds of $|\alpMz|$ (for $\alpMz < 0$), the present-day amplitude
%of GW spectrum is degenerate with respect to the value of $\alpMz$
%AB: inserted "the" for spectrum
of the GW spectrum is degenerate with respect to the value of $\alpMz$
as well as the initial GW amplitude, the latter of which contains a fair
amount of uncertainty itself.
In addition to the GW amplitude, its spectral peaks and slopes are also crucial factors that affects the detectability.
The discussions here are done rather qualitatively, where we inspect the $f^{-8/3}$ behavior (dashed lines) beyond the peak that, in reality, depends on the nature of the turbulence source.
The peak itself, which is assumed to occur at $k_* = 10$ initially, also varies depending on the physics of turbulence after the phase transitions.
%Therefore in practice, constraining $|\alpMz|$ via its amplification
%AB: inserted comma
Therefore, in practice, constraining $|\alpMz|$ via its amplification
effects can be rather challenging.

On the other hand, $\alpMz$ also introduces a spectral enhancement such that $f^3\map f^1$ in the low-frequency end.
Assuming QCDPT as the GW production time, this feature lies below the peak around nHz and thus outside the detectors in that frequency range.
Assuming EWPT as the initial time is more interesting,
as the spectral enhancement could potentially lift the low-frequency tails up into the nHz range of detectors such as NANOGrav and SKA.
Combining the low-frequency spectral enhancement and the overall amplitude boost,
certain values of $\alpMz$ would be particularly interesting to explore,
where the modified GW spectrum could potentially be seen across nHz, mHz, and dHz detection ranges.
For example, it could happen for $\alpMz\lesssim-0.3$ if $\EEGW^*\approx10^{-9}$ or $\alpMz\lesssim-0.1$ if $\EEGW^*\approx10^{-5}$.
This possibility provides a way to reduce the degeneracy of parameters, especially with regards to the initial GW energy $\EEGW^*$.

Finally, we note that values such as $\alpMz = -0.3$ or $-0.5$ could potentially dominate the sensitivity budget of various detectors, mostly LISA and DECIGO but also SKA and NANOGrav,
especially if the GW amplitude $\EEGW^*$ is large to begin with.
If this actually happens, then the said detectors would have extreme difficulties extracting subdominant GW signals from other sources such as astrophysical objects,
as they do not receive the significant boost accumulated over longer periods of time as cosmological ones do.
Therefore, an upper bound could also be placed on $|\alpMz|$ in the event that GWs from closer-by astrophysical sources are identified by, for example, LISA.

\section{Conclusions}
\label{sec:conclusions}

GWs in modified gravity theories exhibit features in their
energy spectra different from what can be expected from GR.
We explored the significance of such features in terms of the spectral
slopes and amplitudes in different frequency ranges, under different
functional forms and values of two modification parameters -- the running
of effective Planck mass $\alpM$ and the tensor speed excess $\alpT$.
We also compared the analytical WKB approximation and numerical solutions
and found that the former is not always in agreement with the latter,
especially in the low-frequency regimes at late times.

The current-day value of $\alpT$ is already tightly constrained.
In this paper we further found the relatively insignificant modifications introduced by $\alpT$, even if taken to be unphysically large values compared to its constraint.
On the other hand, $\alpM$ acts as an (anti-)damping term in the GW equation 
and can impart changes to the energy spectrum of GWs produced by MHD sources 
at EWPT and QCDPT with potential implications to future observations, 
especially in the nHz and mHz frequency ranges.

In principle, the values of $\alpM$ could be constrained by the spectral
shape and amplitude of the GW energy spectrum, especially since
$\alpM < 0$ flattens the low-frequency tail of the spectrum by a
factor of $f^2$ and boosts the overall amplitude of the spectrum.
In practice, this is a rather complex problem with parameter degeneracies,
as the un-modified GW spectrum could also have larger or smaller
amplitude due to a range of possible source amplitudes to begin with,
in addition to the uncertainties of the GW spectral slopes that depend
on the physics of turbulence sources after the phase transitions.
Combining multiple detectors at different frequency ranges could help reducing the degeneracy.
Identifying astrophysical sources could also put constraints on $\alpM$.
But overall, we provide a rather qualitative discussion here and the
%AB: Is "rather qualitative" good? It sounds negative.
utilization of GW spectrum from early-universe sources as a probe of the
parameter $\alpM$ requires more detailed studies on various aspects of
the subject.

\paragraph{Data availability.} 
The source code used for the numerical solutions of this study,
the {\sc Pencil Code}, is freely available~\cite{pencil};
see also ref.~\cite{DATA} for the numerical data.

\acknowledgments
Support through the grant 2019-04234 from the Swedish Research Council
(Vetenskapsr{\aa}det) is gratefully acknowledged.
Nordita is sponsored by Nordforsk.
We acknowledge the allocation of computing resources provided by the
Swedish National Allocations Committee at the Center for Parallel
Computers at the Royal Institute of Technology in Stockholm.

\appendix

\section{Numerical scheme}
\label{sec:numerical_scheme}

Similar to what has been done in refs.~\cite{RoperPol+18,BS21},
we implement a matrix solver for the modified GW equation.
We start by simplifying the notations in \Eq{eqn:GW_mod} as
\begin{equation}
h_{\ij}'' + \sig h_{\ij}' + \om^2h_{\ij} = \ggg T_{\ij},
\label{eqn:GW_mod_rewrite}
\end{equation}
where $\ggg\equiv6/a(\eta)$, and
\begin{equation}
\sig\equiv\alpM\hhh\,,\,\,
\om^2\equiv\cT^2k^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}
<<<<<<< .mine
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.
||||||| .r48
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.
=======
>>>>>>> .r103
For $\del\eta\ll1$, 
<<<<<<< .mine
assume $T_{\ij}(\eta)\approx T_{\ij}(\eta + \del\eta)$,
||||||| .r48
assume $T_\ij(\eta)\approx T_\ij(\eta + \del\eta)$,
=======
we assume $T_\ij(\eta)\approx T_\ij(\eta + \del\eta)$,
which is justified when considering MHD sources as their characteristic time scales are typically much shorter than those of GWs,
>>>>>>> .r103
and $\lambda_{1,2}' \approx \lambda_{1,2}''
\approx 0$.
<<<<<<< .mine
Then \Eq{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
||||||| .r48
Then \Eq{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
=======
Then solutions for the strain quantities in \Eq{eqn:GW_mod_rewrite} take the form of
>>>>>>> .r103
\begin{align}
<<<<<<< .mine
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}.
||||||| .r48
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}.
=======
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},
>>>>>>> .r103
\end{align}
<<<<<<< .mine
Equivalently in matrix form for $h_{\ij}$ and $h_{\ij}'$,
||||||| .r48
Equivalently in matrix form for $h_\ij$ and $h_\ij'$,
=======
where $\ccc_\ij$ and $\ddd_\ij$ are constant amplitude coefficients evaluated at $\eta$.
Equivalently in matrix form for $h_\ij$ and $h_\ij'$,
>>>>>>> .r103
the solutions can be rewritten as
\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{Numerical accuracy}
\label{sec:numerical_accuracy}

\subsection*{Accuracy as a function of time intervals}

In this section we attach a few case studies on numerical accuracy.
Although the GW solver in the {\sc Pencil Code} has been made accurate to second order in time,
this is only true when the coefficients are constant in time.
When the coefficients vary only slowly, the error is also small, but depends on the length of the time step only to first order.
This is demonstrated in \Fig{fig:dEEGW_dn_incr},
where we show on the left panel the relative errors in the GW energy during RD and MD,
and on the right panel that during $\Lam$D.
Recall that during RD and MD, we use increasing time steps, represented by $n_\incr$ here.
During $\Lam$D, we revert back to linear time steps, with the interval shown as $\eta_\incr$ in \Fig{fig:dEEGW_dn_incr}.
\begin{figure}[t]
\centering
\includegraphics[scale = 0.48]{Figures/dEEGW_dn_incr.png}
\caption{Relative errors as a function of time intervals.
\textit{Left panel}: The accuracy during RD (round dots) and MD (square
dots) for $\alpMz = -0.5$ (blue) and $\alpMz = -0.01$ (green), where
$\del n_\incr$ is the difference in the time steps.
\textit{Right panel}: $\del\eta_\incr$ indicates the linear time interval
difference during $\Lam$D.}
\label{fig:dEEGW_dn_incr}
\end{figure}
Another way of seeing the time accuracy is via \Fig{fig:EGW_deta},
where we show the final energy spectrum (left panel) and its averaged results (right panel) for run M0A.
Different colors indicate the linear time intervals during $\Lam$D and
the black curves are the spectrum obtained by running M0A
entirely with increasing time steps using $n_\incr = 0.01$.
In other words, the colored curves only differ from the black one from the linear time stepping during $\Lam$D.
We observe that all chosen values of $\del\eta$ produce converging
energy spectra in the small-$k$ regime after averaging,
although larger values of $\del\eta$ result in more fluctuations
in the un-averaged spectra.
We also note that, compared to results from hybrid time steps, the
entirely nonuniform time steps underestimate the final energy spectrum.
\begin{figure}[t]
\centering
\includegraphics[scale = 0.48]{Figures/EGW_deta.png}
\caption{Final energy spectrum $\EGW(k)$ (left) and the averaged results
(right) for $\alpMz = -0.5$ (run M0A).
All runs are done with increasing time steps using $n_\incr = 0.01$ during
RD and MD but continue in linear time intervals (indicated by the colors)
during $\Lam$D.
The black dashed curves being the run M0A done entirely with increasing
time steps with $n_\incr = 0.01$.}
\label{fig:EGW_deta}
\end{figure}

%YH: suppress (and can probably discard eventually) the subsection below, as now it's discussed in the energy spectrum section in the main text
%\subsection*{Accuracy in the small-$k$ range}

%\begin{figure}[t]
%\centering
%\includegraphics[scale = 0.48]{Figures/EGW_dk1.png}
%\caption{Runs M0A (blue), M0C (orange), and M0D (purple) extended further into the small-$k$ regime.}
%\label{fig:EGW_dk1}
%\end{figure}

%Another aspect of numerical accuracy lies in the small-$k$ range of
%the energy spectra, where the most significant deviations from GR
%solutions occur.
%For this we re-run M0A, M0C, and M0D but with its wave numbers extended four orders of magnitude smaller.
%The resulting final energy spectra and the averaged values are shown
%in the left and right panels of \Fig{fig:EGW_dk1}.
%For $\alpMz = -0.01$ (purple), the spectrum remains flat for all wave numbers below $k_\crit$.
%For $\alpMz = -0.1$ (orange), below $k_\crit$, the spectrum starts out flat as $k^0$ but changes to $k^{0.1}$ as $k$ becomes smaller.
%This steepening of the modified (flat) spectra is seen more clearly
%with $\alpMz = -0.5$ (blue), where the slope first becomes $k^{0.4}$
%between $k = 10^{-3}$ and $k = 10^{-2}$, then $k^{0.6}$ across four
%decades between $k = 10^{-7}$ to around $k = 10^{-3}$.
%However, this only becomes important when considering the extreme IR regime,
%which, after converting to the physical frequency today, lies
%completely outside the lowest detection range around nHz considered in
%\Sec{sec:obs_implications}.

\bibliographystyle{JCAP}
\bibliography{references}

\end{document}
