# Multi-agent reinforcement learning for wall modeling in LES of flow over
periodic hills

D. Zhou111Graduate Aerospace Laboratories, California Institute of Technology
M. P. Whitmore  K. P. Griffin  H. J. Bae$\dagger$

We develop a wall model for large-eddy simulation (LES) that takes into
account various pressure-gradient effects using multi-agent reinforcement
learning (MARL). The model is trained using low-Reynolds-number flow over
periodic hills with agents distributed on the wall along the computational
grid points. The model utilizes a wall eddy-viscosity formulation as the
boundary condition, which is shown to provide better predictions of the mean
velocity field, rather than the typical wall-shear stress formulation. Each
agent receives states based on local instantaneous flow quantities at an off-
wall location, computes a reward based on the estimated wall-shear stress, and
provides an action to update the wall eddy viscosity at each time step. The
trained wall model is validated in wall-modeled LES (WMLES) of flow over
periodic hills at higher Reynolds numbers, and the results show the
effectiveness of the model on flow with pressure gradients. The analysis of
the trained model indicates that the model is capable of distinguishing
between the various pressure gradient regimes present in the flow.

##  1 Introduction

LES is an essential technology for the simulation of turbulent flows. The
basic premise of LES is that energy-containing and dynamically important
eddies must be resolved everywhere in the domain. This requirement is hard to
meet in the near-wall region, as the stress-producing eddies become
progressively smaller. The cost involved in resolving the near-wall region has
motivated the development of WMLES, which utilizes a wall model to account for
the effect of the energetic near-wall eddies. Because of such characteristics,
wall modeling has been accepted as the next step to enable the increased use
of high-fidelity LES in realistic engineering and geophysical applications.

Recently, Bae & Koumoutsakos (2022) demonstrated the efficacy of MARL as a
development tool for wall models in zero-pressure-gradient flat-plate flow. In
that work, a series of learning agents are distributed along the computational
grid points, with each agent receiving local states and rewards then providing
local actions at each time step. The MARL model is able to achieve these
results by training in situ on moderate-Reynolds-number flows with a reward
function only based on the mean wall-shear stress. Therefore, the trained
models do not require any higher-fidelity simulation data and are data
efficient, unlike supervised learning. Furthermore, MARL can develop novel
models that are optimized to accurately reproduce the quantities of interest
by discovering dominant patterns in flow physics rather than using _a-priori_
parameterization.

In the present study, we extend the methodology of Bae & Koumoutsakos (2022)
for pressure-gradient flows and apply it to the flow over periodic hills,
training on low-Reynolds-number cases and testing on higher-Reynolds-number
flows. Our objective is to develop a wall model for LES based on MARL that is
robust to pressure-gradient effects in a data-efficient way. The remainder of
this report is organized as follows. In Section 2, the setup of the flow
simulation and reinforcement learning are introduced. In addition, an
investigation of two candidate wall boundary conditions is presented. In
Section 3, tests of the wall model are presented for flow over periodic hills
at various Reynolds numbers. Finally, conclusions are drawn in Section 4.

##  2 Methodology

###  2.1 Flow simulation

For the flow solver, we utilize a finite volume, unstructured-mesh
incompressible LES code (You et al., 2008). The subgrid-scale (SGS) stress is
modeled using the dynamic Smagorinsky model (Germano et al., 1991; Lilly,
1992).

For training the wall model, a flow configuration that (i) has widely
available wall-shear stress profiles for several Reynolds numbers and (ii)
does not require tuning of the inlet profile or other boundary conditions is
chosen. The flow over periodically arranged hills in a channel (Mellen et al.,
2000) has well-defined boundary conditions, can be computed at affordable
costs, and constitutes a canonical smooth-body separation and reattachment,
which is representative of many aerodynamic applications. Furthermore, the
flow does not require reconfiguration of the inlet boundary condition for
different grid resolutions and wall models, which may be necessary for
nonperiodic flows. This configuration has become a popular benchmark test case
for validating computational fluid dynamics codes, since numerous experimental
and high-fidelity numerical references exist (e.g., Krank et al., 2018;
Gloerfelt & Cinnella, 2019) that provide extensive data on a wide range of
Reynolds numbers covering $700\leq Re_{H}\leq 37000$, where
$Re_{H}=U_{b}H/\nu$, $H$ is the hill height, and $U_{b}$ is the bulk velocity
at the top of the hill. The periodic-hill configuration has the dimensions of
$9H\times 3.035H\times 4.5H$ in the streamwise ($x$), vertical ($y$), and
spanwise ($z$) directions, respectively. In the simulations of the present
study, periodic boundary conditions are applied in streamwise and spanwise
boundaries, and the equilibrium stress-balance wall model (EQWM; Kawai &
Larsson, 2012) is employed at the top wall. To maintain constant bulk velocity
in time, the flow is driven by time-varying body force (Balakumar et al.,
2014). A set of three meshes with increasing resolution are used in the
present study. More details of the grids utilized are listed in Table 1. The
grid is evenly spaced in the $z$ direction and is slightly non-uniform in both
the $x$ and $y$ directions. A CFL number of 1 is used for all simulations.

Case | Mesh size ($N_{x}\times N_{y}\times N_{z}$) | $x_{sep}/H$ | ${x_{rea}/H}$
---|---|---|---
RLWM | $128\times 64\times 64$ | -0.03 | 3.91
RLWM, coarse mesh | $64\times 32\times 32$ | 0.36 | 4.16
RLWM, fine mesh | $128\times 128\times 72$ | 0.21 | 3.85
EQWM | $128\times 64\times 64$ | 0.57 | 3.05
DNS (Krank et al., 2018) | $896\times 448\times 448$ | 0.20 | 4.51
WRLES (Gloerfelt & Cinnella, 2019) | $512\times 256\times 256$ | -0.11 | 4.31
Table 1: Simulation cases in comparison to reference data, including mesh
size, mean separation location $x_{sep}$, and mean reattachment location
$x_{rea}$ at $Re_{H}=10595$.

###  2.2 Reinforcement learning architecture

The MARL architecture of wall-model training in the study is based on Bae &
Koumoutsakos (2022). During model training, the agents distributed above the
wall receive states based on local instantaneous flow information and a reward
based on the estimated wall-shear stress, then provide local actions to update
the wall boundary condition at each time step. The agents infer a single
optimized policy through their repeated interactions with the flow field to
maximize its cumulative long-term reward. To utilize MARL as a model
development tool for wall models, an open-source reinforcement learning (RL)
toolbox smarties (Novati & Koumoutsakos, 2019) is coupled with the flow
solver. The coupling is validated using the same training configuration as the
preliminary study of Bae & Koumoutsakos (2022).

###  2.3 Boundary conditions for wall modeling

Previous studies show that how the wall-shear stress is applied affects the
mean quantities (Balakumar et al., 2014). We test two formulations of the
boundary condition, namely the wall-shear stress and the wall eddy-viscosity
formulations, to determine the appropriate action for the wall model. Two
simulations for the flow over periodic hills at $Re_{H}=10595$ are conducted
using the same mesh ($N_{x}\times N_{y}\times N_{z}=128\times 64\times 64$)
but with the two different boundary conditions. For the wall-shear stress
formulation, the wall-shear stress $\tau_{w}$ is set to the mean wall-shear
stress from WRLES (Gloerfelt & Cinnella, 2019). For the wall eddy-viscosity
formulation, the wall-shear stress is calculated as
$\tau_{w}=(\nu+\nu_{t,w})(\partial{u_{s}}/\partial{n})_{w}$, where ${u_{s}}$
is the wall-parallel velocity, $n$ denotes the wall-normal direction pointing
toward the interior, and the subscript $w$ denotes quantities evaluated at the
wall. The wall eddy viscosity $\nu_{t,w}(x)$ is obtained from the interpolated
mean-flow field and mean skin friction of WRLES (Gloerfelt & Cinnella, 2019).
Similar boundary conditions for wall modeling have been used in the past (Bae
& Lozano-Durán, 2021).

The mean skin friction coefficients of the two simulations are compared in
Figure 1(a). The skin friction coefficient is defined as
$C_{f}=\langle\tau_{w}\rangle/(0.5\rho U_{b}^{2})$, where $\rho$ is the
density and $\langle\cdot\rangle$ denotes the averaging in both time and
spanwise direction. The $C_{f}$ of the wall-shear stress formulation is equal
to that from WRLES (Gloerfelt & Cinnella, 2019) by design, whereas the wall
eddy-viscosity formulation deviates slightly, with some differences observed
near the hill peak and on the leeward side of the hill. Despite the similar
predictions of $C_{f}$, the velocity fields of the two simulations are
different [Figure 1(b)]. The size of the separation bubble behind the hill is
largely underpredicted in the wall-shear stress formulation, whereas the
prediction from the wall eddy-viscosity formulation agrees reasonably well
with the WRLES results (see Figure 2).

Figure 1: (a) Mean skin friction coefficient along the bottom wall at
$Re_{H}=10595$;  , LES using wall eddy-viscosity formulation;  , LES using
wall-shear stress formulation. (b) Contours of the mean streamwise velocity
and the streamlines at $Re_{H}=10595$: (left) LES using wall eddy-viscosity
formulation and (right) LES using wall-shear stress formulation.

The underprediction of the separation bubble in the wall-shear stress
formulation can be explained by the directional inconsistency of the $u_{s}$
at the wall-adjacent cell and $\tau_{w}$. When the direction of the $\tau_{w}$
is different from $u_{s}$ at the wall-adjacent cell, the forces acting on the
control volume of the wall-adjacent cell will push the $u_{s}$ in the
incorrect direction in a positive feedback loop, leading to an inaccurate
prediction of the separation bubble. We find that the wall eddy-viscosity
formulation is a more robust method for WMLES of such flows and will utilize
it for the remainder of the paper.

###  2.4 Training of the wall model

The model training is carried out based on the LES of periodic-hill channel
flow at $Re_{H}=10595$ using the baseline mesh ($N_{x}\times N_{y}\times
N_{z}=128\times 64\times 64$). A total of 512 agents are uniformly distributed
along the bottom wall, and the wall-normal locations of the agents $h_{m}$ are
randomly selected between $0.01H$ and $0.09H$ at each agent location. For each
agent, we set three local flow quantities as states: (i) the instantaneous
wall-parallel velocity $u_{s}$; (ii) the wall-normal location of the agent
$y_{n}=h_{m}$; and (iii) the turbulence strain rate
$S_{12}=(\partial{u_{s}}/\partial{n}+\partial{u_{n}}/\partial{s})/2$, where
$s$ denotes the wall-parallel direction pointing toward the positive
$x$-direction. Moreover, to increase the applicability of the wall model for a
wide range of flow parameters, the states are nondimensionalized using
kinematic viscosity $\nu$ and the instantaneous composite friction velocity
$u_{\tau p}=({u_{\tau}^{2}+u_{p}^{2}})^{1/2}$, where
$u_{p}=|(\nu/\rho)(\partial{p_{w}}/\partial{s})|^{1/3}$, $p_{w}$ is the
pressure on the bottom wall, and $u_{\tau}$ is the friction velocity based on
the instantaneous modeled wall-shear stress. Nondimensionalized quantities are
denoted by superscript $*$. Each agent acts to adjust the local wall eddy
viscosity $\nu_{t,w}$ at each time step through a multiplication factor
$\nu_{t,w}(t_{i+1})=a\nu_{t,w}(t_{i})$, where $a\in\left[1-\alpha\Delta
TU_{b}/\Delta x,1+\alpha\Delta TU_{b}/\Delta x\right]$, $\Delta T$ is the
action time step, $\Delta x$ is the grid size in $x$ direction, and $\alpha$
is a constant that is selected to be 0.001. The local wall-shear stress can be
calculated by the formula
$\tau_{w}=\nu(1+\nu_{t,w}^{*})(\partial{u_{s}}/\partial{n})_{w}$. The reward
$r$ is calculated based on
$r(t_{i})=[|\tau_{w}^{ref}-\tau_{w}(t_{i})|-|\tau_{w}^{ref}-\tau_{w}(t_{i-1})|]/\tau_{w,rms}^{ref}$
at each location, where $\tau_{w}^{ref}$ and $\tau_{w,rms}^{ref}$ are the mean
and root-mean-square wall-shear stress from reference simulation,
respectively. The reward is proportional to the improvement in the modeled
wall-shear stress compared to the one obtained in the previous time step, and
an extra reward of 0.1 is added when the modeled $\tau_{w}$ is within $10\%$
of the reference value.

During model training, each iteration is initialized with the normalized wall
eddy viscosity $\nu_{t,w}^{*}$ that is randomly selected from $(0,6]$. To
generate the initial condition for training, the simulation is started from a
flow field generated by the EQWM and run with the given initial $\nu_{t,w}$
for 20 flow-through times (FTTs) to remove numerical artifacts. Each iteration
of the model training simulation is conducted for 5 FTTs, and to increase the
data efficiency, the $\nu_{t,w}$ is updated every 100 time steps.

The model policy is learned based on a neural network with two hidden layers
of 128 units each with a soft sign activation function. The parameters of the
neural network are identical to those used in Bae & Koumoutsakos (2022). The
model training is advanced for 1.6 million policy gradient steps.

##  3 Results

###  3.1 Testing for flow over periodic hills at
$\textnormal{Re}_{\textnormal{H}}=10595$

To evaluate the performance of the trained wall model (RLWM), three
simulations for the periodic-hill channel flow at $Re_{H}=10595$ are carried
out using meshes with different resolutions. The details of the simulation
cases are listed in Table 1. The number of agents on the bottom wall is
consistent with the number of mesh cells on the wall. The wall-normal matching
location of the agents was chosen to be within the first off-wall cell. The
wall eddy viscosity $\nu_{t,w}$ is updated based on the model action at every
time step. All simulations are run for about 60 FTTs after transients. The
flow statistics of all simulations are averaged over spanwise direction and
time. For comparison, an LES using EQWM is conducted and results from two
high-fidelity reference databases for this flow (Krank et al., 2018; Gloerfelt
& Cinnella, 2019) are included.

Figure 2: Contours of the velocity in $x$ direction and the streamlines at
$Re_{H}=10595$: (top left) RLWM; (top right) RLWM, coarse mesh; (bottom left)
EQWM; and (bottom right) WRLES (Gloerfelt & Cinnella, 2019).

Figure 2 shows the contours of the mean velocity in $x$ direction and the
mean-flow streamlines. The flow separates on the leeward side of the hill due
to a strong adverse pressure gradient (APG), and a shear layer is generated
near the top of the hill. The flow reattaches in the middle section of the
channel, and as the flow approaches the windward side of the downstream hill,
it is subjected to a strong favorable pressure gradient (FPG) and accelerates
rapidly. The simulations with RLWM successfully capture the separation bubble
on the leeward side of the hill and yield more accurate results than the EQWM
(see Table 1 for quantitative comparison).

Figure 3: (a) Mean skin friction coefficient and (b) mean pressure coefficient
along the bottom wall at $Re_{H}=10595$. Lines indicate  , RLWM;  , RLWM,
coarse mesh;  , RLWM, fine mesh;  , EQWM; $\circ$, DNS (Krank et al., 2018).

The predictions of $C_{f}$ and the mean pressure coefficient $C_{p}$ are shown
in Figure 3. The mean pressure coefficient is defined as $C_{p}=(\langle
p_{w}\rangle-\langle p_{ref}\rangle)/(0.5\rho U_{b}^{2})$, where the pressure
at $x/H=0$ on the top wall is chosen as reference pressure $p_{ref}$ (Krank et
al., 2018). Regarding $C_{f}$, the results from RLWM simulations are in
reasonable agreement with the direct numerical simulation (DNS) data, with
large deviations found only near the top of the hill where the skin friction
reaches its maximum value and then rapidly decreases to a negative value. As
the mesh resolution increases, it yields a more accurate prediction of skin
friction on the leeward side, but the maximum value of the skin friction on
the windward side of the hill still is overpredicted. However, the results are
better than the EQWM simulation, which largely underpredicts the skin friction
on the windward side and overpredicts the skin friction on the leeward side.
Furthermore, the mean locations of the separation and reattachment points
(listed in Table 1) are better predicted by the RLWM, consistent with the
streamline shown in Figure 2. Additionally, note that the separation bubble
size changes as the mesh resolution increases or decreases from the baseline,
and more details about the velocity field will be discussed later in this
section.

Figure 4:  Mean velocity and Reynolds stress components at $Re_{H}=10595$: (a)
streamwise velocity, (b) vertical velocity, (c) streamwise Reynolds stress,
(d) vertical Reynolds stress, and (e) Reynolds shear stress. Lines indicate  ,
RLWM;  , RLWM, coarse mesh;  , RLWM, fine mesh;  , EQWM; $\circ$, DNS (Krank
et al., 2018).

All simulations capture the qualitative trend of the mean $C_{p}$ on the
bottom wall including the APG and FPG regimes, but relatively larger
deviations among the simulation cases are visible near the top of the hill
($x/H\geq 8.5$ or $x/H\leq 0.5$) where the pressure sees a sudden change from
strong FPG to strong APG and the flow separation emerges. Overall, the RLWM
provides more accurate predictions of $C_{f}$ and $C_{p}$ than the EQWM.

Quantitative comparison of mean velocity and Reynolds stress components at
five streamwise locations ($x/H=0.05,2,4,6,\text{ and }8$) are shown in Figure
4. The results of the mean velocity profiles and Reynolds stress profiles from
the RLWM simulations agree reasonably well with each other and with the
reference DNS data. The discrepancies are visible for the regions near channel
center. Note that the prediction of the velocity field depends not only on the
wall boundary conditions but also on the SGS model, and the varying
performance of the SGS model in different meshes may contribute to the
inconsistency of those three RLWM simulations. Furthermore, the prediction
from EQWM simulation is less ideal, particularly for the mean streamwise
velocity $\langle u_{x}\rangle/U_{b}$.

To better understand the mechanism of the trained model, we examine the state-
action maps, which are the probability density functions (PDFs) of the
likelihood that the model takes a particular action conditioned on the
occurrence of positive rewards. Figure 5 shows the maps based on instantaneous
states and actions at three streamwise positions $x/H=0.025,2,\text{ and }8.3$
which are located at the top of the hill, within the separation bubble, and
the windward side of the hill, respectively. Overall, the action contour lines
for increasing and decreasing $\nu_{t,w}$ are well separated with respect to
the mean quantities obtained from the WRLES simulation (Gloerfelt & Cinnella,
2019), which illustrates that the model is able to distinguish flow states and
provide appropriate actions. Some minor overlaps of the contour lines not
clearly separated by the mean relation from WRLES are visible at the top of
the hill [see Figures 5(a,d)] where the APG is strong. This implies that the
capability of the trained model to robustly adjust the wall-shear stress
through appropriate actions is less ideal within this regime.

Figure 5: PDFs of states for RLWM conditioned to events with $r>0.1$ and
$a<0.9996$ (blue) or $a>1.0004$ (red) at (a,d) $x/H=0.025$, (b,e) $x/H=2$, and
(c,f) $x/H=8.3$ at $Re_{H}=10595$: PDFs of states (a–c) $u_{s}^{*}$ and
$y_{n}^{*}$ and (d–f) $S_{12}^{*}$ and $u_{s}^{*}$. The black line is the mean
relation of the states from WRLES (Gloerfelt & Cinnella, 2019). Contour levels
are $25$, $50$, and $75\%$ of the maximum value.

###  3.2 Testing for flow over periodic hills at higher Reynolds numbers

In this section, the RLWM is applied to a WMLES of periodic-hill channel flow
at $Re_{H}=19000$ and $37000$. The simulations are conducted by using the
baseline mesh ($128\times 64\times 64$) and the coarse mesh ($64\times
32\times 32$), and the implementation of RLWM is similar to the simulations at
$Re_{H}=10595$. All simulations are run for about 60 FTTs after transients.
The results from EQWM with baseline mesh and the WRLES (Gloerfelt & Cinnella,
2019) are included for comparison.

The mean skin friction coefficients along the bottom wall at $Re_{H}=19000$
and $37000$ are shown in Figure 6(a,b). The distributions of $C_{f}$ at higher
Reynolds numbers have a similar shape as the one shown in Figure 3(a). As the
Reynolds number increases, the maximum value of the $C_{f}$ on the windward
side of the hill decreases. The results from the two RLWM simulations are in
reasonable agreement with the WRLES results. Regarding the separation point,
the predicted locations from the RLWM simulations are further downstream than
that of the WRLES. The predicted reattachment location from the RLWM
simulation with baseline mesh is slightly downstream from that of the WRLES,
but the prediction from the RLWM simulation with coarse mesh is further
upstream. Compared to the EQWM simulation, where the separation bubble size is
largely underpredicted and the maximum value of $C_{f}$ on the windward side
of the hill is much smaller than the value of WRLES, RLWM is more accurate.

Figure 6: Mean skin friction coefficient along the bottom wall at (a)
$Re_{H}=19000$ and (b) $Re_{H}=37000$. Lines indicate  , RLWM;  , RLWM, coarse
mesh;  , EQWM; , WRLES (Gloerfelt & Cinnella, 2019).

Figure 7 shows the profiles of streamwise components of mean velocity and
Reynolds stress at five streamwise stations ($x/H=0.05,2,4,6,\text{ and }8$)
for $Re_{H}=19000$ and $37000$. The deviations among the profiles from RLWM
simulations and the WRLES profiles grow as the Reynolds number increases,
albeit slightly, which implies the degradation of performance in the RLWM.
Particularly for the Reynolds stress, the EQWM outperforms the RLWM even
though the prediction of $C_{f}$ is inaccurate, as shown in Figure 6.
Furthermore, at $Re_{H}=19000$ the coarse-mesh RLWM simulation results agree
better with the WRLES results than the baseline results. However, the
performance of the SGS model affects the velocity field away from the wall
more than the wall model, and this could be an indication that a better SGS
model is necessary.

Figure 7: Mean velocity and Reynolds stress components at (a,b) $Re_{H}=19000$
and (c,d) $Re_{H}=37000$: (a,c) mean streamwise velocity and (b,d) streamwise
Reynolds stress. Lines indicate  , RLWM;  , RLWM, coarse mesh;  , EQWM; ,
WRLES (Gloerfelt & Cinnella, 2019).

##  4 Conclusions

In this work, a wall model that can adapt to various pressure-gradient effects
is developed for turbulent flow over periodic hills using multi-agent
reinforcement learning. The model behaves as a control policy for wall eddy
viscosity to predict the correct wall-shear stress. The optimized policy of
the wall model is learned through the training process based on LES of low-
Reynolds-number flow over periodic hills with cooperating agents using the
recovery of the correct wall-shear stress as a reward.

The developed wall model is first validated in the LES of the periodic hill
configuration at the same Reynolds number of model training. The wall model
provides good predictions of mean wall-shear stress, mean wall pressure, and
the mean velocity as well as Reynolds stress in the flow field. The test
results also show that the developed model outperforms the EQWM. The
performance of the developed model is further evaluated at $Re_{H}=19000$ and
$37000$. Good predictions are obtained for the mean wall-shear stress.

Overall, we have extended the MARL-based wall modeling paradigm and
demonstrated that it is robust in flows with pressure gradients. Additional
investigations into the generalizability of the developed wall model for
different geometries and flow configurations need to be carried out, and the
grid sensitivity of the model needs to be further analyzed.

### Acknowledgments

The authors acknowledge the XSEDE computing resource from PSC under project
PHY210119. The authors would like to thank Dr. Xiang Yang for their insightful
remarks on the report.

## References

  * Bae & Koumoutsakos (2022) Bae, H. J. & Koumoutsakos, P. 2022 Scientific multi-agent reinforcement learning for wall-models of turbulent flows. Nat. Commun. 13, 1443.
  * Bae & Lozano-Durán (2021) Bae, H. J. & Lozano-Durán, A. 2021 Effect of wall boundary conditions on a wall-modeled large-eddy simulation in a finite-difference framework. Fluids 6, 112.
  * Balakumar et al. (2014) Balakumar, P., Park, G. & Pierce, B. 2014 DNS, LES, and wall-modeled LES of separating flow over periodic hills. Proceedings of the Summer Program, Center for Turbulence Research, Stanford University, pp. 407–415.
  * Germano et al. (1991) Germano, M., Piomelli, U., Moin, P. & Cabot, W. H. 1991 A dynamic subgrid-scale eddy viscosity model. Phys. Fluids 3, 1760–1765.
  * Gloerfelt & Cinnella (2019) Gloerfelt, X. & Cinnella, P. 2019 Large eddy simulation requirements for the flow over periodic hills. Flow Turbul. Combust. 103, 55–91.
  * Kawai & Larsson (2012) Kawai, S. & Larsson, J. 2012 Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy. Phys. Fluids 24, 015105.
  * Krank et al. (2018) Krank, B., Kronbichler, M. & Wall, W. A. 2018 Direct numerical simulation of flow over periodic hills up to $\text{Re}_{H}=10{,}595$. Flow Turbul. Combust. 101, 521–551.
  * Lilly (1992) Lilly, D. K. 1992 A proposed modification of the Germano subgrid-scale closure method. Phys. Fluids A 4, 633–635.
  * Mellen et al. (2000) Mellen, C. P., Fröhlich, J. & Rodi, W. 2000 Large eddy simulation of the flow over periodic hills. In Proceedings of the 16th IMACS World Congress, pp. 21–25.
  * Novati & Koumoutsakos (2019) Novati, G. & Koumoutsakos, P. 2019 Remember and forget for experience replay. In Proceedings of the 36th International Conference on Machine Learning, pp. 4851–4860.
  * You et al. (2008) You, D., Ham, F. & Moin, P. 2008 Discrete conservation principles in large-eddy simulation with application to separation control over an airfoil. Phys. Fluids 20, 101515.


#  Anisotropic Cosmology in the Local Limit of Nonlocal Gravity

Javad Tabatabaei1<EMAIL_ADDRESS>Abdolali Banihashemi1
<EMAIL_ADDRESS>Shant Baghram1<EMAIL_ADDRESS>Bahram
Mashhoon2,3<EMAIL_ADDRESS>1Department of Physics, Sharif University
of Technology, Tehran 11155-9161, Iran
2School of Astronomy, Institute for Research in Fundamental Sciences (IPM),
Tehran 19395-5531, Iran
3Department of Physics and Astronomy, University of Missouri, Columbia,
Missouri 65211, USA

###### Abstract

Within the framework of the local limit of nonlocal gravity (NLG), we
investigate a class of Bianchi type I spatially homogeneous but anisotropic
cosmological models. The modified field equations are presented in this case
and some special solutions are discussed in detail. This modified gravity
theory contains a susceptibility function $S(x)$ such that general relativity
(GR) is recovered for $S=0$. In the modified anisotropic cosmological models,
we explore the contribution of $S(t)$ and its temporal derivative to the local
anisotropic cosmic acceleration. The implications of our results for
observational cosmology are briefly discussed.

Gravitation

###### pacs:

04.20.Cv

##  I Introduction

In the current $\Lambda$CDM model of cosmology, the energy content of the
universe consists of about 70% dark energy, about 25% dark matter, and about
5% visible matter. The model is based on the standard spatially homogeneous
and isotropic Friedmann–Lemaître–Robertson–Walker (FLRW) cosmological
solutions of Einstein’s general relativity theory Einstein . The _dark_
features, whose nature and origin are unknown, provide the motivation to
modify and extend general relativity (GR) on galactic scales and beyond in
order to account for observational data purely on the basis of the
gravitational physics of the extended GR without recourse to dark ingredients.

To modify the current benchmark model of cosmology, we consider nonlocal
gravity (NLG) theory Hehl:2008eu ; Hehl:2009es ; BMB , a classical history-
dependent generalization of GR that bears a formal resemblance to the nonlocal
electrodynamics of media Hop ; Poi ; Jackson ; L+L ; HeOb ;
VanDenHoogen:2017nyy . It is important to digress here and mention that there
are indeed various other approaches to nonlocal gravitation. For the sake of
brevity, we only refer to some examples and their cosmological implications.
Nonlocally modified extensions of GR can be generated by the addition of
functions of $\Box$, as in infinite derivative theories, or functions of
$\Box^{-1}$ to the Einstein-Hilbert action. Here, $\Box$ is the d’Alembert-
Beltrami operator. Cosmological implications of such theories have been
investigated by a number of authors in connection with dynamic dark energy and
accelerated expansion of the universe; for instance, see Maggiore:2014sia ;
Capozziello:2021krv  and the references cited therein. Moreover, cosmological
solutions of nonlocal infinite derivative theories have been studied that
involve anisotropic bouncing models Kumar:2021mgc  or an interplay between
dark matter and dark energy Dimitrijevic:2019pct . Quantum field theory
provides the motivation for a different class of nonlocal theories of
gravitation. Higher curvature nonlocal gravity theories have recently been
reviewed in Koshelev:2023elc  and a generalized nonlocal quantum gravity
theory has been formulated within the framework of inflationary cosmology. On
the other hand, the phenomenological approach of Deser and Woodard has been
based on an effective quantum gravitational action and has been designed to
explain cosmic acceleration without dark energy; see Deser:2019lmm  and the
references cited therein. Furthermore, primordial bouncing cosmology and
anisotropy have been investigated within the framework of the Deser-Woodard
nonlocal gravity model in Chen:2019wlu .

We now return to our classical model of nonlocal gravity that is patterned
after the nonlocal electrodynamics of media. Nonlocal gravity (NLG) is a
tetrad theory, where the gravitational potentials are given by the 16
components of a preferred orthonormal tetrad frame field. The extended
geometric framework of NLG is based on the Weitzenböck connection We , which
renders the spacetime a parallelizable manifold. Within the framework of
teleparallelism Maluf:2013gaa ; Aldrovandi:2013wha ; BlHe ; Itin:2018dru , it
is possible to express GR using the Weitzenböck torsion tensor. This
teleparallel equivalent of general relativity (TEGR) is a gauge theory of the
Abelian group of spacetime translations Cho . The formal similarity between
TEGR and electrodynamics can be employed to introduce nonlocality into GR via
constitutive kernels Hehl:2008eu ; Hehl:2009es . In NLG, the gravitational
field is local, but the theory involves an average of the field over past
events resulting in 16 partial integro-differential field equations
Puetzfeld:2019wwo ; Mashhoon:2022ynk . No exact nontrivial solution of NLG is
known at present Bini:2016phe ; however, the linear regime of the theory has
been extensively studied. Nonlocal gravity, in its Newtonian regime, simulates
dark matter. It is therefore possible to account for the gravitational effects
in the solar system as well as in nearby galaxies and clusters of galaxies
Rahvar:2014yta ; Chicone:2015coa ; Roshan:2021ljs ; Roshan:2022zov ;
Roshan:2022ypk . A comprehensive account of these studies is contained in BMB
.

NLG is rather intricate and to study its cosmological implications, we resort
to its local limit, which is easier to analyze. In Section II, we present a
brief account of the modified GR field equations in the local limit of NLG.
For a more detailed treatment of this limiting situation, see
Tabatabaei:2022tbq ; Tabatabaei:2023qxw , where spatially homogeneous and
isotropic (FLRW) cosmological models were investigated in this modified TEGR
scheme. To explore anisotropy in the Hubble flow, we present in Section III
modified gravitational field equations for a Bianchi type I class of time-
dependent spatially homogeneous but anisotropic spacetimes within the
framework of the local limit of NLG. The field equations contain a
susceptibility function $S(t)$ with $1+S>0$ and $dS/dt\neq 0$ that is
characteristic of the dynamic spacetime background. For $S=0$, we recover the
GR field equations. We show that de Sitter and Kasner spacetimes are not
solutions of the modified field equations unless $S(t)$ is independent of
time, which is not physically reasonable. Explicit solutions of the modified
field equations are studied in the next two sections. A well-known class of
dynamic solutions of GR for dust with vanishing cosmological constant is
extended to the local limit of NLG in Section IV. The new solutions contain
the time-dependent function $S(t)$ and allow the possibility of exploring the
dependence of anisotropic acceleration on $S(t)$. These modified cosmological
models are locally anisotropic but tend to the isotropic modified Einstein-de
Sitter model at late times ($t\to\infty$). Similarly, we study the solution of
the modified field equations for a spacetime dominated by dark energy in
Section V and explore anisotropic cosmic acceleration in this cosmological
model, which eventually becomes isotropic as well. The presence of $S(t)$
could be responsible for certain new “dark” features of accelerating bulk
flows in the local universe.

Anisotropy of the Hubble flow would indicate a significant departure from the
presumed large-scale spatial homogeneity and isotropy of the standard FLRW
cosmology. On the other hand, there is recent observational evidence in
support of _local_ anisotropic cosmic acceleration Colin:2019opb ;
Secrest:2020has ; Secrest:2022uvx ; Solanki:2023yoa ; Perivolaropoulos:2023tdt
. The purpose of the present paper is to study theoretically the possible
contribution of the susceptibility function $S(t)$ to anisotropic features of
the Hubble flow.

##  II Local Limit of NLG

We consider a spacetime manifold as in general relativity (GR). In an
admissible system of coordinates $x^{\mu}$, the spacetime metric can be
written as

$ds^{2}=g_{\mu\nu}(x)\,dx^{\mu}\,dx^{\nu}\,.$ (1)

Here, Greek indices run from 0 to 3, Latin indices run from 1 to 3, and the
signature of the metric is +2; moreover, we employ units such that $c=1$. As
in GR, the world lines of free test particles and null rays are geodesics of
the spacetime manifold. We assume the existence of a preferred set of
observers in this gravitational field. The observers have adapted orthonormal
tetrads $e^{\mu}{}_{\hat{\alpha}}(x)$,

$g_{\mu\nu}(x)\,e^{\mu}{}_{\hat{\alpha}}(x)\,e^{\nu}{}_{\hat{\beta}}(x)=\eta_{\hat{\alpha}\hat{\beta}}\,,$
(2)

where $\eta_{\alpha\beta}={\rm diag}(-1,1,1,1)$ is the Minkowski metric
tensor. In our convention, indices without hats are normal spacetime indices,
while hatted indices indicate the tetrad axes in the local tangent space.

We employ the tetrad frame field to define the curvature-free Weitzenböck
connection,

$\Gamma^{\mu}_{\alpha\beta}=e^{\mu}{}_{\hat{\rho}}~{}\partial_{\alpha}\,e_{\beta}{}^{\hat{\rho}}\,.$
(3)

Let $\nabla$ denote covariant differentiation with respect to the Weitzenböck
connection; then, $\nabla_{\nu}\,e_{\mu}{}^{\hat{\alpha}}=0$, thereby the
preferred tetrad frames are parallel throughout the gravitational field and
provide a natural scaffolding for the spacetime manifold. The spacetime is
thus a parallelizable manifold by the Weitzenböck connection. In this
framework of teleparallelism, two distant vectors are considered parallel if
they have the same local components relative to their preferred tetrad frames.
Moreover, it follows from the tetrad orthonormality relation that the
Weitzenböck connection is metric compatible, namely,
$\nabla_{\mu}\,g_{\alpha\beta}=0$.

The difference between two connections on the same manifold is a tensor. We
define the _torsion_ tensor that corresponds to the Weitzenböck connection by

$C_{\mu\nu}{}^{\alpha}=\Gamma^{\alpha}_{\mu\nu}-\Gamma^{\alpha}_{\nu\mu}=e^{\alpha}{}_{\hat{\beta}}\Big{(}\partial_{\mu}e_{\nu}{}^{\hat{\beta}}-\partial_{\nu}e_{\mu}{}^{\hat{\beta}}\Big{)}\,.$
(4)

In the extended GR framework, we have the Weitzenböck connection as well as
the symmetric Levi-Civita connection,

$^{0}\Gamma^{\mu}_{\alpha\beta}=\frac{1}{2}g^{\mu\nu}\,(g_{\nu\alpha,\beta}+g_{\nu\beta,\alpha}-g_{\alpha\beta,\nu})\,.$
(5)

We use a left superscript “0” to refer to geometric quantities directly
derived from the Levi-Civita connection. The _contorsion_ tensor is then
defined by

$K_{\mu\nu}{}^{\alpha}=\,^{0}\Gamma^{\alpha}_{\mu\nu}-\Gamma^{\alpha}_{\mu\nu}\,,$
(6)

which is related to the torsion tensor through the metric compatibility of the
Weitzenböck connection. In fact,

$K_{\mu\nu\rho}=\frac{1}{2}\,(C_{\mu\rho\nu}+C_{\nu\rho\mu}-C_{\mu\nu\rho})\,.$
(7)

The Levi-Civita connection given by the Christoffel symbol is the sum of the
Weitzenböck connection and the contorsion tensor. One can therefore express
the Einstein tensor ${}^{0}G_{\mu\nu}$ and the gravitational field equations
of GR in terms of the teleparallelism framework resulting in the teleparallel
equivalent of GR, namely, TEGR BMB . Indeed, we find

$\displaystyle{{}^{0}}G_{\mu\nu}=\frac{\kappa}{\sqrt{-g}}\Big{[}e_{\mu}{}^{\hat{\gamma}}\,g_{\nu\alpha}\,\frac{\partial}{\partial
x^{\beta}}\,\mathfrak{H}^{\alpha\beta}{}_{\hat{\gamma}}-\Big{(}C_{\mu}{}^{\rho\sigma}\,\mathfrak{H}_{\nu\rho\sigma}-\frac{1}{4}\,g_{\mu\nu}\,C^{\alpha\beta\gamma}\,\mathfrak{H}_{\alpha\beta\gamma}\Big{)}\Big{]}\,$
(8)

and Einstein’s field equations expressed in terms of torsion thus become the
TEGR field equations

$\frac{\partial}{\partial
x^{\nu}}\,\mathfrak{H}^{\mu\nu}{}_{\hat{\alpha}}+\frac{\sqrt{-g}}{\kappa}\,\Lambda\,e^{\mu}{}_{\hat{\alpha}}=\sqrt{-g}\,(T_{\hat{\alpha}}{}^{\mu}+\mathbb{T}_{\hat{\alpha}}{}^{\mu})\,,$
(9)

where $\Lambda$ is the cosmological constant and $\kappa:=8\pi G$. Here, we
define the auxiliary torsion field $\mathfrak{H}_{\mu\nu\rho}$ by means of the
auxiliary torsion tensor $\mathfrak{C}_{\alpha\beta\gamma}$, namely,

$\mathfrak{H}_{\mu\nu\rho}:=\frac{\sqrt{-g}}{\kappa}\,\mathfrak{C}_{\mu\nu\rho}\,,\qquad\mathfrak{C}_{\alpha\beta\gamma}:=C_{\alpha}\,g_{\beta\gamma}-C_{\beta}\,g_{\alpha\gamma}+K_{\gamma\alpha\beta}\,.$
(10)

Moreover, $C_{\mu}:=C^{\alpha}{}_{\mu\alpha}=-C_{\mu}{}^{\alpha}{}_{\alpha}$
is the torsion vector. As in GR, $T_{\mu\nu}$ is the symmetric energy-momentum
tensor of matter. We interpret $\mathbb{T}_{\mu\nu}$ in Equation (9) to be the
traceless energy-momentum tensor of the gravitational field given by

$\mathbb{T}_{\mu\nu}:=(\sqrt{-g})^{-1}\,(C_{\mu\rho\sigma}\,\mathfrak{H}_{\nu}{}^{\rho\sigma}-\tfrac{1}{4}g_{\mu\nu}\,C_{\rho\sigma\delta}\,\mathfrak{H}^{\rho\sigma\delta})\,.$
(11)

This version of GR, namely, TEGR, is the gauge theory of the 4-parameter
Abelian group of spacetime translations Cho ; therefore, though nonlinear, it
bears a certain resemblance to Maxwell’s electrodynamics.

In analogy with the electrodynamics of media, we can consider the torsion
tensor in the form
$C_{\mu\nu}{}^{\hat{\alpha}}=\partial_{\mu}e_{\nu}{}^{\hat{\alpha}}-\partial_{\nu}e_{\mu}{}^{\hat{\alpha}}$
to be similar to the Faraday tensor, while the relationship between
$\mathfrak{H}_{\mu\nu\rho}$ and the torsion tensor in Equation (10) can be
viewed as the local constitutive relation of TEGR. Let us recall that in
Maxwell’s electrodynamics, the constitutive relation may change, but the field
equations remain the same. We adopt the same approach for the purpose of
modifying Einstein’s theory. That is, we modify TEGR by introducing a tensor
$N_{\mu\nu\rho}=-N_{\nu\mu\rho}$ that changes the constitutive relation of
TEGR as follows

$\mathcal{H}_{\mu\nu\rho}=\frac{\sqrt{-g}}{\kappa}(\mathfrak{C}_{\mu\nu\rho}+N_{\mu\nu\rho})\,.$
(12)

To obtain the field equations of modified TEGR, we simply replace
$\mathfrak{H}$ in Equations (9) and (11) by $\mathcal{H}$. The gravitational
field equations of extended GR based on the new tensor field $N_{\mu\nu\rho}$
now take the form

$\frac{\partial}{\partial
x^{\nu}}\,\mathcal{H}^{\mu\nu}{}_{\hat{\alpha}}+\frac{\sqrt{-g}}{\kappa}\,\Lambda\,e^{\mu}{}_{\hat{\alpha}}=\sqrt{-g}\,(T_{\hat{\alpha}}{}^{\mu}+\mathcal{T}_{\hat{\alpha}}{}^{\mu})\,,$
(13)

where $\mathcal{T}_{\mu\nu}$ is the traceless energy-momentum tensor of the
gravitational field. The gravitational energy-momentum tensor is modified by
the presence of $N_{\mu\nu\rho}$; hence, we introduce a traceless tensor
$Q_{\mu\nu}$ that indicates this difference, namely,

$\kappa\,\mathcal{T}_{\mu\nu}=\kappa\,\mathbb{T}_{\mu\nu}+Q_{\mu\nu}\,,$ (14)

where

$Q_{\mu\nu}:=C_{\mu\rho\sigma}N_{\nu}{}^{\rho\sigma}-\frac{1}{4}\,g_{\mu\nu}\,C_{\delta\rho\sigma}N^{\delta\rho\sigma}\,.$
(15)

The total energy-momentum conservation law takes the form

$\frac{\partial}{\partial
x^{\mu}}\,\Big{[}\sqrt{-g}\,(T_{\hat{\alpha}}{}^{\mu}+\mathcal{T}_{\hat{\alpha}}{}^{\mu}-\frac{\Lambda}{\kappa}\,e^{\mu}{}_{\hat{\alpha}})\Big{]}=0\,.$
(16)

It is interesting to see how $N_{\mu\nu\rho}$ modifies GR field equations; to
this end, we substitute

$\mathfrak{H}_{\mu\nu\rho}=\mathcal{H}_{\mu\nu\rho}-\frac{\sqrt{-g}}{\kappa}N_{\mu\nu\rho}\,$
(17)

in the Einstein tensor (8) and employ modified TEGR field Equation (13) to get

$^{0}G_{\mu\nu}+\Lambda g_{\mu\nu}=\kappa
T_{\mu\nu}+Q_{\mu\nu}-\mathcal{N}_{\mu\nu}\,,$ (18)

where $\mathcal{N}_{\mu\nu}$ is a tensor defined by

$\mathcal{N}_{\mu\nu}:=g_{\nu\alpha}e_{\mu}{}^{\hat{\gamma}}\frac{1}{\sqrt{-g}}\frac{\partial}{\partial
x^{\beta}}\,(\sqrt{-g}N^{\alpha\beta}{}_{\hat{\gamma}})\,.$ (19)

Therefore, to find the field equations of modified GR, we must add
$Q_{\mu\nu}-\mathcal{N}_{\mu\nu}$ to the right-hand side of Einstein’s field
equations of GR.

Finally, we have to relate $N_{\mu\nu\rho}=-N_{\nu\mu\rho}$ to the torsion
tensor. In NLG, the components of $N_{\mu\nu\rho}$ measured by the preferred
observers of the theory with adapted tetrads $e^{\mu}{}_{\hat{\alpha}}$ are
associated with the corresponding measured components of $X_{\mu\nu\rho}$ that
are directly connected to the torsion tensor, and its expression has been
discussed in detail in BMB . That is Puetzfeld:2019wwo ; Mashhoon:2022ynk ,

$N_{\hat{\mu}\hat{\nu}\hat{\rho}}(x)=\int\mathcal{K}(x,x^{\prime})\,X_{\hat{\mu}\hat{\nu}\hat{\rho}}(x^{\prime})\sqrt{-g(x^{\prime})}\,d^{4}x^{\prime}\,,$
(20)

where

$X_{\hat{\mu}\hat{\nu}\hat{\rho}}=\mathfrak{C}_{\hat{\mu}\hat{\nu}\hat{\rho}}+\check{p}\,(\check{C}_{\hat{\mu}}\,\eta_{\hat{\nu}\hat{\rho}}-\check{C}_{\hat{\nu}}\,\eta_{\hat{\mu}\hat{\rho}})\,.$
(21)

Here, $\mathcal{K}(x,x^{\prime})$ is the basic causal kernel of NLG that in
essence must be determined via observation BMB , $\check{p}\neq 0$ is a
constant dimensionless parameter, and $\check{C}^{\mu}$ is the torsion
pseudovector,

$\check{C}_{\mu}=\frac{1}{3!}C^{\alpha\beta\gamma}\,E_{\alpha\beta\gamma\mu}\,,$
(22)

where $E_{\alpha\beta\gamma\delta}$ is the Levi-Civita tensor.

Nonlocal gravity (NLG) is thus a classical extension of GR that is highly
nonlinear as well. Linearized NLG has been investigated in detail BMB . Within
the Newtonian regime of NLG, it appears possible to account for the rotation
curves of nearby spiral galaxies as well as for the solar system data
Rahvar:2014yta ; Chicone:2015coa ; Roshan:2021ljs ; Roshan:2022zov ;
Roshan:2022ypk . Beyond the linear domain, no exact solution is known except
for the trivial result that in the absence of gravity we have Minkowski
spacetime Bini:2016phe . On the other hand, it is possible that certain
nonlinear features of NLG that belong to the strong-field regimes such as
those involving black holes or cosmological models may indeed survive in the
local limit of the theory. It is therefore interesting to explore this
limiting case of NLG.

To come up with the local limit of NLG, let us assume that the kernel in
Equation (20) is proportional to the 4D Dirac delta function, namely,

$\mathcal{K}(x,x^{\prime}):=\frac{S(x)}{\sqrt{-g(x)}}\,\delta(x-x^{\prime})\,;$
(23)

then, $N_{\mu\nu\rho}(x)=S(x)X_{\mu\nu\rho}$, where $S(x)$ is a dimensionless
scalar function. Therefore,

$N_{\mu\nu\rho}(x)=S(x)\,[\mathfrak{C}_{\mu\nu\rho}(x)+\check{p}\,(\check{C}_{\mu}\,g_{\nu\rho}-\check{C}_{\nu}\,g_{\mu\rho})]\,$
(24)

and the constitutive relation takes the form

$\mathcal{H}_{\mu\nu\rho}=\frac{\sqrt{-g}}{\kappa}[(1+S)\,\mathfrak{C}_{\mu\nu\rho}+S\,\check{p}\,(\check{C}_{\mu}\,g_{\nu\rho}-\check{C}_{\nu}\,g_{\mu\rho})]\,.$
(25)

Here, the susceptibility function $S(x)$ is a characteristic of the background
spacetime just as $\epsilon(x)$ and $\mu(x)$ are features of the medium in
electrodynamics. In general, the local electric permittivity $\epsilon(x)$ and
magnetic permeability $\mu(x)$ functions are expected to preserve significant
features of the electrodynamics of media such as spatial symmetries and
temporal dependence. Similarly, $S(x)$ is expected to preserve the
characteristics of the background spacetime. Ultimately, $S(x)$ must be
determined based on observational data Tabatabaei:2022tbq ; Tabatabaei:2023qxw
.

For $S(x)=0$, we recover TEGR; otherwise, we have a natural generalization of
GR that contains a new function $S(x)$. Indeed, Equation (25) implies that to
have GR as a limit, we must impose the requirement that $1+S>0$. In this local
limit of nonlocal gravity, explicit deviations from locality have vanished;
however, nontrivial aspects of NLG may have survived through $S(x)$ that would
be interesting to study. Consequently, we explore the cosmological
implications of this local limit of NLG. Spatially homogeneous and isotropic
(FLRW) cosmological models have been treated in Tabatabaei:2022tbq ;
Tabatabaei:2023qxw  in connection with $H_{0}$ tension. Therefore, we
concentrate here on a class of spatially homogeneous but anisotropic
spacetimes.

##  III Anisotropic Models

Let us consider a Bianchi type I model with a metric of the form

$\displaystyle ds^{2}=-dt^{2}+X^{2}dx^{2}+Y^{2}dy^{2}+Z^{2}dz^{2},$ (26)

where $X$, $Y$ and $Z$ are functions of time $t$. This spacetime is spatially
homogeneous with three spacelike commuting Killing vector fields
$\partial_{x}$, $\partial_{y}$, and $\partial_{z}$. A detailed discussion of
such spacetimes is contained in Chapter 13 of SKM ; for a recent discussion
within the context of teleparallelism, see Coley:2023ibm .

Einstein’s gravitational field equations are

$\displaystyle^{0}G_{\mu\nu}+\Lambda\,g_{\mu\nu}=8\pi
G\,T_{\mu\nu}\,,\qquad^{0}G_{\mu\nu}=~{}^{0}R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\,^{0}R\,,$
(27)

where $T_{\mu\nu}$ is assumed to be due to the presence of a comoving perfect
fluid with density $\rho(t)$ and pressure $P(t)$,

$\displaystyle T_{\mu\nu}=(\rho+P)U_{\mu}U_{\nu}+Pg_{\mu\nu},$ (28)

and $\Lambda$ is the cosmological constant. With respect to the system of
coordinates $x^{\mu}=(t,x,y,z)$, the perfect fluid is comoving with
$U^{\mu}=\delta^{\mu}{}_{0}$; hence,

$\displaystyle T_{\mu\nu}={\rm diag}(\rho,PX^{2},PY^{2},PZ^{2})\,.$ (29)

Moreover, the Einstein tensor ${}^{0}G_{\mu\nu}$ is diagonal as well with
components

$\displaystyle^{0}G_{00}$
$\displaystyle=\frac{\dot{X}\dot{Y}}{XY}+\frac{\dot{Y}\dot{Z}}{YZ}+\frac{\dot{Z}\dot{X}}{ZX}\,,$
(30) $\displaystyle^{0}G_{11}$
$\displaystyle=-X^{2}\left(\frac{\ddot{Y}}{Y}+\frac{\ddot{Z}}{Z}+\frac{\dot{Y}\dot{Z}}{YZ}\right)\,,$
(31) $\displaystyle^{0}G_{22}$
$\displaystyle=-Y^{2}\left(\frac{\ddot{Z}}{Z}+\frac{\ddot{X}}{X}+\frac{\dot{Z}\dot{X}}{ZX}\right)\,,$
(32) $\displaystyle^{0}G_{33}$
$\displaystyle=-Z^{2}\left(\frac{\ddot{X}}{X}+\frac{\ddot{Y}}{Y}+\frac{\dot{X}\dot{Y}}{XY}\right)\,.$
(33)

It is interesting to work out the Kretschmann scalar $\mathcal{K}$,

$\displaystyle\mathcal{K}=~{}^{0}R_{\mu\nu\rho\sigma}\,^{0}R^{\mu\nu\rho\sigma}\,,$
(34)

for metric (26). The result is

$\displaystyle\frac{1}{4}\,\mathcal{K}=\left(\frac{\ddot{X}}{X}\right)^{2}+\left(\frac{\ddot{Y}}{Y}\right)^{2}+\left(\frac{\ddot{Z}}{Z}\right)^{2}+\left(\frac{\dot{X}\dot{Y}}{XY}\right)^{2}+\left(\frac{\dot{Y}\dot{Z}}{YZ}\right)^{2}+\left(\frac{\dot{Z}\dot{X}}{ZX}\right)^{2}\,.$
(35)

Detailed discussions of the GR solutions of these models with $\Lambda=0$ for
dust ($P=0$) can be found, for instance, in HS1 ; HS2 , Section 5.4 of Ref. HE
, and Section 12.15 of Ref. PK . We give a brief description of these
solutions in Section IV in connection with cosmic deceleration.

We are interested in the extended GR framework. Therefore, consider the class
of observers that are spatially at rest with adapted tetrad
$e^{\mu}{}_{\hat{\alpha}}$ field given by

$\displaystyle e^{\mu}{}_{\hat{0}}=(1,0,0,0)\,,\quad
e^{\mu}{}_{\hat{1}}=(0,\frac{1}{X},0,0)\,,\quad
e^{\mu}{}_{\hat{2}}=(0,0,\frac{1}{Y},0)\,,\quad
e^{\mu}{}_{\hat{3}}=(0,0,0,\frac{1}{Z})\,,$ (36)

where the spatial axes point along the Cartesian coordinate directions. We
have

$\displaystyle e_{\mu}{}^{\hat{0}}=(1,0,0,0)\,,\quad
e_{\mu}{}^{\hat{1}}=(0,X,0,0)\,,\quad e_{\mu}{}^{\hat{2}}=(0,0,Y,0)\,,\quad
e_{\mu}{}^{\hat{3}}=(0,0,0,Z)\,.$ (37)

We compute the Weitzenböck torsion tensor (4) in this case and we find
$C_{\mu\nu}{}^{0}=0$, $C_{ij}{}^{k}=0$, and the only nonzero components can be
obtained from

$\displaystyle C_{01}{}^{1}=\frac{\dot{X}}{X}\,,\qquad
C_{02}{}^{2}=\frac{\dot{Y}}{Y}\,,\qquad C_{03}{}^{3}=\frac{\dot{Z}}{Z}\,.$
(38)

Similarly, we have $C_{\mu\nu 0}=0$, $C_{ijk}=0$, and the only nonzero
components of $C_{\mu\nu\rho}$ can be obtained from

$\displaystyle C_{011}=\dot{X}\,X\,,\qquad C_{022}=\dot{Y}\,Y\,,\qquad
C_{033}=\dot{Z}\,Z\,.$ (39)

It follows from these results that the torsion vector is given by

$\displaystyle
C_{0}=-\left(\frac{\dot{X}}{X}+\frac{\dot{Y}}{Y}+\frac{\dot{Z}}{Z}\right)\,,\qquad
C_{i}=0\,,$ (40)

while the torsion pseudovector $\check{C}_{\mu}=0$ in this case.

The calculations of contorsion (7) and the auxiliary torsion (9) tensors
produce similar results. That is, $K_{0\mu\nu}=0$, $K_{ijk}=0$, and the only
nonzero components of $K_{\mu\nu\rho}$ can be obtained from

$\displaystyle K_{101}=\dot{X}\,X\,,\qquad K_{202}=\dot{Y}\,Y\,,\qquad
K_{303}=\dot{Z}\,Z\,.$ (41)

Moreover, $\mathfrak{C}_{\mu\nu 0}=0$, $\mathfrak{C}_{ijk}=0$, and the only
nonzero components of $\mathfrak{C}_{\mu\nu\rho}$ can be obtained from

$\displaystyle\mathfrak{C}_{101}=X^{2}\,\left(\frac{\dot{Y}}{Y}+\frac{\dot{Z}}{Z}\right)\,,\quad\mathfrak{C}_{202}=Y^{2}\,\left(\frac{\dot{X}}{X}+\frac{\dot{Z}}{Z}\right)\,,\quad\mathfrak{C}_{303}=Z^{2}\,\left(\frac{\dot{X}}{X}+\frac{\dot{Y}}{Y}\right)\,.$
(42)

In the local limit of NLG, the constitutive relation of modified TEGR is given
by $N_{\mu\nu\rho}(x)=S(x)\mathfrak{C}_{\mu\nu\rho}(x)$, where the
gravitational susceptibility $S$ is a property of the background spacetime. In
the case of the homogeneous time-dependent background (26), we assume that $S$
is a function of time $t$. Therefore, $N_{\mu\nu 0}=0$, $N_{ijk}=0$ and the
only nonzero components of $N_{\mu\nu\rho}$ can be obtained from

$\displaystyle
N_{101}=S(t)\,X^{2}\,\left(\frac{\dot{Y}}{Y}+\frac{\dot{Z}}{Z}\right)\,,\quad
N_{202}=S(t)\,Y^{2}\,\left(\frac{\dot{X}}{X}+\frac{\dot{Z}}{Z}\right)\,,\quad
N_{303}=S(t)\,Z^{2}\,\left(\frac{\dot{X}}{X}+\frac{\dot{Y}}{Y}\right)\,.$ (43)

We can now compute $Q_{\mu\nu}$ given in Equation (14) and
$\mathcal{N}_{\mu\nu}$ given in Equation (18). The results are that these
quantities are diagonal with elements

$\displaystyle
Q_{00}=-S(t)\left(\frac{\dot{X}\dot{Y}}{XY}+\frac{\dot{Y}\dot{Z}}{YZ}+\frac{\dot{Z}\dot{X}}{ZX}\right)\,,$
(44) $\displaystyle Q_{11}=-S(t)X^{2}\,\frac{\dot{Y}\dot{Z}}{YZ}\,,\qquad
Q_{22}=-S(t)Y^{2}\,\frac{\dot{Z}\dot{X}}{ZX}\,,\qquad
Q_{33}=-S(t)Z^{2}\,\frac{\dot{X}\dot{Y}}{XY}\,.$ (45)

For $\mathcal{N}_{\mu\nu}$, however, we find $\mathcal{N}_{00}=0$ and

$\displaystyle\mathcal{N}_{11}$
$\displaystyle=-\frac{X^{2}}{YZ}\,\frac{d}{dt}[S(Y\dot{Z}+\dot{Y}Z)]\,,$ (46)
$\displaystyle\mathcal{N}_{22}$
$\displaystyle=-\frac{Y^{2}}{XZ}\,\frac{d}{dt}[S(X\dot{Z}+\dot{X}Z)]\,,$ (47)
$\displaystyle\mathcal{N}_{33}$
$\displaystyle=-\frac{Z^{2}}{XY}\,\frac{d}{dt}[S(X\dot{Y}+\dot{X}Y)]\,.$ (48)

Collecting everything, the modified GR field Equation (18) can be expressed as

$\displaystyle(1+S)\left(\frac{\dot{X}\dot{Y}}{XY}+\frac{\dot{Y}\dot{Z}}{YZ}+\frac{\dot{Z}\dot{X}}{ZX}\right)=\Lambda+8\pi
G\rho\,$ (49)

and

$\displaystyle(1+S)\left(\frac{\ddot{Y}}{Y}+\frac{\ddot{Z}}{Z}+\frac{\dot{Y}\dot{Z}}{YZ}\right)$
$\displaystyle=\Lambda-8\pi
GP-\frac{dS}{dt}\,\left(\frac{\dot{Y}}{Y}+\frac{\dot{Z}}{Z}\right)\,,$ (50)
$\displaystyle(1+S)\left(\frac{\ddot{X}}{X}+\frac{\ddot{Z}}{Z}+\frac{\dot{X}\dot{Z}}{XZ}\right)$
$\displaystyle=\Lambda-8\pi
GP-\frac{dS}{dt}\,\left(\frac{\dot{X}}{X}+\frac{\dot{Z}}{Z}\right)\,,$ (51)
$\displaystyle(1+S)\left(\frac{\ddot{X}}{X}+\frac{\ddot{Y}}{Y}+\frac{\dot{X}\dot{Y}}{XY}\right)$
$\displaystyle=\Lambda-8\pi
GP-\frac{dS}{dt}\,\left(\frac{\dot{X}}{X}+\frac{\dot{Y}}{Y}\right)\,.$ (52)

###  III.1 Field Equations

To express the gravitational field equations for the anisotropic models under
consideration in a more tractable form, it is useful to consider

$\displaystyle V(t)=XYZ\,,\qquad
W(t):=\frac{\dot{X}\dot{Y}}{XY}+\frac{\dot{Y}\dot{Z}}{YZ}+\frac{\dot{Z}\dot{X}}{ZX}\,,$
(53)

where $|V(t)|=\sqrt{-g}$ and note that

$\displaystyle\frac{\dot{V}}{V}=\frac{\dot{X}}{X}+\frac{\dot{Y}}{Y}+\frac{\dot{Z}}{Z}\,,\qquad\frac{\ddot{V}}{V}=\frac{\ddot{X}}{X}+\frac{\ddot{Y}}{Y}+\frac{\ddot{Z}}{Z}+2\,W\,.$
(54)

Let us add Equations (50)–(52) to get

$\displaystyle(1+S)\left(2\,\frac{\ddot{V}}{V}-3W\right)=3(\Lambda-8\pi
GP)-2\frac{dS}{dt}\,\frac{\dot{V}}{V}\,.$ (55)

Using Equation (49), we can write

$\displaystyle(1+S)\,\frac{\ddot{V}}{V}=3[\Lambda+4\pi
G(\rho-P)]-\frac{dS}{dt}\,\frac{\dot{V}}{V}\,$ (56)

or

$\displaystyle\frac{d}{dt}[(1+S)\,\dot{V}]=3V\,[\Lambda+4\pi G(\rho-P)]\,.$
(57)

Another interesting result is obtained by writing Equation (49) as
$VW=(\Lambda+8\pi G\rho)V/(1+S)$ and taking the time derivative of both sides.
From the relation

$\displaystyle\frac{1}{V}\,\frac{d(VW)}{dt}=\frac{\dot{X}}{X}\,\left(\frac{\ddot{Y}}{Y}+\frac{\ddot{Z}}{Z}+\frac{\dot{Y}\dot{Z}}{YZ}\right)+\frac{\dot{Y}}{Y}\,\left(\frac{\ddot{X}}{X}+\frac{\ddot{Z}}{Z}+\frac{\dot{X}\dot{Z}}{XZ}\right)+\frac{\dot{Z}}{Z}\,\left(\frac{\ddot{X}}{X}+\frac{\ddot{Y}}{Y}+\frac{\dot{X}\dot{Y}}{XY}\right)\,,$
(58)

we find

$\displaystyle(1+S)\frac{d(VW)}{dt}=(\Lambda-8\pi
GP)\dot{V}-2\,\frac{dS}{dt}VW\,.$ (59)

Using Equation (49), we finally get

$\displaystyle\frac{d\rho}{dt}=-(\rho+P)\frac{\dot{V}}{V}-\frac{\frac{dS}{dt}}{(1+S)}\,\left(\rho+\frac{\Lambda}{8\pi
G}\right)\,.$ (60)

Equations (57) and (60) are important consequences of the modified field
equations.

Furthermore, let us define a new temporal variable $\tau$ by

$\displaystyle\tau:=\int_{0}^{t}\frac{dt^{\prime}}{1+S(t^{\prime})}\,,\qquad
dt=(1+S)d\tau\,;$ (61)

hence, the spacetime metric in $(\tau,x,y,z)$ coordinates is

$\displaystyle
ds^{2}=-(1+S)^{2}d\tau^{2}+X^{2}(\tau)dx^{2}+Y^{2}(\tau)dy^{2}+Z^{2}(\tau)dz^{2}\,,$
(62)

where $S$ is now considered, by an abuse of notation, a function of $\tau$.
For instance, let us suppose $S(t)=t$; then, $\tau=\ln(1+t)$ and in the above
metric we have in this case $S=-1+e^{\tau}$. For metric (62), the Kretschmann
scalar $\mathcal{K}$ given by Equation (35) can be expressed in terms of the
new temporal variable $\tau$ using

$\displaystyle\frac{\dot{X}}{X}=(1+S)^{-1}\frac{1}{X}\,\frac{dX}{d\tau}\,$
(63)

and

$\displaystyle\frac{\ddot{X}}{X}=(1+S)^{-2}\,\left(\frac{1}{X}\,\frac{d^{2}X}{d\tau^{2}}-\frac{\frac{dS}{d\tau}}{1+S}\,\frac{1}{X}\,\frac{dX}{d\tau}\right)\,.$
(64)

The gravitational field equations can now be written in terms of the temporal
variable $\tau$ as

$\displaystyle\frac{1}{XY}\,\frac{dX}{d\tau}\,\frac{dY}{d\tau}+\frac{1}{YZ}\,\frac{dY}{d\tau}\,\frac{dZ}{d\tau}+\frac{1}{ZX}\,\frac{dZ}{d\tau}\,\frac{dX}{d\tau}=(1+S)(\Lambda+8\pi
G\rho)\,,$ (65)
$\displaystyle\frac{1}{Y}\,\frac{d^{2}Y}{d\tau^{2}}+\frac{1}{Z}\,\frac{d^{2}Z}{d\tau^{2}}+\frac{1}{YZ}\,\frac{dY}{d\tau}\,\frac{dZ}{d\tau}=(1+S)(\Lambda-8\pi
GP)\,,$ (66)
$\displaystyle\frac{1}{X}\,\frac{d^{2}X}{d\tau^{2}}+\frac{1}{Z}\,\frac{d^{2}Z}{d\tau^{2}}+\frac{1}{XZ}\,\frac{dX}{d\tau}\,\frac{dZ}{d\tau}=(1+S)(\Lambda-8\pi
GP)\,,$ (67)
$\displaystyle\frac{1}{X}\,\frac{d^{2}X}{d\tau^{2}}+\frac{1}{Y}\,\frac{d^{2}Y}{d\tau^{2}}+\frac{1}{XY}\,\frac{dX}{d\tau}\,\frac{dY}{d\tau}=(1+S)(\Lambda-8\pi
GP)\,.$ (68)

To solve Equations (66)–(68), let us subtract, for instance, Equation (66)
from Equation (67) to get

$\displaystyle\frac{1}{X}\,\frac{d^{2}X}{d\tau^{2}}-\frac{1}{Y}\,\frac{d^{2}Y}{d\tau^{2}}+\frac{1}{Z}\,\frac{dZ}{d\tau}\left(\frac{1}{X}\,\frac{dX}{d\tau}-\frac{1}{Y}\,\frac{dY}{d\tau}\right)=0\,,$
(69)

which with $V=XYZ$ can be written as

$\displaystyle\frac{d}{d\tau}\left(\frac{1}{X}\,\frac{dX}{d\tau}-\frac{1}{Y}\,\frac{dY}{d\tau}\right)+\frac{1}{V}\,\frac{dV}{d\tau}\,\left(\frac{1}{X}\,\frac{dX}{d\tau}-\frac{1}{Y}\,\frac{dY}{d\tau}\right)=0\,.$
(70)

Hence,

$\displaystyle
V\,\left(\frac{1}{X}\,\frac{dX}{d\tau}-\frac{1}{Y}\,\frac{dY}{d\tau}\right)=\mathcal{C}_{12}\,,$
(71)

where $\mathcal{C}_{12}$ is a constant of integration and similar results hold
for the other metric functions.

Finally, in terms of temporal variable $\tau$, Equation (57) can be written as

$\displaystyle\frac{1}{V}\,\frac{d^{2}V}{d\tau^{2}}=3(1+S)[\Lambda+4\pi
G(\rho-P)]\,.$ (72)

###  III.2 Special Solutions

We now explore some cases of particular interest.

####  III.2.1 de Sitter

Let us first consider de Sitter’s solution with

$\displaystyle\rho=P=0\,,\qquad X=Y=Z=e^{\lambda t}\,,$ (73)

where $\lambda$ is a nonzero constant. The field equations imply

$\displaystyle 3\lambda^{2}(1+S)=\Lambda\,,\qquad
3\lambda^{2}(1+S)=\Lambda-2\lambda\frac{dS}{dt}\,.$ (74)

Therefore, $dS/dt=0$ and $S$ must be constant. It follows that de Sitter
spacetime is not a solution of the modified TEGR since the susceptibility
function is independent of time while the background spacetime is dynamic.
This is in agreement with the fact that de Sitter spacetime is not a solution
of NLG Mashhoon:2022ynk .

####  III.2.2 Kasner

Let us next consider the standard Kasner metric Kas ; Chicone:2011ie

$ds^{2}=-dt^{2}+t^{2p_{1}}dx^{2}+t^{2p_{2}}dy^{2}+t^{2p_{3}}dz^{2}\,,$ (75)
$p_{1}+p_{2}+p_{3}=p_{1}^{2}+p_{2}^{2}+p_{3}^{2}=1\,.$ (76)

In GR, this empty universe model is a solution of the field equations with
$\Lambda=0$ and $\rho=P=0$. Note that with $p_{1}=1$, say, and hence
$p_{2}=p_{3}=0$, we recover flat spacetime. Therefore, we can assume
$p_{1}<p_{2}<p_{3}$; that is,

$-\frac{1}{3}\leq p_{1}\leq 0,\qquad 0\leq
p_{2}\leq\frac{2}{3},\qquad\frac{2}{3}\leq p_{3}\leq 1\,.$ (77)

It follows from field equations (49) and (50) that

$\frac{(1+S)}{t^{2}}(p_{1}p_{2}+p_{2}p_{3}+p_{1}p_{3})=\Lambda+8\pi G\rho\,,$
(78)
$\frac{(1+S)}{t^{2}}(p_{2}^{2}-p_{2}+p_{3}^{2}-p_{3}+p_{2}p_{3})=\Lambda-8\pi
GP-\frac{dS}{dt}(p_{2}+p_{3})\,,$ (79)

etc. Because $p_{1}p_{2}+p_{2}p_{3}+p_{1}p_{3}=0$ and
$p_{2}^{2}-p_{2}+p_{3}^{2}-p_{3}+p_{2}p_{3}=0$ together with two other
cyclically related terms that vanish, we find that $\rho$, $P$, and $S$ must
be constants such that

$-8\pi G\rho=\Lambda\leq 0\,,\qquad\rho+P=0\,,\qquad S={\rm constant}\,.$ (80)

Therefore, Kasner’s spacetime with constant $S$ is not a solution of the
modified TEGR.

####  III.2.3 Flat FLRW Model

We now consider $X=Y=Z=a(t)$. Then, field equations (49) and (50) imply

$3(1+S)\left(\frac{\dot{a}}{a}\right)^{2}=\Lambda+8\pi G\rho\,,$ (81)
$(1+S)\left[2\frac{\ddot{a}}{a}+\left(\frac{\dot{a}}{a}\right)^{2}\right]=\Lambda-8\pi
GP-2\frac{dS}{dt}\,\frac{\dot{a}}{a}\,.$ (82)

A thorough treatment of these equations is contained in a recent paper
Tabatabaei:2022tbq , where they were employed with $\Lambda=0$ in a detailed
discussion of the implications of the modified Cartesian flat cosmology in
connection with $H_{0}$ tension.

It is important to emphasize that in these time-dependent solutions considered
thus far, $S$ must be dependent upon time as well; otherwise, we do not have a
physically meaningful solution of the theory.

##  IV Solution for Dust with $\Lambda=0$

###  IV.1 $S=0$

When $S=0$, we get the gravitational field equations in GR. A well-known class
of anisotropic solutions corresponds to $\Lambda=0$ and $P=0$ HS1 ; HS2 ; HE ;
PK . In this case, Equations (57) and (60) imply

$\displaystyle\rho\,V=\frac{\ell_{0}}{6\pi G}\,,\qquad
V=\ell_{0}\,t^{2}(1+\Sigma/t)\,,$ (83)

where $\ell_{0}>0$ and $\Sigma>0$ are integration constants with dimensions of
length (or time, since $c=1$) and we have assumed that $V=0$ at $t=0$. Using

$\displaystyle\int\frac{dt}{V}=-\frac{1}{\ell_{0}\Sigma}\ln(1+\Sigma/t)\,,$
(84)

Equation (71) and similar ones can be integrated with the result that ratios
such as $X/Y$, etc., up to constant coefficients are given by $(1+\Sigma/t)$
to some constant powers. In the absence of anisotropy ($\Sigma=0$), we must
recover the standard Einstein-de Sitter solution; that is, $\Sigma$ is the
anisotropy parameter for finite $t>0$. Therefore, we look for the solutions of
the field equations such that the metric coefficients are given by

$\displaystyle X=\ell_{0}{}^{q_{1}}\,t^{2/3}(1+\Sigma/t)^{q_{1}}\,,\quad
Y=\ell_{0}{}^{q_{2}}\,t^{2/3}(1+\Sigma/t)^{q_{2}}\,,\quad
Z=\ell_{0}{}^{q_{3}}\,t^{2/3}(1+\Sigma/t)^{q_{3}}\,,$ (85)

where $q_{i}$, $i=1,2,3$, are constants that must add up to unity in order to
satisfy Equation (83). A detailed examination reveals that the field equations
are all satisfied in the case of dust with $\Lambda=0$ provided

$\displaystyle q_{1}+q_{2}+q_{3}=1\,,\qquad
q_{1}^{2}+q_{2}^{2}+q_{3}^{2}=1\,.$ (86)

It is possible to choose $q_{1}<q_{2}<q_{3}$; that is, $-\frac{1}{3}\leq
q_{1}\leq 0$, $0\leq q_{2}\leq\frac{2}{3}$, and $\frac{2}{3}\leq q_{3}\leq 1$.
Then, a convenient representation of these constants is given by

$\displaystyle q_{i}=\frac{1}{3}-\frac{2}{3}\sin[\theta-2(i-1)\pi/3]\,,\qquad
i=1,2,3\,,\qquad\frac{\pi}{6}\leq\theta\leq\frac{\pi}{2}\,.$ (87)

It is important to note that at late times, $t\to\infty$, we recover the
standard Einstein-de Sitter solution.

###  IV.2 $S\neq 0$

Let us note that a constant $S\neq 0$ is equivalent to a constant rescaling of
the time coordinate in GR, namely, $t\to(1+S)^{1/2}t$. Therefore, new results
are obtained only when $S$ depends upon time, which is necessary for a proper
physical interpretation of the susceptibility function $S$.

We now consider the case $dS/dt\neq 0$. For $P=0$ and $\Lambda=0$, it follows
from Equations (60) and (72) that

$\displaystyle(1+S)\rho\,V=\frac{\ell_{0}}{6\pi G}\,,\qquad
V=\ell_{0}\,\tau^{2}(1+\Sigma/\tau)\,.$ (88)

Moreover, the metric coefficients are given by

$\displaystyle X=\ell_{0}{}^{q_{1}}\,\tau^{2/3}(1+\Sigma/\tau)^{q_{1}}\,,\quad
Y=\ell_{0}{}^{q_{2}}\,\tau^{2/3}(1+\Sigma/\tau)^{q_{2}}\,,\quad
Z=\ell_{0}{}^{q_{3}}\,\tau^{2/3}(1+\Sigma/\tau)^{q_{3}}\,,$ (89)

where $q_{1}+q_{2}+q_{3}=1$ and $q_{1}^{2}+q_{2}^{2}+q_{3}^{2}=1$. For
$\Sigma=0$ or $\tau\to\infty$, we recover the isotropic solution corresponding
to the modified Einstein-de Sitter model.

It is interesting to explore the implications of these solutions for the
dimensionless deceleration parameter defined for the $x$ direction, say, by
$Q_{x}=-(\ddot{X}/X)/(\dot{X}/X)^{2}$, etc.

Let us first recall that $\Sigma>0$ is the anisotropy parameter; in fact, for
$\Sigma=0$, we get the standard Einstein-de Sitter result that $Q=1/2$ in
every direction. Let anisotropic expansion occur in the $z$ direction; then,
Equation (86) allows two possibilities, namely, $(q_{1},q_{2},q_{3})=(0,0,1)$
and $(q_{1},q_{2},q_{3})=(2/3,2/3,-1/3)$. We choose the former case for the
sake of simplicity. We are interested in the nature of deceleration parameter
$Q_{z}$ in the $z$ direction. For $S=0$, we have $Q_{x}=Q_{y}=1/2$ and

$\displaystyle
Q_{z}=2\frac{(t+\Sigma)(t-2\Sigma)}{(2t-\Sigma)^{2}}<\frac{1}{2}\,.$ (90)

Indeed, for $t<2\,\Sigma$ we have $Q_{z}<0$ and hence acceleration in the $z$
direction that later turns into deceleration for $t>2\,\Sigma$ with magnitude
$<1/2$. As $t\to\infty$, $Q_{z}\to 1/2$; that is, isotropy is recovered at
late times.

How do these results change in the presence of $S(t)$? In terms of temporal
parameter $\tau$, we find for the $x$ direction, say,

$\displaystyle
Q_{x}=-X\frac{d^{2}X}{d\tau^{2}}\left(\frac{dX}{d\tau}\right)^{-2}+\frac{dS}{dt}X\left(\frac{dX}{d\tau}\right)^{-1}\,.$
(91)

Therefore, with

$\displaystyle X=Y=\tau^{2/3}\,,\qquad
Z=\ell_{0}\,\tau^{2/3}(1+\Sigma/\tau)\,,$ (92)

the results are

$\displaystyle Q_{x}=Q_{y}=\frac{1}{2}+\frac{3}{2}\tau\,\frac{dS}{dt}\,$ (93)

and

$\displaystyle
Q_{z}=2\frac{(\tau+\Sigma)(\tau-2\Sigma)}{(2\tau-\Sigma)^{2}}+3\left(\frac{\tau+\Sigma}{2\tau-\Sigma}\right)\tau\,\frac{dS}{dt}\,.$
(94)

Let us suppose, for instance, that $dS/dt>0$. Then, the deceleration increases
in the $x$ and $y$ directions. The same is true in the $z$ direction for
$\tau>2\,\Sigma$; however, for $\tau<\Sigma/2$, the presence of $dS/dt>0$
causes extra acceleration in the $z$ direction. For $2\,\Sigma>\tau>\Sigma/2$,
the sign of $Q_{z}$ depends upon the magnitude of $dS/dt$. We note that
isotropy is restored at late times $\tau\gg\Sigma$.

##  V Solution for Dynamic Dark Energy with $\Lambda=0$

Let us imagine a universe that in the absence of the cosmological constant
($\Lambda=0$) is dominated by dynamic dark energy with $P_{de}+\rho_{de}=0$.
Therefore, Equation (60) implies

$\displaystyle\rho_{de}=\rho_{de}(t_{0})\,\frac{1+S_{0}}{1+S}\,,\qquad\rho_{de}(t_{0})>0\,,\quad
S_{0}=S(t_{0})\,,$ (95)

where $t_{0}$ refers to some fiducial epoch in the expansion of the universe.
The pressure of dark energy is always negative, $P_{de}=-\rho_{de}$, and we
assume on the basis of Equation (95) that

$\displaystyle 6\pi G\rho_{de}(1+S)=6\pi
G\rho_{de}(t_{0})\,(1+S_{0}):=\frac{1}{\tau_{0}^{2}}\,,\qquad\eta:=\frac{\tau}{\tau_{0}}\,,$
(96)

where $\tau_{0}>0$ is a constant with the dimensions of time and henceforth we
employ $\eta$ as the new dimensionless temporal variable. To explore the
anisotropic acceleration of this dark energy universe model, we assume for the
sake of simplicity that $X=Y$. Thus we have a cylindrical model with the $z$
direction as the main direction of anisotropy. Let us note here that for $S=0$
our dynamic dark energy source in effect reduces to a cosmological constant.

With $X=Y$, Equation (68) reduces to

$\displaystyle\frac{2}{X}\,\frac{d^{2}X}{d\eta^{2}}+\frac{1}{X^{2}}\,\left(\frac{dX}{d\eta}\right)^{2}=\frac{4}{3}\,.$
(97)

This nonlinear equation can be easily solved with

$\displaystyle
X:=\alpha^{2/3}\,,\qquad\frac{d^{2}\alpha}{d\eta^{2}}-\alpha=0\,,\qquad\alpha=C_{+}\,e^{\eta}+C_{-}\,e^{-\eta}\,,$
(98)

where $C_{\pm}$ are integration constants. Next, Equation (65) implies

$\displaystyle\frac{1}{X^{2}}\,\left(\frac{dX}{d\eta}\right)^{2}+\frac{2}{Z}\,\frac{dZ}{d\eta}\frac{1}{X}\frac{dX}{d\eta}=\frac{4}{3}\,,$
(99)

or

$\displaystyle\frac{1}{Z}\,\frac{dZ}{d\eta}=\frac{\alpha}{\beta}-\frac{1}{3}\frac{\beta}{\alpha}\,,$
(100)

where $\alpha=d\beta/d\eta$ and

$\displaystyle\beta=\frac{d\alpha}{d\eta}=C_{+}\,e^{\eta}-C_{-}\,e^{-\eta}\,.$
(101)

Hence, $Z=C_{0}\,\alpha^{-1/3}\,\beta$, where $C_{0}$ is an integration
constant. Therefore,

$\displaystyle X=Y=\alpha^{2/3}\,,\quad Z=C_{0}\,\alpha^{-1/3}\,\beta\,,\quad
V=X^{2}Z=C_{0}\,\alpha\,\beta=C_{0}(C^{2}_{+}\,e^{2\eta}-C^{2}_{-}\,e^{-2\eta})\,.$
(102)

One can check that this is the general solution of the field equations in the
present case. For these solutions, the Kretschmann scalar $\mathcal{K}$ can be
written as

$\displaystyle\frac{1}{4}(1+S)^{4}\,\tau_{0}^{4}\,\mathcal{K}=\frac{16}{9}K_{1}-\frac{8}{3}\frac{\beta}{\alpha}\left(1-\frac{\beta^{2}}{3\alpha^{2}}\right)\tau_{0}\,\frac{dS}{dt}+K_{2}\,\frac{\alpha^{2}}{\beta^{2}}\left(\tau_{0}\,\frac{dS}{dt}\right)^{2}\,,$
(103)

where

$\displaystyle
K_{1}=1-\frac{2}{3}\frac{\beta^{2}}{\alpha^{2}}+\frac{1}{3}\frac{\beta^{4}}{\alpha^{4}}\,,\qquad
K_{2}=1-\frac{2}{3}\frac{\beta^{2}}{\alpha^{2}}+\frac{\beta^{4}}{\alpha^{4}}\,.$
(104)

We are interested in the behavior of the acceleration of these anisotropic
solutions.

For the deceleration parameters of this model, we find

$\displaystyle
Q_{x}=Q_{y}=\frac{1}{2}\left(1-3\frac{\alpha^{2}}{\beta^{2}}\right)+\frac{3}{2}\tau_{0}\frac{dS}{dt}\frac{\alpha}{\beta}\,,$
(105) $\displaystyle
Q_{z}=-4\left(1-3\frac{\alpha^{2}}{\beta^{2}}\right)^{-2}-3\tau_{0}\frac{dS}{dt}\frac{\alpha}{\beta}\left(1-3\frac{\alpha^{2}}{\beta^{2}}\right)^{-1}\,.$
(106)

At late times $\eta\to\infty$, $\alpha/\beta\to 1$, and these expressions
reduce to $Q_{x}=Q_{y}=Q_{z}=-1+(3\tau_{0}/2)\,dS/dt$. However, for finite
$\eta$ there could be anomalous behavior when $\beta=0$ or
$3\alpha^{2}=\beta^{2}$. If $\beta=0$, then $Q_{x}=Q_{y}=-\infty$, while
$Q_{z}=0$. On the other hand, if $3\alpha^{2}=\beta^{2}$, then
$Q_{x}=Q_{y}=\pm(\sqrt{3}\tau_{0}/2)dS/dt$, while $Q_{z}=-\infty$. To give an
example of the former situation, consider, for instance, the case where
$C_{+}=C_{-}=C>0$ with $\alpha=2C\,\cosh\eta$ and $\beta=2C\,\sinh\eta$. In
this case, $3\alpha^{2}>\beta^{2}$. Then, $X=Y=(2C)^{2/3}\cosh^{2/3}\eta$ and
$Z=C_{0}(2C)^{2/3}\cosh^{-1/3}\eta\,\sinh\eta$. The universe starts from a
singular pancake state at $\eta=0$ and expands with infinite acceleration in
the $x$ and $y$ directions but with zero acceleration in the $z$ direction.
Indeed, inspection of the Kretschmann scalar $\mathcal{K}$ given by Equation
(103) reveals that $\mathcal{K}$ diverges at $\tau=0$ provided $dS/d\tau\neq
0$.

For an example of the case where $3\alpha^{2}=\beta^{2}$, let us consider, for
instance, $C_{+}=1$ and $C_{-}=-(2+\sqrt{3})$. Then, the universe starts from
a cigar state at $\eta=0$ with infinite acceleration in the $z$ direction but
with finite acceleration $-(\sqrt{3}\tau_{0}/2)dS/dt$ in the $x$ and $y$
directions assuming that $dS/dt>0$. The Kretschmann scalar (103) remains
finite in this case.

These results could be interesting in connection with recent observational
evidence in favor of anomalous anisotropic acceleration of bulk flow in the
local universe Colin:2019opb ; Secrest:2020has ; Secrest:2022uvx ;
Solanki:2023yoa ; Perivolaropoulos:2023tdt .

##  VI DISCUSSION

The cosmic acceleration in the standard cosmological model is assumed to be
isotropic. Large-scale anisotropy in the Hubble flow has been the subject of
numerous studies. We have explored anisotropy in cosmic acceleration within
the theoretical framework of a recent teleparallel extension of general
relativity that corresponds to the local limit of nonlocal gravity. This
modified theory involves a function $S(x)$ such that GR is recovered for
$S=0$. We are interested in the cosmological significance of the extra
function $S(x)$. In particular, in a dynamic Bianchi type I model that is
consistent with the modified GR field equations, we theoretically investigate
anisotropy in cosmic acceleration and determine the contribution of $S(t)$ to
local anisotropic acceleration. Our results in Sections IV and V for possible
local anisotropic cosmic acceleration depend explicitly upon $dS(t)/dt$.

The function $S(t)$ may also possibly contribute to the resolution of an
anomaly in the quadrupole anisotropy of CMB temperature. It has been reported
that the CMB temperature angular power spectrum suffers from a deficit in its
quadrupole moment WMAP:2012fli ; wmap ; Planck:2019evm ; Planck:2018vyg . The
anomaly has to do with the low amplitude of the quadrupole anisotropy compared
to the prediction of the standard $\Lambda$CDM model. On the other hand, in an
anisotropic universe described by metric (26), one would in general expect a
change in the quadrupole moment of the CMB temperature due to different
amounts of redshift suffered by photons traveling after recombination along
different directions toward the observer Akarsu:2019pwn ; Campanelli:2006vb ;
Campanelli:2007qn ; Buniy:2005qm ; Sachs:1967er . However, Big Bang
nucleosynthesis (BBN) puts a tight constraint on the anisotropy of the
expansion rate Barrow:1976rda ; hence, it is rather difficult for the
anisotropic expansion to justify the CMB quadrupole deficit Akarsu:2019pwn . A
suitable susceptibility function $S(t)$ of the modified gravity theory may be
able to ameliorate the situation. That is, $S(t)$ may allow the quadrupole
deficit anomaly be alleviated while respecting the BBN constraint. A more
complete discussion is beyond the scope of the present paper.

## ACKNOWLEDGMENTS

The work of AB has been supported financially by Iran Science Elites
Federation.

## References

  * (1) Einstein, A. _The Meaning of Relativity_ ; Princeton University Press: Princeton, NJ, USA, 1955.
  * (2) Hehl, F.W.; Mashhoon, B. Nonlocal Gravity Simulates Dark Matter. _Phys. Lett. B_ 2009, _673_ , 279–282. [arXiv:0812.1059 [gr-qc]]
  * (3) Hehl, F.W.; Mashhoon, B. Formal framework for a nonlocal generalization of Einstein’s theory of gravitation. _Phys. Rev. D_ 2009, _79_ , 064028. [arXiv:0902.0560 [gr-qc]]
  * (4) Mashhoon, B. _Nonlocal Gravity_ ; Oxford University Press: Oxford, UK, 2017.
  * (5) Hopkinson, J. Residual Charge of the Leyden Jar.—Dielectric Properties of different Glasses. _Phil. Trans. Roy. Soc. London_ 1877, _167_ , 599–626.
  * (6) Poisson, S. Mémoire sur la théorie du magnétisme en mouvement. _Mém. Acad. Sci. France_ 1823, _6_ , 441–570.
  * (7) Jackson, J.D. _Classical Electrodynamics_ , 3rd edn; Wiley: Hoboken, NJ, USA, 1999.
  * (8) Landau, L.D.; Lifshitz, E.M. _Electrodynamics of Continuous Media_ ; Pergamon: Oxford, UK, 1960.
  * (9) Hehl, F.W.; Obukhov, Y.N. _Foundations of Classical Electrodynamics: Charge, Flux, and Metric_ ; Birkhäuser: Boston, MA, USA, 2003.
  * (10) van den Hoogen, R.J. Towards a covariant smoothing procedure for gravitational theories. _J. Math. Phys._ 2017, _58_ , no.12, 122501.
  * (11) Maggiore, M.; Mancarella, M. Nonlocal gravity and dark energy. _Phys. Rev. D_ 2014, _90_ , no.2, 023005. [arXiv:1402.0448 [hep-th]]
  * (12) Capozziello, S.; Bajardi, F. Non-Local Gravity Cosmology: an Overview. _Int. J. Mod. Phys. D_ 2022, _31_ , no.06, 2230009. [arXiv:2201.04512 [gr-qc]]
  * (13) Kumar, K.S.; Maheshwari, S.; Mazumdar, A.; Peng, J. An anisotropic bouncing universe in non-local gravity. _JCAP_ 2021, _07_ , 025. [arXiv:2103.13980 [gr-qc]]
  * (14) Dimitrijevic, I.; Dragovich, B.; Koshelev, A.S.; Rakic, Z.; Stankovic, J. Cosmological Solutions of a Nonlocal Square Root Gravity. _Phys. Lett. B_ 2019, _797_ , 134848. [arXiv:1906.07560 [gr-qc]]
  * (15) Koshelev, A.S.; Kumar, K.S.; Starobinsky, A.A. Cosmology in nonlocal gravity, [arXiv:2305.18716 [hep-th]].
  * (16) Deser, S.; Woodard, R.P. Nonlocal Cosmology II — Cosmic acceleration without fine tuning or dark energy. _JCAP_ 2019, _06_ , 034. [arXiv:1902.08075 [gr-qc]]
  * (17) Chen, C.Y.; Chen, P.; Park, S. Primordial bouncing cosmology in the Deser-Woodard nonlocal gravity. _Phys. Lett. B_ 2019, _796_ , 112–116. [arXiv:1905.04557 [gr-qc]]
  * (18) Weitzenböck, R. Invariantentheorie; Noordhoff: Groningen, The Netherlands, 1923.
  * (19) Maluf, J.W. The teleparallel equivalent of general relativity. _Ann. Phys. (Berlin)_ 2013, _525_ , 339–357. [arXiv:1303.3897 [gr-qc]]
  * (20) Aldrovandi, R.; Pereira, J.G. _Teleparallel Gravity: An Introduction_ ; Springer: New York, 2013.
  * (21) Blagojević, M.; Hehl, F.W., Eds. _Gauge Theories of Gravitation_ ; Imperial College Press: London, UK, 2013.
  * (22) Itin, Y.; Obukhov, Y.N.; Boos, J.; Hehl, F.W. Premetric teleparallel theory of gravity and its local and linear constitutive law. _Eur. Phys. J. C_ 2018, _78_ , 907. [arXiv:1808.08048 [gr-qc]]
  * (23) Cho, Y.M. Einstein Lagrangian as the translational Yang-Mills Lagrangian. _Phys. Rev. D_ 1976, _14_ , 2521–2525.
  * (24) Puetzfeld, D.; Obukhov, Y.N.; Hehl, F.W. Constitutive law of nonlocal gravity. _Phys. Rev. D_ 2019, _99_ , no. 10, 104013. [arXiv:1903.04023 [gr-qc]]
  * (25) Mashhoon, B. Nonlocal Gravity: Fundamental Tetrads and Constitutive Relations. _Symmetry_ 2022, _14_ , 2116. [arXiv:2209.05817 [gr-qc]]
  * (26) Bini, D.; Mashhoon, B. Nonlocal gravity: Conformally flat spacetimes. _Int. J. Geom. Meth. Mod. Phys._ 2016, _13_ , 1650081. [arXiv:1603.09477 [gr-qc]]
  * (27) Rahvar, S.; Mashhoon, B. Observational Tests of Nonlocal Gravity: Galaxy Rotation Curves and Clusters of Galaxies. _Phys. Rev. D_ 2014, _89_ , 104011. [arXiv:1401.4819 [astro-ph.GA]]
  * (28) Chicone, C.; Mashhoon, B. Nonlocal Gravity in the Solar System. _Classical Quantum Gravity_ 2016, _33_ , 075005. [arXiv:1508.01508 [gr-qc]]
  * (29) Roshan, M.; Mashhoon, B. Dynamical Friction in Nonlocal Gravity. _Astrophys. J._ 2021, _922_ , 9. [arXiv:2107.05841 [gr-qc]]
  * (30) Roshan, M.; Mashhoon, B. Characteristics of Effective Dark Matter in Nonlocal Gravity. _Astrophys. J._ 2022, _934_ , 9. [arXiv:2201.12852 [astro-ph.GA]]
  * (31) Roshan, M.; Mashhoon, B. Nonlocal Gravity: Modification of Newtonian Gravitational Force in the Solar System. _Universe_ 2022, _8_ , 470. [arXiv:2205.13276 [gr-qc]]
  * (32) Tabatabaei, J.; Baghram, S.; Mashhoon, B. Local Limit of Nonlocal Gravity: A Teleparallel Extension of General Relativity, [arXiv:2212.05536 [gr-qc]].
  * (33) Tabatabaei, J.; Banihashemi, A.; Baghram, S.; Mashhoon, B. Dynamic Dark Energy from the Local Limit of Nonlocal Gravity. _Int. J. Mod. Phys. D_ 2023, in press. [arXiv:2305.07630 [gr-qc]]
  * (34) Colin, J.; Mohayaee, R.; Rameez, M.; Sarkar, S. Evidence for anisotropy of cosmic acceleration. _Astron. Astrophys._ 2019, _631_ , L13. [arXiv:1808.04597 [astro-ph.CO]]
  * (35) Secrest, N.J.; von Hausegger, S.; Rameez, M.; Mohayaee, R.; Sarkar, S.; Colin, J. A Test of the Cosmological Principle with Quasars. _Astrophys. J. Lett._ 2021, _908_ , no.2, L51. [arXiv:2009.14826 [astro-ph.CO]]
  * (36) Secrest, N.J.; von Hausegger, S.; Rameez, M.; Mohayaee, R.; Sarkar, S. A Challenge to the Standard Cosmological Model. _Astrophys. J. Lett._ 2022, _937_ , no.2, L31. [arXiv:2206.05624 [astro-ph.CO]]
  * (37) Solanki, R.; Patel, B.; Jaybhaye, L.V.; Sahoo, P.K. Cosmic acceleration with bulk viscosity in an anisotropic $f(R,L_{m})$ background. _Commun. Theor. Phys._ 2023, _75_ , 075401. [arXiv:2305.07683 [gr-qc]]
  * (38) Perivolaropoulos, L. On the isotropy of SnIa absolute magnitudes in the Pantheon+ and SH0ES samples. [arXiv:2305.12819 [astro-ph.CO]]
  * (39) Stephani, H.; Kramer, D.; MacCallum, M.A.H.; Hoenselaers, C.; Herlt, E. _Exact Solutions of Einstein’s Field Equations_ , 2nd edn; Cambridge University Press: Cambridge, UK, 2003.
  * (40) Coley, A.A.; van den Hoogen, R.J. Spatially Homogeneous Teleparallel Gravity: Bianchi I. [arXiv:2305.12168 [gr-qc]]
  * (41) Heckmann, O.; Schücking, E. World models. In _La Structure et l’evolution de l’univers_ , 149–162; Editions Stoops: Brussels, 1958.
  * (42) Heckmann, O.; Schücking, E. Relativistic Cosmology. In _Gravitation: an Introduction to Current Research_ , ed. L. Witten, 438–469; Wiley: New York, 1962.
  * (43) Hawking, S.W.; Ellis, G.F.R. _The Large Scale Structure of Space-Time_ ; Cambridge University Press: Cambridge, UK, 1973.
  * (44) Plebanski, J.; Krasinski, A. _An Introduction to General Relativity and Cosmology_ ; Cambridge University Press: Cambridge, UK, 2006.
  * (45) Kasner, E. Geometrical Theorems on Einstein’s Cosmological Equations. _Am. J. Math._ 1921, _43_ , 217–221.
  * (46) Chicone, C.; Mashhoon, B.; Rosquist, K. Double-Kasner Spacetime: Peculiar Velocities and Cosmic Jets. _Phys. Rev. D_ 2011, _83_ , 124029. [arXiv:1104.5058 [gr-qc]]
  * (47) Akrami, Y. et al. [Planck Collaboration]. Planck 2018 results. VII. Isotropy and Statistics of the CMB. _Astron. Astrophys._ 2020, _641_ , A7. [arXiv:1906.02552 [astro-ph.CO]]
  * (48) Aghanim, N. et al. [Planck Collaboration]. Planck 2018 results. VI. Cosmological parameters. _Astron. Astrophys._ 2020, _641_ , A6. [erratum: _Astron. Astrophys._ 2021, _652_ , C4.] [arXiv:1807.06209 [astro-ph.CO]]
  * (49) Bennett, C.L. et al. [WMAP Collaboration]. Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Are There Cosmic Microwave Background Anomalies? _Astrophys. J. Suppl._ 2011, _192_ , 17. [arXiv:1001.4758 [astro-ph.CO]]
  * (50) Bennett, C.L. et al. [WMAP Collaboration]. Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Final Maps and Results. _Astrophys. J. Suppl._ 2013, _208_ , 20. [arXiv:1212.5225 [astro-ph.CO]]
  * (51) Akarsu, Ö.; Kumar, S.; Sharma, S.; Tedesco, L. Constraints on a Bianchi type I spacetime extension of the standard $\Lambda$CDM model. _Phys. Rev. D_ 2019, _100_ , no.2, 023532. [arXiv:1905.06949 [astro-ph.CO]]
  * (52) Campanelli, L.; Cea, P.; Tedesco, L. Ellipsoidal Universe Can Solve The CMB Quadrupole Problem. _Phys. Rev. Lett._ 2006, _97_ , 131302. [erratum: _Phys. Rev. Lett._ 2006, _97_ , 209903.] [arXiv:astro-ph/0606266 [astro-ph]]
  * (53) Campanelli, L.; Cea, P.; Tedesco, L. Cosmic Microwave Background Quadrupole and Ellipsoidal Universe. _Phys. Rev. D_ 2007, _76_ , 063007. [arXiv:0706.3802 [astro-ph]]
  * (54) Buniy, R.V.; Berera, A.; Kephart, T.W. Asymmetric inflation: Exact solutions. _Phys. Rev. D_ 2006, _73_ , 063529. [arXiv:hep-th/0511115 [hep-th]]
  * (55) Sachs, R.K.; Wolfe, A.M. Perturbations of a cosmological model and angular variations of the microwave background. _Astrophys. J._ 1967, _147_ , 73–90.
  * (56) Barrow, J. Light elements and the isotropy of the Universe. _Mon. Not. Roy. Astron. Soc._ 1976, _175_ , 359–370.