USTC-ICTS/PCFT-22-27

# Irregular universe in the Nieh-Yan modified teleparallel gravity

Mingzhe Li  Interdisciplinary Center for Theoretical Study, University of
Science and Technology of China, Hefei, Anhui 230026, China  Peng Huanwu
Center for Fundamental Theory, University of Science and Technology of China,
Hefei, Anhui 230026, China  Haomin Rao  School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  University of Chinese Academy of Sciences, 100190 Beijing,
China

###### Abstract

The Nieh-Yan modified teleparallel gravity is a model which modifies the
general relativity equivalent teleparallel gravity by a coupling between the
Nieh-Yan density and an axion-like field. This model predicts parity
violations in the gravitational waves if the axion-like field has a non-
trivial background, and more importantly it is ghost free and avoids the
pathologies presented in other parity-violating gravity models. The
cosmological dynamics and perturbations of the Nieh-Yan modified teleparallel
gravity have been investigated in detail, but all these previous
investigations rely on the symmetry requirement that in the background
universe both the metric and affine connection are homogeneous and isotropic.
In this paper we relax the symmetry constraint on the connection and leave it
arbitrary at the beginning, after all the cosmological principle only needs
the metric of the background spacetime to meet the symmetry requirement. We
find a new flat universe solution for the Nieh-Yan modified teleparallel
gravity, for which the background dynamics itself is unchanged but the
perturbations around it present a new feature that the scalar and tensor
perturbations are coupled together at the linear level. The implications of
this peculiar feature in primordial perturbations from inflation are also
discussed.

##  I introduction

Stimulated by the experimental detections of gravitational waves (GWs) ligo1 ;
ligo2  and the developments in the cosmic microwave background radiation (CMB)
experiments CMB1 ; CMB2 , parity violating gravities attracted lots of
interests in recent years. A famous and frequently studied parity violating
gravity model is the so-called Chern-Simons modified gravity CSgravity1 ;
CSgravity2  which within the framework of Riemannian geometry modifies general
relativity (GR) by a gravitational Chern-Simons term. The Chern-Simons
modified gravity predicts the difference between the amplitudes of the left-
and right-handed polarized components of gravitational waves, i.e., the so-
called amplitude birefringence phenomenon. However, this model was found to
suffer from the problem of vacuum instability because one of the circularly
polarized components of GWs becomes a ghost at high frequencies CSgravity3 .
Further extensions Crisostomi:2017ugk ; Gao:2019liu ; Zhao:2019xmm  to this
model did not circumvent this difficulty because in these extended models the
pathological behavior still appear at high energy scales, as shown in Ref.
Bartolo:2020gsh . It is very difficult to have a ghost-free parity violating
gravity model within the framework of Riemannian geometry.

Successful parity violating gravity models are available if we go beyond the
Riemannian geometry. For example, the Nieh-Yan modified teleparallel gravity
(NYTG) PVtele1 ; PVtele2  is constructed within the framework of the
teleparallel gravity (TG) Tele ; tele2021 , where the gravity is identified
with the spacetime torsion in stead of the curvature. One may have a GR
equivalent TG model Maluf:2013gaa  (we may call it TGR). The NYTG model
PVtele1 ; PVtele2  modifies TGR slightly by the anomalous coupling
$\theta\mathcal{T}\widetilde{\mathcal{T}}$ between an axion-like field
$\theta(x)$ and the Nieh-Yan density Nieh:1981ww :
$\mathcal{T}\widetilde{\mathcal{T}}\equiv(1/2)\varepsilon^{\mu\nu\rho\sigma}\mathcal{T}^{\lambda}_{~{}\mu\nu}\mathcal{T}_{\lambda\rho\sigma}$,
where $\mathcal{T}^{\lambda}_{~{}\mu\nu}$ is the torsion tensor,
$\varepsilon^{\mu\nu\rho\sigma}$ is Levi-Civita tensor which relates the
totally antisymmetric symbol $\epsilon^{\mu\nu\rho\sigma}$ and the determinant
of the metric $g$ through the equation
$\varepsilon^{\mu\nu\rho\sigma}=\epsilon^{\mu\nu\rho\sigma}/\sqrt{-g}$. The
Nieh-Yan density is parity-odd, so at the background with
$\partial_{\mu}\theta\neq 0$, the Nieh-Yan coupling term
$\theta\mathcal{T}\widetilde{\mathcal{T}}$ violates the parity symmetry
spontaneously. The NYTG model has been applied to cosmology in Refs. PVtele1 ;
PVtele2 , where it was found that this model predicts a difference between the
propagating velocities of the left- and right-handed polarized components of
GWs, i.e., the so-called velocity birefringence phenomenon. More importantly,
through detailed investigations on the cosmological perturbations, it was
shown in Refs. PVtele1 ; PVtele2  that the NYTG model is ghost-free. Recently,
this model was found to be compatible with the results of most local tests in
the Solar System at the post-Newtonian order Rao:2021azn ; Qiao:2021fwi , the
upper limit on its model parameters by the GWs data of LIGO/Virgo
Collaboration was obtained in Ref. Wu:2021ndf , and the enhancement of
primordial GWs during inflation due to the velocity birefringence of NYTG
model and its implications in the air-based GWs experiments were studied in
Ref. Cai:2021uup . Other recent studies on parity violating gravities can be
found in Refs. Li:2021mdp ; Gong:2021jgg ; Hohmann:2022wrk ; Tong:2022cdz ;
Li:2022vtn ; Zhang:2022xmm ; Zhu:2022dfq ; Zhu:2022uoq ; Filho:2022yrk ;
Qiao:2022mln ; Cai:2022lec ; Chen:2022wtz .

In all the previous studies of the cosmological applications of the NYTG
model, both the metric and the affine connection of the background universe
are required to be homogeneous and isotropic at the beginning. The spacetime
under this strong symmetry constraint is called the regular universe in this
paper. The background solutions of the regular universe have been well studied
within the TG framework Hohmann:2019nat ; Hohmann:2020zre ; Coley:2022qug ,
and are universally applicable to almost all TG models. In fact these
solutions have been frequently adopted by different authors, e.g.,
Myrzakulov:2010vz ; Cai:2011tc ; Cai:2015emx ; Hohmann:2022wrk  111 Actually,
the cosmological background solution whose tetrad is
$e^{A}_{~{}\mu}=\mathrm{diag}(1,a,a,a)$ or
$e^{A}_{~{}\mu}=\mathrm{diag}(a,a,a,a)$ under the Weitzenböck gauge is the
regular flat universe. However, most of the earlier literature did not clearly
point out that the selection of such a tetrad under the Weitzenböck gauge
actually requires the connection to satisfy the same symmetry of the metric..

However, the cosmological principle only needs the metric of the background
universe to meet the high symmetry requirement. In the Riemannian geometry,
once we impose this symmetry requirement on the metric, the connection (i.e.,
the Christoffel symbol) satisfies the same symmetry requirement automatically.
In TG models, the symmetry constraint on the affine connection is independent
of the one on the metric. If one drops this extra constraint on the connection
and leaves it arbitrary at the beginning, there will be final solutions for
which the connection is neither homogeneous nor isotropic. We call the
universe which has a homogeneous and isotropic metric and a non-homogeneous
and non-isotropic affine connection the irregular universe. So far the
irregular universe has rarely aroused research interest, only a few
literatures have initially studied the flat irregular universe (or irregular
Minkowski spacetime) solutions in $f(\mathbb{T})$ gravity models
Bejarano:2017akj ; Bejarano:2019fii ; BeltranJimenez:2020fvy ;
Golovnev:2020nln . The irregular universe does not violate the cosmological
principle, but questions are in coming: What features and new physical
phenomena could exist in the irregular universe? Or might the irregular
universe have properties that are clearly contradictory to experiments so that
only the regular universe is physically feasible? These questions deserve
detailed studies for any TG models.

In this paper, we will study the irregular universe in the NYTG model.
Firstly, we will obtain a more general flat universe solution than those in
Refs. PVtele1 ; PVtele2  by solving the equations of motion of the NYTG model
directly under the condition that only the metric is required to be
homogeneous and isotropic. By analyzing the symmetry of the connection, we
will show that the flat universe we obtain is generally an irregular flat
universe, and in special cases it reduces back to a regular universe. We will
also show that even in the irregular flat universe, the background equations
in the NYTG model are exactly the same as those in GR. Secondly, we will study
the linear cosmological perturbations around the irregular flat universe. We
will find that tensor perturbations and scalar perturbations are coupled at
the linear perturbation level. This is a peculiar feature that distinguishes
the irregular universe from the regular universe in the NYTG model. We
speculate that this peculiar feature is caused by the fact that the interior
space does not satisfy the homogeneity and isotropy in the irregular universe.
Finally, we will study the primordial fluctuations generated by slow-roll
inflation in the regular and irregular flat universes. We will show that the
primordial fluctuations of left- and right-handed GWs are different whether in
the regular universe or in the irregular universe. We will also show that
there is a strong statistical correlation between primordial scalar
fluctuations and primordial tensor fluctuations generated by slow-roll
inflation in the irregular universe.

This paper is organized as follows. In Sec. II, we briefly introduce the TG
theory and the NYTG model. In Sec. III, we study spatially flat cosmological
background solutions that only requires the metric to be homogeneous and
isotropic in the NYTG model. In Sec. IV, through the quadratic actions for
scalar, vector, and tensor perturbations, we investigate linear perturbations
around the regular and irregular flat universes. In Sec. V, we apply our
result to the early universe and discuss briefly the primordial perturbations
generated by slow-roll inflation.

In this paper, we adopt the unit $8\pi G=1$, and use the signature $(+,-,-,-)$
for the metric. The tensor indices of the interior space are denoted by
$A,B,C,...=0,1,2,3$ and by $a,b,c,...=1,2,3$ when limiting to spatial
components. They are lowered and raised by the Minkowski metric $\eta_{AB}$
and its inverse $\eta^{AB}$. The spacetime tensor indices are denoted by Greek
$\mu,\nu,\rho,...=0,1,2,3$ and by Latin $i,j,k,...=1,2,3$ when limiting to
spatial components. They are lowered and raised by the spacetime metric
$g_{\mu\nu}$ and its inverse $g^{\mu\nu}$. The antisymmetric symbol
$\epsilon^{\mu\nu\rho\sigma}$ has the properties:
$\epsilon^{0ijk}=\epsilon^{ijk}\equiv\epsilon_{ijk}$, and $\epsilon^{123}=1$.
In addition, we distinguish the spacetime affine connection
$\hat{\Gamma}^{\rho}_{~{}\mu\nu}$ and its associated covariant derivative
$\hat{\nabla}$ from the Levi-Civita connection ${\Gamma}^{\rho}_{~{}\mu\nu}$
and its associated covariant derivative ${\nabla}$ respectively.

##  II TG theory and the NYTG model

The TG theory can be considered as a constrained metric-affine theory. It is
formulated in a spacetime endowed with a metric $g_{\mu\nu}$ and an affine
connection $\hat{\Gamma}^{\rho}_{~{}\mu\nu}$, which is curvature free and
metric compatible,

$\hat{R}^{\sigma}_{~{}\rho\mu\nu}=\partial_{\mu}\hat{\Gamma}^{\sigma}_{~{}\nu\rho}-\partial_{\nu}\hat{\Gamma}^{\sigma}_{~{}\mu\rho}+\hat{\Gamma}^{\sigma}_{~{}\mu\lambda}\hat{\Gamma}^{\lambda}_{~{}\nu\rho}-\hat{\Gamma}^{\sigma}_{~{}\nu\lambda}\hat{\Gamma}^{\lambda}_{~{}\mu\rho}=0~{},~{}\hat{\nabla}_{\rho}g_{\mu\nu}=\partial_{\rho}g_{\mu\nu}-\hat{\Gamma}^{\lambda}_{~{}\rho\mu}g_{\lambda\nu}-\hat{\Gamma}^{\lambda}_{~{}\rho\nu}g_{\mu\lambda}=0~{}.$
(1)

Without curvature and nonmetricity, in the TG theory the gravity is identified
with spacetime torsion
$\mathcal{T}^{\rho}_{~{}\mu\nu}=2\hat{\Gamma}_{\,[\mu\nu]}^{\rho}$. One can
also describe the TG theory using the language of the tetrad $e^{A}_{~{}\mu}$
and the spin connection $\omega^{A}_{~{}B\mu}$. They relates the metric
$g_{\mu\nu}$ and the affine connection $\hat{\Gamma}^{\rho}_{~{}\mu\nu}$
through the following relations

$g_{\mu\nu}=\eta_{AB}e^{A}_{~{}\mu}e^{B}_{~{}\nu}~{}~{},~{}~{}\hat{\Gamma}^{\rho}_{~{}\mu\nu}=e_{A}^{~{}\,\rho}(\partial_{\mu}e^{A}_{~{}\nu}+\omega^{A}_{~{}B\mu}e^{B}_{~{}\nu})~{}.$
(2)

The torsion tensor is written as

$\mathcal{T}^{\rho}_{~{}\mu\nu}=2e_{A}^{~{}\,\rho}(\partial_{[\mu}e^{A}_{~{}\nu]}+\omega^{A}_{~{}B[\mu}e^{B}_{~{}\nu]})~{}.$
(3)

The teleparallel constraints (1) dictate that the spin connection can be in
general expressed as

$\omega_{~{}B\mu}^{A}=(\Lambda^{-1})^{A}_{~{}C}\partial_{\mu}\Lambda_{~{}B}^{C}~{},$
(4)

where $\Lambda^{A}_{~{}B}$ is arbitrary element of Lorentz transformation
matrix which is position dependent and satisfies the relation
$\eta_{AB}\Lambda^{A}_{~{}C}\Lambda^{B}_{~{}D}=\eta_{CD}$ at any spacetime
point. Therefore, the tetrad $e^{A}_{~{}\mu}$ and the Lorentz matrix
$\Lambda^{A}_{~{}B}$ can be regarded as the basic variables of the TG theory.
In this way, the teleparallel constraints (1) are automatically satisfied.

The TGR model, as the GR equivalent TG model, has the following action,

$\displaystyle
S_{TGR}=\frac{1}{2}\int\mathrm{d}^{4}x~{}{|e|}\,\mathbb{T}\equiv\int\mathrm{d}^{4}x~{}{|e|}\left(-\frac{1}{2}\mathcal{T}_{\mu}\mathcal{T}^{\mu}+\frac{1}{8}\mathcal{T}_{\alpha\beta\mu}\mathcal{T}^{\alpha\beta\mu}+\frac{1}{4}\mathcal{T}_{\alpha\beta\mu}\mathcal{T}^{\beta\alpha\mu}\right)~{},$
(5)

where ${|e|}=\sqrt{-g}$ is the determinant of the tetrad, $\mathbb{T}$ is the
torsion scalar, and $\mathcal{T}_{\mu}=\mathcal{T}^{\alpha}_{~{}~{}\mu\alpha}$
is the torsion vector. Since we have the identity
$-{R}(e)=\mathbb{T}+2{\nabla}_{\mu}\mathcal{T}^{\mu}$, the action (5) is
identical to the Einstein-Hilbert action up to a surface term, where the
curvature scalar $R(e)$ is defined by the Levi-Civita connection and
considered as being fully constructed from the metric, and in turn from the
tetrad. Since the surface term in the action does not affect the equations of
motion, we say that the TGR is equivalent to GR at the level of the equations
of motion.

The NYTG model PVtele1 ; PVtele2  modifies the TGR model by introducing the
coupling

$S_{NY}=\frac{c}{4}\int\mathrm{d}^{4}x~{}{|e|}~{}\theta\,\mathcal{T}\widetilde{\mathcal{T}}~{},$
(6)

between an axion-like field $\theta$ and the Nieh-Yan density
$\mathcal{T}\widetilde{\mathcal{T}}$. The coupling constant $c$ is
dimensionless. Generally we should also consider its own dynamics of the
axion-like field and take other matter into account, so the full action of the
NYTG model is

$\displaystyle
S_{NYTG}=\int\mathrm{d}^{4}x~{}{|e|}\left[\frac{1}{2}\mathbb{T}+\frac{c}{4}\,\theta\,\mathcal{T}\widetilde{\mathcal{T}}+\frac{1}{2}\nabla_{\mu}\theta\nabla^{\mu}\theta-V(\theta)\right]+S_{m}~{}.$
(7)

Other matter with the action $S_{m}$ is assumed to be coupled to spacetime
minimally through the tetrad. At the background in which the axion-like field
has non-zero spacetime derivatives, the Nieh-Yan coupling term breaks parity
spontaneously. Because only the first-order derivatives of the basic variables
appears in the action, the NYTG model can avoid the Ostrogradski ghost mode,
which is expected to be originated from higher-order derivatives in the action
Woodard:2006nt .

As with most modified TG theories, the NYTG model apparently has two kinds of
gauge symmetries: diffeomorphism invariance and local Lorentz invariance. The
latter transformation makes the following change:

$e^{A}_{~{}\mu}\rightarrow(L^{-1})^{A}_{~{}B}e^{B}_{~{}\mu}~{},~{}\Lambda^{A}_{~{}B}\rightarrow\Lambda^{A}_{~{}C}L^{C}_{~{}B}~{},$
(8)

where $L^{A}_{~{}B}(x)$ are the element of Lorentz matrix. We would like to
use different notations to distinguish two kinds of Lorentz matrices:
$\Lambda^{A}_{~{}B}(x)$ is used to express the spin connection as in Eq. (4),
but $L^{A}_{~{}B}(x)$ represents the local transformation that makes a shift
from one local frame to another. Transformation (8) can be expressed in terms
of tetrad and spin connections as

$e^{A}_{~{}\mu}\rightarrow(L^{-1})^{A}_{~{}B}e^{B}_{~{}\mu}~{},~{}\omega^{A}_{~{}B\mu}\rightarrow(L^{-1})^{A}_{~{}C}\omega^{C}_{~{}D\mu}L^{D}_{~{}B}+(L^{-1})^{A}_{~{}C}\partial_{\mu}L^{C}_{~{}B}~{}.$
(9)

It is easy to prove that the metric $g_{\mu\nu}$ and torsion tensor
$\mathcal{T}_{~{}\mu\nu}^{\rho}$ are invariant under the local Lorentz
transformation (8), as is the action (7). Due to the local Lorentz invariance,
one can choose the gauge $\Lambda^{A}_{~{}B}=\delta^{A}_{~{}B}$, i.e.,
$\omega^{A}_{~{}B\mu}=0$. This is the Weitzenböck connection, which has been
frequently adopted in the literature. In addition, there is another symmetry
hidden in the NYTG model. The Nieh-Yan term (6) can be integrated by parts as

$S_{NY}=-\frac{c}{2}\int\mathrm{d}^{4}x~{}\eta_{AB}\epsilon^{\mu\nu\rho\sigma}(\partial_{\mu}\theta)(\Lambda^{A}_{~{}C}e^{C}_{~{}\nu})\partial_{\rho}(\Lambda^{B}_{~{}D}e^{D}_{~{}\sigma})~{}.$
(10)

It can be seen that the Nieh-Yan term (6) is invariant under the following
transformation

$(\Lambda^{A}_{~{}C}e^{C}_{~{}\mu})\rightarrow
L^{A}_{~{}B}(\theta)(\Lambda^{B}_{~{}C}e^{C}_{~{}\mu})\,,$ (11)

where $L^{A}_{~{}B}(\theta)$ is Lorentz matrix that depends only on axion-like
field $\theta$. Note that $\Lambda^{A}_{~{}C}e^{C}_{~{}\mu}$ is invariant
under transformation (8). Due to the Lorentz symmetry (8), the transformation
(11) can always be attributed to the fact that the tetrad $e^{A}_{~{}\mu}$
remains unchanged while the Lorentz matrix $\Lambda^{A}_{~{}B}$ undergoes a
Lorentz transformation. Obviously the metric and the action of TGR model are
invariant under such a transformation. So the total action of the NYTG model
is invariant under the transformation (11).

The equations of motion follow from the variation of the action (7) with
respect to $e^{A}_{~{}\mu}$ and $\Lambda^{A}_{~{}B}$ separately

$\displaystyle G^{\mu\nu}+N^{\mu\nu}$ $\displaystyle=$ $\displaystyle
T^{\mu\nu}+T^{\mu\nu}_{\theta}~{},$ (12) $\displaystyle N^{[\mu\nu]}$
$\displaystyle=$ $\displaystyle 0~{},$ (13)

where
$N^{\mu\nu}=(c/2)\varepsilon^{\mu\lambda\rho\sigma}\partial_{\lambda}\theta\,\mathcal{T}^{\nu}_{~{}\rho\sigma}$,
$G^{\mu\nu}$ is the Einstein tensor, $T^{\mu\nu}=-(2/\sqrt{-g})(\delta
S_{m}/\delta g_{\mu\nu})$ and
$T^{\mu\nu}_{\theta}=[V(\theta)-\nabla_{\alpha}\theta\nabla^{\alpha}\theta/2]g^{\mu\nu}+\nabla^{\mu}\theta\nabla^{\nu}\theta$
are the energy-momentum tensors for the matter and the axion-like field
$\theta$ respectively. Similar to most modified TG models, the equation of
motion (13) from the variation of $\Lambda^{A}_{~{}B}$ is not independent of
Eq. (12), it is just the antisymmetric part of the latter. As explained in
Ref. PVtele2 , this is due to the local Lorentz invariance of the action, any
change caused by $\delta\Lambda^{A}_{~{}B}$ can always be equivalent to the
change caused by $\delta e^{A}_{~{}\mu}$, so requiring the action to take the
extremum under $\delta e^{A}_{~{}\mu}$ already includes the case where the
action takes the extremum under $\delta\Lambda^{A}_{~{}B}$. There is another
equation following from the variation of the action (7) with respect to
$\theta$,

${\square}\theta+V^{(1)}-\frac{c}{4}\mathcal{T}\widetilde{\mathcal{T}}=0~{},$
(14)

where ${\square}=g^{\mu\nu}{\nabla}_{\mu}{\nabla}_{\nu}$ and
$V^{(n)}=\mathrm{d}^{n}V(\theta)/\mathrm{d}\theta^{n}$. All of these equations
of motion are consistent with the Bianchi identity $\nabla_{\mu}G^{\mu\nu}=0$
and the covariant conservation law $\nabla_{\mu}T^{\mu\nu}=0$.

Also in Refs. PVtele1 ; PVtele2 , the cosmological perturbations of the NYTG
model were analyzed in detail. It was found that the NYTG model makes a
difference between the propagating velocities of the left- and right-handed
polarized components of GWs, but makes no difference between their amplitudes.
This phenomenon is called velocity birefringence, which is a clear physical
signal of parity violation. More importantly, the NYTG model was confirmed to
be ghost free through the quadratic action of cosmological perturbations.

It is worth mentioning that the Nieh-Yan density
$\mathcal{T}\widetilde{\mathcal{T}}$ is not the only parity-odd term within
the TG framework. A more general model including all the parity-odd terms
which are quadratic in the torsion tensor was considered in Ref. PVtele3 . But
then it was found in Ref. PVtele4  that this more general model suffers from
the problem of ghost instability again, unless it completely reduces to the
NYTG model. Therefore, within the TG framework, for all parity-odd terms which
are quadratic in the torsion tensor, only the Nieh-Yan density
$\mathcal{T}\widetilde{\mathcal{T}}$ can avoid the ghost instability. This
means the NYTG model is robust to some extent.

##  III Irregular flat universe in the NYTG model

So far all the studies on the cosmological applications of the NYTG model only
considered the regular universe as the background, that means both the metric
and the affine connection are constrained to be homogeneous and isotropic.
This constraint may be too strong, after all the cosmological principle which
is supported by current observations only needs the metric of the background
spacetime to meet the high symmetry requirement. In this paper, we will drop
the symmetry requirement on the connection and leave it arbitrary at the
beginning. After this relaxation, it is expected that the NYTG model will have
more interesting cosmological background solutions. We are interested in the
irregular universe solutions in which the metric homogeneous and isotropic but
the connection is neither homogeneous nor isotropic. For simplicity, we will
only consider the spatially flat universe.

In flat universe, the metric can be expressed in rectangular coordinate as

$\mathrm{d}s^{2}=g_{\mu\nu}\mathrm{d}x^{\mu}\mathrm{d}x^{\nu}=a^{2}\left(\mathrm{d}\eta^{2}-\delta_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j}\right)~{},$
(15)

where $a=a(\eta)$ is the scale factor of the universe, $\eta$ is the conformal
time. This is the Friedmann-Robertson-Walker (FRW) metric. There are 6 Killing
vector fields $\\{\xi_{I}^{\mu},I=1,2...6\\}$ in flat universe, which can be
expressed as

$\xi^{\mu}_{I}=\delta^{\,\mu}_{~{}I}~{},~{}\xi^{\mu}_{I+3}=\epsilon_{Iij}\delta^{\mu}_{~{}i}x^{j}~{},~{}~{}~{}~{}I=1,2,3$
(16)

where $\xi^{\mu}_{1},\xi^{\mu}_{2},\xi^{\mu}_{3}$ are Killing vector fields
representing the symmetry of spatial translation, and
$\xi^{\mu}_{4},\xi^{\mu}_{5},\xi^{\mu}_{6}$ are Killing vector fields
representing the symmetry of spatial rotation. One can prove that the FRW
metric satisfies the condition: $\mathcal{L}_{\xi_{I}}g_{\mu\nu}=0$, where
$\mathcal{L}_{\xi_{I}}$ is the Lie derivative along the Killing vector field
$\xi^{\mu}_{I}$. This reflects the result that the metric is homogeneous and
isotropic. One can also prove that
$\mathcal{L}_{\xi_{I}}\Gamma^{\rho}_{~{}\mu\nu}=0$ for the Levi-Civita
connection $\Gamma^{\rho}_{~{}\mu\nu}$, which is automatically homogeneous and
isotropic. This is why we do not need to pay extra attention to the symmetry
of the connection within the framework of Riemannian geometry.

###  III.1 Regular flat universe

For TG models, even the metric is determined, the affine connection is still
arbitrary to some extent. Usually, as suggested in Refs Hohmann:2019nat ;
Hohmann:2020zre ; Coley:2022qug , a further constraint was imposed that
requires the connection is also homogeneous and isotropic, that is,

$\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{\rho}_{~{}\mu\nu}=\hat{\nabla}_{\mu}\hat{\nabla}_{\nu}\,\xi^{\rho}_{I}-\hat{\nabla}_{\mu}(\mathcal{T}^{\rho}_{~{}\nu\sigma}\xi^{\sigma}_{I})=0~{}.$
(17)

Although $\hat{\Gamma}^{\rho}_{~{}\mu\nu}$ is coordinate dependent, the Lie
derivative of $\hat{\Gamma}^{\rho}_{~{}\mu\nu}$ does not depend on the
coordinate. Hence the condition (17) is unambiguous. Combining Eqs. (15) and
(17) selected the regular flat universe solution in which the tetrad
$e^{A}_{\mu}$ and Lorentz matrix $\Lambda^{A}_{~{}B}$ have the following
forms:

$e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}~{},~{}\Lambda^{A}_{~{}B}=\mathring{\Lambda}^{A}_{~{}B}~{},$
(18)

where $\mathring{\Lambda}^{A}_{~{}B}$ is a global Lorentz matrix, which does
not depend on spacetime. All other solutions satisfying Eqs. (15) and (17)
differ from the solution (18) only by Lorentz transformation (8), so they are
physically equivalent to the solution (18). The above process does not depend
on a specific TG theory, so the solution (18) is generally applicable to most
TG theories.

For the NYTG model, the solution (18) can automatically satisfy the constraint
$N^{[\mu\nu]}=0$, so the solution (18) is compatible with the NYTG model.
Furthermore, solution (18) leads to $N^{\mu\nu}=0$ and
$\mathcal{T}\widetilde{\mathcal{T}}=0$, which means that the Nieh-Yan term has
no effect on the regular flat universe background. Therefore, the background
equations of the regular flat universe are exactly the same as those of GR
PVtele1 ; PVtele2 .

###  III.2 Irregular flat universe

To look for the irregular universe solution, we should give up the constraint
(17) on the connection. After this relaxation, the connection is left to be
determined by the equation of motion.

In a flat universe, we can always simply find the non-zero components of
$G^{\mu\nu}$, $T^{\mu\nu}$ and $T^{\mu\nu}_{\theta}$ as

$G^{00}=\frac{3\mathcal{H}^{2}}{a^{4}}~{},~{}T^{00}=\frac{\rho}{a^{2}}~{},~{}T_{\theta}^{00}=\frac{\rho_{\theta}}{a^{2}}~{},~{}G^{ij}=-\frac{2\mathcal{H}^{\prime}+\mathcal{H}^{2}}{a^{4}}\delta_{ij}~{},~{}T^{ij}=\frac{p}{a^{2}}\delta_{ij}~{},~{}T_{\theta}^{ij}=\frac{p_{\theta}}{a^{2}}\delta_{ij}~{},$
(19)

where $\mathcal{H}=a^{\prime}/a$ is the conformal Hubble rate, prime
represents the derivative with respect to the conformal time $\eta$,
$\rho_{\theta}=\theta^{\prime 2}/\left(2a^{2}\right)+V$ and
$p_{\theta}=\theta^{\prime 2}/\left(2a^{2}\right)-V$ are the energy density
and pressure of the $\theta$ field, and $\rho$ and $p$ denote the energy
density and pressure of other matter. Thanks to the Lorentz symmetry (8), we
can always reduce the tetrad to the simple form
$e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}$ in flat universe. In order to facilitate
further analysis, we decompose the independent non-zero components of spin
connections $\omega^{A}_{~{}B\mu}$ as follows

$\displaystyle\delta^{a}_{~{}i}\omega^{0}_{~{}a0}=\mathcal{U}_{i}~{}~{},~{}~{}\delta^{i}_{~{}a}\delta^{b}_{~{}j}\omega^{a}_{~{}bk}=\Sigma\epsilon_{ijk}+\epsilon_{ijl}\Sigma_{kl}+\Sigma_{i}\delta_{jk}-\Sigma_{j}\delta_{ik}~{},$
$\displaystyle\delta^{i}_{~{}a}\delta^{b}_{~{}j}\omega^{a}_{~{}b0}=\epsilon_{ijk}\mathcal{V}_{k}~{}~{},~{}~{}\delta^{a}_{~{}i}\omega^{0}_{~{}aj}=\sigma\delta_{ij}+\sigma_{ij}+\epsilon_{ijk}\sigma_{k}~{},$
(20)

where $\Sigma_{ij}$ and $\sigma_{ij}$ are symmetric and traceless spatial
tensors. In the above decomposition we have exploited the property
$\omega_{AB\mu}=-\omega_{BA\mu}$ due to $\hat{\nabla}_{\rho}g_{\mu\nu}=0$.
Note that the variables
$\sigma,\Sigma,\mathcal{U}_{i},\mathcal{V}_{i},\sigma_{i},\Sigma_{i},\sigma_{ij},\Sigma_{ij}$
are not completely independent because we have not yet imposed
$\hat{R}^{\sigma}_{~{}\rho\mu\nu}=0$ on the spin connection. Combining
$e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}$ and Eq. (III.2), $N^{\mu\nu}$ can be
obtained as

$N^{00}=0~{}~{},~{}~{}N^{0i}=0~{}~{},~{}~{}N^{i0}=\frac{2c\theta^{\prime}}{a^{4}}\sigma_{i}~{}~{},~{}~{}N^{ij}=\frac{c\theta^{\prime}}{a^{4}}\left(2\Sigma\delta_{ij}-\Sigma_{ij}+\epsilon_{ijk}\Sigma_{k}\right)~{}.$
(21)

In order for Eqs. (12) and (13) to hold, there must be

$\sigma_{i}=0~{}~{},~{}~{}\Sigma_{i}=0~{}~{},~{}~{}\Sigma_{ij}=0~{}~{},~{}~{}\Sigma=\Sigma(\eta)~{}.$
(22)

Combining $e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}$, Eqs. (III.2) and (22), Nieh-
Yan density can be obtained as

$\mathcal{T}\widetilde{\mathcal{T}}=\frac{24\Sigma}{a^{2}}(\mathcal{H}-\sigma)~{}.$
(23)

In order for Eq. (14) to hold, the Nieh-Yan density
$\mathcal{T}\widetilde{\mathcal{T}}$ can only be a function of time $\eta$, so
$\sigma=\sigma(\eta)$ when $\Sigma\neq 0$.

Combining Eqs. (III.2) and (22), $\hat{R}^{\sigma}_{~{}\rho\mu\nu}=0$ gives

$\displaystyle\mathcal{S}^{\prime}_{ij}-\mathcal{U}_{i,j}+\epsilon_{ijk}\Sigma\,\mathcal{U}_{k}+\epsilon_{ikl}\mathcal{S}_{jk}\mathcal{V}_{l}=0~{},$
(24)
$\displaystyle\Sigma^{\prime}\delta_{ij}-\mathcal{V}_{i,j}+\epsilon_{ijk}\Sigma\,\mathcal{U}_{k}-\epsilon_{ikl}\mathcal{S}_{jk}\mathcal{U}_{l}=0~{},$
(25)
$\displaystyle\epsilon_{ikl}\mathcal{S}_{lj,k}+\Sigma(\mathcal{S}_{ij}-\mathcal{S}_{kk}\delta_{ij})=0~{},$
(26)
$\displaystyle\epsilon_{inm}\mathcal{S}_{jn}\mathcal{S}_{km}-\Sigma^{2}\epsilon_{ijk}=0~{},$
(27)

where $\mathcal{S}_{ij}=\sigma\delta_{ij}+\sigma_{ij}$ and the subscript
“$,i$” represents a derivative with respect to $x^{i}$. The trace of Eq. (26)
gives

$\sigma\Sigma=0~{}.$ (28)

This means that at least one of $\sigma$ and $\Sigma$ is zero. If $\sigma=0$,
the equation after the Hodge duality of the ”$j,k$” index in Eq. (27) can be
decomposed as follows according to the trace part and the traceless part:

$6\,\Sigma^{2}+\sigma_{ij}\sigma_{ij}=0~{}~{},~{}~{}\sigma_{ik}\sigma_{jk}-\frac{1}{3}(\sigma_{kl}\sigma_{kl})\delta_{ij}=0~{}.$
(29)

The solution of Eq. (29) is $\Sigma=0,\sigma_{ij}=0$. This means that Eqs.
(27) and (28) must give

$\Sigma=0~{}.$ (30)

Combining Eqs. (22) and (30) gives $N^{\mu\nu}=0$ and
$\mathcal{T}\widetilde{\mathcal{T}}=0$, which means that the Nieh-Yan term has
no effect even on the irregular flat universe background. Therefore, the
background equations of the irregular flat universe are exactly the same as
those of GR. This is a somewhat unexpected result. But the fact that Nieh-Yan
term has no effect on the background does not mean that it has no effect on
the perturbations. In order to analyze the perturbations, we need to first
find the background solution of the irregular flat universe.

Substituting Eq. (30) into Eqs. (24), (25), (26) and (27), we get

$\displaystyle\mathcal{S}^{\prime}_{ij}-\mathcal{U}_{i,j}+\epsilon_{ikl}\mathcal{S}_{jk}\mathcal{V}_{l}=0~{},$
(31)
$\displaystyle\mathcal{V}_{i,j}+\epsilon_{ikl}\mathcal{S}_{jk}\mathcal{U}_{l}=0~{},$
(32) $\displaystyle\epsilon_{ikl}\mathcal{S}_{lj,k}=0~{},$ (33)
$\displaystyle\epsilon_{inm}\mathcal{S}_{jn}\mathcal{S}_{km}=0~{},$ (34)

Although there are more equations than variables, this does not mean that Eqs.
(31), (32), (33) and (34) have no solution. It can be verified that the
following are the solution of Eqs. (31), (32), (33) and (34)

$\displaystyle\mathcal{S}_{ij}=v_{i}v_{j}f(\eta)F^{(1)}(\vec{v}\cdot\vec{x})~{},$
$\displaystyle\mathcal{V}_{i}=g_{a}(\eta)\alpha^{a}_{i}(\eta,\vec{x})-h_{a}(\eta)\beta^{a}_{i}(\eta,\vec{x})~{},$
$\displaystyle\mathcal{U}_{i}=h_{a}(\eta)\alpha^{a}_{i}(\eta,\vec{x})+g_{a}(\eta)\beta^{a}_{i}(\eta,\vec{x})+v_{i}f^{(1)}(\eta)F(\vec{v}\cdot\vec{x})~{},$
(35)

where

$\displaystyle\alpha^{a}_{i}(\eta,\vec{x})=\cosh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]\delta_{ai}+\frac{v_{a}v_{i}}{v^{2}}\Big{(}1-\cosh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]\Big{)}~{},$
$\displaystyle\beta^{a}_{i}(\eta,\vec{x})=\epsilon_{aij}\frac{v_{j}}{v}\sinh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]~{},$

where $v_{1},v_{2},v_{3}$ are constant parameters,
$v=\sqrt{\delta^{ij}v_{i}v_{j}}$, $\vec{v}\cdot\vec{x}=v_{i}x^{i}$,
$f(\eta),g_{a}(\eta),h_{a}(\eta)$ are arbitrary smooth function of conformal
time $\eta$, $F(\vec{v}\cdot\vec{x})$ is arbitrary smooth function of
$\vec{v}\cdot\vec{x}$, $f^{(n)}(\eta)$ is the $n$ derivative of $f(\eta)$ with
respect to conformal time $\eta$, and $F^{(n)}(\vec{v}\cdot\vec{x})$ is the
$n$ derivative of $F(\vec{v}\cdot\vec{x})$ with respect to
$\vec{v}\cdot\vec{x}$.

Putting solutions (22), (30) and (III.2) into the decomposition (III.2), the
spin connection $\omega^{A}_{~{}B\mu}$ when the tetrad is
$e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}$ can be obtained as

$\displaystyle\omega^{a}_{~{}00}=\omega^{0}_{~{}a0}=h_{c}(\eta)\alpha^{c}_{a}(\eta,\vec{x})+g_{c}(\eta)\beta^{c}_{a}(\eta,\vec{x})+v_{a}f^{(1)}(\eta)F(\vec{v}\cdot\vec{x})~{},$
$\displaystyle\omega^{a}_{~{}b0}=\epsilon_{abi}\left[g_{c}(\eta)\alpha^{c}_{i}(\eta,\vec{x})-h_{c}(\eta)\beta^{c}_{i}(\eta,\vec{x})\right]~{},$
$\displaystyle\omega^{0}_{~{}ai}=\omega^{a}_{~{}0i}=v_{a}v_{i}f(\eta)F^{(1)}(\vec{v}\cdot\vec{x})~{},~{}~{}\omega^{a}_{~{}bi}=0~{}.$
(36)

It can be verified that the spin connection (III.2) does satisfy the
teleparallel constraints (1). Due to the symmetry (11), not every
$h_{I}(\eta)$ and $g_{I}(\eta)$ represent a physically inequivalent solution.
In order to see this better, we perform a Lorentz transformation (9) on the
above solution. The transformation matrix $L^{A}_{~{}B}$ is

$\displaystyle
L^{0}_{~{}0}=\cosh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]~{},~{}L^{0}_{~{}a}=L^{a}_{~{}0}=\frac{v_{a}}{v}\sinh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]~{},$
$\displaystyle
L^{a}_{~{}b}=\delta_{ab}+\frac{v_{a}v_{b}}{v^{2}}\Big{(}\cosh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]-1\Big{)}~{},$
(37)

Then, the tetrad $\tilde{e}^{A}_{~{}\mu}=L^{A}_{~{}B}e^{B}_{~{}\mu}$ and the
corresponding spin connection $\tilde{\omega}^{A}_{~{}B\mu}$ are

$\displaystyle\tilde{e}^{0}_{~{}0}=a\cosh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]~{},\tilde{e}^{a}_{~{}0}=\delta^{ai}\tilde{e}^{0}_{~{}i}=a\frac{v_{a}}{v}\sinh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]~{},$
$\displaystyle\tilde{e}^{a}_{~{}i}=a\left[\delta_{ai}+\frac{v_{a}v_{i}}{v^{2}}\big{(}\cosh\left[vf(\eta)F(\vec{v}\cdot\vec{x})\right]-1\big{)}\right]~{},$
$\displaystyle\tilde{\omega}^{a}_{~{}00}=\tilde{\omega}^{0}_{~{}a0}=h_{a}(\eta)~{},~{}\tilde{\omega}^{a}_{~{}b0}=\epsilon_{abc}g_{b}(\eta)~{},~{}\tilde{\omega}^{A}_{~{}Bi}=0~{}.$
(38)

It can be verified that the metric $g_{\mu\nu}$ and connection
$\hat{\Gamma}^{\rho}_{~{}\mu\nu}$ given by solution (III.2) are the same as
those given by the tetrad $e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}$ and the spin
connection (III.2). Since the solution (III.2) satisfies the teleparallel
constraints (1), the spin connection $\tilde{\omega}^{A}_{~{}B\mu}$ in the
solution (III.2) can be expressed by a Lorentz matrix
$\tilde{\Lambda}^{A}_{~{}B}(\eta,\vec{x})$. And $\tilde{\omega}^{A}_{~{}Bi}=0$
means that
$\tilde{\Lambda}^{A}_{~{}B}(\eta,\vec{x})=\tilde{\Lambda}^{A}_{~{}B}(\eta)$.
So taking different $h_{a}(\eta)$ and $g_{a}(\eta)$ is actually taking
different $\tilde{\Lambda}^{A}_{~{}B}(\eta)$. Since $\theta=\theta(\eta)$ in
the cosmological background, different $\tilde{\Lambda}^{A}_{~{}B}(\eta)$ can
be converted to each other through the Lorentz transformation
$\tilde{\Lambda}^{A}_{~{}B}(\eta)\rightarrow
L^{A}_{~{}C}(\theta)\tilde{\Lambda}^{C}_{~{}B}(\eta)$. Therefore, the
solutions with different $h_{a}(\eta)$ and $g_{a}(\eta)$ can be transformed
into each other by transformation (11), so they are physically equivalent. In
this case, we only need to consider the simplest case below, that is, the case
where $h_{a}(\eta)=g_{a}(\eta)=0$, so that the solution (III.2) can be
simplified to

$\displaystyle e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}~{},$
$\displaystyle\omega^{a}_{~{}00}=\omega^{0}_{a0}=v_{a}f^{(1)}(\eta)F(\vec{v}\cdot\vec{x})~{},~{}\omega^{a}_{~{}b0}=0~{},$
$\displaystyle\omega^{a}_{~{}0i}=\omega^{0}_{ai}=v_{a}v_{i}f(\eta)F^{(1)}(\vec{v}\cdot\vec{x})~{},~{}\omega^{a}_{~{}bi}=0~{}.$
(39)

The solution (III.2) can be expressed by the tetrad $e^{A}_{~{}\mu}$ and the
Lorentz matrix $\Lambda^{A}_{~{}B}$ as

$e^{A}_{~{}\mu}=a\delta^{A}_{~{}\mu}~{},~{}~{}\Lambda=\mathring{\Lambda}\cdot\exp\Big{[}f(\eta)F(\vec{v}\cdot\vec{x})\,v_{a}\mathbb{K}^{a}\Big{]}~{},$
(40)

where $\mathring{\Lambda}$ is a spacetime independent Lorentz matrix, and
$\mathbb{K}^{1},\mathbb{K}^{2},\mathbb{K}^{3}$ are the boost matrices whose
expression are

$\mathbb{K}^{1}=\left(\begin{array}[]{cccc}0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{array}\right)~{},~{}~{}\mathbb{K}^{2}=\left(\begin{array}[]{cccc}0&0&1&0\\\
0&0&0&0\\\ 1&0&0&0\\\
0&0&0&0\end{array}\right)~{},~{}~{}\mathbb{K}^{3}=\left(\begin{array}[]{cccc}0&0&0&1\\\
0&0&0&0\\\ 0&0&0&0\\\ 1&0&0&0\end{array}\right)~{}.$

Regardless of the functional form of $f(\eta)$ and $F(\vec{v}\cdot\vec{x})$,
it can be verified that the solution (40) always satisfies the teleparallel
constraints (1) and makes Eqs. (12) and (14) self-consistent. Note that the
boosted tetrad solutions similar to the solution (40) are also discussed in
Ref. Bejarano:2019fii ; BeltranJimenez:2020fvy ; Golovnev:2020nln . Putting
solution (40) into Eqs. (12) and (14), we can get

$3\mathcal{H}^{2}=a^{2}\left(\rho_{\theta}+\rho\right)~{},~{}~{}2\mathcal{H}^{\prime}+\mathcal{H}^{2}=-a^{2}\left(p_{\theta}+p\right)~{},~{}~{}\theta^{\prime\prime}+2\mathcal{H}\theta^{\prime}+a^{2}V^{(1)}=0~{}.$
(41)

The background equations are exactly the same as those of GR. This means that
the Nieh-Yan term has no effect even on the irregular flat universe
background. This is consistent with our analysis above.

Finally, let’s focus on the symmetry of the connection given by the solution
(40). The non-zero components of
$\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{\rho}_{~{}\mu\nu}$ given by the solution
(40) are

$\displaystyle\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{0}_{~{}0i}=\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{i}_{~{}00}=v_{I}v_{i}f^{(1)}(\eta)F^{(1)}(\vec{v}\cdot\vec{x})~{},$
$\displaystyle\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{0}_{~{}ij}=\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{i}_{~{}j0}=v_{I}v_{i}v_{j}f(\eta)F^{(2)}(\vec{v}\cdot\vec{x})~{},$
$\displaystyle\mathcal{L}_{\xi_{I+3}}\hat{\Gamma}^{0}_{~{}0i}=\mathcal{L}_{\xi_{I+3}}\hat{\Gamma}^{i}_{~{}00}=-\epsilon_{Iij}v_{j}f^{(1)}(\eta)F(\vec{v}\cdot\vec{x})+v_{i}\epsilon_{Ijk}v_{j}x^{k}f^{(1)}(\eta)F^{(1)}(\vec{v}\cdot\vec{x})~{},$
$\displaystyle\mathcal{L}_{\xi_{I+3}}\hat{\Gamma}^{0}_{~{}ij}=\mathcal{L}_{\xi_{I+3}}\hat{\Gamma}^{i}_{~{}j0}=2v_{(i}\epsilon_{j)Ik}v_{k}f(\eta)F^{(1)}(\vec{v}\cdot\vec{x})+v_{i}v_{j}\epsilon_{Ikl}v_{k}x^{l}f(\eta)F^{(2)}(\vec{v}\cdot\vec{x})~{},$
(42)

where $I=1,2,3$ in Eq. (III.2), and the subscript parentheses denotes the
symmetrization. The fact that
$\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{\rho}_{~{}\mu\nu}\neq 0$ indicates that
the spacetime connection given by the solution (40) is neither homogeneous nor
isotropic. So the solution (40) does represent an irregular flat universe.
When $v_{i}=0$ or $f(\eta)=0$ or $F(\vec{v}\cdot\vec{x})=0$, there is
$\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{\rho}_{~{}\mu\nu}=0$, and the solution
(40) dose reduce to the regular flat universe solution (18).

##  IV Perturbations around the irregular flat universe

In the previous section we studied the flat universe solution of the NYTG
model that only requires the metric to be homogeneous and isotropic. We found
that the Nieh-Yan term has no effect even on the irregular flat universe
background. In order to explore the effect of the Nieh-Yan term on the
irregular flat universe, we study the linear cosmological perturbations around
the irregular flat universe (40) in this section. For simplicity, we only
consider the case of $F(\vec{v}\cdot\vec{x})=\vec{v}\cdot\vec{x}$, which is
equivalent to requiring that the coefficients of the equations of linear
perturbations do not depend on the spatial coordinates $\vec{x}$ (see below
for details). And we also ignore other matter so that $S_{m}=0$.

We use the following parametrization for perturbed tetrad Izumi:2012qj :

$\displaystyle e^{0}_{\ 0}=a(1+A)~{},~{}e^{0}_{\
i}=a(\beta_{,i}+\beta_{i}^{V})~{},~{}e^{c}_{\
0}=a\delta_{ci}(\chi_{,i}+\chi_{i}^{V})~{},$ $\displaystyle e^{c}_{\
i}=a\delta_{cj}[(1-\psi)\delta_{ij}+\alpha_{,ij}+\alpha_{j,i}^{V}-\epsilon_{ijk}(\lambda_{,k}+\lambda_{k}^{V})+\frac{1}{2}h^{T}_{ij}]~{},$
(43)

So the perturbed metric components have the familiar forms:

$\displaystyle g_{00}=a^{2}(1+2A)~{},~{}g_{0i}=-a^{2}(B_{,i}+B_{i}^{V})~{},$
$\displaystyle
g_{ij}=-a^{2}[(1-2\psi)\delta_{ij}+2\alpha_{,ij}+\alpha_{i,j}^{V}+\alpha_{j,i}^{V}+h^{T}_{ij}]~{},$
(44)

where $B=\chi-\beta$ and $B^{V}_{i}=\chi^{V}_{i}-\beta^{V}_{i}$. Besides the
familiar scalar perturbations ($A,B,\psi,\alpha$), vector perturbations
($B_{i}^{V},\alpha^{V}_{i}$), and tensor perturbations $h^{T}_{ij}$ in the
metric, the parametrization of tetrad brings six extra variables, which are
scalar perturbation $\lambda,\chi+\beta$ and vector perturbation
$\lambda_{i}^{V},\chi_{i}^{V}+\beta_{i}^{V}$. All the vector perturbations are
transverse and denoted by the superscript $V$, both the tensor perturbations
are transverse and traceless and denoted by the superscript $T$. In addition,
the scalar field $\theta$ is decomposed as
$\theta(\eta,\vec{x})=\bar{\theta}(\eta)+\delta\theta(\eta,\vec{x})$.

Although we can perform a similar decomposition on the Lorentz matrix
$\Lambda^{A}_{~{}B}$ following the parametrization in Ref. PVtele2 , we do not
need to do so in this paper. Because we can always transform the perturbed
Lorentz matrix into the background Lorentz matrix in Eq. (40) through the
infinitesimal Lorentz transformation (8). In other words, we can always absorb
the perturbations of the Lorentz matrix $\Lambda^{A}_{~{}B}$ into the
perturbations of the tetrad $e^{A}_{~{}\mu}$ through the infinitesimal Lorentz
transformation (8), so that we only need to deal with the perturbations of the
the tetrad.

Due to the diffeomorphism invariance, it is safe to take the unitary gauge
$\delta\theta=0,~{}\alpha=0,~{}\alpha_{i}^{V}=0$. This simplifies the
calculations, for example, the gauge invariant scalar perturbation
$\zeta=-(\psi+\mathcal{H}\delta\theta/\theta^{\prime})$ representing the
curvature perturbation of the hypersurfaces of constant $\theta$ reduces to
$-\psi$ under the unitary gauge. Since both $\alpha$ and $\alpha_{i}^{V}$ are
perturbations which enter the metric, the perturbations $\alpha$,
$\alpha_{i}^{V}$ and $\delta\theta$ are invariant under the infinitesimal
Lorentz transformation (8). Therefore, the unitary gauge is compatible with
the operation of absorbing the perturbations of the Lorentz matrix into the
perturbations of the tetrad.

The non-isotropic nature of the background connection may lead to coupling of
scalar, vector and tensor perturbations. Therefore, when studying linear
perturbations around the irregular flat universe (40), we should not deal with
scalar, vector, or tensor perturbations individually, but should deal with all
perturbation variables simultaneously. In the following we choose $A$,
$\zeta$, $B$, $B_{i}^{V}$, $\beta_{i}=\beta_{,i}+\beta_{i}^{V}$,
$\lambda_{i}=\lambda_{,i}+\lambda_{i}^{V}$ and $h^{T}_{ij}$ as independent
variables, and we study the linear perturbations around the irregular flat
universe by means of quadratic action.

For the NYTG model (7) with $S_{m}=0$, one can directly obtain the quadratic
action as

$\displaystyle
S^{(2)}=\int\mathrm{d}^{4}x~{}a^{2}\bigg{\\{}6\mathcal{H}\zeta^{\prime}A-3{\zeta^{\prime}}^{2}-(2A+\zeta)\zeta_{,ii}-a^{2}VA^{2}+2(\zeta^{\prime}-\mathcal{H}A)B_{,ii}+\frac{1}{8}\left(h^{T\prime}_{ij}h^{T\prime}_{ij}-h^{T}_{ij,k}h^{T}_{ij,k}\right)$
$\displaystyle\quad\quad\quad\quad-\frac{1}{4}B_{i}^{V}B_{i,jj}^{V}+c\theta^{\prime}\Big{[}2\lambda_{i}\zeta_{,i}+\frac{1}{2}\epsilon^{ijk}(\beta_{i}\beta_{j,k}-\lambda_{i}\lambda_{j,k})+\hat{\mathcal{S}}_{ij}\lambda_{i}\beta_{j}-\frac{1}{2}\epsilon^{ijk}\mathcal{S}_{il}h^{T}_{jl}\beta_{k}-\frac{1}{8}\epsilon^{ijk}h^{T}_{il}h^{T}_{jl,k}\Big{]}\bigg{\\}}.~{}~{}~{}$
(45)

where $\mathcal{S}_{ij}=v_{i}v_{j}f(\eta)F^{(1)}(\vec{v}\cdot\vec{x})$ and
$\hat{\mathcal{S}}_{ij}=(v_{i}v_{j}-v^{2}\delta_{ij})f(\eta)F^{(1)}(\vec{v}\cdot\vec{x})$.
In general, the coefficients $\mathcal{S}_{ij}$ and $\hat{\mathcal{S}}_{ij}$
are explicitly dependent on the spatial coordinate $\vec{x}$ unless $F$ is a
linear function of $\vec{v}\cdot\vec{x}$. It means that generally the
evolution equations for the linear perturbations are neither homogeneous nor
isotropic. This is reasonable since we are studying the perturbations around
an irregular background. In this paper we only consider the case where
$F(\vec{v}\cdot\vec{x})=\vec{v}\cdot\vec{x}$ for simplicity 222The expression
of $F(\vec{v}\cdot\vec{x})$ can differ by a constant term, which does not
change the coefficients $\mathcal{S}_{ij}$ and $\hat{\mathcal{S}}_{ij}$. And a
constant factor of the difference of $F(\vec{v}\cdot\vec{x})$ can be absorbed
into $f(\eta)$.. In this way, $\mathcal{S}_{ij}$ and $\hat{\mathcal{S}}_{ij}$
are constant coefficients. So the evolution equations for the linear
perturbations are homogeneous. We will see later that when all constraints are
lifted, the quadratic action (IV) will become homogeneous and isotropic. In
addition, the terms $\hat{\mathcal{S}}_{ij}\lambda_{i}\beta_{j}$ and
$\epsilon^{ijk}\mathcal{S}_{il}h^{T}_{jl}\beta_{k}$ in the action (IV) show
that there is a coupling of scalar, vector and tensor perturbations. But such
coupling may be eliminated by the constraints imposed by the action (IV)
itself. Therefore, only after the constraints are lifted can we know whether
there is really a coupling of scalar, vector and tensor perturbations.

To further simplify the quadratic action, we change to the momentum space in
terms of Fourier transformations,

$\zeta(\eta,\vec{x})=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{\frac{3}{2}}}\,\zeta(\eta,\vec{k})\,\mathrm{e}^{\mathrm{i}\vec{k}\cdot\vec{x}}~{},$
(46)

and we also expand the variables $A$, $B$, ${\lambda}_{i}$, ${\beta}_{i}$ and
$h^{T}_{ij}$ in the same way. Note that we use the normal letter $\mathrm{i}$
for imaginary unit to distinguish it from the italic letter $i$ used for
spatial indices. The tensor perturbation $h^{T}_{ij}$ can be further expanded
as

$h^{T}_{ij}(\eta,\vec{k})=\sum_{\mathrm{A}}h_{\mathrm{A}}(\eta,\vec{k})\,\hat{e}_{ij}^{\mathrm{A}}(\vec{k})~{},$
(47)

where
$\\{\hat{e}^{\mathrm{A}}_{ij}(\vec{k}),~{}\mathrm{A}=\mathrm{L},\mathrm{R}\\}$
are circular polarization bases 333 Note that the choice of circular
polarization bases is not unique, $\hat{e}^{\mathrm{A}}_{ij}(\vec{k})$ can be
rotated along the $\vec{k}$-axis while maintaining all the properties of the
circular polarization bases. For the case where there is a constant vector
$\vec{v}\neq 0$ on the background, we can always choose the circular
polarization bases to satisfy
$v_{i}v_{j}\hat{e}^{\mathrm{A}}_{ij}(\vec{k})=(v^{2}/\sqrt{2})\sin^{2}\vartheta$,
where $\vartheta$ is the angle between $\vec{k}$ and $\vec{v}$. This choice
maximally simplifies the quadratic action (57), so we adopt this choice in
this paper.  satisfying
$\hat{k}^{l}\epsilon_{lik}\hat{e}_{jk}^{\mathrm{A}}(\vec{k})=\mathrm{i}\mathfrak{p}_{\mathrm{A}}\hat{e}_{ij}^{\mathrm{A}}(\vec{k})$,
where $\hat{k}$ is the unit vector of $\vec{k}$,
$\mathfrak{p}_{\mathrm{L}}=-1$ and $\mathfrak{p}_{\mathrm{R}}=1$. Note that we
use the normal letter $\mathrm{A}$ for the left- and right- hand indices to
distinguish it from the italic letter $A$ used to represent the tetrad
indices. The quadratic action in the momentum space can be expressed as

$\displaystyle
S^{(2)}=\int\mathrm{d}\eta\int\mathrm{d}^{3}k~{}a^{2}\bigg{\\{}6\mathcal{H}\zeta^{\prime}A^{*}-3\zeta^{*\prime}\zeta^{\prime}+k^{2}(2A+\zeta)\zeta^{*}+2k^{2}(\mathcal{H}A-\zeta^{\prime})B^{*}$
$\displaystyle\quad\quad\quad\quad-a^{2}VA^{*}A+\frac{1}{4}k^{2}B_{i}^{V*}B_{i}^{V}+\frac{1}{4}\sum_{\mathcal{A}}\left[h_{\mathrm{A}}^{*\prime}h_{\mathrm{A}}^{\prime}-(k^{2}-c\theta^{\prime}\mathfrak{p}_{\mathrm{A}}k)h^{*}_{\mathrm{A}}h_{\mathrm{A}}\right]$
$\displaystyle\quad\quad\quad\quad+c\theta^{\prime}\Big{[}2\mathrm{i}k^{i}\lambda^{*}_{i}\zeta+\frac{\mathrm{i}}{2}\epsilon^{ijk}k^{i}(\beta^{*}_{j}\beta_{k}-\lambda^{*}_{j}\lambda_{k})+\hat{\mathcal{S}}_{ij}\lambda^{*}_{i}\beta_{j}-\frac{1}{2}\beta^{*}_{i}\Big{(}\sum_{\mathrm{A}}\mathcal{S}^{\mathrm{A}}_{i}h_{\mathrm{A}}\Big{)}\Big{]}\bigg{\\}},$
(48)

where
$\mathcal{S}^{\mathrm{A}}_{i}(\vec{k})=\epsilon_{ijk}{\mathcal{S}}_{jl}\hat{e}^{\mathrm{A}}_{kl}(\vec{k})$.
It can be seen that $A$, $B$, $B_{i}^{V}$, $\lambda_{i}$ and $\beta_{i}$ are
all non-dynamical fields and the variations of the action (IV) with them lead
to the following constraints:

$\displaystyle B^{V}_{i}=0~{},$ (49)
$\displaystyle\mathcal{H}A-\zeta^{\prime}=0~{},$ (50) $\displaystyle
3\mathcal{H}\zeta^{\prime}+k^{2}\zeta-a^{2}VA+\mathcal{H}k^{2}B=0~{},$ (51)
$\displaystyle\epsilon_{ijk}k^{j}\lambda_{k}-\mathrm{i}\hat{S}_{ij}\beta_{j}+2k^{i}\zeta=0~{},$
(52)
$\displaystyle-\hat{S}_{ij}\lambda_{j}+\mathrm{i}\epsilon_{ijk}k^{j}\beta_{k}+\frac{1}{2}\sum_{A}\mathcal{S}^{\mathrm{A}}_{i}h_{\mathrm{A}}=0~{}.$
(53)

For the regular flat universe case with $v_{i}=0$ or $f(\eta)=0$, there are
$\hat{S}_{ij}=0$ and $\mathcal{S}^{\mathrm{A}}_{i}=0$, so the solution of Eqs.
(49), (50), (51), (52) and (53) is

$\zeta=0~{},~{}A=0~{},~{}B=0~{},~{}B_{i}^{V}=0~{},~{}\lambda_{i}=\mathrm{i}k^{i}\lambda~{},~{}\beta_{i}=\mathrm{i}k^{i}\beta~{},$
(54)

where $\lambda$ and $\beta$ are arbitrary scalar perturbations. Substituting
the Eq. (54) back into the action (IV), the action (IV) can be simplified as

$S^{(2)}=\int\mathrm{d}\eta\int\mathrm{d}^{3}k\,\frac{a^{2}}{4}\sum_{\mathrm{A}}\left(|h_{\mathrm{A}}^{\prime}|^{2}-\omega_{\mathrm{A}}^{2}|h_{\mathrm{A}}|^{2}\right)~{},$
(55)

where
$\omega^{2}_{\mathrm{A}}=k^{2}-c\theta^{\prime}\mathfrak{p}_{\mathrm{A}}k$. It
can be seen that there is no scalar dynamical degree of freedom at the linear
perturbation level. This is a bit strange because the action (7) clearly shows
that there is a scalar dynamical degree of freedom. Further research in Ref.
PVtele2  shows that the missing scalar dynamical degree of freedom reappears
in the regular curved universe. The phenomenon of degrees of freedom being
hidden under special background also appears in $f(\mathbb{T})$ gravity
Ong:2013qja  and massive gravity DeFelice:2012mx . This implies that such a
special background is likely to suffer from strong coupling issue
Delhom:2022vae . It can also be seen that the modified dispersion relation
$\omega^{2}_{\mathrm{A}}$ is helicity dependent. This means that GWs with
different helicities will have different propagation velocities. This
phenomenon is called velocity birefringence, which is a direct reflection of
the parity violation in the NYTG model. These results are consistent with the
results in Refs. PVtele1 ; PVtele2  444The subtle difference in the dispersion
relation $\omega^{2}_{\mathrm{A}}$ is due to the difference between expanding
by $e^{\mathrm{i}\vec{k}\cdot\vec{x}}$ and expanding by
$e^{-\mathrm{i}\vec{k}\cdot\vec{x}}$ in the Fourier transformation..

For the irregular flat universe case with $v_{i}\neq 0$ and $f(\eta)\neq 0$,
the solution of Eqs. (49), (50), (51), (52) and (53) is

$\displaystyle A=\zeta^{\prime}/\mathcal{H}~{},~{}B=-\left[\theta^{\prime
2}\zeta^{\prime}+2k^{2}\mathcal{H}\zeta\right]/2k^{2}\mathcal{H}^{2}~{},~{}B^{V}_{i}=0~{},~{}$
$\displaystyle\lambda_{i}=\Big{(}\frac{2\cos\vartheta}{kv\sin^{2}\vartheta}\epsilon_{ijk}k^{j}v_{k}\Big{)}\zeta-\frac{\mathrm{i}}{2\sqrt{2}k}k^{i}\Big{(}\sum_{\mathrm{A}}\mathfrak{p}_{\mathrm{A}}h_{\mathrm{A}}\Big{)}~{},$
$\displaystyle\beta_{i}=\Big{(}\frac{2\mathrm{i}}{v^{2}f(\eta)\sin^{2}\vartheta}k^{i}+\frac{2\mathrm{i}vf(\eta)\cos\vartheta}{k\sin^{2}\vartheta}v_{i}\Big{)}\zeta+\frac{\mathrm{i}vf(\eta)\cos\vartheta}{2\sqrt{2}k}v_{i}\Big{(}\sum_{\mathrm{A}}h_{\mathrm{A}}\Big{)}~{},$
(56)

where $\vartheta$ is the angle between $\vec{k}$ and $\vec{v}$. Substituting
the above results back into the action (IV), the action (IV) can be simplified
as

$S^{(2)}=\int\mathrm{d}\eta\int\mathrm{d}^{3}k\,\bigg{\\{}\frac{z^{2}}{2}\left(|\zeta^{\prime}|^{2}-k^{2}|\zeta|^{2}\right)+\frac{a^{2}}{4}\sum_{\mathrm{A}}\left(|h_{\mathrm{A}}^{\prime}|^{2}-\omega_{\mathrm{A}}^{2}|h_{\mathrm{A}}|^{2}\right)-\frac{ca^{2}\theta^{\prime}k}{\sqrt{2}}\zeta^{*}\Big{(}\sum_{\mathrm{A}}\mathfrak{p}_{\mathrm{A}}h_{\mathrm{A}}\Big{)}\bigg{\\}}~{},$
(57)

where $z^{2}=a^{2}\theta^{\prime 2}/\mathcal{H}^{2}$. For the action (57), the
following points need to be emphasized. Firstly, it can be seen that there is
indeed a scalar dynamical degree of freedom, which again verifies that there
is a scalar dynamical degree of freedom hidden under the regular flat universe
at the linear perturbation level. Secondly, there are two tensor dynamics
degrees of freedom and the dispersion relation $\omega_{\mathrm{A}}^{2}$ is
helicity dependent, as is the case for the regular universe. This means that
the velocity birefringence phenomenon of GWs also exists in the irregular
universe. Thirdly, it is surprising that $v_{i}$ and $f(\eta)$ are completely
cancelled in the step of lifting the constraints, so that the action (57) no
longer depends on $v_{i}$ and $f(\eta)$. This makes the case of
$v_{i}=0,f(\eta)=0$ not the limit of the case of $v_{i}\rightarrow
0,f(\eta)\rightarrow 0$. This is somewhat analogous to the case where a
massless photon is not the limit of a photon with mass tends to zero. Fourth,
it can be seen that the coefficients in the action (57) are homogeneous and
isotropic. This means that the evolution equations of the scalar perturbation
$\zeta$ and the tensor perturbations $h_{\mathrm{A}}$ are homogeneous and
isotropic. It should be emphasized that such a property depends on the fact
that we only consider the simple case with
$F(\vec{v}\cdot\vec{x})=\vec{v}\cdot\vec{x}$. In the next section we will see
that this results in isotropic power spectra for perturbations. This is
consistent with observations since the current data from the CMB found no
statistically significant evidence for direction dependence of the power
spectra Planck:2018jri . More general considerations of the function $F$ and
the constrains on it from observational data through detailed analysis deserve
further studies. These considerations are model dependent and beyond the scope
of this paper. Finally, it can be seen that even after the constraints are
lifted, there is still a coupling of scalar and tensor degrees of freedom.
This is a feature that neither in the regular flat universe nor in the regular
curved universe. This means that scalar perturbations and tensor perturbations
can influence each other at the linear perturbation level. This can be seen
more clearly from the perspective of the equations of motion. From the action
(57), the linear equations of $\zeta$ and $h_{\mathrm{A}}$ can be obtained as

$\displaystyle\zeta^{\prime\prime}+2\frac{z^{\prime}}{z}\zeta^{\prime}+k^{2}\zeta+\frac{ca^{2}\theta^{\prime}k}{\sqrt{2}z^{2}}\Big{(}\sum_{\mathrm{A}}\mathfrak{p}_{\mathrm{A}}h_{\mathrm{A}}\Big{)}=0~{},$
(58) $\displaystyle
h_{\mathrm{A}}^{\prime\prime}+2\mathcal{H}h_{\mathrm{A}}^{\prime}+\omega_{\mathrm{A}}^{2}h_{\mathrm{A}}+\sqrt{2}c\theta^{\prime}\mathfrak{p}_{\mathrm{A}}k\zeta=0~{}.$
(59)

Eq. (58) shows that the tensor perturbations $h_{\mathrm{A}}$ can be used as a
source of the scalar perturbation $\zeta$. The scalar perturbation $\zeta$ can
be excited when left- and right- handed GWs have different amplitudes or
phases. And Eq. (59) shows that the scalar perturbation $\zeta$ can be used as
a source of the tensor perturbations $h_{\mathrm{A}}$. It is worth noting that
the source of the tensor perturbations $h_{\mathrm{A}}$ caused by $\zeta$ is
helicity-dependent, that is, the excitation effects caused by $\zeta$ on the
left- and right-handed GWs are different.

##  V Primordial fluctuations generated by inflation

In the previous section, we preliminarily studied the the linear perturbations
around the regular and irregular flat universe, and obtained the quadratic
action after the constraints was lifted. In this section, we will
preliminarily study the primordial fluctuations generated by slow-roll
inflation in the regular and irregular flat universe.

###  V.1 The case of the regular universe

For the case of regular universe, the quadratic action (55) can be expressed
as

$S^{(2)}=\int\mathrm{d}\eta\int\mathrm{d}^{3}k\,\frac{a^{2}}{2}\sum_{\mathrm{A}}\left[\Big{|}\frac{1}{\sqrt{2}}h_{\mathrm{A}}^{\prime}\Big{|}^{2}-\left(k^{2}-c\theta^{\prime}\mathfrak{p}_{\mathrm{A}}k\right)\Big{|}\frac{1}{\sqrt{2}}h_{\mathrm{A}}\Big{|}^{2}\right]~{}.$
(60)

Note that since there are only tensor degrees of freedom in the regular flat
universe at the linear perturbation level, a scalar field other than $\theta$
needs to be introduced to generate the primordial scalar perturbation PVtele1
; Cai:2021uup . In this subsection we do not consider the case of introducing
additional scalar fields, and we only focus on the tensor perturbations.

Next we consider the case of slow-roll inflation dominated by the axion-like
field $\theta$. Since the background equations of the regular flat universe
are exactly the same as those in GR, the background evolution during inflation
will be exactly the same as the case of slow-roll inflation in GR
Baumann:2009ds ; Wang:2013zva . So we don’t need to repeat the analysis of the
details of single scalar field inflation. We introduce two commonly used slow-
roll parameters

$\varepsilon\equiv-\frac{\dot{H}}{H^{2}}~{},~{}\delta\equiv\frac{\ddot{\theta}}{H\dot{\theta}}~{},$
(61)

where $H=\dot{a}/{a}=\mathcal{H}/a$ is the Hubble rate, the upper dot
represents the derivative with respect to the physical time $t$. We assume
$\varepsilon\sim|\delta|\ll 1$, $|\dot{\varepsilon}/H|\ll|\varepsilon|$ and
$|\dot{\delta}/H|\ll|\delta|$ during inflation. Under the slow-roll
approximation,

$\mathcal{H}\approx-\frac{1+\varepsilon}{\eta}~{}~{},~{}~{}\theta^{\prime}\approx\frac{\sqrt{2\varepsilon}}{\eta}~{}.$
(62)

Without loss of generality, in Eq. (62) we have assumed that the value of
$\theta$ decreases during inflation.

Next, by combining Eqs (60) and (62), the correlation function of
$h_{\mathrm{A}}$ can be obtained through the process in Appendix C:

$\langle h_{\mathrm{A}}^{\dagger}h_{\mathrm{A}}\rangle\approx
H^{2}\mathrm{e}^{-\mathfrak{p}_{\mathrm{A}}\sqrt{\varepsilon/2}c\pi}k^{-(3+2\varepsilon)}~{},$
(63)

and $\langle h_{\mathrm{L}}^{\dagger}h_{\mathrm{R}}\rangle=0$. Through the
correlation functions (63), the power spectrum of the left- and right-handed
GWs can be obtained as

$\mathcal{P}_{\mathrm{A}}(k)=\frac{k^{3}}{\pi^{2}}\langle
h_{\mathrm{A}}^{\dagger}h_{\mathrm{A}}\rangle\approx\frac{H^{2}}{\pi^{2}}\mathrm{e}^{-\mathfrak{p}_{\mathrm{A}}\sqrt{\varepsilon/2}c\pi}k^{-2\varepsilon}~{}.$
(64)

The power spectrum of the tensor perturbations can be obtained as

$\mathcal{P}_{T}(k)=\mathcal{P}_{\mathrm{L}}(k)+\mathcal{P}_{\mathrm{R}}(k)\approx\frac{H^{2}}{\pi^{2}}\left[1+\cosh\left(\sqrt{\frac{\varepsilon}{2}}c\pi\right)\right]k^{-2\varepsilon}~{}.$
(65)

The relative different between the power spectrum of the left- and right-
handed GWs can be obtained as

$\Pi\equiv\frac{\mathcal{P}_{\mathrm{R}}-\mathcal{P}_{\mathrm{L}}}{\mathcal{P}_{\mathrm{R}}+\mathcal{P}_{\mathrm{L}}}\approx-\tanh\left(\sqrt{\frac{\varepsilon}{2}}c\pi\right)\approx-\sqrt{\frac{\varepsilon}{2}}c\pi~{}.$
(66)

$\Pi\neq 0$ means that the magnitudes of the primordial fluctuations of left-
and right-handed GWs are different. This is a clear physical signal of parity
violation. But this seems to contradict the conclusion in Refs. PVtele1 ;
PVtele2  that there is only velocity birefringence of GWs but no amplitude
birefringence of GWs in the NYTG model. The reason for this contradiction is
that $\theta^{\prime}$ is approximated as a constant in the analysis of the
evolution of GWs in Refs. PVtele1 ; PVtele2 . Of course, this approximation is
valid when studying the propagation of GWs in a slowly expanding universe.
However, $\theta^{\prime}=a\dot{\theta}\propto 1/\eta$ cannot be approximated
as a constant during the slow-roll inflation dominated by $\theta$. We know
that for a harmonic oscillator (the equation of motion is
$\ddot{x}+\omega^{2}x=0$), the amplitude of the harmonic oscillator can be
changed when the frequency $\omega$ is time-dependent. And when the time
dependence of $\theta^{\prime}$ is not negligible, the time dependence of
$\omega_{L}$ and $\omega_{R}$ will be different, resulting in different
effects on the amplitudes of left- and right-hand GWs. This is why the
magnitudes of the primordial fluctuations of left- and right-handed GWs
generated by slow-roll inflation in the regular flat universe are different.
If $\varepsilon\rightarrow 0$, it can be seen from Eq. (62) that
$\theta^{\prime}\approx 0$ can be approximated as a constant, and from Eq.
(66), it can be seen that $\Pi\rightarrow 0$ too, that is, the magnitudes of
the primordial fluctuation of the left- and right-handed GWs are the same.

Finally, let’s look at the case when the coupling constant $c\rightarrow 0$,
then

$\mathcal{P}_{T}(k)\approx\frac{2H^{2}}{\pi^{2}}k^{-2\varepsilon}~{}~{},~{}~{}\Pi\approx
0~{}.$ (67)

This is exactly the result of the slow-roll inflation of single scalar field
in GR.

###  V.2 The case of the irregular universe

For the case of irregular universe, since the coupling of $\zeta$ and
$h_{\mathrm{A}}$ in the action (57) makes it difficult to analyze the quantum
fluctuations, we first diagonalize the variables $\zeta$ and $h_{\mathrm{A}}$
below. Firstly, for the convenience of analysis, we introduce new variables
$\xi_{1}=(z/a)\zeta$, $\xi_{2}=(1/\sqrt{2})h_{\mathrm{L}}$ and
$\xi_{3}=(1/\sqrt{2})h_{\mathrm{R}}$, so that the action (57) can be
simplified as

$S^{(2)}=\int\mathrm{d}\eta\int\mathrm{d}^{3}k\,\frac{a^{2}}{2}\bigg{\\{}\sum_{\mathrm{s}=1}^{3}\xi^{*\prime}_{\mathrm{s}}\xi_{\mathrm{s}}-\sum_{\mathrm{s}_{1}=1}^{3}\sum_{\mathrm{s}_{2}=1}^{3}\mathbf{M}_{\mathrm{s}_{1}\mathrm{s}_{2}}\xi^{*}_{\mathrm{s}_{1}}\xi_{\mathrm{s}_{2}}\bigg{\\}}~{},~{}~{}\text{with}~{}~{}\mathbf{M}=\left(\begin{array}[]{ccc}k^{2}-\Omega&-\kappa&\kappa\\\
-\kappa&k^{2}-\sigma&0\\\ \kappa&0&k^{2}+\sigma\\\ \end{array}\right)~{},$
(68)

where $\Omega=z^{\prime\prime}/z-a^{\prime\prime}/a$,
$\sigma=-c\theta^{\prime}k$ and $\kappa=c\mathcal{H}k$ are background
quantities. Secondly, we introduce an orthogonal matrix $\mathbf{T}$ that can
diagonalize the matrix $\mathbf{M}$, and its expression is

$\mathbf{T}=\left(\begin{array}[]{c}\mathbf{t}_{1}^{\mathrm{T}}\\\
\mathbf{t}_{2}^{\mathrm{T}}\\\ \mathbf{t}_{3}^{\mathrm{T}}\\\
\end{array}\right)~{},~{}\text{with}~{}~{}\mathbf{t}_{\mathrm{s}}=\frac{-\mathrm{s}^{2}+5\mathrm{s}-5}{\sqrt{1+\frac{(\tau_{\mathrm{s}}-\sigma)^{2}}{\kappa^{2}}+\left(1-\frac{(\tau_{\mathrm{s}}-\sigma)(\tau_{\mathrm{s}}+\Omega)}{\kappa^{2}}\right)^{2}}}\left(\begin{array}[]{c}(\tau_{\mathrm{s}}-\sigma)/\kappa\\\
1-(\tau_{\mathrm{s}}-\sigma)(\tau_{\mathrm{s}}+\Omega)/\kappa^{2}\\\ 1\\\
\end{array}\right)~{},$ (69)

where the superscript $\mathrm{T}$ means transpose, and
$\\{\tau_{\mathrm{s}},~{}\mathrm{s}=1,2,3\\}$ are the solutions of the cubic
equation

$\tau^{3}+\Omega\tau^{2}-(2\kappa^{2}+\sigma^{2})\tau-\sigma^{2}\Omega=0~{}.$
(70)

The specific expressions of $\\{\tau_{\mathrm{s}},~{}\mathrm{s}=1,2,3\\}$ are
in Appendix A. Finally, we introduce new variables
$\\{q_{\mathrm{s}},~{}\mathrm{s}=1,2,3\\}$, which are defined as

$\left(\begin{array}[]{c}q_{1}\\\ q_{2}\\\ q_{3}\\\
\end{array}\right)=\mathbf{T}\left(\begin{array}[]{c}\xi_{1}\\\ \xi_{2}\\\
\xi_{3}\\\ \end{array}\right).$ (71)

Thus, the action (68) can be further simplified as

$S^{(2)}=\sum_{\mathrm{s}=1}^{3}\int\mathrm{d}\eta\int\mathrm{d}^{3}k\,\frac{a^{2}}{2}\Big{\\{}|q_{\mathrm{s}}^{\prime}|^{2}-(k^{2}+\tau_{\mathrm{s}})|q_{\mathrm{s}}|^{2}\Big{\\}}~{}.$
(72)

So far, we have simplified the action (57) with coupling between variables to
the action (72) without coupling between variables. The latter form makes it
easier to calculate the primordial fluctuations generated by inflation.

Next we consider the case of slow-roll inflation dominated by the axion-like
field $\theta$. Since in Sec. III we proved that the background equations of
the irregular flat universe are exactly the same as those in GR, the
background evolution during inflation will be exactly the same as the case of
slow-roll inflation in GR. Under the slow-roll approximation, the background
quantities $\Omega$, $\sigma$ and $\kappa$ can be approximately expressed as

$\Omega\approx\frac{3(\varepsilon+\delta)}{2\eta^{2}}~{},~{}\sigma\approx-\frac{\sqrt{2\varepsilon}ck}{\eta}~{},~{}\kappa\approx-\frac{(1+\varepsilon)ck}{\eta}~{}.$
(73)

In this section, we also assume that the coupling constant $c\sim 1$ (it can
also be seen as a requirement of naturalness), so that
$c\gg\sqrt{\varepsilon}$. Ignoring high-order small quantities such as
$\varepsilon^{2}$, $\\{\tau_{\mathrm{s}},~{}\mathrm{s}=1,2,3\\}$ in Eq. (A)
can be approximated as

$\tau_{1}\approx\frac{(2+3\varepsilon)ck}{\sqrt{2}\eta}-\frac{3(\varepsilon+\delta)}{2\eta^{2}}~{},~{}\tau_{2}\approx
0~{},~{}\tau_{3}\approx-\frac{(2+3\varepsilon)ck}{\sqrt{2}\eta}-\frac{3(\varepsilon+\delta)}{2\eta^{2}}~{}.$
(74)

If only up to the order of $\sqrt{\varepsilon}$ is retained, the orthogonal
matrix $\mathbf{T}$ can be approximated as

$\mathbf{T}\approx\left(\begin{array}[]{ccc}\frac{1}{\sqrt{2}}&\frac{1+\sqrt{\varepsilon}}{2}&-\frac{1-\sqrt{\varepsilon}}{2}\\\
-\sqrt{\varepsilon}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\\
\frac{1}{\sqrt{2}}&-\frac{1-\sqrt{\varepsilon}}{2}&\frac{1+\sqrt{\varepsilon}}{2}\\\
\end{array}\right)$ (75)

Regarding the approximate expression (75), there are two points that need
additional explanation. First, the order $\sqrt{\varepsilon}$ is the lowest
order approximation required to preserve the difference in the power spectrum
of left- and right-handed GWs. If we further ignore the contribution of
$\sqrt{\varepsilon}$ in $\mathbf{T}$, the difference in the power spectrum of
left- and right-handed GWs disappears. And if we keep the higher-order terms,
it brings only more complex but less important corrections in the power
spectrum. Second, it can be seen that the matrix $\mathbf{T}$ does not tend to
the identity matrix as $c\rightarrow 0$ in the approximate expression (75).
This is confusing because the three variables are all decoupled as
$c\rightarrow 0$ in the action (68). The reason for this confusing phenomenon
is that we have used the approximation $c\gg\sqrt{\varepsilon}$ in Eqs. (74)
and (75). If $c$ is too small, neither the Eq. (74) nor Eq. (75) hold. See
Appendix B for the approximate behavior of orthogonal matrix $\mathbf{T}$ when
$c\rightarrow 0$.

Next, by combining Eqs (72) and (74), the correlation function between
variables $q_{\mathrm{s}}$ can be obtained through the process in Appendix C:

$\langle
q_{1}^{\dagger}q_{1}\rangle\approx\frac{H^{2}}{2}\mathrm{e}^{\frac{c\pi}{\sqrt{2}}}k^{-(3+3\varepsilon+\delta)}~{},~{}\langle
q_{2}^{\dagger}q_{2}\rangle\approx\frac{H^{2}}{2}k^{-(3+2\varepsilon)}~{},~{}\langle
q_{3}^{\dagger}q_{3}\rangle\approx\frac{H^{2}}{2}\mathrm{e}^{-\frac{c\pi}{\sqrt{2}}}k^{-(3+3\varepsilon+\delta)}~{},$
(76)

and $\langle q_{\mathrm{s}_{1}}^{\dagger}q_{\mathrm{s}_{2}}\rangle=0$ when
$\mathrm{s}_{1}\neq\mathrm{s}_{2}$. Then, using the approximation techniques
in Appendix D and combining Eqs (71), (75) and (76), the correlation functions
for the variables $\zeta$ and $h_{\mathrm{A}}$ can be obtained as

$\displaystyle\langle\zeta^{\dagger}\zeta\rangle\approx\frac{1}{2\varepsilon}\cosh\left(\frac{c\pi}{\sqrt{2}}\right)H^{2}k^{n_{S}-4}~{},$
$\displaystyle\langle
h_{\mathrm{A}}^{\dagger}h_{\mathrm{A}}\rangle\approx\left[\frac{1}{2}+\frac{1}{2}\cosh\left(\frac{c\pi}{\sqrt{2}}\right)-\mathfrak{p}_{\mathrm{A}}\sqrt{\varepsilon}\sinh\left(\frac{c\pi}{\sqrt{2}}\right)\right]H^{2}k^{n_{T}-3}~{},$
$\displaystyle\langle\zeta^{\dagger}h_{\mathrm{A}}\rangle\approx-\frac{\mathfrak{p}_{\mathrm{A}}}{2\sqrt{2\varepsilon}}\sinh\left(\frac{c\pi}{\sqrt{2}}\right)H^{2}k^{-(3+3\varepsilon+\delta)}~{},$
$\displaystyle\langle
h_{\mathrm{L}}^{\dagger}h_{\mathrm{R}}\rangle\approx\frac{1}{2}\left[1-\cosh\left(\frac{c\pi}{\sqrt{2}}\right)\right]H^{2}k^{-(3+3\varepsilon+\delta)-\frac{1}{2}\,\mathrm{csch}^{2}\left(\frac{c\pi}{2\sqrt{2}}\right)(\varepsilon+\delta)}~{},$
(77)

where

$n_{S}\approx
1-(\delta+3\varepsilon)~{}~{},~{}~{}n_{T}\approx-(3\varepsilon+\delta)+\frac{1}{2}\,\mathrm{sech}^{2}\left(\frac{c\pi}{2\sqrt{2}}\right)(\varepsilon+\delta)~{}.$
(78)

It should be noted that since Eqs. (74) and (75) are approximately true only
when $c\gg\sqrt{\varepsilon}$, Eqs. (V.2) and (78) are also approximately true
only when $c\gg\sqrt{\varepsilon}$.

Through the correlation functions (V.2), the power spectrum of the scalar
perturbation $\zeta$ can be obtained as

$\mathcal{P}_{S}(k)=\frac{k^{3}}{2\pi^{2}}\langle\zeta^{\dagger}\zeta\rangle\approx\frac{H^{2}}{8\pi^{2}\varepsilon}\cosh\left(\frac{c\pi}{\sqrt{2}}\right)k^{n_{S}-1}~{}.$
(79)

The power spectrum of the left- and right-handed GWs can be obtained as

$\mathcal{P}_{\mathrm{A}}(k)=\frac{k^{3}}{\pi^{2}}\langle
h_{\mathrm{A}}^{\dagger}h_{\mathrm{A}}\rangle\approx\frac{H^{2}}{2\pi^{2}}\left[1+\cosh\left(\frac{c\pi}{\sqrt{2}}\right)-2\mathfrak{p}_{\mathrm{A}}\sqrt{\varepsilon}\sinh\left(\frac{c\pi}{\sqrt{2}}\right)\right]k^{n_{T}}.$
(80)

The power spectrum of the tensor perturbations can be obtained as

$\mathcal{P}_{T}(k)=\mathcal{P}_{\mathrm{L}}(k)+\mathcal{P}_{\mathrm{R}}(k)\approx\frac{H^{2}}{\pi^{2}}\left[1+\cosh\left(\frac{c\pi}{\sqrt{2}}\right)\right]k^{n_{T}}~{}.$
(81)

The tensor-to-scalar ratio $r$ can be obtained as

$r\equiv\frac{\mathcal{P}_{T}}{\mathcal{P}_{S}}=8\left[1+\mathrm{sech}\left(\frac{c\pi}{\sqrt{2}}\right)\right]\varepsilon~{}.$
(82)

The relative different between the power spectrum of the left- and right-
handed GWs can be obtained as

$\Pi\equiv\frac{\mathcal{P}_{\mathrm{R}}-\mathcal{P}_{\mathrm{L}}}{\mathcal{P}_{\mathrm{R}}+\mathcal{P}_{\mathrm{L}}}\approx-2\sqrt{\varepsilon}\tanh\left(\frac{c\pi}{\sqrt{2}}\right)~{}.$
(83)

Strictly speaking, since Eqs. (V.2) and (78) are only approximately true when
$c\gg\sqrt{\varepsilon}$, Eqs. (79)-(83) are also approximately true only when
$c\gg\sqrt{\varepsilon}$. But If we ignore this fact and force $c\rightarrow
0$, then

$\mathcal{P}_{S}\approx\frac{H^{2}}{8\pi^{2}\varepsilon}k^{n_{S}-1}~{},~{}\mathcal{P}_{T}\approx\frac{2H^{2}}{\pi^{2}}k^{n_{T}}~{},~{}r\approx
16\varepsilon~{},~{}\Pi\approx 0~{}.$ (84)

It can be seen that except for the spectral indices $n_{S}$ and $n_{T}$, Eq.
(84) is the result of the slow-roll inflation in GR.

Once the spectral index $n_{S},n_{T}$ and the tensor-to-scalar ratio $r$ are
measured experimentally, we can obtain the slow-roll parameters
$\varepsilon,\delta$ and the coupling constant $c$ as

$\displaystyle\varepsilon=\frac{1}{16}\Big{[}4(1-n_{S})+r+\sqrt{(r+4-4n_{S})^{2}+16r\,n_{T}}\Big{]}~{},$
(85)
$\displaystyle\delta=\frac{1}{16}\Big{[}4(1-n_{S})-3r+3\sqrt{(r+4-4n_{S})^{2}+16r\,n_{T}}\Big{]}~{},$
(86) $\displaystyle
c=\frac{\sqrt{2}}{\pi}\mathrm{arccosh}\bigg{[}\frac{\sqrt{(r+4-4n_{S})^{2}+16r\,n_{T}}-4(1-n_{S})-8n_{T}}{8(1-n_{S}+n_{T})}\bigg{]}~{}.$
(87)

In order to ensure that $c$ given by Eq. (87) is real, we need the values of
$n_{S},n_{T}$ and $r$ to fall in the blue region in FIG. 1.

Figure 1:  The blue region is the region where the coupling constant $c$ is
real. Since _Planck_ 2018 Planck:2018nkj  gives $r<0.101$ and $n_{S}\approx
0.966$, we only considered the interval $0<r/(1-n_{S})<3$.

If the values of $n_{S},n_{T}$ and $r$ fall outside the blue region in FIG. 1,
it means that $c$ given by Eq. (87) is not real, and this irregular universe
model cannot be consistent with the experiment. If the values of the
$n_{S},n_{T}$ and $r$ fall in the blue region in FIG. 1, then we can obtain
the values of $\varepsilon,\delta$ and $c$, so that we can calculate the value
of the relative different $\Pi$ between the left- and right-handed primordial
GWs through Eq. (83). On this basis, if the relative different $\Pi$ can be
measured experimentally, we can test whether this irregular universe model is
consistent with the experiment by comparing the theoretical and experimental
values of $\Pi$.

Unfortunately, primordial GWs have not yet been observed experimentally, so we
cannot know the values of tensor spectral index $n_{T}$, the tensor-to-scalar
ratio $r$ and the relative different $\Pi$ between the left- and right-handed
primordial GWs. From the _Planck_ 2018 Planck:2018nkj , we only know that the
scalar spectral index $n_{S}\approx 0.966$ and the tensor-to-scalar ratio
$r<0.101$. Therefore, all we can know is that the allowable value range of the
slow-roll parameters $\varepsilon,\delta$ is

$0<\varepsilon<\frac{0.101}{8\left[1+\mathrm{sech}\left(c\pi/\sqrt{2}\right)\right]}<0.012625~{}~{},~{}~{}\delta\approx
0.034-3\varepsilon~{}.$ (88)

The maximum value of $\varepsilon$ depends on the coupling constant $c$, but
will not exceed $0.012625$ (the upper limit of $\varepsilon$ when
$c\rightarrow\infty$). The allowable value of $\delta$ is determined by
$\varepsilon$. FIG. 2 shows the allowable value range of slow-roll parameters
$\varepsilon,\delta$ when $c=1$.

Figure 2: In the $\varepsilon$-$\delta$ plane, the blue line is the allowable
value range when $c=1$.

Although by comparing the results in subsections V.1 and V.2, we can find that
the power spectrum of the left- and right-handed GWs given by the irregular
universe is different from that of the regular universe. But this is not the
main difference between irregular and regular universes for primordial
fluctuations. For primordial fluctuations, the most important feature of the
irregular universe compared to the regular universe is that the correlation
function of scalar perturbation and tensor perturbations
$\langle\zeta^{\dagger}h_{\mathrm{A}}\rangle\neq 0$ at the linear perturbation
level. This means that there is a strong statistical correlation between
primordial scalar fluctuations and primordial tensor fluctuations generated by
slow-roll inflation in the irregular universe. The apparent reason for this
phenomenon is that the quadratic action contains the coupling of scalar
perturbations and tensor perturbations in the irregular universe, as exhibited
by the action (57). The deeper reason may be that the condition
$\mathcal{L}_{\xi_{I}}\hat{\Gamma}^{\rho}_{~{}\mu\nu}\neq 0$ destroys the
homogeneity and isotropy of the interior space, so that the scalar
fluctuations and the tensor fluctuations can interact with each other in the
irregular universe.

##  VI Conclusion

As a step towards exploring the irregular universe within the TG framework, in
this paper, we studied the irregular flat universe of the NYTG model. Firstly,
we obtained the irregular flat universe solution of the NYTG model under the
condition that only the symmetry of the metric is required. We found that the
cosmological background equations of the NYTG model are exactly the same as
those of GR in both the regular flat universe and the irregular flat universe.
Secondly, we studied the linear cosmological perturbations around the
irregular flat universes. We found a peculiar feature of the irregular flat
universe: the tensor and scalar perturbations are coupled together at the
linear perturbation level. We speculate that this peculiar feature is caused
by the fact that the interior space does not satisfy the homogeneity and
isotropy in the irregular universe. Finally, we applied the NYTG model to the
early universe and studied the primordial perturbations generated by slow-roll
inflation in the regular and irregular flat universes. We found that the left-
and right-handed primordial GWs are different in both the regular flat
universe and the irregular flat universe. We also found that there is a strong
statistical correlation between the primordial scalar and tensor perturbations
generated by slow-roll inflation in the case of irregular universe, this is a
direct consequence of the direct coupling between the scalar and tensor
perturbations at linear order.

Acknowledgement: This work is supported in part by National Key R&D Program of
China Grant No. 2021YFC2203102, and by NSFC under Grant No. 12075231 and
12047502.

##  Appendix A Solutions of the cubic equation

Consider a cubic equation with respect to the variable $\tau$ as

$\mathfrak{a}\tau^{3}+\mathfrak{b}\tau^{2}+\mathfrak{c}\tau+\mathfrak{d}=0~{},$
(89)

where $\mathfrak{a}$, $\mathfrak{b}$, $\mathfrak{c}$ and $\mathfrak{d}$ are
real coefficients. In order to express the solution of Eq. (89) conveniently,
we introduce the following parameters

$\mathfrak{A}=\mathfrak{b}^{2}-3\mathfrak{a}\mathfrak{c}~{},~{}\mathfrak{B}=\mathfrak{b}\mathfrak{c}-9\mathfrak{a}\mathfrak{d}~{},~{}\mathfrak{C}=\mathfrak{c}^{2}-3\mathfrak{b}\mathfrak{d}~{},~{}\Delta=\mathfrak{B}^{2}-4\mathfrak{A}\mathfrak{C}~{},~{}\Theta=\frac{1}{3}\arccos\left(\frac{2\mathfrak{A}\mathfrak{b}-3\mathfrak{B}\mathfrak{a}}{2\mathfrak{A}^{3/2}}\right)~{}.$
(90)

When $\Delta<0$, Eq. (89) has three real solutions, which are

$\displaystyle\tau_{1}=-\frac{1}{3\mathfrak{a}}\left(\mathfrak{b}+2\sqrt{\mathfrak{A}}\cos\Theta\right)~{},$
$\displaystyle\tau_{2}=\frac{1}{3\mathfrak{a}}\left[-\mathfrak{b}+\sqrt{\mathfrak{A}}\left(\cos\Theta-\sqrt{3}\sin\Theta\right)\right]~{},$
$\displaystyle\tau_{3}=\frac{1}{3\mathfrak{a}}\left[-\mathfrak{b}+\sqrt{\mathfrak{A}}\left(\cos\Theta+\sqrt{3}\sin\Theta\right)\right]~{},$
(91)

The Eq. (70) in the main text is the result of taking $\mathfrak{a}=1$,
$\mathfrak{b}=\Omega$, $\mathfrak{c}=-(2\kappa^{2}+\sigma^{2})$ and
$\mathfrak{d}=-\sigma^{2}\Omega$ in Eq. (89). In this case, there are always
$\mathfrak{A}\geq 0$ and $\Delta\leq 0$, where the equal sign holds if and
only if $\kappa=\sigma=\Omega=0$. And when $\kappa=\sigma=\Omega=0$, obviously
the three solutions of Eq. (70) are $\tau_{1}=\tau_{2}=\tau_{3}=0$, and the
orthogonal matrix $\mathbf{T}$ is the identity matrix.

##  Appendix B The orthogonal matrix $\mathbf{T}$ when $c\rightarrow 0$

In this appendix, we discuss the approximate behavior of the orthogonal matrix
$\mathbf{T}$ in Eq. (69) as $c\rightarrow 0$ in a more general background (not
only during inflation). Since $\sigma\propto c$ and $\kappa\propto c$, then

$\frac{\kappa}{\Omega}\propto c~{},~{}\frac{\sigma}{\Omega}\propto
c~{},~{}\frac{\kappa^{2}}{2\sigma\Omega}\propto c~{}.$ (92)

When $c$ is much smaller than any other background quantities such as
$\sqrt{\varepsilon}$, $\dot{\theta}$ and $\mathcal{H}^{-1}$, ignoring the
quadratic and higher terms of $c$, the solutions of Eq. (70) can be
approximately expressed as

$\tau_{1}\approx\Omega~{},~{}\tau_{2}\approx\sigma~{},~{}\tau_{3}\approx\sigma~{}.$
(93)

So the orthogonal matrix $\mathbf{T}$ in Eq. (69) can be approximately
expressed as

$\mathbf{T}=\left(\begin{array}[]{ccc}1&\frac{\kappa}{\Omega}&-\frac{\kappa}{\Omega}\\\
-\frac{\kappa}{\Omega}&1&\frac{\kappa^{2}}{2\sigma\Omega}\\\
\frac{\kappa}{\Omega}&-\frac{\kappa^{2}}{2\sigma\Omega}&1\\\
\end{array}\right)\xrightarrow{\text{when}~{}c\rightarrow
0}\left(\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\ 0&0&1\\\ \end{array}\right)$
(94)

It can be easily seen from Eqs. (92) and (94) that when $c\rightarrow 0$, the
orthogonal matrix $\mathbf{T}$ does tend to the identity matrix. This is
consistent with the fact that all variables in the action (68) tend to be
decoupled when $c\rightarrow 0$.

##  Appendix C Correlation function generated by inflation

The purpose of this appendix is to show how to calculate the correlation
function generated by inflation. Consider a univariate system whose effective
action during inflation is

$S=\frac{1}{2}\int\mathrm{d}\eta\,\mathrm{d}^{3}k\,a^{2}\left[|q_{\vec{k}}^{\prime}|^{2}-\left(k^{2}-\frac{2\mathfrak{a}k}{\eta}-\frac{3\mathfrak{b}}{\eta^{2}}\right)|q_{\vec{k}}|^{2}\right]~{},$
(95)

where $\mathfrak{a}$ and $\mathfrak{b}$ are real parameters, and
$\mathfrak{b}$ has the same order of magnitude as the slow-roll parameter
$\varepsilon$. Here $q(\eta,\vec{x})$ is the variable and we have changed to
the Fourier space $q_{\vec{k}}(\eta)$. After quantization, the variable
$q_{\vec{k}}(\eta)$ can be expanded as

$q_{\vec{k}}(\eta)=\frac{1}{a(\eta)}\left(v_{k}(\eta)\hat{a}_{\vec{k}}+v_{k}^{*}(\eta)\hat{a}_{\vec{k}}^{\dagger}\right)~{},$
(96)

where $\hat{a}_{\vec{k}}^{\dagger}$ and $\hat{a}_{\vec{k}}$ are the generation
and annihilation operators that satisfy the following commutation relations

$[\hat{a}_{\vec{k}}~{}\hat{a}_{\vec{k}^{\prime}}^{\dagger}]=\updelta^{(3)}(\vec{k}-\vec{k}^{\prime})~{}~{},~{}~{}[\hat{a}_{\vec{k}}~{}\hat{a}_{\vec{k}^{\prime}}]=[\hat{a}_{\vec{k}}^{\dagger}~{}\hat{a}_{\vec{k}^{\prime}}^{\dagger}]=0~{},$
(97)

and $v_{k}(\eta)$ satisfies the following equation

$v_{k}^{\prime\prime}+\left(k^{2}-\frac{2\mathfrak{a}k}{\eta}-\frac{\mu^{2}-1/4}{\eta^{2}}\right)v_{k}=0~{},$
(98)

where $\mu\approx 3/2+\varepsilon+\mathfrak{b}$. Note that in Eq. (98), we
used the approximation
$a^{\prime\prime}/a\approx[(3/2+\varepsilon)^{2}-1/4]/\eta$, and we ignored
the higher-order terms of $\varepsilon$ and $\mathfrak{b}$. Next we choose the
Bunch-Davies vacuum at $\eta\rightarrow-\infty$, that is,

$\lim_{\eta\rightarrow-\infty}v_{k}=\frac{1}{\sqrt{2k}}\mathrm{e}^{-\mathrm{i}k\eta}~{}.$
(99)

Under this condition, the solution for Eq. (98) is (for more detail, see
Satoh:2010ep )

$v_{k}(\eta)=\mathrm{e}^{-\mathrm{i}k\eta}(-2k\eta)^{\mu}(-\eta)^{\frac{1}{2}}\mathrm{e}^{-\mathrm{i}\pi(\frac{1}{4}+\frac{\mu}{2})}\mathbf{U}\left(1/2+\mu-\mathrm{i}\mathfrak{a},1+2\mu;~{}2\mathrm{i}k\eta\right)\mathrm{e}^{-\frac{\mathfrak{a}\pi}{2}}~{},$
(100)

where $\mathbf{U}(c_{1},c_{2};\,\mathrm{z})$ is the confluent hypergeometric
function. The $|v_{k}|$ has the following asymptotic form when
$k\eta\rightarrow 0^{-}$ (super-horizon scale)

$|v_{k}|\approx
2^{\mu-1}\pi^{-\frac{1}{2}}\Gamma(\mu)k^{-\mu}(-\eta)^{\frac{1}{2}-\mu}\mathrm{e}^{-\frac{\mathfrak{a}\pi}{2}}\approx
2^{-\frac{1}{2}}\mathrm{e}^{-\frac{\mathfrak{a}\pi}{2}}aHk^{-\mu}$ (101)

where $\Gamma(\mathrm{z})$ is the Gamma function. In the last approximately
equal sign in Eq. (101), we used the approximations $\mu\approx 3/2$ and
$(-\eta)^{-1}\approx aH$. Combining Eqs. (96), (97) and (101), we can obtain
the correlation function on the super-horizon scale as

$\langle
0|q^{\dagger}_{\vec{k}}q_{\vec{k}^{\prime}}|0\rangle\approx\frac{H^{2}}{2}\mathrm{e}^{-\mathfrak{a}\pi}k^{-(3+2\varepsilon+2\mathfrak{b})}\updelta^{(3)}(\vec{k}+\vec{k}^{\prime})~{}.$
(102)

where $|0\rangle$ is the vacuum state, which satisfies
$\hat{a}_{\vec{k}}|0\rangle=0$. For the sake of convenience, we can omit the
subscript $\vec{k}$ and throw away the annoying delta function
$\updelta^{(3)}(\vec{k}+\vec{k}^{\prime})$, so that the correlation function
(102) can be abbreviated as

$\langle
q^{\dagger}q\rangle\approx\frac{H^{2}}{2}\mathrm{e}^{-\mathfrak{a}\pi}k^{-(3+2\varepsilon+2\mathfrak{b})}~{}.$
(103)

##  Appendix D Summation of nearly scale-invariant functions

Consider there are $N$ nearly scale-invariant functions
$\\{\mathrm{f}_{\mathfrak{i}}(k)=C_{\mathfrak{i}}k^{n_{\mathfrak{i}}},~{}\mathfrak{i}=1,2,...,N\\}$,
where $|n_{\mathfrak{i}}|\ll 1$. Then the sum of these functions should also
be a nearly scale-invariant function, so it can be approximated as

$\mathrm{f}(k)=\sum_{\mathfrak{i}=1}^{N}\mathrm{f}_{\mathfrak{i}}(k)=\sum_{\mathfrak{i}=1}^{N}C_{\mathfrak{i}}k^{n_{\mathfrak{i}}}\approx
Ck^{n}~{},~{}\text{with}~{}|n|\ll 1~{}.$ (104)

Next we need to find the coefficient $C$ and the exponent $n$ in Eq. (104).
Since $n_{\mathfrak{i}}\approx 0$ and $n\approx 0$, we can approximately let
$n_{\mathfrak{i}}=n=0$ in Eq. (104), so that Eq. (104) becomes

$C\approx\sum_{\mathfrak{i}=1}^{N}C_{\mathfrak{i}}~{}.$ (105)

Next, let Eq. (104) take the derivative of $k$ and then let
$n_{\mathfrak{i}}=n=0$ on the exponent of $k$. Then the approximate expression
of $n$ can be obtained as

$n\approx\frac{1}{C}\sum_{\mathfrak{i}=1}^{N}C_{\mathfrak{i}}n_{\mathfrak{i}}~{}.$
(106)

## References

  * (1) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, no.6, 061102 (2016) doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * (2) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 119, no.16, 161101 (2017) doi:10.1103/PhysRevLett.119.161101 [arXiv:1710.05832 [gr-qc]].
  * (3) H. Li, S. Y. Li, Y. Liu, Y. P. Li, Y. Cai, M. Li, G. B. Zhao, C. Z. Liu, Z. W. Li and H. Xu, et al. “Probing Primordial Gravitational Waves: Ali CMB Polarization Telescope,” Natl. Sci. Rev. 6, no.1, 145-154 (2019) doi:10.1093/nsr/nwy019 [arXiv:1710.03047 [astro-ph.CO]].
  * (4) K. Abazajian et al. [CMB-S4], Astrophys. J. 926, no.1, 54 (2022) doi:10.3847/1538-4357/ac1596 [arXiv:2008.12619 [astro-ph.CO]].
  * (5) R. Jackiw and S. Y. Pi, Phys. Rev. D 68, 104012 (2003) doi:10.1103/PhysRevD.68.104012 [arXiv:gr-qc/0308071 [gr-qc]].
  * (6) S. Alexander and N. Yunes, Phys. Rept. 480, 1-55 (2009) doi:10.1016/j.physrep.2009.07.002 [arXiv:0907.2562 [hep-th]].
  * (7) S. Dyda, E. E. Flanagan and M. Kamionkowski, Phys. Rev. D 86, 124031 (2012) doi:10.1103/PhysRevD.86.124031 [arXiv:1208.4871 [gr-qc]].
  * (8) M. Crisostomi, K. Noui, C. Charmousis and D. Langlois, Phys. Rev. D 97, no.4, 044034 (2018) doi:10.1103/PhysRevD.97.044034 [arXiv:1710.04531 [hep-th]].
  * (9) X. Gao and X. Y. Hong, Phys. Rev. D 101, no.6, 064057 (2020) doi:10.1103/PhysRevD.101.064057 [arXiv:1906.07131 [gr-qc]].
  * (10) W. Zhao, T. Zhu, J. Qiao and A. Wang, “Waveform of gravitational waves in the general parity-violating gravities,” Phys. Rev. D 101, no.2, 024002 (2020) doi:10.1103/PhysRevD.101.024002 [arXiv:1909.10887 [gr-qc]].
  * (11) N. Bartolo, L. Caloni, G. Orlando and A. Ricciardone, JCAP 03, 073 (2021) doi:10.1088/1475-7516/2021/03/073 [arXiv:2008.01715 [astro-ph.CO]].
  * (12) M. Li, H. Rao and D. Zhao, JCAP 11, 023 (2020) doi:10.1088/1475-7516/2020/11/023 [arXiv:2007.08038 [gr-qc]].
  * (13) M. Li, H. Rao and Y. Tong, Phys. Rev. D 104, no.8, 084077 (2021) doi:10.1103/PhysRevD.104.084077 [arXiv:2104.05917 [gr-qc]].
  * (14) R. Aldrovandi and J. G. Pereira, Teleparallel Gravity, Vol. 173. Springer, 23 Dordrecht, (2013).
  * (15) S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. L. Said, J. Mifsud and E. Di Valentino, [arXiv:2106.13793 [gr-qc]].
  * (16) J. W. Maluf, Annalen Phys. 525, 339-357 (2013) doi:10.1002/andp.201200272 [arXiv:1303.3897 [gr-qc]].
  * (17) H. T. Nieh and M. L. Yan, J. Math. Phys. 23, 373 (1982) doi:10.1063/1.525379
  * (18) H. Rao, Phys. Rev. D 104, no.12, 124084 (2021) doi:10.1103/PhysRevD.104.124084 [arXiv:2107.08597 [gr-qc]].
  * (19) J. Qiao, T. Zhu, G. Li and W. Zhao, JCAP 04, no.04, 054 (2022) doi:10.1088/1475-7516/2022/04/054 [arXiv:2110.09033 [gr-qc]].
  * (20) Q. Wu, T. Zhu, R. Niu, W. Zhao and A. Wang, Phys. Rev. D 105, no.2, 024035 (2022) doi:10.1103/PhysRevD.105.024035 [arXiv:2110.13870 [gr-qc]].
  * (21) R. G. Cai, C. Fu and W. W. Yu, Phys. Rev. D 105, no.10, 103520 (2022) doi:10.1103/PhysRevD.105.103520 [arXiv:2112.04794 [astro-ph.CO]].
  * (22) M. Li and D. Zhao, Phys. Lett. B 827 (2022), 136968 doi:10.1016/j.physletb.2022.136968 [arXiv:2108.01337 [gr-qc]].
  * (23) C. Gong, T. Zhu, R. Niu, Q. Wu, J. L. Cui, X. Zhang, W. Zhao and A. Wang, Phys. Rev. D 105 (2022) no.4, 044034 doi:10.1103/PhysRevD.105.044034 [arXiv:2112.06446 [gr-qc]].
  * (24) M. Hohmann and C. Pfeifer, Phys. Lett. B 834 (2022), 137437 doi:10.1016/j.physletb.2022.137437 [arXiv:2203.01856 [gr-qc]].
  * (25) X. Tong and Z. Z. Xianyu, JHEP 10 (2022), 194 doi:10.1007/JHEP10(2022)194 [arXiv:2203.06349 [hep-ph]].
  * (26) M. Li, Y. Tong and D. Zhao, Phys. Rev. D 105 (2022) no.10, 104002 doi:10.1103/PhysRevD.105.104002 [arXiv:2203.06912 [gr-qc]].
  * (27) F. Zhang, J. X. Feng and X. Gao, JCAP 10 (2022), 054 doi:10.1088/1475-7516/2022/10/054 [arXiv:2205.12045 [gr-qc]].
  * (28) T. Zhu, W. Zhao and A. Wang, Phys. Rev. D 107, no.2, 024031 (2023) doi:10.1103/PhysRevD.107.024031 [arXiv:2210.05259 [gr-qc]].
  * (29) T. Zhu, W. Zhao and A. Wang, Phys. Rev. D 107, no.4, 044051 (2023) doi:10.1103/PhysRevD.107.044051 [arXiv:2211.04711 [gr-qc]].
  * (30) A. A. A. Filho, J. R. Nascimento, A. Y. Petrov and P. J. Porfírio, [arXiv:2211.11821 [gr-qc]].
  * (31) J. Qiao, Z. Li, T. Zhu, R. Ji, G. Li and W. Zhao, Front. Astron. Space Sci. 9, 1109086 (2023) doi:10.3389/fspas.2022.1109086 [arXiv:2211.16825 [gr-qc]].
  * (32) Y. Cai, Phys. Rev. D 107, no.6, 063512 (2023) doi:10.1103/PhysRevD.107.063512 [arXiv:2212.10893 [gr-qc]].
  * (33) Z. Chen, Y. Yu and X. Gao, [arXiv:2212.14362 [gr-qc]].
  * (34) M. Hohmann, L. Järv, M. Krššák and C. Pfeifer, Phys. Rev. D 100, no.8, 084002 (2019) doi:10.1103/PhysRevD.100.084002 [arXiv:1901.05472 [gr-qc]].
  * (35) M. Hohmann, Int. J. Geom. Meth. Mod. Phys. 18, no.supp01, 2140005 (2021) doi:10.1142/S0219887821400053 [arXiv:2008.12186 [gr-qc]].
  * (36) A. A. Coley, R. J. v. Hoogen and D. D. McNutt, [arXiv:2205.10719 [gr-qc]].
  * (37) R. Myrzakulov, Eur. Phys. J. C 71, 1752 (2011) doi:10.1140/epjc/s10052-011-1752-9 [arXiv:1006.1120 [gr-qc]].
  * (38) Y. F. Cai, S. H. Chen, J. B. Dent, S. Dutta and E. N. Saridakis, Class. Quant. Grav. 28, 215011 (2011) doi:10.1088/0264-9381/28/21/215011 [arXiv:1104.4349 [astro-ph.CO]].
  * (39) Y. F. Cai, S. Capozziello, M. De Laurentis and E. N. Saridakis, Rept. Prog. Phys. 79, no.10, 106901 (2016) doi:10.1088/0034-4885/79/10/106901 [arXiv:1511.07586 [gr-qc]].
  * (40) C. Bejarano, R. Ferraro and M. J. Guzmán, Eur. Phys. J. C 77, no.12, 825 (2017) doi:10.1140/epjc/s10052-017-5394-4 [arXiv:1707.06637 [gr-qc]].
  * (41) C. Bejarano, R. Ferraro, F. Fiorini and M. J. Guzmán, Universe 5, 158 (2019) doi:10.3390/universe5060158 [arXiv:1905.09913 [gr-qc]].
  * (42) J. Beltrán Jiménez, A. Golovnev, T. Koivisto and H. Veermäe, Phys. Rev. D 103, no.2, 024054 (2021) doi:10.1103/PhysRevD.103.024054 [arXiv:2004.07536 [gr-qc]].
  * (43) A. Golovnev and M. J. Guzman, Phys. Rev. D 103, no.4, 044009 (2021) doi:10.1103/PhysRevD.103.044009 [arXiv:2012.00696 [gr-qc]].
  * (44) R. P. Woodard, Lect. Notes Phys. 720, 403-433 (2007) doi:10.1007/978-3-540-71013-4_14 [arXiv:astro-ph/0601672 [astro-ph]].
  * (45) M. Hohmann and C. Pfeifer, Eur. Phys. J. C 81, no.4, 376 (2021) doi:10.1140/epjc/s10052-021-09165-x [arXiv:2012.14423 [gr-qc]].
  * (46) M. Li, Z. Li and H. Rao, Phys. Lett. B 834, 137395 (2022) doi:10.1016/j.physletb.2022.137395 [arXiv:2201.02357 [gr-qc]].
  * (47) K. Izumi and Y. C. Ong, JCAP 06, 029 (2013) doi:10.1088/1475-7516/2013/06/029 [arXiv:1212.5774 [gr-qc]].
  * (48) Y. C. Ong, K. Izumi, J. M. Nester and P. Chen, Phys. Rev. D 88, 024019 (2013) doi:10.1103/PhysRevD.88.024019 [arXiv:1303.0993 [gr-qc]].
  * (49) A. De Felice, A. E. Gumrukcuoglu and S. Mukohyama, Phys. Rev. Lett. 109, 171101 (2012) doi:10.1103/PhysRevLett.109.171101 [arXiv:1206.2080 [hep-th]].
  * (50) A. Delhom, A. Jiménez-Cano and F. J. Maldonado Torralba, [arXiv:2207.13431 [gr-qc]].
  * (51) Y. Akrami et al. [Planck], Astron. Astrophys. 641 (2020), A10 doi:10.1051/0004-6361/201833887 [arXiv:1807.06211 [astro-ph.CO]].
  * (52) D. Baumann, doi:10.1142/9789814327183_0010 [arXiv:0907.5424 [hep-th]].
  * (53) Y. Wang, Commun. Theor. Phys. 62, 109-166 (2014) doi:10.1088/0253-6102/62/1/19 [arXiv:1303.1523 [hep-th]].
  * (54) N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A1 doi:10.1051/0004-6361/201833880 [arXiv:1807.06205 [astro-ph.CO]].
  * (55) M. Satoh, JCAP 11, 024 (2010) doi:10.1088/1475-7516/2010/11/024 [arXiv:1008.2724 [astro-ph.CO]].