# Production of Primordial Gravitational Waves in Teleparallel Gravity

Geovanny A. Rave-Franco<EMAIL_ADDRESS>Instituto de Ciencias
Nucleares, Universidad Nacional Autónoma de México, Circuito Exterior C.U.,
A.P. 70-543, México D.F. 04510, México.  Celia Escamilla-Rivera
<EMAIL_ADDRESS>Instituto de Ciencias Nucleares, Universidad
Nacional Autónoma de México, Circuito Exterior C.U., A.P. 70-543, México D.F.
04510, México.  Jackson Levi Said<EMAIL_ADDRESS>Institute of Space
Sciences and Astronomy, University of Malta, Malta, MSD 2080. Department of
Physics, University of Malta, Malta, MSD 2080.

###### Abstract

We study the production of primordial gravitational waves in the context of
extended teleparallel gravity models and compare them with those of general
relativity. Teleparallel gravity has been widely studied in the context of the
late universe but not much in the early universe. Two sources of primordial
gravitational waves are considered, vacuum fluctuations and tensor
anisotropies within two inflation-compatible backgrounds: a perfect de Sitter
and a quasi de Sitter background. We find that in the vacuum case with a
perfect de Sitter background, the gravitational waves propagation equation is
the same as that of general relativity, however, if the background is promoted
to a quasi de Sitter background, the propagation equations are different
resulting in an important difference on the tensor spectral index. When tensor
anisotropies are included, we compute the most general solution for
gravitational waves in terms of a retarded Green’s function and analyze the
contributions to the power spectrum from these anisotropies. Finally, we
investigate the energy density of these gravitational waves.

##  1 Introduction

General Relativity (GR) has been a successful theory to explain different
cosmological phenomena at different cosmological scales in the light of
current observational data within the framework of the flat $\Lambda$CDM
concordance model [1, 2], which can explain the cosmological dynamics at
different stages, particularly, the formation of large-scale structure [3],
thermal history [4], the current accelerated expansion [5], and a wide range
of other cosmological observations [6, 7, 8]. However, despite its success,
this model faces several theoretical and observational issues, namely, the
fine-tuning problem [9, 10], coincidence problem [11, 12], and currently, a
growing tension between direct and indirect measurements of some prominent
cosmological parameters [13, 14].

Several approaches to solve these problems have been made, some of them
include scalar quintessence models [15, 16], high dimensional theories [17], a
variable cosmological constant [18], and extensions of gravity [19, 20] as
well as a reexamination of the cosmological principle [21, 22], possible early
Universe dark energy [23], and extra degrees of freedom for neutrino species
in the early Universe [24, 25].

In the direction of the latter idea, Teleparallel Gravity (TG) is a gauge
theory of gravity where the gravitational field is mediated by torsion instead
of curvature [26]. Here, a teleparallel equivalent to general relativity
(TEGR) can be defined in which all classical phenomena are equivalently
predicted, making them dynamically equivalent at this level [27, 28, 29].
Naturally, this can be extended in similar ways that we do in the regular
curvature-based setting. The look for $f(T)$ and $f(T,B)$ extended models in
TG have been done to address these issues in the context of the late universe,
showing promising in explaining the late time acceleration as a consequence of
the geometry of each model [30, 31, 32], alleviating to some extent the
$H_{0}$ tension [33, 34], the structure formation and inhomogeneities [35,
36], among other topics (see Ref. [27]).

On the other hand, there have been a few studies made in the context of the
early universe, some of them include the study of primordial black holes [37],
constraints from inflation [38] and the study of inflation in the context of
TG [39, 40]. An important effect of the early universe is primordial
gravitational waves (GW) since they represent a window to explore the physics
of the inflationary era and the first stages of the universe [41]. The problem
of primordial GW in TG was partially addressed in [37] but no investigation on
the power spectrum was performed in this research.

In this work, we will give a first look into the generation of GW as linear
tensor perturbations around a flat FLRW background, in the early universe,
that is primordial GW. Particularly, we will obtain the solution in TG for GW
coming from vacuum fluctuations and the observational implications on the
power spectrum and the tensor spectral index. We will also derive the general
solution for the GW including an arbitrary possible source of tensor
anisotropic stress and see how the power spectrum is changed by the presence
of extra geometry given by the non-linearity of the extended $f(T)$ and
$f(T,B)$ models, in comparison with GR. As a further step, we will discuss the
observational and theoretical implications of GW energy density.

This work is divided as follows: In Sec. 2 we will discuss the basic elements
of TG, its geometrical setup, and the background cosmology on both $f(T)$ and
$f(T,B)$ models. In Sec. 3 we discuss the linear perturbation theory in TG
from a transverse and traceless tensor perturbation, we will study the
production of GW coming from vacuum fluctuations and tensor anisotropies.
Finally, in Sec. 5 we discuss the principal results of this work and its
observational implications, and in Sec. 6 we present our conclusions. In this
work, Greek indices $\mu,\nu,\ldots$ refer to spacetime indices, and capital
Latin indices $A,B,\ldots$ refer to a tangent or Minkowski spacetime indices.

##  2 Basics of teleparallel gravity and its cosmology

TG is a gauge theory of translations invariant under local Lorentz
transformations [26]. The geometrical setup of TG consists of two fundamental
objects [42], the tetrad fields
$\bm{e}\in\Omega^{1}(\mathcal{M},\mathbb{R}^{1,3})$, a set differential
$1$-forms on $M$ assuming values on the Minkowski space
$\bm{e}=\\{e^{A}\\}^{3}_{A=0}$, and the spin connection
$\bm{\omega}\in\Omega^{1}(\mathcal{M},\mathfrak{so}(1,3))$ a differential
$1$-form assuming values on the Lie Algebra of the Lorentz Group. On a local
coordinate system $U\subset\mathcal{M}$ of the base manifold $\mathcal{M}$,
these objects read as

$e^{A}_{x}=e^{A}_{\
\mu}(x)dx^{\mu}\,,\quad\bm{\omega}=\frac{1}{2}\omega^{AB}_{\ \
\mu}dx^{\mu}\otimes
S^{AB}\in\Gamma\left(T^{*}U\otimes\mathfrak{so}(1,3)\right)\,,\quad
x\mapsto(x,\bm{\omega}_{x})\,,$ (1)

where $S^{AB}$ are the generators of the Lorentz group, $\mathfrak{so}(1,3)$
the Lie Algebra of the Lorentz group, $\Gamma(T^{*}U)$ the sections over the
cotangent bundle and $\bm{\omega}_{x}=\frac{1}{2}\omega^{AB}_{\ \
\mu}(x)dx^{\mu}\otimes S^{AB}\in T^{*}U\otimes\mathfrak{so}(1,3)$ [43, 44].
However, in TG the tetrad fields and the spin connection are not any arbitrary
fields satisfying the aforementioned definitions but particular fields as we
shall see later on in this section. The tetrad fields $e^{A}=e^{A}_{\
\mu}dx^{\mu}$ and $E_{A}=E_{A}^{\ \mu}\partial_{\mu}$ constitute the non-
coordinate basis of the sections over the tangent bundle and cotangent bundles
$\Gamma(TU)$ and $\Gamma(T^{*}U)$ in every local coordinate system $U$ of the
manifold $\mathcal{M}$ such that

$e^{A}_{\ \mu}E_{A}^{\ \nu}=\delta_{\mu}^{\nu}\,,\quad\text{and}\quad e^{A}_{\
\mu}E_{B}^{\ \mu}=\delta_{A}^{B}\,,$ (2)

and at every point of the base manifold $p\in\mathcal{M}$ they diagonalize the
metric [45]

$g_{AB}(p)=\eta_{AB}(p)=g_{\mu\nu}(p)E^{\mu}_{A}(p)E^{\nu}_{B}(p)\,,$ (3)

or analogously,

$g_{\mu\nu}(p)=e^{A}_{\ \mu}(p)e^{B}_{\ \nu}(p)\eta_{AB}(p)\,.$ (4)

Under diffeomorphism over the base manifold $\mathcal{M}$ and local Lorentz
transformation $\Lambda^{A}_{\ B}(\bm{x})\in SO^{+}(1,3)$ over the Minkowski
space, the tetrad fields transform as

$e^{A}_{\ \mu^{\prime}}=\frac{\partial x^{\nu}}{\partial
x^{\mu^{\prime}}}e^{A}_{\ \nu}\,,\quad\text{and}\quad e^{A}_{\
\mu}=\Lambda^{A}_{\ B}e^{B}_{\ \mu}\,.$ (5)

In the case of the spin connection, we aim to describe a globally flat
manifold but with non-trivial geometry given by the torsion, hence, we require
the curvature $2$-form to satisfy

$\mathbf{R}=\frac{1}{4}R^{A}_{\text{ }B\nu\mu}S_{A}^{B}dx^{\nu}\wedge
dx^{\mu}\quad\text{with}\quad R^{A}_{\text{
}B\nu\mu}=\partial_{\nu}\omega^{A}_{\text{
}B\mu}-\partial_{\mu}\omega^{A}_{\text{ }B\nu}+\omega^{A}_{\text{
}D\nu}\omega^{D}_{\text{ }B\mu}-\omega^{A}_{\text{ }D\mu}\omega^{D}_{\text{
}B\nu}\equiv 0\,,$ (6)

which is achieved by

$\omega^{A}_{\ B\mu}=\Lambda^{A}_{\ C}(x)\partial_{\mu}\Lambda^{C}_{\
B}(x)\,,$ (7)

called the purely inertial spin connection, accounting for the inertial
effects of the theory. The purely inertial spin connection is the portrayal of
a zero spin connection in an arbitrary Lorentz frame

${\omega^{A}_{\ B\mu}}(\bm{x})=\Lambda^{A}_{\ C}(\bm{x}){\omega^{C}_{\
D\mu}}^{\prime}(\bm{x})\Lambda_{B}^{\ D}(\bm{x})+\Lambda^{A}_{\
C}(\bm{x})\partial_{\mu}\Lambda^{C}_{\
B}(\bm{x})\quad\text{with}\quad{\omega^{C}_{\ D\mu}}^{\prime}=0\,,$ (8)

and thus can be transformed into a zero spin connection using the inverse
Lorentz transformation [42]. The process of choosing a vanishing spin
connection ${\omega^{A}_{\ B\mu}}=0$ is called the Weitzenböck gauge, and can
only be applied when the spin connection field equations are identically
satisfied. Finally, TG requires a pair of tetrad and spin connections where
the torsion $2$ form satisfies

$\mathbf{T}=\frac{1}{2}T^{A}_{\nu\mu}P_{A}dx^{\nu}\wedge
dx^{\mu}\,,\quad\text{with}\quad T^{A}_{\text{
}\nu\mu}=\partial_{\nu}e^{A}_{\mu}-\partial_{\mu}e^{A}_{\nu}+\omega^{A}_{\text{
}C\nu}e^{C}_{\mu}-\omega^{A}_{\text{ }C\mu}e^{C}_{\nu}\neq 0\,,$ (9)

where $P_{A}$ are the generators of the translational group. Hence, TG is a
gravitational theory where gravity is a manifestation of the torsion and not
of curvature of spacetime. The relation between the torsion and curvature
tensors with the torsion and curvature $2$-forms are given by

$R^{\rho}_{\ \gamma\nu\mu}=E^{\rho}_{A}e^{B}_{\gamma}R^{A}_{\ B\nu\mu},\quad
T^{\rho}_{\ \nu\mu}=E_{A}^{\ \rho}T^{A}_{\
\nu\mu}=-2\Gamma^{\rho}_{\left[\nu\mu\right]}\,,$ (10)

where the TG connection is given by

$\Gamma_{\nu\mu}^{\rho}\equiv
E_{A}^{\rho}\partial_{\mu}e^{A}_{\nu}+E_{A}^{\rho}\omega^{A}_{\text{
}B\mu}e^{B}_{\nu}=E_{A}^{\rho}\mathscr{D}_{\mu}e^{A}_{\nu}\,.$ (11)

Throughout this work, we will work on the Weitzenböck gauge, since for the
case of cosmological symmetry, the field equations associated to variations
with respect to the spin connection are satisfied [42].

Under this scheme, we can write the TG scalars111We can extend this to any
object, for instance, the Riemann tensor can be written w.r.t the Levi-Civita
connection using the contorsion tensor. in terms of the Levi-Civita connection
with the use of the contorsion tensor222All geometric objects with an over
circle are computed using the Levi-Civita connection.

$\Gamma^{\rho}_{\mu\nu}=\accentset{\circ}{\Gamma}^{\rho}_{\mu\nu}+K^{\rho}_{\mu\nu}\,.$
(12)

Particularly, if we compute the Ricci scalar, we arrive at [46]

$R=\accentset{\circ}{R}+T-B\,,$ (13)

with

$\displaystyle T$ $\displaystyle=T^{\alpha}_{\text{
}\sigma\rho}S_{\alpha}^{\text{
}\sigma\rho}=\frac{1}{4}T^{\mu\nu\lambda}T_{\mu\nu\lambda}+\frac{1}{2}T^{\mu\nu\lambda}T_{\nu\mu\lambda}-T^{\mu}T_{\mu}\,,$
(14) $\displaystyle S_{\alpha}^{\text{ }\sigma\rho}$
$\displaystyle=\frac{1}{4}(T_{\alpha}^{\text{
}\sigma\rho}+T^{\rho\sigma}_{\hskip 8.5359pt\alpha}-T^{\sigma\rho}_{\hskip
8.5359pt\alpha}-2T^{\lambda\sigma}_{\hskip
8.5359pt\lambda}\delta_{\alpha}^{\rho}+2T^{\lambda\rho}_{\hskip
8.5359pt\lambda}\delta_{\alpha}^{\sigma})=\frac{1}{2}\left(K_{\hskip
8.5359pt\alpha}^{\sigma\rho}+T^{\sigma}\delta_{\alpha}^{\rho}-T^{\rho}\delta_{\alpha}^{\sigma}\right)\,,$
(15) $\displaystyle B$
$\displaystyle=\frac{2}{e}\partial_{\mu}\left(eT^{\mu}\right)=2\nabla_{\mu}T^{\mu}\,.$
(16)

In TG, we arrive at the following equation

$\accentset{\circ}{R}=-T+B\,,$ (17)

which is arrived at due to the vanishing of curvature in this setting. Through
this expression and substituting the tetrad determinant $\sqrt{-g}=e$ into the
Einstein-Hilbert action [47], considering that $B$ is a boundary term as a
vanishing object at linear level, we obtain the action of the TEGR

$S_{\text{TEGR}}=-\frac{1}{2\kappa}\int eTd^{4}x+\int
e\mathcal{L}_{\text{m}}d^{4}x\,,$ (18)

with $\kappa=8\pi G$. The field equations of the TEGR action (18) recover the
same equations as in GR at the level of the equations of motion for particular
systems, leaving us with the same phenomenology of GR. In the same spirit as
$f(\accentset{\circ}{R})$ [48, 49, 50] extensions of TEGR have been gaining
interest in recent years [51, 52, 53], particularly, we look for extensions of
TEGR in the form of an arbitrary functional $f(T,B)$ into the action

$S_{f(T,B)}=-\frac{1}{2\kappa}\int ef(T,B)d^{4}x+\int
e\mathcal{L}_{\text{m}}d^{4}x\,,$ (19)

such that variations concerning the tetrad gives the following field equations

$\displaystyle W^{\lambda}_{\ \nu}\equiv$
$\displaystyle\delta^{\lambda}_{\nu}\accentset{\circ}{\square}f_{B}-\accentset{\circ}{\nabla}^{\lambda}\accentset{\circ}{\nabla}_{\nu}f_{B}+\frac{1}{2}f_{B}B\delta^{\lambda}_{\nu}+2\left[\partial_{\mu}f_{B}+\partial_{\mu}f_{T}\right]S_{\nu}^{\text{
}\mu\lambda}+\frac{2}{e}e^{A}_{\nu}f_{T}\partial_{\mu}(eS_{A}^{\text{
}\mu\lambda})$ (20) $\displaystyle-2f_{T}T^{\alpha}_{\text{
}\mu\nu}S_{\alpha}^{\text{
}\lambda\mu}-\frac{1}{2}f\delta^{\lambda}_{\nu}=\kappa\Theta^{\lambda}_{\
\nu}\,,$

where the energy-momentum tensor is defined as

$\Theta^{\lambda}_{\ \nu}=e^{A}_{\ \nu}\theta^{\lambda}_{\ A},\quad\theta_{\
A}^{\lambda}=\frac{1}{e}\frac{\delta(e\mathcal{L}_{\text{m}})}{\delta e^{A}_{\
\lambda}}\,.$ (21)

If the function is not dependent on the boundary term $f(T,B)=f(T)$, we
recover the field equations of $f(T)$ gravity, and if $f(T,B)=f(T)=-T$ we
recover the TEGR field equations [27]. The spin connection field equations
turn out to be written as the anti-symmetry of these equations, given by

$W_{[\mu\nu]}=0\,,$ (22)

where square brackets denote the antisymmetry operator.

In the case of cosmological symmetry, where the spacetime manifold is foliated
with maximally symmetric space-like slices
$\mathcal{M}=\mathbb{R}\times\Sigma$ [47, 54] showing isotropy and
homogeneity, we work with the FLRW metric given by

$ds^{2}=dt^{2}-a(t)^{2}\left(dx^{2}+dy^{2}+dz^{2}\right),$ (23)

in the spatially flat case $k=0$, with Cartesian coordinates The corresponding
tetrad field which exhibits the cosmological symmetry in the Weitzenböck gauge
is given by

$e^{A}_{\ \mu}=\text{diag}(1,a(t),a(t),a(t))\,,$ (24)

which turns out to be compatible with the Weitzenböck gauge. When working with
cosmological symmetry, the antisymmetric part of the field equations is
completely satisfied [42], leaving us the task to solve only the symmetric
part of the field equations. Using the background tetrad (24), we compute the
torsion and superpotential tensors to be [55]

$T^{i}_{\ 0j}=H\delta^{i}_{j}\,,\quad S_{i}^{\ 0j}=-H\delta^{i}_{j}\,,$ (25)

with the torsion scalar and boundary terms are

$T=-6H^{2}\,,\quad\text{and}\quad B=-6(3H^{2}+\dot{H})\,,$ (26)

and $H=\dot{a}/a$ is the Hubble parameter in cosmic time. The background field
equations for $f(T,B)$ are then given by

$\displaystyle
3H\dot{f}_{B}-3H^{2}(3f_{B}+2f_{T})-3f_{B}\dot{H}-\frac{1}{2}f(T,B)=\kappa\rho\,,$
(27)
$\displaystyle-(3H^{2}+\dot{H})(2f_{T}+3f_{B})-2H\dot{f}_{T}+\ddot{f}_{B}-\frac{1}{2}f(T,B)=-\kappa
P\,.$ (28)

The energy density $\rho$ and pressure $P$ correspond to those of a perfect
fluid.

From now on, we will focus only on these latter equations

for simplicity of the discussion, taking into account that $f(T,B)=f(T)$ and
$f(T,B)=-T$ recover the appropriate equations of each theory. We have been
working so far with cosmic time, however, in the case of linear perturbation
theory, we will switch to conformal time [56]

$ad\eta=dt\quad\rightarrow\quad\eta-\eta_{i}=\int_{t_{i}}^{t}\frac{dt^{\prime}}{a(t^{\prime})}\,,$
(29)

where the background tetrad becomes

$e^{A}_{\ \mu}=a(\eta)\text{diag}(1,1,1,1)\,,$ (30)

and the primes will indicate derivatives concerning conformal time.

##  3 Tensor perturbations in TG

In this section, we explore the tensor perturbations around a flat FLRW
background in order to probe the production of primordial GWs in TG. We break
the background cosmological symmetry by considering the first order deviation
of the background tetrad333The quantities with an overbar are background
quantities, which in this context exhibit the cosmological symmetry.

$e^{A}_{\ \mu}=\bar{e}^{A}_{\ \mu}+\delta e^{A}_{\ \mu}\,,$ (31)

such that $|\delta e^{A}_{\ \mu}|\ll 1$, such that the corresponding metric
perturbation is given by [57]

$\delta
g_{\mu\nu}=2\tau_{(\mu\nu)}\quad\text{where}\quad\tau_{\mu\nu}=\bar{e}^{A}_{\
\mu}\delta e^{B}_{\ \nu}\eta_{AB}\,.$ (32)

We consider the linear perturbation to be a symmetric transverse and traceless
tensor perturbation, in conformal time, given by

$\delta e^{A}_{\ \mu}=\frac{a(\eta)}{2}\begin{pmatrix}0&0\\\ 0&h_{ij}\\\
\end{pmatrix}\,,$ (33)

such that the linear perturbation of the torsion and superpotential tensors
are given by

$\delta T^{i}_{\ 0j}=\frac{1}{2}h_{ij}^{\prime},\quad\delta T^{i}_{\
jk}=\frac{1}{2}(\partial_{j}h_{ik}-\partial_{k}h_{ij}),\delta S_{0}^{\
0i}=0,\quad\delta S_{i}^{\ 0j}=\frac{1}{4a^{2}}h_{ij}^{\prime},\quad\delta
S_{i}^{\ jk}=-\frac{1}{4a^{2}}(\partial_{j}h_{ik}-\partial_{k}h_{ij})\,.$ (34)

The torsion scalar and boundary terms are given by

$\delta T=\delta B=0\,,$ (35)

from where we observe the Lagrangian functional form to have dependencies

$f(T,B)=f(\bar{T},\bar{B})\,,$ (36)

that is, they only depend on the background torsion scalar and boundary term.
Thus, it turns out that every time the function or its derivatives appear in
calculations, we should keep in mind that they are dependent only on
background quantities. Finally, the perturbed part on the left-hand side of
the field equations is given by

$\delta W^{i}_{\
j}=\frac{f_{T}}{2a^{2}}\delta^{ik}\left(h^{\prime\prime}_{kj}+[2+\nu]\mathcal{H}h^{\prime}_{kj}-\nabla^{2}h_{kj}\right)\,,$
(37)

from where, including tensor anisotropic effects, we arrive at the following
equation in Fourier space

$h^{\prime\prime}_{kj}+[2+\nu]\mathcal{H}h^{\prime}_{kj}+k^{2}h_{kj}=\frac{16\pi
Ga^{2}}{f_{T}}\pi_{ij}^{T}\,,$ (38)

where $\nu=\frac{1}{\mathcal{H}}\frac{f^{\prime}_{T}}{f_{T}}$ and
$\mathcal{H}=a^{\prime}/a$ [58, 59] Eq. (38) is the propagation equation for a
massless GW that propagates at the speed of light $c$, with Planck mass run
rate $\nu$, and which agrees with the literature form of this evolution
equation [60]. Notice that since we are working with a symmetric transverse
and traceless tensor perturbation, only two degrees of freedom are
propagating, namely, the so-called $+$ and $\times$ polarizations. We now
proceed to investigate the production of these GW coming from vacuum
fluctuations and possible tensor anisotropies. We will study the effects of
the extensions on the power spectrum and the energy density of GW. In Refs.
[55, 27], the condition $f_{T}<0$ is required for stability of the tensor
perturbation.

We begin by considering an inflationary scenario, an exponentially accelerated
expansion in the early universe driven by a scalar field that solves several
theoretical problems on the $\Lambda$CDM model [61]. By quantizing the GW
field as

$h_{ij}(\eta,\mathbf{x})=\sum_{\lambda=\pm
2}\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\left[h(\eta,k)e^{i\mathbf{k}\cdot\mathbf{x}}a(\mathbf{k},\lambda)e^{\lambda}_{ij}(\hat{k})+h^{*}(\eta,k)e^{-i\mathbf{k}\cdot\mathbf{x}}a^{\dagger}(\mathbf{k},\lambda){e^{\lambda}}^{*}_{ij}(\hat{k})\right]\;\,,$
(39)

with $e_{ij}^{\lambda}$ the polarization tensors, and we have defined
$h^{\lambda}(\eta,\mathbf{k})=h(\eta,k)a(\mathbf{k},\lambda)$ such that [62]

$h_{ij}=\sum_{\lambda}e_{ij}^{\lambda}h^{\lambda},\quad\sum_{i,j}e_{ij}^{\lambda}{e_{ij}^{\lambda^{\prime}}}^{*}=\delta^{\lambda\lambda^{\prime}}\,.$
(40)

The operators $a(\bm{k},\lambda),a^{\dagger}(\bm{k},\lambda)$ are quantum
operators that annihilate and create gravitons of momentun $k$ and helicity
$\lambda$, and satisfy the canonical commutation relations [56]

$\left[a(\mathbf{k},\lambda),a(\mathbf{k}^{\prime},\lambda^{\prime})\right]=0\;\,,\qquad\left[a(\mathbf{k},\lambda),a^{\dagger}(\mathbf{k}^{\prime},\lambda^{\prime})\right]=(2\pi)^{3}\delta^{(3)}(\mathbf{k}-\mathbf{k}^{\prime})\delta_{\lambda\lambda^{\prime}}\;\,.$
(41)

We will assume the state of the universe to be in a vacuum quantum state
$\lvert 0\rangle$ such that

$a(\bm{k},\lambda)\lvert 0\rangle=0\,,$ (42)

such that in the infinite past $\eta\to-\infty$, the modes $h(\eta,k)$
resemble the quantisation of a free field in Minkowski space [56][63]. This
state is called the Bunch-Davies vacuum state.

The correlator of the quantized tensor field is given by

$\langle
0|h_{ij}(\eta,\mathbf{x})h_{lm}(\eta,\mathbf{x}^{\prime})|0\rangle=\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}|h(\eta,k)|^{2}e^{i\mathbf{k}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\Pi_{ij,lm}(\hat{k})\;\,,$
(43)

where the sum over helicities is given by [64]

$\displaystyle\Pi_{ij,lm}(\hat{k})$ $\displaystyle\equiv\sum_{\lambda=\pm
2}e_{ij}(\hat{k},\lambda)e^{*}_{lm}(\hat{k},\lambda)\;$ (44)
$\displaystyle=\delta_{im}\delta_{jn}+\delta_{in}\delta_{jm}-\delta_{ij}\delta_{mn}+\delta_{ij}\hat{k}_{m}\hat{k}_{n}+\delta_{mn}\hat{k}_{i}\hat{k}_{j}-\delta_{im}\hat{k}_{j}\hat{k}_{n}-\delta_{jn}\hat{k}_{i}\hat{k}_{m}-\delta_{in}\hat{k}_{j}\hat{k}_{m}-\delta_{jm}\hat{k}_{i}\hat{k}_{n}$
$\displaystyle+\hat{k}_{i}\hat{k}_{j}\hat{k}_{m}\hat{k}_{n}\;\,.$ (45)

In superhorizon limits $k|\eta|\to 0$, particularly when the $k$-mode re-
enters the horizon, the quantum field can be seen as a classical field, hence,
we expect the correlator to become a classical average ensemble, from where
assuming gaussianity, we identify the power spectrum to be

$P_{h}(\eta,k)=|h(\eta,k)|^{2}\,,$ (46)

and the dimensionless power spectrum [65]

$\Delta_{h}^{2}(\eta,k)=\frac{k^{3}}{2\pi^{2}}|h(\eta,k)|^{2}\,.$ (47)

Now, we proceed to analyze the power spectrum for GW in TG coming from vacuum
fluctuations and also the contributions from tensor anisotropies.

###  3.1 Power spectrum from vacuum fluctuations

Vacuum fluctuations are generated where the inflaton field dominates. Since
the inflaton is a scalar field, no anisotropic stress is present during this
epoch [64], thus, the propagation equation for GW in TG is given by

$h_{ij}^{\prime\prime}+[2+\nu]\mathcal{H}h_{ij}^{\prime}+k^{2}h_{ij}=0\,.$
(48)

If we multiply by the polarization tensor and sum over helicities, we arrive
at

${h^{\lambda}}^{\prime\prime}+[2+\nu]\mathcal{H}{h^{\lambda}}^{\prime}+k^{2}h^{\lambda}=0\,,$
(49)

from which we observe that both polarization states satisfy the same
propagation equation in a vacuum, allowing us to work with $h^{\lambda}\equiv
h$ to simplify notation. We will now consider two cases of inflation, one with
a perfect de Sitter expansion and another with a quasi de Sitter expansion
modeled by the slow-roll parameter $\epsilon$.

From now on, we will be working with canonically normalized tensor
perturbations [56, 62]

$\hat{h}_{ij}=\frac{h_{ij}}{\sqrt{32\pi
G}}\equiv\frac{M_{\text{pl}}}{2}h_{ij}\,,$ (50)

for convenience

  * 1.

De Sitter background. A perfect de-sitter background is described by

$H=H_{\Lambda}=\text{const}\rightarrow\dot{H}=0\,,$ (51)

which has important implications for GWs in TG. Recalling that for tensor
perturbation the functional is only dependent on background quantities
$f(T,B)=f(\bar{T},\bar{B})$, we observe that, since $\bar{T}=-6H^{2}$ and
$\bar{B}=-6(3H^{2}+\dot{H})$, we have

$f^{\prime}_{T}=\frac{dt}{d\eta}\dot{f}_{T}=\frac{dt}{d\eta}\left(f_{TT}\dot{\bar{T}}+f_{TB}\dot{\bar{B}}\right)=0\,,$
(52)

from where the $\nu$ parameter vanishes, thus, the equation becomes

${\hat{h}}^{\prime\prime}+2\mathcal{H}{\hat{h}}^{\prime}+k^{2}{\hat{h}}=0\,,$
(53)

which is the equation for GW in GR [66]. Therefore, in a vacuum with a perfect
de Sitter background TG and GR produce the same gravitational waves. To solve
Eq. (53), we need to recast it as an equation with no damping term [56], which
is achieved by performing the change of variable

$\hat{h}=\frac{g}{a}\,,$ (54)

from where Eq. (53) becomes

$g^{\prime\prime}+\left(k^{2}-\frac{2}{\eta^{2}}\right)g=0\,,$ (55)

where we have used the fact that $a(\eta)=-\frac{1}{H_{\Lambda}\eta}$ in a
perfect de Sitter background. The solutions of Eq. (55) are given by

$g\subset\left\\{C(k)k\left(1-\frac{i}{k\eta}\right)e^{-ik\eta},\quad
C^{*}(k)k\left(1+\frac{i}{k\eta}\right)e^{ik\eta}\right\\}\,,$ (56)

where $C(k)$ is an arbitrary integration function of $k$. However, to satisfy
the initial conditions in the infinite past, we need to choose the first
solution

$g(\eta,k)=C(k)k\left(1-\frac{i}{k\eta}\right)e^{-ik\eta}\,,$ (57)

and evaluate it in the infinite past $k|\eta_{i}|\to\infty$, or
$k|\eta_{i}|\gg 1$, imposing the condition of the field to be a free field

$g(\eta_{i},k)=C(k)ke^{-ik\eta_{i}}=\frac{1}{\sqrt{2k}}e^{-ik\eta_{i}}\rightarrow
C(k)k=\frac{1}{\sqrt{2k}}\,,$ (58)

and then we can just return to the original $h$ field given by

$\hat{h}=\frac{1}{a}\frac{e^{-ik\eta}}{\sqrt{2k}}\left(1-\frac{i}{k\eta}\right)\,,$
(59)

from where the power spectrum, taking into account the canonical
normalization, is

$P_{h}(\eta,k)=\frac{16\pi
G}{k}H^{2}_{\Lambda}\eta^{2}\left(1+\frac{1}{k^{2}\eta^{2}}\right)\,,$ (60)

and the dimensionless power spectrum, considering that there are two
polarization states [66], is given by

$\Delta^{2}_{h}(\eta,k)\equiv\frac{d\langle 0\rvert\hat{h}^{2}_{ik}\lvert
0\rangle}{d\ln k}=64\pi
G\frac{k^{3}}{2\pi^{2}}|\hat{h}(\eta,k)|^{2}=\frac{2H^{2}_{\Lambda}}{\pi^{2}M_{\text{pl}}^{2}}[1+k^{2}\eta^{2}]\,.$
(61)

On the super horizon limit $k|\eta|\to 0$, we recover the scale-invariant
dimensionless power spectrum

$\Delta^{2}_{h}(k)=\frac{2}{\pi^{2}}\frac{H^{2}_{\Lambda}}{M_{\text{pl}}^{2}}\,,$
(62)

which is the same as GR as stated before [65]. As we have seen, in the case of
a perfect de Sitter background, there is no difference between the GW of GR
and TG, resulting in the same power spectrum. this prompts us to work in a
slightly different scenario, a quasi-de Sitter background.

  * 2.

Quasi de Sitter background. The quasi de Sitter background is defined by a
slowly varying derivative of the Hubble factor

$\dot{H}=-\epsilon H^{2}_{\Lambda}\,,\quad
a(\eta)=\frac{1}{H_{\Lambda}}\frac{1}{|\eta|^{1+\varepsilon}}\,,\quad\text{where}\quad|\eta|=-\eta\,,$
(63)

with $|\epsilon|\ll 1$ is the first slow roll parameter. Since the slow roll
parameter is small, all calculations will be performed in the first order in
$\epsilon$. Now, we need to compute the $\nu$ parameter for a general $f(T,B)$
functional. The idea is to Taylor expand at first order in $\epsilon$ the
$\nu$ parameter

$\nu(\epsilon)=\nu(\epsilon=0)+\epsilon\partial_{\epsilon}\nu|_{\epsilon=0}\,,$
(64)

such that $\nu(\epsilon=0)=0$ as we have argued previously since at
$\epsilon\to 0$, $f^{\prime}_{T}=0$. Using the same argument we obtain

$\partial_{\epsilon}\nu|_{\epsilon=0}=-\frac{\eta}{f_{T_{\Lambda}}}\partial_{\epsilon}f^{\prime}_{T}|_{\epsilon=0}\,,$
(65)

where introducing the notation $T_{\Lambda}=-6H_{\Lambda}^{2}$ and
$B_{\Lambda}=3T_{\Lambda}$, we arrive at

$\partial_{\epsilon}\nu|_{\epsilon=0}=2\left(\frac{f_{T_{\Lambda}T_{\Lambda}}|T_{\Lambda}|+f_{T_{\Lambda}B_{\Lambda}}|B_{\Lambda}|}{f_{T_{\Lambda}}}\right)\,,$
(66)

from where we obtain the $\nu$ parameter to be

$\nu=2\gamma\epsilon\,,\quad\text{with}\quad\gamma=\left(\frac{f_{T_{\Lambda}T_{\Lambda}}|T_{\Lambda}|+f_{T_{\Lambda}B_{\Lambda}}|B_{\Lambda}|}{f_{T_{\Lambda}}}\right)\,.$
(67)

Nevertheless, the stability condition $f_{T}<0$ needs to be fulfilled and it
depends on the particular form of the function. In the literature there have
been several models that can achieve these conditions and have been tested
against observations [67, 68, 69], we will expose some of these models next
444Recall that $\delta T=\delta B=0$..

    1. (a)

Power law models. There are models of $f(T,B)$ and $f(T)$ in the form of power
law as

$f(T)=-T+f_{0}(-T)^{m}\,,\quad f(T,B)=T+f_{0}(-T)^{m}+f_{1}(-B)^{m}\,.$ (68)

These models produce the same type of GW since both of them imply

$f_{T}=-1-mf_{0}(-T)^{m-1}\,,$ (69)

and then have associated the same $\nu$ parameter. From observations [70, 71]
we have $m,f_{0}>0$ such that $f_{T}<-1$, fulfilling the stability condition.

    2. (b)

Mixed power law model. In this case, the function is given by

$f(T,B)=-T+f_{0}(-T)^{m}(-B)^{n}\,,$ (70)

where the partial derivative w.r.t the torsion scalar is

$f_{T}=-1-mf_{0}(-T)^{m-1}(-B)^{n}\,,$ (71)

and according to observations [71] we have $f_{0},m>0$ implying $f_{T}<-1$,
also fulfilling the stability condition.

    3. (c)

Exponential model. In this case, we deal with a $f(T)$ model given by

$f(T)=-T+\beta T_{\Lambda}(1-e^{-qT/T_{\Lambda}})\,.$ (72)

In the limit, $q\to\infty$ the model is reduced to be GR [72]. Using this
model we have

$f_{T}=-1+\beta qe^{-qT/T_{\Lambda}}\,,$ (73)

where the stability condition $f_{T}<0$ requires $\beta q<0$ which is
fulfilled, since observations requires $q>0$ and $\beta<0$ [73].

The gamma parameter for these models is given by

$\displaystyle\Large{\gamma=\left\\{\begin{array}[]{ll}\frac{f_{0}m(1-m)|T_{\Lambda}|^{m-1}}{1+mf_{0}|T_{\Lambda}|^{m-1}}\,,&\quad\text{\normalsize{for
the power law models}},\\\
\frac{f_{0}m(1-m-n)|T_{\Lambda}|^{m-1}|B_{\Lambda}|^{n}}{1+mf_{0}|T_{\Lambda}|^{m-1}|B_{\Lambda}|^{n}}\,,&\quad\text{\normalsize{for
the mixed power law model,}}\\\ \frac{\beta q^{2}}{\beta
q-e^{q}}\,,&\quad\text{\normalsize{for the exponential
model.}}\end{array}\right.}$ (77)

With this in mind, we are now in a position to compute the power spectrum for
GW coming from vacuum fluctuations. The equation for GW for vacuum
fluctuations is

${\hat{h}}^{\prime\prime}+[2+\nu]\mathcal{H}{\hat{h}}^{\prime}+k^{2}{\hat{h}}=0\,,$
(78)

where we again aim to recast it as the equation of a harmonic oscillator
without the damping effect, however, the change of variable in (54) is not
suitable for that purpose. For an arbitrary $\nu$ parameter, the appropriate
change of variable is given by

$h(\eta,k)=f(\eta)g(\eta,k),\quad\text{with}\quad
f(\eta)=\frac{1}{a(\eta)}\text{Exp}\left(-\frac{1}{2}\int\mathcal{H}\nu
d\eta\right)\,.$ (79)

With this change of variable, the equation for GW in TG in the vacuum case
reduces to

$g^{\prime\prime}+g\left(k^{2}+\frac{f^{\prime\prime}}{f}+[2+\nu]\mathcal{H}\frac{f^{\prime}}{f}\right)=0\,,$
(80)

where in the case of a quasi-de Sitter background reads as

$g^{\prime\prime}+g\left(k^{2}-\frac{2+3(1+\gamma)\epsilon}{\eta^{2}}\right)=0\,.$
(81)

The solutions of equation (81) are given in terms of Hankel functions

$g(\eta,k)\subset\left\\{C(k)\sqrt{|\eta|}H_{\alpha}^{(1)}(k|\eta|),\quad
C^{*}(k)\sqrt{|\eta|}H_{\alpha}^{(2)}(k|\eta|)\right\\},\quad\alpha=\frac{3}{2}+\epsilon(1+\gamma)\,,$
(82)

such that the initial condition requires the choosing of the first solution,
wherein the infinite past implies the integration constant to be

$C(k)=\frac{\sqrt{\pi}}{2}e^{i\alpha\frac{\pi}{2}+i\frac{\pi}{4}}\,,$ (83)

where $\alpha$ is given in Eq. 82. The form of the solutions is completely
analogous to those of GR [56] but with a different $\alpha$ for the Hankel
functions, due to the presence of $\gamma$. If $\gamma=0$ we recover the same
solution as GR, which is expected. In the case of a quasi de Sitter, we have

$f(\eta)=\frac{1}{a(\eta)}\text{Exp}\left(-\frac{1}{2}\int\mathcal{H}\nu
d\eta\right)=a(\eta)^{-(1+\epsilon\gamma)}\,,$ (84)

from where the original $h$ field is

$\hat{h}=a^{-(1+\epsilon\gamma)}\frac{\sqrt{\pi}}{2}e^{i\frac{\alpha\pi}{2}+i\frac{\pi}{4}}\sqrt{|\eta|}H_{\alpha}^{(1)}(k|\eta|)\,.$
(85)

The power spectrum is then given by

$P_{h}(\eta,k)=32\pi
Ga^{-2(1+\epsilon\gamma)}\frac{\pi}{4}|\eta||H_{\alpha}^{(1)}(k|\eta|)|^{2}\,,$
(86)

and the dimensionless power spectrum is given by

$\Delta^{2}_{h}(\eta,k)=\frac{k^{3}|\eta|}{\pi
M_{\text{pl}}^{2}}a^{-2(1+\epsilon\gamma)}|H_{\alpha}^{(1)}(k|\eta|)|^{2}\,.$
(87)

In the superhorizon limit, the Hankel function has the following asymptotic
behavior

$H_{\alpha}^{(1)}(x)\underset{x\to
0}{\sim}-i\frac{\Gamma(\alpha)}{\pi}\left(\frac{x}{2}\right)^{-\alpha}\,,$
(88)

from where the dimensionless power spectrum, using the explicit form of the
scale factor, becomes

$\Delta^{2}_{h}(\eta,k)=\frac{k^{3}|\eta|\Gamma(\alpha)^{2}}{\pi^{3}M_{\text{pl}}^{2}}H_{\Lambda}^{2(1+\gamma\epsilon)}|\eta|^{2+2\epsilon(1+\gamma)}\left(\frac{k|\eta|}{2}\right)^{-3-2\epsilon(1+\gamma)}=\frac{\Gamma(\alpha)^{2}2^{3+2\epsilon(1+\gamma)}}{\pi^{3}M_{\text{pl}}^{2}}H_{\Lambda}^{2(1+\gamma\epsilon)}k^{-2\epsilon(1+\gamma)}\,,$
(89)

such that taking the $k$-dependence at first order in $\epsilon$ and zero-
order everywhere else, we obtain

$\Delta^{2}_{h}(\eta,k)=\frac{2H_{\Lambda}^{2}}{\pi^{2}M_{\text{pl}}^{2}}k^{-2\epsilon(1+\gamma)}\,.$
(90)

The power law dependence on the scale is called tensor spectral index, where
in the case of TG reads as

$n_{T}=-2\epsilon(1+\gamma)\,,$ (91)

which has an extra contribution $\gamma$ due to the presence of extensions of
TEGR. This is a key result since if the GW background is measured and a large
tensor spectral index is measured, it would suggest the need for extensions of
gravity, since in GR, $n_{T}=-2\epsilon$, and then must be small consequently.

###  3.2 Power spectrum including tensor anisotropies

For this section, we will follow a similar approach as in Ref. [74] but for a
general source of tensor anisotropic stress. The idea is to decompose the
tensor field into two modes, a short wavelength mode, and a slowly variant
long wavelength mode

$h(\eta,\bm{x})=h_{<}(\eta,\bm{x})+h_{>}(\eta,\bm{x})\,.$ (92)

The short wavelength part contains all the information of quantum vacuum
fluctuations and its effects are manifested through the parameters of the long
wavelength perturbation, in this case, through the power spectrum. The short
wavelength mode represents high momentum modes with wavelength much smaller
than the horizon, We then write the short wavelength mode as a quantized field
but with a filter or windows function

$h_{<}(\eta,\bm{x})=\int\frac{d^{3}\bm{k}}{\sqrt{(2\pi)^{3}}}W(\eta,k)\left[h(\eta,k)e^{i\mathbf{k}\cdot\mathbf{x}}a(\mathbf{k},\lambda)+h^{*}(\eta,k)e^{-i\mathbf{k}\cdot\mathbf{x}}a^{\dagger}(\mathbf{k},\lambda)\right]\,,$
(93)

where the filter function is given by

$W(\eta,k)=\theta(k-\varepsilon\mathcal{H})\,,$ (94)

and a suitable parameter $\varepsilon$ has been introduced such that it allows
us to guarantee that $h_{<}$ is a short wavelength contribution. Taking the GW
equations of TG in a vacuum and the physical space for a quasi-de Sitter
background, we have

$h^{\prime\prime}-\frac{2}{\eta}\left[1+\epsilon(1+\gamma)\right]h^{\prime}-\nabla^{2}h=0\,,$
(95)

such that, if we perform the split (92), we arrive at

$h^{\prime\prime}_{>}-\frac{2}{\eta}\left[1+\epsilon(1+\gamma)\right]h^{\prime}_{>}-\nabla^{2}h_{>}=\varrho_{Q},\quad\varrho_{Q}=-\left(h^{\prime\prime}_{<}-\frac{2}{\eta}\left[1+\epsilon(1+\gamma)\right]h^{\prime}_{<}-\nabla^{2}h_{<}\right)\,,$
(96)

with $\varrho_{Q}$ is called the quantum noise. In the case of TG in a quasi-
de Sitter background, the quantum noise has the following form

$\varrho_{Q}(\eta,\bm{x})=-\int\frac{d^{3}\bm{k}}{\sqrt{(2\pi)^{3}}}\left[a(\mathbf{k},\lambda)f^{\text{TG}}_{k}(\eta)e^{i\mathbf{k}\cdot\mathbf{x}}+a^{\dagger}(\mathbf{k},\lambda){f^{\text{TG}}_{k}}^{*}(\eta)e^{-i\mathbf{k}\cdot\mathbf{x}}\right]\,,$
(97)

with

$f^{\text{TG}}_{k}(\eta)=\left(W^{\prime\prime}-\frac{2}{\eta}\left[1+\epsilon(1+\gamma)\right]W^{\prime}\right)h+2W^{\prime}h^{\prime}\,.$
(98)

In the limit $\epsilon\to 0$ we recover the same quantum noise found in ref.
[62]. The correlator for the quantum noise is

$\langle
0\rvert\varrho_{Q}(\eta_{1},\bm{x}_{1})\varrho_{Q}(\eta_{2},\bm{x}_{2})\lvert
0\rangle=\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}e^{i\bm{k}\cdot(\bm{x}_{1}-\bm{x}_{2})}f_{k}(\eta_{1})f^{*}_{k}(\eta_{2})\,,$
(99)

or in Fourier space is

$\langle 0\rvert\varrho_{Q}(\eta_{1},\bm{k})\varrho_{Q}(\eta_{2},\bm{q})\lvert
0\rangle=\delta_{D}(\bm{k}-\bm{q})f_{k}(\eta_{1})f^{*}_{k}(\eta_{2})\,.$ (100)

However, we want to include explicitly the contributions of tensor anisotropic
stress, for that, we consider the propagation equation for GW in TG

$\hat{h}^{\prime\prime}_{kj}+[2+\nu]\mathcal{H}\hat{h}^{\prime}_{kj}+k^{2}\hat{h}_{kj}=\frac{16\pi
Ga^{2}}{f_{T}}\pi_{ij}^{T}\,,$ (101)

and then multiply both sides by $\epsilon_{ij}^{\lambda*}$, sum over the
spatial indices and use the orthonormality (40) of the polarization tensors to
arrive at

$\left(\partial^{2}_{\eta}+[2+\nu]\mathcal{H}\partial_{\eta}+k^{2}\right)\hat{h}^{\lambda}(\eta,\bm{k})=\frac{16\pi
Ga^{2}}{f_{T}}\sum_{ij}\epsilon_{ij}^{\lambda*}\pi_{ij}^{T}\equiv-\frac{1}{f_{T}}\varrho_{T}^{\lambda}(\eta,\bm{k})\,,$
(102)

with $\varrho_{T}^{\lambda}(\eta,\bm{k})$ is the tensor anisotropic noise. If
$f_{T}=-1$ we obtain the same equation as in ref. [62]. Now, we perform the
split into short and long wavelength modes (92), leaving us with the following
equation

$\left(\partial^{2}_{\eta}+[2+\nu]\mathcal{H}\partial_{\eta}+k^{2}\right)\hat{h}^{\lambda}_{>}(\eta,\bm{k})=\hat{\varrho}_{Q}^{\lambda}(\eta,\bm{k})-\frac{1}{f_{T}}\varrho_{T}^{\lambda}(\eta,\bm{k})\,,$
(103)

where the hat in $\hat{\varrho}_{Q}^{\lambda}(\eta,\bm{k})$ is to indicate
that it is computed using the canonical normalization of the tensor
perturbation. Now, we perform the same change of variable as before

$\hat{h}_{>}^{\lambda}(\eta,\bm{k})=f(\eta)g^{\lambda}_{>}(\eta,\bm{k})\quad\text{with}\quad
f(\eta)=a(\eta)^{-(1+\gamma\epsilon)}\,,$ (104)

such the equation for GW in a quasi de Sitter becomes

$\left(\partial^{2}_{\eta}+k^{2}-\frac{2+3\epsilon(1+\gamma)}{\eta^{2}}\right)g^{\lambda}_{>}(\eta,\bm{k})=a(\eta)^{1+\gamma\epsilon}\left(\hat{\varrho}_{Q}^{\lambda}(\eta,\bm{k})-\frac{1}{f_{T}}\varrho_{T}^{\lambda}(\eta,\bm{k})\right)\,.$
(105)

We will solve this equation using the method of Green’s function with retarded
boundary conditions, such that

$\left(\partial^{2}_{\eta}+k^{2}-\frac{2+3\epsilon(1+\gamma)}{\eta^{2}}\right)G_{R}(\eta,\eta_{i},k)=\delta_{D}(\eta-\eta_{i})\,,$
(106)

and $G_{R}(\eta,\eta_{i},k)=0$ if $\eta<\eta_{i}$. The Green’s function is
computed from the linear independent solution from the vacuum case in (82)
from where it reads as

$G_{\text{TG}}^{R}(\eta,\eta_{i},k)=\theta(\eta-\eta_{i})\frac{\pi}{2}\sqrt{|\eta|}\sqrt{|\eta_{i}|}\text{
Im}\left[H_{\alpha}^{(1)}(k|\eta_{i}|)H_{\alpha}^{(2)}(k|\eta|)\right],\quad\alpha=\frac{3}{2}+\epsilon(1+\gamma)\,.$
(107)

The solution of the equation is then

$g^{\lambda}_{>}(\eta,\bm{k})=\int_{-\infty}^{\eta}d\eta_{i}a(\eta_{i})^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{i},k)\hat{\varrho}_{Q}^{\lambda}(\eta_{i},\bm{k})-\int_{-\infty}^{\eta}d\eta_{i}\frac{a(\eta_{i})^{1+\gamma\epsilon}}{f_{T}(\eta_{i})}G_{\text{TG}}^{R}(\eta,\eta_{i},k)\varrho_{T}^{\lambda}(\eta_{i},\bm{k})\,.$
(108)

If we choose $\eta_{i}$ to be the time when the mode enters the horizon, we
can interpret the retarded Green’s function as the Green’s function with
boundary conditions such that when the mode has not re-entered the horizon
yet, the long wavelength mode is zero and we end up with only the short-
wavelength mode. Returning to the $h$ field, we obtain the general solution
for GW in TG in a quasi-de Sitter background that includes the effects of
tensor anisotropies

$\displaystyle\hat{h}^{\lambda}_{>}(\eta,\bm{k})=$
$\displaystyle\int_{-\infty}^{\eta}d\eta_{i}\left[\frac{a(\eta_{i})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{i},k)\hat{\varrho}_{Q}^{\lambda}(\eta_{i},\bm{k})-\int_{-\infty}^{\eta}d\eta_{i}\frac{1}{f_{T}(\eta_{i})}\left[\frac{a(\eta_{i})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{i},k)\varrho_{T}^{\lambda}(\eta_{i},\bm{k})\,.$
(109)

This result is of high importance since it opens the door to studying widely
different possible scenarios of generations of GW in the context of TG. if
$f_{T}=-1$, then $\gamma=0$, and then we have also the general solution in GR
for a quasi-de Sitter background. In the limit $\epsilon\to 0$, the Green’s
function satisfy

$\frac{a(\eta_{i})}{a(\eta)}G^{R}(\eta,\eta_{i},k)=\frac{\theta(\eta-\eta_{i})}{\eta_{i}^{2}k^{3}}\text{Im}\left[e^{ik(\eta-\eta_{i})}(1-ik\eta)(1+ik\eta_{i})\right]\,,$
(110)

which is the same Green’s function found in [62]. We are now in a position to
see the form of the power spectrum using this result. Taking the equal time
correlator of the general solution, we have

$\displaystyle\left\langle\hat{h}^{\lambda}_{>}(\eta,\bm{k})\hat{h}^{\lambda^{\prime}}_{>}(\eta,\bm{q})\right\rangle=$
(111)
$\displaystyle\int_{-\infty}^{\eta}d\eta_{1}\left[\frac{a(\eta_{1})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{1},k)\int_{-\infty}^{\eta}d\eta_{2}\left[\frac{a(\eta_{2})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{2},k)\langle
0\rvert\hat{\varrho}_{Q}^{\lambda}(\eta_{1},\bm{k})\hat{\varrho}_{Q}^{\lambda^{\prime}}(\eta_{2},\bm{q})\lvert
0\rangle$
$\displaystyle+\int_{-\infty}^{\eta}d\eta_{1}\frac{1}{f_{T}(\eta_{1})}\left[\frac{a(\eta_{1})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{1},k)\int_{-\infty}^{\eta}d\eta_{2}\frac{1}{f_{T}(\eta_{2})}\left[\frac{a(\eta_{2})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{2},k)$
$\displaystyle\langle
0\rvert\varrho_{T}^{\lambda}(\eta_{1},\bm{k})\varrho_{T}^{\lambda^{\prime}}(\eta_{2},\bm{q})\lvert
0\rangle\,,$

where recalling the correlator of the quantum noise

$\langle
0\rvert\hat{\varrho}_{Q}^{\lambda}(\eta_{1},\bm{k})\hat{\varrho}_{Q}^{\lambda^{\prime}}(\eta_{2},\bm{q})\lvert
0\rangle=32\pi
G\delta_{D}(\bm{k}-\bm{q})\delta^{\lambda\lambda^{\prime}}f_{k}^{\text{TG}}(\eta_{1}){f_{k}^{\text{TG}}}^{*}(\eta_{2})\,,$
(112)

and taking into account the two polarization states implies that

$\displaystyle\left\langle\hat{h}_{>}(\eta,\bm{k})\hat{h}_{>}(\eta,\bm{q})\right\rangle=64\pi
G\delta_{D}(\bm{k}-\bm{q})\left\lvert\int_{-\infty}^{\eta}d\eta_{i}\left[\frac{a(\eta_{i})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{i},k)f_{k}^{\text{TG}}(\eta_{i})\right\rvert^{2}$
(113)
$\displaystyle+\int_{-\infty}^{\eta}d\eta_{1}\frac{1}{f_{T}(\eta_{1})}\left[\frac{a(\eta_{1})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{1},k)\int_{-\infty}^{\eta}d\eta_{2}\frac{1}{f_{T}(\eta_{2})}\left[\frac{a(\eta_{2})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{2},k)\sum_{\lambda}\langle
0\rvert\varrho_{T}^{\lambda}(\eta_{1},\bm{k})\varrho_{T}^{\lambda}(\eta_{2},\bm{q})\lvert
0\rangle\,.$

By construction, before the mode enter the horizon $-\eta<-\eta_{*}$, with
$\eta_{*}$ the time of horizon crossing, the term

$P_{\Lambda}(\eta,k)=64\pi
G\left\lvert\int_{-\infty}^{\eta}d\eta_{i}\left[\frac{a(\eta_{i})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{i},k)f_{k}^{\text{TG}}(\eta_{i})\right\rvert^{2}=\frac{2\pi|\eta|}{M_{\text{pl}}^{2}}a^{-2(1+\epsilon\gamma)}|H_{\alpha}^{(1)}(k|\eta|)|^{2}\,,$
(114)

which is the contribution from the vacuum computed previously in (86)
considering both polarization states. Analogously, we define the contribution
to the power spectrum from tensor anisotropic stress as

$\displaystyle\delta_{D}(\bm{k}-\bm{q})P_{T}(\eta,k)=$
$\displaystyle\int_{-\infty}^{\eta}d\eta_{1}\frac{1}{f_{T}(\eta_{1})}\left[\frac{a(\eta_{1})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{1},k)\int_{-\infty}^{\eta}d\eta_{2}\frac{1}{f_{T}(\eta_{2})}\left[\frac{a(\eta_{2})}{a(\eta)}\right]^{1+\gamma\epsilon}G_{\text{TG}}^{R}(\eta,\eta_{2},k)$
(115) $\displaystyle\sum_{\lambda}\langle
0\rvert\varrho_{T}^{\lambda}(\eta_{1},\bm{k})\varrho_{T}^{\lambda}(\eta_{2},\bm{q})\lvert
0\rangle\,.$

There exist possible sources for anisotropic stress studied in the context of
GR, for instance, local thermal fluctuations of the primordial plasma [62],
first-order cosmological phase transitions [75, 76, 77], and magnetic fields
[78], to name a few. Whether or not those physical mechanisms sourcing GW are
present in TG is a matter of future research, however, at this point, it is
possible to discuss the implication of the extended models on those peaks on
the power spectrum, as we shall see later on.

##  4 Energy density of gravitational waves in teleparallel gravity

The energy density spectrum is defined as the gravitational-wave energy
density per logarithmic wave number interval is given by [66]

$\Omega_{\text{GW}}(\eta,k)=\frac{1}{\rho_{\text{crit}}(\eta)}\frac{\langle
0\rvert\rho_{\text{GW}}(\eta)\lvert 0\rangle}{d\ln k}\,,$ (116)

where the critical density is given by

$\rho_{\text{crit}}(\eta)=\frac{3H^{2}(\eta)}{8\pi G}\,.$ (117)

To compute $\rho_{\text{GW}}$ we will compute the energy-momentum tensor
associated with the action for GW in TG, which is given by [27]

$S=\int d^{3}\bm{x}d\eta(-a^{2}f_{T})\left[h^{\prime}_{ij}h^{\prime
ij}-\partial_{k}h_{ij}\partial^{k}h^{ij}\right]\,,$ (118)

or in terms of the background tetrad

$S=\int d^{3}\bm{x}d\eta\bar{e}(-f_{T})\left[\bar{E}^{\mu}_{\
A}\bar{E}^{\nu}_{\
B}\eta^{AB}\partial_{\mu}h_{ij}\partial_{\nu}h_{ij}\right]\,,$ (119)

where it is possible to identify the Lagrangian density to be

$\mathcal{L}=(-f_{T})\left[\bar{E}^{\mu}_{\ A}\bar{E}^{\nu}_{\
B}\eta^{AB}\partial_{\mu}h_{ij}\partial_{\nu}h_{ij}\right]\,.$ (120)

With this Lagrangian density, it is possible to compute the associated energy-
momentum tensor given by eq. (21), taking into account the signature
$(+,-,-,-)$, whose result is

$\mathcal{T}^{\lambda}_{\ \nu}=2(-f_{T})\bar{E}^{\lambda}_{\
A}\bar{E}^{\mu}_{\
B}\eta^{AB}\partial_{\mu}h_{ij}\partial_{\nu}h_{ij}+\left[\bar{E}^{\mu}_{\
C}\bar{E}^{\alpha}_{\
B}\eta^{CB}\partial_{\mu}h_{ij}\partial_{\alpha}h_{ij}\right]\bar{e}^{A}_{\
\nu}\frac{1}{\bar{e}}W^{\lambda}_{\ A}[f_{T},\bar{e}]\,,$ (121)

where we have defined

$\displaystyle W^{\lambda}_{\
A}[f_{T},\bar{e}]\equiv\frac{\delta(\bar{e}f_{T}(T,B))}{\delta\bar{e}^{A}_{\
\lambda}}=$
$\displaystyle-\bar{e}f_{TB}\bar{B}\bar{E}^{\lambda}_{A}-2\bar{e}\bar{E}^{\lambda}_{A}\accentset{\circ}{\square}f_{TB}+2\bar{e}\bar{E}^{\nu}_{A}\accentset{\circ}{\bar{\nabla}}^{\lambda}\accentset{\circ}{\bar{\nabla}}_{\nu}f_{TB}-4\bar{e}(\partial_{\mu}f_{TB})\bar{S}_{A}^{\text{
}\mu\lambda}$ (122)
$\displaystyle-4\partial_{\mu}(f_{TT})\bar{e}\bar{S}_{A}^{\text{
}\mu\lambda}-4f_{TT}\partial_{\mu}\left(\bar{e}\bar{S}_{A}^{\mu\lambda}\right)+4\bar{e}f_{TT}\bar{T}^{\alpha}_{\text{
}\mu A}\bar{S}_{\alpha}^{\text{
}\lambda\mu}+\bar{e}f_{T}\bar{E}^{\lambda}_{A}\,,$

to keep the same notation as the case of the field equations and we have
neglected total derivatives. Thus, we identify

$\rho_{\text{GW}}^{\text{TG}}=\Theta^{0}_{\
0}=\frac{2}{a^{2}(\eta)}\left[-(h^{\prime}_{ij})^{2}(f_{T}+\kappa\bar{\rho}_{f_{T}})+\kappa\bar{\rho}_{f_{T}}(\nabla
h_{ij})^{2}\right]\,,$ (123)

with

$a^{2}(\eta)\kappa\bar{\rho}_{f_{T}}\equiv
3\mathcal{H}f^{\prime}_{TB}-6\mathcal{H}^{2}(f_{TB}+f_{TT})-3f_{TB}\mathcal{H}^{\prime}-\frac{a^{2}(\eta)}{2}f_{T}\,.$
(124)

Notice that $f_{T}=-1$ implies $\kappa\bar{\rho}_{f_{T}}=\frac{1}{2}$ and we
recover the GR limit [66]. Using the condition at horizon crossing [66]

$|h^{\prime}(\eta,k)|^{2}=k^{2}|h(\eta,k)|^{2}\,,$ (125)

and going to Fourier space, we obtain

$\Omega(\eta,k)=\frac{1}{12}\frac{k^{2}\Delta_{h}^{2}(\eta,k)}{a^{2}(\eta)H^{2}(\eta)}\left(-f_{T}\right)\,,$
(126)

where the GR result is trivially recovered by $f_{T}=-1$.

##  5 Discussions

There are several implications of the main results of this work. On the one
hand, when working in a vacuum for a perfect de Sitter background, both GW in
GR and TG satisfy the same propagation equation (53), and therefore, have the
same analytical solution. This implies that the associated scale-invariant
power spectrum is the same for both theories. Thus, if the GW background is
measured with zero or nearly zero tensor spectral index (scale-invariant power
spectrum), both theories will survive, and there are no direct implications of
the plausibility of TG as the underlying gravitational theory.

On the other hand, in the case when the background is promoted to a quasi de
Sitter background, driven by the first slow roll parameter $\epsilon$, the
situation is quite different. As we have seen, the $\gamma$ parameter encodes
all the information about the presence of extensions of gravity. In the case
of a quasi de Sitter background, the tensor spectral index is
$n_{T}=-2\epsilon(1+\gamma)$. Hence, if a future GW detector is capable of
directly measuring the tensor spectral index and observes a value way
different from zero, this would indicate a need of extensions of gravity since
TG allows for high values of the tensor scalar index. This is an important
result and relies on the efficacy of future GW detectors for measures of the
GW background, since the only way we have to infer the tensor spectral index
is through the tensor-to-scalar ratio, and this ratio is highly dependent on
the gravitational theory.

Therefore, a direct measurement of zero tensor spectral index does not provide
enough evidence to discard this expression of TG, but a measurement of a
tensor spectral index way higher than zero is a direct implication of the need
for an extended gravitational model, where we have shown that such result is
perfectly explained by TG, providing a significant prediction and not only a
validation.

Furthermore, if we take the limit $\epsilon\to 0$, the $f(T,B)$ functional
will not depend on time $f^{\prime}_{T}=0$, and Green’s function is given in
Eq. (107) will recover that of GR

$\frac{a(\eta_{i})}{a(\eta)}G^{R}(\eta,\eta_{i},k)=\frac{\theta(\eta-\eta_{i})}{\eta_{i}^{2}k^{3}}\text{Im}\left[e^{ik(\eta-\eta_{i})}(1-ik\eta)(1+ik\eta_{i})\right]\equiv
G_{R}(\eta,\eta_{i},k)\,,$ (127)

using the same notation as in [62]. This implies that the tensor contribution
to the power spectrum can be written as

$\displaystyle\delta_{D}(\bm{k}-\bm{q})P_{T}(\eta,k)=$
$\displaystyle\frac{1}{f_{T}^{2}}\int_{-\infty}^{\eta}d\eta_{i}G_{R}^{2}(\eta,\eta_{i},k)\sum_{\lambda}\langle
0\rvert\varrho_{T}^{\lambda}(\eta_{i},\bm{k})\varrho_{T}^{\lambda}(\eta_{i},\bm{q})\lvert
0\rangle\equiv\frac{1}{f_{T}^{2}}\delta_{D}(\bm{k}-\bm{q})P_{T}^{GR}(\eta,k)\,,$
(128)

with

$\delta_{D}(\bm{k}-\bm{q})P_{T}^{GR}(\eta,k)=\int_{-\infty}^{\eta}d\eta_{i}G_{R}^{2}(\eta,\eta_{i},k)\sum_{\lambda}\langle
0\rvert\varrho_{T}^{\lambda}(\eta_{i},\bm{k})\varrho_{T}^{\lambda}(\eta_{i},\bm{q})\lvert
0\rangle\,,$ (129)

the contribution from tensor anisotropies on the power spectrum. Therefore, we
will have either an amplification on the peaks of the power spectrum, if
$-1<f_{T}<0$, or a reduction, if $f_{T}<-1$, coming from tensor anisotropies.
For particular models studied in the literature, like the ones we presented
previously in Sec. 3.1, the functional satisfies $f_{T}=-1+F_{T}$ with
$F_{T}<0$, hence $f_{T}<-1$. This would imply that the peaks from tensor
anisotropies are squeezed in those particular models, giving us another
important observational prediction. If the peaks on the power spectrum, from
tensor anisotropies, of primordial GW are decreased compared to what should be
expected from GR, it is a hint of the need for TG as a possible gravitational
theory, where the amount of squeezing will depend highly on the particular
choice of the $(T,B)$ functional. Finally, on the same $\epsilon\to 0$ limit
and $f_{T}=-1+F_{T}$ with $F_{T}<0$, it is possible to observe that

$\Omega(\eta,k)=\frac{1}{12}\frac{k^{2}\Delta_{h}^{2}(\eta,k)}{a^{2}(\eta)H^{2}(\eta)}-F_{T}\frac{1}{12}\frac{k^{2}\Delta_{h}^{2}(\eta,k)}{a^{2}(\eta)H^{2}(\eta)}\,,$
(130)

and identifying the contribution from GR and the contribution from extensions
of gravity

$\Omega_{\text{GR}}\equiv\frac{1}{12}\frac{k^{2}\Delta_{h}^{2}(\eta,k)}{a^{2}(\eta)H^{2}(\eta)}\,,\quad\Omega_{F_{T}}=-F_{T}\frac{1}{12}\frac{k^{2}\Delta_{h}^{2}(\eta,k)}{a^{2}(\eta)H^{2}(\eta)}\,,$
(131)

we have

$\Omega(\eta,k)=\Omega_{\text{GR}}+\Omega_{F_{T}}>\Omega_{GR}\,,$ (132)

this implies that for the particular cosmological viable models of $f(T,B)$
presented in subsection 3.1, the amount of energy density in TG is greater
than that of GR, which is also an important prediction for TG since if the
energy spectrum observed is greater than the expected from GR, the need for an
extension of gravity will be manifest with $f(T,B)$ or $f(T)$ being a
candidate to explain that extra energy density.

The bigger energy density of GW in TG, along with the results of reduction in
the peaks of the power spectrum, implies that the vacuum fluctuations have a
bigger contribution to the energy density of GW in TG compared with its
contribution in GR, and the contributions from tensor anisotropies have less
contribution.

##  6 Conclusions

TG is a gravitational theory built as a gauge theory of translations and
invariance under local Lorentz transformations. In its equivalent extension
limit to TEGR, TG is completely equivalent to GR, hence showing the same
advantages and disadvantages in explaining the cosmological phenomena. The
need for extended models is an important subject of research in TG to try to
solve the GR problems which are manifest in the TEGR limit. These models have
been explored widely in the context of the late universe to address the
cosmological constant issue, cosmography, and other problems, and some
contributions in the context of the early universe have been performed,
including studies of primordial black holes [37] and inflation [69, 38, 40].
As a further contribution, in this paper, we study the production of
primordial GW in the context of TG. We showed that in a perfect de Sitter
background the propagation equation in TG is the same as in the GR case, then,
the same signal as GR of the power spectrum is recovered in TG since it is
given from the squared modulus of the solution of the same equation for GW. In
the quasi de Sitter limit we obtained that TG allows for a high value of the
tensor spectral index, providing a crucial test of TG if a direct measurement
of the power spectrum of GW from vacuum fluctuations is done and a tensor
spectral index obtained is considerably different from zero, this would imply
the need of TG. It is interesting to notice that our results could be tested
via future data from pulsar timing arrays, e.g. PTAs. In particular, such
experiments have constrained some effects related to stochastic gravitational
wave backgrounds as sources of cosmic inflation [79]. A successful program
through these experiments could determine if we will observe a blue tilt on
the tensor power spectrum. Finally, we obtained that although the vacuum
contribution in a perfect de Sitter background is the same for both GR and TG.
In both cases, the peaks of the contributions from tensor anisotropies are
different, particularly, we explored these differences for some particular
cosmological viable models, showing that the peaks are squeezed in TG compared
to GR.

Other works on the topic of gravitational waves induced from primordial black
holes include Refs. [80, 81] which tackle different aspects of the topic. In
the present work, we focus on the generation of primordial gravitational
waves, in the context of linear perturbation theory, coming from vacuum
fluctuations and tensor anisotropic stress which is a first order perturbation
quantity of the energy momentum tensor and do not address the generation of
gravitational waves induced from second order scalar perturbation. Secondly,
and the most important one, primordial black holes are formed in the radiation
dominated era, and in this work, we go far earlier and address the generation
during the inflaton dominated era where the universe stays in a quantum state.

Finally, we showed that we expect a greater amount of energy density in TG
compared to GR. This result implies that along with a reduction in the peaks
of the power spectrum, the vacuum fluctuations have a greater contribution to
the energy density in TG and the tensor anisotropies have a less contribution
when compared to GR. These results are the basis for the study of primordial
GW in TG, but to perform a complete numerical analysis that involves
cosmological constraints, we need to compute the transfer functions in TG and
study the characteristic signals that should be expected within the redshift
range of the GW experiment, and then use the data to test both theories.

Within the same latter scheme, we included tensor anisotropies sourcing GW,
and the direct observational implications on the peaks of the power spectrum
were computed. This is an important step in studying TG in the context of the
early universe. In the case of only vacuum contributions, the result for the
power spectrum needs to be taken in the superhorizon limit and evaluated at
the time of horizon crossing when the $k$ mode re-enters the horizon. In the
case of including tensor anisotropies, it should be evaluated at a time
$\tilde{\eta}$ shortly after the end of inflation, and it will depend on the
phenomenology of the tensor source. Nevertheless, these expressions need to be
studied at cosmic times (redshifts) achievable by future experiments, and that
will require the study of transfer functions in TG which will allow us to
evaluate these results at the required scales. However such a task lies beyond
the scope of this paper, but in the light of future GW detectors, like LISA or
the Einstein Telescope, this should bring the attention of the TG community to
explore these topics with the general solutions provided here.

## Acknowledgmenents

GRF is supported by CONACyT National Grant. CE-R acknowledges the Royal
Astronomical Society as FRAS 10147 and is supported by PAPIIT UNAM Project
TA100122. This article is based upon work from COST Action CA21136 Addressing
observational tensions in cosmology with systematics and fundamental physics
(CosmoVerse) supported by COST (European Cooperation in Science and
Technology).

##  Appendix A Tensor perturbations

The idea of tensor perturbation consists in splitting the tetrad field into
two arts, a background tetrad, referring to the background manifold exhibiting
cosmological symmetry $\bar{e}^{A}_{\ \mu}$, and a small deviation from this
background $\delta e^{A}_{\ \mu}$ such that

$e^{A}_{\ \mu}=\bar{e}^{A}_{\ \mu}+\delta e^{A}_{\ \mu}.$ (133)

We can identify the geometrical part of field equations (20) as

$\displaystyle W^{\lambda}_{\nu}=\delta^{\lambda}_{\nu}A-B^{\lambda}_{\
\nu}+\frac{1}{2}\delta^{\lambda}_{\nu}C+2D^{\lambda}_{\ \nu}+2F^{\lambda}_{\
\nu}-2H^{\lambda}_{\ \nu}-\frac{1}{2}f\delta^{\lambda}_{\nu},$ (134)

with

$\displaystyle A=\accentset{\circ}{\square}f_{B},\quad B^{\lambda}_{\
\nu}=\accentset{\circ}{\nabla}^{\lambda}\accentset{\circ}{\nabla}_{\nu}f_{B},\quad
C=f_{B}B,\quad D^{\lambda}_{\
\nu}=\left[\partial_{\mu}f_{B}+\partial_{\mu}f_{T}\right]S_{\nu}^{\text{
}\mu\lambda},$ (135) $\displaystyle F^{\lambda}_{\ \nu}=\frac{1}{e}e^{A}_{\
\nu}f_{T}\partial_{\mu}(eS_{a}^{\ \mu\lambda}),\quad H^{\lambda}_{\
\nu}=f_{T}T^{\alpha}_{\ \mu\nu}S_{\alpha}^{\ \lambda\mu}.$ (136)

The linear perturbation of the geometrical part of field equations (134) due
to the linear perturbation (133) is

$\displaystyle\delta W^{\lambda}_{\nu}=\delta^{\lambda}_{\nu}\delta A-\delta
B^{\lambda}_{\ \nu}+\frac{1}{2}\delta^{\lambda}_{\nu}\delta C+2\delta
D^{\lambda}_{\ \nu}+2\delta F^{\lambda}_{\ \nu}-2\delta H^{\lambda}_{\
\nu}-\frac{1}{2}\delta f\delta^{\lambda}_{\nu},$ (137)

with

$\displaystyle\delta A$ $\displaystyle=\delta
g^{\mu\beta}\left(\partial_{\mu}\partial_{\beta}\bar{f}_{\bar{B}}-\accentset{\circ}{\bar{\Gamma}}^{\rho}_{\beta\mu}\partial_{\rho}\bar{f}_{\bar{B}}\right)+\bar{g}^{\mu\beta}\left(\partial_{\mu}\partial_{\beta}\delta
f_{B}-\delta\accentset{\circ}{\Gamma}^{\rho}_{\beta\mu}\partial_{\rho}\bar{f}_{\bar{B}}-\accentset{\circ}{\bar{\Gamma}}^{\rho}_{\beta\mu}\partial_{\rho}\delta
f_{B}\right),$ (138) $\displaystyle\delta B^{\lambda}_{\ \nu}$
$\displaystyle=\delta
g^{\lambda\beta}\accentset{\circ}{\bar{\nabla}}_{\beta}\accentset{\circ}{\bar{\nabla}}_{\nu}\bar{f}_{\bar{B}}+\bar{g}^{\lambda\beta}\left(\accentset{\circ}{\bar{\nabla}}_{\beta}\accentset{\circ}{\bar{\nabla}}_{\nu}\delta
f_{B}-\delta\accentset{\circ}{\Gamma}^{\rho}_{\nu\beta}\partial{\rho}\bar{f}_{\bar{B}}\right),$
(139) $\displaystyle\delta C$ $\displaystyle=\bar{f}_{\bar{B}}\delta B+\delta
f_{B}\bar{B},$ (140) $\displaystyle\delta D^{\lambda}_{\ \nu}$
$\displaystyle=\partial_{\mu}\left(\delta f_{B}+\delta
f_{T}\right)\bar{S}_{\nu}^{\
\mu\lambda}+\partial_{\mu}\left(\bar{f}_{\bar{B}}+\bar{f}_{\bar{T}}\right)\delta
S_{\nu}^{\ \mu\lambda},$ (141) $\displaystyle\delta e$
$\displaystyle=\bar{e}\bar{E}^{\mu}_{\ B}\delta e^{B}_{\ \mu},\quad\delta
F^{\lambda}_{\ \nu}=\frac{1}{\bar{e}}\bar{e}^{A}_{\
\nu}\bar{f}_{\bar{T}}\partial_{\mu}\left(\bar{e}\delta S_{A}^{\
\mu\lambda}+\delta e\bar{S}_{A}^{\
\mu\lambda}\right)+\left[\frac{1}{\bar{e}}\bar{e}^{A}_{\ \nu}\delta
f_{T}+\left(\frac{1}{\bar{e}}\delta e^{A}_{\ \nu}-\frac{\delta
e}{\bar{e}^{2}}\bar{e}^{A}_{\
\nu}\right)\bar{f}_{\bar{T}}\right]\partial_{\mu}\left(\bar{e}\bar{S}_{A}^{\
\mu\lambda}\right),$ (142) $\displaystyle\delta H^{\lambda}_{\ \nu}$
$\displaystyle=\bar{S}_{\alpha}^{\ \lambda\mu}\left[\bar{f}_{\bar{T}}\delta
T^{\alpha}_{\mu\nu}+\delta
f_{T}\bar{T}^{\alpha}_{\mu\nu}\right]+\bar{f}_{\bar{T}}\bar{T}^{\alpha}_{\mu\nu}\delta
S_{\alpha}^{\ \lambda\mu},$ (143)

where

$\displaystyle\delta\accentset{\circ}{\Gamma}^{\rho}{}_{\mu\nu}$
$\displaystyle=\frac{1}{2}\bar{g}^{\rho\sigma}\left(\accentset{\circ}{\bar{\nabla}}_{\mu}\delta
g_{\sigma\nu}+\accentset{\circ}{\bar{\nabla}}_{\nu}\delta
g_{\mu\sigma}-\accentset{\circ}{\bar{\nabla}}_{\sigma}\delta
g_{\mu\nu}\right)=\bar{g}^{\rho\sigma}\left(\accentset{\circ}{\bar{\nabla}}_{\mu}\tau_{(\sigma\nu)}+\accentset{\circ}{\bar{\nabla}}_{\nu}\tau_{(\mu\sigma)}-\accentset{\circ}{\bar{\nabla}}_{\sigma}\tau_{(\mu\nu)}\right)$
(144)
$\displaystyle=\bar{g}^{\rho\sigma}\left(\tau_{(\sigma\mu),\nu}+\tau_{(\sigma\nu),\mu}-\tau_{(\mu\nu),\sigma}-2\tau_{(\sigma\alpha)}\accentset{\circ}{\bar{\Gamma}}^{\alpha}_{\
\mu\nu}\right),$ (145)

$\tau_{\mu\nu}$ given by (32), the bar indicates that the corresponding
geometrical is object is computed with the background tetrad, and the
perturbations for the functional are

$\displaystyle f(T,B)=f(\bar{T}+\delta T,\bar{B}+\delta
B)=\bar{f}+\bar{f}_{\bar{T}}\delta T+\bar{f}_{\bar{B}}\delta
B\equiv\bar{f}+\delta f,$ (146)

where we have defined $\bar{f}=f(\bar{T},\bar{B})$ and $\delta
f=\bar{f}_{\bar{T}}\delta T+\bar{f}_{\bar{B}}\delta B$, and for $f_{T}$ and
$f_{B}$ we have

$\displaystyle f_{T}(T,B)=\bar{f}_{\bar{T}}+\delta f_{T},\quad
f_{B}(T,B)=\bar{f}_{\bar{B}}+\delta f_{B}$ (147)

where

$\displaystyle\delta f_{T}=\bar{f}_{\bar{T}\bar{T}}\delta
T+\bar{f}_{\bar{T}\bar{B}}\delta B\quad\text{and}\quad\delta
f_{B}=\bar{f}_{\bar{B}\bar{T}}\delta T+\bar{f}_{\bar{B}\bar{B}}\delta B.$
(148)

For the tensor perturbation given by eq. (33), we obtain

$\displaystyle\delta T^{i}_{\ 0j}=\frac{1}{2}h_{ij}^{\prime},\quad\delta
T^{i}_{\ jk}=\frac{1}{2}(\partial_{j}h_{ik}-\partial_{k}h_{ij}),\delta
S_{0}^{\ 0i}=0,\quad\delta S_{i}^{\ 0j}=\frac{1}{4a^{2}}h_{ij}^{\prime},$
(149) $\displaystyle\delta S_{i}^{\
jk}=-\frac{1}{4a^{2}}(\partial_{j}h_{ik}-\partial_{k}h_{ij}),\quad\delta
T=\delta B=0,$ (150)

and then using the aforementioned perturbed quantities, the corresponding
perturbed geometrical part of the field equations is

$\displaystyle\delta W^{i}_{\
j}=\frac{f_{T}}{2a^{2}}\delta^{ik}\left(h^{\prime\prime}_{kj}+[2+\nu]\mathcal{H}h^{\prime}_{kj}-\nabla^{2}h_{kj}\right),$
(151)

giving us the GW propagation equation in Fourier space as shown in eq. (38).

## References

  * Clifton et al. [2012] T. Clifton, P. G. Ferreira, A. Padilla, C. Skordis,  Modified Gravity and Cosmology,  Phys. Rept. 513 (2012) 1–189.
  * Aghanim et al. [2020] N. Aghanim, et al. (Planck),  Planck 2018 results. VI. Cosmological parameters,  Astron. Astrophys. 641 (2020) A6.
  * Tröster et al. [2020] T. Tröster, et al.,  Cosmology from large-scale structure,  Astronomy $\&$ Astrophysics 633 (2020) L10.
  * Dodelson and Schmidt [2021] S. Dodelson, F. Schmidt, Modern Cosmology, Elsevier, 2021\.
  * Sahni [2004] V. Sahni,  5 dark matter and dark energy,  in: The Physics of the Early Universe, Springer Berlin Heidelberg, 2004, pp. 141–179. URL: https://doi.org/10.1007/978-3-540-31535-3_5. doi:10.1007/978-3-540-31535-3_5.
  * Riess et al. [1998] A. G. Riess, et al.,  Observational evidence from supernovae for an accelerating universe and a cosmological constant,  The Astronomical Journal 116 (1998) 1009–1038.
  * Bull et al. [2016] P. Bull, et al.,  Beyond $\Lambda$CDM: Problems, solutions, and the road ahead,  Physics of the Dark Universe 12 (2016) 56–99.
  * Fazlollah Hajkarim and Jürgen Schaffner-Bielich [2020] Fazlollah Hajkarim and Jürgen Schaffner-Bielich,  Thermal history of the early universe and primordial gravitational waves from induced scalar perturbations,  Physical Review D 101 (2020).
  * Perivolaropoulos and Skara [2022] L. Perivolaropoulos, F. Skara,  Challenges for $\Lambda$CDM: An update,  New Astronomy Reviews 95 (2022) 101659.
  * Martin [2012] J. Martin,  Everything you always wanted to know about the cosmological constant problem (but were afraid to ask),  Comptes Rendus Physique 13 (2012) 566–665.
  * Velten et al. [2014] H. E. S. Velten, R. F. vom Marttens, W. Zimdahl,  Aspects of the cosmological “coincidence problem”,  The European Physical Journal C 74 (2014).
  * Weinberg [1989] S. Weinberg,  The cosmological constant problem,  Reviews of Modern Physics 61 (1989) 1–23.
  * Valentino et al. [2021] E. D. Valentino, et al.,  Snowmass2021 - letter of interest cosmology intertwined II: The hubble constant tension,  Astroparticle Physics 131 (2021) 102605.
  * Abdalla et al. [2022] E. Abdalla, et al.,  Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies,  Journal of High Energy Astrophysics 34 (2022) 49–211.
  * Tián [2020] S. X. Tián,  Cosmological consequences of a scalar field with oscillating equation of state: A possible solution to the fine-tuning and coincidence problems,  Physical Review D 101 (2020).
  * Chimento et al. [2003] L. P. Chimento, A. S. Jakubi, D. Pavón, W. Zimdahl,  Interacting quintessence solution to the coincidence problem,  Physical Review D 67 (2003).
  * Kang et al. [2020] G.-Z. Kang, D.-S. Zhang, L. Jun, H.-S. Zong,  Fine tuning problem of the cosmological constant in a generalized randall-sundrum model,  Chinese Physics C 44 (2020) 125102.
  * Feng and Li [2014] C.-J. Feng, X.-Z. Li,  Towards a realistic solution of the cosmological constant fine-tuning problem by higgs inflation,  Physical Review D 90 (2014).
  * Bisabr [2010] Y. Bisabr,  Coincidence problem in $f(R)$ gravity models,  Physical Review D 82 (2010).
  * Rudra [2015] P. Rudra,  Towards a possible solution for the coincidence problem: f(g) gravity as background,  International Journal of Modern Physics D 24 (2015) 1550013.
  * Krishnan et al. [2022] C. Krishnan, R. Mohayaee, E. O. Colgáin, M. M. Sheikh-Jabbari, L. Yin,  Hints of FLRW breakdown from supernovae,  Phys. Rev. D 105 (2022) 063514.
  * Krishnan et al. [2021] C. Krishnan, R. Mohayaee, E. O. Colgáin, M. M. Sheikh-Jabbari, L. Yin,  Does Hubble tension signal a breakdown in FLRW cosmology?,  Class. Quant. Grav. 38 (2021) 184001.
  * Poulin et al. [2023] V. Poulin, T. L. Smith, T. Karwal,  The Ups and Downs of Early Dark Energy solutions to the Hubble tension: a review of models, hints and constraints circa 2023 (2023).
  * Di Valentino and Melchiorri [2022] E. Di Valentino, A. Melchiorri,  Neutrino Mass Bounds in the Era of Tension Cosmology,  Astrophys. J. Lett. 931 (2022) L18.
  * Di Valentino et al. [2022] E. Di Valentino, S. Gariazzo, C. Giunti, O. Mena, S. Pan, W. Yang,  Minimal dark energy: Key to sterile neutrino and Hubble constant tensions?,  Phys. Rev. D 105 (2022) 103511.
  * Aldrovandi and Pereira [2013] R. Aldrovandi, J. G. Pereira, Teleparallel Gravity, volume 173, Springer, Dordrecht, 2013\. doi:10.1007/978-94-007-5143-9.
  * Bahamonde et al. [2021] S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. L. Said, J. Mifsud, E. Di Valentino,  Teleparallel Gravity: From Theory to Cosmology (2021).
  * Krssak et al. [2019] M. Krssak, R. van den Hoogen, J. Pereira, C. Böhmer, A. Coley,  Teleparallel theories of gravity: illuminating a fully invariant approach,  Class. Quant. Grav. 36 (2019) 183001.
  * Cai et al. [2016] Y.-F. Cai, S. Capozziello, M. De Laurentis, E. N. Saridakis,  $f(T)$ teleparallel gravity and cosmology,  Rept. Prog. Phys. 79 (2016) 106901.
  * Kadam et al. [2022] S. A. Kadam, J. L. Said, B. Mishra,  Accelerating cosmological models in $f(T,B)$ gravitational theory (2022).
  * Franco et al. [2020] G. A. R. Franco, C. Escamilla-Rivera, J. Levi Said,  Stability analysis for cosmological models in $f(T,B)$ gravity,  Eur. Phys. J. C 80 (2020) 677.
  * Mirza and Oboudiat [2017] B. Mirza, F. Oboudiat,  Constraining f(t) gravity by dynamical system analysis,  Journal of Cosmology and Astroparticle Physics 2017 (2017) 011–011.
  * Rave-Franco et al. [2021] G. A. Rave-Franco, C. Escamilla-Rivera, J. L. Said,  Dynamical complexity of the teleparallel gravity cosmology,  Physical Review D 103 (2021).
  * Briffa et al. [2022] R. Briffa, C. Escamilla-Rivera, J. L. Said, J. Mifsud, N. L. Pullicino,  Impact of $H_{0}$ priors on f(t) late time cosmology,  The European Physical Journal Plus 137 (2022).
  * Nunes [2018] R. C. Nunes,  Structure formation in $f(T)$ gravity and a solution for $H_{0}$ tension,  Journal of Cosmology and Astroparticle Physics 2018 (2018) 052–052.
  * Nájera et al. [2022] S. Nájera, A. Aguilar, G. A. Rave-Franco, C. Escamilla-Rivera, R. A. Sussman,  Inhomogeneous solutions in $f(T,B)$ gravity,  International Journal of Geometric Methods in Modern Physics 19 (2022).
  * Bourakadi et al. [2022] K. E. Bourakadi, B. Asfour, Z. Sakhi, M. Bennai, T. Ouali,  Primordial black holes and gravitational waves in teleparallel gravity,  The European Physical Journal C 82 (2022).
  * Sahlu et al. [2020] S. Sahlu, J. Ntahompagaze, A. Abebe, D. F. Mota,  Inflationary constraints in teleparallel gravity theory,  International Journal of Geometric Methods in Modern Physics 18 (2020) 2150027\.
  * Raatikainen and Räsänen [2019] S. Raatikainen, S. Räsänen,  Higgs inflation and teleparallel gravity,  Journal of Cosmology and Astroparticle Physics 2019 (2019) 021–021.
  * Chakrabortty et al. [2021] M. Chakrabortty, N. Sk, S. Sanyal, A. K. Sanyal,  Inflation with f(t) teleparallel gravity,  The European Physical Journal Plus 136 (2021).
  * Ricciardone [2017] A. Ricciardone,  Primordial gravitational waves with LISA,  Journal of Physics: Conference Series 840 (2017) 012030.
  * Hohmann et al. [2019] M. Hohmann, L. Järv, M. Krššák, C. Pfeifer,  Modified teleparallel theories of gravity in symmetric spacetimes,  Physical Review D 100 (2019).
  * Socolovsky [2012a] M. Socolovsky,  Fiber bundles, connections, general relativity, and the einstein-cartan theory – part i,  Advances in Applied Clifford Algebras 22 (2012a) 837–872.
  * Socolovsky [2012b] M. Socolovsky,  Fiber bundles, connections, general relativity, and the einstein-cartan theory – part II,  Advances in Applied Clifford Algebras 22 (2012b) 873–909.
  * Blau [2011] M. Blau, Lecture Notes on General Relativity, Technical Report, Albert Einstein Center for Fundamental Physics, University of Bern, 2011. URL: http://www.blau.itp.unibe.ch/newlecturesGR.pdf.
  * Golovnev et al. [2017] A. Golovnev, T. Koivisto, M. Sandstad,  On the covariance of teleparallel gravity theories,  Classical and Quantum Gravity 34 (2017) 145013.
  * Carroll [2004] S. Carroll, Spacetime and Geometry: An Introduction to General Relativity, Addison Wesley, 2004\.
  * Sotiriou and Faraoni [2010] T. P. Sotiriou, V. Faraoni,  f(R) Theories Of Gravity,  Rev. Mod. Phys. 82 (2010) 451–497.
  * Faraoni [2008] V. Faraoni,  f(R) gravity: Successes and challenges,  in: 18th SIGRAV Conference, 2008\. arXiv:0810.2602.
  * Capozziello et al. [2018] S. Capozziello, M. Capriolo, M. Transirico,  The gravitational energy-momentum pseudotensor: the cases of $f(R)$ and $f(T)$ gravity,  Int. J. Geom. Meth. Mod. Phys. 15 (2018) 1850164.
  * Farrugia et al. [2020] G. Farrugia, J. Levi Said, A. Finch,  Gravitoelectromagnetism, Solar System Test and Weak-Field Solutions in $f(T,B)$ Gravity with Observational Constraints,  Universe 6 (2020) 34.
  * Farrugia et al. [2018] G. Farrugia, J. Levi Said, V. Gakis, E. N. Saridakis,  Gravitational Waves in Modified Teleparallel Theories,  Phys. Rev. D 97 (2018) 124064.
  * Bahamonde et al. [2015] S. Bahamonde, C. G. Böhmer, M. Wright,  Modified teleparallel theories of gravity,  Phys. Rev. D92 (2015) 104042.
  * Wald [1984] R. M. Wald, General Relativity, University of Chicago Pr., 1984.
  * Bahamonde et al. [2021] S. Bahamonde, V. Gakis, S. Kiorpelidi, T. Koivisto, J. L. Said, E. N. Saridakis,  Cosmological perturbations in modified teleparallel gravity models: boundary term extension,  The European Physical Journal C 81 (2021).
  * Piattella [2018] O. Piattella, Lecture Notes in Cosmology, Springer International Publishing, 2018. URL: https://doi.org/10.1007/978-3-319-95570-4. doi:10.1007/978-3-319-95570-4.
  * Hohmann [2021] M. Hohmann,  General cosmological perturbations in teleparallel gravity,  The European Physical Journal Plus 136 (2021).
  * Chen et al. [2011] S.-H. Chen, J. B. Dent, S. Dutta, E. N. Saridakis,  Cosmological perturbations in f(T) gravity,  Phys. Rev. D83 (2011) 023508.
  * Bahamonde et al. [2023] S. Bahamonde, K. F. Dialektopoulos, M. Hohmann, J. Levi Said, C. Pfeifer, E. N. Saridakis,  Perturbations in non-flat cosmology for f(T) gravity,  Eur. Phys. J. C 83 (2023) 193.
  * Nunes et al. [2018] R. C. Nunes, S. Pan, E. N. Saridakis,  New observational constraints on $f(T)$ gravity through gravitational-wave astronomy,  Phys. Rev. D 98 (2018) 104055.
  * Ryden [2017] B. Ryden, Introduction to Cosmology, Cambridge University Press, 2017.
  * Klose et al. [2022] P. Klose, M. Laine, S. Procacci,  Gravitational wave background from non-abelian reheating after axion-like inflation,  Journal of Cosmology and Astroparticle Physics 2022 (2022) 021.
  * Kundu [2012] S. Kundu,  Inflation with general initial conditions for scalar perturbations,  Journal of Cosmology and Astroparticle Physics 2012 (2012) 005--005.
  * Weinberg [2008] S. Weinberg, Cosmology, Cosmology, OUP Oxford, 2008. URL: https://books.google.com.mt/books?id=nqQZdg020fsC.
  * Baumann [2022] D. Baumann, Cosmology, Cambridge University Press, Cambridge, England, 2022.
  * Boyle and Steinhardt [2008] L. A. Boyle, P. J. Steinhardt,  Probing the early universe with inflationary gravitational waves,  Physical Review D 77 (2008).
  * Bahamonde and Capozziello [2017] S. Bahamonde, S. Capozziello,  Noether Symmetry Approach in $f(T,B)$ teleparallel cosmology,  Eur. Phys. J. C77 (2017) 107.
  * Bahamonde et al. [2019] S. Bahamonde, U. Camci, S. Capozziello,  Noether symmetries and boundary terms in extended teleparallel gravity cosmology,  Classical and Quantum Gravity 36 (2019) 065013.
  * Rezazadeh et al. [2016] K. Rezazadeh, A. Abdolmaleki, K. Karami,  Power-law and intermediate inflationary models in f(t)-gravity,  Journal of High Energy Physics 2016 (2016).
  * dos Santos et al. [2022] F. B. M. dos Santos, J. E. Gonzalez, R. Silva,  Observational constraints on f(t) gravity from model-independent data,  The European Physical Journal C 82 (2022).
  * Escamilla-Rivera and Levi Said [2020] C. Escamilla-Rivera, J. Levi Said,  Cosmological viable models in $f(T,B)$ theory as solutions to the $H_{0}$ tension,  Class. Quant. Grav. 37 (2020) 165002.
  * Li et al. [2018] C. Li, Y. Cai, Y.-F. Cai, E. N. Saridakis,  The effective field theory approach of teleparallel gravity, $f(t)$ gravity and beyond,  Journal of Cosmology and Astroparticle Physics 2018 (2018) 001--001.
  * Nesseris et al. [2013] S. Nesseris, S. Basilakos, E. N. Saridakis, L. Perivolaropoulos,  Viable $f(t)$ models are practically indistinguishable from lcdm,  Physical Review D 88 (2013).
  * Klose et al. [2022] P. Klose, M. Laine, S. Procacci,  Gravitational wave background from vacuum and thermal fluctuations during axion-like inflation,  Journal of Cosmology and Astroparticle Physics 2022 (2022) 020.
  * Hogan [1986] C. J. Hogan,  Gravitational radiation from cosmological phase transitions,  Monthly Notices of the Royal Astronomical Society 218 (1986) 629--636.
  * Durrer [2010] R. Durrer,  Gravitational waves from cosmological phase transitions,  Journal of Physics: Conference Series 222 (2010) 012021.
  * Grojean and Servant [2007] C. Grojean, G. Servant,  Gravitational waves from phase transitions at the electroweak scale and beyond,  Physical Review D 75 (2007).
  * Caprini and Durrer [2006] C. Caprini, R. Durrer,  Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields,  Physical Review D 74 (2006).
  * Benetti et al. [2022] M. Benetti, L. L. Graef, S. Vagnozzi,  Primordial gravitational waves from NANOGrav: A broken power-law approach,  Physical Review D 105 (2022).
  * Papanikolaou et al. [2023] T. Papanikolaou, C. Tzerefos, S. Basilakos, E. N. Saridakis,  No constraints for f(T) gravity from gravitational waves induced from primordial black hole fluctuations,  Eur. Phys. J. C 83 (2023) 31.
  * Tzerefos et al. [2023] C. Tzerefos, T. Papanikolaou, E. N. Saridakis, S. Basilakos,  Scalar induced gravitational waves in modified teleparallel gravity theories (2023).