# Inertial Effects on the Stress Generation of Active Fluids

S.C. Takatori111Present address: Miller Institute for Basic Research in
Science, University of California, Berkeley, 2536 Channing Way, D-104,
Berkeley, California 94720, USA; Department of Bioengineering, University of
California, Berkeley, California 94720, USA and J.F. Brady  Division of
Chemistry and Chemical Engineering, California Institute of Technology,
Pasadena, California 91125, USA

###### Abstract

Suspensions of self-propelled bodies generate a unique mechanical stress owing
to their motility that impacts their large-scale collective behavior. For
microswimmers suspended in a fluid with negligible particle inertia, we have
shown that the virial ‘swim stress’ is a useful quantity to understand the
rheology and nonequilibrium behaviors of active soft matter systems. For
larger self-propelled organisms like fish, it is unclear how particle inertia
impacts their stress generation and collective movement. Here, we analyze the
effects of finite particle inertia on the mechanical pressure (or stress)
generated by a suspension of self-propelled bodies. We find that swimmers of
all scales generate a unique ‘swim stress’ and ‘Reynolds stress’ that impacts
their collective motion. We discover that particle inertia plays a similar
role as confinement in overdamped active Brownian systems, where the reduced
run length of the swimmers decreases the swim stress and affects the phase
behavior. Although the swim and Reynolds stresses vary individually with the
magnitude of particle inertia, the sum of the two contributions is independent
of particle inertia. This points to an important concept when computing
stresses in computer simulations of nonequilibrium systems—the Reynolds and
the virial stresses must both be calculated to obtain the overall stress
generated by a system.

##  I Introduction

Active matter systems constitute an intriguing class of materials whose
constituents have the ability to self-propel, generate internal stress, and
drive the system far from equilibrium. Because classical concepts of
thermodynamics do not apply to nonequilibrium active matter, recent work has
focused on invoking the mechanical pressure (or stress) as a framework to
understand the complex dynamic behaviors of active systems (Takatori14, ;
Takatori16a, ; Ginot15, ; Yang14, ; Mallory14, ; Solon15, ). In particular,
active swimmers exert a unique “swim pressure” as a result of their self
propulsion (Takatori14, ; Yang14, ). A physical interpretation of the swim
pressure is the pressure (or stress) exerted by the swimmers as they interact
with the surrounding boundaries that confine them, similar to molecular or
colloidal solutes that collide into the container walls to exert an osmotic
pressure. The swim pressure is a purely entropic, destabilizing quantity
(Takatori15a, ) that can explain the self-assembly and phase separation of a
suspension of self-propelled colloids into dilute and dense phases that
resemble an equilibrium gas-liquid coexistence (Theurkauff12, ; Palacci13, ;
Redner13, ; Bialke13, ; Stenhammar13, ; Cates10, ).

Existing studies of the pressure of active systems (Takatori14, ; Takatori14b,
; Takatori15a, ; Yang14, ; Mallory14, ; Ginot15, ; Solon15, ) have focused on
overdamped systems where swimmer inertia is neglected (i.e, the particle
Reynolds number is small). However, the swim pressure has no explicit
dependence on the body size and may exist at all scales including larger
swimmers (e.g., fish, birds) where particle inertia is not negligible
(Takatori14, ). The importance of particle inertia is characterized by the
Stokes number, $St_{R}\equiv(M/\zeta)/\tau_{R}$, where $M$ is the particle
mass, $\zeta$ is the hydrodynamic drag factor, and $\tau_{R}$ is the
reorientation time of the active swimmer. Here we analyze the role of nonzero
$St_{R}$ on the mechanical stress exerted by a system of self-propelled bodies
and provide a natural extension of existing pressure concepts to swimmers with
finite inertia. We maintain negligible fluid inertia so that the fluid motion
satisfies the steady Stokes equation.

We consider a suspension of self-propelled spheres of radii $a$ that translate
with an intrinsic swim velocity $\boldsymbol{U}_{0}$ and tumble with a time
$\tau_{R}$ in a continuous Newtonian fluid. The random tumbling results in a
diffusive process for time $t\gg\tau_{R}$ with $D^{swim}=U_{0}^{2}\tau_{R}/6$
in 3D. An isolated active particle generates a swim pressure
$\Pi^{swim}=n\zeta D^{swim}=n\zeta U_{0}^{2}\tau_{R}/6$, where $n$ is the
number density of active particles. We do not include the effects of
hydrodynamic interactions, and there is no macroscopic polar order of the
swimmers or any large-scale collective motion (e.g., bioconvection).

From previous work (Takatori16b, ; Takatori16a, ; Yan15b, ; Ezhilan15, ;
Tailleur09, ), we know that geometric confinement of overdamped active
particles plays a significant role in their dynamics and behavior. Confinement
from potential traps, physical boundaries, and collective clustering can
reduce the average run length and swim pressure of the particles. We have
shown experimentally (Takatori16b, ) that active Brownian particles trapped
inside a harmonic well modifies the swim pressure to $\Pi^{swim}=(n\zeta
U_{0}^{2}\tau_{R}/2)\left(1+\alpha\right)^{-1}$ in 2D, where $\alpha\equiv
U_{0}\tau_{R}/R_{c}$ is a ratio of the run length, $U_{0}\tau_{R}$, to the
size of the trap, $R_{c}$. For weak confinement, $\alpha\ll 1$, we obtain the
‘ideal-gas’ swim pressure of an isolated swimmer. For strong confinement,
$\alpha\gg 1$, the swim pressure decreases as $\Pi^{swim}/(n\zeta
U_{0}^{2}\tau_{R}/2)\sim 1/\alpha$. Confinement reduces the average distance
the swimmers travel between reorientation events, which results in a decreased
swim pressure.

In this work, we find that particle inertia plays a similar role as
confinement by reducing the correlation between the position and self-
propulsive swim force of the swimmers. In addition to the swim pressure,
active swimmers exert the usual kinetic or Reynolds pressure contribution
associated with their average translational kinetic energy. For swimmers with
finite particle inertia we find that the sum of the swim and Reynolds
pressures is the relevant quantity measured by confinement experiments and
computer simulations. We also study systems at finite swimmer concentrations
and extend our existing mechanical pressure theory to active matter of any
size or mass.

An important implication of this work pertains to the computation of
mechanical stresses of colloidal suspensions at the appropriate level of
analysis. Consider the Brownian osmotic pressure of molecular fluids and
Brownian colloidal systems, $\Pi^{B}=nk_{B}T$, where $n$ is the number density
of particles and $k_{B}T$ is the thermal energy. At the Langevin-level
analysis where mass (or inertia) is explicitly included, the Reynolds stress
is the source of the Brownian osmotic pressure,
$-\langle\rho\boldsymbol{U}^{\prime}\boldsymbol{U}^{\prime}\rangle$, where
$\rho$ is the density and $\boldsymbol{U}^{\prime}$ is the velocity
fluctuation. The virial stress or the moment of the Brownian force,
$-n\langle\boldsymbol{x}\boldsymbol{F}^{B}\rangle$, is identically equal to
zero—the position and Brownian force are uncorrelated at the Langevin level.
However, at the overdamped Fokker-Planck or Smoluchowski level where inertia
is not explicitly included, the position and Brownian force are correlated,
and the virial stress $-n\langle\boldsymbol{x}\boldsymbol{F}^{B}\rangle$ is
the source of the Brownian osmotic pressureBrady93 , whereas the Reynolds
stress is zero. The important point is that the sum of the Reynolds and virial
stresses gives the correct Brownian osmotic pressure at both levels of
analysis, as it must be since the osmotic pressure of colloidal Brownian
suspensions is $\Pi^{B}=nk_{B}T$ whether or not inertia is explicitly included
in the analysis.

We report in this work that an identical concept applies for active swimmers.
The virial stress arising from the correlation between the particle position
and its ‘internal’ force, $-n\langle\boldsymbol{x}\boldsymbol{F}\rangle$, is a
term that is separate and in addition to the Reynolds stress associated with
their average translational kinetic energy. Interestingly, the mechanical
stress generated by active swimmers has a nonzero contribution from both the
Reynolds and virial stresses because the internal force associated with self-
propulsion has an autocorrelation that is not instantaneous in time and
instead decays over a finite timescale modulated by the reorientation time of
the active swimmer, $\tau_{R}$:
$\langle\boldsymbol{F}^{swim}(t)\boldsymbol{F}^{swim}(0)\rangle\sim
e^{-2t/\tau_{R}}$. A distinguishing feature of active swimmers compared with
passive Brownian particles is that their direction of self-propulsion can
relax over large timescales, and that the ‘internal’ swim force
autocorrelation cannot in general be described by a delta-function in time.

##  II Swim stress

All self-propelled bodies exert a swim pressure, a unique pressure associated
with the confinement of the active body inside a bounded domain. The swim
pressure is the trace of the swim stress, which is defined as the symmetric
first moment of the self-propulsive force,
$\boldsymbol{\sigma}^{swim}=-n\langle\boldsymbol{x}\boldsymbol{F}^{swim}\rangle^{sym}$,
where $n$ is the number density of swimmers, $\boldsymbol{x}$ is the position,
$\boldsymbol{F}^{swim}$ is the swimmer’s self-propulsive swim force, and
$\langle\cdot\rangle^{sym}$ denotes the symmetric part of the tensor. As we
noted previouslyTakatori17a , the swim stress is properly defined as the
symmetric force moment since the force arises from the fluid which can only
generate symmetric stresses. For the active Brownian particle model, the swim
force can be written as $\boldsymbol{F}^{swim}=\zeta U_{0}\boldsymbol{q}$
where $\boldsymbol{q}$ is a unit vector specifying the swimmer’s direction of
self-propulsion. For a dilute suspension of active particles with negligible
particle inertia the “ideal-gas” swim stress is given by
$\boldsymbol{\sigma}^{swim}=-n\zeta
U_{0}^{2}\tau_{R}\boldsymbol{I}/6=-nk_{s}T_{s}\boldsymbol{I}$, where we define
$k_{s}T_{s}\equiv\zeta U_{0}^{2}\tau_{R}/6$ as the swimmer’s “energy scale”
(force ($\zeta U_{0}$) $\times$ distance ($U_{0}\tau_{R}$)). The swim pressure
(or stress) is entropic in origin and is the principle destabilizing term that
facilitates a phase transition in active systems (Takatori15a, ).

In the absence of any external forces, the motion of an active Brownian
particle is governed by the Langevin equations:

$\displaystyle M\frac{\partial\boldsymbol{U}}{\partial t}$
$\displaystyle=-\zeta(\boldsymbol{U}-\boldsymbol{U}_{0})+\sqrt{2\zeta^{2}D_{0}}\boldsymbol{\Lambda}_{T},$
(1) $\displaystyle I\frac{\partial\boldsymbol{\Omega}}{\partial t}$
$\displaystyle=-\zeta_{R}\boldsymbol{\Omega}+\sqrt{\frac{2\zeta_{R}^{2}}{\tau_{R}}}\boldsymbol{\Lambda}_{R},$
(2)

where $M$ and $I$ are the particle mass and moment-of-inertia,
$\boldsymbol{U}$ and $\boldsymbol{\Omega}$ are the translational and angular
velocities, $\zeta_{R}$ is the hydrodynamic drag factor coupling angular
velocity to torque, $\sqrt{2\zeta^{2}D_{0}}\boldsymbol{\Lambda}_{T}$ and
$\sqrt{2\zeta_{R}^{2}/\tau_{R}}\boldsymbol{\Lambda}_{R}$ are the Brownian
translational force and rotational torque, respectively,
$\boldsymbol{\Lambda}_{T}$ and $\boldsymbol{\Lambda}_{R}$ are unit random
normal deviates, $\tau_{R}\sim 1/D_{R}$ is the reorientation timescale set by
rotational Brownian motion, and $D_{0}$ is the Stokes-Einstein-Sutherland
translational diffusivity. The translational diffusivity and the reorientation
dynamics are modeled with the usual white noise statistics,
$\langle\Lambda_{i}(t)\rangle=0$ and
$\langle\Lambda_{i}(t)\Lambda_{j}(0)\rangle=\delta(t)\delta_{ij}$. The swimmer
orientation $\boldsymbol{q}(t)$ is related to the angular velocity by the
kinematic relation
$\boldsymbol{\Omega}\times\boldsymbol{q}=d\boldsymbol{q}/dt$. The
translational and angular velocities may be combined into a single vector,
$\mathcal{U}=(\boldsymbol{U},\boldsymbol{\Omega})^{T}$, and similarly for the
force and torque, $\mathcal{F}=(\boldsymbol{F},\boldsymbol{L})^{T}$, to obtain
a general solution to the system of ordinary differential equations (Hinch75,
). Although a general solution is available for any particle mass and moment-
of-inertia, inclusion of nonzero moment-of-inertia leads to calculations that
are analytically involved. For convenience and to make analytical progress,
here we summarize the case of zero moment-of-inertia ($I=0$) and focus on
finite mass ($St_{R}\equiv(M/\zeta)/\tau_{R}\neq 0$) to elucidate the effects
of inertia on the dynamics of active matter.

We can solve Eqs 1 and 2 for the swimmer configuration
($\boldsymbol{x}(t),\boldsymbol{q}(t)$), and calculate the swim stress. As
shown in the Appendix, the swim stress for arbitrary particle inertia is

$\boldsymbol{\sigma}^{swim}=-nk_{s}T_{s}\left(\frac{1}{1+2St_{R}}\right)\boldsymbol{I},$
(3)

where we have taken times $t>\tau_{M}$ and $t>\tau_{R}$, $\tau_{M}\equiv
M/\zeta$ is the swimmer momentum relation time, energy scale
$k_{s}T_{s}\equiv\zeta U_{0}^{2}\tau_{R}/6$, and
$St_{R}\equiv\tau_{M}/\tau_{R}=(M/\zeta)/\tau_{R}$ is the Stokes number. For
$St_{R}=0$ we recover the “ideal-gas” swim pressure for an overdamped system:
$\Pi^{swim}=-\text{tr}\boldsymbol{\sigma}^{swim}/3=n\zeta
U_{0}^{2}\tau_{R}/6=nk_{s}T_{s}$ (Takatori14, ). This is precisely the
mechanical force per unit area that a dilute system of confined active micro-
swimmers exert on its surrounding container (Takatori14, ; Yang14, ;
Mallory14, ).

Notice in the other limit as $St_{R}\rightarrow\infty$,
$\boldsymbol{\sigma}^{swim}$ vanishes. Physically, the magnitude of the swim
stress decreases because inertia may translate the swimmer in a trajectory
that is different from the direction of the swim force, reducing the
correlation $\langle\boldsymbol{x}\boldsymbol{F}^{swim}\rangle$ between the
moment arm $\boldsymbol{x}$ and the orientation-dependent swim force
$\boldsymbol{F}^{swim}=\zeta U_{0}\boldsymbol{q}$. Our earlier work
(Takatori16b, ) showed that active particles confined by an acoustic trap
exert a swim pressure that is reduced by a factor of $(1+\alpha)^{-1}$, where
$\alpha\equiv U_{0}\tau_{R}/R_{c}$ is the degree of confinement of the swimmer
run length relative to the size of the trap $R_{c}$. Equation 3 has a
reduction in the swim pressure by a similar factor, $(1+2St_{R})^{-1}$,
suggesting that particle inertia may play a similar role as confinement by
reducing the correlation between the position and self-propulsive direction of
the swimmers.

Particle inertia may be interpreted as imposing a confinement effect on the
swimmers because their effective run length between reorientation events
decreases. The average run length of the swimmers with inertia reduces to
$\sim(U_{0}\tau_{R}-\Delta x^{\prime})$, where $\Delta x^{\prime}=U_{0}\Delta
t$ is the distance over which inertia translates the swimmer along a
trajectory that is independent of the direction of its swim force, and the
time over which this occurs scales with the inertial relaxation time, $\Delta
t\sim M/\zeta$. Substituting these terms into the virial expression for the
swim pressure, we obtain $\Pi^{swim}\sim
n\langle\boldsymbol{x}\cdot\boldsymbol{F}^{swim}\rangle\sim n\zeta
U_{0}^{2}\tau_{R}(1-(M/\zeta)/\tau_{R})$. Using the definition of the Stokes
number $St_{R}\equiv(M/\zeta)/\tau_{R}$ and for small $St_{R}$, we can rewrite
the swim pressure as $\Pi^{swim}\sim n\zeta U_{0}^{2}\tau_{R}(1+St_{R})^{-1}$.
Aside from the factor of $2$ in the denominator (which arises from spatial
dimensionality), this scaling argument agrees with Eq 3 and shows that
particle inertia plays a confining role in the swim pressure, analogous to the
physical confinement of swimmers in a potential well(Takatori16b, ).

Figure 1: Swim and Reynolds pressures of a dilute system of swimmers with
finite inertia, where $\Pi=-\text{tr}\boldsymbol{\sigma}/3$. The blue
($\Pi^{swim}$) and red ($\Pi^{Rey}$) curves and symbols are the analytical
theory of Eqs 3 and 4 and simulation data, respectively. The solid black line
is the sum of the swim and Reynolds stresses. The Brownian osmotic pressure
$\Pi^{B}=nk_{B}T$ has been subtracted from $\Pi^{Rey}$.

##  III Reynolds stress

For systems with finite particle inertia, an additional stress contribution
arises owing to particle acceleration: the Reynolds stress. This term is seen
in Bernoulli’s equation and is associated with the average translational
kinetic energy of a particle,
$\boldsymbol{\sigma}^{Rey}=-nM\langle\boldsymbol{U}^{\prime}\boldsymbol{U}^{\prime}\rangle$.
In atomic or molecular systems this is often referred to as the kinetic
stress. This contribution was not included in previous studies since
overdamped active systems have no particle mass, $M=0$ (i.e.,
$St_{R}\equiv(M/\zeta)/\tau_{R}=0$). As shown in the Appendix, we can use the
solution to Eqs 1 and 2 to obtain the Reynolds stress for arbitrary $St_{R}$,
given by

$\boldsymbol{\sigma}^{Rey}=-nk_{B}T\boldsymbol{I}-nk_{s}T_{s}\left(\frac{1}{1+1/(2St_{R})}\right)\boldsymbol{I},$
(4)

which is a sum of the Brownian osmotic stress
$\boldsymbol{\sigma}^{B}=-nk_{B}T\boldsymbol{I}$ and a self-propulsive
contribution that depends on $St_{R}$. For an overdamped system where
$St_{R}=0$, the self-propulsive contribution to the Reynolds stress vanishes,
justifying the neglect of this term in previous studies of overdamped systems.
Mallory et al (Mallory14, ) analytically calculated the expression of the
Reynolds stress, Eq 4, but not the swim stress, Eq 3; the mechanical pressure
that is measured from the walls of an enclosing container is the sum of the
Reynolds and swim pressures, however.

Notice that the Brownian osmotic stress,
$\boldsymbol{\sigma}^{B}=-nk_{B}T\boldsymbol{I}$, arises solely from the
Reynolds stress and not from taking the virial moment of the Brownian force,
$\langle\boldsymbol{x}\boldsymbol{F}^{B}\rangle$. As stated earlier, this is
precisely because of the delta-function statistics imposed for the Brownian
force in the Langevin-level analysis where mass is explicitly included:
$\langle\boldsymbol{F}^{B}(t)\boldsymbol{F}^{B}(0)\rangle=2k_{B}T\zeta\delta(t)\boldsymbol{I}$.
In contrast, the active stress has a nonzero contribution from both the
Reynolds and swim stresses because the swim force autocorrelation is a
decaying exponential modulated by the reorientation timescale $\tau_{R}$:
$\langle\boldsymbol{F}^{swim}(t)\boldsymbol{F}^{swim}(0)\rangle\sim
e^{-2t/\tau_{R}}$. If one were to model the Brownian force autocorrelation as
one that relaxes over a finite solvent relaxation timescale,
$\langle\boldsymbol{F}^{B}(t)\boldsymbol{F}^{B}(0)\rangle\sim
e^{-2t/\tau_{S}}$, then we would obtain nonzero contributions from both the
Reynolds and virial stresses, with their sum equal to
$\boldsymbol{\sigma}^{B}=-nk_{B}T\boldsymbol{I}$ for all $\tau_{S}$.

For a dilute suspension, the stress exerted by an active swimmer is the sum of
the swim and Reynolds stresses. Adding Eqs 3 and 4, we find

$\boldsymbol{\sigma}^{swim}+\boldsymbol{\sigma}^{Rey}=-n(k_{B}T+k_{s}T_{s})\boldsymbol{I}.$
(5)

Remarkably, the Stokes number $St_{R}$ disappears. The magnitude of the swim
pressure that decreases with increasing $St_{R}$ cancels exactly the increase
in the magnitude of the Reynolds stress. This verifies that swimmers of all
scales exert the pressure, $nk_{s}T_{s}$, regardless of their mass and
inertia. We conducted simulations where the dynamics of active Brownian
particles were evolved following Eqs 1 and 2 using the velocity verlet
algorithm (Allen89, ), and the results are shown in Fig 1. Results from the
simulations agree with our theoretical predictions in Eqs 3 \- 5. In an
experiment or simulation, the average mechanical pressure exerted on a
confining boundary gives the sum of the swim and Reynolds pressures, and not
their separate values.

We have shown previously that random motion gives rise to a micromechanical
stress via the relationship $\boldsymbol{\sigma}=-n\zeta\boldsymbol{D}$, where
$\boldsymbol{D}$ is the effective translational diffusivityTakatori14 . For
overdamped active systems, the self-propulsive contribution to the diffusivity
is described solely through the swim stress,
$\boldsymbol{\sigma}^{swim}=-n\zeta\boldsymbol{D}^{swim}=-n\zeta
U_{0}^{2}\tau_{R}\boldsymbol{I}/6$. With finite particle inertia, the particle
diffusivity is unaltered and the stress-diffusivity relationship still applies
but the stress is now a sum of two independent contributions, as shown in Eq
5.

In the presence of a nonzero moment-of-inertia, there is another dimensionless
Stokes number, $St_{I}\equiv(I/\zeta_{R})/\tau_{R}$, which is a ratio of the
inertial reorientation timescale, $\tau_{I}=I/\zeta_{R}$, and the swimmers’
intrinsic reorientation timescale, $\tau_{R}$. Similar to the translational
Stokes number $St_{R}$ that does not appear in Eq 5, the moment-of-inertia is
not expected to appear explicitly in the total stress generated by an active
swimmer; at long times the random walk diffusive motion is unaffected by the
particle inertia, whether that be translational or rotational.

Figure 2: (A) Swim pressure, $\Pi^{swim}$, and (B) Reynolds pressure,
$\Pi^{Rey}$, as a function of volume fraction of particles, $\phi$, for
different values of $St_{R}\equiv(M/\zeta)/\tau_{R}$ and a fixed reorientation
Péclet number, $Pe_{R}\equiv a/(U_{0}\tau_{R})=0.01$. The symbols and solid
curves are the simulation data and analytical theory, respectively. The
Brownian osmotic pressure $\Pi^{B}=nk_{B}T$ has been subtracted from the
Reynolds pressure.

##  IV Finite concentrations

The results thus far are for a dilute suspension of active swimmers. At finite
concentrations, experiments and computer simulations have observed unique
phase behavior and self-assembly in active matter (Theurkauff12, ; Palacci13,
; Redner13, ; Bialke13, ; Stenhammar13, ; Cates10, ). Recently a new
mechanical pressure theory was developed to provide a phase diagram and a
natural extension of the chemical potential and other thermodynamic quantities
to nonequilibrium active matter (Takatori15a, ).

At finite concentrations of swimmers, dimensional analysis shows that the
nondimensional swim and Reynolds stresses depend in general on
$(St_{R},\phi,Pe_{R},k_{s}T_{s}/(k_{B}T))$, where $\phi=4\pi a^{3}n/3$ is the
volume fraction of active swimmers and $Pe_{R}\equiv a/(U_{0}\tau_{R})$ is the
reorientation Péclet number—the ratio of the swimmer size $a$ to its run
length $U_{0}\tau_{R}$. The ratio $k_{s}T_{s}/(k_{B}T)$ quantifies the
magnitude of the swimmers’ activity ($k_{s}T_{s}\sim\zeta U_{0}^{2}\tau_{R}$)
relative to the thermal energy $k_{B}T$; this ratio can be a large quantity
for typical micro-swimmers.

From previous work on overdamped active systems with negligible particle
inertia (Takatori14, ), we know that $Pe_{R}$ is a key parameter controlling
the phase behavior of active systems. For large $Pe_{R}$ the swimmers reorient
rapidly and take small swim steps, behaving as Brownian walkers; the swimmers
thus do not clump together to form clusters and the system remains
homogeneous. For small $Pe_{R}$ the swimmers obstruct each others’ paths when
they collide for a time $\tau_{R}$ until they reorient. This decreases the run
length of the swimmers between reorientation events and causes the system to
self-assemble into dense and dilute phases resembling an equilibrium liquid-
gas coexistence.

As reported previously for $St_{R}=0$ (Takatori14, ; Takatori15a, ), for small
$Pe_{R}$ the swim pressure decreases with increasing swimmer concentration. To
verify how finite particle inertia affects the swim pressure at larger
concentrations, we conducted simulations by evolving the motion of active
particles following Eqs 1 and 2, with an additional hard-sphere interparticle
force $\boldsymbol{F}^{P}$ that prevents particle overlaps using a potential-
free algorithm (Heyes93, ). Care was taken to ensure that the simulation time
step was small enough to preclude unwanted numerical errors associated with
resolution of particle collisions. We varied the simulation time step from
$dt/\tau_{R}=10^{-5}-10^{-3}$ and found a negligible difference in our
results. As shown in Fig 2A, for finite $St_{R}$ the data from our simulations
are well described by the expression
$\Pi^{swim}=nk_{s}T_{s}(1-\phi-\phi^{2})/(1+2St_{R})$, which is simply a
product of a volume fraction dependence and a Stokes number dependence of Eq
3. The volume fraction dependence $(1-\phi-\phi^{2})$ was used previously to
model the phase behavior of active matter (Takatori15a, ).

The clustering of swimmers reduces their translational velocity
autocorrelation,
$\langle\boldsymbol{U}^{\prime}\boldsymbol{U}^{\prime}\rangle$, and hence
decreases the Reynolds pressure. As shown in Fig 2B, our simulations show that
the Reynolds pressure decreases with concentration, increases with $St_{R}$,
and is well described by the expression
$\Pi^{Rey}=nk_{s}T_{s}(1-\phi-\phi^{2})/(1+1/(2St_{R}))+nk_{B}T$.

The sum of the swim and Reynolds stresses is given by

$\Pi^{swim}+\Pi^{Rey}=nk_{B}T+nk_{s}T_{s}(1-\phi-\phi^{2}),$ (6)

which again has no dependence on $St_{R}$ (nor $Pe_{R}$ for $Pe_{R}<1$
considered here), even at finite $\phi$. Equation 6 is corroborated by our
simulations as shown in Fig 3. This result implies that the existing
mechanical pressure theory (Takatori15a, ) developed for overdamped systems
can be used directly for swimmers with finite Stokes numbers, as long as we
include the Reynolds stress contribution into the active pressure. Inclusion
of the Reynolds stress is critical, as confinement experiments and computer
simulations measure the total active pressure, including both the swim and
Reynolds contributions.

Figure 3: Sum of swim and Reynolds pressures,
$\Pi^{swim}+\Pi^{Rey}=-\text{tr}(\boldsymbol{\sigma}^{swim}+\boldsymbol{\sigma}^{Rey})/3$,
as a function of volume fraction of particles $\phi$ for different values of
$St_{R}\equiv(M/\zeta)/\tau_{R}$ and a fixed reorientation Péclet number,
$Pe_{R}\equiv a/(U_{0}\tau_{R})=0.01$. The symbols and solid curve are the
simulation data and analytical theory of Eq 6, respectively. The Brownian
osmotic pressure $\Pi^{B}=nk_{B}T$ has been subtracted from the total
pressure.

In addition to the swim and Reynolds stresses, interparticle interactions
between the swimmers at finite concentrations give rise to an interparticle
stress, $\boldsymbol{\sigma}^{P}(St_{R},\phi,Pe_{R},k_{s}T_{s}/(k_{B}T))$. For
repulsive interactions, the interparticle (or collisional) pressure,
$\Pi^{P}=-\text{tr}\boldsymbol{\sigma}^{P}/3$, increases monotonically with
concentration and helps to stabilize the system. As shown in Fig 4, we find
that the expression $\Pi^{P}/(nk_{s}T_{s})=3Pe_{R}\phi g(\phi)/(1+0.5St_{R})$
agrees with the simulation data for a fixed value of $Pe_{R}=0.01$, where
$g(\phi)=\left(1-\phi/\phi_{0}\right)^{-1}$ is the pair distribution function
at particle contact, and $\phi_{0}=0.65$ is a parameter obtained from the
interparticle pressure of hard-sphere molecular fluids (Takatori15a, ). We can
add $\Pi^{P}$ to Eq 6 to construct phase diagrams of a system of inertial
swimmers, which are qualitatively similar to those presented in (Takatori15a,
) for small values of $St_{R}$. Adding particle inertia shifts the
stablization to larger $\phi$.

Figure 4: Interparticle collisional pressure, $\Pi^{P}$, as a function of
volume fraction of particles, $\phi$, for different values of
$St_{R}\equiv(M/\zeta)/\tau_{R}$ and a fixed reorientation Péclet number,
$Pe_{R}\equiv a/(U_{0}\tau_{R})=0.01$. The symbols and solid curves are the
simulation data and analytical expression, respectively.

As shown by Batchelor (Batchelor70, ), there may be an additional contribution
to the particle stress arising from local fluctuations in acceleration,
$\boldsymbol{f}^{\prime}$, and is given by
$-(1/V)\sum\int_{V_{p}}\rho\boldsymbol{f}^{\prime}\boldsymbol{r}dV$, where $V$
is the volume of the suspension (fluid plus particles), $V_{p}$ is the volume
of an individual particle, $\rho$ is the uniform density of the particle,
$\boldsymbol{r}$ is a position vector (or the moment arm) from the particle
center, and the summation is over the number of particles in the volume $V$.
For a dilute system of rigid particles, this term arises only from solid body
rotation of the particle and takes the form
$\int_{V_{p}}\boldsymbol{f}^{\prime}\boldsymbol{r}dV=(4\pi
a^{5}/15)(\boldsymbol{\Omega}_{P}\boldsymbol{\Omega}_{P}-\boldsymbol{\Omega}_{P}\cdot\boldsymbol{\Omega}_{P}\boldsymbol{I})$,
where $\boldsymbol{\Omega}_{P}$ is the average angular velocity of the rigid
particle of size $a$. Here the active swimmers have no average angular
velocity, so there is no stress arising from local fluctuations in
acceleration for dilute active systems of rigid particles.

In addition to using a potential-free algorithm to model hard-sphere
particles, we have also tested a short-ranged, repulsive Weeks-Chandler-
Andersen potential with an upper cut-off at particle separation distances of
$r=2^{1/6}(2a)$. Using this softer potential, the swim pressure does not
exhibit a concentration dependence of $(1-\phi-\phi^{2})$ because the
effective radius of the particles decreases as the system becomes denser,
meaning that colliding particles exhibit increasingly large overlaps.
Increasing particle inertia (i.e., larger $St_{R}$) also changed the effective
particle size. As previously stated (Stenhammar14, ), one must use care when
soft potentials are used to model hard-sphere particle collisions because the
effective particle size may depend on system parameters (like $Pe_{R}$ and
$St_{R}$).

##  V Conclusion

Here we presented a mechanical pressure theory for active Brownian particles
with finite inertia. We neglected hydrodynamic interactions between the
swimmers, which may contribute additional terms (like the “hydrodynamic
stresslet” (Saintillan08a, )) to the active pressure. The ratio of the
magnitudes of the hydrodynamic stress to the swim stress is
$\sigma^{H}/\sigma^{swim}\sim(n\zeta U_{0}a)/(n\zeta
U_{0}^{2}\tau_{R})=a/(U_{0}\tau_{R})\equiv Pe_{R}$. The hydrodynamic stress
contribution becomes negligible when phase-separation occurs at low $Pe_{R}$.

We assumed that the surrounding fluid obeys the steady Stokes equations, which
may not be true for larger swimmers that propel themselves using fluid
inertia. However, the concepts of the swim and Reynolds stress apply for
swimmers with a nonlinear hydrodynamic drag factor,
$\zeta_{inertia}(|\boldsymbol{U}|)$, where $|\boldsymbol{U}|$ is the magnitude
of the swimmer velocity. For example, a self-propelled body may experience a
fluid drag that is quadratic in the velocity
$\boldsymbol{F}^{drag}\sim\zeta_{inertia}(|\boldsymbol{U}|)\boldsymbol{U}\sim(\rho_{s}a^{2}|\boldsymbol{U}|)\boldsymbol{U}$,
where $\rho_{s}$ is the fluid density and $a$ is the characteristic size of
the body. The nondimensional Langevin equations would become
$d\boldsymbol{U}/dt=-A|\boldsymbol{U}|(\boldsymbol{U}-\boldsymbol{q})$, where
$\boldsymbol{q}$ is the orientation vector of the swimmer and
$A=\rho_{s}a^{2}U_{0}\tau_{R}/M\sim(1/St_{R})(\zeta_{inertia}/(\eta a))$ is
the relevant quantity that must be varied.

##  VI acknowledgments

The authors thank Luis Nieves-Rosado for his contributions to this work as
part of the Undergraduate Research Program at the California Institute of
Technology. S.C.T. acknowledges support by the Gates Millennium Scholars
fellowship and the National Science Foundation Graduate Research Fellowship
under Grant No. DGE-1144469. This work was also supported by NSF Grant CBET
1437570.

##  VII Appendix

Integrating Eq 1 twice in time, we obtain the position of the swimmer,

$\boldsymbol{x}(t)=\boldsymbol{x}(0)+\boldsymbol{U}(0)\tau_{M}\left(1-e^{-t/\tau_{M}}\right)+\\\
\int_{0}^{t}\left(U_{0}\boldsymbol{q}(t^{\prime})+\sqrt{2D_{0}}\boldsymbol{\Lambda}_{T}(t^{\prime})\right)(1-e^{-(t-t^{\prime})/\tau_{M}})\
dt^{\prime},$ (7)

where $\boldsymbol{x}(0)$ and $\boldsymbol{U}(0)$ are the arbitrary initial
position and velocity, respectively, $\tau_{M}\equiv M/\zeta$ is the momentum
relaxation time, $U_{0}$ is the intrinsic swimmer velocity, $\boldsymbol{q}$
is the unit orientation vector of the swimmer, and $\boldsymbol{\Lambda}_{T}$
is a unit random deviate. Integrating Eq 2 and using the kinematic relation,
$\boldsymbol{\Omega}\times\boldsymbol{q}=d\boldsymbol{q}/dt$, we obtain

$\frac{d\boldsymbol{q}}{dt}=\boldsymbol{\Omega}(0)\times\boldsymbol{q}(t)e^{-t/\tau_{I}}+\frac{1}{\tau_{I}}\sqrt{\frac{2}{\tau_{R}}}\int_{0}^{t}\boldsymbol{\Lambda}_{R}(t^{\prime})\times\boldsymbol{q}(t^{\prime})e^{-(t-t^{\prime})/\tau_{I}}\
dt^{\prime},$ (8)

where $\boldsymbol{\Omega}(0)$ is the initial angular velocity,
$\tau_{I}=I/\zeta_{R}$ is the angular momentum relaxation time, and
$\boldsymbol{\Lambda}_{R}$ is a unit random deviate. Equation 8 is of the form
$dq_{i}/dt=A_{ik}(t)q_{k}$, where $A_{ik}(t)$ is a coefficient matrix. The
general solution of Eq 8 is
$q_{i}(t)=q_{k}(0)e^{\int_{0}^{t}A_{ik}(t^{\prime})\ dt^{\prime}}$, where
$q_{k}(0)$ is an arbitrary initial orientation of the swimmer.

We are interested in the orientation autocorrelation

$\langle q_{i}(t)q_{n}(t^{\prime})\rangle=\frac{1}{3}\delta_{nk}\langle
e^{\int_{t^{\prime}}^{t}A_{ik}(t^{\prime\prime})dt^{\prime\prime}}\rangle,$
(9)

where
$A_{ik}(t)=\epsilon_{ijk}\left(\Omega_{j}(0)e^{-t/\tau_{I}}+\sqrt{2/(\tau_{R}\tau_{I}^{2})}\int_{0}^{t}\Lambda_{j}(t^{\prime})e^{-(t-t^{\prime})/\tau_{I}}dt^{\prime}\right)$
is the coefficient matrix, and $\boldsymbol{\epsilon}$ is the unit alternating
tensor. In the limit of small angular momentum relaxation time,
$\tau_{I}\rightarrow 0$, we obtain
$A_{ik}(t)=\sqrt{2/\tau_{R}}\epsilon_{ijk}\Lambda_{j}(t)$, and

$\langle
q_{i}(t)q_{n}(t^{\prime})\rangle=\frac{1}{3}\delta_{in}e^{-2(t-t^{\prime})/\tau_{R}}.$
(10)

Notice that as $\tau_{R}\rightarrow 0$ the autocorrelation becomes a delta
function and the swimmer reorients rapidly and behaves as a Brownian walker.

Using Eqs 7, 10, and the swim force $\boldsymbol{F}^{swim}\equiv\zeta
U_{0}\boldsymbol{q}$, the swim stress is given by

$\boldsymbol{\sigma}^{swim}=-n\langle\boldsymbol{x}\boldsymbol{F}^{swim}\rangle^{sym}=\frac{n}{3}\zeta
U_{0}^{2}\boldsymbol{I}\left[\frac{\tau_{R}}{2}\left(1-e^{-2t/\tau_{R}}\right)-\right.\\\
\left.\frac{1}{2/\tau_{R}+1/\tau_{M}}\left(1-e^{-(2/\tau_{R}+1/\tau_{M})t}\right)\right],$
(11)

where we have used that
$\langle\boldsymbol{x}(0)\boldsymbol{F}^{swim}(t)\rangle=\langle\boldsymbol{U}(0)\boldsymbol{F}^{swim}(t)\rangle=\langle\boldsymbol{F}^{swim}(t^{\prime})\boldsymbol{\Lambda}_{T}(t)\rangle=\boldsymbol{0}$.
Taking times $t>\tau_{M}$ and $t>\tau_{R}$, we obtain Eq 3 of the main text.

Following a similar procedure, the Reynolds stress is given by

$\boldsymbol{\sigma}^{Rey}=-nM\langle\boldsymbol{U}^{\prime}\boldsymbol{U}^{\prime}\rangle=-nM\langle\boldsymbol{U}(0)\boldsymbol{U}(0)\rangle
e^{-2t/\tau_{M}}-\\\
\frac{nM}{3}\left(\frac{U_{0}}{\tau_{M}}\right)^{2}\boldsymbol{I}\left\\{\frac{1}{2/\tau_{R}+1/\tau_{M}}\left[\frac{\tau_{M}}{2}\left(1-e^{-2t/\tau_{M}}\right)-\right.\right.\\\
\left.\left.\frac{1}{-2/\tau_{R}+1/\tau_{M}}\left(e^{-(2/\tau_{R}+1/\tau_{M})t}-e^{-2t/\tau_{M}}\right)\right]+\right.\\\
\left.\frac{1}{-2/\tau_{R}+1/\tau_{M}}\left[-\frac{\tau_{M}}{2}\left(1-e^{-2t/\tau_{M}}\right)+\frac{1}{2/\tau_{R}+1/\tau_{M}}\left(1-e^{-(2/\tau_{R}+1/\tau_{M})t}\right)\right]\right\\}-\\\
nk_{B}T\left(1-e^{-2t/\tau_{M}}\right)\boldsymbol{I}.$ (12)

Taking times $t>\tau_{M}$ and $t>\tau_{R}$, we obtain the Reynolds stress as
given in Eq 4 of the main text.

## References

  * [1] S. C. Takatori, W. Yan, and J. F. Brady.  Swim pressure: Stress generation in active matter.  Phys Rev Lett, 113(2):028103, 2014.
  * [2] S. C. Takatori and J. F. Brady.  Forces, stresses and the (thermo?) dynamics of active matter.  Curr Opin Colloid Interface Sci, 21:24–33, 2016.
  * [3] F. Ginot, I. Theurkauff, D. Levis, C. Ybert, L. Bocquet, L. Berthier, and C. Cottin-Bizonne.  Nonequilibrium equation of state in suspensions of active colloids.  Phys Rev X, 5(1):011004, 2015.
  * [4] X. Yang, M. L. Manning, and M. C. Marchetti.  Aggregation and segregation of confined active particles.  Soft Matter, 10(34):6477–6484, 2014.
  * [5] S. A. Mallory, A. Šarić, C. Valeriani, and A. Cacciuto.  Anomalous thermomechanical properties of a self-propelled colloidal fluid.  Phys Rev E, 89(5):052303, 2014.
  * [6] A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates, and J. Tailleur.  Pressure and phase equilibria in interacting active Brownian spheres.  Phys Rev Lett, 114(19):198301, 2015.
  * [7] S. C. Takatori and J. F. Brady.  Towards a thermodynamics of active matter.  Phys Rev E, 91(3):032117, 2015.
  * [8] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet.  Dynamic clustering in active colloidal suspensions with chemical signaling.  Phys Rev Lett, 108(26):268303, 2012.
  * [9] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and P. M. Chaikin.  Living crystals of light-activated colloidal surfers.  Science, 339(6122):936–940, 2013.
  * [10] G. S. Redner, M. F. Hagan, and A. Baskaran.  Structure and dynamics of a phase-separating active colloidal fluid.  Phys Rev Lett, 110(5):055701, 2013.
  * [11] J. Bialké, H. Löwen, and T. Speck.  Microscopic theory for the phase separation of self-propelled repulsive disks.  Europhys Lett, 103(3):30008, 2013.
  * [12] J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo, and M. E. Cates.  Continuum theory of phase separation kinetics for active Brownian particles.  Phys Rev Lett, 111(14):145702, 2013.
  * [13] M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and J. Tailleur.  Arrested phase separation in reproducing bacteria creates a generic route to pattern formation.  Proc Natl Acad Sci U.S.A, 107(26):11715–11720, 2010.
  * [14] S. C. Takatori and J. F. Brady.  Swim stress, motion, and deformation of active matter: Effect of an external field.  Soft Matter, 10(47):9433–9445, 2014.
  * [15] S. C. Takatori, R. De Dier, J. Vermant, and J. F. Brady.  Acoustic trapping of active matter.  Nat Commun, 7, 2016.
  * [16] W. Yan and J. F. Brady.  The force on a boundary in active matter.  J Fluid Mech, 785:R1, 2015.
  * [17] B. Ezhilan, R. Alonso-Matilla, and D. Saintillan.  On the distribution and swim pressure of run-and-tumble particles in confinement.  J Fluid Mech, 781:R4, 2015.
  * [18] J. Tailleur and M. E. Cates.  Sedimentation, trapping, and rectification of dilute bacteria.  Europhys Lett, 86(6):60002, 2009.
  * [19] J. F. Brady.  Brownian motion, hydrodynamics, and the osmotic pressure.  J Chem Phys, 98(4):3335–3341, 1993.
  * [20] S. C. Takatori and J. F. Brady.  Superfluid behavior of active suspensions from diffusive stretching.  Phys Rev Lett, 118(1):018003, 2017.
  * [21] E. J. Hinch.  Application of the langevin equation to fluid suspensions.  J Fluid Mech, 72(3):499–511, 1975.
  * [22] M. P. Allen and D. J. Tildesley.  Computer simulation of liquids.  Clarendon Press, Oxford, UK, 1989.
  * [23] D. M. Heyes and J. R. Melrose.  Brownian dynamics simulations of model hard-sphere suspensions.  J Non-Newtonian Fluid Mech, 46(1):1–28, 1993.
  * [24] G. K. Batchelor.  The stress system in a suspension of force-free particles.  J Fluid Mech, 41(03):545–570, 1970.
  * [25] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E. Cates.  Phase behaviour of active Brownian particles: The role of dimensionality.  Soft Matter, 10(10):1489–1499, 2014.
  * [26] D. Saintillan and M. J. Shelley.  Instabilities, pattern formation, and mixing in active suspensions.  Phys Fluids, 20(12):123304–15, 2008.


# HOW A PROJECTIVELY FLAT GEOMETRY REGULATES $F(R)$-GRAVITY THEORY?

Tee-How Loo 0000-0003-4099-9843<EMAIL_ADDRESS>Institute of Mathematical
Sciences, University of Malaya, 50603 Kuala Lumpur, Malaysia Avik De
0000-0001-6475-3085<EMAIL_ADDRESS>Department of Mathematical and Actuarial
Sciences, Universiti Tunku Abdul Rahman,
Jalan Sungai Long, 43000 Cheras, Malaysia Sanjay Mandal 0000-0003-2570-2335
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute of
Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. P.K. Sahoo 0000-0003-2130-8832
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute
of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India.

###### Abstract

In the present paper we examine a projectively flat spacetime solution of
$F(R)$-gravity theory. It is seen that once we deploy projective flatness in
the geometry of the spacetime, the matter field has constant energy density
and isotropic pressure. We then make the condition weaker and discuss the
effects of projectively harmonic spacetime geometry in $F(R)$-gravity theory
and show that the spacetime in this case reduces to a generalised Robertson-
Walker spacetime with a shear, vorticity, acceleration free perfect fluid with
a specific form of expansion scalar presented in terms of the scale factor.
Role of conharmonic curvature tensor in the spacetime geometry is also briefly
discussed. Some analysis of the obtained results are conducted in terms of
couple of $F(R)$-gravity models.

##  1 Introduction

In 1918, H. Weyl constructed, what today is known as the Weyl geometry [1] and
in 1921, in a subsequent paper [2], he defined a projective curvature tensor
$P^{l}{}_{ijk}$ as a projectively invariant component of the Riemann curvature
tensor. Suppose in an $n$-dimensional (pseudo-)Riemannian manifold $M^{n}$, we
consider a torsion-less affine connection $\Gamma^{l}_{\>\>jk}$, with the
corresponding curvature tensor $R^{l}{}_{ijk}$ and Ricci tensor $R_{ik}$
respectively. The projective curvature tensor is given by

$\displaystyle
P^{l}{}_{ijk}=R^{l}{}_{ijk}+L_{ji}\delta^{l}_{k}-L_{ki}\delta^{l}_{j}+(L_{jk}+L_{kj})\delta^{l}_{j},$
(1)

where

$L_{ik}=\frac{nR_{ki}+R_{ik}}{n^{2}-1}.$

A curve $x^{l}(t)$ in $M$ is said to be path induced by the affine connection
$\Gamma^{l}{}_{jk}$ if it satisfies the following system of differential
equations

$\ddot{x}^{l}+\dot{x}^{j}\dot{x}^{k}\Gamma^{l}{}_{jk}=\alpha\dot{x}^{l},$

where $\alpha(t)$ is a scalar tensor and $\dot{x}^{l}=dx^{l}/dt$. Suppose
there is another torsion-less affine connection $\tilde{\Gamma}^{l}{}_{jk}$ on
$M$ such that both $\Gamma^{l}{}_{jk}$ and $\tilde{\Gamma}^{l}{}_{jk}$ give
the same system of paths. Then we must have the relation

$\tilde{\Gamma}^{l}_{\>\>jk}-\Gamma^{l}_{\>\>jk}=p_{k}\delta^{l}_{j}+p_{j}\delta^{l}_{k},$

for some 1-form $p_{i}$. We say that $\Gamma^{l}{}_{jk}$ and
$\tilde{\Gamma}^{l}{}_{jk}$ are projectively related. Furthermore, if
$\tilde{P}^{l}_{ijk}$ denote the projective curvature tensor corresponding to
$\tilde{\Gamma}^{l}{}_{jk}$, then the following relation hold.

$\tilde{P}^{l}{}_{ijk}=P^{l}{}_{ijk}.$

This means that the project curvature tensor is invariant under a projective
change of $\Gamma^{l}{}_{jk}$. In particular, if $\Gamma^{l}{}_{jk}$ is the
Levi-Civita connection, then (1) becomes

$P_{hijk}=R_{hijk}+\frac{R_{ij}g_{hk}-R_{ik}g_{jh}}{n-1}.$ (2)

It is known that the projective curvature tensor coincides with the Weyl
conformal curvature tensor if and only if the underline metric is Einstein
[3]. Projective curvature symmetry which preserves the Weyl projective
curvature tensor carries significant interest in the standard theory of
gravity [4] governed by Einstein’s field equations (EFE):

$R_{ij}-\frac{R}{2}g_{ij}=\kappa^{2}T_{ij}.$ (3)

With the 1998 discovery of the accelerating expansion of the universe using
type Ia supernovae, and due to the limited knowledge about the character of
dark energy, an alternative theory of gravity was resurfaced by modifying the
Einstein-Hilbert action by replacing the Ricci scalar $R$ with an arbitrary
function $F(R)$. We should mention here that the Weyl geometry mentioned in
the beginning of the present article was primarily proposed as an alternative
of Einstein’s gravity theory and a unified theory of gravitation and
electromagnetism [5]. Some of the initial $F(R)$-models, for example,
$F(R)=R+\alpha R^{2}$, $(\alpha>0)$ proposed by Starobinsky [6] and
$F(R)=R-\mu^{4}/R$, $(\mu>0)$ proposed by Carroll et al. [7] were aimed at
explaining respectively the cosmic inflation without using any scalar field
and late-time acceleration as a pure gravitational effect without using any
dark energy and these models with all their limitations, still managed to
popularise the $F(R)$-models, in general. But it also showed how non-trivial
the process of finding a viable, stable $F(R)$-model is. For a detailed survey
on $F(R)$-gravity theory see [8] and the references therein.

The action term

$S=\frac{1}{2\kappa^{2}}\int F(R)\sqrt{-g}d^{4}x+\int L_{m}\sqrt{-g}d^{4}x,$

gives us the well-known field equations of the $F(R)$-gravity as

$F_{R}(R)R_{ij}-\frac{1}{2}F(R)g_{ij}+(g_{ij}\Box-\nabla_{i}\nabla_{j})F_{R}(R)=\kappa^{2}T_{ij},$
(4)

where $\Box=\nabla^{k}\nabla_{k}$, $L_{m}$ is the matter Lagrangian, and

$T_{ij}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta g^{ij}},$

is the stress-energy tensor. We consider a perfect fluid type stress-energy
tensor given by

$T_{ij}=ph_{ij}+\rho u_{i}u_{j},$ (5)

where $p$ and $\rho$ are the isotropic pressure and energy density,
respectively related by an equation of state (EoS) $\omega$ as $p=\omega\rho$.
$h_{ij}=g_{ij}+u_{i}u_{i}$ denote the orthogonl projector.

The field equations (4) can be rewritten as an effective theory

$R_{ij}-\frac{R}{2}g_{ij}=\frac{\kappa^{2}}{F_{R}(R)}T_{ij}^{\text{eff}}$ (6)

where

$\displaystyle\kappa^{2}T_{ij}^{\text{eff}}=\kappa^{2}T_{ij}+\frac{F(R)-RF_{R}(R)}{2}g_{ij}+\left(\nabla_{i}\nabla_{j}-g_{ij}\Box\right)F_{R}(R).$

We often use the term

$\kappa^{2}T_{ij}^{\text{curv}}=\frac{F(R)-RF_{R}(R)}{2}g_{ij}+\left(\nabla_{i}\nabla_{j}-g_{ij}\Box\right)F_{R}(R).$

Different type of symmetries of the spacetimes, although initiated probably
for plain mathematical purpose, are useful in the gravity theories to find
exact solutions of the field equations. Finding a solution of Einstein’s field
equations or of any modified gravity theories like the present $F(R)$-theory
means finding a metric $g_{ij}$ given in terms of elementary functions (the
geometry) where the physics part, the $T_{ij}$ is guided by physical laws such
as energy conditions etc. This is already very complicated considering the
higher order non-linear differentiation of the metric tensor involved. And
simply a local solution of the set of differential equations is not enough,
either, the global solution or topological properties are to be taken into
account. So traditionally we always impose some kind of symmetry into the
metric, the Riemannian tensor, the Ricci tensor and alike. Then once the
solution is found, we study the physical interpretation. Spherical symmetric,
axisymmetric static, non-static spherically symmetric, plane symmetric etc are
some of the popular outcome of this. In the present paper we impose two
constraints on the spacetime, namely, the projectively flatness and null
divergence of the projective curvature tensor to investigate the role of the
projective curvature tensor in $F(R)$-gravity theories.

After the introduction, we discuss the field equation of $F(R)$-gravity in a
projectively flat spacetime. In the next section we analyse our findings using
two popular models of $F(R)$-gravity. In the next section, we investigate the
effect of divergence free projective curvature tensor followed by a
discussion.

##  2 Projectively flat spacetime

We begin this section with a lemma which is useful in dealing with the
modified $F(R)$-gravity as an effective counterpart of the standard theory of
gravity, in general.

######  Lemma 2.1.

In a perfect fluid spacetime solution of the $F(R)$-gravity theories with
constant Ricci scalar $R$, the effective pressure $p^{\text{eff}}$ and energy
density $\rho^{\text{eff}}$ satisfy the following relation:

$\displaystyle p^{\text{eff}}=p+\frac{F(R)-RF_{R}(R)}{2\kappa^{2}},$
$\displaystyle\rho^{\text{eff}}=\rho-\frac{F(R)-RF_{R}(R)}{2\kappa^{2}}.$

###### Proof.

The field equations (4) of $F(R)$-theories of gravity in case of a constant
$R$ reduces to the form

$R_{ij}-\frac{R}{2}g_{ij}=\frac{\kappa^{2}}{F_{R}(R)}T_{ij}+\frac{F(R)-RF_{R}(R)}{2F_{R}(R)}g_{ij}.$
(7)

Using the perfect fluid form (5), simple calculations produce the result. ∎

######  Theorem 2.1.

In any projectively flat perfect fluid spacetime solution of $F(R)$-gravity
theory, the energy density $\rho$ and isotropic pressure $p$ are constants and
separately satisfy certain relations with the gravity sector as follows:

$p=\frac{RF_{R}(R)}{4\kappa^{2}}-\frac{F(R)}{2\kappa^{2}},\qquad\rho=\frac{F(R)}{2\kappa^{2}}-\frac{RF_{R}(R)}{4\kappa^{2}}.$

###### Proof.

In a 4-dimensioanl projectively flat spacetime, the projective curvature
tensor vanishes, therefore from (2) we can express the Riemannian curvature
tensor as

$3R_{hijk}+R_{ij}g_{hk}-R_{ik}g_{jh}=0.$ (8)

Using the identities of Riemannian curvature tensor, $R_{hijk}+R_{ihjk}=0$,
from (8) we can readily obtain

$\displaystyle R_{ij}g_{hk}-R_{ik}g_{jh}=R_{hk}g_{ij}-R_{hj}g_{ik}$

which on contraction of $h$ and $k$ gives

$R_{ij}=\frac{R}{4}g_{ij}.$ (9)

Again, covariantly differentiating (8) we get

$3\nabla^{h}R_{kijh}=\nabla_{j}R_{ik}-\nabla_{k}R_{ij}$ (10)

Using the contracted Bianchi’s second identity
$\nabla^{m}R_{jklm}=\nabla_{j}R_{kl}-\nabla_{k}R_{jl}$ in (10) we conclude
that the Ricci curvature is of Codazzi type

$\nabla_{i}R_{jk}-\nabla_{j}R_{ik}=0.$ (11)

Contracting $j$ and $k$ in (11) and using the identity
$\nabla^{k}R_{ik}=\frac{1}{2}\nabla_{i}R$ we obtain

$\nabla_{i}R=0.$ (12)

(9) and (12) together imply

$\nabla_{k}R_{ij}=0.$ (13)

Due to (9), for constant $R$ the field equations (7) reduces to

$-\frac{R}{4}g_{ij}=\left[\frac{\kappa^{2}}{F_{R}(R)}p+\frac{F(R)-RF_{R}(R)}{2F_{R}(R)}\right]g_{ij}+\frac{\kappa^{2}}{F_{R}(R)}(p+\rho)u_{i}u_{j}.$
(14)

Contracting $i$ and $j$ in (14) we get

$R=\frac{3\kappa^{2}p}{F_{R}(R)}-\frac{\kappa^{2}\rho}{F_{R}(R)}+\frac{2F(R)}{F_{R}(R)}.$
(15)

Again transvecting (14) by $u^{j}$ and after some calculations we get

$R=\frac{2F(R)}{F_{R}(R)}-\frac{4\kappa^{2}\rho}{F_{R}(R)}.$ (16)

(15) and (16) together imply

$\displaystyle p=\frac{RF_{R}(R)}{4\kappa^{2}}-\frac{F(R)}{2\kappa^{2}},$

and

$\displaystyle\rho=\frac{F(R)}{2\kappa^{2}}-\frac{RF_{R}(R)}{4\kappa^{2}}.$

∎

######  Remark 2.1.

In a projectively flat spacetime, Ricci scalar $R$ is constant.

Using Lemma 2.1 and Theorem 2.1 it is straightforward to obtain the following
result:

######  Corollary 2.1.

In any projectively flat perfect fluid spacetime solution of $F(R)$-gravity
theory, the effective energy density $\rho^{\text{eff}}$ and isotropic
pressure $p^{\text{eff}}$ are constants and given by the expressions:

$\displaystyle
p^{\text{eff}}=-\frac{RF_{R}(R)}{4\kappa^{2}},\qquad\rho^{\text{eff}}=\frac{RF_{R}(R)}{4\kappa^{2}}.$

##  3 Graphical analysis of the scenario

Energy conditions are one of the greatest tools to examine the cosmological
models’ self-stability, which are generally derived from the well-known
Raychaudhuri equation [9, 10, 11]. It also helps us describe the space-time
curve’s geometrical behavior such as spacelike, timelike and lightlike, and
gives new insights into dreadful singularities [12, 13]. The Raychaudhuri
equations can be written if the following forms

$\frac{d\theta}{d\tau}=-\frac{1}{3}\theta^{2}-\sigma_{\mu\nu}\sigma^{\mu\nu}+\omega_{\mu\nu}\omega^{\mu\nu}-R_{\mu\nu}u^{\mu}u^{\nu}\,,$
(17)
$\frac{d\theta}{d\tau}=-\frac{1}{2}\theta^{2}-\sigma_{\mu\nu}\sigma^{\mu\nu}+\omega_{\mu\nu}\omega^{\mu\nu}-R_{\mu\nu}n^{\mu}n^{\nu}\,,$
(18)

where $\theta$ is the expansion factor, $n^{\mu}$ is the null vector, and
$\sigma^{\mu\nu}$ and $\omega_{\mu\nu}$ are, respectively, the shear and the
rotation associated with the vector field $u^{\mu}$. For attractive gravity,
equations (17), and (18) satisfy the following conditions

$\displaystyle R_{\mu\nu}u^{\mu}u^{\nu}\geq 0\,,$ $\displaystyle
R_{\mu\nu}n^{\mu}n^{\nu}\geq 0\,.$

Therefore, if we are working with a perfect fluid matter distribution and
collaborating with the work from Capozziello et al., [14], the energy
conditions recovered from standard GR are

  * •

Strong energy conditions (SEC) if $\rho^{eff}+3p^{eff}\geq 0\,$;

  * •

Weak energy conditions (WEC) if $\rho^{eff}\geq 0,\rho^{eff}+p^{eff}\geq 0\,$;

  * •

Null energy condition (NEC) if $\rho^{eff}+p^{eff}\geq 0\,$;

  * •

Dominant energy conditions (DEC) if $\rho^{eff}\geq 0,|p^{eff}|\leq\rho\,$.

Using Lemma 2.1, we have rewritten the energy conditions as follows

  * •

Strong energy conditions (SEC) if $\rho+3p\geq 0\,$ with $RF_{R}(R)-F(R)\geq
0$;

  * •

Weak energy conditions (WEC) if $\rho\geq 0,\rho+p\geq 0\,$ with
$-RF_{R}(R)+F(R)\geq 0$;

  * •

Null energy condition (NEC) if $\rho+p\geq 0\,$;

  * •

Dominant energy conditions (DEC) if $\rho\geq 0,|p|\leq\rho\,$.

Now, in the following subsections, we discuss the energy conditions for two
specific functions of Ricci scalar $R$ under this formulation. Furthermore, we
would like to note here that under projectively flat spacetime, $R$ is
constant. Therefore, the continuity equation follows directly i.e.,
$\nabla^{i}T_{ij}=0$.

###  3.1 Model-1: $F(R)=R-\alpha\left(1-e^{-\frac{R}{\alpha}}\right)$

Here we presume $F(R)$ as a exponential function with free parameter $\alpha$.
Therefore, the energy density and pressure can be written as

$\rho=\frac{1}{4}e^{-\frac{R}{\alpha}}\left(2\alpha+Re^{R/\alpha}-2\alpha
e^{R/\alpha}+R\right),$ (19)

$p=-\frac{1}{4}e^{-\frac{R}{\alpha}}\left(2\alpha+Re^{R/\alpha}-2\alpha
e^{R/\alpha}+R\right).$ (20)

Now, using Eqn. (19) and (20), one can discuss the energy conditions with this
setup. The profiles of all energy conditions are depicted in Fig 1, 2, 3
except WEC. In this setup, $\rho+p$ reduces to zero. The energy density is
positive throughout the evolution of Ricci Scalar $R>1$ and $\alpha>0$. From
Figure 1, one can observe that the energy density is high for higher values of
$R$, which suggests the energy density of our universe was high in the early
time, and later it decreases as Ricci scalar $R$ decreases with time. Also,
NEC is part of WEC. Therefore, WEC and NEC are satisfied. The profile of DEC
is presented in Fig. 2, and it takes its value in the positive range. SEC is
violated, and this result indicates the late-time acceleration of the universe
[15, 16, 17, 18, 19, 20]. In addition, the equation of state parameter,
$\omega=p/\rho=-1$ for this formulation. Moreover, All of the results are
compatible with the $\Lambda$CDM model [21].

Figure 1: Evolution of energy density with respect to $R$ and $\alpha$ for
model-1. Figure 2: Evolution of DEC with respect to $R$ and $\alpha$ for
model-1. Figure 3: Evolution of SEC with respect to $R$ and $\alpha$ for
model-1.

###  3.2 Model-2: $F(R)=R-\log(\beta R)$

In this subsection, we consider $F(R)$ as a logarithmic function of the Ricci
scalar $R$ with free parameter $\beta>0$. For this $F(R)$, we can write the
energy density and pressure as follows

$\rho=\frac{1}{4}(-2\log(\beta R)+R+1),$ (21)

$p=\frac{1}{4}(2\log(\beta R)-R-1).$ (22)

Using Eqn. (21) and (22), we present the profiles of energy density $\rho$,
DEC, and SEC in Fig. 4, 5, and 6, respectively. We observe that the energy
density and DEC are satisfied from those figures, whereas SEC is violated.
Moreover, WEC and NEC are also satisfied as $\rho+p$ reduces to zero for this
formulation. Nonetheless, the equation of state parameter $\omega$ takes its
value as $-1$. These results indicate the accelerated expansion of the
universe and compatible with the $\Lambda$CDM model [21].

Figure 4: Evolution of energy density with respect to $R$ and $\beta$ for
model-2. Figure 5: Evolution of DEC with respect to $R$ and $\beta$ for
model-2. Figure 6: Evolution of SEC with respect to $R$ and $\beta$ for
model-2.

##  4 Projectively harmonic spacetime

In this section, instead of projective flatness, we impose only vanishing
projective divergence, namely, $\nabla^{h}P_{hijk}=0$, into the spacetime
geometry to see how this weaker restriction affects the theories of gravity.

Putting $n=4$ in (2) and then taking divergence we obtain same equation as
(10) for a divergence free projective curvature. Thereafter, proceeding
similarly we get the results as in equations (11) and (12). Since the Ricci
scalar $R$ is constant, we take covariant derivative of (7) and use (11) to
produce

$\nabla_{i}T_{jk}-\nabla_{j}T_{ik}=0.$ (23)

Hence we get the following result:

######  Theorem 4.1.

In a projectively harmonic spacetime solution of $F(R)$-gravity theories, the
stress-energy tensor is of Codazzi type.

Now, covariant derivative of (5) gives

$\displaystyle\nabla_{k}T_{ij}=$
$\displaystyle\nabla_{k}pg_{ij}+(\nabla_{k}\rho+\nabla_{k}p)u_{i}u_{j}+(\rho+p)\\{\nabla_{k}u_{i}u_{j}+u_{i}\nabla_{k}u_{j}\\}$
$\displaystyle=$ $\displaystyle\nabla_{k}ph_{ij}+\nabla_{k}\rho
u_{i}u_{j}+(\rho+p)\\{\nabla_{k}u_{i}u_{j}+u_{i}\nabla_{k}u_{j}\\},$

which due to (23) produces

$\displaystyle\nabla_{k}ph_{ij}+$ $\displaystyle\nabla_{k}\rho
u_{i}u_{j}+(\rho+p)\\{\nabla_{k}u_{i}u_{j}+u_{i}\nabla_{k}u_{j}\\}$
$\displaystyle=$ $\displaystyle\nabla_{j}ph_{ik}+\nabla_{j}\rho
u_{i}u_{k}+(\rho+p)\\{\nabla_{j}u_{i}u_{k}+u_{i}\nabla_{j}u_{k}\\}.$ (24)

Transvecting (4) with $u^{i}$ gives

$\displaystyle\nabla_{k}\rho u_{j}+(\rho+p)\nabla_{k}u_{j}=$
$\displaystyle\nabla_{j}\rho u_{k}+(\rho+p)\nabla_{j}u_{k}.$ (25)

Substituting this back into (4) gives

$\displaystyle\nabla_{k}ph_{ij}+(\rho+p)\nabla_{k}u_{i}u_{j}=$
$\displaystyle\nabla_{j}ph_{ik}+(\rho+p)\nabla_{j}u_{i}u_{k}.$ (26)

Transvecting (26) with $h^{ij}$ and $h^{kl}$ yields $2h^{lk}\nabla_{k}p=0$,
implying

$\displaystyle\nabla_{k}p=-\dot{p}u_{k};\qquad(\dot{p}=u^{i}\nabla_{i}p).$
(27)

Transvecting (26) with $u^{k}$, with the help of (27), we obtain

$\displaystyle(\rho+p)\nabla_{j}u_{i}=-\dot{p}h_{ji}-(\rho+p)u_{j}\dot{u}_{i};\qquad(\dot{u}_{i}=u^{k}\nabla_{k}u_{i}).$
(28)

By using (25) and (28), we obtain

$\displaystyle\nabla_{k}\rho u_{j}-(\rho+p)u_{k}\dot{u}_{j}=$
$\displaystyle\nabla_{j}\rho u_{k}-(\rho+p)u_{j}\dot{u}_{k}.$

Transvecting this equation with $u^{j}$ and $h^{lk}$ gives

$\displaystyle
h^{lk}\nabla_{k}\rho=-(\rho+p)\dot{u}_{k}h^{lk}=-(\rho+p)\dot{u}^{l}.$ (29)

The trace of (7)

$\displaystyle F_{R}(R)R=2F(R)+\kappa^{2}(3p-\rho),$ (30)

from (27) and (29) gives

$\kappa^{2}(\rho+p)\dot{u}^{l}=-\kappa^{2}h^{lk}\nabla_{k}\rho=h^{lk}\nabla_{k}(F_{R}(R)R-2F(R)-3\kappa^{2}p)=0.$

If $\rho+p=0$, then NEC satisfied. Furthermore, NEC describes the null
geodesic congruences, and it is convergent in sufficiently small neighborhoods
of every point of spacetime. The physical interpretation of NEC is that
particles following null geodesics will observe that gravity tends locally to
be attractive (or at least not repulsive) when acting on nearby particles also
following null geodesics.

If $\rho+p\neq 0$, then it follows from (28) and (30) that

$\displaystyle\nabla_{j}u_{i}=-\frac{\dot{p}}{\rho+p}h_{ji}=\frac{-\kappa^{2}\dot{\rho}}{F_{R}(R)-2F(R)+4\kappa^{2}\rho}h_{ji}.$
(31)

Since $h^{lk}\nabla_{k}\rho=0$, the spacetime is a GRW-spacetime with a scale
factor $a(t)$. Since every GRW-spacetime had the property

$\nabla_{j}u_{i}=\frac{\dot{a}}{a}h_{ji}.$

comparing with (31) gives

$\frac{4\dot{a}}{a}=\frac{-4\kappa^{2}\dot{\rho}}{F_{R}(R)-2F(R)+4\kappa^{2}\rho},$

which means that

$a=\left(\frac{A}{F_{R}(R)-2F(R)+4\kappa^{2}\rho}\right)^{1/4},$

where $A$ is a constant. Hence we get the following results.

######  Theorem 4.2.

A projectively harmonic perfect fluid spacetime solution of $F(R)$-gravity
theories is either a GRW-spacetime with scale factor

$a=\left(\frac{A}{F_{R}(R)-2F(R)+4\kappa^{2}\rho}\right)^{1/4},$

or $p+\rho=0$. We have discussed previously that $p+\rho$ is represented as
NEC and describes the null energy congruence of the spacetime. Moreover, it
describes the nature of gravity in the neighbourhoods of each spacetime point.
Primarily, it describes the attractive nature of gravity in the locality of
each point.

The profile of scale factor for two models are discussed below. For model-1,
the scale factor reads

$a=\sqrt[4]{-\frac{Ae^{R/\alpha}}{(R-1)\left(e^{R/\alpha}-1\right)}},\,\ R\neq
1,$

and, for model-2, the scale factor reads

$a=\sqrt[4]{-\frac{AR}{(R-1)^{2}}},\,\ R\neq 1.$

From Fig. 7 and 8, one can see that the value of scale factor decreases as the
Ricci scalar increases, which indicates that for higher curvature the scale
factor of our universe is small and vice-verse. Also, it suggests that
curvature reduces as time goes on and scale factor increases.

Figure 7: Evolution of scale factor with respect to $R$ and $\alpha$ with
$A<0$ for model-1. Figure 8: Evolution of scale factor with respect to $R$ and
$A<0$ for model-2.

######  Remark 4.1.

In [22], the authors rightly showed that conharmonic flatness in a spacetime
constrains the standard theory of gravity to a trivial vacuum case with
$R_{ij}=0$ and $R=0$ case. Even in modified gravity theories the situation
does not improve much, it ends up resulting into vanishing effective pressure
and energy density.

We further can add to the findings of the article [22] that instead of
conharmonic flatness if one imposes divergence free conharmonic curvature into
the geometry, then the Ricci scalar $R$ of the spacetime is constant and the
Ricci curvature tensor is of Codazzi type. Therefore, same conclusions as in
Theorem 4.1 and Theorem 4.2 follow from similar arguments.

######  Remark 4.2.

All the results of this section are also true for the standard theory of
gravity governed by Einstein’s field equations (3).

##  5 Discussion

This manuscript has focused on studying the projectively flat spacetime in the
framework of $F(R)$ gravity. We have analyzed our study both analytically and
graphically. We have used the analytical method to construct our formulation,
and graphical analysis has been used to check the stability of two
cosmological toy models for a better understanding, such as
$F(R)=R-\alpha\left(1-e^{-\frac{R}{\alpha}}\right)$ and $F(R)=R-\log(\beta
R)$.

Moreover, we have examined the stability analysis of the cosmological models
using energy conditions. For model-1, we have presented the profiles of energy
conditions in Figures 1,2 and 3. It has been shown that the energy density is
positive throughout the evolution of parameters $R>1$ and $\alpha>0$.
Furthermore, DEC, NEC and WEC have been satisfied, whereas SEC violated. All
of the above profiles of energy conditions are in agreement with the
accelerated expansion of the universe. In addition, the equation of state
parameter $\omega$ takes its value as $-1$, which also indicates the dark
energy era. Nevertheless, these results are also compatible with the
$\Lambda$CDM model. Similarly, for the second model, we have depicted all the
energy conditions in Figures 4, 5, and 6. The results we have obtained for
model-2 are aligned with the results of model-1.

Finally, we have studied the projectively harmonic spacetime in the framework
of $F(R)$ gravity. We have formulated the scale factor in terms of $F(R)$ and
alternately found $\rho+p=0$. Then we have presented the evolution of the
scale factor for two models concerning Ricci scalar $R$ and model parameter in
Fig. 7 and 8, respectively. From those figures, we have observed that for
large values of $R$, scale factor is small, which suggests that the scale
factor of the universe is small in the early phase due to the high curvature,
and later it increases. Alternately, NEC has been satisfied, that is,
$\rho+p=0$. The advantage of this setup is that it does not allow hypothetical
astronomical objects like wormholes (where NEC should be violated) and allow
the Big-Bang singularity.

## Acknowledgments

A.D. and L.T.H. are supported by the FRGS research grant
(FRGS/1/2019/STG06/UM/02/6). S.M. acknowledges Department of Science &
Technology (DST), Govt. of India, New Delhi, for awarding Junior Research
Fellowship (File No. DST/INSPIRE Fellowship/2018/IF180676). We are very much
grateful to the honorable referee and to the editor for the illuminating
suggestions that have significantly improved our work in terms of research
quality, and presentation.

## References

  * [1] H. Weyl, Reine Infinitesimalgeometrie, Math. Z., 2 (1918) 384. Reprinted in Gesammelte Abhandlungen ed. K.Chandrasekharan Vol.2 1–28.
  * [2] H. Weyl, Zur Infinitesimalgeometrie: Einordnung der projektiven und der konformen Auffassung, Göttingen Nachr. (1921) 99. Reprinted in esammelte Abhandlungen ed. K.Chandrasekharan Vol.2 195-207
  * [3] C. Lübbe, A note on the coincidence of the projective and conformal Weyl tensors, arXiv:1301.5659.
  * [4] G. Shabbir, Weyl Projective Curvature Symmetry in FRW $k=0$ Model, arXiv:0705.1882.
  * [5] Hans-Jurgen Schmidt, Fourth order gravity: equations, history, and applications to cosmology, Int. J. Geom. Methods Mod. Phys., 04, (2007) 209.
  * [6] A. A. Starobinsky, A new type of isotropic cosmological models without singularity, Phys. Lett. B, 91 (1980) 99.
  * [7] S. M. Carroll, V. Duvvuri, M. Trodden and M. S. Turner, Is cosmic speed-up due to new gravitational physics?, Phys. Rev. D, 70 (2004) 043528.
  * [8] A. De Felice, S. Tsujikawa, f(R) Theories. Living Rev. Relativ. 13 (2010) 3.
  * [9] A. Raychaudhuri, Relativistic Cosmology Phys. Rev., 98 (1955) 1123.
  * [10] J. Ehlers, A. K. Raychaudhuri and his equations Int. J. Mod. Phys. D, 15 (2006) 1573.
  * [11] S. Nojiri, S. D. Odintsov, Introduction to Modified Gravity and Gravitational Alternative for Dark Energy Int. J. Geom. Methods Mod. Phys., 4 (2007) 115 .
  * [12] P. H. R. S. Moraes, P. K. Sahoo, The simples non-minimal matter-geometry coupling in the $f(R,T)$ cosmology, Eur. Phys. J. C., 77, (2017) 480.
  * [13] R. M. Wald, General relativity (University of Chicago Press, Chicago, 1984.
  * [14] S. Capozzeillo, S. Nojiri, S. D. Odinstov, The role of energy conditions in $f(R)$ cosmology, Phys. Lett. B, 781 (2018) 99.
  * [15] S. Mandal, P.K. Sahoo, J.R.L. Santos, Energy Conditions in f(Q) gravity, Phys. Rev. D, 102 (2020) 024057.
  * [16] S. Arora, J. R. L. Santos, P.K. Sahoo, Phys. Dark Univ., 31 (2021) 100790.
  * [17] S. Mandal, P. K. Sahoo, A Complete Cosmological Scenario in Teleparallel Gravity, Eur. Phys. J. Plus 135(2020) 706.
  * [18] S. Mandal et al., Accelerating Universe in Hybrid and Logarithmic Teleparallel Gravity, Phys. Dark Univ., 28 (2020) 100551.
  * [19] A. De et al., Energy conditions for a (WRS)4 spacetime in $F(R)$ gravity. Eur. Phys. J. Plus, 136 (2021) 218.
  * [20] S. Mandal et al., Almost-Pseudo-Ricci Symmetric FRW Universe with a Dynamic Cosmological Term and Equation of State, Universe, 7 (2021) 205.
  * [21] Planck 2018 results. VI. Cosmological parameters, arXiv:1807.06209.
  * [22] A. Beesham and F. Makhanya, How Extra Symmetries Affect Solutions in General Relativity, Universe, 6 (2020) 170.