# A complete study of conformally flat pseudo-symmetric spacetimes in the
theory of $F(R)$-gravity

Avik De  0000-0001-6475-3085<EMAIL_ADDRESS>Department of Mathematical and
Actuarial Sciences, Universiti Tunku Abdul Rahman, Jalan Sungai Long, 43000
Cheras, Malaysia  Simran Arora  0000-0003-0326-8945<EMAIL_ADDRESS>Department of Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India.  Uday Chand De  0000-0002-8990-4609
<EMAIL_ADDRESS>Department of Pure Mathematics, University of Calcutta, West
Bengal, India  P.K. Sahoo  0000-0003-2130-8832<EMAIL_ADDRESS>pilani.ac.in Department of Mathematics, Birla Institute of Technology and
Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India.

###### Abstract

The motive in this article is twofold. First we investigate the geometrical
structures of a pseudo-symmetric spacetime $(PS)_{4}$ with a timelike vector
under the condition of conformal flatness. We classify it into two possible
types: constant Ricci scalar and closed velocity vector. Then we further study
this spacetime as a solution of $F(R)$-gravity theory and show that the
pressure and energy density of the effective cosmological perfect fluid
separately possess certain typical relations with the geometry and the gravity
sector of the theory. Based on this result, some observational and
cosmological analyses are done for the different $F(R)$-gravity models and the
energy conditions are investigated accordingly, supporting the accelerated
expansion of the universe.

##  I Introduction

To efficiently reason the late time accelerated expansion of the universe
ignoring the assumption of the yet undetected existence of dark energy,
researchers tried to modify the Einstein’s field equations (EFE). One of the
most popular and established such modified theories of gravity is the
$F(R)$-theory. By replacing the Ricci scalar $R$ in the Einstein-Hilbert
action with an arbitrary function $F(R)$ of $R$ and then varying the action
term

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

with respect to $g^{\mu\nu}$ gives the required field equations of
$F(R)$-gravity

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

where $\Box=\nabla^{l}\nabla_{l}$. The stress-energy tensor is obtained from
the matter Lagrangian $L_{m}$ by
$T_{ij}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta g^{ij}}$.
Viability of $F(R)$-gravity models are extensively studied in the literature
[2, 1, 3]. The freedom in developing different functional forms of $F(R)$
gives rise to the challenge of how to confine these numerous possibilities
from theoretical or observational features, even though these theories provide
an alternate way to explain the cosmic speed up. Additional constraints to
$F(R)$ theories may also come by imposing the so-called energy conditions.
Recently, some solutions of $F(R)$-gravity has been given with the spacetime
geometry restricted to almost pseudo-Ricci symmetric [4], weakly Ricci
symmetric [5], generalized Ricci recurrent [6] type, among others. In fact,
imposing these kind of symmetries into the spacetime geometry is always used
in general relativity due to the complications involved to solve the system of
non-linear partial differential equations. Study of Ricci calculus has a rich
history perhaps as old as the theory of general relativity. The study
presumably started with locally symmetric manifolds ($\nabla_{i}R_{hjkl}=0$)
and continues to countless newer structures with the development in Physics,
specially in the theory of gravity. Chaki [7] introduced the notion of pseudo-
symmetric manifold $(PS)_{n}$ as a non-flat Riemannian manifold $M^{n},\,n\geq
2$, whose Riemannian curvature tensor satisfies the relation

$\nabla_{l}R_{hijk}=2E_{l}R_{hijk}+E_{h}R_{lijk}+E_{i}R_{hljk}+E_{j}R_{hilk}+E_{k}R_{hijl}$
(2)

where $E_{i}$ is a non-zero 1-form. A pseudo-Riemannian manifold with
Lorentzian type metric tensor of signature $(-,+,+,+)$ is a pseudo-symmetric
spacetime if $R_{ijkl}$ satisfies (2).

Incidentally, in the literature, there are two different notions of pseudo-
symmetric manifolds, the one we considered, by Chaki [7] and the other by
Deszcz [8, 9], and there is a necessary and sufficient condition for these two
concepts to be equivalent [10]. Researchers studied the curvature properties
and physical implications of $(PS)_{n}$ in different settings, including the
vanishing Weyl conformal curvature tensor. In [11] a $(PS)_{4}$ spacetime was
studied as the solution of the standard theory of gravity and it was shown
that if the Ricci scalar of such a spacetime is nonzero constant, then the
expansion scalar and the acceleration vector of the perfect fluid type matter
content vanish. Very recently, it is showed that a pseudo-symmetric
generalized Robertson-Walker spacetime is a perfect fluid spacetime and a
conformally flat $(PS)_{n}$ is a generalized Robertson-Walker spacetime [12].
This rejuvenates the research interests in $(PS)_{n}$ spacetimes with its
cosmological implications and this time within the realm of the standard
theory of gravity to study its interactions with modified theories of gravity.

The perfect fluid type stress-energy tensor is given by
$T_{ij}=(p+\rho)u_{i}u_{j}+pg_{ij}$, where $p$ and $\rho$ denote the isotropic
pressure and energy density, respectively and $u^{i}$ is the four-velocity
vector of the fluid. In the present discussion we consider a perfect fluid
with a four-velocity vector $u^{i}$ identical to $E^{i}$. The isotropic
pressure $p$ and the energy density $\rho$ of the perfect fluid satisfies a
simplified relation $p=w\rho$, where $w$ is called the equation of state (EoS)
parameter. The evolution of the energy density and the expansion of the
universe is significantly related to the EoS parameter. We can observe
different phases of the universe using EoS [13]. The dust phase, at
$\omega=0$, the radiation-dominated era when $\omega=\frac{1}{3}$, while
$\omega=-1$ corresponds to the vacuum energy, i.e. $\Lambda$CDM model.
Besides, the accelerating phase of the universe, which has recently been in
discussion is portrayed when $\omega<-\frac{1}{3}$. It includes the
quintessence phase when $-1<\omega<0$ and phantom regime $\omega<-1$.

In modified gravity, energy conditions define the attractive behavior,
spacetime geodesic, and casual structure. Further in physical scenario, energy
conditions are also fundamental tools to study black holes and wormholes in
different modified gravity [14, 15, 16, 17] The well-known Raychaudhuri
equations are used to explain an attractive aspect of gravity and positive
energy density. The null energy condition (NEC), weak energy condition (WEC),
dominant energy condition (DEC), and strong energy condition (SEC) are the
four primary conditions. These are obtained from Raychaudhuri equations, which
allow us to study the complete spacetime structure crucial for understanding
cosmological gravitational interactions. The NEC is helpful in discussing the
second law of black hole thermodynamics and its violation causes the universe
to collapse into a Big-Rip singularity [18]. Under modified theories of
gravity, SEC is good at characterizing repulsive or attractive aspects of
gravity. The observed accelerated expansion of the universe is supported by
the violation of SEC. Some of the works on energy conditions in $F(R)$ gravity
are studied in [19, 20]. Also, the study of energy conditions using the
mathematical framework known as cosmography enables us to bound the free
parameters constraining the families of different gravity models compatible
with the accelerated expansion of the universe [21, 22]. The energy conditions
enable us to specify our free parameters, limiting the families of of $F(R)$
models that are compatible with the accelerated expansion of the universe.

The present paper is organized as follows: We consider a 4-dimensional pseudo-
symmetric spacetime and examine the equation of motion in section II after
introducing the necessary ingredients of the present study in the first
section. In section III, we observe the behavior of cosmological parameters
such as density parameter, the equation of state parameter. Further, various
energy conditions are investigated in section IV. The discussions and
conclusion are presented in section V.

Throughout the article we use the notations
$\dot{R}=E^{k}\nabla_{k}R,\,\ddot{R}=E^{k}\nabla_{k}\dot{R},\,\dot{E_{i}}=E^{k}\nabla_{k}E_{i}$.

##  II The equations of motion

In this section we consider a 4-dimensional pseudo-symmetric spacetime, hence
its $(0,4)$\- type Riemannian curvature tensor $R_{ijkl}$ satisfies (2).
Contracting simultaneously $h$, $k$ and $j$, $l$ in (2) we respectively obtain

$\nabla_{l}R_{ij}=2E_{l}R_{ij}+E_{i}R_{lj}+E_{j}R_{il}+E^{k}R_{lijk}+E^{k}R_{ljik}$
(3)

and

$\nabla_{i}R=4E^{l}R_{il}+2RE_{i}.$ (4)

Transvecting (4) by $E^{i}$ we get

$\dot{R}+2R=4E^{i}E^{l}R_{il}.$ (5)

Now, from the contracted Bianchi’s second identity we have the well-known
relation [23]

$\nabla^{m}R_{jklm}=\nabla_{j}R_{kl}-\nabla_{k}R_{jl}.$ (6)

Contracting $k$ and $l$ in (2) and using (6) we obtain

$\nabla_{h}R_{ij}-\nabla_{i}R_{hj}=3E^{k}R_{hijk}+E_{h}R_{ij}-E_{i}R_{hj}.$
(7)

Since we consider a flat conformal curvature tensor, we also have vanishing
conformal divergence in this spacetime, which gives the relation [23]

$6[\nabla_{i}R_{jk}-\nabla_{j}R_{ik}]=g_{jk}\nabla_{i}R-g_{ik}\nabla_{j}R.$
(8)

Using (4), (7) and (8) we get

$9E^{k}R_{hijk}+3E_{h}R_{ij}-3E_{i}R_{hj}\\\
=g_{ij}[RE_{h}+2E^{l}R_{hl}]-g_{hj}[RE_{i}+2E^{l}R_{il}].$ (9)

Now a conformally flat spacetime secures a special format of the Riemannian
curvature tensor expressed in terms of the Ricci curvature tensor and the
Ricci scalar [23]

$6R_{hijk}=(3R_{ij}-Rg_{ij})g_{hk}-(3R_{hj}-Rg_{hj})g_{ik}\\\
+3(g_{ij}R_{hk}-g_{hj}R_{ik}).$ (10)

Transvecting (10) by $E^{k}$ and using (9) we obtain

$E_{h}(9R_{ij}-3Rg_{ij})-E_{i}(9R_{hj}-3Rg_{hj})+9E^{k}(g_{ij}R_{hk}-g_{hj}R_{ik})\\\
+6E_{h}R_{ij}-6E_{i}R_{hj}=g_{ij}[2RE_{h}+4E^{l}R_{lh}]-\\\
g_{hj}[2RE_{i}+4E^{l}R_{il}]$ (11)

which simplifies to

$3E_{h}R_{ij}-RE_{h}g_{ij}-3E_{i}R_{hj}+RE_{i}g_{hj}+g_{ij}E^{l}R_{lh}-g_{hj}E^{l}R_{il}=0.$
(12)

Transvecting (12) by $E^{h}$ we obtain

$3R_{ij}=Rg_{ij}+RE_{i}E_{j}+g_{ij}E^{k}E^{l}R_{kl}-3E_{i}E^{l}R_{jl}-E_{j}E^{l}R_{il}.$
(13)

On the other hand, transvecting (9) by $E^{j}$ we obtain

$E_{h}E^{l}R_{il}=E_{i}E^{l}R_{hl}.$ (14)

And hence by (5) we have

$4E^{l}R_{il}=-4E_{i}E^{h}E^{l}R_{hl}=-E_{i}[\dot{R}+2R].$ (15)

(14) and (15) reduces (13) to

$R_{ij}=\frac{[6R+\dot{R}]}{12}g_{ij}+\frac{[3R+\dot{R}]}{3}E_{i}E_{j}.$ (16)

####  II.0.1 $\dot{R}+3R=0$

In this particular case, the discussed conformally flat spacetime reduces to
an Einstein manifold by the virtue of (16) and thus turns out to be a
spaceform itself. Also, on integration of $\dot{R}+3R=0$ we obtain
$R=R_{0}e^{-3t}$ where by $R_{0}$ we denote the initial value of the Ricci
scalar.

####  II.0.2 $\dot{R}+3R\neq 0$

Covariantly differentiating (16) and then transvecting by $E^{j}$ we get

$4[\dot{R}+3R]\nabla_{k}E_{i}=-3[\nabla_{k}\dot{R}+2\nabla_{k}R]E_{i}-12E^{j}\nabla_{k}R_{ij}$
(17)

Transvecting (17) by $E^{k}$ we obtain

$4(\dot{R}+3R)\dot{E_{i}}=-3(\ddot{R}+2\dot{R})E_{i}-12E^{j}E^{k}\nabla_{k}R_{ij}.$
(18)

From (3) and (15) we obtain

$E^{j}E^{k}\nabla_{k}R_{ij}=(\dot{R}+2R)E_{i}.$ (19)

From (18) and (19) we get

$4(\dot{R}+3R)\dot{E_{i}}=(-3\ddot{R}-18\dot{R}-24R)E_{i}$ (20)

which by $E_{i}\dot{E^{i}}=0$, readily gives us the acceleration vector

$\dot{E_{i}}=0.$ (21)

Whereas using (15) from (4) we obtain the correspondence

$\nabla_{i}R-2RE_{i}=-\dot{R}E_{i}-2RE_{i}$ (22)

which provides

$\nabla_{i}R=-\dot{R}E_{i}.$ (23)

Differentiate (23) covariantly we obtain

$\nabla_{j}\nabla_{i}R=-\dot{R}\nabla_{j}E_{i}-E_{i}\nabla_{j}\dot{R}.$ (24)

Swapping $i$ and $j$ in (24) we obtain

$\nabla_{i}\nabla_{j}R=-\dot{R}\nabla_{i}E_{j}-E_{j}\nabla_{i}\dot{R}.$ (25)

$R$ being a scalar field, this combination produces

$\dot{R}(\nabla_{j}E_{i}-\nabla_{i}E_{j})=E_{j}\nabla_{i}\dot{R}-E_{i}\nabla_{j}\dot{R}.$
(26)

Transvecting the above equation by $E^{j}$ we get

$\nabla_{i}\dot{R}=-E_{i}\ddot{R}$ (27)

which transforms (24) and (25) respectively into

$\nabla_{j}\nabla_{i}R=-\dot{R}\nabla_{j}E_{i}+\ddot{R}E_{i}E_{j}$ (28)
$\nabla_{i}\nabla_{j}R=-\dot{R}\nabla_{i}E_{j}+\ddot{R}E_{i}E_{j}.$ (29)

But these two equations together imply

$\dot{R}[(\nabla_{j}E_{i}-\nabla_{i}E_{j}))]=0$ (30)

So, either $\dot{R}=0$ which by (23) implies a constantor $R$, or
$(\nabla_{j}E_{i}-\nabla_{i}E_{j})=0$, that is, $E$ is a closed 1-form. We
investigate these two cases in two separate subsections.

###  II.1 Case I: Closed 1-form $E$

In this subsection we consider a non-constant $R$. (3) by (5), (15) and (16)
gives

$E^{j}\nabla_{k}R_{ij}=-\left[\frac{13}{12}\dot{R}+\frac{5}{2}R\right]E_{i}E_{k}-\\\
\frac{1}{12}(6R+\dot{R})g_{ik}+E^{l}E^{j}R_{kjil}$ (31)

and, (10) by (5), (15) and (16) gives

$6E^{l}E^{j}R_{kjil}=-(R+\frac{\dot{R}}{2})g_{ik}-(R+\frac{\dot{R}}{2})E_{i}E_{k}.$
(32)

Combining (17),(31),(32), we obtain

$\nabla_{k}E_{i}=\frac{\dot{R}+4R}{2(\dot{R}+3R)}g_{ik}+\frac{3\ddot{R}+20\dot{R}+32R}{4(\dot{R}+3R)}E_{i}E_{k}.$
(33)

Therefore, using (33), from (29) we finally conclude that

$\nabla_{i}\nabla_{j}R=\frac{\dot{R}\ddot{R}+12R\ddot{R}-32R\dot{R}-20\dot{R}^{2}}{4(\dot{R}+3R)}E_{i}E_{j}-\frac{\dot{R}^{2}+4R\dot{R}}{2(\dot{R}+3R)}g_{ij}$
(34)

and

$\Box R=\frac{12\dot{R}^{2}-12R\ddot{R}-\dot{R}\ddot{R}}{4(\dot{R}+3R)}.$ (35)

We summarise these results as:

Theorem 2.1

In a conformally flat $(PS)_{4}$ spacetime with non-constant Ricci scalar, if
$\dot{R}+3R\neq 0$, we have the followings:

  1. 1.

$R_{ij}=\frac{[6R+\dot{R}]}{12}g_{ij}+\frac{[3R+\dot{R}]}{3}E_{i}E_{j}.$

  2. 2.

$E^{l}\nabla_{l}E_{i}=0.$

  3. 3.

$\nabla_{i}R=-\dot{R}E_{i}.$

  4. 4.

$\nabla_{k}E_{i}=\frac{\dot{R}+4R}{2(\dot{R}+3R)}g_{ik}+\frac{3\ddot{R}+20\dot{R}+32R}{4(\dot{R}+3R)}E_{i}E_{k}.$

  5. 5.

$\nabla_{i}\nabla_{j}R=\frac{\dot{R}\ddot{R}+12R\ddot{R}-32R\dot{R}-20\dot{R}^{2}}{4(\dot{R}+3R)}E_{i}E_{j}-\frac{\dot{R}^{2}+4R\dot{R}}{2(\dot{R}+3R)}g_{ij}.$

  6. 6.

$\Box R=\frac{12\dot{R}^{2}-12R\ddot{R}-\dot{R}\ddot{R}}{4(\dot{R}+3R)}.$

Now, for any analytic function $F(R)$ we have

$\nabla_{i}\nabla_{j}F_{R}(R)=F_{RR}(R)\nabla_{i}\nabla_{j}R+F_{RRR}(R)(\nabla_{i}R)(\nabla_{j}R).$
(36)

Using the curvature properties in Theorem 2.1 we conclude that

Theorem 2.2

In a conformally flat $(PS)_{4}$ spacetime with non-constant Ricci scalar, if
$\dot{R}+3R\neq 0$, for any analytic function $F(R)$ we have

$\nabla_{i}\nabla_{j}F_{R}(R)=\left[\left(\frac{\dot{R}\ddot{R}+12R\ddot{R}-32R\dot{R}-20\dot{R}^{2}}{4(\dot{R}+3R)}\right)F_{RR}(R)+\dot{R}^{2}F_{RRR}(R)\right]E_{i}E_{j}-\frac{\dot{R}^{2}+4R\dot{R}}{2(\dot{R}+3R)}F_{RR}(R)g_{ij}$
(37)

and

$\Box
F_{R}(R)=\frac{12\dot{R}^{2}-12R\ddot{R}-\dot{R}\ddot{R}}{4(\dot{R}+3R)}F_{RR}(R)-\dot{R}^{2}F_{RRR}(R)$
(38)

If a conformally flat $(PS)_{4}$ is a solution of the $F(R)$-gravity equations
(1), using Theorem 2.2 we can express

$\kappa
T_{ij}=\left[-\frac{F(R)}{2}+\frac{\dot{R}+6R}{12}F_{R}(R)+\frac{14\dot{R}^{2}+8R\dot{R}-12R\ddot{R}-\dot{R}\ddot{R}}{4(\dot{R}+3R)}F_{RR}(R)-\dot{R}^{2}F_{RRR}(R)\right]g_{ij}\\\
+\left[\frac{\dot{R}+3R}{3}F_{R}(R)-\frac{\dot{R}\ddot{R}+12R\ddot{R}-32R\dot{R}-20\dot{R}^{2}}{4(\dot{R}+3R)}F_{RR}(R)-\dot{R}^{2}F_{RRR}(R)\right]E_{i}E_{j},$
(39)

which implies that $T_{ij}$ denotes a perfect fluid energy-momentum tensor
with isotropic pressure given by

$\kappa
p=-\frac{F(R)}{2}+\frac{\dot{R}+6R}{12}F_{R}(R)+\frac{14\dot{R}^{2}+8R\dot{R}-12R\ddot{R}-\dot{R}\ddot{R}}{4(\dot{R}+3R)}F_{RR}(R)-\dot{R}^{2}F_{RRR}(R)$
(40)

and the energy density given by

$\kappa\rho=\frac{F(R)}{2}+\frac{\dot{R}+2R}{4}F_{R}(R)-\frac{3\dot{R}(\dot{R}+4R)}{2(\dot{R}+3R)}F_{RR}(R).$
(41)

This leads to one of the main results of the present study:

Theorem 2.3

In a conformally flat $(PS)_{4}$ ($\dot{R}+3R\neq 0$) satisfying the theory of
$F(R)$-gravity with a non-rotating matter content is a perfect fluid with
four-velocity vector $E^{i}$; energy density

$\rho=\frac{\frac{F(R)}{2}+\frac{\dot{R}+2R}{4}F_{R}(R)-\frac{3\dot{R}(\dot{R}+4R)}{2(\dot{R}+3R)}F_{RR}(R)}{\kappa}$
(42)

and isotropic pressure

$p=\frac{-\frac{F(R)}{2}+\frac{\dot{R}+6R}{12}F_{R}(R)+\frac{14\dot{R}^{2}+8R\dot{R}-12R\ddot{R}-\dot{R}\ddot{R}}{4(\dot{R}+3R)}F_{RR}(R)-\dot{R}^{2}F_{RRR}(R)}{\kappa}$
(43)

Remark 2.1

In general relativity, $F(R)=R$, so the perfect fluid in that case has
isotropic pressure $p=\frac{\dot{R}}{12\kappa}$ and the energy density
$\rho=\frac{R+\dot{R}/4}{\kappa}$.

It is known that a Lorentzian manifold $M$ of dimension $n>3$ with Ricci
curvature tensor $R_{kl}=Ag_{kl}+Bu_{k}u_{l}$, where $A$ and $B$ are scalar
fields, $u$ is a closed timelike unit vector field and with a conformally
harmonic curvature tensor is a generalized Robertson-Walker (GRW) spacetime
with a Einstein fibre [24]. Moreover, a GRW spacetime $M$ is conformally flat
if and only if its fibre is a space of constant curvature, that is $M$ is
Robertson-Walker (RW) [25]. Using (16) therefore we assert that

Theorem 2.4

A conformally flat $(PS)_{4}$ with a non-constant Ricci scalar $R$ is a
Robertson-Walker spacetime.

###  II.2 Case II: constant $R$

In this case, from (4) we get

$2E^{l}R_{il}=-RE_{i}.$ (44)

For a constant $R$, the field equations of $F(R)$-gravity (1) reduces to

$R_{ij}-\frac{R}{2}g_{ij}=\frac{\kappa}{F_{R}(R)}T^{\text{eff}}_{ij},$ (45)

where $T^{\text{eff}}_{ij}=T_{ij}+\frac{1}{2k}[F(R)-RF_{R}(R)]g_{ij}$.
Transvecting (45) by $E^{j}$ and using (44) we get

$R=\frac{2\kappa\rho-F(R)}{F_{R}(R)}$ (46)

which gives

$\rho=\frac{RF_{R}(R)+F(R)}{2\kappa}.$ (47)

Again, the trace of (45) by the expressions of $R$ and $\rho$ found above
gives us

$p=\frac{RF_{R}(R)-F(R)}{2\kappa}.$ (48)

We state that

Theorem 2.4

In a conformally flat $(PS)_{4}$ (with $\dot{R}=0$) satisfying the theory of
$F(R)$-gravity, the matter content is a perfect fluid with four-velocity
vector $E^{i}$; isotropic pressure $p=\frac{RF_{R}(R)-F(R)}{2\kappa}$ and
energy density $\rho=\frac{RF_{R}(R)+F(R)}{2\kappa}$.
In such case, the EoS $w$ is given by

$w=\frac{p}{\rho}=1-\frac{2F(R)}{F(R)+RF_{R}(R)}.$ (49)

Remark 2.2

In a conformally flat $(PS)_{4}$ with constant $R$ satisfying the field
equations of $F(R)$-gravity with $F_{R}(R)>0$, we readily have the followings:
(i) For a positive Ricci scalar $R>0$, $w>-1$.
(ii)For a negative Ricci scalar $R<0$, $w<-1$, a hypothetical phantom energy
state.

##  III Cosmological Study

The Hubble parameter $H=\frac{\dot{a}}{a}$, where $a$ is the scale factor and
$(\cdot)$ represent the derivative with respect to time, is a cosmological
parameter extensively used to study the rate of expansion of the universe.
Besides, the expansion of scale factor with respect to cosmic time yields the
higher order derivatives. The first four derivatives are known as velocity,
acceleration $(q)$, jerk $(j)$ and snap $(s)$ parameters respectively [26].
This cosmography is a mathematical tool for the description of the universe.
The current values of these parameters with subscript $0$ represent the
evolution of the universe. That is, $q_{0}<0$ signifies an accelerated
expansion, whereas $j_{0}$ and $s_{0}$ distinguish between different
accelerating models.

We first note that the Ricci scalar and its derivatives for a spatially flat
FLRW geometry can be expressed in terms of the deceleration ($q$), jerk ($j$),
and snap ($s$) parameters [28],

$\displaystyle R$ $\displaystyle=-6H^{2}(1-q)$ (50) $\displaystyle\dot{R}$
$\displaystyle=-6H^{3}(j-q-2)$ (51) $\displaystyle\ddot{R}$
$\displaystyle=-6H^{4}(s+q^{2}+8q+6)$ (52)

where

$\displaystyle q$ $\displaystyle=-\frac{1}{H^{2}}\frac{\ddot{a}}{a}$ (53)
$\displaystyle j$ $\displaystyle=\frac{1}{H^{3}}\frac{\dddot{a}}{a}$ (54)
$\displaystyle s$ $\displaystyle=\frac{1}{H^{4}}\frac{\ddddot{a}}{a}$ (55)

Case 2.1.1 : $\dot{R}+3R=0$
As studied in section 2, we have a solution in this case as $R=R_{0}e^{-3t}$.
This form itself describes the space. Since, there is no functional form $f$
in the solution, therefore it can be treated as the General Relativity (GR)
case.

Case 2.1.2 : $\dot{R}+3R\neq 0$

We obtain density and pressure in terms of deceleration parameter, jerk and
snap parameters using eqs.(50),(51),(52) as

$\rho=\frac{1}{\kappa}\left(\frac{F}{2}-\frac{3}{2}F_{R}\left(H^{3}(j-q-2)+2H^{2}(1-q)\right)+\frac{3F_{RR}H^{3}(j-q-2)\left(-6H^{3}(j-q-2)-24H^{2}(1-q)\right)}{-2H^{3}(j-q-2)-6H^{2}(1-q)}\right)$
(56)

$p=\frac{1}{\kappa}\left(-\frac{F}{2}-\frac{1}{2}F_{R}\left(H^{3}(j-q-2)+6H^{2}(1-q)\right)-36F_{RRR}H^{6}(j-q-2)^{2}+\right.\\\
\left.\frac{F_{RR}(-36H^{7}(j-q-2)\left(q^{2}+8q+s+6\right)+504H^{6}(j-q-2)^{2}-432H^{6}(1-q)\left(q^{2}+8q+s+6\right)+288H^{5}(1-q)(j-q-2))}{4\left(-6H^{3}(j-q-2)-18H^{2}(1-q)\right)}\right)$
(57)

###  III.1 Model-I

We assume the functional form of $F(R)=R+\alpha log(\beta R)$ [30], where
$\alpha$ and $\beta$ are constants. It is seen that Logarithm function is
continuous and differential when $\beta R>0$. Therefore, we choose $\beta\neq
0$. For $\alpha=0$, the model reduces to the well motivated GR case. We
consider the present values of $q_{0}=-0.55$, $j_{0}=1$, $s_{0}=-0.35$ and
$H_{0}=67.9$ km/s/Mpc [31, 32]. Then we can have expressions for density
$\rho$ and pressure $p$ as,

$\displaystyle\rho=\frac{1}{4k}\left(\frac{(H(-j+q+2)+2(q-1))(\alpha+6H^{2}(q-1))}{q-1}+2\left(\alpha\log\left(6\beta
H^{2}(q-1)\right)+6H^{2}(q-1)\right)\right.\\\
\left.+\frac{\alpha(-j+q+2)(H(-j+q+2)+4(q-1))}{H(q-1)^{2}(H(-j+q+2)+3(q-1))}\right)$
(58)

and

$\displaystyle
p=-\frac{1}{24\kappa}\left(\frac{2(H(j-q-2)-6q+6)\left(\alpha+6H^{2}(q-1)\right)}{q-1}+12\left(\alpha\log\left(6\beta
H^{2}(q-1)\right)+6H^{2}(q-1)\right)\right.\\\
\left.+\frac{\alpha\left(H^{2}(-(j-q-2))\left(q^{2}+8q+s+6\right)+2H(7j^{2}-14j(q+2)+6q^{3}+49q^{2}+6qs+16q-6s-8)+8(-jq+j+q^{2}+q-2)\right)}{H(q-1)^{2}(H(-j+q+2)+3(q-1))}\right.\\\
\left.+\frac{8\alpha(-j+q+2)^{2}}{(q-1)^{3}}\right)$ (59)

Since $F(R)=R+\alpha log(\beta R)$ is defined for $\beta R>0$ and $R$ is
negative with the value of $q_{0}$, therefore we consider $\beta<0$. Also,
$\alpha=0$ corresponds to $F(R)=R$, which reduces to the case of GR. We obtain
the bounds of $\alpha$ through the behavior of energy density $\rho$, which is
always positive, as shown in Fig. 1(a) and negative pressure depicted in Fig.
1(b). It is readily seen that $\beta$ contributes more to the energy density
and $\rho$ increases for more negative $\beta$ and decreases for less negative
$\beta$. Further, the equation of state parameter(EoS) is the relationship
between $p$ and $\rho$ i.e., $\omega=\frac{p}{\rho}$. It is used to classify
various epochs, such as decelerated to accelerated phases of the universe. So,
$\omega=0$ represent the matter-dominated phase, $\omega=\frac{1}{3}$ shows
the radiation whereas the accelerated phases are depicted as the quintessence
when $-1<\omega<0$, the $\Lambda$CDM as $\omega=-1$ and $\omega<-1$ is the
phantom era. As a result, the EoS parameter is thought to be a good choice for
comparing our models to $\Lambda$CDM. So, the EoS parameter is used to
constrain model parameters to investigate various energy conditions. The
behavior of the EoS parameter for our considered model is shown in Fig. 1(c).
The plot of EoS shows a transition from positive to negative supporting an
acceleration in the universe and lies in the quintessence phase $-1<\omega<0$.

(a) The density parameter.

(b) The pressure parameter.

(c) The EoS parameter.

Figure 1: The behavior of the density, pressure and the EoS parameters
considering $F(R)=R+\alpha log(\beta R)$ for Case 2.1.2 with $-1\leq\alpha\leq
50000$ and $-40\leq\beta\leq-10$.

Case II : $\dot{R}=0$

In this case $R$ is constant. We obtain $\rho$ and $p$ for this case in
section II. Using equations (47), (48), (50) and the functional form
$F(R)=R+\alpha log(\beta R)$, where $\alpha$, $\beta$ are constants. Again,
for $\alpha=0$, the model reduces to GR. Further from (66), we get $p=0$ for
$\alpha=0$, and according to it $\rho$ is constant which is instant the case
of cosmological constant. We assumed $\alpha$ different from zero to analyze
the derived model to study different cases rather than the cosmological
constant.

$\rho=\frac{\alpha+\alpha\log\left(6\beta
H^{2}(q-1)\right)+12H^{2}(q-1)}{2\kappa}.$ (60)

$p=\frac{\alpha-\alpha\log\left(6\beta H^{2}(q-1)\right)}{2\kappa}.$ (61)

The Eos parameter $\omega=\frac{p}{\rho}$ is obtained as

$\omega=1-\frac{2\left(\alpha\log\left(-6\beta
H^{2}(1-q)\right)-6H^{2}(1-q)\right)}{\alpha\log\left(-6\beta
H^{2}(1-q)\right)-6H^{2}(1-q)\left(1-\frac{\alpha}{6H^{2}(1-q)}\right)-6H^{2}(1-q)}$
(62)

We consider the present values $q_{0}=-0.55$, $j_{0}=1$, $s_{0}=-0.35$ and
$H_{0}=67.9$ km/s/Mpc to study the behavior of different parameters. It is
observed from Figs. 2(a) and 2(b) that density is positive whereas pressure is
negative in the range $40000\leq\alpha\leq 70000$ and $-2\leq\beta\leq-1$. As
a result, the EoS parameter is used to constrain model parameters to
investigate various energy conditions. In Fig. 2(c), it is seen that the EoS
parameter lies in the quintessence phase, and $\omega$ is close to -1, which
shows an accelerated phase of the universe.

(a) The density parameter.

(b) The pressure parameter.

(c) The EoS parameter.

Figure 2: The behavior of the density, pressure and the EoS parameters
considering $F(R)=R+\alpha log(\beta R)$ for Case II with $40000\leq\alpha\leq
70000$ and $-2\leq\beta\leq-1$.

###  III.2 Model-II

Case 2.1.2 : $\dot{R}+3R\neq 0$

We assume the functional form of $F(R)=\gamma Exp(\frac{\zeta}{R})-R$ [30],
where $\gamma$ and $\zeta$ are constants. In the limiting case when $\gamma$
tends to $0$, the model reduces to the case of GR. We consider the present
values of $q_{0}=-0.55$, $j_{0}=1$, $s_{0}=-0.35$ and $H_{0}=67.9$ km/s/Mpc
[31, 32]. Then the expressions for density $\rho$ and pressure $p$ read as,

$\rho=\frac{\frac{1}{2}\gamma
e^{\frac{\zeta}{6H^{2}(q-1)}}-3H^{2}(q-1)-\frac{\gamma\zeta(-j+q+2)(H(-j+q+2)+4(q-1))e^{\frac{\zeta}{6H^{2}(q-1)}}\left(\zeta+12H^{2}(q-1)\right)}{144H^{5}(q-1)^{4}(H(-j+q+2)+3(q-1))}+\frac{(H(j-q-2)-2q+2)\left(36H^{4}(q-1)^{2}+\gamma\zeta
e^{\frac{\zeta}{6H^{2}(q-1)}}\right)}{24H^{2}(q-1)^{2}}}{\kappa}$ (63)

$p=\frac{1}{\kappa}\left(-\frac{1}{2}\gamma
e^{\frac{\zeta}{6H^{2}(q-1)}}+3H^{2}(q-1)+\frac{\begin{multlined}\gamma\zeta
e^{\frac{\zeta}{6H^{2}(q-1)}}\left(\zeta+12H^{2}(q-1)\right)(H^{2}(-(j-q-2))(q(q+8)+s+6)+2H(7j^{2}\\\
-14j(q+2)+q(q(6q+49)+6s+16)-6s-8)+8(-jq+j+q^{2}+q-2))\end{multlined}\gamma\zeta
e^{\frac{\zeta}{6H^{2}(q-1)}}\left(\zeta+12H^{2}(q-1)\right)(H^{2}(-(j-q-2))(q(q+8)+s+6)+2H(7j^{2}\\\
-14j(q+2)+q(q(6q+49)+6s+16)-6s-8)+8(-jq+j+q^{2}+q-2))}{864H^{5}(q-1)^{4}(H(-j+q+2)+3(q-1))}\right.\\\
\left.+\frac{(H(j-q-2)-6q+6)\left(36H^{4}(q-1)^{2}+\gamma\zeta
e^{\frac{\zeta}{6H^{2}(q-1)}}\right)}{72H^{2}(q-1)^{2}}+\frac{\gamma\zeta(-j+q+2)^{2}e^{\frac{\zeta}{6H^{2}(q-1)}}\left(\zeta^{2}+216H^{4}(q-1)^{2}+36\zeta
H^{2}(q-1)\right)}{1296H^{6}(q-1)^{6}}\right)$ (64)

(a) The density parameter.

(b) The pressure parameter.

(c) The EoS parameter.

Figure 3: The behavior of the density, pressure and the EoS parameters
considering $F(R)=\gamma Exp(\zeta/R)-R$ for Case 2.1.2 with $9\times
10^{6}\leq\gamma\leq 5\times 10^{7}$ and $0.1\leq\zeta\leq 0.5$.

We obtained the bound for $\gamma$ and $\zeta$ through the density parameter
in 3(a), pressure 3(b) and the EoS parameter in 3(c). According to the
observations, the density behaves positively, whereas the EoS parameter should
be negative and close to -1. We notice that $\omega\approx-1$ from lower
values in, revealing an unexpected behavior consistent with a phantom era for
dark energy. Furthermore, when dealing with a standard explanation of gravity
plus an essential phantom field, a phantom period for dark energy could lead
to a problematic description for the universe because the DEC would be
violated [33]. Moreover, the equation of state parameter is bound to be
greater than -1, implying that DEC is validated [34]. However, given the
existing experimental bounds for $\omega=-1.03^{+0.03}_{-0.03}$ [32] show that
a phantom description of dark energy is still possible.

Case II : $\dot{R}=0$

In this case $R$ is constant. We obtain $\rho$ and $p$ for this case in
section II. Using equations (47), (48), (50) and the functional form
$F(R)=\gamma Exp(\frac{\zeta}{R})-R$, where $\gamma$, $\zeta$ are constants.
So, the density and pressure are as follows.

$\rho=\frac{6H^{2}(q-1)\left(\gamma\zeta e^{6\zeta H^{2}(q-1)}-2\right)+\gamma
e^{6\zeta H^{2}(q-1)}}{2\kappa}.$ (65)

$p=\frac{\gamma e^{6\zeta H^{2}(q-1)}\left(6\zeta
H^{2}(q-1)-1\right)}{2\kappa}.$ (66)

The Eos parameter $\omega=\frac{p}{\rho}$ is obtained as

$\omega=\frac{\gamma e^{6\zeta H^{2}(q-1)}\left(6\zeta
H^{2}(q-1)-1\right)}{6H^{2}(q-1)\left(\gamma\zeta e^{6\zeta
H^{2}(q-1)}-2\right)+\gamma e^{6\zeta H^{2}(q-1)}}.$ (67)

We consider the present values $q_{0}=-0.55$, $j_{0}=1$, $s_{0}=-0.35$ and
$H_{0}=67.9$ km/s/Mpc to study the behavior of different parameters. It is
observed from Figs. 4(a) and 4(b) that density is positive whereas pressure is
negative in the range $7\times 10^{5}\leq\gamma\leq 8\times 10^{5}$ and
$0.0001\leq\zeta\leq 0.0002$. Here, we have considered the bounds of
parameters large due to higher powers of $H$. In Fig. 4(c), it is seen that
the EoS parameter lies in the quintessence phase, with a transition from
positive to negative and is close to -1, depicting an accelerated phase of the
universe.

(a) The density parameter.

(b) The pressure parameter.

(c) The EoS parameter.

Figure 4: The behavior of the density, pressure and the EoS parameters
considering $F(R)=\gamma Exp(\frac{\zeta}{R})-R$ for Case II with $7\times
10^{5}\leq\gamma\leq 8\times 10^{5}$ and $0.0001\leq\zeta\leq 0.0002$.

##  IV Energy conditions

Raychaudhuri equations, which play a fundamental part in any consideration of
the congruence of null and timelike geodesics with the criterion that only
gravity is attractive but also the energy density is positive, give rise to
the concept of energy conditions. Let $u^{i}$ be the tangent vector field to a
congruence of timelike geodesics in a spacetime manifold furnished with a
metric $g_{ij}$ to illustrate this concept. The Raychaudhuri equations
describe the temporal evolution of the expansion scalar. For a brief review
see[27].

$\frac{d\theta}{d\tau}=-\frac{\theta^{2}}{3}+\omega_{ij}\omega^{ij}-\sigma_{ij}\sigma^{ij}-R_{ij}u^{i}u^{j},$
(68)

$\frac{d\theta}{d\tau}=-\frac{\theta^{2}}{2}+\omega_{ij}\omega^{ij}-\sigma_{ij}\sigma^{ij}-R_{ij}k^{i}k^{j}.$
(69)

where $\omega_{ij}$, $\sigma_{ij}$, $u^{i}$ and $k^{i}$ are the rotation,
shear tensor, timelike and null tangent vector fields respectively. Since, we
are considering the case $\omega_{ij}=0$ and neglecting the small distortions,
we obtained $\theta=-\tau R_{ij}u^{i}u^{j}$ and $\theta=-\tau
R_{ij}k^{i}k^{j}$. According to attractive nature of gravity,
$R_{ij}u^{i}u^{j}\geq 0$ and $R_{ij}k^{i}k^{j}\geq 0$. Further according to
Einstein field equations, we get $(T_{ij}-\frac{1}{2}g_{ij}T)u^{i}u^{j}\geq 0$
and $(T_{ij}-\frac{1}{2}g_{ij}T)k^{i}k^{j}\geq 0$. In case of perfect fluid
matter distribution, the energy conditions are defined as

  * •

Null energy condition (NEC): $\rho+p\geq 0$

  * •

Weak energy condition (WEC): $\rho+p\geq 0$ and $\rho\geq 0$

  * •

Strong energy conditions (SEC): $\rho+p\geq 0$ and $\rho+3p\geq 0$

  * •

Dominant energy condition (DEC): $\rho\pm p\geq 0$ and $\rho\geq 0$.

The concept of energy boundaries in modified theories of gravity can be
expanded with the premise that the whole cosmic matter behaves like a perfect
fluid due to the utterly geometric nature of the Raychaudhuri equations. So,
we can consider $p^{eff}$ as $p$ and $\rho^{eff}$ as $\rho$. The classical
energy conditions and their cosmological implications were studied extensively
in $F(R)$-gravity theory in FLRW setting [28, 29].

###  IV.1 Model I

Case 2.1.2: $\dot{R}+3R\neq 0$

The energy conditions for the functional form of $F(R)=R+\alpha log(\beta R)$
are the following.

$\displaystyle\rho+p=\frac{H^{2}}{2\kappa}\left(\frac{\begin{multlined}3F_{RR}H(H^{2}(j-q-2)(q(q+8)+s+6)-4H(2j^{2}-4j(q+2)+q(q(3q+23)+3s+2)-3s-10)-\\\
16(q-1)(j-q-2))\end{multlined}3F_{RR}H(H^{2}(j-q-2)(q(q+8)+s+6)-4H(2j^{2}-4j(q+2)+q(q(3q+23)+3s+2)-3s-10)-\\\
16(q-1)(j-q-2))}{H(j-q-2)-3q+3}+\right.\\\
\left.4F_{R}(H(-j+q+2)+3(q-1))-72F_{RRR}H^{4}(-j+q+2)^{2}\right)$ (70)

$\displaystyle\rho-p=\frac{1}{2\kappa}\left(\frac{\begin{multlined}3F_{RR}(H^{2}(-(j-q-2))(q(q+8)+s+6)+4H(5j^{2}-10j(q+2)+q(q(3q+26)+3s+14)-3s+2)-\\\
32(q-1)(j-q-2))\end{multlined}3F_{RR}(H^{2}(-(j-q-2))(q(q+8)+s+6)+4H(5j^{2}-10j(q+2)+q(q(3q+26)+3s+14)-3s+2)-\\\
32(q-1)(j-q-2))}{H(j-q-2)-3q+3}+\right.\\\
\left.2F+H^{3}(2F_{R}(-j+q+2)+72F_{RRR}H^{3}(-j+q+2)^{2}\right)$ (71)

$\displaystyle\rho+3p=\frac{1}{2\kappa}\left(\frac{9F_{RR}H^{4}(H(j-q-2)(q(q+8)+s+6)-12(j^{2}-2j(q+2)+q(q(q+8)+s+2)-s-2))}{H(j-q-2)-3q+3}\right.\\\
\left.-2F+6F_{R}H^{2}(H(-j+q+2)+4(q-1))-216F_{RRR}H^{6}(-j+q+2)^{2}\right)$
(72)

As, WEC ( $\rho\geq 0$ ) depends on the observational values of $q_{0}$ and
$j_{0}$, i.e. independent of $s_{0}$, so we obtain the restrictions on the
model parameters $\alpha$ and $\beta$. Since no accurate measurement of the
snap parameter has been reported, we can concentrate on the WEC requirement in
confronting energy conditions. The condition reduces to

$0.5\alpha\log(-42876.8\beta)-4.42758\alpha+168430\geq 0.$ (73)

It is seen that $\beta$ can take negative values and $\alpha$ can vary
negative to positive satisfying the above condition.

(a) $\rho$ \+ p.

(b) $\rho$ \- p.

(c) $\rho$ \+ 3 p.

Figure 5: The behavior of energy conditions considering $F(R)=R+\alpha
log(\beta R)$ for Case 2.1.2 with $-1\leq\alpha\leq 50000$ and
$-40\leq\beta\leq-10$.

There may be a violation of energy conditions in some instances without the
system becoming unacceptable. The success of the hypothesis of inflation and
recent observations of cosmic acceleration point to a violation of SEC in the
universe. So, we can observe from the plots 5(a), 5(b) that the WEC, DEC are
satisfying their conditions (showing positive behavior) in the range of
$-1\leq\alpha\leq 50000$ and $-40\leq\beta\leq-10$. On the other, SEC in Fig.
5(c) exhibits negative behavior violating its condition and hence supports an
accelerated expansion of the universe.

Case II: $\dot{R}=0$

The energy conditions for this case should satisfy the following:

$\rho+p=\frac{\alpha+\alpha\log\left(6\beta
H^{2}(q-1)\right)+12H^{2}(q-1)}{2\kappa}+\frac{\alpha-\alpha\log\left(6\beta
H^{2}(q-1)\right)}{2\kappa}\geq 0$ (74)

$\rho-p=\frac{\alpha+\alpha\log\left(6\beta
H^{2}(q-1)\right)+12H^{2}(q-1)}{2\kappa}-\frac{\alpha-\alpha\log\left(6\beta
H^{2}(q-1)\right)}{2\kappa}\geq 0$ (75)

$\rho+3p=\frac{\alpha+\alpha\log\left(6\beta
H^{2}(q-1)\right)+12H^{2}(q-1)}{2\kappa}-\frac{\alpha-\alpha\log\left(6\beta
H^{2}(q-1)\right)}{2\kappa}\geq 0$ (76)

(a) $\rho$ \+ p.

(b) $\rho$ \- p.

(c) $\rho$ \+ 3 p.

Figure 6: The behavior of energy conditions considering $F(R)=R+\alpha
log(\beta R)$ for Case II with $40000\leq\alpha\leq 70000$ and
$-2\leq\beta\leq-1$ .

WEC is the combination of $\rho+p\geq 0$ and the positive density. The
condition on WEC gives

$0.5(\alpha\log(-42876.8\beta)+\alpha-85753.6)\geq 0.$ (77)

which confronts the negative value of $\beta$ and the variation of $\alpha$.
Here we can observe from Figs. 6(a) & 6(b) that WEC, DEC are validating their
respective conditions with $40000\leq\alpha\leq 70000$ and
$-2\leq\beta\leq-1$. On the other hand, SEC in Fig. 6(c) shows negative
behavior violating its condition in the same range. Hence, the violation of
the SEC supports an accelerated expansion of the universe.

###  IV.2 Model-II

Case 2.1.2: $\dot{R}+3R\neq 0$

The energy conditions for the functional form $F(R)=\gamma
Exp(\frac{\zeta}{R})-R$ are obtained and the behavior is shown in the
following plots.

(a) $\rho$ \+ p.

(b) $\rho$ \- p.

(c) $\rho$ \+ 3 p.

Figure 7: The behavior of energy conditions considering $F(R)=\gamma
Exp(\zeta/R)-R$ for Case 2.1.2 with $9\times 10^{6}\leq\gamma\leq 5\times
10^{7}$ and $0.1\leq\zeta\leq 0.5$.

Furthermore, the plot 7 reveals that DEC in 7(b) is satisfying its conditions
whereas WEC in 7(a) is partially obeyed(i.e. satisfying only $\rho>0$).
However, we can also see that SEC in 7(c) is again violated confirming that
our universe experiences an accelerated expansion of the universe. Hence, the
WEC violation along with the positive density acts as a scalar-tensor gravity
model in this case [35], and such a violation can be understood as natural
contributions from quantum processes to classical gravity [36].

Case II: $\dot{R}=0$

The energy conditions for the functional form $F(R)=\gamma
Exp(\frac{\zeta}{R})-R$ are obtained. Since equations are little big. We have
shown the behavior in the following plots.

(a) $\rho$ \+ p.

(b) $\rho$ \- p.

(c) $\rho$ \+ 3 p.

Figure 8: The behavior of energy conditions considering $F(R)=\gamma
Exp(\zeta/R)-R$ for Case II with $7\times 10^{5}\leq\gamma\leq 8\times 10^{5}$
and $0.0001\leq\zeta\leq 0.0002$.

Here we can observe from Figs. 8(a) & 8(b) that WEC, DEC are validating their
respective conditions with $7\times 10^{5}\leq\gamma\leq 8\times 10^{5}$ and
$0.0001\leq\zeta\leq 0.0002$. We have considered the bounds of parameters
large due to higher powers of $H$. On the other hand, SEC in Fig. 8(c) shows
negative behavior violating its condition in the same range. Hence, the
violation of the SEC supports an accelerated expansion of the universe.

##  V Discussion

There is a need to investigate some fundamental questions about the
unexplained dark energy, the underlying acceleration of the universe. What is
causing this sudden acceleration? Is it due to the addition of a dark energy
component such as scalar field or quintessence, or is it due to a modification
in Einstein’s General Theory of Relativity? In this work, we attempt to
investigate the accelerated expansion in the realm of modified $F(R)$ gravity.
We assume the functional form of $F(R)=R+\alpha log(\beta R)$ where $\alpha$
and $\beta$ are constants and $F(R)=\gamma Exp(\zeta/R)-R$ with $\gamma$ and
$\zeta$ are constants. The spacetime is assumed to be pseudo-symmetric, and a
detailed analysis of it has been carried out, providing some insight into the
behavior of density and equation of state parameters. The density shows
positive behavior, whereas the EoS parameter shows a negative behavior, lying
in the quintessence phase in model I whereas a phantom regime in model II,
indicating the acceleration in the universe. Further energy conditions are
investigated using the current estimated values of deceleration, jerk, and
snap parameters. We obtain the bounds of model parameters $\alpha$, $\beta$,
$\gamma$ and $\zeta$ showing the viability of the considered models. However,
it is seen that all the energy conditions NEC, WEC, and DEC are satisfying
their derived conditions while SEC shows a violation supporting the
accelerated expansion of the universe in the Model I. Whereas in the Model II,
DEC is met but WEC and SEC are not obeyed, depicting the behavior similar to
the scalar tensor gravity models and the accelerating phase, respectively.

## Acknowledgments

SA acknowledges CSIR, New Delhi, India for JRF.

## References

  * [1] I. De Martino et al., Constraining $f(R)$ gravity with Planck data on galaxy cluster profiles, MNRAS, 442 (2014) 921.
  * [2] S. Capozziello, V. F. Cardone, V. Salzano, Cosmography of $f(R)$ gravity, Phys. Rev. D, 78 (2008) 063504.
  * [3] T. P. Sotiriou, V. Faraoni, $f(R)$ theories of gravity, Rev. Mod. Phys., 82 (2010) 451.
  * [4] A. De and L. T. How. Almost pseudo-Ricci symmetric spacetime solutions in F(R)-gravity. Gen Relativ Gravit, 53, 5 (2021).
  * [5] Avik. De et al., Energy conditions for a $(WRS)_{4}$ spacetime in $F(R)$ gravity, Eur. Phys. J. Plus, 136 (2021) 218.
  * [6] Avik De et al., A conformally flat generalized Ricci recurrent spacetime in $F(R)$ gravity, Phys. Scr. 96 (2021) 085001.
  * [7] M. C. Chaki and U. C. De,On pseudo symmetric spaces, Acta Math.Hungarica 54 (1989), 185-190.
  * [8] F. Defever, R. Deszcz, L. Verstraelen and L. Vrancken, On pseudosymmetric spacetimes, Journal of Mathematical Physics, 35, (1994) 5908.
  * [9] S. Haesen and L. Verstraelen, Classification of the pseudosymmetric spacetimes, Journal of Mathematical Physics, 45, (2004) 2343.
  * [10] A.A. Shaikh, R. Deszcz, M. Hotlos, J. Jelwicki and H. Kundu, On pseudosymmetric manifolds, Publicationes Mathematicae, 86/3-4 (2015), 433-456.
  * [11] B. Chaki, S. Guha and A. Konar, On Chaki Pseudo Symmetric Perfect Fluid Space-Time, Bulgarian Journal of Physics, 26 (1999) 204-209.
  * [12] P. Zhao, U. C. De, B. Unal and K. De, Sufficient conditions for a pseudosymmetric spacetime to be a perfect fluid spacetime, arXiv:2105.03597.
  * [13] T. D. Saini, S. Raychaudhury, V. Sahni, A. A. Starobinsky, Reconstructing the Cosmic equation of state from supernovae distances, Phys. Rev. Lett., 85 (2000) 1162.
  * [14] K. Bamba, M. Ilyas, M.Z. Bhatti, Z. Yousaf, Energy conditions in modified $f(G)$ gravity, Gen. Relativ. Gravit, 49 (2017) 112.
  * [15] Z. Yousaf, M. Ilyas, M.Z. ul-Haq Bhatti, Static spherical wormhole models in $f(R,T)$ gravity, Eur. Phys. J. Plus, 132 268 (2017).
  * [16] M. Halilsoy, A. Ovgun, Particle acceleration by static black holes in a model of $f(R)$ gravity, Canadian Journal of Physics, 95(11) (2017) 1037.
  * [17] A. Ovgun, M. Halilsoy, Existense of traversable wormholes in the spherical stellar systems, Astrophys. Space Sci., 361, 214 (2016).
  * [18] S. Carroll, Spacetime and Geometry: An Introduction to General Relativity (Addison Wesley, Boston, 2004).
  * [19] K. Atazadeh et al., Energy conditions in $f(R)$ gravity and Brans-Dicke theories, Int. J. Mod. Phys. D, 18 (2009) 1101.
  * [20] J. Santos, M. J. Reboucas, J. S. Alcaniz, Energy conditions on a class of $f(R)$ gravity, Int. J. Mod. Phys. D, 19 (2010) 1315.
  * [21] S. Mandal, P.K. Sahoo, J.R.L. Santos, Energy conditions in $f(Q)$ gravity, Phys. Rev. D, 102(2) (2020) 024057.
  * [22] S. Arora, J.R.L. Santos, P.K. Sahoo, Constraining $f(Q,T)$ gravity from energy conditions, Phys. Dark Universe, 31 (2021) 100790.
  * [23] S. W. Hawking and G. F. R. Ellis, The Large Scale Structure of Space-Time, (Cambridge University Press, 1973) doi:10.1017/CBO9780511524646.
  * [24] C. A. Mantica, L.G. Molinari, U.C. De, A condition for a perfect fluid space-time to be a generalized Robertson-Walker space-time, J. Math. Phys., 57 (2016) 022508.
  * [25] M. Brozos-Vazquez, E. Garcia-Rio, and R. Vazquez-Lorenzo, Some remarks on locally conformally flat static spacetimes, J. Math. Phys., 46 (2005) 022501.
  * [26] M. Visser, Jerk, snap and the cosmological equation of state, Class. Quantum Grav., 21 (2004) 2603.
  * [27] S. Kar, S. Sengupta, The Raychaudhuri equations: A brief review, Pramana, 69 (2007) 49-76.
  * [28] J. Santos et al., Energy conditions in $f(R)$ gravity, Phys. Rev. D, 76 (2007) 083513.
  * [29] J. Wang et al., Energy conditions and stability in generalized gravity with arbitrary coupling between matter and geometry, Phys. Lett. B, 689 (2010) 133.
  * [30] Z. Gironés et al., Cosmological data analysis of $f(R)$ gravity, JCAP, 11 (2010) 004.
  * [31] S. Capozziello, R. D Agostino, O. Luongo, Extended gravity cosmography, Int. J. Mod. Phys. D, 28 (2019) 1930016.
  * [32] N. Aghanim et al., Planck 2018 results-VI. Cosmological parameters, Astronomy & Astrophysics, 641 (2020) A6.
  * [33] S. M. Carroll, A. De Felice, M. Trodden, Can we be tricked into thinking that $\omega$ is less than -1?, Phys.Rev.D, 71 (2005) 023525.
  * [34] S. M. Carroll, M. Hoffman, M. Trodden, Can the dark energy equation of state parameter $\omega$ be less than -1?, Phys.Rev.D, 68 (2003) 023509.
  * [35] A. W. Whinnett, D. F. Torres, A new strong-field effect in scalar tensor gravity: spontaneous violation of the energy conditions, The Astrophys. J, 603 (2004) L133.
  * [36] G. Calcagni, Classical and Quantum Cosmology, Springer Switzerland, 2017, doi: 10.1007/978-3-319-41127-9


# Electron-ion heating partition in imbalanced solar-wind turbulence

Jonathan Squire Physics Department, University of Otago, Dunedin 9010, New
Zealand<EMAIL_ADDRESS>Romain Meyrand Physics Department,
University of Otago, Dunedin 9010, New Zealand  Matthew W. Kunz Department of
Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ
08544, USA  Princeton Plasma Physics Laboratory, PO Box 451, Princeton, NJ
08543, USA

###### Abstract

A likely candidate mechanism to heat the solar corona and solar wind is low-
frequency “Alfvénic” turbulence sourced by magnetic fluctuations near the
solar surface. Depending on its properties, such turbulence can heat different
species via different mechanisms, and the comparison of theoretical
predictions to observed temperatures, wind speeds, anisotropies, and their
variation with heliocentric radius provides a sensitive test of this physics.
Here we explore the importance of normalized cross helicity, or imbalance, for
controlling solar-wind heating, since it is a key parameter of magnetized
turbulence and varies systematically with wind speed and radius. Based on a
hybrid-kinetic simulation in which the forcing’s imbalance decreases with
time—a crude model for a plasma parcel entrained in the outflowing wind—we
demonstrate how significant changes to the turbulence and heating result from
the “helicity barrier” effect. Its dissolution at low imbalance causes its
characteristic features—strong perpendicular ion heating with a steep
“transition-range” drop in electromagnetic fluctuation spectra—to disappear,
driving more energy into electrons and parallel ion heat, and halting the
emission of ion-scale waves. These predictions seem to agree with a diverse
array of solar-wind observations, offering to explain a variety of complex
correlations and features within a single theoretical framework.

††journal: Astrophys. J. Lett.††software: Pegasus++

##  1 Introduction

The solar corona and its extended outflow, the solar wind, provide us with an
unparalleled laboratory for studying the physics of magnetized collisionless
plasmas. Decades of observations have revealed a highly complex system filled
with electromagnetic fluctuations across a vast range of scales, whose
properties (e.g., total power or spectra) are correlated in surprising and
nontrivial ways with the plasma’s flow speed, temperatures, anisotropies, and
elemental abundances (Marsch, 2006; Horbury et al., 2012; Bruno & Carbone,
2013). These correlations, as well as the extended heating at large
heliocentric distances needed to explain the high speed of fast-wind streams
(Parker, 1965), suggest that the solar wind is shaped by both properties of
its low-coronal source and turbulent heating at larger altitudes. Of
particular interest are the decades of observations that hint at the role
played by the fluctuations’ _imbalance_ (i.e., normalized cross helicity), a
key parameter in the theory of magnetized turbulence (Dobrowolny et al., 1980;
Schekochihin, 2022) that is observed to correlate with wind speed $U$ and
decrease with increasing heliocentric radius $R$ (e.g., Roberts et al., 1987;
Marsch, 2006; D’Amicis et al., 2021).

In this Letter, we argue for the importance of imbalance in shaping turbulence
and heating in the low-$\beta$ solar wind. We focus on the physics of the
“helicity barrier” (Meyrand et al., 2021), which dramatically alters
imbalanced turbulence in $\beta\lesssim 1$ plasmas by restricting the
turbulent energy cascade at perpendicular scales $\lambda_{\perp}$ below the
ion gyroradius $\rho_{i}$. Previous kinetic simulations have shown how
numerous features of helicity-barrier-mediated turbulence match measurements
of the low-$\beta$ solar wind well (Squire et al., 2022; hereafter S+22),
including properties of the steep “transition-range” drop in electromagnetic
field spectra around $\rho_{i}$ scales and the ion velocity distribution
function (VDF; see, e.g., Duan et al., 2021; Bowen et al., 2022, 2023).

We study how the turbulence evolves as imbalance is slowly decreased in time,
transitioning between the highly imbalanced and balanced regimes. This physics
is of interest for two reasons: first, it is a crude model for the expected
response of the turbulent heating as the plasma flies outwards from the Sun
and its imbalance decreases; second, it captures how the helicity barrier
becomes weaker and eventually disappears at low imbalance, thus probing its
robustness. The results, which are based on a kinetic simulation of forced
Alfvénic turbulence, can explain the observed correlation of wind speed with
ion temperatures, as well as the switch from negative to positive correlation
of wind speed with electron temperature at increasing $R$ (e.g., Burlaga &
Ogilvie, 1973; Marsch et al., 1989; Shi et al., 2023; Bandyopadhyay et al.,
2023). Correlations with other properties such as the spectral transition
range, proton VDFs, and plasma-wave/instability activity seem to explain a
diverse array of observations within a single theoretical framework.

##  2 Formulation of the problem and method of solution

###  2.1 The helicity barrier

The helicity barrier arises due to the conservation of a “generalized
helicity” $\mathcal{H}$ in low-$\beta$ gyrokinetics (Schekochihin et al.,
2019; Meyrand et al., 2021). At large scales $\lambda_{\perp}\gg\rho_{i}$,
$\mathcal{H}$ is the cross helicity, which, along with the energy, cascades
forward to smaller scales; at small scales $\lambda_{\perp}\ll\rho_{i}$,
$\mathcal{H}$ is a form of magnetic helicity, which cannot cascade to smaller
$\lambda_{\perp}$ without violating energy conservation (Fjørtoft, 1953).
Thus, the sub-$\rho_{i}$ cascade must have $\mathcal{H}=0$, restricting its
energy flux to be at most the “balanced portion” of that injected,
$\varepsilon-\varepsilon_{H}$, where $\varepsilon$ and $\varepsilon_{H}$ are
the large-scale injection rates of energy and helicity, respectively. The
remaining energy input,
$\varepsilon-(\varepsilon-\varepsilon_{H})=\varepsilon_{H}$, is “trapped” at
$\lambda_{\perp}\lesssim\rho_{i}$ scales, causing the turbulence to reach
large amplitudes and (through critical balance) small parallel scales
approaching the ion inertial length $d_{i}$. The fluctuations then take on the
character of oblique ion-cyclotron waves (ICWs) with frequencies $\omega$ that
approach the ion gyrofrequency $\Omega_{i}$. This activates strong wave-
particle interactions and quasi-linear resonant heating (Kennel & Engelmann,
1966), which can absorb the remaining energy input, allowing the turbulence to
saturate (Li et al., 1999). It also causes the plasma to become unstable and
emit small-scale parallel ICWs (Chandran et al., 2010b). By diverting a large
fraction of the energy flux that would have reached electron scales, and thus
electron heat, into ions instead, the helicity barrier causes large ion-to-
electron heating ratios
$Q_{i}/Q_{e}\simeq\varepsilon_{H}/(\varepsilon-\varepsilon_{H})$ in
anisotropic low-$\beta$ turbulence.

###  2.2 Numerical method

We use the hybrid-kinetic method implemented in the Pegasus++ code (Kunz et
al., 2014; Arzamasskiy et al., 2023). This approach treats the ion (proton)
dynamics fully kinetically using a particle-in-cell (PIC) approach, while the
electrons constitute a massless, neutralizing, isothermal fluid. The ion
macro-particle positions ($\bm{r}$) and velocities ($\bm{v}$) are drawn from
an initially Maxwellian $f_{i}(\mbox{\boldmath{$v$}})$ and evolved via

$\frac{{\rm d}\bm{r}}{{\rm d}t}=\bm{v},\quad\frac{{\rm d}\bm{v}}{{\rm
d}t}=\frac{e}{m_{i}}\left[\bm{E}(\bm{r},t)+\frac{\bm{v}}{c}\,\mbox{\boldmath{$\times$}}\,\bm{B}(\bm{r},t)\right]+\frac{1}{m_{i}}\bm{F}_{\perp}^{U},$
(1)

where $E$ and $B$ are the electric and magnetic fields and
$\bm{F}_{\perp}^{U}$ stirs turbulence by injecting incompressible motions
perpendicular (“$\perp$”) to a mean (“guide”) magnetic field $\bm{B}_{0}$. The
magnetic field satisfies a modified version of Faraday’s law,

$\frac{\partial\bm{B}}{\partial
t}=-c\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,(\bm{E}+\bm{F}_{\perp}^{B})+\eta_{4}\nabla^{4}\bm{B},$
(2)

where $\bm{F}_{\perp}^{B}$ forces solenoidal magnetic fluctuations
perpendicular to $\bm{B}_{0}$ and $\eta_{4}$ is a hyper-resistivity that
absorbs small-scale magnetic energy. The electric field is

$\bm{E}=-\frac{\bm{u}}{c}\,\mbox{\boldmath{$\times$}}\,\bm{B}-\frac{T_{e}}{en}\mbox{\boldmath{$\nabla$}}{n}+\frac{(\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\bm{B})\,\mbox{\boldmath{$\times$}}\,\bm{B}}{4\pi
en},$ (3)

where $n$ is the ion (and electron) density and $\bm{u}$ is the ion flow
velocity, both of which are computed via a weighted sum of the macro-particles
in the relevant region of space; $T_{e}$ is the electron temperature, which is
a parameter of the model; and $e$, $m_{i}$, and $c$ are the electron/ion
charge, the ion mass, and the speed of light, respectively. We also define the
Alfvén speed $v_{\rm A}\equiv B_{0}/\sqrt{4\pi nm_{i}}$, $\Omega_{i}\equiv
eB_{0}/(m_{i}c)$, the ion thermal speed $v_{\rm th}$, and the “Elsässer”
fields $\bm{z}^{\pm}\equiv\bm{u}_{\perp}\pm\bm{B}_{\perp}/\sqrt{4\pi
nm_{i}}\equiv\bm{u}_{\perp}\pm\bm{b}_{\perp}$. Volume averages are denoted by
$\langle\,\cdot\,\rangle$ and $f_{i}(w_{\perp},w_{\|})$ is the gyro-averaged
VDF, with $w_{\perp}$ and $w_{\|}$ the particle velocities perpendicular and
parallel to the local magnetic field in the frame of the plasma (i.e.,
$\bm{w}\equiv\bm{v}-\bm{u}$).

###  2.3 Problem set up

Our basic simulation set up follows S+22. The simulation represents a small
co-moving patch of plasma, capturing realistic solar-wind turbulence
amplitudes and anisotropies around $k_{\perp}\rho_{i}\sim 1$ scales, where
$k_{\perp}\sim 1/\lambda_{\perp}$ denotes the perpendicular wavenumber. The
domain is Cartesian and periodic, with coordinates $\\{x,y,z\\}$, and is
elongated along $\bm{B}_{0}=-B_{0}\hat{\mbox{\boldmath{$z$}}}$, with
$L_{z}=6L_{\perp}$ and $L_{\perp}=67.5d_{i}$ ($d_{i}=v_{\rm A}/\Omega_{i}$).
With these parameters, the box shape approximately matches the (statistical)
shape of turbulent eddies measured at similar scales in the solar wind (Chen
et al., 2016; S+22). The grid resolution is $N_{\perp}^{2}\times
N_{z}=392^{2}\times 2352$, so that the smallest resolved scales are
$k_{\perp,{\rm max}}d_{i}\simeq\pi N_{\perp}d_{i}/L_{\perp}\approx 18$. We use
$N_{\rm ppc}=216$ ion-macro-particles per cell. The hyper-resistivity was
increased from $\eta_{4}\approx 2.4\times 10^{-5}d_{i}^{4}\Omega_{i}$ to
$\eta_{4}=5\times 10^{-5}d_{i}^{4}\Omega_{i}$ at $t=4.3\tau_{\rm A}$ in
response to the strengthening kinetic-range cascade.

The simulation is initialized at $t=0$ using the final snapshot from S+22, and
thus represents saturated highly imbalanced turbulence with strong
perpendicular ion heating and a helicity barrier, similar to turbulence
observed in the near-Sun solar wind by PSP (Bowen et al., 2023). The initial
ion temperature is such that $\beta_{i}=8\pi\langle nT_{i}\rangle/\langle
B^{2}\rangle\approx 0.33$; electrons are isothermal with
$\beta_{e0}=8\pi\langle n\rangle T_{e}/B_{0}^{2}=0.3$. The simulation is run
for an additional ${\approx}18\tau_{\rm A}$, where $\tau_{\rm A}=L_{z}/v_{\rm
A}$ is the outer-scale Alfvén time, which is also comparable to the turnover
time of the turbulence $\tau_{\rm turb}\sim L_{\perp}/u_{{\rm rms}}$ due to
our choice of forcing parameters (see below; $u_{{\rm rms}}=\langle
u^{2}\rangle^{1/2}$ is the root-mean-square velocity). Because the ions heat
up, $\beta_{i}$ increases over the course of the simulation to ${\simeq}0.45$
and $\rho_{i}=\sqrt{\beta_{i}}d_{i}$ changes modestly. The domain resolves
scales between $k_{\perp 0}\rho_{i0}\equiv 2\pi\rho_{i0}/L_{\perp}\approx
0.05$ and $k_{\perp{\rm max}}\rho_{i0}\approx 10$, where
$\rho_{i0}=\rho_{i}(t=0)$.

###  2.4 Decreasing imbalance with distance from the Sun

A novel feature of this work is the forcing, which changes from imbalanced to
balanced over the simulation, heuristically mimicking the radial evolution of
the solar wind. We inject energy and cross helicity at the rates $\varepsilon$
and $\varepsilon_{H}$, respectively, with forcing functions
$\bm{F}_{\perp}^{U}$ and $\bm{F}_{\perp}^{B}$ that are intended to capture the
effect of stirring due to turbulent eddies above the box scale. These consist
of random combinations of large-scale Fourier modes with wavenumbers $k_{j}$
satisfying $2\pi/L_{j}\leq k_{j}\leq 4\pi/L_{j}$ for $j=\\{x,y,z\\}$. They are
computed as $\bm{F}_{\perp}^{U}=f^{U}\bm{F}_{0}$ and
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\bm{F}_{\perp}^{B}=f^{B}\bm{F}_{0}$,
where $\bm{F}_{0}$ is divergence-free, perpendicular to $\bm{B}_{0}$, and
evolved in time via an Ornstein–Uhlenbeck process with correlation time
$\tau_{\rm A}/2$. We fix $\varepsilon$ and $\varepsilon_{H}$ at each time step
by adjusting $f^{U}$ and $f^{B}$ so that
$n\bm{u}\,\mbox{\boldmath{$\cdot$}}\,\bm{F}_{\perp}^{U}$ and
$\bm{B}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\bm{F}_{\perp}^{B}$
take the values needed to inject the required energies into $\bm{z}^{\pm}$.
$\bm{F}_{\perp}^{B}$ is computed from $\bm{F}_{0}$ via a fast Fourier
transform and used in the standard Pegasus++ constrained-transport algorithm
to evolve $\bm{B}$, ensuring that
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,\bm{B}=0$ to machine
precision.

We fix $\varepsilon=C_{\rm A}m_{i}n(L_{\perp}/L_{z})^{2}v_{\rm
A}^{2}V/\tau_{\rm A}\approx 37m_{i}nv_{\rm A}^{2}\Omega_{i}d_{i}^{3}$, where
$C_{\rm A}=0.29$ is the Kolmogorov constant and $V$ is the simulation volume;
the factor $(L_{\perp}/L_{z})^{2}v_{\rm A}^{2}/\tau_{\rm A}$ guarantees
critically balanced fluctuations at the outer scale with $u_{\rm
rms}\sim(L_{\perp}/L_{z})v_{\rm A}$ and $\tau_{\rm turb}\sim\tau_{\rm A}$. In
contrast, $\varepsilon_{H}$ decreases in time during the simulation so that
the “injection imbalance” $\varepsilon_{H}/\varepsilon$ starts at $0.9$ (as in
S+22) then decreases linearly in time at a rate of $0.1$ every $1.5\tau_{\rm
A}$, reaching $\varepsilon_{H}=0$ (balanced forcing) at $t=13.5\tau_{\rm A}$,
where it remains for the rest of the simulation (see Fig. 1). Because the
cross helicity and energy are approximately conserved at $k_{\perp}\rho_{i}\ll
1$ scales, this evolution is intended to mimic crudely the effect of larger
scales becoming more balanced with radius, thereby driving the smaller scales
with decreasing $\varepsilon_{H}/\varepsilon$. However, in the actual solar
wind, the timescale over which the imbalance decreases is comparable to both
the turbulent decay timescale and the expansion timescale $\tau_{\rm
exp}\approx R/U$ (Roberts et al., 1987; Meyrand et al., 2023). By keeping
$\varepsilon$ fixed, we effectively assume that the turbulent decay and
expansion are slow compared to the simulation’s duration, which is appropriate
since $\tau_{\rm exp}/\tau_{\rm A}\simeq 1300\,(B/80~{}{\rm nT})(R/35~{}{\rm
R}_{\odot})(U/350~{}{\rm km~{}s}^{-1})^{-1}$ (see e.g., Bale et al. 2019).
However, this also implies that our adopted imbalance-decrease timescale of
${\approx}15\tau_{\rm A}$ is too short. This trade-off is unavoidable given
the extreme computational expense of kinetic-turbulence simulations, but it
should be kept mind that the simulation can only qualitatively capture
realistic features of this transition in the solar wind.

##  3 Results

Figure 1: (Top) Time evolution of $\varepsilon_{H}/\varepsilon$ (purple),
energy imbalance $\sigma_{c}$ (orange), and total fluctuation energy per unit
volume $E_{\perp}(t)$ (blue). (Bottom) Electron heating rate
$Q_{e}=\varepsilon_{\eta}$ (green) and ion heating rate $Q_{i}$ (red), with
simple helicity-barrier expectations based on the forcing shown with dotted
lines.

Figure 2: (_Top_) Perpendicular spectra of the electric field
($\mathcal{E}_{E_{\perp}}$; solid lines) and magnetic field
($\mathcal{E}_{B_{\perp}}$; dashed lines) at a selection of times. The sub-
panel illustrates the local spectral slopes $\alpha$, computed via a fit over
range $\pm\log_{10}(k_{\perp})\approx 0.1$ around each point; the dotted lines
indicate $k_{\perp}\rho_{i}=1$ and the representative power laws
$k_{\perp}^{-5/3}$ and $k_{\perp}^{-0.8}$. (_Bottom_) Evolution of the break
scale $k_{\perp}^{*}$ with imbalance $1-\sigma_{c}$, with the line color
indicating the time as in the top panel. The dotted line shows the empirical
scaling $k_{\perp}^{*}\rho_{i}=2(1-\sigma_{c})^{1/4}$, which provides a good
fit to the simulation until $t\approx 11\tau_{\rm A}$.

###  3.1 Evolution of the helicity barrier

The time evolution as the injection imbalance decreases is shown in Fig. 1.
Properties of the saturated imbalanced turbulence, from which the system is
initialized, are discussed in S+22; a portion of this phase is shown at $t<0$
for comparison. The imbalance of the fluctuations,

$\sigma_{c}\equiv\frac{\langle|\bm{z}^{+}|^{2}\rangle-\langle|\bm{z}^{-}|^{2}\rangle}{\langle|\bm{z}^{+}|^{2}\rangle+\langle|\bm{z}^{-}|^{2}\rangle}=\frac{2\langle\bm{u}_{\perp}\,\mbox{\boldmath{$\cdot$}}\,\bm{b}_{\perp}\rangle}{\langle|\bm{u}_{\perp}|^{2}\rangle+\langle|\bm{b}_{\perp}|^{2}\rangle},$
(4)

starts at $\sigma_{c}\approx 0.98$ for $t<0$, then decreases as the system
adjusts to the changing forcing ($\varepsilon_{H}/\varepsilon$). We halt the
simulation once $\sigma_{c}\simeq 0.4$, a value similar to that in the solar
wind around $1~{}{\rm au}$; as shown below, the heating properties have
changed significantly by this point even though $\sigma_{c}$ is still nonzero.
The blue line shows the fluctuation energy density $E_{\perp}=\langle
nm_{i}(|\bm{u}_{\perp}|^{2}+|\bm{b}_{\perp}|^{2})\rangle/2$, which decreases
over the simulation, as expected because the turbulence can dissipate more
effectively at lower $\varepsilon_{H}/\varepsilon$.

The lower panel of Fig. 1 compares the measured heating rates, $Q_{i}$ and
$\varepsilon_{\eta}$, to helicity-barrier predictions. Here
$\varepsilon_{\eta}$, the hyper-resistive dissipation, is taken as a proxy for
the electron-heating rate $Q_{e}$, which assumes that ions cannot interact
with $k_{\perp}\rho_{i}\gg 1$ eddies and that there is minimal electron
heating via ICWs or other $k_{\perp}\rho_{i}\lesssim 1$ fields (Chandran et
al., 2010b). For most of the simulation, $\varepsilon_{\eta}$ (green line)
tracks well the helicity-barrier prediction
$Q_{e}=\varepsilon-\varepsilon_{H}$ (dotted-green line), aside from a delay of
${\sim}\tau_{\rm A}$ as the injected energy cascades towards smaller scales.
$Q_{i}$ increases initially because the energy $E_{\perp}(t)$ decreases while
the sub-$\rho_{i}$ energy flux is fixed by the helicity barrier. Its large
value results from the unrealistically fast decrease of
$\varepsilon_{H}/\varepsilon$: in reality $E_{\perp}(t)$ would change more
slowly, implying that $Q_{i}=\varepsilon-
Q_{e}-\partial_{t}E_{\perp}\approx\varepsilon_{H}-\partial_{t}E_{\perp}$
should more closely track $\varepsilon_{H}$. As $\varepsilon_{H}/\varepsilon$
nears zero ($t/\tau_{\rm A}\gtrsim 10$) and the helicity barrier erodes, the
heating departs from $Q_{e}\approx\varepsilon-\varepsilon_{H}$ and approaches
that found in kinetic simulations of balanced turbulence at $\beta\approx 0.3$
(Kawazura et al., 2019; Cerri et al., 2021), in which ion heating via Landau
damping and stochastic heating absorbs a portion of the energy flux before it
reaches the smallest scales. At around the same time, $Q_{i}$ drops
significantly.

Figure 2 shows the evolution of the perpendicular spectra
$\mathcal{E}(k_{\perp})$ of $\bm{E}_{\perp}$ and $\bm{B}_{\perp}$. The flatter
sub-$\rho_{i}$ range in the former helps to highlight the double-kinked
“transition-range” power law; for $k_{\perp}\rho_{i}\lesssim 1.5$, including
in the transition range,
$\mathcal{E}_{E_{\perp}}\approx\mathcal{E}_{B_{\perp}}$. As the turbulence
becomes balanced, the spectral break smoothly moves towards smaller scales,
creating a less pronounced transition range that is narrower in $k_{\perp}$
and less steep, eventually transitioning directly into the kinetic range with
$\mathcal{E}_{E_{\perp}}\sim k_{\perp}^{-0.8}$. At late times, $\delta
B_{\perp}$ is resistively damped by the larger resistivity and the magnetic
spectrum does not exhibit a clear ${\sim}k_{\perp}^{-2.8}$ kinetic range.

The lower panel of Fig. 2 provides the evolution of the break scale
$k^{*}_{\perp}$, obtained by fitting a broken power-law to
$\mathcal{E}_{B_{\perp}}(k_{\perp})$.111 The functional form of the fit is
$\mathcal{E}_{\rm
fit}=[(k_{\perp}/k_{\perp}^{*})^{n\alpha_{1}}+(k_{\perp}/k_{\perp}^{*})^{n\alpha_{2}}]^{-1/n}$,
where $\alpha_{1}$ ($\alpha_{2}$) is the power-law index for
$k_{\perp}<k_{\perp}^{*}$ ($k_{\perp}>k_{\perp}^{*}$), and $n$ controls the
break’s sharpness; we fit $\mathcal{E}_{B_{\perp}}(k_{\perp})$ in the range
$0.14<k_{\perp}\rho_{i}<2$. The measured proportionality between
$k_{\perp}^{*}\rho_{i}$ and $(1-\sigma_{c})^{1/4}$ varies with fitting
choices, but the power-law exponent ($1/4$) is robust. We find a clear power-
law dependence on the energy imbalance,
$k_{\perp}^{*}\rho_{i}\propto(1-\sigma_{c})^{1/4}\sim(z^{-}_{\rm
rms}/z^{+}_{\rm rms})^{1/2}$, for $\sigma_{c}\gtrsim 0.8$ ($t\lesssim
11\tau_{\rm A}$) when the barrier is active
($\varepsilon_{H}/\varepsilon\gtrsim 0$). Although this scaling remains
unexplained theoretically, it matches the results of low-$\beta$ gyrokinetic
simulations for various choices of $\varepsilon_{H}/\varepsilon$ (Meyrand et
al. 2021, figure 7), suggesting it is a robust consequence of the helicity
barrier and independent of the energy-dissipation mechanism (gyrokinetics is
ignorant of ICW kinetic physics). The persistence of the transition-range
drop, as well as the evolution of $\varepsilon_{\eta}$ and $Q_{i}$ in Fig. 2,
provide good evidence that a helicity barrier actively mediates the turbulence
until $t\approx 11\tau_{\rm A}$, at which point it transitions towards the
balanced regime.

Figure 3: (_Left_) Time evolution of the perpendicular diffusion energy
coefficient $D_{\perp\perp}^{E}$, showing the transition from the quasi-linear
expectation $D_{\perp\perp}^{E}\propto w_{\perp}^{2}$ to a flatter profile
around $t\gtrsim 11\tau_{\rm A}$. The inset shows the normalized perpendicular
(orange) and parallel (purple) ion heating rates, illustrating a sharp drop in
$Q_{\perp,i}$ associated with the change in $D_{\perp\perp}^{E}$, while
$Q_{\|,i}$ remains almost constant. (_Right_) Differential heating rates,
$v_{\rm th0}^{2}\varepsilon^{-1}\partial^{2}Q_{\perp,\|,i}/\partial
w_{\perp}\partial w_{\|}$, computed from
$\langle\mbox{\boldmath{$E$}}_{\perp,\|}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$w$}}_{\perp,\|}\rangle$.
The $Q_{\perp}$ profile (panels b and c) exhibits resonant structure at early
times (top); the shaded region shows $w_{\|{\rm res}}$ for oblique ICWs with
$d_{i}^{-1}<k_{\|}<2d_{i}^{-1}$, around where there is a sharp dropoff in wave
power (Fig. 4). By late times it drops significantly and an ill-defined peak
around $w_{\perp}\simeq v_{\rm th}$ appears, reminiscent of stochastic
heating. The $Q_{\|}$ profile (panels d and e) maintains its structure,
although spreads out modestly (accounting partially for its lower values).
Figure 4: (a)–(c) 2D spectra of $\delta B_{\perp}$,
$\mathcal{E}_{B_{\perp}}(k_{\perp},k_{\|})$, at $t\approx\tau_{\rm A}$ (left),
$t\approx 11\tau_{\rm A}$ (around the time of the heating transition; middle),
and $t\approx 18\tau_{\rm A}$ (right). The dotted white lines show the
critical-balance condition $k_{\|}v_{\rm A}=2k_{\perp}z^{+}_{\rm
rms}(k_{\perp}/k_{\perp 0})^{-1/3}$. (d) $E_{\rm ICW}$ for different ICW
obliquities, which is the energy contained in the shaded regions of panel (b)
(see text). (e)–(f) Frequency spectra at different times, using the same
colors as Fig. 3a. A sharp drop in high-frequency $\delta n$ fluctuations,
taken as a signature of oblique ICW activity, matches the drop in
$D_{\perp\perp}^{E}$. In contrast, $\delta B_{\perp}$ fluctuations, present in
both parallel and oblique ICWs, drop less abruptly and at later times.

###  3.2 Consequences for ion heating

This drop in ion heating associated with the balanced transition is diagnosed
in Fig. 3. We show the perpendicular energy-diffusion coefficient
$D_{\perp\perp}^{E}$ (panel a), which is computed from the measured evolution
of $f_{i}$ via

$D_{\perp\perp}^{E}=\left(\frac{\partial f_{i}}{\partial
e_{\perp}}\right)^{-1}\int_{0}^{e_{\perp}}{\rm
d}e_{\perp}^{\prime}\frac{\partial f_{i}(e_{\perp}^{\prime})}{\partial t},$
(5)

where $e_{\perp}\equiv w_{\perp}^{2}/2$ and $f_{i}(e_{\perp})=\int{\rm
d}w_{\|}w_{\perp}\,f_{i}(w_{\perp},w_{\|})$. Equation (5) is taken directly
from the $D_{\perp\perp}^{E}$ definition, $\partial f_{i}/\partial
t=\partial/\partial e_{\perp}(D_{\perp\perp}^{E}\partial f_{i}/\partial
e_{\perp})$ (Vasquez et al., 2020), averaging over the $w_{\|}$ variation of
$f_{i}$ and assuming that the heating is predominantly perpendicular, as
appropriate for $t\lesssim 12\tau_{\rm A}$ (Fig. 3a inset). In quasi-linear
theory, waves strongly scatter “resonant” particles with $w_{\|}=w_{\|{\rm
res}}\equiv\omega(\bm{k})/k_{\|}-n\Omega_{i}/k_{\|}$ for $n\in\mathbb{Z}$,
flattening $f_{i}$ along contours of constant energy in the wave’s frame (the
“resonance contours”; we consider only $n=1$ resonances). When parallel waves
dominate the spectrum, $D_{\perp\perp}^{E}\propto w_{\perp}^{2}$, because
$w_{\|{\rm res}}$ and various wave-intensity/polarization factors are all
independent of $w_{\perp}$ (Kennel & Engelmann, 1966). Although the wave
spectrum here is more complex, involving a mix of oblique and parallel ICW
modes (see below), this prediction is nonetheless well satisfied until
$t\approx 11\tau_{\rm A}$.222The effect of obliquity is to flatten
$D_{\perp\perp}^{E}$ at larger $w_{\perp}$, particularly for shorter-
wavelength waves (Isenberg & Vasquez, 2011). The highest-$k_{\|}$ power here
is predominantly parallel (Fig. 4) so these effects are likely minor. After
this, $D_{\perp\perp}^{E}$ drops and flattens, causing a large drop in the
perpendicular ion heating $Q_{\perp,i}$, even while the parallel heating
$Q_{\|,i}$ remains almost constant (see inset). At later times, the form of
$D_{\perp\perp}^{E}$ does not clearly indicate a particular heating mechanism,
but is plausibly consistent with stochastic heating (Chandran et al., 2010a):
computing the $D_{\perp\perp}^{E}$ predicted for stochastic heating using the
method of Cerri et al. (2021) gives a similar shape (not shown). As a
complementary analysis, panels (b) and (c) display
$\partial^{2}Q_{\perp,i}/\partial w_{\perp}\partial w_{\|}$, the differential
perpendicular heating per unit velocity, which is computed directly from
$\langle\bm{E}_{\perp}\,\mbox{\boldmath{$\cdot$}}\,\bm{w}_{\perp}\rangle$
evaluated along particle trajectories (the so-called field-particle
correlation technique; Klein & Howes, 2016; Arzamasskiy et al., 2019). At
early times a resonant feature is centered around the $w_{\|}$ that resonates
with the smallest-($d_{i}$)-scale oblique ICWs with significant power (the
shaded region shows $w_{\|{\rm res}}$ for $d_{i}^{-1}\lesssim k_{\|}\lesssim
2d_{i}^{-1}$; see Fig. 4); at later times, the magnitude drops significantly
into a diffuse peak around $w_{\perp}\simeq v_{\rm th}$, consistent with
stochastic heating. In contrast, $\partial^{2}Q_{\|,i}/\partial
w_{\perp}\partial w_{\|}$ (the differential parallel heating; panels d and e),
maintains a form consistent with Landau damping of kinetic Alfvén waves (KAWs)
throughout the simulation.

In Fig. 4 we provide evidence that the shut off in ion heating occurs because
the turbulence amplitude decreases to the point where it can no longer drive
quasi-linear heating by oblique ICWs. In the top panels, we show 2D spectra of
$\delta B_{\perp}$ fluctuations, $\mathcal{E}_{B_{\perp}}(k_{\perp},k_{\|})$,
which are computed by filtering each field into $k_{\perp}$ bins, then
interpolating these onto the exact magnetic-field lines to take a $k_{\|}$
spectrum (see Methods in S+22). Energy is concentrated at $k_{\perp}>k_{\|}$
(the turbulence) and in a bump at $k_{\perp}\ll k_{\|}\simeq 0.8d_{i}$
(parallel ICWs). Examining the time evolution from left to right, we see that
the energy migrates to lower $k_{\|}$ with time (it moves downwards), with the
dotted white lines indicating how this follows the critical-balance scaling
$k_{\|}v_{\rm A}=Ak_{\perp}z^{+}_{\rm rms}(k_{\perp}/k_{\perp 0})^{-1/3}$ (the
coefficient $A=2$ is chosen to align with the “peak” on the cone and is
consistent across time). A corollary is that larger-amplitude turbulence feeds
more power into smaller-$k_{\|}$ oblique ICWs with smaller $w_{\|{\rm res}}$,
driving more quasi-linear heating.

Figure 5: (_Top_) Contours of the ion VDF $f_{i}$ at $t\approx 5\tau_{\rm A}$
(white lines and colormap) and the $t\approx 18\tau_{\rm A}$ (red lines).
Dotted vertical lines highlight $w_{\|}=0$ and $w_{\|}=v_{\rm A}$. The red
contour levels are chosen to emphasize changes in $f_{i}$ and do not
correspond to the white contours (see text). (_Bottom_) $f_{i}$ versus
distance $s$ along the oblique-ICW resonant contour shown with the black line
in the top panel ($s=0$ at $w_{\|}=0$). Colors are as in Fig. 3a. The thin-
black line is flat ($\mathcal{G}_{\rm res}[f_{i}]=0$), while the dashed line
shows an $f_{i}$ that is independent of $w_{\perp}$ (as expected from
stochastic heating).

We quantify this effect in panel (d) by computing the energy of ICWs ($E_{\rm
ICW}$), defined as fluctuations with $0.7\leq k_{\|}d_{i}\leq 2$, as a
function of wavevector obliquity $\theta=\tan^{-1}(k_{\perp}/k_{z})$. The
minimum wavenumber $k_{\|}d_{i}=0.7$ is chosen to capture fluctuations with
$w_{\|{\rm res}}\lesssim v_{\rm th}$, which interact with the VDF core
(Isenberg & Vasquez, 2011), while the angle ranges capture quasi-parallel
($0<\theta<45^{\circ}$), moderately oblique ($45^{\circ}<\theta<60^{\circ}$),
and highly oblique ($60^{\circ}<\theta<75^{\circ}$) populations; these ranges
are indicated by the shaded regions in panel (b). Despite $E_{\perp}$
decreasing continuously, the energies of the oblique-ICW populations increase
slightly for $t\lesssim 12\tau_{\rm A}$ before dropping rapidly. The
similarity of this functional form with that of $Q_{\perp,i}(t)$ provides good
evidence for oblique ICWs being the primary driver of heating. The energy of
parallel ICWs ($\theta\lesssim 45^{\circ}$) follows a similar trend, but they
are not driven directly by the turbulence, which has little power at
$k_{\|}\gtrsim k_{\perp}$. Instead they arise because oblique-ICW quasi-linear
heating causes $f_{i}$ to increase along the resonance contours of parallel
ICWs, thus emitting waves (causing instability; Kennel & Wong, 1967; Chandran
et al., 2010b). Their energy is driven and undamped, building up substantially
for $t\lesssim 12\tau_{\rm A}$, then decreasing at later times as $f_{i}$
changes shape and renders parallel ICWs stable. Further evidence for this
scenario is presented in panels (e) and (f), which show frequency spectra
$\mathcal{E}(\omega)$ of the fluctuations in density $\delta n$ (e) and
$\delta B_{\perp}$ (f), computed from high-cadence time data at $100$ points,
windowed around different times matching Fig. 3 (we highlight the ICW
frequency range where $k_{\|}d_{i}\gtrsim 0.7$,
$0.3\lesssim\omega/\Omega_{i}\lesssim 1$). While both parallel and oblique
ICWs involve $\delta B_{\perp}$ fluctuations, only oblique modes involve
$\delta n$. Indeed, we see a significant drop in $\mathcal{E}_{\delta n}$ at
the same time as the drop in $Q_{\perp}$, while $\mathcal{E}_{\delta
B_{\perp}}$ drops more smoothly and only for $t\gtrsim 13\tau_{\rm A}$ (orange
to red lines). The behavior of $\mathcal{E}_{\delta B_{\perp}}$ is consistent
with the $0<\theta\leq 45^{\circ}$ $E_{\rm ICW}$ line in panel (d), which only
drops below its initial value at $t\approx 13.5\tau_{\rm A}$.

The interpretation described above predicts that $Q_{i}$ scales with the
oblique-ICW power, which itself scales with $k_{\|,{\rm max}}d_{i}$, where
$k_{\|,{\rm max}}$ corresponds to the smallest parallel scale accessible by
the imbalanced portion of the cascade. An independent estimate of $k_{\|,{\rm
max}}d_{i}$ can be obtained from the spectrum using $k_{\|,{\rm max}}v_{\rm
A}\sim k_{\perp}^{*}z^{+}_{\rm rms}(k_{\perp}^{*}/k_{\perp 0})^{-1/3}$, where
$k_{\perp}^{*}\rho_{i}\propto(z^{-}_{\rm rms}/z^{+}_{\rm rms})^{1/2}$ is the
break scale (Fig. 2). This gives $k_{\|,{\rm
max}}d_{i}\propto\beta^{-1/2}(k_{\perp 0}\rho_{i})^{1/3}(z^{+}_{\rm
rms})^{2/3}(z^{-}_{\rm rms})^{1/3}/v_{\rm A}$. In the simulation, $(z^{+}_{\rm
rms})^{2}z^{-}_{\rm rms}$ increases modestly for $t\lesssim 11\tau_{\rm A}$
before dropping significantly (not shown), remaining well correlated with
$Q_{i}$ and $E_{\rm ICW}$ at all times, including during pre-saturated phase
of S+22. This agreement demonstrates the internal consistency of the results
across diverse diagnostics ($E_{\rm ICW}$, $Q_{i}$, and $k_{\perp}^{*}$ via
$\mathcal{E}_{B_{\perp}}(k_{\perp})$), and could prove useful for developing a
closure model for helicity-barrier-mediated turbulence.

The effect of these dynamics on $f_{i}$ is illustrated in Fig. 5, whose top
panel superimposes isocontours of $f_{i}(w_{\perp},w_{\|})$ at $t\approx
18\tau_{\rm A}$ over $f_{i}(w_{\perp},w_{\|})$ at $t\approx 5\tau_{\rm A}$.
Particles are scattered to render $f_{i}$ constant along the resonant
contours, which are computed from $\mathcal{G}_{\rm res}[f_{i}]=0$, where
$\mathcal{G}_{\rm res}\equiv(1-w_{\|}/v_{\rm ph})\partial/\partial
w_{\perp}+(w_{\perp}/v_{\rm ph})\partial/\partial w_{\|}$ with $v_{\rm
ph}(w_{\|})=\omega/k_{\|}$ for resonant waves that satisfy
$\omega/k_{\|}-\Omega_{i}/k_{\|}=w_{\|}$. With heating being driven by oblique
ICWs but exciting parallel ICWs for $t\lesssim 11\tau_{\rm A}$, we expect
$f_{i}$ to decrease along the oblique-ICW contours (which approximately
satisfy the $\sin\theta\approx 1$ cold-plasma relation $\omega=k_{\|}v_{\rm
A}/\sqrt{1+k_{\|}^{2}d_{i}^{2}}$; Stix, 1992) and to increase along the
parallel-ICW contours ($\omega=k_{\|}v_{\rm A}\sqrt{1-\omega/\Omega_{i}}$; see
Isenberg 2012 for further discussion of cold-plasma relations in this
context). Its evolution towards a flatter core at $t\gtrsim 11\tau_{\rm A}$ is
hard to discern on $f_{i}(w_{\perp},w_{\|})$, so in the lower panel we plot
$f_{i}$ along an oblique-ICW resonant contour (thin black line). It changes
from decreasing slightly to increasing along the contour, with the latter a
signature of $f_{i}(w_{\perp},w_{\|})$ becoming flatter in $w_{\perp}$, likely
due to stochastic heating (a $w_{\perp}$-independent $f_{i}$ is shown with the
dashed line; Klein & Chandran, 2016; Cerri et al., 2021). Various other
modifications to $f_{i}$ over the simulation are clear in the top panel,
including the smoothing of the sharp “ridge” bordering the quasi-linearly
heated region of phase space at $w_{\|}/v_{\rm th0}\approx-0.5$, and the
flattening of $f_{i}(w_{\perp})$ at $w_{\|}>0$.

Another feature of $f_{i}$ is the strong parallel plateau, or beam, which
results from the Landau damping of $k_{\perp}\gg k_{\|}$ fluctuations at
$k_{\perp}\rho_{i}\lesssim 1$ (see also S+22; Li et al. 2010). Such
fluctuations, which lie in the $k_{\perp}$ range between Alfvén waves and
KAWs, propagate with phase speed ${>}v_{\rm A}$ (Howes et al., 2006), growing
a modestly super-Alfvénic plateau in $f_{i}$. As the turbulence becomes
balanced, such fluctuations would likely flatten $f_{i}$ at $w_{\|}\approx-
v_{\rm A}$ also (see figure 7 of Arzamasskiy et al. 2019), but this feature is
not yet observable. The total parallel heating $Q_{\|,i}$, which captures this
beam formation, remains almost constant throughout the full simulation (Fig.
3a inset). This suggests that the $k_{\perp}\rho_{i}\sim 1$ Alfvénic
fluctuations are Landau damping some fixed portion of their energy before
dissipating into oblique ICWs (at earlier times) or KAW turbulence (at later
times). A corollary is that there is no fundamental difference between beam
formation and standard resonant parallel heating, aside from the fluctuations’
imbalance.

##  4 Discussion

A complete theory for the helicity barrier would be able to predict the
parameters ($\varepsilon_{H}/\varepsilon$, $\beta$, etc.) at which it is
operative. Low-$\beta$ gyrokinetics states only that the barrier occurs unless
the generalized helicity can be destroyed at least as fast as
${\sim}\varepsilon_{H}$ before $k_{\perp}\rho_{i}\sim 1$ scales are reached
(Meyrand et al., 2021). Although generalized helicity is not a true invariant
of hybrid kinetics, we can infer from our simulation, based on the agreement
between $Q_{e}$ and $\varepsilon-\varepsilon_{H}$ (Fig. 1) and the continuous
evolution of the transition-range spectrum (Fig. 2) for $t\lesssim 11\tau_{\rm
A}$, that the critical $\varepsilon_{H}/\varepsilon$ at $\beta\approx 0.4$ is
${\approx}0.2$. The changes that occur once
$\varepsilon_{H}/\varepsilon\lesssim 0.2$ also coincide with the sharp drop in
$Q_{\perp,i}$ (Fig. 3), suggesting that breaking the helicity barrier curtails
ion heating. These results highlight the surprising robustness of the helicity
barrier in the face of the additional complexities of a true kinetic system.
However, $\varepsilon_{H}/\varepsilon$ decreases unrealistically fast in our
simulation; understanding how the barrier evolves when
$\varepsilon_{H}/\varepsilon$ changes on timescales comparable to the
turbulent decay time should be a priority for future work.

Our results provide a helpful roadmap for understanding the radially dependent
interplay between turbulence, heating, and instabilities in the $\beta\lesssim
1$ solar wind. At smaller $R$ and in faster streams, which have higher
observed imbalance, we predict strong perpendicular ion heating via quasi-
linear resonance, continual emission (instability) of parallel ICWs, a steep
and wide ion-Larmor-scale transition range, proton beam formation, and little
electron heating. At larger $R$ and in slower streams, which have lower
imbalance ($\sigma_{c}\lesssim 0.8$ in our simulation), electron heating
dominates, parallel ICWs are absorbed/damped, there is no transition-range
spectrum, and sharp features in the VDF are smoothed out. These correlations
and features match those measured in the low-$\beta$ solar wind by _in situ_
spacecraft (e.g., Marsch, 2006; Bruno et al., 2014; Zhao et al., 2021; Shi et
al., 2023).

Under the assumption that, far from the Sun, faster wind is heated more than
slower wind (Hansteen & Leer, 1995; Totten et al., 1995; Halekas et al.,
2023), our results explain qualitatively various interesting properties of
observed temperature profiles. Close to the Sun, proton and minor-ion
temperatures correlate positively with wind speed (Burlaga & Ogilvie, 1973),
while the electron temperature is negatively correlated (Marsch et al., 1989);
this is natural if the highly imbalanced turbulence of faster streams has low
$Q_{e}/Q_{i}$ (Shi et al., 2023). At larger distances (${\sim}1~{}{\rm au}$),
the electron-temperature correlation flips to positive for $U\lesssim
500~{}{\rm km~{}s}^{-1}$ (Shi et al., 2023), as would occur if $Q_{e}/Q_{i}$
increased to ${\gtrsim}1$ as the turbulence becomes balanced at larger radii.
More directly, Abraham et al. (2022) measure $Q_{e}\propto R^{-2}$ for
$R\lesssim 0.3~{}{\rm au}$, followed by a rapid drop in $Q_{e}$ at larger $R$.
This profile, $Q_{e}\propto R^{-2}$, is much flatter than the “standard”
total-heating profile $Q\propto R^{-4}$ (e.g., Totten et al., 1995)—a natural
explanation is that $Q=Q_{i}+Q_{e}$ is steeper than $Q_{e}$, because
$Q_{e}/Q_{i}$ increases as $\varepsilon_{H}/\varepsilon$ decreases with $R$.
Then, once $\varepsilon_{H}/\varepsilon\ll 1$ (around $R\simeq 0.3$ in this
scenario), $Q_{e}/Q_{i}$ saturates and $Q_{e}$ drops rapidly.

Together, these observations suggest that imbalanced turbulence and the
helicity barrier are actively shaping global coronal and solar-wind dynamics.
More generally, they highlight the crucial role of imbalance in controlling
collisionless plasma thermodynamics.

We thank S. S. Cerri, C. H. K. Chen, B. D. G. Chandran, E. Quataert, and A. A.
Schekochihin for useful discussions. JS and RM acknowledge the support of the
Royal Society Te Apārangi, through Marsden-Fund grants MFP-UOO2221 (JS) and
MFP U0020 (RM), as well as through the Rutherford Discovery Fellowship
RDF-U001804 (JS). This research was part of the Frontera computing project at
the Texas Advanced Computing Center, which is made possible by National
Science Foundation award OAC-1818253. Further computational support was
provided by the New Zealand eScience Infrastructure (NeSI) high performance
computing facilities, funded jointly by NeSI’s collaborator institutions and
through the NZ MBIE, and through PICSciE-OIT TIGRESS High Performance
Computing Center and Visualization Laboratory at Princeton University.

## References

  * Abraham et al. (2022) Abraham, J. B., Verscharen, D., Wicks, R. T., et al. 2022, Astrophys. J., 941, 145, doi: 10.3847/1538-4357/ac9fd8
  * Arzamasskiy et al. (2019) Arzamasskiy, L., Kunz, M. W., Chandran, B. D. G., & Quataert, E. 2019, Astrophys. J., 879, 53, doi: 10.3847/1538-4357/ab20cc
  * Arzamasskiy et al. (2023) Arzamasskiy, L., Kunz, M. W., Squire, J., Quataert, E., & Schekochihin, A. A. 2023, Phys. Rev. X, 13, 021014, doi: 10.1103/PhysRevX.13.021014
  * Bale et al. (2019) Bale, S. D., Badman, S. T., Bonnell, J. W., et al. 2019, Nature, 576, 237, doi: 10.1038/s41586-019-1818-7
  * Bandyopadhyay et al. (2023) Bandyopadhyay, R., Meyer, C. M., Matthaeus, W. H., et al. 2023, Astrophys. J. Lett., 955, L28, doi: 10.3847/2041-8213/acf85e
  * Bowen et al. (2022) Bowen, T. A., Chandran, B. D. G., Squire, J., et al. 2022, Phys. Rev. Lett., 129, 165101, doi: 10.1103/PhysRevLett.129.165101
  * Bowen et al. (2023) Bowen, T. A., Bale, S. D., Chandran, B. D. G., et al. 2023, arXiv:2306.04881, doi: 10.48550/arXiv.2306.04881
  * Bruno & Carbone (2013) Bruno, R., & Carbone, V. 2013, Living Rev. Solar Phys., 10, 2, doi: 10.12942/lrsp-2013-2
  * Bruno et al. (2014) Bruno, R., Trenchi, L., & Telloni, D. 2014, Astrophys. J. Lett., 793, L15, doi: 10.1088/2041-8205/793/1/L15
  * Burlaga & Ogilvie (1973) Burlaga, L. F., & Ogilvie, K. W. 1973, J. Geophys. Res., 78, 2028, doi: 10.1029/JA078i013p02028
  * Cerri et al. (2021) Cerri, S. S., Arzamasskiy, L., & Kunz, M. W. 2021, Astrophys. J., 916, 120, doi: 10.3847/1538-4357/abfbde
  * Chandran et al. (2010a) Chandran, B. D. G., Li, B., Rogers, B. N., Quataert, E., & Germaschewski, K. 2010a, Astrophys. J., 720, 503, doi: 10.1088/0004-637X/720/1/503
  * Chandran et al. (2010b) Chandran, B. D. G., Pongkitiwanichakul, P., Isenberg, P. A., et al. 2010b, Astrophys. J., 722, 710, doi: 10.1088/0004-637X/722/1/710
  * Chen et al. (2016) Chen, C. H. K., Matteini, L., Schekochihin, A. A., et al. 2016, Astrophys. J., 825, 1, doi: 10.3847/2041-8205/825/2/l26
  * D’Amicis et al. (2021) D’Amicis, R., Alielden, K., Perrone, D., et al. 2021, Astron. Astro., 654, A111, doi: 10.1051/0004-6361/202140600
  * Dobrowolny et al. (1980) Dobrowolny, M., Mangeney, A., & Veltri, P. 1980, Phys. Rev. Lett., 45, 144, doi: 10.1103/PhysRevLett.45.144
  * Duan et al. (2021) Duan, D., He, J., Bowen, T. A., et al. 2021, Astrophys. J. Lett., 915, L8, doi: 10.3847/2041-8213/ac07ac
  * Fjørtoft (1953) Fjørtoft, R. 1953, Tellus, 5, 225
  * Halekas et al. (2023) Halekas, J. S., Bale, S. D., Berthomier, M., et al. 2023, arXiv:2305.13424, doi: 10.48550/arXiv.2305.13424
  * Hansteen & Leer (1995) Hansteen, V. H., & Leer, E. 1995, J. Geophys. Res., 100, 21577, doi: 10.1029/95JA02300
  * Horbury et al. (2012) Horbury, T. S., Wicks, R. T., & Chen, C. H. K. 2012, Space Sci. Rev., 172, 325, doi: 10.1007/s11214-011-9821-9
  * Howes et al. (2006) Howes, G. G., Cowley, S. C., Dorland, W., et al. 2006, Astrophys. J., 651, 590, doi: 10.1086/506172
  * Isenberg (2012) Isenberg, P. A. 2012, Phys. Plasmas, 19, 032116, doi: 10.1063/1.3697721
  * Isenberg & Vasquez (2011) Isenberg, P. A., & Vasquez, B. J. 2011, Astrophys. J., 731, 88, doi: 10.1088/0004-637X/731/2/88
  * Kawazura et al. (2019) Kawazura, Y., Barnes, M., & Schekochihin, A. A. 2019, Proc. Nat. Acad. Sci., 116, 771, doi: 10.1073/pnas.1812491116
  * Kennel & Engelmann (1966) Kennel, C. F., & Engelmann, F. 1966, Phys. Fluids, 9, 2377, doi: 10.1063/1.1761629
  * Kennel & Wong (1967) Kennel, C. F., & Wong, H. V. 1967, J. Plasma Phys., 1, 75, doi: 10.1017/S002237780000310X
  * Klein & Chandran (2016) Klein, K. G., & Chandran, B. D. G. 2016, Astrophys. J., 820, 47, doi: 10.3847/0004-637X/820/1/47
  * Klein & Howes (2016) Klein, K. G., & Howes, G. G. 2016, Astrophys. J. Lett., 826, L30, doi: 10.3847/2041-8205/826/2/L30
  * Kunz et al. (2014) Kunz, M. W., Stone, J. M., & Bai, X.-N. 2014, J. Comp. Phys., 259, 154, doi: 10.1016/j.jcp.2013.11.035
  * Li et al. (1999) Li, X., Habbal, S. R., Hollweg, J. V., & Esser, R. 1999, J. Geophys. Res., 104, 2521, doi: 10.1029/1998JA900126
  * Li et al. (2010) Li, X., Lu, Q., Chen, Y., Li, B., & Xia, L. 2010, Astrophys. J. Lett., 719, L190, doi: 10.1088/2041-8205/719/2/L190
  * Marsch (2006) Marsch, E. 2006, Living Rev. Solar Phys., 3, 1, doi: 10.12942/lrsp-2006-1
  * Marsch et al. (1989) Marsch, E., Pilipp, W. G., Thieme, K. M., & Rosenbauer, H. 1989, J. Geophys. Res., 94, 6893, doi: 10.1029/JA094iA06p06893
  * Meyrand et al. (2023) Meyrand, R., Squire, J., Mallet, A., & Chandran, B. D. G. 2023, arXiv:2308.10389, doi: 10.48550/arXiv.2308.10389
  * Meyrand et al. (2021) Meyrand, R., Squire, J., Schekochihin, A. A., & Dorland, W. 2021, J. Plasma Phys., 87, 535870301, doi: 10.1017/S0022377821000489
  * Parker (1965) Parker, E. N. 1965, Space Sci. Rev., 4, 666, doi: 10.1007/BF00216273
  * Roberts et al. (1987) Roberts, D. A., Goldstein, M. L., Klein, L. W., & Matthaeus, W. H. 1987, J. Geophys. Res., 92, 12023, doi: 10.1029/JA092iA11p12023
  * Schekochihin (2022) Schekochihin, A. A. 2022, J. Plasma Phys., 88, 155880501.  https://arxiv.org/abs/2010.00699
  * Schekochihin et al. (2019) Schekochihin, A. A., Kawazura, Y., & Barnes, M. A. 2019, J. Plasma Phys., 85, 905850303, doi: 10.1017/S0022377819000345
  * Shi et al. (2023) Shi, C., Velli, M., Lionello, R., et al. 2023, Astrophys. J., 944, 82, doi: 10.3847/1538-4357/acb341
  * Squire et al. (2022) Squire, J., Meyrand, R., Kunz, M. W., et al. 2022, Nature Astron., 6, 715, doi: 10.1038/s41550-022-01624-z
  * Stix (1992) Stix, T. H. 1992, Waves in plasmas (American Institute of Physics, Melville, NY)
  * Totten et al. (1995) Totten, T. L., Freeman, J. W., & Arya, S. 1995, J. Geophys. Res., 100, 13, doi: 10.1029/94JA02420
  * Vasquez et al. (2020) Vasquez, B. J., Isenberg, P. A., & Markovskii, S. A. 2020, Astrophys. J., 893, 71, doi: 10.3847/1538-4357/ab7e2b
  * Zhao et al. (2021) Zhao, G. Q., Lin, Y., Wang, X. Y., et al. 2021, Astrophys. J., 906, 123, doi: 10.3847/1538-4357/abca3b