# Modified propagation of gravitational waves from the early radiation era

Yutong He  Alberto Roper Pol  and Axel Brandenburg

###### Abstract

We study the propagation of cosmological gravitational wave (GW) backgrounds
from the early radiation era until the present day in modified theories of
gravity. Comparing to general relativity (GR), we study the effects that
modified gravity parameters, such as the GW friction $\alpha_{\rm M}$ and the
tensor speed excess $\alpha_{\rm T}$, have on the present-day GW spectrum. We
use both the WKB estimate, which provides an analytical description but fails
at superhorizon scales, and numerical simulations that allow us to go beyond
the WKB approximation. We show that a constant $\alpha_{\rm T}$ makes
relatively insignificant changes to the GR solution, especially taking into
account the constraints on its value from GW observations by the LIGO–Virgo
collaboration, while $\alpha_{\rm M}$ can introduce modifications to the
spectral slopes of the GW energy spectrum in the low-frequency regime
depending on the considered time evolution of $\alpha_{\rm M}$. The latter
effect is additional to the damping or growth occurring equally at all scales
that can be predicted by the WKB approximation. In light of the recent
observations by pulsar timing array (PTA) collaborations, and the potential
observations by future detectors such as SKA, LISA, DECIGO, BBO, or ET, we
show that, in most of the cases, constraints cannot be placed on the effects
of $\alpha_{\rm M}$ and the initial GW energy density $\mathcal{E}_{\rm
GW}^{*}$ separately, but only on the combined effects of the two, unless the
signal is observed at different frequency ranges. In particular, we provide
some constraints on the combined effects from the reported PTA observations.

##  1 Introduction

The present day Hubble constant, $H_{0}$, is measured to be around $H_{0}\sim
74\,{\rm km/s/Mpc}$ by astrophysical tests using type Ia supernovae [1, 2, 3],
lensed quasars [4], and megamaser-hosting galaxies [5]. On the other hand,
cosmological tests of $H_{0}$ constrain its value to $H_{0}\sim 67\,{\rm
km/s/Mpc}$ from cosmic microwave background (CMB) [6] and baryon acoustic
oscillations (BAO) [7] experiments, assuming the standard $\Lambda$CDM
($\Lambda$ cold dark matter) model of cosmology. These measurements suggest
the presence of a $4$–$5\,\sigma$ deviation discrepancy between the early- and
late-universe measurements of $H_{0}$ [8], known as the Hubble tension. The
observation of gravitational waves (GWs) from “standard sirens” provides an
independent measurement of the Hubble rate [9, 10]. Current observations from
LIGO–Virgo present large uncertainties and cannot yet solve the tension [11,
12, 13, 14]. However, the measurement of $H_{0}$ is expected to drastically
improve with the space-based GW detector LISA [15] or the next-generation
ground-based detectors like the Einstein Telescope (ET) [16]. Besides the
well-known $H_{0}$ tension, a number of other observational discrepancies
within the $\Lambda$CDM model have been reported, such as the $S_{8}$ tension,
where $S_{8}\propto\sigma_{8}\,\Omega_{{\rm mat},0}^{1/2}$ characterizes the
structure growth rate $\sigma_{8}$ and today’s matter density $\Omega_{{\rm
mat},0}$ [17]. A recent summary of over a dozen cosmological tensions and
their varied significance can be found in ref. [18].

In view of these tensions, $\Lambda$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 $\Lambda$CDM, modifying $\Lambda$CDM often requires
modifying GR.

In search for modified gravity (MG) theories, certain systematic approaches
have been developed to go beyond GR, such as relaxing one or some of the
assumptions made in Lovelock’s theorem that proves the uniqueness of GR [19,
20]. Overarching families of MG theories have been constructed to represent a
vast number of individual models, for example Horndeski’s theory, which
corresponds to all four-dimensional scalar-tensor models with second-order
derivatives [21, 22, 23], beyond Horndeski [24, 25, 26], Degenerate Higher-
Order Scalar-Tensor (DHOST) [27, 28, 29], Hořava-Lifshitz [30, 31],
Generalized Einstein-Aether gravity [32, 33], standard [34, 35] and beyond
[36, 37, 38] generalized Proca, bimetric gravity [39, 40, 41], and tensor-
scalar-vector (TeVeS) [42] theories. Some specific examples of MG theories
are: Brans-Dicke [43], quintessence [44], $f(R)$ [45], $f(G)$ [46],
$k$-essence [47], kinetic gravity braiding [48], and galileon [49]. For a
comprehensive review see, e.g., ref. [50].

At the level of tensor linear perturbations, MG theories introduce three
parameters to the standard gravitational wave (GW) equation that can be
conveniently used as a phenomenological probe in the parameter space of MG:

  * •

GW friction $\alpha_{\rm M}$, usually111This is not necessarily the only
interpretation of the GW friction; see, e.g., ref. [51]. associated to the
running of the effective Planck mass $\alpha_{\rm M}={\rm\,d}\ln M_{\rm
eff}^{2}/{\rm\,d}\ln a$. The friction term yields a GW luminosity distance
different than the electromagnetic (EM) counterpart $d_{\rm GW}\neq d_{\rm
EM}$ that allows us to put constraints on $|\alpha_{\rm
M}|\sim\mathcal{O}(10^{0})$ [52, 53, 54, 55, 56, 51]. Signals propagating over
longer distances could provide tighter constraints. Such “standard siren”
measurements have so far yielded mild constraints on $\alpha_{\rm M}$, as
LIGO–Virgo covers a relatively small distance of $z\sim\mathcal{O}(10^{-1})$
[57, 58]. Cosmological tests involving the CMB provide a tighter constraint of
$|\alpha_{\rm M}|\lesssim\mathcal{O}(10^{-1})$ [6] at the present day; see
section 4 for details. This aspect can also be improved with LISA [59, 60, 15]
and ET [56, 51, 61, 62, 63, 16], which might probe distances up to
$z\sim\mathcal{O}(10^{1})$.

  * •

Tensor speed excess $\alpha_{\rm T}$, related to the GW speed $c_{\rm T}$ via
$c_{\rm T}^{2}=\alpha_{\rm T}+1$. It has been tightly bounded to $\alpha_{\rm
T}\lesssim\mathcal{O}(10^{-15})$ by the multi-messenger detection of GW170817
and GRB 170817A [64]. This observation in turn constrains the theory space of
MG [65, 66, 67, 68]. However, this constraint only applies to the frequency
range probed by LIGO–Virgo and might not necessarily apply at different
frequencies, as proposed in ref. [69].

  * •

Graviton mass $m_{\rm g}$, which modifies the otherwise linear dispersion
relation in GR.

In this paper, we focus on the study of the modifications on the expected
relic cosmological GW spectrum under MG for different time-dependent GW
friction models, and for constant values of $\alpha_{\rm T}\neq 0$. We present
the general theoretical framework for time-dependent $\alpha_{\rm M}$ and
$\alpha_{\rm T}$, but defer the study of time and/or frequency dependent
$\alpha_{\rm T}$ to future work. For a similar study concerning the effects of
$m_{\rm g}$, see ref. [70]. Our focus on the cosmological effects of
$\alpha_{\rm M}$ is motivated by its less stringent constraints compared to
$\alpha_{\rm T}$, and its accumulated effects over cosmic history. We ignore
the effects of scalar perturbations in MG that can lead to gravitational slip
and focus only on the phenomenological effects on the tensor sector. For the
combined study of scalar and tensor perturbations in MG see refs. [53, 63].
Similarly, we ignore modifications on the dark energy equation of state and
set $w_{\rm DE}=-1$ (i.e., $w_{0}=-1$ and $w_{a}=0$) [71, 6], such that the
$\Lambda$CDM background expansion is unmodified.

To this end, we initialize a GW spectrum based on previous studies of GWs
sourced by primordial magnetohydrodynamic (MHD) turbulent fields in the early
radiation era—e.g., at the electroweak or QCD phase transitions (EWPT or
QCDPT)—and then propagate it through cosmic history to obtain the present-day
relic spectrum. Actually, we take a generic double broken power law GW
spectrum and we show that the obtained results are independent of the initial
spectral shape. Hence, our results can be applied to predict the expected GW
spectrum from a generic source in MG. We then analytically and numerically
compare the relic spectra in GR and MG, and discuss their potential
observational implications. The numerical solutions in this study are obtained
using the Pencil Code [72], which has been used as a tool for simulations of
GWs from primordial turbulent sources since the implementation of a GW solver
[73]. The cosmological history of the universe is modelled by numerically
solving Friedmann equations. We provide a cosmological solver and a tutorial
within the associated online material [74]. It also contains the routines that
have been used to process the data obtained from the Pencil Code and to
generate the results and plots of the present work.

We introduce the propagation of tensor-mode perturbations described by the GW
equation in MG in section 2. In section 3, we introduce the WKB approximation
commonly used to describe approximate solutions to the modified GW equation.
We derive the WKB solution in section 3.1, study its limitations and range of
validity in section 3.2, and compute the GW spectrum in section 3.3. Then, in
section 4, we present some common parameterizations of the parameter
$\alpha_{\rm M}$ through the 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 Pencil Code in section 5. In particular, we
study the effects on the GW spectrum that can not be predicted using the WKB
In particular, we study the effects on the GW spectrum that cannot be
predicted using the WKB approximation. Finally, we discuss potential
observational implications of MG compared to GR for the different
parameterizations of $\alpha_{\rm M}$ in section 6 and conclude in section 7.

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

##  2 Modified GWs on FLRW background

For simplicity, we assume the background evolution to follow $\Lambda$CDM,
while allowing a generic modification of the GW evolution. Therefore, we start
with the homogeneous and isotropic background described by the Friedmann-
Lemaître-Robertson-Walker (FLRW) metric. Including tensor perturbations, the
FLRW line element reads

${\rm\,d}s^{2}=a^{2}\bigl{(}-{\rm\,d}\eta^{2}+[\delta_{ij}+a^{-1}h_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle
x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle
x$}}{\mbox{\boldmath$\scriptscriptstyle
x$}}{},\eta)]{\rm\,d}x^{i}{\rm\,d}x^{j}\bigr{)},$ (2.1)

where $h_{ij}=ah_{ij}^{\rm phys}$ are the strains obtained by scaling the
physical strains $h_{ij}^{\rm phys}$ with the scale factor $a$. In GR, the GW
equation in Fourier space222We use a tilde for variables in Fourier space and
the following convention  $\tilde{h}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}{})=\int
h(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle
x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle
x$}}{})\,e^{-i\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{}\cdot\mathchoice{\mbox{\boldmath$\displaystyle
x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle
x$}}{\mbox{\boldmath$\scriptscriptstyle
x$}}{}}{\rm\,d}^{3}\mathchoice{\mbox{\boldmath$\displaystyle
x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle
x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},\quad
h(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle
x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle
x$}}{})=\frac{1}{(2\pi)^{3}}\int\tilde{h}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{})\,e^{i\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{}\cdot\mathchoice{\mbox{\boldmath$\displaystyle
x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle
x$}}{\mbox{\boldmath$\scriptscriptstyle
x$}}{}}{\rm\,d}^{3}\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}{}.$ reads [75, 73]

${\tilde{h}}_{ij}^{\prime\prime}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)+\biggl{(}k^{2}-\frac{a^{\prime\prime}}{a}\biggr{)}\tilde{h}_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)=\frac{6}{a}\tilde{T}_{ij}^{\rm
TT}(\mathchoice{\mbox{\boldmath$\displaystyle k$}}{\mbox{\boldmath$\textstyle
k$}}{\mbox{\boldmath$\scriptstyle k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta),$ (2.2)

where $\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 ${\cal E}_{\rm rad}=3H^{2}/\kappa$,
being $H\equiv\dot{a}/a$ the Hubble rate at time $\eta$. The projection
operator is $\Lambda_{ijlm}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}}{})=P_{il}P_{jm}-{\textstyle{1\over 2}}P_{ij}P_{lm}$, where
$P_{ij}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}}{})=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$ and
$\hat{\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}}{}=\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}{}/k$. The scale factor at
$\eta_{*}$ is set to unity, and the wave numbers $k$ are also normalized by
the (normalized) conformal Hubble rate ${\cal H}\equiv a^{\prime}/a$. We use
an asterisk to refer to the time in the early universe at which the GWs were
generated.

Including the additional parameters $\alpha_{\rm M}$ and $\alpha_{\rm T}$
mentioned in section 1, equation (2.2) generalizes to [76, 77]

$\tilde{h}_{ij}^{\prime\prime}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}{},\eta)+\alpha_{\rm M}\,{\cal
H}\,\tilde{h}_{ij}^{\prime}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}{},\eta)+\biggl{(}c_{\rm
T}^{2}k^{2}-\alpha_{\rm M}{\cal
H}^{2}-\frac{a^{\prime\prime}}{a}\biggr{)}\,\tilde{h}_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)=\frac{6}{a}\tilde{T}_{ij}^{\rm
TT}(\mathchoice{\mbox{\boldmath$\displaystyle k$}}{\mbox{\boldmath$\textstyle
k$}}{\mbox{\boldmath$\scriptstyle k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta).$ (2.3)

The effect of massive gravitons could be included by adding a $m_{\rm g}\neq
0$ term multiplying to $h_{ij}$ in equation (2.3), as done in ref. [70], but
we omit it in the present work.

Equation (2.3) can be recast in a simplified equation by introducing the new
variable $\chi_{ij}=e^{\mathcal{D}}h_{ij}$ [56, 51, 60], with

$\mathcal{D}(\eta)=\frac{1}{2}\int_{1}^{\eta{\cal H}_{*}}\alpha_{\rm
M}\,\mathcal{H}{\rm\,d}\tau,$ (2.4)

where the integral is performed over normalized time
$\tau=\eta^{\prime}\,{\cal H}_{*}$. This yields

$\tilde{\chi}_{ij}^{\prime\prime}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}{},\eta)+\biggl{(}c_{\rm
T}^{2}k^{2}-\frac{\tilde{a}^{\prime\prime}}{\tilde{a}}\biggr{)}\tilde{\chi}_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)=\frac{6}{a}\tilde{T}_{ij}^{\rm
TT}(\mathchoice{\mbox{\boldmath$\displaystyle k$}}{\mbox{\boldmath$\textstyle
k$}}{\mbox{\boldmath$\scriptstyle k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta),$ (2.5)

where $\tilde{a}=ae^{\mathcal{D}}$, and $\tilde{a}^{\prime\prime}/\tilde{a}$
can be assembled from $a$ and $\alpha_{\rm M}$ as

$\frac{\tilde{a}^{\prime\prime}}{\tilde{a}}=\frac{a^{\prime\prime}}{a}\biggl{(}1+\frac{1}{2}\alpha_{\rm
M}\biggr{)}+\frac{1}{2}\alpha_{\rm M}{\cal
H}^{2}\biggl{(}1+\frac{1}{2}\alpha_{\rm M}\biggr{)}+\frac{1}{2}\alpha_{\rm
M}^{\prime}{\cal H}.$ (2.6)

The variable $\chi$ is related to the difference between GW and EM luminosity
distances mentioned in section 1. We first express the latter in $\Lambda$CDM,

$d_{\rm
EM}=\frac{a_{0}}{aH_{0}}\int_{a}^{a_{0}}\frac{a_{0}{\rm\,d}a^{\prime}}{a^{\prime
2}\sqrt{\Omega(a)}},$ (2.7)

where $H_{0}$ is today’s Hubble rate and $\Omega(a)$ is the total energy
density, which includes matter, radiation, and dark energy contributions; see
appendix A. The GW luminosity distance is

$d_{\rm GW}=\frac{a}{\tilde{a}}d_{\rm EM}=e^{-\mathcal{D}}d_{\rm EM},$ (2.8)

which reduces to $d_{\rm GW}=d_{\rm EM}$ in GR; see ref. [78] for a detailed
calculation.

Previous numerical work was usually restricted to specific eras, e.g.,
radiation (RD) [73, 79, 80, 81, 82, 83, 84, 85] or matter domination (MD) [86,
87, 88], 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. The dynamical thermal history
of the universe during RD, represented by the relativistic $g_{*}$ and
adiabatic $g_{\rm S}$ degrees of freedom, was also ignored. In GR, this
formulation is justified since the evolution of the physical strains
$h_{ij}^{\rm 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 degrees of freedom and transitions between
radiation, matter, and dark energy dominations ($\Lambda$D) [89]. In the
present work, we focus on including such effects and solve the GW equation in
modified gravity, which present richer dynamics even when the source is
inactive, from the time of GW generation up to present time.

Assuming a piece-wise EOS, such that $w=1/3$, $0$, and $-1$ during RD, MD, and
$\Lambda$D, respectively, would create discontinuities in the time evolution
of $a(\eta)$ and its derivatives that appear in equations (2.3) and (2.5).
Therefore, to find a smooth $a(\eta)$, we directly solve the Friedmann
equations (see appendix A for details),

${\cal H}=\frac{H_{0}}{{\cal
H}_{*}}\,a\sqrt{\Omega(a)},\quad\frac{a^{\prime\prime}}{a}=\frac{1}{2}{\cal
H}^{2}\bigl{[}1-3w(a)\bigr{]},$ (2.9)

where $\Omega(a)$ denotes the total energy density, and the term ${\cal
H}_{*}$ appears due to our definition of the normalized ${\cal
H}=a^{\prime}/a$. We set $a_{*}=1$ for consistency with our GW equation; see
the discussion below equation (2.2). Using equation (2.9), we solve for the
time evolution of the scale factor $a(\eta)$ and finally obtain the time
evolution of ${\cal H}$ and $a^{\prime\prime}$ that are required to solve the
GW equation, either when they appear explicitly in equation (2.3) for $h_{ij}$
or to compute the term $\tilde{a}^{\prime\prime}$, defined in equation (2.6),
that appears in equation (2.5) for $\chi_{ij}$.

##  3 WKB approximation

###  3.1 Solution of the GW equation

Modifications of the GW propagation, in the absence of sources333In the
absence of sources, the GW equation [see equation (2.3)] does not depend on
the wave vector $\textstyle k$ 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 [55, 59]. The WKB
solution can be obtained using the following ansatz

$\tilde{\chi}_{ij}(k,\eta)=\chi(k,\eta)\,e_{ij}=Ae_{ij}\,e^{iB},$ (3.1)

where $A$ and $B$ are generic coefficients, and $e_{ij}$ is the polarization
tensor.444 In the absence of sources, the GW propagation of any polarization
mode is the same, so we can just call $\chi$ 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 equation (2.5), one
gets

$\displaystyle 0$
$\displaystyle=2\frac{A^{\prime}}{A}+\frac{B^{\prime\prime}}{B^{\prime}}\quad\quad\quad\quad\quad\quad\Rightarrow
A(\eta)=e^{-\tilde{\mathcal{D}}},$ (3.2) $\displaystyle B^{\prime 2}$
$\displaystyle=\frac{A^{\prime\prime}}{A}+k^{2}c_{\rm
T}^{2}-\frac{\tilde{a}^{\prime\prime}}{\tilde{a}}\approx k^{2}c_{\rm
T}^{2}\Rightarrow B(k,\eta)=\pm k(\eta\mathcal{H}_{*}-\Delta T),$ (3.3)

where we have neglected $A^{\prime\prime}/A$ and
$\tilde{a}^{\prime\prime}/\tilde{a}$ compared to $k^{2}c_{\rm T}^{2}$,
following the WKB approximation. We have defined the additional damping factor
$\tilde{\mathcal{D}}$, and the time delay due to the effective GW speed
$\Delta T$

$\tilde{\mathcal{D}}=\frac{1}{2}\int_{1}^{\eta\mathcal{H}_{*}}\frac{B^{\prime\prime}}{B^{\prime}}{\rm\,d}\tau=\frac{1}{2}\int_{1}^{\eta\mathcal{H}_{*}}\frac{c_{\rm
T}^{\prime}}{c_{\rm T}}{\rm\,d}\tau,\quad\Delta
T=\int_{1}^{\eta\mathcal{H}_{*}}(1-c_{\rm T}){\rm\,d}\tau.$ (3.4)

Hence, the particular solution is

$\chi(k,\eta)=e^{-\tilde{\mathcal{D}}}e^{\pm ik(\eta\mathcal{H}_{*}-\Delta
T)}=e^{-\tilde{\mathcal{D}}\mp ik\Delta T}\chi^{\rm GR}(k,\eta).$ (3.5)

where $\tilde{\mathcal{D}}=\Delta T=0$ in GR. For the initial conditions555In
the current work, we focus on the propagation of a GW background after it has
already been generated and reached a stationary solution, when 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{\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}}{}$, for each polarization mode.
$\chi(k,\eta_{*})=\chi_{*}(k)$ and
$\chi^{\prime}(k,\eta_{*})=\chi_{*}^{\prime}(k)$, the general solution reads

$\displaystyle\chi(k,\eta)$
$\displaystyle=e^{-\tilde{\mathcal{D}}}\Bigg{[}\chi_{*}(k)\cos k\tilde{c}_{\rm
T}(\eta\mathcal{H}_{*}-1)+\frac{c_{\rm
T}^{*}\chi_{*}^{\prime}(k)+{\textstyle{1\over 2}}{c_{\rm
T}^{\prime}}^{*}\,\chi_{*}(k)}{k{c_{\rm T}^{*}}^{2}}\sin k\tilde{c}_{\rm
T}(\eta\mathcal{H}_{*}-1)\Bigg{]},$ (3.6)

where $\tilde{c}_{\rm T}(\eta)=\int c_{\rm T}{\rm\,d}\eta^{\prime}/(\eta{\cal
H}_{*}-1)$, and $c_{\rm T}^{*}$ and ${c_{\rm T}^{\prime}}^{*}$ are the values
of the GW speed (and its derivative) at the initial time of GW production.
Note that the specific choice of initial conditions does not allow to give a
linear relation between $\chi$ and $\chi^{\rm GR}$ as done in ref. [55] and,
in general, $\chi\neq e^{-\tilde{\mathcal{D}}-ik\Delta T}\chi^{\rm GR}$, which
is only true when referring to the particular solution of the ODE, given in
equation (3.5). In terms of the initial conditions
$h(k,\eta_{*})=h_{*}(k)=\chi_{*}(k)$ and
$h^{\prime}(k,\eta_{*})=h^{\prime}_{*}(k)=\chi^{\prime}_{*}(k)-{\textstyle{1\over
2}}\alpha_{\rm M}^{*}\,\chi_{*}(k)$, the solution (3.6) can be rewritten as

$\displaystyle\chi(k,\eta)=e^{-\tilde{\mathcal{D}}}\biggl{[}$
$\displaystyle\,h_{*}(k)\cos k\tilde{c}_{\rm T}(\eta{\cal H}_{*}-1)$
$\displaystyle\,+\frac{c_{\rm T}^{*}h_{*}^{\prime}(k)+{\textstyle{1\over
2}}\bigl{(}c_{\rm T}^{*}\,\alpha_{\rm M}^{*}+{c_{\rm
T}^{\prime}}^{*}\bigr{)}h_{*}(k)}{k{c_{\rm T}^{*}}^{2}}\sin k\tilde{c}_{\rm
T}(\eta{\cal H}_{*}-1)\biggr{]}.$ (3.7)

The general solution reduces to the following when $c_{\rm T}=1$ at all times,

$\chi(k,\eta)=h_{*}(k)\cos
k(\eta\mathcal{H}_{*}-1)+\frac{h_{*}^{\prime}(k)+{\textstyle{1\over
2}}\alpha_{\rm M}^{*}\,h_{*}(k)}{k}\sin k(\eta\mathcal{H}_{*}-1).$ (3.8)

This expression holds both in GR with $h(k,\eta)=\chi(k,\eta)$ and
$\alpha_{\rm M}=0$, and in MG theories with
$h(k,\eta)=e^{-\mathcal{D}}\chi(k,\eta)$ and $\alpha_{\rm M}\neq 0$.

###  3.2 Limitations of the WKB approximation

To obtain the WKB solution, we have neglected two terms in equation (3.3),
which can be expressed as

$\displaystyle\frac{A^{\prime\prime}}{A}-\frac{\tilde{a}^{\prime\prime}}{\tilde{a}}$
$\displaystyle=-\frac{c_{\rm T}^{\prime\prime}}{2c_{\rm T}}+\frac{3c_{\rm
T}^{\prime 2}}{4c_{\rm T}^{2}}-\frac{1}{2}\alpha_{\rm
M}\Big{(}1+\frac{1}{2}\alpha_{\rm
M}\Big{)}\mathcal{H}^{2}-\frac{1}{2}\alpha_{\rm
M}^{\prime}\mathcal{H}-\frac{a^{\prime\prime}}{a}\Big{(}1+\frac{1}{2}\alpha_{\rm
M}\Big{)},$ (3.9)

where we have used ${\cal H}^{\prime}=a^{\prime\prime}/a-{\cal H}^{2}$. Hence,
the WKB assumption could break down when at least one of these terms is not
negligible when compared to $k^{2}c_{\rm T}^{2}$. First, we obtain the
limiting values of $k$ due to $c_{\rm T}$ and its derivatives in the first two
terms on the right-hand side of equation (3.9)

$k_{\rm lim,c_{\rm T}^{\prime\prime}}\lesssim\sqrt{\frac{|c_{\rm
T}^{\prime\prime}|}{2c_{\rm T}^{3}}},\quad k_{\rm lim,c_{\rm
T}^{\prime}}\lesssim\frac{\sqrt{3}|c_{\rm T}^{\prime}|}{2c_{\rm T}^{2}}.$
(3.10)

Note that in this paper we only consider $c_{\rm T}$ to be constant in time,
i.e., $k_{\rm lim,c_{\rm T}^{\prime\prime}}$ and $k_{\rm lim,c_{\rm
T}^{\prime}}$ are zero. We defer the study of time-varying $c_{\rm T}$ and its
effects to future work.

In addition, we have one term that depends on $\alpha_{\rm M}$, one term that
depends on the time evolution of $\alpha_{\rm M}$, and one term that depends
on $a^{\prime\prime}/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 equation (3.8). This term is given in equation (2.9) and its upper
bound can be found using Friedmann equations

$\frac{a^{\prime\prime}}{a}\lesssim\frac{1}{2}\frac{\Omega_{{\rm
mat},0}}{\Omega_{{\rm rad},0}}\frac{g_{*}}{g_{*}^{0}}\biggl{(}\frac{g_{\rm
S}}{g_{\rm S}^{0}}\biggr{)}^{-{4\over 3}}{\cal H}\,\frac{a_{*}}{a_{0}}\lesssim
10^{4}\,h^{2}\,{\cal H}\,\frac{a_{*}}{a_{0}},$ (3.11)

where the first inequality is an approximation 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 limiting wave numbers (for
$h=0.67$), at which the WKB approximation might break down,

$\displaystyle\,k_{{\rm lim},\,a^{\prime\prime}}^{\rm EW}\sim 1.87\times
10^{-6}\,{\cal H}\,\sqrt{\frac{\bigl{|}1+{\textstyle{1\over 2}}\alpha_{\rm
M}\bigr{|}}{c_{\rm T}}}\leq 1.87\times
10^{-6}\,\sqrt{\frac{\bigl{|}1+{\textstyle{1\over 2}}\alpha_{\rm
M}\bigr{|}}{c_{\rm T}}},$ (3.12) $\displaystyle\,k_{{\rm
lim},\,a^{\prime\prime}}^{\rm QCD}\sim 5.90\times 10^{-5}\,{\cal
H}\,\sqrt{\frac{\bigl{|}1+{\textstyle{1\over 2}}\alpha_{\rm M}\bigr{|}}{c_{\rm
T}}}\leq 5.90\times 10^{-5}\,\sqrt{\frac{\bigl{|}1+{\textstyle{1\over
2}}\alpha_{\rm M}\bigr{|}}{c_{\rm T}}}.$ (3.13)

In both cases, this term is subdominant up to very large superhorizon scales
and it is bounded by the values in equations (3.12) and (3.13) since ${\cal
H}=\eta_{*}/\eta$ during RD era as can be seen using equation (2.9). This
limit can be modified by the inclusion of $\alpha_{\rm M}$ and $\alpha_{\rm
T}$ 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 $\alpha_{\rm M}$ parameter,

$\displaystyle k_{{\rm lim},\,\alpha_{\rm M}}$ $\displaystyle\,\sim\frac{{\cal
H}}{c_{\rm T}}\sqrt{\biggl{|}\frac{\alpha_{\rm
M}}{2}\Bigl{(}1+{\textstyle{1\over 2}}\alpha_{\rm
M}\Bigr{)}\biggr{|}}\lesssim\frac{1}{c_{\rm T}}\sqrt{\frac{|\alpha_{\rm
M}^{*}|}{2}},$ (3.14) $\displaystyle k_{{\rm lim},\,\alpha_{\rm M}^{\prime}}$
$\displaystyle\,\sim\frac{1}{c_{\rm T}}\sqrt{\frac{{\cal H}\,|\alpha_{\rm
M}^{\prime}|}{2}}\leq\frac{1}{c_{\rm
T}}\sqrt{\frac{|\alpha^{{}^{\prime}*}_{\rm M}|}{2}},$ (3.15)

where we can neglect the term ${\textstyle{1\over 2}}\alpha_{\rm M}$ in front
of 1 in the first limit for small values of $\alpha_{\rm M}$. Hence, the WKB
limit can break down in MG around the horizon or at larger scales. In general,
we expect the limit from $\alpha_{\rm M}$ to be more restrictive than that
from $\alpha_{\rm M}^{\prime}$ and to dominate at the initial time when ${\cal
H}=1$. However, this depends on the parameterization of $\alpha_{\rm M}$,
which can give different results for both limits; see section 4.

###  3.3 GW spectrum using the WKB approximation

The spectrum of GW energy density can be expressed as

$\Omega_{\rm GW}(k,\eta)=\frac{1}{\rho_{{\rm crit},0}}\frac{{\rm\,d}\rho_{\rm
GW}}{{\rm\,d}\ln
k}=\frac{1}{6}\biggl{(}\frac{H_{*}}{H_{0}}\biggr{)}^{2}\biggl{(}\frac{a_{*}}{a_{0}}\biggr{)}^{4}kS_{h^{\prime}}(k,\eta),$
(3.16)

where $2S_{h^{\prime}}$ is the spectrum666 Following ref. [73], we define the
spectrum $S_{h^{\prime}}(k,\eta)$ from the $+$ and $\times$ polarization
modes, giving an extra factor of 2 due to the property
$\tilde{h}^{\prime}_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)\,\tilde{h}^{{}^{\prime}*}_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{}^{\prime},\eta)=2\bigl{[}h^{\prime}_{+}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)\,h^{{}^{\prime}*}_{+}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{}^{\prime},\eta)+h^{\prime}_{\times}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)\,h^{{}^{\prime}*}_{\times}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle k$}}{}^{\prime},\eta)\bigr{]}.$ of the
scaled strain derivatives,

$\langle\tilde{h}_{ij}^{\prime}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{},\eta)\,\tilde{h}_{ij}^{{}^{\prime}*}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{}^{\prime},\eta)\rangle=(2\pi)^{6}\delta^{3}(\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{}-\mathchoice{\mbox{\boldmath$\displaystyle
k$}}{\mbox{\boldmath$\textstyle k$}}{\mbox{\boldmath$\scriptstyle
k$}}{\mbox{\boldmath$\scriptscriptstyle
k$}}{}^{\prime})\frac{2S_{h^{\prime}}(k,\eta)}{4\pi k^{2}}.$ (3.17)

Taking $\mathcal{D}^{\prime}=\frac{1}{2}\alpha_{\rm M}\mathcal{H}$, the
modified GW energy spectrum can be obtained using the WKB solution for MG,
given in equation (3.6), and assuming that at time $\eta_{*}$ the GWs are in
the free-propagation regime following GR, such that
$|h^{\prime}_{*}(k)|=k\,|h_{*}(k)|$,

$\displaystyle\Omega_{\rm GW}$
$\displaystyle\,(k,\eta)=e^{-2\mathcal{D}}\Omega_{\rm GW}^{\rm
GR}(k,\eta)\,\xi(k,\eta),$ (3.18)

where $\xi(k,\eta)$ corresponds to the spectral modifications with respect to
the GR spectrum after taking into account the amplification (or depletion)
produced by the damping factor $\mathcal{D}$ and assuming constant-in-time
$c_{\rm T}$ (see appendix B for details and, in particular, equation (B.9) for
generic time-dependent $c_{\rm T}$),

$\xi(k,\eta)=1+\frac{1}{2}\alpha_{\rm T}+\frac{\alpha_{\rm
M}^{*}}{2k}+\frac{1}{8k^{2}}\Biggl{[}\alpha_{\rm M}^{2}{\cal
H}^{2}\biggl{(}1+\frac{1}{c_{\rm T}^{2}}\biggr{)}+\alpha_{\rm
M}^{*2}\Biggr{]}+\frac{\alpha_{\rm M}^{*}\alpha_{\rm M}^{2}{\cal
H}^{2}}{8k^{3}c_{\rm T}^{2}}+\frac{\alpha_{\rm M}^{*2}\alpha_{\rm M}^{2}{\cal
H}^{2}}{32k^{4}c_{\rm T}^{2}}.$ (3.19)

Note that we have averaged over oscillations in time $\eta$, such that the
expression is valid at $k\,\tilde{c}_{\rm T}\gtrsim(\eta{\cal H}_{*}-1)^{-1}$.
At late times, ${\cal H}^{2}=(\eta_{*}/\eta)^{2}$ becomes small, so we can
neglect the terms that dilute with ${\cal H}^{2}$,

$\xi(k,\eta\gg\eta_{*})={\textstyle{1\over 2}}\bigl{(}1+\alpha_{\rm
T}\bigr{)}+{\textstyle{1\over 2}}\Bigl{(}1+\frac{\alpha_{\rm
M}^{*}}{2k}\Bigr{)}^{2}.$ (3.20)

From equation (3.20), we can identify the IR and UV limiting ranges of the
spectrum,

$\displaystyle\Omega_{\rm GW}^{\rm IR}=\Omega_{\rm
GW}\Bigl{(}k\ll{\textstyle{1\over 2}}\alpha_{\rm
M}^{*},\eta\gg\eta_{*}\Bigr{)}=$ $\displaystyle\,\frac{{\alpha_{\rm
M}^{*}}^{2}}{8k^{2}}\,e^{-2\mathcal{D}}\,\Omega_{\rm GW}^{\rm GR}(k,\eta),$
(3.21) $\displaystyle\Omega_{\rm GW}^{\rm UV}=\Omega_{\rm
GW}\Bigl{(}k\gg{\textstyle{1\over 2}}\alpha_{\rm
M}^{*},\eta\gg\eta_{*}\Bigr{)}=$ $\displaystyle\,\Bigl{(}1+{\textstyle{1\over
2}}\alpha_{\rm T}\Bigr{)}\,e^{-2\mathcal{D}}\,\Omega_{\rm GW}^{\rm
GR}(k,\eta).$ (3.22)

Note that the factor ${\alpha_{\rm M}^{*}}^{2}$ is part of the IR limit,
indicating that the IR enhancement of the form $k^{-2}$ holds regardless of
the sign of $\alpha_{\rm M}^{*}$. The critical $k_{\rm
crit}={\textstyle{1\over 2}}|\alpha_{\rm M}^{*}|$ indicates where the IR
regime begins, i.e., where the $k^{-2}$ term in equation (3.20) becomes
dominant.

Equation (3.19) shows that the GW spectrum in the IR regime can present up to
$k^{-4}$ modifications to the GR spectrum but this and other terms vanish as
time evolves since they are proportional to ${\cal H}^{2}$. Hence, at late
times, we end up with the GR spectrum amplified (or diluted) for negative (or
positive) $\mathcal{D}$ by $e^{-2\mathcal{D}}$ at all wave numbers with a
$k^{-2}$ enhancement in the IR regime, proportional to the additional
${\textstyle{1\over 8}}{\alpha_{\rm M}^{*}}^{2}$ factor, as shown in equation
(3.21). In the UV regime, we find the $e^{-2\mathcal{D}}$ enhancement (or
depletion) and an additional factor $1+{\textstyle{1\over 2}}\alpha_{\rm T}$,
which is, in general, a function of $\eta$. The sign of $\alpha_{\rm T}$
determines if the GW spectrum is amplified or decreased with respect to that
obtained from GR.

We have found that, according to the WKB approximation, the parameter
$\alpha_{\rm M}$ introduces changes in the spectral shape at $k\leq k_{\rm
crit}={\textstyle{1\over 2}}|\alpha_{\rm M}^{*}|$ that do not dilute as time
evolves and depend on the value of $\alpha_{\rm M}$ only at the time of GW
generation. On the other hand, note that, using equations (3.14) and (3.15),
we can estimate the WKB approximation to break down around $k\sim\max(k_{{\rm
lim},\,\alpha_{\rm M}},k_{{\rm lim},\,\alpha_{\rm M}^{\prime}})$. This means
that when ${\cal H}\sim 1$, the spectral changes occur at wave numbers around
and below the critical $k_{\rm crit}$, where the WKB approximation might not
be valid (for $c_{\rm T}\sim 1$).

We investigate the resulting spectra in MG using numerical simulations in
section 5 with the objective to test the validity of the WKB approximation and
its potential limitations and to confirm the resulting GW spectra when MG
parameters are introduced.

##  4 Phenomenological parameterizations

In section 3.3, we noted that, according to the WKB approximation, a constant
$\alpha_{\rm T}\neq 0$ induces modifications in the total GW energy density,
but not in its spectral shape as long as $\alpha_{\rm T}\ll 1$. Realistically,
the tensor speed excess at the present day is observationally constrained to
be $\alpha_{\rm T,0}\lesssim\mathcal{O}(10^{-15})$ by the binary neutron star
merger GW170817 and its gamma-ray burst GRB 170817A [64]. This constraint can
be circumvented if $\alpha_{\rm T}$ is either frequency- or time-dependent,
such that larger deviations could hide outside the LIGO–Virgo frequency band
[69, 90] or in the past [91, 92]. In this study, we show that even using
larger values of $\alpha_{\rm T}$, 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: (i) constant $\alpha_{\rm T}$ with
$\alpha_{\rm M}=0$, and (ii) $\alpha_{\rm T}=0$ allowing $\alpha_{\rm M}$ to
take different values as a function of time.

In particular, we are interested in the perspectives of negative values of
$\alpha_{\rm M}$, as they would enhance the amplitude of the GW spectrum.
Various constraints on the present-day value $\alpha_{\rm M,0}$ exist in the
literature, but are significantly less stringent than those on $\alpha_{\rm
T}$. Assuming the cosmic acceleration due to scalar-tensor theories of
gravity, it can be shown that $|\alpha_{\rm M}|\lesssim 0.5$ at late times
(see figure 4 in ref. [54]). This is compatible with the numerical results of
ref. [59] (see their figure 2), giving $|\alpha_{\rm M,0}|\lesssim 1$. On the
other hand, the Planck Collaboration reports a slightly tighter constraint on
$\alpha_{\rm M,0}\gtrsim-0.1$ (see table 8 of ref. [6]). Note that they
parameterize $\alpha_{\rm M}=\alpha_{\rm M,0}\,(a/a_{0})^{n}$ with
$n\in(0.5,1)$. In addition, to avoid having an infinite number of GW sources
in the early Universe, $d_{\rm GW}$ is required to increase monotonically at
earlier times, i.e., ${\rm\,d}(d_{\rm GW})/{\rm\,d}a\leq 0$, as argued in ref.
[57]. Hence, using equations (2.7)–(2.8), one gets a bound for $\alpha_{\rm
M}$,

$\alpha_{\rm
M}(a)\geq\frac{a_{0}}{a\sqrt{\Omega(a)}}\Bigg{[}\int_{a}^{a_{0}}\frac{a_{0}\,da^{\prime}}{a^{\prime
2}\sqrt{\Omega(a^{\prime})}}\Bigg{]}^{-1}-1,$ (4.1)

which is satisfied by all positive values of $\alpha_{\rm M}$ and for negative
values with $|\alpha_{\rm M}|\leq|\alpha_{\rm M}(a_{*})|\lesssim 1$. Taking
these different constraints into account, we explore a range of $\alpha_{\rm
M,0}\in[-0.5,0.3]$.

Besides the simplest consideration of a constant $\alpha_{\rm M}$ in time,
which we call choice 0, its time-dependent forms can be written in accordance
with specific gravity models [93]. Following refs. [94, 77, 95], we choose
phenomenological forms of $\alpha_{\rm M}$ as simplified parameterizations
motivated by effective descriptions of modified gravity [96]. Specifically:

$\alpha_{\rm M}(\eta)=\begin{dcases}\alpha_{\rm
M,0}\biggl{[}\frac{a(\eta)}{a_{0}}\biggr{]}^{n}&(\text{choice I}),\\\
\alpha_{\rm M,0}\frac{1}{\Omega(\eta)}&(\text{choice II}),\\\ \alpha_{\rm
M,0}\frac{1-\Omega_{\rm mat}(\eta)/\Omega(\eta)}{1-\Omega_{{\rm
mat},0}}&(\text{choice III}),\end{dcases}$ (4.2)

where choices II and III give a value of $\alpha_{\rm M}$ 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 the present time. The exponent $n$ in choice I needs to satisfy the
following stability conditions [97]:

$\mbox{stability for:}\;\begin{dcases}0<n<{\textstyle{3\over 2}}\,\Omega_{{\rm
mat},0}\simeq{\textstyle{1\over 2}},\ &\mbox{if}\;\alpha_{\rm M,0}>0,\\\
n>{\textstyle{3\over 2}},\ &\mbox{if}\;\alpha_{\rm M,0}<0.\end{dcases}$ (4.3)

To provide an intuition on the aforementioned $\alpha_{\rm M}$
parameterizations, we show in figure 1 (left panel), for the different time
dependencies given in equation (4.2), the evolution of $\alpha_{\rm
M}\mathcal{H}$, which characterizes the growth or damping produced by
$\alpha_{\rm M}$, according to the WKB approximation [see equation (3.4)]. We
also show the time evolutions of $k_{{\rm lim},\,\alpha_{\rm M}}\,c_{\rm T}$
and $k_{{\rm lim},\,\alpha_{\rm M}^{\prime}}\,c_{\rm T}$ (right panel), which
correspond to the wave numbers at which the WKB approximation might break
down, according to equations (3.14) and (3.15). Figure 1 shows the time
evolution from the time of generation (e.g., EWPT) up to present time. We see
that for choices 0 and III, $\alpha_{\rm M}\mathcal{H}$ converges to the same
values during RD and $\Lambda$D, and only becomes different during MD, as
expected. On the other hand, for choices I and II, $\alpha_{\rm M}{\cal H}$ is
negligible for all of RD and most of MD and rapidly increases later on,
especially during $\Lambda$D. For all parameterizations, the GW friction
converges to the value $\alpha_{\rm M,0}$ at present time.

Figure 1:  Time evolution of $(\alpha_{\rm M}/\alpha_{\rm M,0})\mathcal{H}$
(left panel), which contributes to the change of amplitude over time,
$h_{ij}\sim e^{-\mathcal{D}}$ with $\mathcal{D}={\textstyle{1\over
2}}\int^{\eta}\alpha_{\rm M}{\cal H}{\rm\,d}\eta^{\prime}$, and of $k_{{\rm
lim},\,\alpha_{\rm M}}\,c_{\rm T}$ and $k_{{\rm lim},\,\alpha_{\rm
M}^{\prime}}\,c_{\rm T}$ (right panel), which are the terms neglected when
compared to $kc_{\rm T}$ under the WKB approximation; see equations (3.14) and
(3.15). 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 MD. 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.4\times 10^{13}$.

We see in figure 1 (right panel) that for choices 0 and III, the limiting $k$
is around the horizon scale at the time of generation, i.e., $k_{{\rm
lim},\,\alpha_{\rm M}}\,c_{\rm T}\sim\sqrt{\alpha_{\rm M}^{*}/2}$ and then it
decreases into superhorizon scales since it is proportional to ${\cal
H}=\eta_{*}/\eta$ (this is because for these choices, $\alpha_{\rm
M}\simeq\alpha_{\rm M}^{*}$ during the RD era). The term $\alpha_{\rm
M}^{\prime}=0$ for choice 0 and $\alpha_{\rm M}^{\prime}\ll\alpha_{\rm M}$ for
choice III, as can be seen in figure 1. The ratio between $\alpha_{\rm
M}^{\prime}$ and $\alpha_{\rm M}$ at the time of GW generation, within RD, can
be found to be (see the associated online material [74])

$\biggl{(}\frac{|\alpha_{\rm M}^{\prime}|}{\alpha_{\rm M}}\biggr{)}^{1\over
2}_{\rm III}\lesssim\biggl{[}\frac{\Omega_{{\rm mat},0}}{\Omega_{{\rm
rad},0}}\frac{g_{*}}{g_{*}^{0}}\biggl{(}\frac{g_{\rm S}}{g_{\rm
S}^{0}}\biggr{)}^{-{4\over 3}}\frac{a_{*}}{a_{0}}\biggr{]}^{1\over
2}\lesssim\sqrt{2}\times
10^{2}\,h\biggl{(}\frac{a_{*}}{a_{0}}\biggr{)}^{1\over 2},$ (4.4)

which only becomes of order $1$ when the time of GW generation is around the
end of the RD era (see right panel of figure 1). For choices I and II,
although the two terms are of the same order, ${\alpha}^{{}^{\prime}*}_{\rm
M}\sim\alpha_{\rm M}^{*}$, their value at the time of GW generation is much
smaller than $\alpha_{\rm M,0}$ (for all $n>0$ in the choice I). Therefore,
for choices 0 and III, the WKB estimate might break down around the horizon
and at superhorizon scales, while for choices I and II, the approximation is
estimated to be valid also at superhorizon scales (down to the limiting wave
number). Hence, for the latter, we do not expect to observe any relevant
spectral change, since the IR enhancement is determined by the value of
$\alpha_{\rm M}^{*}$, as seen in equation (3.20).

Note that, since many efforts of modifying gravity are aimed at addressing the
late-time acceleration of the universe, the $\alpha_{\rm M}$ 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
$\Lambda$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 $\alpha_{\rm M}\mathcal{H}$
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 MG during inflation [77,
98] and around recombination [99, 100], where $\alpha_{\rm M}$ is essentially
a free function of the scalar field, there seems to be a relative lack of
numerical studies on the effects of $\alpha_{\rm M}$ during RD (note, however,
a brief discussion of $\alpha_{\rm M}=-1$ during RD in ref. [101]). Hence, we
emphasize the aim of this work to provide an understanding of the
phenomenological behavior of the GW spectrum due to $\alpha_{\rm M}$ for
signals produced during RD, although the GW background computed in section 3.3
for MG is relevant in the most general case.

##  5 Numerical solutions

To explore the limits and validity of the WKB approximation, we use the Pencil
Code to numerically solve the GW equation under MG, given in equation (2.3).
The Pencil Code is a highly parallelized modular code that can be used to
solve various differential equations [72]. In the context of cosmological GWs,
it has previously been used to study GWs generated by hydrodynamic and MHD
stresses in the early universe [79]. 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 [73]. Previous
numerical work solved the GW equation under GR while, in the present work, we
have extended the code to solve equation (2.3); see appendix C for further
details on the numerical scheme.

###  5.1 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 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. [70,
84]). Hence, equation (2.3) is solved in one dimension (i.e., in $k>0$) in
order to improve the efficiency of the code when studying the propagation of a
GW background along the history of the universe. We take the initial spectrum
for the time derivative of the strains, defined in equation (3.17), to be a
smoothed double broken power law described by

$S_{h^{\prime}}(k,\eta_{*})=S_{h^{\prime}}^{*}\frac{\biggl{[}1+\Bigl{(}\frac{k_{\ast}}{k_{\rm
b}}\Bigr{)}^{a-b}\biggr{]}\Big{(}\frac{k}{k_{*}}\Big{)}^{a}}{\bigg{[}1+\Big{(}\frac{k}{k_{\rm
b}}\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}}},$
(5.1)

where $S_{h^{\prime}}^{*}$ and $k_{*}$ are approximately the initial amplitude
and position of the peak, respectively, $\eta_{*}\approx{\cal H}_{*}^{-1}$
corresponds to the time of GW production during the RD era, and
$\alpha_{1}=\alpha_{2}=2$ are fixed smoothness parameters. We similarly
initialize the spectrum of the strains
$S_{h}(k,\eta_{*})=S_{h^{\prime}}(k,\eta_{*})/k^{2}$. We choose the slope in
the IR range to be $k^{2}$, set by $a=2$ (as is expected for causal sources of
GWs like the ones produced during a phase transition) up to the break wave
number $k_{\rm b}=1$, which corresponds to the horizon scale. At intermediate
wave numbers, $k_{\rm b}\leq k\leq k_{\ast}$, the slope becomes $b=0$ as
found, for example, for MHD turbulence in refs. [79, 84, 85], followed by the
power law $k^{-{11\over 3}}$ in the UV range, set by $c={11\over 3}$. This
corresponds to the spectrum obtained for Kolmogorov-like MHD turbulence [79].
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 $\alpha_{\rm M}^{*}$, the resulting GW spectral shape might be modified by
the inclusion of an additional IR branch $\xi^{\rm WKB}(k)\sim k^{-2}$ as
predicted in section 3.3 by the WKB approximation. Finally, the resulting
spectrum is expected to be enhanced at all wave numbers by a factor
$e^{-2\mathcal{D}}(1+{\textstyle{1\over 2}}\alpha_{\rm T})$ owing to the
presence of non-zero MG parameters $\alpha_{\rm M}$ and $\alpha_{\rm T}$.
Figure 2 shows the resulting GW spectrum estimated using the WKB approximation
for different values of $\alpha_{\rm M}^{*}$. 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.

Figure 2: Left panel: Expected WKB spectra [see equation (3.18)] of the strain derivatives $S_{h^{\prime}}(k)$, compensated by $e^{-2\mathcal{D}}$, 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 equation (5.1). The horizontal line shows the value at the peak $S_{h^{\prime}}^{*}$. Right panel: Expected spectral modifications $\xi^{\rm WKB}(k)\equiv e^{2\mathcal{D}}S_{h^{\prime}}^{\rm WKB}(k)/S_{h^{\prime}}^{\rm GR}(k)$, given in equation (3.20). The black curves indicate the GR solution and the colored curves correspond to different values of $\alpha_{\rm M}^{*}$, with the solid and dashed ones being positive and negative values, respectively. Vertical lines indicate the estimated value of $k_{{\rm lim},\,\alpha_{\rm M}}$, below which the WKB approximation might fail. The inset zooms in on the behaviors around the horizon at $k\in(1,10)$. | run | $\alpha_{\rm M,0}$ | $n$ | $\alpha_{\rm T,0}$ | $\mathcal{E}_{\rm GW}/\mathcal{E}_{\rm GW}^{*}$
---|---|---|---|---|---
Choice 0 | T0A | $0$ | – | $-0.5$ | $\ \,7.50\times 10^{-1}$
($\alpha_{\rm T}$) | T0B | $0$ | – | $-0.2$ | $\ \,9.00\times 10^{-1}$
| T0C | $0$ | – | $0.2$ | $1.10\times 10^{0}$
| T0D | $0$ | – | $0.5$ | $1.25\times 10^{0}$
Choice 0 | M0A | $-0.5$ | – | 0 | $3.44\times 10^{7}$
($\alpha_{\rm M}$) | M0B | $-0.3$ | – | 0 | $3.32\times 10^{4}$
| M0C | $-0.1$ | – | 0 | $3.21\times 10^{1}$
| M0D | $-0.01$ | – | 0 | $1.41\times 10^{0}$
| M0E | $0.1$ | – | 0 | $\ \,3.11\times 10^{-2}$
| M0F | $0.3$ | – | 0 | $\ \,3.03\times 10^{-5}$
Choice I | M1A | $-0.5$ | $2$ | 0 | $1.28\times 10^{0}$
($\alpha_{\rm M}$) | M1B | $-0.3$ | $2$ | 0 | $1.16\times 10^{0}$
| M1C | $-0.1$ | $2$ | 0 | $1.05\times 10^{0}$
| M1D | $0.1$ | $0.4$ | 0 | $\ \,7.79\times 10^{-1}$
| M1E | $0.3$ | $0.4$ | 0 | $\ \,4.73\times 10^{-1}$
Choice II | M2A | $-0.5$ | – | 0 | $1.32\times 10^{0}$
($\alpha_{\rm M}$) | M2B | $-0.3$ | – | 0 | $1.18\times 10^{0}$
| M2C | $-0.1$ | – | 0 | $1.06\times 10^{0}$
| M2D | $0.1$ | – | 0 | $\ \,9.46\times 10^{-1}$
| M2E | $0.3$ | – | 0 | $\ \,8.47\times 10^{-1}$
Choice III | M3A | $-0.5$ | – | 0 | $3.57\times 10^{8}$
($\alpha_{\rm M}$) | M3B | $-0.3$ | – | 0 | $1.35\times 10^{5}$
| M3C | $-0.1$ | – | 0 | $5.12\times 10^{1}$
| M3D | $-0.01$ | – | 0 | $1.48\times 10^{0}$
| M3E | $0.1$ | – | 0 | $\ \,1.95\times 10^{-2}$
| M3F | $0.3$ | – | 0 | $\ \,7.48\times 10^{-6}$
Table 1: 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. $\mathcal{E}_{\rm GW}/\mathcal{E}_{\rm GW}^{*}$
indicates the present-day values of the energy density due to MG with respect
to the corresponding GR values. For each of the runs in series M0 and M3, we
have performed an additional run with $k_{1}=10^{-7}$ to study the spectral
modifications at even smaller $k$.

Table 1 summarizes the parameters of the runs for the numerical studies. The
values of $\alpha_{\rm T,0}$ are chosen to be much larger than the current
constraints in order to later show the relative insignificance of $\alpha_{\rm
T,0}$ even when it takes large values. The choices of $\alpha_{\rm M,0}$ are
made in line with the limits discussed in section 4. For the choice II of
$\alpha_{\rm M}$ parameterization, we take $n=2$ for negative values of
$\alpha_{\rm M,0}$ and $n=0.4$ for positive values to ensure that the
stability conditions of equation (4.3) hold.

For series T0 (T0A–T0D), we evolve the solution entirely with increasing time
steps, such that $\eta_{\rm next}=\eta_{\rm current}(1+\delta n_{\rm incr})$
with $\delta n_{\rm 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 $\Lambda$D such that $\eta_{\rm next}=\eta_{\rm
current}+\delta\eta$ with $\delta\eta\,{\cal H}_{*}=5\times 10^{9}$. We choose
such time evolution to improve the accuracy of late-time results for the cases
of time-dependent $\alpha_{\rm M}$, especially the modifications that they
present in the IR limit. We show in appendix D that decreasing the time step
below $\delta\eta\,{\cal H}_{*}=5\times 10^{9}$ does not affect the IR range
of the spectra, which indicates that the observed modifications are not due to
numerical inaccuracy. 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
appendix D.

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 degrees of freedom are $g_{\rm S}\approx g_{*}\sim 100$, which
yields the values of $a_{*}$ and ${\cal H}_{*}$ in equations (A.8) and (A.9).
With these values, we can then solve equation (2.3) in units of the normalized
time $\eta{\cal H}_{*}$, by mapping the parameterizations in equation (4.2)
$\alpha_{\rm M}(\eta)\rightarrow\alpha_{\rm M}(\eta/\eta_{*})$ and using the
results from Friedmann equations for ${\cal H}$ and $a^{\prime\prime}/a$
(which are already normalized since $a^{\prime}$ is computed as the derivative
with respect to $\eta{\cal H}_{*}$), given in equation (2.9). The effects on
$S_{h^{\prime}}(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^{\prime\prime}/a$ and $\alpha_{\rm M}^{\prime}$ that involve time
derivatives compared to $k$. These terms have been parameterized in equations
(3.11) and (4.4), respectively, and their magnitude has been discussed. We
have shown that the term $a^{\prime\prime}/a$ can only induce modifications to
the solution at scales several orders of magnitude above the horizon scale,
while $\alpha_{\rm M}^{\prime}$ is only of the order of $\alpha_{\rm M}$ for
choice III (see section 4) at very late times within the RD era. For other
choices of $\alpha_{\rm M}(\eta)$, either $\alpha_{\rm M}^{\prime}=0$ (choice
0) or $\alpha_{\rm M}$ itself is orders of magnitude below its present-time
value $\alpha_{\rm M,0}$ (choices I and II). We find that the WKB
approximation is expected to be valid around $k\gg k_{{\rm lim},\,\alpha_{\rm
M}}\sim\sqrt{|\alpha_{\rm M}^{*}|/2}$ and, on the other hand, it predicts an
enhancement $k^{-2}$ in the IR regime $k\lesssim k_{\rm crit}\leq k_{{\rm
lim},\,\alpha_{\rm M}}$, which is not on the range of validity of the WKB
estimate. Using the results of the numerical simulations, we find in section
5.3 an IR enhancement $\xi(k)\propto k^{-\beta_{0}}$ with $\beta_{0}\sim 2$
and $\beta_{0}\gtrsim 0$ for negative and positive values of $\alpha_{\rm
M,0}$, respectively, 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 $\alpha_{\rm M}^{*}$, 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 equation (3.16).

The amplitude $S_{h^{\prime}}^{*}$ and the position of the peak $k_{*}$ can be
chosen to represent a specific model. For example, the values
$S_{h^{\prime}}^{*}=5\times 10^{-10}$ and $k_{*}=10$ are used in our runs (see
table 1) and they would produce a normalized initial total energy density777
The turbulent and GW energy densities are normalized to the radiation energy
density, such that $\mathcal{E}_{\rm GW}$ during RD can be computed as
$\mathcal{E}_{\rm
GW}(\eta)=\frac{1}{6}\int_{0}^{\infty}S_{h^{\prime}}(k,\eta){\rm\,d}k.$ For
the spectral shape $S_{h^{\prime}}(k)$ given in equation (5.1) with
$k_{*}=10$, one finds $\mathcal{E}_{\rm GW}\simeq 2S_{h^{\prime}}^{*}$ (see
the associated online material [74]).  $\mathcal{E}_{\rm GW}^{*}\simeq
10^{-9}$, which corresponds to a vortically turbulent source with an energy
density of roughly $q\mathcal{E}_{\rm turb}^{*}\simeq 1.74\times 10^{-4}$.
This is related via

$\mathcal{E}_{\rm GW}^{*}=(q\mathcal{E}_{\rm turb}^{*}/k_{\rm s})^{2},$ (5.2)

where $q$ is an empirically determined coefficient for a specific type of
turbulence source $\mathcal{E}_{\rm 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 [79, 81, 82], and $k_{\rm s}$ is the characteristic
scale of the turbulence sourcing, which is related to $k_{*}$ in equation
(5.1) as $k_{*}\simeq 1.143\times 1.6\,k_{\rm s}\simeq 1.83\,k_{\rm s}$, where
the factor 1.143 relates the position of the maximum of $\Omega_{\rm
GW}(k)\propto kS_{h^{\prime}}(k)$ and $k_{*}$ (see the associated online
material [74]) and $1.6$ the relation between $k_{\rm s}$ and the maximum
[84].

###  5.2 Time evolution

For constant $\alpha_{\rm T}\neq 0$ and $\alpha_{\rm M}=0$, the GW energy
density stays constant in time and its magnitude is modified by the specific
value of $\alpha_{\rm T}$, as it can be predicted from the analytical solution
to equation (2.3) in the absence of sources during RD, given in equation
(3.20) under the WKB approximation. We find excellent agreement between the
WKB estimate and the numerical solution of the GW energy density. Note that
this corresponds to a boost of energy for $\alpha_{\rm T}>0$ and a depletion
for $\alpha_{\rm T}<0$, given by $1+{\textstyle{1\over 2}}\alpha_{\rm T}$.
Hence, the relative changes upon the GR solutions are of order $\alpha_{\rm
T}\sim\mathcal{O}(10^{-1})$ (see table 1), which would become much smaller if
we restrict $\alpha_{\rm T}$ to the present-time constraint at the LIGO–Virgo
frequency band $\alpha_{\rm T,0}\lesssim\mathcal{O}(10^{-15})$.

Figure 3 presents the time evolution of the total GW energy density
$\mathcal{E}_{\rm GW}(\eta)$ of the runs in series M0 (upper left), M1 (upper
right), M2 (lower left), and M3 (lower right). In these runs, $\alpha_{\rm M}$
follows each of the parameterizations given in equation (4.2) and $\alpha_{\rm
T}=0$. Regardless of the specific parameterization, the numerical solutions
and the WKB approximations agree on the total GW energy density time
dependence, enhanced or depleted by a factor $e^{-2\mathcal{D}}$ for negative
or positive values of $\alpha_{\rm M,0}$, respectively.

Figure 3: Time evolution of the normalized total energy density
$\mathcal{E}_{\rm GW}/\mathcal{E}_{\rm GW}^{*}$ for choices 0 to III and
$\alpha_{\rm M,0}\in\\{-0.5,-0.3,-0.1,-0.01,0.1,0.3\\}$. In all panels, the
WKB estimates are shown as solid lines and numerical solutions as dashed
lines, being almost indistinguishable. For choices I and II, the inset
corresponds to the $\Lambda$D era, when there are relevant modifications to
the GR solution $\mathcal{E}_{\rm GW}^{*}$. All runs are initialized at the
EWPT such that the present time is $\eta_{0}/\eta_{*}\simeq 2.4\times
10^{13}$.

Figure 4: Relative error on the time series in figure 3 of the WKB
approximation compared to the numerical simulations, $\varepsilon^{\rm
WKB}({\cal E}_{\rm GW})$.

Choices 0 and III yield similar results, with a time evolution
$\mathcal{E}_{\rm GW}(\eta)\sim\eta^{-\alpha_{\rm M,0}}$ and $\mathcal{E}_{\rm
GW}(\eta)\sim\eta^{-1.5\,\alpha_{\rm M,0}}$ during RD era, and
$\mathcal{E}_{\rm GW}\sim\eta^{-2\alpha_{\rm M,0}}$ and $\mathcal{E}_{\rm
GW}(\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-\Omega_{{\rm mat},0})\simeq 1.46$, leading to the different scalings. In
other words, the evolution of $\mathcal{E}_{\rm GW}$, proportional to
$a^{-\alpha_{\rm M}}$ when $\alpha_{\rm M}$ 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 density that heavily depends on the values of $\alpha_{\rm
M,0}$, i.e., $\mathcal{E}_{\rm GW}/\mathcal{E}_{\rm GW}^{*}$ can range from
$\mathcal{O}(10^{-5})$ (choice 0) and $\mathcal{O}(10^{-6})$ (choice III) for
$\alpha_{\rm M,0}=0.3$ to $\mathcal{O}(10^{7})$ (choice 0) and
$\mathcal{O}(10^{8})$ (choice III) for $\alpha_{\rm M,0}=-0.5$. In general,
choice III leads to a larger modification due to larger values of $\alpha_{\rm
M}$ during RD than choice 0. The potential implications of such a large GW
energy density enhancement are discussed in section 6.

For both choices I and II (see the upper right and lower left panels of figure
3), the modified GW solutions remain close to their GR counterparts for most
of the time, and rapidly depart from GR as $\eta$ enters $\Lambda$D and
approaches the present day. This is expected since $\alpha_{\rm M}$ is
proportional to the scale factor and the dark energy density for choices I and
II, respectively, and hence, the values of $\alpha_{\rm M}$ are negligibly
small for most of the cosmic history until $\Lambda$D era. For this reason,
for the same values of $\alpha_{\rm M,0}$, the final values of
$\mathcal{E}_{\rm GW}$ in both of these cases are significantly lower than
those in choices 0 and III, where the modifications are accumulated from RD
onward.

In figure 3, the differences between the WKB and the numerical solutions are
indistinguishable. We quantify in figure 4 the relative error between the two,
defined to be $\varepsilon^{\rm WKB}(\mathcal{E}_{\rm
GW})\equiv[\mathcal{E}_{\rm GW}^{\rm WKB}(\eta)-\mathcal{E}_{\rm
GW}(\eta)]/\mathcal{E}_{\rm GW}(\eta)$. In all panels, there exists a brief
but relatively large error region, amounting to $\lesssim 15\%$, around the
initial time $\eta_{*}$. This is due to the sinusoidal oscillations of each
$k$ mode that are present in GR and MG alike [see equation (3.8)]. It is a
consequence of setting an already saturated spectra $S_{h^{\prime}}(k)$ as the
initial condition at time $\eta_{*}$ and hence, it is not due to the WKB
estimate. After the initial oscillations settle down, $\varepsilon^{\rm
WKB}(\mathcal{E}_{\rm GW})$ for choices 0 and III decrease over the majority
of time in RD era. Although they both increase somewhat later on, the maximum
error at the final time is only $\varepsilon^{\rm WKB}(\mathcal{E}_{\rm
GW})\sim 3\%$. For choices I and II, due to the negligible impact of
$\alpha_{\rm M}$ during RD, the GW solutions settle down to the same magnitude
after the initial oscillations. Therefore, the relative errors during RD
remain close to zero. During $\Lambda$D, the errors grow as the effects of
$\alpha_{\rm M}$ become more significant. But this is also limited, since at
most $\varepsilon^{\rm WKB}(\mathcal{E}_{\rm GW})\sim 1\%$ is found at the
present day, lower than in choices 0 and III.

###  5.3 Energy spectrum

In figure 5, we show the final energy spectrum averaged over oscillations and
compensated by the factor $e^{-2\mathcal{D}}$ to study the changes on the
spectral shape. The saturated spectra are shown for the choices 0 and III of
$\alpha_{\rm M}$ parameterization, since the rest of runs ($\alpha_{\rm T}$
choice 0, $\alpha_{\rm M}$ choices I and II) exhibit the same spectral shape
as in GR. To study the spectral changes, we again define
$\xi(k)=e^{2\mathcal{D}}S_{h^{\prime}}(k)/S_{h^{\prime}}^{\rm GR}(k)$, in
analogy to equation (3.18). As a reminder of the potential limitations of the
WKB approximation, we mark with thin vertical lines the values of $k_{{\rm
lim},\,a^{\prime\prime}}^{\rm EW}$ and $k_{{\rm lim},\,\alpha_{\rm M}}$, and
note that the latter occurs at wave numbers larger than the IR regime
characterized by $k_{\rm crit}={\textstyle{1\over 2}}|\alpha_{\rm M}^{*}|$.

Figure 5: Saturated final energy spectra compensated by the change on the
total GW energy density $e^{2\mathcal{D}}S_{h^{\prime}}^{\rm
num}(k)/S_{h^{\prime}}^{*}$ for the numerical runs with the choices 0 and III
of $\alpha_{\rm M}$ parameterization. Left and right panels represent negative
and positive values of $\alpha_{\rm M,0}$, respectively. The GR solution and
WKB approximation are shown in black solid and color dotted lines,
respectively. Numerical solutions are shown in solid (choice 0) and dash-
dotted (choice III) lines. Thin vertical lines indicate the corresponding
limits to the WKB approximation $k_{{\rm lim},\,a^{\prime\prime}}^{\rm EW}$
and $k_{{\rm lim},\,\alpha_{\rm M}}$. The horizontal lines show the values at
the peak $S_{h^{\prime}}^{*}$.

Comparing the numerical results to the expected spectra obtained using the WKB
approximation (see figure 2), we note two main differences. In the first
place, we find an IR enhancement $\xi(k)\sim\\!k^{-\beta_{0}}$ with
$\beta_{0}\in(1,2)$, where $\beta_{0}=2$ corresponds to the WKB approximation,
for negative values of $\alpha_{\rm M,0}$, whereas for positive values there
are more moderate modifications with $\beta_{0}\lesssim 0.5$. Moreover, even
for negative $\alpha_{\rm M,0}$, the modifications to the spectral shape are
different than those predicted by the WKB estimate. These modifications can be
seen more clearly in figure 6, where the left panel shows the numerically
obtained $\xi(k)$ and the right panel directly presents the spectral slopes of
the compensated spectra $\beta\equiv-\partial\log\xi(k)/\partial\log k$,
measured at each wave number.

Figure 6: Left and right panels respectively show $\xi^{\rm num}(k)$, the
saturated spectra normalized with the GR solution, and $\beta$, its spectral
slope as a function of wave number. We show, for comparison, the slopes
$\beta$ predicted by the WKB estimate for $\alpha_{\rm M,0}=-0.5$ and $0.3$.
The horizontal thin black lines indicate the GR solution in the left panel and
the slopes $0$ and $2$ in the right panel, for reference. Different values of
$\alpha_{\rm M,0}$ and line colors are consistent with the ones used in
figures 3–5.

The slope $\beta$ of the WKB estimate can be computed from equation (3.20) as

$\beta=\frac{2\alpha_{\rm M}^{*}\bigl{(}\alpha_{\rm
M}^{*}+2k\bigr{)}}{8k^{2}+\alpha^{*2}_{\rm M}+4k\alpha_{\rm M}^{*}},$ (5.3)

which has an initial departure $\beta<0$ for $\alpha_{\rm M}^{*}<0$ as $k$
decreases until it reaches a peak $\beta_{\rm min}=1-\sqrt{2}\simeq-0.41$ at
$k_{\beta_{\rm min}}={\textstyle{1\over 4}}\alpha_{\rm
M}^{*}\sqrt{2}/(1-\sqrt{2})\simeq 0.85\,|\alpha_{\rm M}^{*}|$. It then grows
again to zero at $k_{\rm crit}$ and continues to increase up to $\beta_{\rm
max}=1+\sqrt{2}\simeq 2.41$ at $k_{\beta_{\rm max}}=-{\textstyle{1\over
4}}\alpha_{\rm M}^{*}\sqrt{2}/(1+\sqrt{2})\simeq 0.15\,|\alpha_{\rm M}^{*}|$,
crossing the aforementioned IR enhanced slope of $-2$ at the wave number
$k_{\beta=2}={\textstyle{1\over 4}}|\alpha_{\rm M}^{*}|$ (for details, see the
associated online material [74]). Finally, the slope becomes slightly
shallower and asymptotically tends to $-2$ as $k\rightarrow 0$, as can be seen
by taking the limit in equation (5.3). For positive values of $\alpha_{\rm
M}^{*}$, the slope $\beta$ monotonically increases from $0$ to $2$ with
decreasing $k$ so it does not present any local minima or maxima. The slope
$\beta$ for $\alpha_{\rm M}^{*}=-0.5$ and $0.3$ under the WKB estimate is
shown in the right panel of figure 6 for comparison with the numerical
results.

At wave numbers close to where the IR regime starts, $k\lesssim k_{\rm crit}$,
the $k^{-2}$ enhancement is reproduced by the numerical results for negative
$\alpha_{\rm M,0}$, flattening that part of the original GR spectrum to
$k^{0}$ (see figure 5). But as $k$ decreases, the numerical spectra exhibit
slopes steeper than $k^{0}$, closer to the original GR $k^{2}$ spectrum, that
are not predicted by the WKB approximation, which maintains the $k^{0}$ slope
throughout all scales $k\lesssim k_{\rm crit}$. The departure in the IR regime
occurs at different scales depending on the values of $\alpha_{\rm M,0}$. For
$\alpha_{\rm M,0}=-0.01$, this is the least obvious, where the enhanced
spectrum still keeps a roughly $k^{0}$ shape even at $k\ll k_{\rm crit}$. For
$\alpha_{\rm M,0}=-0.5$, however, the difference becomes visible when $k$ is
still relatively close to the start of the IR regime at $k\sim k_{\rm crit}$.
We characterize the differences of the numerical results with those from the
WKB estimate in table 2 by computing the critical $k$ at which $\beta=2$,
$k_{\beta=2}^{\rm crit}$. Under WKB this only occurs at ${\textstyle{1\over
4}}|\alpha_{\rm M}^{*}|$ and in the limit $k\rightarrow 0$ while in the
simulations it occurs at a similar $k$ and then again at a larger but finite
scale $k_{\beta=2}^{{\rm crit},2}$, when the spectra becomes steeper (with
respect to the flat spectrum induced by $\xi^{\rm WKB}(k)\sim k^{-2}$). We
also give the final slope $\beta_{0}$, such that $\xi(k)\sim k^{-\beta_{0}}$
as $k\rightarrow 0$.

On the other hand, for positive $\alpha_{\rm M,0}$, the IR enhancement
predicted by the WKB estimate is not observed in the numerical results. In
this case, the modified spectra present small increases around $k\gtrsim 1$
and slightly shallower slopes than their GR counterparts at all superhorizon
scales, characterized by $\beta_{0}\gtrsim 0$.

Finally, to quantify the departures of the WKB approximation from the
numerical solutions, especially at $k\leq k_{{\rm lim},\,\alpha_{\rm M}}$, we
show the relative errors in the saturated spectra as $\varepsilon^{\rm
WKB}(S_{h^{\prime}})\equiv[S_{h^{\prime}}(k)^{\rm
WKB}-S_{h^{\prime}}(k)]/S_{h^{\prime}}(k)$ in figure 7. We observe that for
positive values of $\alpha_{\rm M,0}$, the errors grow monotonically as $k$
becomes smaller. In fact, $\varepsilon^{\rm WKB}(S_{h^{\prime}})\sim k^{-2}$
is found due to a consistent disagreement between the flat spectra predicted
by the WKB formalism and $k^{-\beta_{0}}$ with $\beta_{0}\gtrsim 0$ obtained
numerically (see table 2). On the other hand, for negative values of
$\alpha_{\rm M,0}$, $\varepsilon^{\rm WKB}(S_{h^{\prime}})$ also increase as
$k$ becomes smaller. But since $\alpha_{\rm M,0}<0$ still produces a similar
IR enhancement than that predicted by the WKB, the slope in the error spectrum
does not become as steep as $k^{-2}$ and depends on the corresponding values
of $\beta_{0}$ in table 2. In all of the runs, the errors are bounded by
$\lesssim 15\%$ at $k\gtrsim k_{{\rm lim},\,\alpha_{\rm M}}$, where the WKB
estimate is expected to be valid.

| $\alpha_{\rm M,0}$ | WKB | Choice 0 ($\alpha_{\rm M}$) | Choice III ($\alpha_{\rm M}$)
---|---|---|---|---
$k^{\rm crit}_{\beta=2}$ | $-0.5$ | $1.25\times 10^{-1}$ | $1.32\times 10^{-1}$ | $1.79\times 10^{-1}$
| $-0.3$ | $7.5\times 10^{-2}$ | $9.04\times 10^{-2}$ | $1.20\times 10^{-1}$
| $-0.1$ | $2.5\times 10^{-2}$ | $3.83\times 10^{-2}$ | $5.16\times 10^{-2}$
| $-0.01$ | $2.5\times 10^{-3}$ | $4.58\times 10^{-3}$ | $6.63\times 10^{-3}$
$k^{{\rm crit},2}_{\beta=2}$ | $-0.5$ | 0 | $6.73\times 10^{-2}$ | $9.93\times 10^{-2}$
| $-0.3$ | 0 | $3.09\times 10^{-2}$ | $5.54\times 10^{-2}$
| $-0.1$ | 0 | $4.42\times 10^{-3}$ | $9.62\times 10^{-3}$
| $-0.01$ | 0 | $5.04\times 10^{-5}$ | $7.44\times 10^{-5}$
$\beta_{0}$ | $-0.5$ | $2$ | $1.38$ | $0.95$
| $-0.3$ | $2$ | $1.66$ | $1.47$
| $-0.1$ | $2$ | $1.89$ | $1.84$
| $-0.01$ | $2$ | $2.00$ | $1.99$
| $0.1$ | $2$ | $0.10$ | $0.15$
| $0.3$ | $2$ | $0.30$ | $0.44$
Table 2: Spectral features in the IR regime for the values of $\alpha_{\rm
M,0}$ used in the numerical runs of table 1 and for choices 0 and III. $k^{\rm
crit}_{\beta=2}$ and $k^{{\rm crit},2}_{\beta=2}$ indicate the wave numbers at
which the spectral slope becomes $-2$, which occurs at
$k_{\beta=2}=\\{0,{\textstyle{1\over 4}}|\alpha_{\rm M}^{*}|\\}$ for the WKB
estimate. $\beta_{0}$ is the asymptotic slope in the numerical runs, such that
$\xi(k)\sim k^{-\beta_{0}}$ at large scales, which corresponds to
$\beta_{0}=2$ for the WKB estimate.

Figure 7:  Relative error $\varepsilon^{\rm WKB}(S_{h^{\prime}})$ in the
saturated spectra at late times of the WKB approximation compared to the
numerical simulations. The vertical lines correspond to the wave number
$k_{{\rm lim},\,\alpha_{\rm M}}$, below which the WKB estimate is expected to
break down. The gray shaded regions in both panels represent relative errors
of $\leq 15\%$. Different values of $\alpha_{\rm M,0}$ and line colors are
consistent with the ones used in figures 3–5.

##  6 Observational implications

To infer the observational prospects of detecting a modified GW background, we
convert the linear energy spectrum $S_{h^{\prime}}(k)$ directly obtained from
the code to the commonly used logarithmic energy spectrum $\Omega_{\rm GW}(f)$
via equation (3.16). As a result, the GR solution in the form of the double
broken power law used in the numerical simulations, i.e.,
$S_{h^{\prime}}(k)\propto k^{2}$, $k^{0}$, and $k^{-{11\over 3}}$ in the low,
intermediate, and high wave number regimes, given by the smoothed double
broken power law in equation (5.1), becomes $\Omega_{\rm GW}(f)\propto f^{3}$,
$f$, and $f^{-{8\over 3}}$, respectively. The normalized wave numbers used in
the code are converted to the present-day physical frequencies via

$f=\frac{c_{\rm T}k\mathcal{H}_{*}}{2\pi}\,\frac{a_{*}}{a_{0}},$ (6.1)

where ${\cal H}_{*}$ and $a_{*}/a_{0}$ are given in equations (A.8) and (A.9)
and depend on the time at which the GWs are generated. The term $c_{\rm T}$
appears in equation (6.1) due to the modified dispersion relation at present
time. We have found that only $\alpha_{\rm M}$ introduces spectral changes, so
we focus on the cases with $c_{\rm T}=1$. However, modified theories of
gravity with $c_{\rm T}\neq 1$ would allow the resulting GW spectra
$\Omega_{\rm GW}(f)$ to be shifted in $f$ by $c_{\rm T}$.

Assuming that the origin of the GWs is the EWPT and that the characteristic
scale of the source is ${\cal H}_{*}\lambda_{\rm s}=10^{-2}$ (given, for
example, by the mean separation of the broken-phase bubbles in a first-order
EWPT [102]), then $\Omega_{\rm GW}(f)$ today peaks around the mHz band, which
corresponds to the peak sensitivity of the Laser Interferometer Space Antenna
(LISA). On the other hand, if the GW signal is produced at the QCDPT with a
characteristic length scale of the order of the Hubble scale ${\cal
H}_{*}\lambda_{\rm s}\approx 1$, then it would be compatible [84] with the
recent observations by the different pulsar timing array (PTA) collaborations
[103, 104, 105, 106] in the nHz band of frequencies, which are, however, not
yet confirmed to correspond to a GW background. In particular, for a signal
produced from MHD turbulence, the peak is estimated to occur at $k_{\rm
GW}\simeq 1.6\times 2\pi/({\cal H}_{*}\lambda_{\rm s})$ [84], which is
obtained by taking $k_{*}\simeq 0.875\,k_{\rm GW}$ in equation (5.1).

Choosing the cases with the most pronounced modifications, we show in figure 8
the resulting $\Omega_{\rm GW}(f)$ of the runs using the choice III of
$\alpha_{\rm M}$ parameterizations, i.e., the runs of series M3; see table 1,
for the cases in which the GW signal is produced at the EWPT and at the QCDPT.
For observational comparisons, we show the power law integrated sensitivity
(PLIS) curves [107] of proposed future detectors such as Square Kilometer
Array (SKA) [108], DECi-hertz Interferometer Gravitational wave Observatory
(DECIGO) [109], Big Bang Observer [110], LISA [111], and Einstein Telescope
(ET) [112].

Figure 8: Present-day modified GW energy spectrum for the numerical runs of
series M3 (see table 1) for GW signals produced at the EWPT (upper panel) and
at the QCDPT (lower panel). GR solution is shown in black. Different values of
$\alpha_{\rm M,0}$ and line colors are consistent with the ones used in
figures 3–5. Several detectors’ PLIS curves for a signal-to-noise ratio of 10
are shown in gray for comparison. The recent PTA reported observations are
shown as shaded regions in blue (NANOGrav), red (PPTA), purple (EPTA), and
black (IPTA) for their $2\sigma$ confidence amplitudes in the range of slopes
$\gamma\geq 0$ for $\Omega(f)\sim f^{\gamma}$. Figure 9: Total amplification
or depletion factor as a function of $\alpha_{\rm M,0}$ for a GW signal
produced at the EWPT (blue) and at the QCDPT (green) for the different choices
of $\alpha_{\rm M}$ parameterization in equation (4.2).

Due to the factor $e^{-2\mathcal{D}}$, negative values of $\alpha_{\rm M,0}$
boost the GW spectral amplitude, which increases the chance of GW detection.
Hence, the (non-)observation of a GW signal could serve as a constraint on the
value of $\alpha_{\rm M,0}$. However, the present-day amplitude of the GW
spectrum is degenerate with respect to the value of $e^{-2\mathcal{D}}$ at
present time (shown in figure 9, it depends on $\alpha_{\rm M,0}$, the time-
dependent parameterization of $\alpha_{\rm M}$, and the time of generation of
the GWs) as well as the GW energy density produced at the time of generation
$\mathcal{E}_{\rm GW}^{*}$. We show in figure 8, for different values of
$\alpha_{\rm M,0}$, the expected spectra at present time obtained for fixed
values $e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}\simeq 2.5\times
10^{-7}$ and $4.44\times 10^{-1}$ at the EWPT with $k_{\rm s}=2\pi/({\cal
H}_{*}\lambda_{\rm s})\simeq 600$, which correspond to the spectral peaks at
present time $h^{2}\,\Omega_{\rm GW}\simeq 2.34\times 10^{-12}$ and
$4.15\times 10^{-6}$ at $2.5\,$mHz. Note that the former amplitude would
correspond to $q\mathcal{E}_{\rm turb}^{*}\simeq 0.3$ in GR while the latter
would be $q\mathcal{E}_{\rm turb}^{*}\simeq 400$. This is completely
unrealistic in GR while in MG, it would again be compatible with
$q\mathcal{E}_{\rm turb}^{*}\simeq 0.3$ if $e^{-2\mathcal{D}}(\eta_{0})\simeq
2\times 10^{6}$ that can be obtained if $\alpha_{\rm M,0}\approx-0.3$ (see
choice III in figure 9). In this case, the GW signal could be amplified by MG
such that the low-frequency tail could become detectable by SKA and the
different PTA collaborations in the future, in a way that would break the
degeneracy by a potential observation over different frequency bands. At the
nHz frequencies, the $f^{3}$ branch is modified to a different slope that
depends on the value of $\alpha_{\rm M}^{*}$ at the time of generation,
yielding a potential indirect observation of the MG parameter during RD era.
The non-observation of such a signal at nHz frequencies would similarly
constrain the value of $\alpha_{\rm M}^{*}$ if a compatible signal is observed
by LISA. Otherwise, only the product
$e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}$ could be constrained
by not detecting such a signal in the LISA band. Similarly, an even larger
amplification $e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}\simeq 50$
could in principle allow the large frequency branch $f^{-{8/3}}$ to fall into
ET sensitivity. Hence, ET could also be used to put upper bounds on the
amplifications, although less stringent than those combining LISA and SKA.

For the QCDPT, we take the fixed value
$e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}\simeq 5\times 10^{-4}$
with $k_{\rm s}\simeq 10$, which corresponds to a spectral peak
$h^{2}\,\Omega_{\rm GW}\simeq 7.4\times 10^{-9}$ at $50.6\,$nHz. In GR, this
corresponds to $q\mathcal{E}_{\rm turb}^{*}\simeq 0.22$ so it is compatible
with, for example, the primordial magnetic fields proposed in refs. [113, 83,
84] in connection to the reported PTA observations. In this case, the horizon
at the QCDPT scale is within the range of frequencies where PTAs are sensitive
and hence, the amplitude can be slightly modified by the value of $\alpha_{\rm
M}^{*}$ at the time of generation. Such modifications are small and at the
moment, the data is not accurate enough for such a precision measurement.
However, this is an interesting feature that occurs around the horizon and the
potential observation of a deviation with respect to the GR spectrum around
the horizon scale could be an indication of non-zero $\alpha_{\rm M}^{*}$
around the QCD scale. On the other hand, the current PTA observations allow us
to constrain the product $e^{-2\mathcal{D}}\mathcal{E}_{\rm GW}^{*}$ for
different values of $k_{\rm s}$ since we can find their values at which the GW
signal would be larger than the reported PTA common-noise observations:
$e^{-2\mathcal{D}}\mathcal{E}_{\rm GW}^{*}\leq 2.5\times 10^{-3}$ for $k_{\rm
s}\approx 10$, $\leq 2.5\times 10^{-2}$ for $k_{\rm s}\approx 60$, $\leq
4\times 10^{-2}$ for $k_{\rm s}\approx 100$, $\leq 2.5\times 10^{-1}$ for
$k_{\rm s}\approx 600$.

In addition, the expected value of $\mathcal{E}_{\rm GW}^{*}$ depends on the
mechanisms leading to GW production, in particular, their efficiency $q$,
which presents a fair amount of uncertainty. Recall equation (5.2), which
states that the initial GW energy is determined by the initial turbulent
source energy $\mathcal{E}_{\rm turb}^{*}$, the production efficiency $q$, and
the characteristic size of turbulent eddies $k_{\rm s}$. Therefore, in
practice, constraining $\alpha_{\rm M,0}$ via its amplification effects can be
rather challenging.

Finally, we note that large amplifications due to MG could potentially lead to
GW signals produced at the EWPT that would dominate the sensitivity budget of
various detectors, mostly LISA, but also DECIGO and BBO, for example if the
product $e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}\sim 10^{-1}$
(see figure 8). Hence, at the moment that GW signals of astrophysical origin
that propagate over a shorter period of cosmic history are detected by LISA,
DECIGO and/or BBO, one can put upper bounds to
$e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}$.

##  7 Conclusions

GWs in modified gravity exhibit features in their energy spectra different
from what can be expected from GR. We have explored the significance of such
features in terms of the spectral slopes and amplitudes in different frequency
ranges, under different time-dependent functional forms and values of two
modification parameters — the GW friction $\alpha_{\rm M}$ and the tensor
speed excess $\alpha_{\rm T}$.

We have computed the expected GW energy density spectra using the WKB
approximation, and have found that the GW energy density depends on the time
evolution of $\alpha_{\rm T}$ and $\alpha_{\rm M}$, while its spectral shape
only depends on $\alpha_{\rm M}^{*}$ at the time of GW generation,and the time
dependence of $\alpha_{\rm T}$. However, a number of limiting scales below
which the WKB approximation might not hold have been identified in equations
(3.12)–(3.15), motivating numerical studies to go beyond the WKB estimate,
especially at superhorizon scales. For this purpose, we have performed
numerical simulations using the Pencil Code to study the propagation of GWs
along the history of the universe from their time of generation, assumed to be
during the RD era, e.g., at the EWPT or at the QCDPT, until the present time.

The current-day value of the tensor speed excess is already tightly
constrained by the precise measurement of the speed of GWs to $\alpha_{\rm
T,0}\lesssim\mathcal{O}(10^{-15})$. In this paper, we have found that the
modifications on the GW energy density spectra introduced by a constant
$\alpha_{\rm T}$ are negligible, even when $\alpha_{\rm T}$ takes larger
values (as it might be allowed where it is not constrained, e.g., in the past
or at frequencies out of the LIGO–Virgo band). In addition, $\alpha_{\rm T}$
can modify the dispersion relation, leading to a shift in $\Omega_{\rm GW}(f)$
by $c_{\rm T}$. On the other hand, $\alpha_{\rm M}$ effectively acts as an
(anti-)damping term in the GW equation and can impart changes to the energy
spectrum of GWs with potential implications to future observations.
Considering sources of GWs at the EWPT and QCDPT, the present-day GW spectrum
is expected to peak in the mHz and nHz frequency bands, respectively.

Compared to the standard GW spectra obtained in GR, the changes due to
$\alpha_{\rm M}$ are characterized by an overall boost or depletion factor
$e^{-2\mathcal{D}}$, where $\mathcal{D}$ is given in equation (3.4). However,
the GW energy density at present time due to modified gravity is degenerate
with the produced density at the time of generation $\mathcal{E}_{\rm
GW}^{*}$. Therefore, potential observations or non-observations of the GW
spectra by PTA, SKA, LISA, BBO, DECIGO, ET, or others, could only place bounds
on the quantity $e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}$ and
not on $\alpha_{\rm M}^{*}$ or $\mathcal{E}_{\rm GW}^{*}$ separately if we
only consider the amplitude at the peak of the GW spectrum. For example, we
have provided some upper bounds on
$e^{-2\mathcal{D}}(\eta_{0})\,\mathcal{E}_{\rm GW}^{*}$ using the reported
observations of a common-process spectrum over several pulsars by the
different PTA collaborations; see section 6.

In addition, the GW spectrum also presents spectral modifications $\xi(k)$
around the horizon and at superhorizon scales. If $\alpha_{\rm M}^{*}<0$, the
modification leads to an enhancement of the GW spectrum $\xi(k)\sim
k^{-\beta_{0}}$ with $\beta_{0}\in(0,2)$, which flattens the low-frequency
part of the spectrum; see table 2. The enhancement found from the numerical
simulations is different than the one predicted by the WKB estimate
($\beta_{0}=2$) but it still increases the detection prospects of the low-
frequency tail of the GW spectra. We also find more moderate modifications to
the spectral shape around the horizon scale for all $\alpha_{\rm M}^{*}$ and
at all scales when $\alpha_{\rm M}^{*}>0$. These modifications allow us to
characterize the resulting spectra due to MG.

The degeneracy between $\alpha_{\rm M}$ and $\mathcal{E}_{\rm GW}^{*}$ can be
broken by observing the signal at multiple frequency bands, as $\alpha_{\rm
M}^{*}$ leads to a low-frequency spectral enhancement, while under GR such
enhancement does not appear. This implies that, for example, a GW spectrum
peaking in the mHz band can potentially have its tail in the nHz band, above
PTA or SKA sensitivities. The observation of such a low-frequency tail would
then be a direct indication of a non-zero $\alpha_{\rm M}$ at the time at
which the GW signal has been produced since the slope and amplitude would only
depend on $\alpha_{\rm M}^{*}$.

Overall, a cosmological GW spectrum spanning many orders of magnitudes in
frequency has the potential to constrain the parameter space of modified
theories of gravity, especially if the modifications occur in the early
universe such as the RD era.

#### Data availability.

The source code used for the numerical solutions of this study, the Pencil
Code, along with the additions included for the present study, is freely
available [72]; see also ref. [114] for the numerical data. The calculations
and the routines generating the plots are publicly available on
GitHub888https://github.com/AlbertoRoper/GW_turbulence/tree/master/horndeski
[74].

## Acknowledgments

Support through the grant 2019-04234 from the Swedish Research Council
(Vetenskapsrådet) is gratefully acknowledged. Nordita is sponsored by
Nordforsk. A.R.P. acknowledges support by the Swiss National Science
Foundation (SNSF Ambizione grant 182044), the French National Research Agency
(ANR) project MMUniverse (ANR-19-CE31-0020), and the Shota Rustaveli National
Science Foundation (SRNSF) of Georgia (grant FR/18-1462). 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, and by the Grand Equipement National de Calcul
Intensif (GENCI) to the project “Opening new windows on Early Universe with
multi-messenger astronomy” (A0110412058).

##  Appendix A Friedmann equations

Friedmann equations can be expressed in terms of the energy density
$\Omega(a)$ and a smooth equation of state (EOS) $w(a)$ as

$\frac{\ddot{a}}{a}=-\frac{1}{2}H_{0}^{2}\,\Omega(a)\bigl{[}1+3w(a)\bigr{]},\quad\biggl{(}\frac{\dot{a}}{a}\biggr{)}^{2}=H_{0}^{2}\,\Omega(a),$
(A.1)

where $\Omega(a)$ is defined to be the ratio of total energy density to the
present-day critical energy density $\rho_{{\rm crit},0}\equiv
3H_{0}^{2}/\kappa$, i.e.,

$\displaystyle\Omega(a)=\frac{\rho(a)}{\rho_{{\rm crit},0}}=$
$\displaystyle\,\Omega_{\rm rad}(a)+\Omega_{\rm mat}(a)+\Omega_{\Lambda,0}$
$\displaystyle=$ $\displaystyle\,\biggl{(}{a\over
a_{0}}\biggr{)}^{-4}\,\frac{g_{*}}{g_{*}^{0}}\biggl{(}\frac{g_{\rm S}}{g_{\rm
S}^{0}}\biggr{)}^{-{4\over 3}}\,\Omega_{{\rm rad},0}+\biggl{(}{a\over
a_{0}}\biggr{)}^{-3}\,\Omega_{{\rm mat},0}+\Omega_{\Lambda,0}.$ (A.2)

Numerically, we take the present-time values of $\Omega_{\Lambda,0}\simeq
0.684$, $\Omega_{{\rm rad},0}\simeq 4.16\times 10^{-5}\,h^{-2}$, and
$\Omega_{{\rm mat},0}=1-\Omega_{{\rm rad},0}-\Omega_{\Lambda,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 by CMB experiments [6], and $g_{*}^{0}\simeq 3.36$
and $g_{\rm S}^{0}\simeq 3.91$ are the reference relativistic and adiabatic
degrees of freedom at the present time.999Note that neutrinos’ contribution to
the radiation energy density is accounted for by taking
$g_{*}^{0}=2\Biggl{[}1+N_{\rm
eff}\frac{7}{8}\biggl{(}\frac{4}{11}\biggr{)}^{4\over 3}\Biggr{]}\simeq
3.363,$ where $N_{\rm eff}\approx 3$, 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 degrees of freedom as a
function of the temperature during RD are taken from ref. [115] and expressed
as a function of $a$ by taking $a^{3}\,T^{3}g_{\rm S}$ to be constant,
following an adiabatic expansion of the universe.

The characteristic EOS, corresponding to the energy density of equation (A.2),
can be computed combining equations (A.1) by taking the time derivative of the
second equation and introducing the first equation. This yields

$\displaystyle\frac{\dot{\Omega}(a)}{\Omega(a)}=-3(1+w)\frac{\dot{a}}{a},$
(A.3)

which allows us to compute $w(a)$ using equation (A.2),

$w(a)=\bigg{(}{1\over 3}\Omega_{{\rm
rad}}(a)-\Omega_{\Lambda,0}\bigg{)}\,\Omega^{-1}(a).$ (A.4)

Equation (A.3) justifies the expected evolution of $\Omega(a)$ used in
equation (A.2): 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 the degrees of freedom induces some modifications with respect to
the $a^{-4}$ evolution. Equation (A.4) yields $w=1/3$ for RD, $w=0$ for MD,
and $w=-1$ for $\Lambda$D, as expected. In the intermediate times, the
functions $\Omega(a)$ and $w(a)$ transition smoothly. We can express equation
(A.1) in terms of the normalized conformal time $\eta{\cal H}_{*}$ as

${\cal H}=\frac{H_{0}}{{\cal
H}_{*}}\,a\sqrt{\Omega(a)},\quad\frac{a^{\prime\prime}}{a}=\frac{1}{2}{\cal
H}^{2}\bigl{[}1-3w(a)\bigr{]},$ (A.5)

which corresponds to equation (2.9). Finally, since equation (2.5) is
expressed in terms of conformal time $\eta$ but equation (A.5) is still
written in terms of $a$, we would like to substitute the variables via $a(t)$
or $a(\eta)$, which can be obtained using equation (A.1):

$\dot{a}=H_{0}\,a\sqrt{\Omega(a)}\Rightarrow{\rm\,d}\bigl{(}H_{0}t\bigr{)}=\frac{{\rm\,d}a}{a\sqrt{\Omega(a)}}\Rightarrow
H_{0}\bigl{(}t-t_{\rm ini}\bigr{)}=\int_{a_{\rm
ini}}^{a}\frac{{\rm\,d}a}{a\sqrt{\Omega(a)}}.$ (A.6)

This allows us to compute $t(a)$ and then invert the relation to obtain
$a(t)$. Similarly, in conformal time, we have ${\rm\,d}t=a{\rm\,d}\eta$, so we
solve

$H_{0}\bigl{(}\eta-\eta_{\rm ini}\bigr{)}=\int_{a_{\rm
ini}}^{a}\frac{{\rm\,d}a}{a^{2}\sqrt{\Omega(a)}}.$ (A.7)

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}$. Further details and numerical results to Friedmann equations can be
found on the GW_turbulence GitHub project [74].

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 degrees of freedom,

$\frac{a_{*}}{a_{0}}=\frac{T_{0}}{T_{*}}\biggl{(}\frac{g_{\rm S}^{0}}{g_{\rm
S}}\biggr{)}^{1\over 3}\simeq 7.97\times
10^{-16}\frac{T_{*}}{100\,{\rm\,GeV}}\biggl{(}\frac{g_{\rm
S}}{100}\biggr{)}^{-{1\over 3}},$ (A.8)

where we take $T_{0}\simeq 2.7255$ K. Setting $a_{*}=1$, such that ${\cal
H}_{*}=H_{*}$, one gets the following value of the Hubble rate (valid during
the RD era),

${\cal H}_{*}=\frac{\pi T_{*}^{2}}{3}\sqrt{\frac{\kappa\,g_{*}}{10}}\simeq
2.066\times
10^{10}{\rm\,Hz}\,\biggl{(}\frac{T_{*}}{100{\rm\,GeV}}\biggr{)}^{2}\biggl{(}\frac{g_{*}}{100}\biggr{)}^{1\over
2},$ (A.9)

with101010During RD with $w={1\over 3}$, $\eta=1/{\cal H}$ holds exactly, but
the dependence of the radiation energy density with $g_{*}$ and $g_{\rm S}$
induce small deviations. $\eta_{*}\approx{\cal H}_{*}^{-1}$. These results
allow us to use the solutions from Friedmann equations and adapt them to
compute the variables that appear in equations (2.3) and (2.5), normalized to
the specific epoch of GW generation.

##  Appendix B GW energy spectrum via WKB approximation

In this section, we present the derivation of the GW spectrum obtained using
the WKB approximation introduced in section 3.3. We start with the WKB
solution for $\chi=e^{\mathcal{D}}h$, given in equation (3.7), and assume that
at the initial time, the solution is in the free-propagating regime and
follows GR111111We neglect here the effects of MG during the GW production and
focus on the MG effects after the GWs start to propagate. Otherwise,
$|h^{\prime}_{*}(k)|\neq k|h_{*}(k)|$ in general and the final GW spectrum
would present some additional terms. This scenario is left for future work.,
such that $|h^{\prime}_{*}(k)|=k|h_{*}(k)|$,

$\chi(k,\eta)=e^{-\tilde{\mathcal{D}}}h_{*}^{\prime}(k)\Bigg{[}\frac{1}{k}\cos
k\theta+\frac{1}{kc_{\rm T}^{*}}\biggl{(}1+\frac{{\cal A}_{*}}{k}\biggr{)}\sin
k\theta\Biggr{]},$ (B.1)

where $\theta=\tilde{c}_{\rm T}(\eta\mathcal{H}_{*}-1)$ and ${\cal
A}_{*}={\textstyle{1\over 2}}(\alpha_{\rm M}^{*}+{{c_{\rm
T}^{\prime}}^{*}\over c_{\rm T}^{*}})$. Its time derivative is

$\displaystyle\chi^{\prime}(k,\eta)=$
$\displaystyle\,e^{-\tilde{\mathcal{D}}}h_{*}^{\prime}(k)\biggl{[}\frac{c_{\rm
T}}{c_{\rm T}^{*}}\biggl{(}1+\frac{{\cal A}_{*}}{k}\biggr{)}\cos k\theta-
c_{\rm T}\sin k\theta\biggr{]}-\frac{c_{\rm T}^{\prime}}{2c_{\rm
T}}\chi(k,\eta).$ (B.2)

We have used that, for any function $f$,
$[f(k\theta)]^{\prime}=k\theta^{\prime}f^{\prime}(k\theta)=c_{\rm
T}\,k\,f(k\theta)$, since $\theta^{\prime}=\tilde{c}_{\rm T}+\tilde{c}_{\rm
T}^{\prime}(\eta{\cal H}_{*}-1)=c_{\rm T}$, and
$\tilde{\mathcal{D}}^{\prime}=c_{\rm T}^{\prime}/(2c_{\rm T})$. The spectrum
of GWs in normalized variables, given in equation (3.16), is proportional to

$\langle h^{\prime
2}(k,\eta)\rangle=e^{-2\mathcal{D}(\eta)}\bigl{\langle}\bigl{[}\chi^{\prime}(k,\eta)-\mathcal{D}^{\prime}(\eta)\chi(k,\eta)\bigr{]}^{2}\bigr{\rangle},$
(B.3)

where the brackets indicate an average over oscillations in time. Taking into
account that $\mathcal{D}^{\prime}={\textstyle{1\over 2}}\alpha_{\rm M}{\cal
H}$ and using equation (B.2), the term in brackets is

$\displaystyle\chi^{\prime}(k,\eta)-{\cal
D}^{\prime}\chi(k,\eta)=e^{-\tilde{\mathcal{D}}}h_{*}^{\prime}(k)\Bigl{[}{\cal
C}_{1}(k)\cos k\theta-{\cal C}_{2}(k)\sin k\theta\Bigr{]},$ (B.4)

where we have defined

${\cal C}_{1}(k)=\frac{c_{\rm T}}{c_{\rm T}^{*}}\biggl{(}1+\frac{{\cal
A}_{*}}{k}\biggr{)}-\frac{1}{2k}\biggl{(}\frac{c_{\rm T}^{\prime}}{c_{\rm
T}}+\alpha_{\rm M}{\cal H}\biggr{)},\quad{\cal C}_{2}(k)=c_{\rm
T}+\frac{1}{2kc_{\rm T}^{*}}\biggl{(}\frac{c_{\rm T}^{\prime}}{c_{\rm
T}}+\alpha_{\rm M}{\cal H}\biggr{)}\biggl{(}1+\frac{{\cal
A}_{*}}{k}\biggr{)}.$ (B.5)

Hence, taking the square of equation (B.4) and averaging over oscillations in
time we find

$\bigl{\langle}\bigl{[}\chi^{\prime}(k,\eta)-{\cal
D}^{\prime}\chi(k,\eta)\bigr{]}^{2}\bigr{\rangle}={\textstyle{1\over
2}}e^{-2\tilde{\mathcal{D}}}\langle
h_{*}^{\prime}(k)\,h_{*}^{\prime}(k)\rangle\Bigl{[}{\cal C}_{1}^{2}(k)+{\cal
C}_{2}^{2}(k)\Bigr{]},$ (B.6)

where

$\displaystyle{\cal C}_{1}^{2}(k)+{\cal C}_{2}^{2}(k)=$ $\displaystyle\,c_{\rm
T}^{2}\biggl{(}1+\frac{1}{{c_{\rm
T}^{*}}^{2}}\biggr{)}+\frac{1}{k^{2}}\biggl{[}\frac{c_{\rm T}^{2}{\cal
A}_{*}^{2}}{{c_{\rm T}^{*}}^{2}}+\frac{1}{4}\biggl{(}1+\frac{1}{{c_{\rm
T}^{*}}^{2}}\biggr{)}\biggl{(}\frac{c_{\rm T}^{\prime}}{c_{\rm T}}+\alpha_{\rm
M}{\cal H}\biggr{)}^{2}\biggr{]}$ $\displaystyle\,+\frac{2c_{\rm T}^{2}{\cal
A}_{*}}{k{c_{\rm T}}^{*2}}+\frac{{\cal A}_{*}}{2k^{3}{c_{\rm
T}^{*}}^{2}}\biggl{(}1+\frac{{\cal A}_{*}}{2k}\biggr{)}\biggl{(}\frac{c_{\rm
T}^{\prime}}{c_{\rm T}}+\alpha_{\rm M}{\cal H}\biggr{)}^{2}.$ (B.7)

This expression allows us to compute the GW spectrum under the MG theory from
the time evolution of $\alpha_{\rm M}$ and $c_{\rm T}$, and, in particular,
from these values at the time of GW production,

$\Omega_{\rm GW}(k,\eta)=e^{-2(\mathcal{D}+\tilde{\mathcal{D}})}\,\Omega_{\rm
GW}^{\rm GR}(k)\,\xi(k,\eta),$ (B.8)

where $\xi(k,\eta)$ characterizes the GW spectrum in MG with respect to the
one expected under GR,

$\displaystyle\xi(k,\eta)=$ $\displaystyle\,\frac{c_{\rm
T}^{2}}{2}\biggl{(}1+\frac{1}{{c_{\rm T}^{*}}^{2}}\biggr{)}+\frac{c_{\rm
T}^{2}}{2k{c_{\rm T}^{*}}^{2}}\biggl{(}\alpha_{\rm M}^{*}+\frac{{c_{\rm
T}^{\prime}}^{*}}{c_{\rm T}^{*}}\biggr{)}+\frac{1}{8k^{3}{c_{\rm
T}^{*}}^{2}}\biggl{(}\alpha_{\rm M}^{*}+\frac{{c_{\rm T}^{\prime}}^{*}}{c_{\rm
T}^{*}}\biggr{)}\biggl{(}\frac{c_{\rm T}^{\prime}}{c_{\rm T}}+\alpha_{\rm
M}{\cal H}\biggr{)}^{2}$ $\displaystyle+$
$\displaystyle\,\frac{1}{8k^{2}}\Biggl{[}\frac{c_{\rm T}^{2}}{{c_{\rm
T}^{*}}^{2}}\biggl{(}\alpha_{\rm M}^{*}+\frac{{c_{\rm T}^{\prime}}^{*}}{c_{\rm
T}^{*}}\biggr{)}^{2}+\biggl{(}1+\frac{1}{{c_{\rm
T}^{*}}^{2}}\biggr{)}\biggl{(}\frac{c_{\rm T}^{\prime}}{c_{\rm T}}+\alpha_{\rm
M}{\cal H}\biggr{)}^{2}\Biggr{]}$ $\displaystyle+$
$\displaystyle\,\frac{1}{32k^{4}{c_{\rm T}^{*}}^{2}}\biggl{(}\alpha_{\rm
M}^{*}+\frac{{c_{\rm T}^{\prime}}^{*}}{c_{\rm
T}^{*}}\biggr{)}^{2}\biggl{(}\frac{c_{\rm T}^{\prime}}{c_{\rm T}}+\alpha_{\rm
M}{\cal H}\biggr{)}^{2}.$ (B.9)

If we set $c_{\rm T}^{*}=c_{\rm T}$ and $c_{\rm T}^{\prime}={c_{\rm
T}^{\prime}}^{*}=0$, we recover equation (3.19)

$\xi(k,\eta)=1+\frac{1}{2}\alpha_{\rm T}+\frac{\alpha_{\rm
M}^{*}}{2k}+\frac{1}{8k^{2}}\Biggl{[}\alpha_{\rm M}^{2}{\cal
H}^{2}\biggl{(}1+\frac{1}{c_{\rm T}^{2}}\biggr{)}+\alpha_{\rm
M}^{*2}\Biggr{]}+\frac{\alpha_{\rm M}^{*}\alpha_{\rm M}^{2}{\cal
H}^{2}}{8k^{3}c_{\rm T}^{2}}+\frac{\alpha_{\rm M}^{*2}\alpha_{\rm M}^{2}{\cal
H}^{2}}{32k^{4}c_{\rm T}^{2}}.$ (B.10)

##  Appendix C Numerical scheme

Following refs. [73, 86], we implement a matrix solver for the modified GW
equation. Equation (2.3) can be expressed as

$h_{ij}^{\prime\prime}+\sigma h_{ij}^{\prime}+\omega^{2}h_{ij}=0,$ (C.1)

where

$\sigma\equiv\alpha_{\rm M}\mathcal{H}\,,\,\,\omega^{2}\equiv c_{\rm
T}^{2}k^{2}-\alpha_{\rm M}\mathcal{H}^{2}-\frac{a^{\prime\prime}}{a}.$ (C.2)

Then, assuming a solution of the type $h_{ij}=e_{ij}Ae^{\lambda\eta}$ gives
the characteristic equation

$\lambda^{2}+\sigma\lambda+\omega^{2}=0.$ (C.3)

The eigenvalues $\lambda$ can be obtained as

$\lambda_{1,2}=-\frac{1}{2}\Big{(}\sigma\mp\sqrt{\sigma^{2}-4\omega^{2}}\Big{)}.$
(C.4)

For $\delta\eta\ll 1$, we approximate the solution by assuming that
$\lambda_{1,2}^{\prime}\approx\lambda_{1,2}^{\prime\prime}\approx 0$ during
the time interval within time steps of the numerical solver. Then, the
solution for the strains in equation (C.1) takes the form

$\displaystyle h_{ij}(\eta+\delta\eta)$
$\displaystyle=\mathcal{C}_{ij}e^{\lambda_{1}\delta\eta}+\mathcal{D}_{ij}e^{\lambda_{2}\delta\eta},$
(C.5) $\displaystyle h_{ij}^{\prime}(\eta+\delta\eta)$
$\displaystyle=\mathcal{C}_{ij}\lambda_{1}e^{\lambda_{1}\delta\eta}+\mathcal{D}_{ij}\lambda_{2}e^{\lambda_{2}\delta\eta},$
(C.6)

where $\mathcal{C}_{ij}$ and $\mathcal{D}_{ij}$ are constant amplitude
coefficients evaluated at $\eta$. Equivalently, in matrix form, the solution
can be rewritten as

$\begin{pmatrix}h\\\
h^{\prime}\end{pmatrix}_{ij}^{\eta+\delta\eta}=\begin{pmatrix}e^{\lambda_{1}\delta\eta}&e^{\lambda_{2}\delta\eta}\\\
\lambda_{1}e^{\lambda_{1}\delta\eta}&\lambda_{2}e^{\lambda_{2}\delta\eta}\end{pmatrix}\begin{pmatrix}\mathcal{C}\\\
\mathcal{D}\end{pmatrix}_{ij}.$ (C.7)

Hence, the amplitude coefficients can be obtained as

$\displaystyle\begin{pmatrix}\mathcal{C}\\\ \mathcal{D}\end{pmatrix}_{ij}$
$\displaystyle=\begin{pmatrix}e^{\lambda_{1}\delta\eta}&e^{\lambda_{2}\delta\eta}\\\
\lambda_{1}e^{\lambda_{1}\delta\eta}&\lambda_{2}e^{\lambda_{2}\delta\eta}\end{pmatrix}^{-1}\begin{pmatrix}h\\\
h^{\prime}\end{pmatrix}_{ij}^{\eta+\delta\eta}$ (C.8)
$\displaystyle=\frac{1}{(\lambda_{2}-\lambda_{1})e^{(\lambda_{1}+\lambda_{2})\delta\eta}}\begin{pmatrix}\lambda_{2}e^{\lambda_{2}\delta\eta}&-e^{\lambda_{2}\delta\eta}\\\
-\lambda_{1}e^{\lambda_{1}\delta\eta}&e^{\lambda_{1}\delta\eta}\end{pmatrix}\begin{pmatrix}h\\\
h^{\prime}\end{pmatrix}_{ij}^{\eta+\delta\eta},$ (C.9)

and, in the limit $\delta\eta\rightarrow 0$, we get

$\displaystyle\lim_{\delta\eta\rightarrow 0}\begin{pmatrix}\mathcal{C}\\\
\mathcal{D}\end{pmatrix}_{ij}$
$\displaystyle=\frac{1}{\lambda_{2}-\lambda_{1}}\begin{pmatrix}\lambda_{2}&-1\\\
-\lambda_{1}&1\end{pmatrix}\begin{pmatrix}h\\\
h^{\prime}\end{pmatrix}_{ij}^{\eta},$ (C.10)

where we assume that $\mathcal{C}_{ij}$ and $\mathcal{D}_{ij}$ are time-
independent within time steps. Therefore, the time evolution of the relevant
quantities can be obtained via a matrix multiplication as

$\begin{pmatrix}h\\\
h^{\prime}\end{pmatrix}_{ij}^{\eta+\delta\eta}=\mathcal{M}\begin{pmatrix}h\\\
h^{\prime}\end{pmatrix}_{ij}^{\eta},$ (C.11)

where

$\mathcal{M}=\frac{1}{\lambda_{1}-\lambda_{2}}\begin{pmatrix}\lambda_{1}e^{\lambda_{2}\delta\eta}-\lambda_{2}e^{\lambda_{1}\delta\eta}&e^{\lambda_{1}\delta\eta}-e^{\lambda_{2}\delta\eta}\\\
\lambda_{1}\lambda_{2}(e^{\lambda_{2}\delta\eta}-e^{\lambda_{1}\delta\eta})&\lambda_{1}e^{\lambda_{1}\delta\eta}-\lambda_{2}e^{\lambda_{2}\delta\eta}\end{pmatrix}.$
(C.12)

In the case of $\sigma=0$ and $\omega^{2}=k^{2}$, the eigenvalues become
$\lambda_{1,2}=\pm ik$ and the evolution matrix reduces to the one relevant
for GR,

$\mathcal{M}=\begin{pmatrix}\cos{k\delta\eta}&\sin{k\delta\eta}\\\
-\sin{k\delta\eta}&\cos{k\delta\eta}\end{pmatrix}.$ (C.13)

##  Appendix D Numerical accuracy

In this section, we study the numerical accuracy of the runs. The GW solver in
the Pencil Code has been made accurate to second order in time in the presence
of a source $T_{ij}$ under GR. However, the GW equation under modified
theories of gravity lead to an equation that can be approximated within small
but finite time steps under the assumption that the coefficients are constant
in time, leading to the numerical solver described in section C. This method
is accurate up to first order in the duration of the time step. As the
coefficient variations in time become smaller, the error is also smaller, but
the accuracy of the solver is still of first order. This is demonstrated in
figure 10, where we show on the left panel the relative errors in the GW
energy during RD and MD, and on the right panel the error during $\Lambda$D.
Recall that, as introduced in section 5.1, during RD and MD, we use increasing
time steps, represented by $n_{\rm incr}$ here. During $\Lambda$D, we use
linear time steps, with the interval shown as $\eta_{\rm incr}$ in figure 10.

Figure 10: Relative errors as a function of the numerical time step. Left
panel: Accuracy during RD (round dots) and MD (square dots) for $\alpha_{\rm
M,0}=-0.5$ (blue) and $\alpha_{\rm M,0}=-0.01$ (green), where $\delta n_{\rm
incr}$ indicates the difference in the time steps such that
$\delta\eta=\eta\,\delta n_{\rm incr}$. Right panel: Accuracy during
$\Lambda$D with $\delta\eta=\delta\eta_{\rm incr}$ directly indicating the
linear time interval. The errors are calculated relative to the values
obtained from using the smallest $\delta n_{\rm incr}=10^{-2}$ (left panel)
and $\delta\eta_{\rm incr}=5\times 10^{9}$ (right panel). Figure 11: Saturated
energy spectrum $e^{2\mathcal{D}}S_{h^{\prime}}(k)$ for $\alpha_{\rm
M,0}=-0.5$ (run M0A) with different time steps $\delta\eta\,\mathcal{H}_{*}$
during $\Lambda$D, where $n_{\rm incr}=0.01$ is used during RD and MD. The
black dashed curves correspond to the run M0A performed entirely, i.e.,
including $\Lambda$D, with increasing time steps with $n_{\rm incr}=0.01$. The
inset directly shows the spectrum at the final time without averaging over
time oscillations.

We also test the accuracy of the spectra in figure 11, where we show the
saturated final energy spectrum for run M0A. Different colors indicate the
linear time intervals during $\Lambda$D and the black curve is the spectrum
obtained by running M0A entirely with increasing time steps using $n_{\rm
incr}=0.01$. In other words, the colored curves only differ from the black one
during $\Lambda$D. We observe that all of the chosen values of $\delta\eta$
produce converging energy spectra after averaging over oscillations in time,
although larger values of $\delta\eta$ result in more fluctuations in the
numerical spectra (see the inset in figure 11). We also note that, compared to
results from hybrid time steps, the entirely nonuniform time steps somewhat
underestimate the final energy spectrum.

These accuracy tests justify our choice to use nonuniform time steps during RD
and MD to improve the efficiency of the simulations, and uniform linear time
steps (with $\delta\eta$ that already gives converged results) during
$\Lambda$D.

## References

  * [1] A.G. Riess, S. Casertano, W. Yuan, L.M. Macri and D. Scolnic, _Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM_, _Astrophys. J._ 876 (2019) 85 [1903.07603].
  * [2] A.G. Riess, S. Casertano, W. Yuan, J.B. Bowers, L. Macri, J.C. Zinn et al., _Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM_, _Astrophys. J. Lett._ 908 (2021) L6 [2012.08534].
  * [3] A.G. Riess et al., _A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s -1 Mpc-1 Uncertainty from the Hubble Space Telescope and the SHOES Team_, _Astrophys. J. Lett._ 934 (2022) L7 [2112.04510].
  * [4] K.C. Wong et al., _H0LiCOW – XIII. A 2.4 per cent measurement of H0 from lensed quasars: 5.3 $\sigma$ tension between early- and late-Universe probes_, _Mon. Not. Roy. Astron. Soc._ 498 (2020) 1420 [1907.04869].
  * [5] D.W. Pesce et al., _The Megamaser Cosmology Project. XIII. Combined Hubble constant constraints_ , _Astrophys. J. Lett._ 891 (2020) L1 [2001.09213].
  * [6] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [7] BOSS collaboration, _The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample_ , _Mon. Not. Roy. Astron. Soc._ 470 (2017) 2617 [1607.03155].
  * [8] L. Verde, T. Treu and A.G. Riess, _Tensions between the Early and the Late Universe_ , _Nature Astron._ 3 (2019) 891 [1907.10625].
  * [9] B.F. Schutz, _Determining the Hubble Constant from Gravitational Wave Observations_ , _Nature_ 323 (1986) 310.
  * [10] D.E. Holz and S.A. Hughes, _Using gravitational-wave standard sirens_ , _Astrophys. J._ 629 (2005) 15 [astro-ph/0504616].
  * [11] S. Nissanke, D.E. Holz, N. Dalal, S.A. Hughes, J.L. Sievers and C.M. Hirata, _Determining the Hubble constant from gravitational wave observations of merging compact binaries_ , 1307.2638.
  * [12] LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER collaboration, _A gravitational-wave standard siren measurement of the Hubble constant_ , _Nature_ 551 (2017) 85 [1710.05835].
  * [13] DES collaboration, _A statistical standard siren measurement of the Hubble constant from the LIGO/Virgo gravitational wave compact object merger GW190814 and Dark Energy Survey galaxies_ , _Astrophys. J. Lett._ 900 (2020) L33 [2006.14961].
  * [14] A. Palmese, C.R. Bom, S. Mucesh and W.G. Hartley, _A Standard Siren Measurement of the Hubble Constant Using Gravitational-wave Events from the First Three LIGO/Virgo Observing Runs and the DESI Legacy Survey_ , _Astrophys. J._ 943 (2023) 56 [2111.06445].
  * [15] LISA Cosmology Working Group collaboration, _Cosmology with the Laser Interferometer Space Antenna_ , 2204.05434.
  * [16] M. Branchesi et al., _Science with the Einstein Telescope: a comparison of different designs_ , 2303.15923.
  * [17] S. Joudaki et al., _KiDS-450 + 2dFLenS: Cosmological parameter constraints from weak gravitational lensing tomography and overlapping redshift-space galaxy clustering_ , _Mon. Not. Roy. Astron. Soc._ 474 (2018) 4894 [1707.06627].
  * [18] L. Perivolaropoulos and F. Skara, _Challenges for $\Lambda$CDM: An update_, _New Astron. Rev._ 95 (2022) 101659 [2105.05208].
  * [19] D. Lovelock, _The Einstein tensor and its generalizations_ , _J. Math. Phys._ 12 (1971) 498.
  * [20] D. Lovelock, _The four-dimensionality of space and the einstein tensor_ , _J. Math. Phys._ 13 (1972) 874.
  * [21] G.W. Horndeski, _Second-order scalar-tensor field equations in a four-dimensional space_ , _Int. J. Theor. Phys._ 10 (1974) 363.
  * [22] C. Deffayet, X. Gao, D.A. Steer and G. Zahariade, _From k-essence to generalised Galileons_ , _Phys. Rev. D_ 84 (2011) 064039 [1103.3260].
  * [23] T. Kobayashi, M. Yamaguchi and J. Yokoyama, _Generalized G-inflation: Inflation with the most general second-order field equations_ , _Prog. Theor. Phys._ 126 (2011) 511 [1105.5723].
  * [24] M. Zumalacárregui and J. García-Bellido, _Transforming gravity: from derivative couplings to matter to second-order scalar-tensor theories beyond the Horndeski Lagrangian_ , _Phys. Rev. D_ 89 (2014) 064046 [1308.4685].
  * [25] J. Gleyzes, D. Langlois, F. Piazza and F. Vernizzi, _Healthy theories beyond Horndeski_ , _Phys. Rev. Lett._ 114 (2015) 211101 [1404.6495].
  * [26] J. Gleyzes, D. Langlois, F. Piazza and F. Vernizzi, _Exploring gravitational theories beyond Horndeski_ , _JCAP_ 02 (2015) 018 [1408.1952].
  * [27] D. Langlois and K. Noui, _Degenerate higher derivative theories beyond Horndeski: evading the Ostrogradski instability_ , _JCAP_ 02 (2016) 034 [1510.06930].
  * [28] J. Ben Achour, M. Crisostomi, K. Koyama, D. Langlois, K. Noui and G. Tasinato, _Degenerate higher order scalar-tensor theories beyond Horndeski up to cubic order_ , _JHEP_ 12 (2016) 100 [1608.08135].
  * [29] D. Langlois, M. Mancarella, K. Noui and F. Vernizzi, _Effective Description of Higher-Order Scalar-Tensor Theories_ , _JCAP_ 05 (2017) 033 [1703.03797].
  * [30] P. Hořava, _Quantum Gravity at a Lifshitz Point_ , _Phys. Rev. D_ 79 (2009) 084008 [0901.3775].
  * [31] D. Blas, O. Pujolas and S. Sibiryakov, _Models of non-relativistic quantum gravity: The Good, the bad and the healthy_ , _JHEP_ 04 (2011) 018 [1007.3503].
  * [32] T. Jacobson and D. Mattingly, _Gravity with a dynamical preferred frame_ , _Phys. Rev. D_ 64 (2001) 024028 [gr-qc/0007031].
  * [33] T.G. Zlosnik, P.G. Ferreira and G.D. Starkman, _Modifying gravity with the Aether: An alternative to Dark Matter_ , _Phys. Rev. D_ 75 (2007) 044017 [astro-ph/0607411].
  * [34] G. Tasinato, _Cosmic Acceleration from Abelian Symmetry Breaking_ , _JHEP_ 04 (2014) 067 [1402.6450].
  * [35] L. Heisenberg, _Generalization of the Proca Action_ , _JCAP_ 05 (2014) 015 [1402.7026].
  * [36] L. Heisenberg, R. Kase and S. Tsujikawa, _Beyond generalized Proca theories_ , _Phys. Lett. B_ 760 (2016) 617 [1605.05565].
  * [37] E. Allys, P. Peter and Y. Rodriguez, _Generalized Proca action for an Abelian vector field_ , _JCAP_ 02 (2016) 004 [1511.03101].
  * [38] J. Beltran Jimenez and L. Heisenberg, _Derivative self-interactions for a massive vector field_ , _Phys. Lett. B_ 757 (2016) 405 [1602.03410].
  * [39] C. de Rham and G. Gabadadze, _Generalization of the Fierz-Pauli Action_ , _Phys. Rev. D_ 82 (2010) 044020 [1007.0443].
  * [40] C. de Rham, G. Gabadadze and A.J. Tolley, _Resummation of Massive Gravity_ , _Phys. Rev. Lett._ 106 (2011) 231101 [1011.1232].
  * [41] S.F. Hassan and R.A. Rosen, _Bimetric Gravity from Ghost-free Massive Gravity_ , _JHEP_ 02 (2012) 126 [1109.3515].
  * [42] E. Sagi, _Propagation of gravitational waves in generalized TeVeS_ , _Phys. Rev. D_ 81 (2010) 064031 [1001.1555].
  * [43] C. Brans and R.H. Dicke, _Mach’s principle and a relativistic theory of gravitation_ , _Phys. Rev._ 124 (1961) 925.
  * [44] B. Ratra and P.J.E. Peebles, _Cosmological Consequences of a Rolling Homogeneous Scalar Field_ , _Phys. Rev. D_ 37 (1988) 3406.
  * [45] A. De Felice and S. Tsujikawa, _f(R) theories_ , _Living Rev. Rel._ 13 (2010) 3 [1002.4928].
  * [46] S.M. Carroll, A. De Felice, V. Duvvuri, D.A. Easson, M. Trodden and M.S. Turner, _The cosmology of generalized modified gravity models_ , _Phys. Rev. D_ 71 (2005) 063513 [astro-ph/0410031].
  * [47] C. Armendariz-Picon, V.F. Mukhanov and P.J. Steinhardt, _Essentials of k essence_ , _Phys. Rev. D_ 63 (2001) 103510 [astro-ph/0006373].
  * [48] C. Deffayet, O. Pujolas, I. Sawicki and A. Vikman, _Imperfect Dark Energy from Kinetic Gravity Braiding_ , _JCAP_ 10 (2010) 026 [1008.0048].
  * [49] A. Nicolis, R. Rattazzi and E. Trincherini, _The Galileon as a local modification of gravity_ , _Phys. Rev. D_ 79 (2009) 064036 [0811.2197].
  * [50] T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, _Modified Gravity and Cosmology_ , _Phys. Rept._ 513 (2012) 1 [1106.2476].
  * [51] E. Belgacem, Y. Dirian, S. Foffa and M. Maggiore, _Modified gravitational-wave propagation and standard sirens_ , _Phys. Rev. D_ 98 (2018) 023510 [1805.08731].
  * [52] C. Deffayet and K. Menou, _Probing Gravity with Spacetime Sirens_ , _Astrophys. J. Lett._ 668 (2007) L143 [0709.0003].
  * [53] I.D. Saltas, I. Sawicki, L. Amendola and M. Kunz, _Anisotropic Stress as a Signature of Nonstandard Propagation of Gravitational Waves_ , _Phys. Rev. Lett._ 113 (2014) 191101 [1406.7139].
  * [54] L. Lombriser and A. Taylor, _Breaking a Dark Degeneracy with Gravitational Waves_ , _JCAP_ 03 (2016) 031 [1509.08458].
  * [55] A. Nishizawa, _Generalized framework for testing gravity with gravitational-wave propagation. I. Formulation_ , _Phys. Rev. D_ 97 (2018) 104037 [1710.04825].
  * [56] E. Belgacem, Y. Dirian, S. Foffa and M. Maggiore, _Gravitational-wave luminosity distance in modified gravity theories_ , _Phys. Rev. D_ 97 (2018) 104066 [1712.08108].
  * [57] S. Mastrogiovanni, D. Steer and M. Barsuglia, _Probing modified gravity theories and cosmology using gravitational-waves and associated electromagnetic counterparts_ , _Phys. Rev. D_ 102 (2020) 044009 [2004.01632].
  * [58] S. Mastrogiovanni, L. Haegel, C. Karathanasis, I.M.n. Hernandez and D.A. Steer, _Gravitational wave friction in light of GW170817 and GW190521_ , _JCAP_ 02 (2021) 043 [2010.04047].
  * [59] S. Arai and A. Nishizawa, _Generalized framework for testing gravity with gravitational-wave propagation. II. Constraints on Horndeski theory_ , _Phys. Rev. D_ 97 (2018) 104038 [1711.03776].
  * [60] LISA Cosmology Working Group collaboration, _Testing modified gravity at cosmological distances with LISA standard sirens_ , _JCAP_ 07 (2019) 024 [1906.01593].
  * [61] R. D’Agostino and R.C. Nunes, _Probing observational bounds on scalar-tensor theories from standard sirens_ , _Phys. Rev. D_ 100 (2019) 044041 [1907.05516].
  * [62] I.S. Matos, M.O. Calvão and I. Waga, _Gravitational wave propagation in $f(R)$ models: New parametrizations and observational constraints_, _Phys. Rev. D_ 103 (2021) 104059 [2104.10305].
  * [63] I.S. Matos, E. Bellini, M.O. Calvão and M. Kunz, _Testing gravity with gravitational wave friction and gravitational slip_ , 2210.12174.
  * [64] LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL collaboration, _Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A_ , _Astrophys. J. Lett._ 848 (2017) L13 [1710.05834].
  * [65] T. Baker, E. Bellini, P.G. Ferreira, M. Lagos, J. Noller and I. Sawicki, _Strong constraints on cosmological gravity from GW170817 and GRB 170817A_ , _Phys. Rev. Lett._ 119 (2017) 251301 [1710.06394].
  * [66] J.M. Ezquiaga and M. Zumalacárregui, _Dark Energy After GW170817: Dead Ends and the Road Ahead_ , _Phys. Rev. Lett._ 119 (2017) 251304 [1710.05901].
  * [67] P. Creminelli and F. Vernizzi, _Dark Energy after GW170817 and GRB170817A_ , _Phys. Rev. Lett._ 119 (2017) 251302 [1710.05877].
  * [68] J. Sakstein and B. Jain, _Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories_ , _Phys. Rev. Lett._ 119 (2017) 251303 [1710.05893].
  * [69] C. de Rham and S. Melville, _Gravitational Rainbows: LIGO and Dark Energy at its Cutoff_ , _Phys. Rev. Lett._ 121 (2018) 221101 [1806.09417].
  * [70] Y. He, A. Brandenburg and A. Sinha, _Tensor spectrum of turbulence-sourced gravitational waves as a constraint on graviton mass_ , _JCAP_ 07 (2021) 015 [2104.03192].
  * [71] M. Chevallier and D. Polarski, _Accelerating universes with scaling dark matter_ , _Int. J. Mod. Phys. D_ 10 (2001) 213 [gr-qc/0009008].
  * [72] Pencil Code collaboration, _The Pencil Code, a modular MPI code for partial differential equations and particles: multipurpose and multiuser-maintained_ , _J. Open Source Softw._ 6 (2021) 2807 [2009.08231].
  * [73] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky and S. Mandal, _The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence_ , _Geophys. Astrophys. Fluid Dynamics_ 114 (2020) 130 [1807.05479].
  * [74] A. Roper Pol, _GitHub project “GW_turbulence”_, _doi: 10.5281/zenodo.6045844_ .
  * [75] C. Caprini and D.G. Figueroa, _Cosmological Backgrounds of Gravitational Waves_ , _Class. Quant. Grav._ 35 (2018) 163001 [1801.04268].
  * [76] J.-C. Hwang and H.-R. Noh, _Gauge ready formulation of the cosmological kinetic theory in generalized gravity theories_ , _Phys. Rev. D_ 65 (2002) 023512 [astro-ph/0102005].
  * [77] R.C. Nunes, M.E.S. Alves and J.C.N. de Araujo, _Primordial gravitational waves in Horndeski gravity_ , _Phys. Rev. D_ 99 (2019) 084022 [1811.12760].
  * [78] G. Tasinato, A. Garoffolo, D. Bertacca and S. Matarrese, _Gravitational-wave cosmological distances in scalar-tensor theories of gravity_ , _JCAP_ 06 (2021) 050 [2103.00155].
  * [79] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, _Numerical simulations of gravitational waves from early-universe turbulence_ , _Phys. Rev. D_ 102 (2020) 083512 [1903.08585].
  * [80] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal and A. Roper Pol, _Circular polarization of gravitational waves from early-Universe helical turbulence_ , _Phys. Rev. Res._ 3 (2021) 013193 [2011.05556].
  * [81] A. Brandenburg, G. Gogoberidze, T. Kahniashvili, S. Mandal, A. Roper Pol and N. Shenoy, _The scalar, vector, and tensor modes in gravitational wave turbulence simulations_ , _Class. Quant. Grav._ 38 (2021) 145002 [2103.01140].
  * [82] A. Roper Pol, S. Mandal, A. Brandenburg and T. Kahniashvili, _Polarization of gravitational waves from helical MHD turbulent sources_ , _JCAP_ 04 (2022) 019 [2107.05356].
  * [83] A. Brandenburg, E. Clarke, Y. He and T. Kahniashvili, _Can we observe the QCD phase transition-generated gravitational waves through pulsar timing arrays?_ , _Phys. Rev. D_ 104 (2021) 043513 [2102.12428].
  * [84] A. Roper Pol, C. Caprini, A. Neronov and D. Semikoz, _Gravitational wave signal from primordial magnetic fields in the Pulsar Timing Array frequency band_ , _Phys. Rev. D_ 105 (2022) 123502 [2201.05630].
  * [85] R. Sharma and A. Brandenburg, _Low frequency tail of gravitational wave spectra from hydromagnetic turbulence_ , _Phys. Rev. D_ 106 (2022) 103536 [2206.00055].
  * [86] A. Brandenburg and R. Sharma, _Simulating Relic Gravitational Waves from Inflationary Magnetogenesis_ , _Astrophys. J._ 920 (2021) 26 [2106.03857].
  * [87] A. Brandenburg, Y. He and R. Sharma, _Simulations of Helical Inflationary Magnetogenesis and Gravitational Waves_ , _Astrophys. J._ 922 (2021) 192 [2107.12333].
  * [88] Y. He, A. Roper Pol and A. Brandenburg, _Leading-order nonlinear gravitational waves from reheating magnetogeneses_ , 2110.14456.
  * [89] K. Saikawa and S. Shirai, _Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model_ , _JCAP_ 05 (2018) 035 [1803.01038].
  * [90] LISA Cosmology Working Group collaboration, _Measuring the propagation speed of gravitational waves with LISA_ , _JCAP_ 08 (2022) 031 [2203.00566].
  * [91] Y. Cai, Y.-T. Wang and Y.-S. Piao, _Is there an effect of a nontrivial $c_{T}$ during inflation?_, _Phys. Rev. D_ 93 (2016) 063005 [1510.08716].
  * [92] Y. Cai, Y.-T. Wang and Y.-S. Piao, _Propagating speed of primordial gravitational waves and inflation_ , _Phys. Rev. D_ 94 (2016) 043002 [1602.05431].
  * [93] S.D. Odintsov, V.K. Oikonomou and R. Myrzakulov, _Spectrum of Primordial Gravitational Waves in Modified Gravities: A Short Overview_ , _Symmetry_ 14 (2022) 729 [2204.00876].
  * [94] E. Bellini and I. Sawicki, _Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity_ , _JCAP_ 07 (2014) 050 [1404.3713].
  * [95] J. Kennedy, L. Lombriser and A. Taylor, _Reconstructing Horndeski theories from phenomenological modified gravity and dark energy models on cosmological scales_ , _Phys. Rev. D_ 98 (2018) 044051 [1804.04582].
  * [96] J. Gleyzes, _Parametrizing modified gravity for cosmological surveys_ , _Phys. Rev. D_ 96 (2017) 063516 [1705.04714].
  * [97] M. Denissenya and E.V. Linder, _Gravity’s Islands: Parametrizing Horndeski Stability_ , _JCAP_ 11 (2018) 010 [1808.00013].
  * [98] M. Yoshimura, _Stronger gravity in the early universe_ , 2204.11384.
  * [99] M. Zumalacarregui, _Gravity in the Era of Equality: Towards solutions to the Hubble problem without fine-tuned initial conditions_ , _Phys. Rev. D_ 102 (2020) 023523 [2003.06396].
  * [100] M. Braglia, M. Ballardini, F. Finelli and K. Koyama, _Early modified gravity in light of the $H_{0}$ tension and LSS data_, _Phys. Rev. D_ 103 (2021) 043528 [2011.12934].
  * [101] E.V. Linder, _Gravitational Wave Distances in Horndeski Cosmology_ , 2108.11526.
  * [102] M.S. Turner, E.J. Weinberg and L.M. Widrow, _Bubble nucleation in first order inflation and other cosmological phase transitions_ , _Phys. Rev. D_ 46 (1992) 2384.
  * [103] NANOGrav collaboration, _The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background_ , _Astrophys. J. Lett._ 905 (2020) L34 [2009.04496].
  * [104] B. Goncharov et al., _On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 917 (2021) L19 [2107.12112].
  * [105] S. Chen et al., _Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search_ , _Mon. Not. Roy. Astron. Soc._ 508 (2021) 4970 [2110.13184].
  * [106] J. Antoniadis et al., _The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background_ , _Mon. Not. Roy. Astron. Soc._ 510 (2022) 4873 [2201.03980].
  * [107] K. Schmitz, _New Sensitivity Curves for Gravitational-Wave Signals from Cosmological Phase Transitions_ , _JHEP_ 01 (2021) 097 [2002.04615].
  * [108] C.J. Moore, R.H. Cole and C.P.L. Berry, _Gravitational-wave sensitivity curves_ , _Class. Quant. Grav._ 32 (2015) 015014 [1408.0740].
  * [109] N. Seto, S. Kawamura and T. Nakamura, _Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space_ , _Phys. Rev. Lett._ 87 (2001) 221103 [astro-ph/0108011].
  * [110] J. Crowder and N.J. Cornish, _Beyond LISA: Exploring future gravitational wave missions_ , _Phys. Rev. D_ 72 (2005) 083005 [gr-qc/0506015].
  * [111] C. Caprini, D.G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni et al., _Reconstructing the spectral shape of a stochastic gravitational wave background with LISA_ , _JCAP_ 11 (2019) 017 [1906.09244].
  * [112] M. Punturo et al., _The Einstein Telescope: A third-generation gravitational wave observatory_ , _Class. Quant. Grav._ 27 (2010) 194002.
  * [113] A. Neronov, A. Roper Pol, C. Caprini and D. Semikoz, _NANOGrav signal from magnetohydrodynamic turbulence at the QCD phase transition in the early Universe_ , _Phys. Rev. D_ 103 (2021) 041302 [2009.14174].