# Aberrant High-Order Dependencies in Schizophrenia Resting-State Functional
MRI Networks

Qiang Li1∗, Vince D. Calhoun1, Adithya Ram Ballem1, Armin Iraji1,2∗
1Tri-Institutional Center for Translational Research in Neuroimaging and Data
Science (TReNDS)
Georgia State University, Georgia Institute of Technology, Emory University,
Atlanta, GA, USA
2Department of Computer Science, Georgia State University, Atlanta, GA, USA
&Shujian Yu3
3Department of Computer Science
Vrije Universiteit Amsterdam, The Netherlands &Jesus Malo4
4Image Processing Laboratory
University of Valencia, Spain

###### Abstract

The human brain has a complex, intricate functional architecture. While many
studies primarily emphasize pairwise interactions, delving into high-order
associations is crucial for a comprehensive understanding of how functional
brain networks intricately interact beyond simple pairwise connections.
Analyzing high-order statistics allows us to explore the nuanced and complex
relationships across the brain, unraveling the heterogeneity and uncovering
patterns of multilevel overlap on the psychosis continuum. Here, we employed
high-order independent component analysis (ICA) plus multivariate information-
theoretical metrics ($O$-information and $S$-information) to estimate high-
order interaction to examine schizophrenia using resting-state fMRI. The
results show that multiple brain regions networks may be altered in
schizophrenia, such as temporal, subcortical, and higher-cognitive brain
regions, and meanwhile, it also shows that revealed synergy gives more
information than redundancy in diagnosing schizophrenia. All in all, we showed
that high-order dependencies were altered in schizophrenia. Identification of
these aberrant patterns will give us a new window to diagnose schizophrenia.

##  1 Introduction

Schizophrenia is a major psychotic disease that severely affects people’s
quality of life, and it is usually characterized by delusions, hallucinations,
disorganization, and unusual behavior. The pathogenesis behind schizophrenia
is complex, multifaceted, and involves multiple brain regions, and previous
studies have shown that pair-wise functional connectomics is altered as one of
the features of schizophrenia FU2021117385 (1, 2, 3, 4). However, pair-wise
functional connectomics will not capture multiple network statistical
relationships and will ignore some high-order functional interactions that may
play an important role in its use as a biomarker for schizophrenia
xie2021constructing (5, 6, 7). The related studies have shown that
multivariate mutual information can be used as a high-order functional
connectivity descriptor and applied to the diagnosis of psychotic disorders
herzog2022genuine (8, 9, 10, 11, 12, 13, 14). Therefore, in this study, we
hypothesize that aberrant high-order dependencies are present in
schizophrenia, and we investigate this using high model-order independent
component analysis (ICA) plus information-theoretical metrics ($O$-information
and $S$-information) to estimate high-order brain network interaction to
examine schizophrenia using resting-state fMRI (rsfMRI).

##  2 Materials and Methods

###  2.1 rsfMRI Dataset

The total dataset of 1004 subjects from the Bipolar and Schizophrenia Network
for Intermediate Phenotypes (BSNIP) consortium tamminga2013clinical (15) (each
half subject in typical controls and schizophrenia) used in this work was
obtained by using 105 intrinsic connectivity network (ICN) (see Appendices, in
Fig. 3) time courses derived from a multi-spatial-scale, spatially constrained
ICA approach, and this allows generalizability and comparability of findings
across studies. The shape of the time course is [T, 105], where T is the
number of time points that are different between subjects Iraji2022CanonicalAR
(16). The rsfMRI preprocessed steps as shown in Fig. 1: (a.) Quality control
was applied to identify high-quality data. (b.) Each subject’s rsfMRI data
were preprocessed using a common procedure, including rigid body motion
correction, slice timing correction, and distortion correction. (c.)
Preprocessed subject data were registered into a common space, resampled to
$3mm^{3}$ isotropic voxels, and spatially smoothed using a Gaussian kernel
with a $6mm$ full width at half-maximum. (d.) We used a multi-spatial-scale
template of 105 ICNs obtained from 100k+ subjects and a constraint ICA
approach to obtain subject-specific ICN time courses, and they were cleaned
using a common standard. (e.) The estimated 105 intrinsic connectivity
networks are categorized into visual (VI), cerebellar (CB), temporal (TP),
subcortical (SC), somatomotor (SM), and higher cognitive (HC) groups for the
next information theory analysis.

###  2.2 Information-Theoretical Analysis

The _Total Correlation (TC)_ and _Dual Total Correlation (DTC)_ describe the
dependence among $n$ variables and can be considered as a non-negative
generalization of the concept of mutual information from two parties to $n$
parties. Let the definition of total correlation due Watanabe
watanabe1960information (17) be denoted as:

$\begin{split}&\mathbf{TC}\left(X^{1},\cdots,X^{n}\right)=\sum_{i=1}^{n}\mathbf{H}\left(X^{i}\right)-\mathbf{H}\left(X^{1},\cdots,X^{n}\right)\\\
\end{split}$ (1)

Analogously with _TC_ , _DTC_ can be defined as Han78 (18),

$\begin{split}&\mathbf{DTC}\left(X^{1},\cdots,X^{n}\right)=\mathbf{H}\left(X^{1},\cdots,X^{n}\right)-\sum_{i=1}^{n}\mathbf{H}\left(X^{i}\mid
X^{[n]\backslash i}\right)\\\
&=\left[\sum_{i=1}^{n}\mathbf{H}\left(X^{[n]\backslash
i}\right)\right]-(n-1)\mathbf{H}\left(X^{1},X^{2},\cdots,X^{n}\right)\end{split}$
(2)

where $X^{[n]\backslash i}=\\{X^{1},\cdots,X^{i-1},X^{i+1},\cdots,X^{n}\\}$,
i.e., the set of all variables excluding $X^{i}$. From these definitions in
Eqs. 1, 2, if all variables are independent, both _TC_ and _DTC_ will be zero.

From the _TC_ and _DTC_ , we get $O$-Information and $S$-Information, which
capture the overall balance between redundant and synergistic high-order
dependencies in multivariate system rosas2019quantifying (19).

The $O$-Information of a set of $n$ random variables $X^{n}$ can be get as,

$\mathbf{\Omega_{n}}\left(X^{1},\cdots,X^{n}\right)=\mathbf{TC}\left(X^{1},\cdots,X^{n}\right)-\mathbf{DTC}\left(X^{1},\cdots,X^{n}\right)$
(3)

The $S$-Information of a set of $n$ random variables can be calculate as,

$\mathbf{\Sigma_{n}}\left(X^{1},\cdots,X^{n}\right)=\mathbf{TC}\left(X^{1},\cdots,X^{n}\right)+\mathbf{DTC}\left(X^{1},\cdots,X^{n}\right)$
(4)

if $\emph{TC}>\emph{DTC}$ or $\mathbf{\Omega}>0$, the system corresponds to
redundancy-dominated dependencies, and in the other way around,
$\emph{DTC}>\emph{TC}$ or $\mathbf{\Omega}<0$, the system corresponds to
synergy-dominated dependencies rosas2019quantifying (19).

Consider the fMRI signal as a Gaussian distribution, as shown in Fig. 1C, and
for a univariate Gaussian random variable $X\sim\mathcal{N}(\mu,\sigma)$, the
entropy (given in nats) will be
$\mathbf{H}^{\mathcal{N}}(X)=\frac{\ln\left(2\pi e\sigma^{2}\right)}{2}$, and
for a multivariate Gaussian distribution, the joint entropy will be
$\mathbf{H}^{\mathcal{N}}(\left(X^{1},\cdots,X^{n}\right)=\frac{\ln\left[(2\pi
e)^{n}|\Sigma|\right]}{2}$, where $|\Sigma|$ refers to the determinant of the
covariance matrix of $\left(X^{1},\cdots,X^{n}\right)$. The for the
multivariate Gaussian case, the mutual information between
$\left(X^{1},\cdots,X^{n}\right)$ and $\left(Y^{1},\cdots,Y^{n}\right)$ is
given by, $\mathbf{I}((X^{1},\cdots,X^{n});(Y^{1},\cdots,Y^{n}))=\frac{1}{2\ln
2}\ln\left[\frac{\left|\Sigma_{(}X^{1},\cdots,X^{n})\right|\left|\Sigma_{(}Y^{1},\cdots,Y^{n})\right|}{\left|\Sigma_{(X^{1},\cdots,X^{n})(Y^{1},\cdots,Y^{n})}\right|}\right]$.
Finally, the Gaussian estimator for total correlation can be:

Figure 1: Flowchart of capturing high-order dependencies among brain
connectivity states. A. Group ICA is performed on two independent healthy
control and schizophrenia datasets, and corresponding 105 spatial network maps
and time courses are extracted from each subject. B. The estimated independent
components are matched with the NeuroMark_2.1 template Iraji2022CanonicalAR
(16) based on spatial correlation, and then all components are categorized
into six groups, i.e., VI, CB, TP, SC, SM, and HC. C. The BOLD signal from
independent components matches a multivariate Gaussian distribution and
related high-order (i.e., interaction order: 3, 4, 5) information measures
with $O$-information and $S$-information.

$\mathbf{TC}^{\mathcal{N}}\left(X^{1},\cdots,X^{n}\right)=\frac{-\ln(|\Sigma|)}{2}$
(5)

From Varley2022MultivariateIT (20), the Eq. 3, 4 will be changed to,

$\mathbf{\Omega_{n}}^{\mathcal{N}}\left(X^{1},\cdots,X^{n}\right)=(2-n)\mathbf{TC}^{\mathcal{N}}\left(X^{1},\cdots,X^{n}\right)+\sum_{i=1}^{n}\mathbf{TC}^{\mathcal{N}}\left(X^{[n]\backslash
i}\right)$ (6)

From Eq. 3, 4, 5, 6, it can calculate all of the metrics described above,
i.e., _TC_ , _DTC_ , _O_ -information and _S_ -information for multivariate
Gaussian variables.

##  3 Results and Discussion

The low-order functional connectivity (LOFC) based on the pair-wise Pearson
correlation in healthy control (HC) subjects and schizophrenia (SC) patients
was present in Fig. 2 and Fig. 4 (see Appendices). We see that major aberrant
brain network interaction happened between visual-somatomotor and temporal-
high cognitive brain areas, and it reflects that three brain networks, such as
visual, somatomotor, and higher cognitive, may be altered in the SC patients.

Figure 2: Comparison of low-order functional connectivity (LOFC) and high-
order functional connectivity (HOFC) between HC subjects and SC patients. The
top row shows that LOFC and the main difference present a pair-wise
interaction between VI-SM and TP-HC. The bottom row presents the HOFC with
redundancy and synergy information-theoretic measures made from HC subjects
and SC patients, and the HOFC not only captures redundancy information but
also synergy information in the human brain, and with interaction order
increasing, we see that more difference is captured by the HOFC.

From the high-order functional connectivity (HOFC), the brain information
networks integration and segregation both present from multivariate
information theory measures (in Fig. 2), and we see that redundancy
information basically mainly distributed at the SC and HC brain regions with
interaction order gradually increases, and it also captures more aberrant
brain networks, such as TP, SC, and HC brain regions, and it may suggest that
these brain regions are involved in the SC patients. Moreover, synergy
information was also uncovered by HOFC, and it showed that synergy information
plays a more important role than redundancy in SC patients, and it may suggest
that synergy is mainly aberrant in SC patients, as it can be a biomarker for
diagnosing SC patients.

However, there are still some limitations, and the next step would be to dig
into specific independent component networks and see which independent
component networks mainly contribute to the SC patients. Second, explaining
the role of redundancy and synergy in the SC patients and how to properly
explain them from a neurophysics perspective.

## References

  * (1) Zening Fu et al.  “Dynamic state with covarying brain activity-connectivity: On the pathophysiology of schizophrenia”  In _NeuroImage_ 224, 2021, pp. 117385
  * (2) Xing Meng et al.  “Multi-model order spatially constrained ICA reveals highly replicable group differences and consistent predictive results from resting data: A large N fMRI schizophrenia study”  In _NeuroImage: Clinical_ 38, 2023, pp. 103434
  * (3) Yuhui Du et al.  “Evidence of shared and distinct functional and structural brain signatures in schizophrenia and autism spectrum disorder”  In _Communications Biology_ 4, 2021
  * (4) Armin Iraji et al.  “Spatial Dynamic Subspaces Encode Sex-Specific Schizophrenia Disruptions in Transient Network Overlap and its Links to Genetic Risk”  In _bioRxiv : the preprint server for biology_ , 2023
  * (5) Qingsong Xie et al.  “Constructing high-order functional connectivity network based on central moment features for diagnosis of autism spectrum disorder”  In _PeerJ_ 9 PeerJ Inc., 2021, pp. e11692
  * (6) Yu Zhang et al.  “Hybrid high-order functional connectivity networks using resting-state functional MRI for mild cognitive impairment diagnosis”  In _Scientific reports_ 7.1 Nature Publishing Group UK London, 2017, pp. 6530
  * (7) Xiaobo Chen et al.  “High-order resting-state functional connectivity network for MCI classification”  In _Human brain mapping_ 37.9 Wiley Online Library, 2016, pp. 3282–3296
  * (8) Rubén Herzog et al.  “Genuine high-order interactions in brain networks and neurodegeneration”  In _Neurobiology of Disease_ 175 Elsevier, 2022, pp. 105918
  * (9) Qiang Li  “Measuring Functional Connectivity of Human Intra-Cortex Regions with Total Correlation”  In _Entropy 2021: The Scientific Tool of the 21st Century_ , 2021, pp. 9797
  * (10) Qiang Li  “Functional connectivity inference from fMRI data using multivariate information measures”  In _Neural Networks_ 146 Elsevier, 2022, pp. 85–97
  * (11) Qiang Li, Greg Ver Steeg, Shujian Yu and Jesus Malo  “Functional Connectome of the Human Brain with Total Correlation”  In _Entropy_ 24.12 Multidisciplinary Digital Publishing Institute, 2022, pp. 1725
  * (12) Qiang Li, Greg Ver Steeg and Jesus Malo  “Functional connectivity in visual areas from Total Correlation”  In _arXiv preprint arXiv:2208.05770_ , 2022
  * (13) Qiang Li et al.  “Higher-Order Organization in the Human Brain From Matrix-Based Rényi’s Entropy”  In _2023 IEEE International Conference on Acoustics, Speech, and Signal Processing Workshops (ICASSPW)_ , 2023, pp. 1–5
  * (14) Marilyn Gatica et al.  “High-order interdependencies in the aging brain”  In _Brain connectivity_ 11.9 Mary Ann Liebert, Inc., publishers 140 Huguenot Street, 3rd Floor New …, 2021, pp. 734–744
  * (15) Carol A Tamminga et al.  “Clinical phenotypes of psychosis in the Bipolar-Schizophrenia Network on Intermediate Phenotypes (B-SNIP)”  In _American Journal of psychiatry_ 170.11 Am Psychiatric Assoc, 2013, pp. 1263–1274
  * (16) Armin Iraji et al.  “Canonical and Replicable Multi-Scale Intrinsic Connectivity Networks in 100k+ Resting-State fMRI Datasets”  In _bioRxiv_ , 2022  DOI: 10.1101/2022.09.03.506487
  * (17) Satosi Watanabe  “Information theoretical analysis of multivariate correlation”  In _IBM Journal of research and development_ 4.1 IBM, 1960, pp. 66–82
  * (18) Te Sun Han  “Nonnegative Entropy Measures of Multivariate Symmetric Correlations”  In _Inf. Control._ 36.2, 1978, pp. 133–156
  * (19) Fernando E Rosas, Pedro AM Mediano, Michael Gastpar and Henrik J Jensen  “Quantifying high-order interdependencies via multivariate extensions of the mutual information”  In _Physical Review E_ 100.3 APS, 2019, pp. 032305
  * (20) Thomas F. Varley, Maria Pope, Joshua Faskowitz and Olaf Sporns  “Multivariate information theory uncovers synergistic subsystems of the human cerebral cortex”  In _Communications Biology_ 6, 2022
  * (21) Ronald Coifman and Stéphane Lafon  “Diffusion maps”  In _Applied and Computational Harmonic Analysis_ 21, 2006, pp. 5–30

##  4 Appendices

###  4.1 Independent Component Networks

Total of 105 independent component networks identified from multivariate
spatially constrained ICA are shown in Fig. 3.

Figure 3: A total of 105 independent component networks are presented on the
overlay brain volume, and different colors refer to each independent component
network. Each component is also displayed separately.

###  4.2 Low-Order Functional Connectogram

The low-order functional connectivity estimated for 105 independent component
networks, and we applied the diffusion embedding algorithm Coifman06acha (21)
to cluster the related independent component networks.

Figure 4: Comparison of clustering low-order functional connectivity (LOFC)
between healthy controls (HC subjects, left side) and schizophrenia (SC
patients, right side) The top row shows that LOFC presents a pair-wise
interaction. The bottom row presents the functional connectogram corresponding
to each correlation map with a threshold set to 0.6.

Here we clearly see that interaction-independent component network patterns
are clumped, and the pattern can be matched to six multi model order ICA brain
template networks (see Fig. 4). The corresponding functional connectogram was
also presented with a threshold of 0.6.

###  4.3 High-Order Functional Connectogram

Given the number of independent component networks, exploring the high-order
functional connectogram with 105 independent component networks along with
increasing interaction order will face some challenges. Furthermore, it will
also face how to present the result with a graph connectogram, and all these
problems will be major problems that need to be solved, but it will be very
important for us to precision diagnose schizophrenia, and the previous related
research li2022functionalen (11) has already provided us with a window to
solve it, and we can apply it to the precision diagnosis of schizophrenia.


# Exact models of pure radiation in $R^{2}$ gravity
for spatially homogeneous
Shapovalov spacetimes type II

Konstantin Osetrin,
Altair Filippov,
Evgeny Osetrin,
Tomsk State Pedagogical University,
Tomsk, 634061, Russia
osetrin@tspu.edu.rualtair<EMAIL_ADDRESS>
###### Abstract

Presented are exactly integrable models with pure radiation in $R^{2}$ gravity
with a cosmological constant, related to wave-like Shapovalov spacetimes type
II. Spatially homogeneous models of Shapovalov spaces were considered. The
obtained solutions belong to spaces of type III according to the Bianchi
classification. For the models under consideration, exact solutions for the
equations of motion of test particles are obtained in the Hamilton-Jacobi
formalism.

##  1 Introduction

The construction of models for the early stages of the of the Universe
expansion is often based on the concept of the radiation dominant era when
deviations from the homogeneity and isotropy of space-time are possible and it
is necessary to take into account quantum corrections to the field equations.
One of the most famous models, in this case is, $R^{2}$ gravity (see f.e.
[1]-[4]). Note that due to the complexity of the field equations, there are
few exact solutions for this theory, usually these are solutions inherited
from Einstein’s theory.

In this work, we obtained exact models in $R^{2}$ gravity, which have a number
of interesting properties.

First, the models considered in this work refer to Shapovalov wave-like spaces
(see [5],[6],[7]), i.e. allow separating the variables in the eikonal equation
and the Hamilton-Jacobi equation for test particles, and one of the separating
variables is the wave one, which indicates the wave nature of the solutions.

Secondly, the considered spaces are homogeneous, but non-isotropic models of
space-time, i.e., are more complex models than the standard Friedman-
Robertson-Walker solution and are adequate as models for the early stages of
the expansion of the Universe.

In the third, the considered models use pure radiation as the energy-momentum
tensor of matter, which can simulate the high-energy component of various
types of radiation (gravitational, electromagnetic, dark energy) and
corresponds to the concept of the radiation dominant epoch of the Universe
expansion at the stage of its isotropization.

As it is known ([1]-[4]), the field equations of $R^{2}$ gravity are obtained
from the action of the form

$S=\frac{1}{2\kappa^{2}}\int d^{4}x\sqrt{-g}\,(2\Lambda+R+\gamma
R^{2})+S_{M},$ (1)

where $\Lambda$ is the cosmology constant, $R$ is the scalar curvature,
$\gamma$ is the constant and $S_{M}$ is the action of the matter fields.

The field equations of the theory under consideration can be written in the
following form:

$(1+2\gamma R)\,R_{ij}-\frac{1}{2}\left(2\Lambda+R+\gamma
R^{2}\right)g_{ij}-2\gamma\left(\nabla_{i}\nabla_{j}-g_{ij}\nabla^{k}\nabla_{k}\right)R=\kappa^{2}T_{ij},$
(2)

where $\nabla_{i}$ is the covariant derivative, $T_{ij}$ is the energy-
momentum tensor of the matter fields.

The energy-momentum tensor of pure radiation has the form

$T_{ij}=\varepsilon L_{i}L_{j}\qquad L^{k}L_{k}=0,$ (3)

where $\varepsilon$ is the radiation energy density, $L^{i}$ is the wave
vector.

The paper deals with Shapovalov spaces of type II. According to the definition
of Shapovalov spaces (see [5, 6, 7]), these spaces admit the privileged
coordinate systems, where the eikonal equation and the Hamilton-Jacobi
equation of test particles can be integrated by the method of complete
separation of variables, and among the non-ignored separated variables, there
are wave variables. Therefore, Shapovalov spaces are called wave-like. The
type of Shapovalov space indicates the number of ignored variables on which
the metric of the space does not depend or, which is the same, the number of
mutually commuting Killing vectors admitted by these spaces.

In this paper, we obtain exact solutions to the field equations of $R^{2}$
theory of gravity (2) for classes of wave-like spatially homogeneous space-
time models that allow separation of isotropic (wave) variables in the eikonal
equation

$g^{ij}\,\partial_{i}\Psi\,\partial_{j}\Psi=0$ (4)

and in the Hamilton-Jacobi equation for a test particle of mass $m$

$g^{ij}\,\partial_{i}S\,\partial_{j}S=m^{2},$ (5)

where $\Psi$ is the eikonal function, $S$ is the test particle action
function, $\partial_{i}$ is the partial derivative.

Spaces allowing the integration of the eikonal equation and the Hamilton-
Jacobi equation of the test particles by the method of complete separation of
variables have found wide application in the theory of gravity and cosmology,
including the study of models with an electromagnetic field ([8], [9],
[10],[11],[12]), dust matter [13], pure radiation ([14],[15],[16]), scalar
fields ([17], [18], [19]), spinor fields ([20], [21]).

In works ([22], [23], [24], [25] ) we have carried out a classification of
these types of spaces from the point of view of the symmetries they admit,
ensuring their spatial homogeneity, i.e. selection of types related to
cosmological problems.

For Shapovalov spaces of type II (two non-ignored variables in a privileged
coordinate system), there are two spatially homogeneous wave-like spacetime
models related to the type III according to Bianchi’s classification [25]. For
the case of two non-ignored variables, there are two types of wave-like
Shapovalov spaces - types B1 and B2 (see [23], [25]), which will be discussed
below.

The methods considered in this paper may be of interest when constructing
cosmological models for f(R) theories of gravity and in models with Gauss –
Bonnet terms (see f.e. [26], [27] ).

##  2 Spatially homogeneous wave-like model type B1

The interval for the space-time model of type B1 has the form (see [25]):

$ds^{2}=\frac{1}{{x^{3}}^{2}}\left(2\,dx^{0}dx^{1}+({x^{0}}-{\alpha})^{1-{\beta}}({x^{0}}+{\alpha})^{1+{\beta}}\,{dx^{2}}^{2}+{dx^{3}}^{2}\right),$
(6)

where $x^{0}$ is an isotropic (wave) variable, $\alpha$ and $\beta$ are
constants ($\alpha\neq 0$, $\beta\neq\pm 1$).

${g=\det
g_{ij}=}-\frac{({x^{0}}-{\alpha})^{1-{\beta}}({\alpha}+{x^{0}})^{{\beta}+1}}{{x^{3}}^{8}},\qquad
x^{0}>|\alpha|.$

Independent Killing vector fields of model type B1 in a privileged coordinate
system can be selected in the form:

$X_{0}=\partial_{1},\qquad X_{1}=\partial_{2},\qquad
X_{2}=2\,x^{1}\partial_{1}+x^{2}\partial_{2}+x^{3}\partial_{3},$
$X_{3}=({x^{0}}^{2}-\alpha^{2})\partial_{0}-\frac{{x^{3}}^{2}}{2}\partial_{1}+{\alpha}\beta
x^{2}\partial_{2}+x^{0}x^{3}\partial_{3}.$ (7)

Killing vectors $X_{1}$, $X_{2}$, $X_{3}$ define a subgroup of spatial
homogeneity of the model. Killing vector commutators of model type B1 have the
form:

$[{X_{0}},{X_{1}}]=0,\qquad[{X_{0}},{X_{2}}]=2{X_{0}},\qquad[{X_{0}},{X_{3}}]=0,$
$[{X_{1}},{X_{2}}]={X_{1}},\qquad[{X_{1}},{X_{3}}]={\alpha}\beta{}{X_{1}},\qquad[{X_{2}},{X_{3}}]=0.$
(8)

For $\alpha\beta=0$, this space admits a third commuting Killing vector and
degenerates into a space with one non-ignored variable only in the privileged
coordinate system.

The Riemann tensor $R_{ijkl}$, the Ricci tensor $R_{ij}$ and the scalar
curvature $R$ have the following nonzero components:

$R_{0313}=-\frac{1}{{x^{3}}^{4}},\quad
R_{0212}=R_{2323}=-\frac{(x^{0}-\alpha)^{1-\beta}(x^{0}+\alpha)^{1+\beta}}{{x^{3}}^{4}},$
(9)
$R_{0202}=\frac{\alpha^{2}(1-\beta^{2})(x^{0}-\alpha)^{-1-\beta}(x^{0}+\alpha)^{-1+\beta}}{{x^{3}}^{2}},$
(10)
$R_{00}=\frac{\alpha^{2}(1-\beta^{2})}{({x^{0}}^{2}-\alpha^{2})^{2}},\qquad
R_{01}=R_{33}=-\frac{3}{{x^{3}}^{2}},$ (11)
$R_{22}=-\frac{3\,(x^{0}-\alpha)^{1-\beta}(x^{0}+\alpha)^{1+\beta}}{{x^{3}}^{2}},\qquad
R=-12.$ (12)

Substituting the metric (6) into the field equations (2) with the energy-
momentum tensor of pure radiation (3), we obtain the following results:

$\Lambda=3,\qquad L_{k}=\left(L_{0},0,0,0\right),\qquad L^{k}=(0,L^{1},0,0),$
(13)

and for the radiation we obtain the condition

${\varepsilon}\kappa^{2}{L^{1}}^{2}={\alpha}^{2}\left(1-{\beta}^{2}\right)\left(1-24{\gamma}\right)\left(\frac{{x^{3}}^{2}}{{x^{0}}^{2}-{\alpha}^{2}}\right)^{2}.$
(14)

If we assume the radiation energy density ${\varepsilon}$ and the factor
$(1-24{\gamma})$ to be positive, then we obtain restrictions on the values of
the constant $\beta$:

$-1<\beta<1.$ (15)

Thus, the metric (6) and the conditions (13) – (15) give us an exact solution
to the field equations of $R^{2}$ gravity with the energy-momentum tensor of
pure radiation.

If $\alpha\neq 0$ and $\beta\neq\pm 1$ the resulting solution cannot be
conformally flat, since the two components of the Weyl tensor $C_{ijkl}$ are
not equal to zero:

$C_{0202}=\frac{{\alpha}^{2}\left(1-{\beta}^{2}\right)({x^{0}}-{\alpha})^{-{\beta}-1}({x^{0}}+{\alpha})^{{\beta}-1}}{2{x^{3}}^{2}},$
(16)
$C_{0303}=\frac{-{\alpha}^{2}\left(1-{\beta}^{2}\right)}{2{x^{3}}^{2}({x^{0}}-\alpha)^{2}({x^{0}}+{\alpha})^{2}}.$
(17)

If ${\alpha}=0$ or $\beta=0,\pm 1$, the metric of the model type B1
degenerates - in a privileged coordinate system it depends on one variable
only.

The model B1 is of type III according to the Bianchi classification and has
the type N according to the Petrov classification.

###  2.1 Solution of the Hamilton-Jacobi equation of test particle for the
model type B1

As mentioned earlier, the model type B1 admits in the used privileged
coordinate system complete separation of variables in the eikonal equation and
the Hamilton-Jacobi equation for test particles. The separation constants
$\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are integrals of motion and are
determined by the initial data.

From the Hamilton-Jacobi equation (5) by the method of complete separation of
variables we have the function $S$ of the action of the test particle for
model type B1 in the form ($\alpha\beta\lambda_{1}\neq 0$):

$S=\frac{\lambda_{3}}{2{\lambda_{1}}}\,x^{0}-\frac{{\lambda_{2}}^{2}}{2\alpha\beta{\lambda_{1}}}({x^{0}}-{\alpha})^{\beta}({x^{0}}+{\alpha})^{-\beta}+\lambda_{1}x^{1}+\lambda_{2}x^{2}+m\log
x^{3}$
$+\sqrt{m^{2}-\lambda_{3}{x^{3}}^{2}}-m\log\left(m+\sqrt{m^{2}-\lambda_{3}{x^{3}}^{2}}\right)+F(\lambda_{1},\lambda_{2},\lambda_{3}),$
(18)

where $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are the independent
constants of the motion of test particles and
$F(\lambda_{1},\lambda_{2},\lambda_{3})$ is an arbitrary function of
parameters.

Using the proper time of the test particle $\tau=S/m$ (we put the mass of the
test particle m = 1)

$\tau=\frac{{\lambda_{3}}{x^{0}}}{2{\lambda_{1}}}-\frac{{\lambda_{2}}^{2}({x^{0}}-\alpha)^{\beta}(\alpha+{x^{0}})^{-\beta}}{2\alpha\beta{\lambda_{1}}}+{\lambda_{1}}{x^{1}}+{\lambda_{2}}{x^{2}}$
$\mbox{}+\sqrt{1-{\lambda_{3}}{x^{3}}^{2}}+\log({x^{3}})-\log\left(\sqrt{1-{\lambda_{3}}{x^{3}}^{2}}+1\right),$
(19)

one can obtain a parametric dependence of the coordinates of a particle on its
proper time.

According to the Hamilton-Jacobi method, the trajectory of a particle is
determined by relations of the form:

$\displaystyle\frac{\partial
S}{\partial\lambda_{1}}=\frac{{\lambda_{2}}^{2}({x^{0}}-\alpha)^{\beta}(\alpha+{x^{0}})^{-\beta}}{2\alpha\beta{\lambda_{1}}^{2}}-\frac{{\lambda_{3}}{x^{0}}}{2{\lambda_{1}}^{2}}+{x^{1}}+{c_{1}}=0,$
(20) $\displaystyle\frac{\partial
S}{\partial\lambda_{2}}=-\frac{{\lambda_{2}}({x^{0}}-\alpha)^{\beta}(\alpha+{x^{0}})^{-\beta}}{\alpha\beta{\lambda_{1}}}+{x^{2}}+{c_{2}}=0,$
(21) $\displaystyle\frac{\partial
S}{\partial\lambda_{3}}=\frac{{x^{0}}}{2{\lambda_{1}}}-\frac{{x^{3}}^{2}}{2\left(m+\sqrt{m^{2}-{\lambda_{3}}{x^{3}}^{2}}\right)}+{c_{3}}=0.$
(22)

where $c_{k}=\partial F/\partial\lambda_{k}$ are additional constants.

If $\lambda_{3}<0$, then the motion of the test particle is infinite. For
$\lambda_{3}>0$ the motion of the test particle along the coordinate $x^{0}$
and $x^{3}$ is finite and there are turning points.

The turning points of the motion of test particles along the coordinate
$x^{0}$ for $\lambda_{3}>0$ are determined by solutions of the equation

${{\lambda_{2}}}{}^{2}\left[({x^{0}}-{\alpha})/({x^{0}}+{\alpha})\right]^{\beta}=\lambda_{3}\left({x^{0}}^{2}-\alpha^{2}\right).$

For example, for $\lambda_{3}>0$ and $\beta=0$ we get

$|\alpha|<x^{0}\leq\sqrt{\alpha^{2}+{\lambda_{2}}^{2}/\lambda_{3}}.$

There are also for $\lambda_{3}>0$ the turning points along the coordinate
$x^{3}$:

$0<\left|x^{3}\right|\leq m/\sqrt{\lambda_{3}}.$

Points $x^{0}=|\alpha|$ and $x^{3}=0$ are special because the model has a
singularity at these points.

Let us choose the variable $x^{0}$ as a parameter $x^{0}=t$, then from (20) we
have $x^{1}=x^{1}(t)$, from (21) we get $x^{2}=x^{2}(t)$, from (22)
$x^{3}=x^{3}(t)$ and from (19)
$\tau=\tau\bigl{(}t,x^{1}(t),x^{2}(t),x^{3}(t)\bigr{)}=\tau(t)$. Thus, we have
determined the parametric dependence of the coordinates of the particle as
functions of its proper time $\tau$.

The synchronous coordinate system associated with the test particle, where the
test particle is at rest, is determined by transformation from $x^{k}$ to
$\tilde{x}^{k}=\left(\tau,\lambda_{1},\lambda_{2},\lambda_{3}\right)$, given
by equations (19), (20), (21), (22).

The eikonal function $\Psi$ from Eq.(4) for the model type B1 has the form:

$\Psi=\frac{{\lambda_{3}}^{2}}{2{\lambda_{1}}}\,x^{0}-\frac{{\lambda_{2}}^{2}}{2\alpha\beta{\lambda_{1}}}\left(\frac{{x^{0}}-{\alpha}}{{x^{0}}+{\alpha}}\right)^{\beta}+\lambda_{1}x^{1}+\lambda_{2}x^{2}+\lambda_{3}x^{3}+F(\lambda_{1},\lambda_{2},\lambda_{3}),$
(23)

where $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are the independent
constants and $F(\lambda_{1},\lambda_{2},\lambda_{3})$ is an arbitrary
function of parameters.

##  3 Spatially homogeneous wave-like model type B2

The space-time interval for a type B2 model can be written as (see [25]):

$ds^{2}=\frac{1}{{x^{3}}^{2}}\,\left(2\,dx^{0}dx^{1}+{x^{0}}^{\alpha}\,{dx^{2}}^{2}+{dx^{3}}^{2}\right),$
(24)

where $x^{0}$ is an isotropic (wave) variable, and ${\alpha}$ is a constant.

${g=\det g_{ij}=}-{x^{0}}^{\alpha}/{x^{3}}^{8},\qquad x^{0}>0.$

Independent Killing vector fields in a privileged coordinate system can be
selected in the form:

$X_{0}=\partial_{1},\qquad X_{1}=\partial_{2},\qquad
X_{2}=2\,x^{1}\partial_{1}+x^{2}\partial_{2}+x^{3}\partial_{3},$
$X_{3}=x^{0}\partial_{0}+\frac{1-\alpha}{2}\,x^{2}\partial_{2}+\frac{x^{3}}{2}\,\partial_{3}.$
(25)

Killing vectors $X_{1}$, $X_{2}$, $X_{3}$ define a subgroup of spatial
homogeneity of the model. Killing vector commutators of model type B2 have the
form:

$[{X_{0}},{X_{1}}]=0,\qquad[{X_{0}},{X_{2}}]=2{X_{0}},\qquad[{X_{0}},{X_{3}}]=0,$
$[{X_{1}},{X_{2}}]={X_{1}},\qquad[{X_{1}},{X_{3}}]=\frac{1-\alpha}{2}\,{X_{1}},\qquad[{X_{2}},{X_{3}}]=0.$
(26)

For ${\alpha}=1$, this space admits an additional commuting Killing vector and
degenerates into a space with one non-ignored variable.

The Riemann tensor $R_{ijkl}$, the Ricci tensor $R_{ij}$ and the scalar
curvature $R$ have the following nonzero components:

$R_{0101}=-R_{0313}=\frac{1}{{x^{3}}^{4}},\ \ \
R_{0212}=R_{2323}=-\frac{{x^{0}}^{\alpha}}{{x^{3}}^{4}},\ \ \
R_{0202}=\frac{\alpha(2-\alpha)}{4{x^{3}}^{2}{x^{0}}^{(2-\alpha)}},$
$R_{01}=R_{33}=-\frac{3}{{x^{3}}^{2}},\ \ \
R_{22}=-\frac{3{x^{0}}^{\alpha}}{{x^{3}}^{2}},\ \ \
R_{00}=\frac{\alpha(2-\alpha)}{4{x^{0}}^{2}},\ \ \ R=-12.$ (27)

Substituting the metric (24) into the field equations (2) with the energy-
momentum tensor of pure radiation (3), we obtain the following results:

$\Lambda=3,\qquad L_{k}=\left(L_{0},0,0,0\right),\quad
L^{k}=\left(0,L^{1},0,0,0\right),\quad L^{1}={x^{3}}^{2}\,L_{0},$ (28)
$\kappa{\varepsilon}{L^{1}}^{2}=\frac{{\alpha}(2-{\alpha})(1-24{\gamma})\,{x^{3}}^{4}}{4\,{x^{0}}^{2}}.$
(29)

If we assume the radiation energy density ${\varepsilon}$ and the factor
$(1-24{\gamma})$ to be positive, then we obtain restrictions on the values of
the constant $\alpha$:

$0\leq\alpha\leq 2,$ (30)

Thus, the metric (24) and the conditions (28) – (30) give us an exact solution
of the field equations of $R^{2}$ gravity with the energy-momentum tensor of
pure radiation.

If $\alpha\neq 0$ or $\alpha\neq 2$ the resulting solution cannot be
conformally flat, since the two components of the Weyl tensor $C_{ijkl}$ are
not equal to zero:

${\rm
C}_{0202}=\frac{{\alpha}\,(2-{\alpha}){x^{0}}^{({\alpha}-2)}}{8{x^{3}}^{2}},\qquad{\rm
C}_{0303}=-\frac{{\alpha}\,(2-{\alpha})}{8{x^{0}}^{2}{x^{3}}^{2}}$ (31)

If $\alpha=0$ or $\alpha=2$, the Weyl tensor vanishes (conformally flat
space), but the Ricci tensor, scalar curvature and the Riemann curvature
tensor does not vanish. If ${\alpha}=0,1,2$, the metric of the model type B2
degenerates - in a privileged coordinate system it depends on one variable
only.

This spatially homogeneous space-time model is of type III according to the
Bianchi classification and has type N according to the Petrov classification.

###  3.1 Solution of the Hamilton-Jacobi equation of test particle for the
model type B2

As mentioned earlier, the model type B2 admits in the privileged coordinate
system complete separation of variables in the eikonal equation and the
Hamilton-Jacobi equation for test particles. The separation constants
$\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are integrals of motion and
are determined by the initial data.

From the Hamilton-Jacobi equation (5) by the method of complete separation of
variables we have the function $S$ of the action of the test particle for the
model type B2 in the form ($\lambda_{1}\neq 0$, $\alpha\neq 1$):

$S=\frac{x^{0}}{2{\lambda_{1}}}\left({\lambda_{3}}+{\lambda_{2}}^{2}{x^{0}}^{-\alpha}/(\alpha-1)\right)+{\lambda_{1}}x^{1}+{\lambda_{2}}x^{2}+m\log
x^{3}$
$\mbox{}+\sqrt{m^{2}-{\lambda_{3}}{x^{3}}^{2}}-m\log{\left(m+\sqrt{m^{2}-{\lambda_{3}}{x^{3}}^{2}}\right)}+F(\lambda_{1},\lambda_{2},\lambda_{3}),$
(32)

where $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are the independent
constants of the motion of test particles and
$F(\lambda_{1},\lambda_{2},\lambda_{3})$ is an arbitrary function of
parameters.

According to the Hamilton-Jacobi method, the trajectory of a test particle is
determined by relations of the form:

$\displaystyle\frac{\partial
S}{\partial\lambda_{1}}=x^{1}-\frac{x^{0}}{2{\lambda_{1}}^{2}}\left({\lambda_{3}}+{\lambda_{2}}^{2}{x^{0}}^{-\alpha}/(\alpha-1)\right)+c_{1}=0,$
(33) $\displaystyle\frac{\partial
S}{\partial\lambda_{2}}=x^{2}+\frac{\lambda_{2}}{{\lambda_{1}}(\alpha-1)}{x^{0}}^{(1-\alpha)}+c_{2}=0,$
(34) $\displaystyle\frac{\partial
S}{\partial{\lambda_{3}}}=\frac{x^{0}}{2{\lambda_{1}}}-\frac{{x^{3}}^{2}}{2\left(m+\sqrt{m^{2}-\lambda_{3}{x^{3}}^{2}}\right)}+c_{3}=0,$
(35)

where $c_{k}=\partial F/\partial\lambda_{k}$ are constants.

The proper time of a test particle $\tau=S/m$ can be written as ($m\to 1$):

$\tau=\frac{1}{2{\lambda_{1}}}\left({\lambda_{3}}x^{0}-{\lambda_{2}}^{2}\log
x^{0}\right)+{\lambda_{1}}x^{1}+{\lambda_{2}}x^{2}$
$\mbox{}+\sqrt{1-{\lambda_{3}}{x^{3}}^{2}}+\log
x^{3}-\log{\left(1+\sqrt{1-{\lambda_{3}}{x^{3}}^{2}}\right)},$ (36)

Then, using the variable $x^{0}=t$ as a parameter, from equations (33)-(36) we
obtain the parametric dependence of the coordinates of the test particle
$x^{k}(t)$ on its proper time
$\tau(t)=\tau\left(t,x^{1}(t),x^{2}(t),x^{3}(t)\right)$.

If $\lambda_{3}<0$, then the motion of the test particle for model type B2 is
infinite. For $\lambda_{3}>0$ the motion of the test particle along the
coordinate $x^{0}$ and $x^{3}$ is finite and there are turning points:

$0<x^{0}\leq\left({\lambda_{2}}^{2}/\lambda_{3}\right)^{1/\alpha},\qquad
0<\left|x^{3}\right|\leq m/\sqrt{\lambda_{3}}.$ (37)

The eikonal function $\Psi$ from Eq.(4) for the model type B2 has the
following form ($\lambda_{1}\neq 0$, $\alpha\neq 1$):

$\Psi=\frac{x^{0}}{2{\lambda_{1}}}\left({\lambda_{3}}^{2}+{\lambda_{2}}^{2}{x^{0}}^{-\alpha}/(\alpha-1)\right)+{\lambda_{1}}x^{1}+{\lambda_{2}}x^{2}+{\lambda_{3}}{x^{3}}+F(\lambda_{1},\lambda_{2},\lambda_{3}),$
(38)

where $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are the independent
constants and $F(\lambda_{1},\lambda_{2},\lambda_{3})$ is an arbitrary
function of parameters.

##  4 Conclusion

Two classes of exact solutions are obtained in $R^{2}$ gravity with a
cosmological constant and pure radiation for spatially homogeneous
cosmological models allowing the integration of the equations of motion of
test particles by the method of separation of variables in the Hamilton-Jacobi
formalism.

The solutions belong to the class II of Shapovalov wave-like spaces and have a
separating wave variable in a privileged coordinate system, which reflects the
wave nature of the solutions obtained. The positiveness of the energy density
of pure radiation imposes restrictions on the parameters of the considered
spacetime models.

For the obtained models, the equations of motion of test particles were
integrated in the Hamilton-Jacobi formalism. Complete integrals are presented
for the action function of test particles and the eikonal function. Parametric
dependence of the coordinates of test particles as functions of proper time is
obtained. The solutions obtained make it possible to go over to a synchronous
coordinate system with respect to which the freely falling observer is at
rest.

The obtained models belong to type III according to Bianchi’s classification
and to type N according to Petrov’s classification.

Wave-like spatially homogeneous models of spacetime can describe aperiodic
primordial gravitational waves of the Universe.

## Acknowledgments

The reported study was funded by RFBR, project number N 20-01-00389 A.

## References

  * [1] S. Capozziello, M. De Laurentis, Extended Theories of Gravity, Physics Reports, 509 (4-5), pp. 167-321, DOI: 10.1016/j.physrep.2011.09.003
  * [2] S. Nojiri and S.D. Odintsov, Introduction to modified gravity and gravitational alternative for dark energy, Int.J.Geom.Meth.Mod.Phys. 4 (2007), 115-146.
  * [3] S. Nojiri and S.D. Odintsov, Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models, Physics Reports, 505 (2011), 59-144.
  * [4] S. Nojiri, S.D. Odintsov and V.K. Oikonomou, Modified gravity theories on a nutshell: Inflation, bounce and late-time evolution, Physics Reports (2017), 692, pp. 1-104, DOI: 10.1016/j.physrep.2017.06.001
  * [5] V.N. Shapovalov, Symmetry and separation of variables in Hamilton-Jacobi equation, Izv. vuzov. Fizika (Sov. Phys. Journ.) 9 (1978) 18.
  * [6] V.N. Shapovalov, The Stäckel spaces, Sib. Math. Journal (Sov. J. of Math.) 20 (1979) 1117.
  * [7] K. Osetrin, E. Osetrin, Shapovalov wave-like spacetimes (2020) Symmetry, 12 (8), 1372. DOI: 10.3390/SYM12081372
  * [8] V.G. Bagrov, V.V Obukhov, Classes of Exact Solutions of the Einstein‐Maxwell Equations (1983) Annalen der Physik, 495 (4-5), pp. 181-188. DOI: 10.1002/andp.19834950402
  * [9] V.G. Bagrov, V.V Obukhov, A.V. Shapovalov, Special Stäckel electrovac spacetimes (1986) Pramana, 26 (2), pp. 93-108. DOI: 10.1007/BF02847629
  * [10] V.G. Bagrov, V.V. Obukhov, K.E. Osetrin, Classification of null-Stäckel electrovac metrics with cosmological constant (1988) General Relativity and Gravitation, 20 (11), pp. 1141-1154. DOI: 10.1007/BF00758935
  * [11] V.V. Obukhov, Separation of variables in Hamilton–Jacobi equation for a charged test particle in the Stäckel spaces of type (2.1), International Journal of Geometric Methods in Modern Physics (2020), 17(14), 2050186. DOI: 10.1142/S0219887820501868
  * [12] V.V. Obukhov, Hamilton–Jacobi equation for a charged test particle in the Stäckel spaces of type (2.0), Symmetry (2020), 12 (8), 1289. DOI: 10.3390/sym12081289
  * [13] K. Osetrin, A Filippov and E. Osetrin, The spacetime models with dust matter that admit separation of variables in Hamilton-Jacobi equations of a test particle. Mod. Phys. Lett. A31 (2016), N 06, 1650027.
  * [14] V.G. Bagrov, A.D. Istomin, V.V. Obukhov and K.E. Osetrin, Classification of conformal Stäckel spaces in the Vaidya problem (1996) Russian Physics Journal, 39 (8), pp. 744-749. DOI: 10.1007/BF02437084
  * [15] V.V Obukhov, K.E. Osetrin, A.E. Filippov, Y.A. Rybalov, The Vaidya problem in conformally flat Stäckel spaces of type (1.1) (2009) Russian Physics Journal, 52 (1), pp. 11-14. DOI: 10.1007/s11182-009-9198-3
  * [16] E. Osetrin and K. Osetrin, Pure radiation in space-time models that admit integration of the eikonal equation by the separation of variables method. J. Math. Phys. 58 (2017), N 11, 112504.
  * [17] K. Osetrin, A. Filippov and E. Osetrin, Wave-like spatially homogeneous models of Stäckel spacetimes (2.1) type in the scalar-tensor theory of gravity (2020) Modern Physics Letters A, 2050275. DOI: 10.1142/S0217732320502752
  * [18] E.K. Osetrin, K.E. Osetrin, A.E. Filippov and I.V. Kirnos, Wave-like spatially homogeneous models of Stäckel spacetimes (3.1) type in the scalar-tensor theory of gravity (2020) International Journal of Geometric Methods in Modern Physics, v 17, N 12, 2050184 DOI: 10.1142/S0219887820501844
  * [19] K.E. Osetrin, A.E. Filippov, E.K. Osetrin, Models of Generalized Scalar-Tensor Gravitation Theories with Radiation Allowing the Separation of Variables in the Eikonal Equation (2018) Russian Physics Journal, 61 (8), pp. 1383-1391. DOI: 10.1007/s11182-018-1546-8
  * [20] V.G. Bagrov, V.V Obukhov, New method of integration for the Dirac equation on a curved space-time (1992) Journal of Mathematical Physics, 33 (6), pp. 2279-2289. DOI: 10.1063/1.529600
  * [21] K.E. Osetrin, Y.A. Rybalov, Cosmological Models with Scalar and Spinor Fields (2013) Russian Physics Journal, 55 (12), pp. 1416-1424. DOI: 10.1007/s11182-013-9975-x
  * [22] V.V Obukhov, K.E. Osetrin, A.E. Filippov, Metrics of homogeneous spaces admitting (3.1)-type complete sets (2002) Russian Physics Journal, 45 (1), pp. 42-48. DOI: 10.1023/A:1016093620137
  * [23] V.V. Obukhov, K.E. Osetrin, A.E. Filippov, Homogeneous spacetimes and separation of variables in the Hamilton-Jacobi equation. Journal of Physics A: Mathematical and General. 39 (2006), N 21, 6641-6647.
  * [24] E.K. Osetrin, K.E. Osetrin, A.E. Filippov, Spatially Homogeneous Conformally Stäckel Spaces of Type (3.1) (2020) Russian Physics Journal, 63 (3), pp. 403-409. DOI: 10.1007/s11182-020-02050-2
  * [25] E.K. Osetrin, K.E. Osetrin, A.E. Filippov, Spatially Homogeneous Models Stäckel Spaces of Type (2.1) (2020) Russian Physics Journal, 63 (3), pp. 410-419. DOI: 10.1007/s11182-020-02051-1
  * [26] V.D. Ivashchuk, A.A. Kobtsev, On exponential cosmological type solutions in the model with Gauss-Bonnet term and variation of gravitational constant (2015) European Physical Journal C, 75 (5), 177, pp. 1-12. DOI: 10.1140/epjc/s10052-015-3394-9
  * [27] V.D. Ivashchuk, On stable exponential solutions in Einstein-Gauss-Bonnet cosmology with zero variation of G (2016) Gravitation and Cosmology, 22 (4), pp. 329-332. DOI: 10.1134/S0202289316040095