# Retrieving Event-related Human Brain Dynamics from Natural Sentence Reading

Xinping Liu
School of ICT
University of Tasmania
Sandy Bay, TAS 7005
<EMAIL_ADDRESS>
&Zehong Cao
School of ICT
University of Tasmania
Sandy Bay, TAS 7005
<EMAIL_ADDRESS>

###### Abstract

Electroencephalography (EEG) signals recordings when people reading natural
languages are commonly used as a cognitive method to interpret human language
understanding in neuroscience and psycholinguistics. Previous studies have
demonstrated that the human fixation and activation in word reading associated
with some brain regions, but it is not clear when and how to measure the brain
dynamics across time and frequency domains. In this study, we propose the
first analysis of event-related brain potentials (ERPs), and event-related
spectral perturbations (ERSPs) on benchmark datasets which consist of
sentence-level simultaneous EEG and related eye-tracking recorded from human
natural reading experiment tasks. Our results showed peaks evoked at around
162 ms after the stimulus (starting to read each sentence) in the occipital
area, indicating the brain retriving lexical and semantic visual information
processing approaching 200 ms from the sentence onset. Furthermore, the
occipital ERP around 200ms presents negative power and positive power in short
and long reaction times. In addition, the occipital ERSP around 200ms
demonstrated increased high gamma and decreased low beta and low gamma power,
relative to the baseline. Our results implied that most of the semantic-
perception responses occurred around the 200ms in alpha, beta and gamma bands
of EEG signals. Our findings also provide potential impacts on promoting
cognitive natural language processing models evaluation from EEG dynamics.

##  1 Introduction

Preserving human languages by text has been influencing the human cognitive
growth, and it is also considered as the primary medium for preservation of
natural language processing (NLP). However, research on how human brains
retrieve and process the meaning from text-based languages is relatively few.
Previous research Mason and Just (2006) has used custom reading materials to
explore particular facets of linguistic material such as event-related
potentials (ERPs) and eye motions. Later on, using Electroencephalography
(EEG) in conjunction with eye-tracking has become significant measures to
study the temporal dynamics of natural sentence reading, because EEG signals
have excellent time resolution and comparatively low costs to be recorded
Dimigen et al. (2011), Loberg et al. (2018), Hollenstein and Zhang (2019). For
example, a recent study Pfeiffer et al. (2020) indicated that text-based
emotion analyze could be significantly enhanced by computing fixation-related
potentials (FRPs) as one of ERPs on a neurocognitive dataset named ZuCo 1.0
Hollenstein et al. (2018) that recorded from human readers during naturalistic
reading. To measure event-related synchronisations, Roach and Mathalon (2008)
quantified the timeline of event-related spectral perturbations (ERSPs) that
could generate brain dynamic changes in amplitude of the broadband EEG
frequency spectrum.

However, the studies above only focused on wording or naive reading tasks
without comprehensive event-related brain dynamic analyze, so it is still
unknown how ERP and ERPS changes in human natural sentence reading. In this
work, we investigated EEG brain dynamics of ERPs and ERSPs through both
channel and component forms for the total of 18 subjects during a recent
natural sentence reading task, ZuCo 2.0 (Hollenstein et al., 2020). This task
supports to analyse the differences in cognitive processing between natural
reading, and also allows us to precisely extract the EEG signals for sentence-
level processing and detecting human retrieving lexical and semantic
information in natural sentence reading.

Figure 1: Natural Sentence Reading Experiment. (a) Natural reading phases. (b)
Each sentence onset and finished state and switching to the next sentence by
control pad. The sentence reading and sentence switching were performed
alternately, and each sentence switching lasted 100 ms, while sentence reading
lasted from 1.2 second to 12 seconds. (c) Location map of multi-channel EEG.

##  2 Related Work

ERPs Frank et al. (2013, 2015) are electrical potential responses or reaction
measured in relation to an event, such as reading each sentence in our study.
For example, N200 is a negative wave in connection with a baseline because it
occurs at around 200 to 350 ms after the onset of the first word in a
sentence, but it may cause some mismatch issues after an infrequent change in
the case of visual stimuli Brattico et al. (2006). The recent study Pfeiffer
et al. (2020) has indicated that text-based emotion analyze can be
significantly enhanced by using FRPs to analyze human readers during
naturalistic reading, but it only focused on the potentials of word-based
reading. Also, Hale et al. (2018), Jat et al. (2019) applied ERPs as a
neuroimage feature to retrieve the lexical and semantic information from the
natural language content for designing cognitive word embedding evaluation in
NLP, but still, there has some significant loss possibly due to not precise
patterns. In terms of ERSPs, one study zhang2017investigation investigated
emotion and working memory analyze to reading, but it did not touch the
natural sentence reading tasks.

##  3 Method

###  3.1 Experiment Paradigm and Data

The natural sentence reading experiment and data was performed and collected
in ZuCo 2.0 (Hollenstein et al., 2020), which recorded from 18 healthy
subjects who are all native English speakers originating from Canada, USA, UK
or Australia. The sentences were displayed one by one at the same position on
the monitor during the natural reading paradigm. The text was displayed on a
light grey backdrop in black with font size 20-point Arial, resulting in a
letter height of 0.8 mm or 0.674 ∘. During the entire reading task, a sum of
80 letters or 13 words were presented per line, and the subjects who
participate in this natural reading experiment read each sentence at their own
speed.

In the task, the 349 sentences were displayed one by one at the same position
on the monitor during the natural reading task for each subject. The recording
devices include both 128-channel EEG and eye-tracking that used to collect
brain signals and eye motions when subjects conducted natural sentence
reading, as shown in Fig. 1. In this work, one trial was considered as reading
each sentence, and there are a total of 6,282 trials (18 subjects multiplied
by 349 sentences in the natural reading paradigm).

###  3.2 Data Processing and analyze

#### EEG Processing

The raw EEG data were loaded by EEGLAB (Delorme and Makeig, 2004) with the
removal of bad electrodes, keeping 105 EEG channels. We then removed
electromyogram (EMG) and electrooculogram (EOG) artifacts and minimised other
artifacts by setting a finite impulse response filter 1-50Hz. Next, we
conducted the EEG electrodes re-referencing by the vertex electrode Cz, and
the re-referenced EEG data were measured by Z Score formula normalisation.
Furthermore, we conducted the data segmentation to extract all trials where
two time-locking events (sentence onset event and sentence finished event)
were selected. What is more, the trial limits were 100 ms before starting
sentence reading (-100 ms, act as the baseline) to 1200 ms after the reading
event (+1200 ms), because the sentence interval sets 100 ms and the minimum
interval between sentence onset and sentence finished is 1223 ms.

In addition, to represent synchronous brain activities, decomposing the pre-
processed data by Independent Component analyze (ICA) is an effective approach
to investigate the diversity of source information typically contained in EEG
data. In this study, we merged a total of 18 pre-processed single-subject EEG
data and conducted the statistical comparisons by running the merged-subjects
ICA decomposition to find the potential similar independent component sources
across 18 subjects.

#### EEG analyze: ERP/ERSP in channel-based and component-based forms

The 1-D event latency scalp maps and 2-D ERP/ERSP images could produce to
intuitively discover the brain dynamics of human natural sentence reading in
both forms of channel-based and component-based EEG patterns, respectively. In
terms of ERPs, to analyse the channel-based ERPs, we chose the most relevant
channel to the corresponding events from the scalp map where the power is
concentrated among a total of 18 subjects. Meanwhile, to analyse component-
based ERPs, we calculated and selected the power spectrum for isolated
components to identify which component contributes the most to the ERPs. For
ERSPs, similarly, we determined the channel which is most relevant to the
corresponding events and analysed spectral perturbations, and then calculated
and selected a component that the source contributes the most to the task.
Please note that the significance level sets 0.01 with removing the baseline
for plotting the ERSPs in this study.

##  4 Results

Our results show the comprehensive group analyze and visualisation of 18
subjects of EEG dynamics for natural sentence reading, including 1-D event
latency potentials plots and 2-D ERP and ERSP images in channel-based and
component-based forms.

#### Channel/Component 1-D event latency potentials

As shown in Fig. 2, from 1-D event latency potentials among all subjects, we
found highest peaks at 162 ms from the onset natural sentence reading, which
may indicate that a significant EEG dynamic (particular in reduced power in
the occipital region) approaching 200 ms (N200) and this finding also
supported by evidence for the early latency of lexical and semantic
information retrieval in visual word recognition (Hauk et al., 2012, Pfeiffer
et al., 2020). Thus, we assume that the human start retrieves lexical and
semantic information near-simultaneously at the central occipital region (Oz
channel) within 200 ms of the sentence reading onset.

Figure 2: Channel-based 1-D event latency potientials with the scalp map among
a total of 18 subjects for natural sentence reading.

Following above, we also extracted 7 most significant component maps from 1-D
event latency potentials, as shown in Fig. 3, which showed the occipital
component contributed most to the task (marked as rank 1). The results in the
section are used for determining that the occipital region is the most
contribution area related to the natural sentence reading task, because the
occipital results displayed the strongest response in the text-based retrieval
of sentiment information than that of the other regions of the brain.

Figure 3: Components-based 1-D event latency potientials with the scalp map
among a total of 18 subjects for natural sentence reading.

#### Channel/Component 2-D ERP Results

As shown in Fig. 4-A, we located the channel Oz as the central occipital
region to plot channel-based ERP results, which present the negative power
(marked as the blue colour) around 200 ms in short reaction time (between
1-2000 trials). However, from the occipital component-based ERP image as shown
in Fig. 4-B, we observed the positive power (marked as the red colour) around
200 ms in medium (between 2000-4000 trials) and long (between 4000-6000
trials) reaction times. This is a new finding that may be beneficial for
retrieving the diverse patterns of brain visual reading processing.

Figure 4: 2-D ERP results for natural sentence reading in channel-based (Oz)
and component-based (occipital) forms.

#### Channel/Component 2-D ERSP Results

In terms of the channel-based and component-based ERSPs in the occipital
region, they demonstrated a similar performance around 200 ms, as shown in
Fig. 5. The Oz-channel ERSP as presented in Fig. 5-A, at the time of around
200ms, demonstrated increased high alpha (10-12Hz), high beta (25-30Hz), and
high gamma (40-50Hz) power, and decreased low beta (13-25Hz) and low gamma
(30-40Hz) power, relative to the baseline (the interval of sentence reading).
In terms of the occipital component-channel ERSP around 200 ms, Fig. 5-B
presented increased high alpha (10-12Hz) and high gamma (40-50Hz) power, and
decreased low beta (13-25Hz) and low gamma (30-40Hz) power, relative to the
baseline (the interval of sentence reading).

Furthermore, ERSP results in the occipital region around 200 ms primarily
consist with the longer latencies covering from 300 ms to 1200 ms (the ending
point of each sentence reading), particularly in the enhanced high gamma
(40-50Hz) and the reduced low beta (13-25Hz) and low gamma (30-40Hz) power
spectral, relative to the baseline, suggesting the patterns for retrieving
lexical and semantic information in the natural sentence reading task.

Figure 5:  2-D ERSP results for natural sentence reading in channel-based (Oz)
and component-based (occipital) forms.

##  5 Conclusion

Our study successfully identified the spatiotemporal neural dynamics of
sentiment processing during the naturalistic reading of English sentences.
Combining high-density EEG and eye-tracking data, we provide clear evidence
that essential processes of lexical and semantic information retrieval occur
around the 200 ms of each sentence onset, by retrieving EEG patterns of
negative power and positive power in short and long reaction times and
increased high gamma and decreased low beta and low gamma power. We believe
that these EEG dynamics found from this study for human language understanding
could benefit for semantic representations of sentence reading and cognitive
evaluation modelling for NLP.

## References

  * Mason and Just (2006) Robert A. Mason and Marcel Adam Just.  Chapter 19 - neuroimaging contributions to the understanding of discourse processes.  In Matthew J. Traxler and Morton A. Gernsbacher, editors, _Handbook of Psycholinguistics (Second Edition)_ , pages 765 – 799. Academic Press, London, second edition edition, 2006.  ISBN 978-0-12-369374-7.  doi: https://doi.org/10.1016/B978-012369374-7/50020-1.  URL http://www.sciencedirect.com/science/article/pii/B9780123693747500201.
  * Dimigen et al. (2011) Olaf Dimigen, Werner Sommer, Annette Hohlfeld, Arthur M Jacobs, and Reinhold Kliegl.  Coregistration of eye movements and eeg in natural reading: analyses and review.  _Journal of experimental psychology. General_ , 140(4):552—572, November 2011.  ISSN 0096-3445.  doi: 10.1037/a0023885.  URL https://doi.org/10.1037/a0023885.
  * Loberg et al. (2018) Otto Loberg, Jarkko Hautala, Jarmo A. Hämäläinen, and Paavo H. T. Leppänen.  Semantic anomaly detection in school-aged children during natural sentence reading – a study of fixation-related brain potentials.  _PLoS ONE_ , 13(12):e0209741, 2018.  doi: 10.1371/journal.pone.0209741.  URL https://app.dimensions.ai/details/publication/pub.1110936695andhttps://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0209741&type=printable.
  * Hollenstein and Zhang (2019) Nora Hollenstein and Ce Zhang.  Entity recognition at first sight improving NER with eye movement information.  _CoRR_ , abs/1902.10068, 2019.  URL http//arxiv.org/abs/1902.10068.
  * Pfeiffer et al. (2020) Christian Pfeiffer, Nora Hollenstein, Ce Zhang, and Nicolas Langer.  Neural dynamics of sentiment processing during naturalistic sentence reading.  _NeuroImage_ , 218:116934, 2020.  ISSN 1053-8119.  doi: https://doi.org/10.1016/j.neuroimage.2020.116934.  URL http://www.sciencedirect.com/science/article/pii/S1053811920304201.
  * Hollenstein et al. (2018) Nora Hollenstein, Jonathan Rotsztejn, Marius Troendle, Andreas Pedroni, Ce Zhang, and Nicolas Langer.  Zuco, a simultaneous eeg and eye-tracking resource for natural sentence reading.  _Scientific Data_ , 5, 2018.
  * Roach and Mathalon (2008) Brian J. Roach and Daniel H. Mathalon.  Event-related eeg time-frequency analysis: an overview of measures and an analysis of early gamma band phase locking in schizophrenia.  _Schizophrenia bulletin_ , 34(5):907–926, Sep 2008.
  * Hollenstein et al. (2020) Nora Hollenstein, Marius Troendle, Ce Zhang, and Nicolas Langer.  ZuCo 2.0: A dataset of physiological recordings during natural reading and annotation.  In _Proceedings of the 12th Language Resources and Evaluation Conference_ , pages 138–146, Marseille, France, May 2020. European Language Resources Association.  ISBN 979-10-95546-34-4.  URL https://www.aclweb.org/anthology/2020.lrec-1.18.
  * Frank et al. (2013) Stefan L. Frank, Irene Fernandez Monsalve, Robin L. Thompson, and Gabriella Vigliocco.  Reading time data for evaluating broad-coverage models of english sentence processing.  _Behavior Research Methods_ , 45(4):1182–1190, 2013.
  * Frank et al. (2015) Stefan L. Frank, Leun J. Otten, Giulia Galli, and Gabriella Vigliocco.  The erp response to the amount of information conveyed by words in sentences.  _Brain and Language_ , 140:1–11, 2015.
  * Brattico et al. (2006) Elvira Brattico, Mari Tervaniemi, Risto Naatanen, and Isabelle Peretz.  Musical scale properties are automatically processed in the human auditory cortex.  _Gene Expression Patterns_ , 1117:162–174, October 2006\.  ISSN 1567-133X.  doi: 10.1016/j.brainres.2006.08.023.
  * Hale et al. (2018) John Hale, Chris Dyer, Adhiguna Kuncoro, and Jonathan Brennan.  Finding syntax in human encephalography with beam search.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 2727–2736, Melbourne, Australia, July 2018. Association for Computational Linguistics.  doi: 10.18653/v1/P18-1254.  URL https://www.aclweb.org/anthology/P18-1254.
  * Jat et al. (2019) Sharmistha Jat, Hao Tang, Partha P. Talukdar, and Tom Michael Mitchell.  Relating simple sentence representations in deep neural networks and the brain.  _CoRR_ , abs/1906.11861, 2019.  URL http//arxiv.org/abs/1906.11861.
  * Delorme and Makeig (2004) Arnaud Delorme and Scott Makeig.  Eeglab: an open source toolbox for analysis of single-trial eeg dynamics including independent component analysis.  _Journal of neuroscience methods_ , 134(1):9—21, March 2004.  ISSN 0165-0270.  doi: 10.1016/j.jneumeth.2003.10.009.  URL https://doi.org/10.1016/j.jneumeth.2003.10.009.
  * Hauk et al. (2012) O. Hauk, C. Coutout, A. Holden, and Y. Chen.  The time-course of single-word reading: Evidence from fast behavioral and brain responses.  _NeuroImage_ , 60(2):1462–1477, 2012.  doi: 10.1016/j.neuroimage.2012.01.061.  URL https://doi.org/10.1016/j.neuroimage.2012.01.061.


# Thermodynamic of the $f(Q)$ universe

Haomin Rao<EMAIL_ADDRESS>School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  University of Chinese Academy of Sciences, 100190 Beijing,
China  Chunhui Liu<EMAIL_ADDRESS>School of Fundamental
Physics and Mathematical Sciences, Hangzhou Institute for Advanced Study,
UCAS, Hangzhou 310024, China  University of Chinese Academy of Sciences,
100190 Beijing, China  Institute of Theoretical Physics, Chinese Academy of
Sciences, Beijing 100190, China  Chao-Qiang Geng<EMAIL_ADDRESS>School of
Fundamental Physics and Mathematical Sciences, Hangzhou Institute for Advanced
Study, UCAS, Hangzhou 310024, China  Department of Physics and Synergetic
Innovation Center for Quantum Effects and Applications, Hunan Normal
University, Changsha, Hunan 410081, China

###### Abstract

nothing

We investigate thermodynamics of apparent horizon in the $f(Q)$ universe with
trivial and nontrivial connections. We first explore the perspectives of the
first law, generalized second law and $P-V$ phase transition with trivial
connection. We show that the lowest-order correction of entropy has the same
form as that in loop quantum gravity, and the critical exponents of the phase
transition caused by the lowest-order correction are consistent with those in
mean field theory. We then examine the thermodynamic implication of nontrivial
connections. We find that nontrivial connections in the $f(Q)$ universe imply
non-equilibrium states from the perspective of thermodynamics.

##  I Introduction

The four laws of black hole mechanics in general relativity (GR) are very
similar to the four laws of thermodynamics [1]. Classically, black holes seem
to have neither temperature nor entropy. However, the area law of black hole
entropy can be derived from information theory [2] and the black hole
temperature can be derived from quantum field theory in curved spacetime [3].
This suggests that the black hole thermodynamics is not just an analogy, but
reveals a hidden relationship between gravity, thermodynamics and quantum
theory [4, 5, 6, 7]. Following this idea, more thermodynamic concepts such as
phase transition have been introduced into the study of black holes [8, 9,
10]. In addition, it has been found that black holes in different gravity
models have different thermodynamic behaviors [11, 12, 13, 14], indicating
that black hole thermodynamics contains information about gravity models.

The universe is another neat physical system closely related to gravity.
Inspired by the black hole thermodynamics, people found that the evolution
equations of the universe can be written in the form of thermodynamic laws
[15, 16]. Since the temperature of the horizon in the universe can also be
derived by quantum methods [17, 18, 19, 20], the thermodynamics of the
universe once again reveals the hidden relationship between gravity and
quantum theory. Even more surprising is that in some gravity models the
horizon entropy of the black hole and the horizon entropy of the universe have
exactly the same form [21, 22, 23]. This reflects that the thermodynamics of
the universe and the thermodynamics of black holes are not independent in a
given gravity model. Clearly, it is necessary and attractive to discuss the
thermodynamics of the universe in various gravity models [24, 25, 26, 27, 28,
29, 30, 31].

Symmetric teleparallel gravity (STG) is a newly popular modified gravity
scheme that identifies gravity as non-metricity rather than curvature or
torsion. The simplest STG model classically is equivalent to GR [32, 33, 34],
which provides us with another way to modify GR, that is, to modify the
simplest STG model. The most eye-catching and interesting modified STG model
is the $f(Q)$ model [35, 36], and its cosmological applications have been
extensively studied in the literature [37, 38, 39, 40, 41, 42]. More recent
studies have shown that in the $f(Q)$ model, there are three different
connections that satisfy cosmological symmetry in the flat universe [43, 44,
45, 46]. And different connections can result in different background
evolutions as well as different perturbation behaviors [47, 48, 49, 50]. This
makes the flat universe in the $f(Q)$ model more complex but with richer
phenomenology.

So far, the thermodynamics of the $f(Q)$ universe and the thermodynamic
implication of different connections remain to be explored. In this paper, we
will investigate the thermodynamics of $f(Q)$ universes with different
connections. Since the case of trivial connection has attracted overwhelming
attention in past studies of the $f(Q)$ universe, the thermodynamics of the
$f(Q)$ universe with trivial connection is also the primary target in this
paper. We will explore it from the perspectives of the first law, generalized
second law and $P-V$ phase transition. After that, we will examine the
thermodynamic implication of nontrivial connections. This will provides the
first glimpse into how different connections affect the thermodynamics of the
universe in the same gravity model.

The present paper is organized as follows. In section II, we briefly introduce
the $f(Q)$ gravity and its flat universe background with trivial and
nontrivial connections. In section III, we analyse the thermodynamics of the
flat universe with trivial connection from the following three aspects. In
subsection III.1, we examine the first law of thermodynamics and derive the
area law of entropy for the most general $f(Q)$ model. Subsequently, we
briefly discuss the restrictions on the $f(Q)$ universe brought by the
generalized second law in subsection III.2. In subsection III.3, we study the
$P-V$ phase transition in a simple $f(Q)$ model and calculate all critical
exponents. Finally, in section IV, we show the thermodynamic implication of
nontrivial connections from the perspective of the first law.

In this paper, we adopt the unit $G=1$ and the signature $(-,+,+,+)$. The
spacetime indices are denoted by Greek indices $\mu,\nu,\rho,...=0,1,2,3$ and
the spatial indices are represented by Latin indices $i,j,k,...=1,2,3$. In
addition, we distinguish the spacetime affine connection
$\Gamma^{\rho}_{\leavevmode\nobreak\ \mu\nu}$ and its associated covariant
derivative $\nabla$ from the Levi-Civita connection
$\mathring{\Gamma}^{\rho}_{\leavevmode\nobreak\ \mu\nu}$ and its associated
covariant derivative $\mathring{\nabla}$, respectively.

##  II $f(Q)$ gravity and its cosmology

The so-called STG theory is formulated in a spacetime endowed with a metric
$g_{\mu\nu}$ and an affine connection $\Gamma^{\rho}_{\leavevmode\nobreak\
\mu\nu}$, which is curvature-free and torsion-free. Without curvature and
torsion, the gravitation is identified with non-metircity

$Q_{\rho\mu\nu}=\nabla_{\rho}g_{\mu\nu}.$ (1)

The connection $\Gamma^{\rho}_{\leavevmode\nobreak\ \mu\nu}$ on such a
geometry can be generally expressed as

$\Gamma^{\rho}_{\leavevmode\nobreak\ \mu\nu}=\frac{\partial x^{\rho}}{\partial
y^{\sigma}}\partial_{\mu}\partial_{\nu}y^{\sigma}.$ (2)

As a result, we can regard $g_{\mu\nu}$ and $y^{\mu}$ as the basic variables
of the STG theory.

The simplest STG model is the so-called symmetric teleparallel equivalent of
general relativity (STEGR) model, whose action is

$I_{\text{STEGR}}=\frac{1}{16\pi}\int\mathrm{d}^{4}x\sqrt{-g}\,Q+I_{m},$ (3)

with the non-metricity scalar $Q$ defined as

$Q=P^{\rho\mu\nu}Q_{\rho\mu\nu}=-\frac{1}{4}Q_{\rho\mu\nu}Q^{\rho\mu\nu}+\frac{1}{2}Q_{\rho\mu\nu}Q^{\mu\nu\rho}+\frac{1}{4}Q_{\mu}Q^{\mu}-\frac{1}{2}Q_{\mu}\tilde{Q}^{\mu},$
(4)

where the non-metricity conjugate is
$P^{\rho\mu\nu}=-\frac{1}{4}Q^{\rho\mu\nu}+\frac{1}{2}Q^{(\mu\nu)\rho}+\frac{1}{4}(Q^{\rho}-\tilde{Q}^{\rho})g^{\mu\nu}-\frac{1}{4}g^{\rho(\mu}Q^{\nu)}$
with $Q_{\mu}=Q_{\mu\nu}{}^{\nu}$, $\tilde{Q}_{\mu}=Q^{\nu}{}_{\nu\mu}$. Since
we have the identity $\mathring{R}\equiv
Q-\mathring{\nabla}_{\mu}(Q^{\mu}-\tilde{Q}^{\mu})$, the action (3) is
identical to the Einstein-Hilbert action up to a surface term, where the
curvature scalar $\mathring{R}$ is defined by the Levi-Civita connection.
Since the surface term in the action does not affect the equations of motion,
we say that STEGR is equivalent to GR at the level of equations of motion.

This equivalence provides us with another way to modify GR, which is to modify
the STEGR model within the STG framework. Along this line, a variety of
modified STG models have been proposed. The most interesting one is the $f(Q)$
model, which generalizes $Q$ in the action (3) to a smooth function $f(Q)$,
given by

$I=\frac{1}{16\pi}\int\mathrm{d}^{4}x\sqrt{-g}\,f(Q)+I_{m},$ (5)

The equations of motion follow from the variations with respect to
$g_{\mu\nu}$ and $y^{\mu}$ are

$\displaystyle
f_{Q}\mathring{G}^{\mu\nu}+\frac{1}{2}(Qf_{Q}-f)g^{\mu\nu}+2\nabla_{\rho}f_{Q}P^{\rho\mu\nu}$
$\displaystyle=$ $\displaystyle 8\pi T^{\mu\nu},$ (6)
$\displaystyle\nabla_{\mu}\nabla_{\nu}\left(\sqrt{-g}f_{Q}P^{\mu\nu}{}_{\rho}\right)$
$\displaystyle=$ $\displaystyle 0,$ (7)

where $f_{Q}=\frac{\mathrm{d}f(Q)}{\mathrm{d}Q}$,
$\mathring{G}^{\mu\nu}=\mathring{R}^{\mu\nu}-\frac{1}{2}\mathring{R}g^{\mu\nu}$
is the Einstein tensor, $T^{\mu\nu}=\frac{2}{\sqrt{-g}}\frac{\delta
I_{m}}{\delta g_{\mu\nu}}$ is the energy-momentum tensor of matter. It can be
verify that if $f(Q)=Q$, Eq. (6) becomes the Einstein field equation, while
Eq. (7) is reduced to an identity.

After clarifying the $f(Q)$ model, let’s briefly introduce its cosmology. In
the flat universe, the metric can be expressed in rectangular coordinate as

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\delta_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j},$
(8)

where $a=a(t)$ is the scale factor. Eq. (8) comes from the symmetry
requirement that the metric is homogeneous and isotropic, i.e.
$\mathcal{L}_{\xi}g_{\mu\nu}=0$, where $\mathcal{L}_{\xi}$ is the Lie
derivative and $\xi$ represents all six Killing vector fields in the flat
universe. In the STG theory, the connection cannot be completely determined by
the metric. As suggested in Refs. [43, 44, 45, 46], it is natural to further
require that the connection is also homogeneous and isotropic, that is,

$\mathcal{L}_{\xi}{\Gamma}^{\rho}_{\leavevmode\nobreak\
\mu\nu}={\nabla}_{\mu}{\nabla}_{\nu}\,\xi^{\rho}=0.$ (9)

It can be obtained from Eq. (9) that the non-zero components of the connection
${\Gamma}^{\rho}_{\leavevmode\nobreak\ \mu\nu}$ are

$\Gamma^{0}_{\leavevmode\nobreak\ 00}=K_{1}\leavevmode\nobreak\
,\leavevmode\nobreak\ \Gamma^{0}_{\leavevmode\nobreak\
ij}=a^{2}K_{2}\,\delta_{ij}\leavevmode\nobreak\ ,\leavevmode\nobreak\
\Gamma^{i}_{\leavevmode\nobreak\ 0j}=\Gamma^{i}_{\leavevmode\nobreak\
j0}=K_{3}\,\delta_{ij},$ (10)

with $\\{K_{1},K_{2},K_{3}\\}$ having three branch solutions, given by

$\displaystyle\text{branch 1 :}\leavevmode\nobreak\
K_{1}=\gamma\leavevmode\nobreak\ ,\leavevmode\nobreak\
K_{2}=0\leavevmode\nobreak\ ,\leavevmode\nobreak\ K_{3}=0,$ (11)
$\displaystyle\text{branch 2 :}\leavevmode\nobreak\
K_{1}=\dot{\gamma}/\gamma+\gamma\leavevmode\nobreak\ ,\leavevmode\nobreak\
K_{2}=0\leavevmode\nobreak\ ,\leavevmode\nobreak\ K_{3}=\gamma,$ (12)
$\displaystyle\text{branch 3 :}\leavevmode\nobreak\
K_{1}=-(\dot{\gamma}/\gamma+2H),\leavevmode\nobreak\
K_{2}=\gamma\leavevmode\nobreak\ ,\leavevmode\nobreak\ K_{3}=0,$ (13)

where $H=\dot{a}/a$ is the Hubble rate, $\gamma=\gamma(t)$ is a function of
time $t$ and the superscript “dot” represents the derivative with respect to
$t$. It can be seen that even if the metric and connection have been required
to be homogeneous and isotropic, the form of the cosmological background is
not unique. This is a unique feature of the STG universe, which does not
appear in $f(R)$ and $f(T)$.

Putting the metric (8) and the connection (10) into Eqs. (6) and (7), we
obtain the background equations as

$\displaystyle
3f_{Q}H^{2}+\frac{1}{2}(f-Qf_{Q})+\frac{3}{2}(K_{3}-K_{2})\dot{Q}f_{QQ}=8\pi\rho_{m},$
(14) $\displaystyle
f_{Q}(2\dot{H}+3H^{2})+\frac{1}{2}(f-Qf_{Q})+\frac{1}{2}(4H-K_{2}-3K_{3})\dot{Q}f_{QQ}=-8\pi
p_{m},$ (15)
$\displaystyle(K_{2}-K_{3})\ddot{f}_{Q}+\left(HK_{2}-3HK_{3}-2K_{1}K_{2}\right)\dot{f}_{Q}=0,$
(16)

where $Q=-6H^{2}$ in branch 1 and $Q=-6H^{2}+9H\gamma+3\dot{\gamma}$ in branch
2 and branch 3, $f_{QQ}=\mathrm{d}f_{Q}/\mathrm{d}Q$, $\rho_{m}$ and $p_{m}$
are the energy density and pressure of matter, respectively. Eqs. (14)-(16)
show that the connection can affect the evolution of the cosmological
background. In fact, this feature has only recently been taken seriously.
Almost all early studies of the $f(Q)$ universe have assumed the trivial
connection $\Gamma^{\rho}_{\leavevmode\nobreak\ \mu\nu}=0$. This is equivalent
to considering only the branch 1 solution with $\gamma=0$, which is only a
special solution of the condition (9). In this case, Eqs. (14) and (15) can be
simplified to

$\displaystyle 3f_{Q}H^{2}+\frac{1}{2}(f-Qf_{Q})=8\pi\rho_{m},$ (17)
$\displaystyle
f_{Q}(2\dot{H}+3H^{2})+\frac{1}{2}(f-Qf_{Q})+2H\dot{f}_{Q}=-8\pi p_{m},$ (18)

while Eq. (16) degenerates into an identity. Eq. (17) and (18) are exactly the
background equations adopted in most studies of the $f(Q)$ universe.

##  III Thermodynamic with trivial connection

After a brief introduction to $f(Q)$ gravity and its flat universe, in this
section, we study the thermodynamics of the $f(Q)$ universe. Since the case of
trivial connection $\Gamma^{\rho}_{\leavevmode\nobreak\ \mu\nu}=0$ has
attracted overwhelming attention in the studies of the $f(Q)$ universe, we
only focus on the thermodynamics of the universe with trivial connection in
this section. We will explore the thermodynamics from the perspectives of the
first law, generalized second law and $P-V$ phase transition. The
thermodynamic effects of nontrivial connections are left to the next section.

###  III.1 First Law of Thermodynamics

Similar to the black hole horizon, there is also thermodynamics at the
cosmological horizon. For the convenience of discussions, we express the
metric (8) as the following form with obvious spherical symmetry

$\mathrm{d}s^{2}=h_{\alpha\beta}\mathrm{d}z^{\alpha}\mathrm{d}z^{\beta}+R^{2}\mathrm{d}\Omega_{2}^{2},$
(19)

where $R=a|\vec{x}|$ is physics radius, $\mathrm{d}\Omega_{2}^{2}$ is the line
element of a 2-dimensional sphere with unit radius,
$z^{\alpha}=(t,|\vec{x}|)$, and $h_{\alpha\beta}=\text{diag}\\{-1,a^{2}\\}$ is
the induced metric on the $z$-plane. Since the event horizon depends on the
whole history of the universe, the apparent horizon is the more natural
horizon, which is a trapped surface with a vanishing expansion of
$h^{\alpha\beta}{\partial_{\alpha}}R{\partial_{\beta}}R=0$. The apparent
horizon of the flat universe can be easily found as

$R_{A}=H^{-1}.$ (20)

which is exactly equal to the Hubble horizon. Then, the surface gravity on the
apparent horizon (20) can be obtained as

$\kappa=\frac{1}{2\sqrt{-h}}\partial_{\alpha}\left(\sqrt{-h}h^{\alpha\beta}\partial_{\beta}R\right)_{R_{A}}=-\frac{1}{R_{A}}(1-\frac{\dot{R}_{A}}{2})$
(21)

Inspired by the thermodynamics of black holes, the temperature of the apparent
horizon in the flat universe is defined by the surface gravity $\kappa$ as

$T=\frac{|\kappa|}{2\pi}=\frac{1}{2\pi
R_{A}}\left(1-\frac{\dot{R}_{A}}{2}\right).$ (22)

The temperature (22) is consistent with the one obtained by the Hamilton-
Jacobi tunneling method for the apparent horizon [19]. Following Refs. [51,
52, 53], we define the work density by
$W=-\frac{1}{2}h_{\alpha\beta}T^{\alpha\beta}$. Using background equations
(17) and (18) we can find the work density in the $f(Q)$ flat universe as

$W=\frac{1}{8\pi}\Big{[}f_{Q}(\dot{H}+3H^{2})+\frac{1}{2}(f-Qf_{Q})+H\dot{f}_{Q}\Big{]}.$
(23)

In addition, the total energy of matter inside the apparent horizon can be
easily obtained from the background equation (17)

$E=\rho_{m}V=\frac{R_{A}f_{Q}}{2}+\frac{R_{A}^{3}}{12}(f-Qf_{Q}),$ (24)

where $V=\frac{4}{3}\pi R_{A}^{3}$ is the thermodynamic volume, which is also
the physical volume within the apparent horizon in the flat universe. Finally,
from the internal energy $U=-E$ and the thermodynamic pressure $P=W$, the
standard first law of thermodynamics can be established

$\mathrm{d}U=T\mathrm{d}S-P\mathrm{d}V,$ (25)

where the entropy $S$ can be integrated as a function of the horizon area

$S(A)=\frac{1}{4}\int\left(f_{Q}+2Qf_{QQ}\right)\mathrm{d}A,$ (26)

with $Q=-24\pi A^{-1}$, where $A=4\pi R_{A}^{2}$ is the thermodynamic area as
well as the physical area of the apparent horizon. This means that the
thermodynamics of the apparent horizon is always the equilibrium
thermodynamics in the $f(Q)$ universe with trivial connection.

In order to explore the horizon entropy (26) more intuitively, let’s take a
look at the entropy (26) in some simple $f(Q)$ models. The simplest $f(Q)$
model is the STEGR model, that is, $f(Q)=Q$. The entropy of the STEGR model
can be obtained by Eq. (26) as

$S(A)=\frac{A}{4},$ (27)

which is exactly the Bekenstein-Hawking entropy. Next, from the perspective of
the Taylor expansion, the $f(Q)$ model retaining the lowest-order correction
is $f(Q)=Q+\lambda\,Q^{2}$. The entropy in this model can be obtained by Eq.
(26) as

$S(A)=\frac{A}{4}-36\pi\lambda{\rm ln}\frac{A}{A_{0}},$ (28)

where $A_{0}$ is an integration constant. Compared with the Bekenstein-Hawking
entropy, the entropy (28) has an additional logarithmic correction. Such a
logarithmic correction can also arise from the lowest-order correction of loop
quantum gravity due to thermal equilibrium fluctuations and quantum
fluctuations [54, 55, 56, 57, 58, 59]. Next, if we continue to follow the
Taylor expansion and keep it to the second-order, then the function $f(Q)$ can
be expressed as $f(Q)=Q+\lambda\,Q^{2}+c\,Q^{3}$. The entropy in this model
can be obtained by Eq. (26) as

$S(A)=\frac{A}{4}-36\pi\lambda\ln\frac{A}{A_{0}}-\frac{2160\pi^{2}c}{A}.$ (29)

Comparing with Eq. (28), it is found that the second-order correction of
entropy is a power law correction. The power law correction of entropy also
appears in the second-order correction of loop quantum gravity [57, 58, 59] as
well as the entanglement of quantum fields between in and out the horizon [60,
61, 62]. A succession of coincidences seems to suggest that the $f(Q)$ model
is linked to the effects of quantum gravity. However, this topic is beyond the
scope of this paper and is left for future research.

###  III.2 Generalized Second Law of Thermodynamics

In the previous subsection, we derive the horizon entropy from the first law
of thermodynamics. It follows that one might wonder what would happen if we
required that the total entropy does not decrease. In this subsection, we
investigate the generalized second law of thermodynamics in the $f(Q)$
universe, which asserts that the entropy of the horizon plus the entropy of
matter within the horizon does not decrease.

Suppose that there are $i=1,2,...,n$ components in the universe, and the
$i$-th component follows continuity equation

$\rho_{i}+3H(\rho_{i}+p_{i})=q_{i},$ (30)

where $\rho_{i}$ and $p_{i}$ are the energy density and pressure of the $i$-th
component, respectively, and $q_{i}$ is the interaction term, which reflects
that energy can be exchanged between different components. The total energy
density $\rho_{m}=\sum\rho_{i}$ and pressure $p_{m}=\sum p_{i}$ must satisfy
the continuity equation, which means

$\dot{\rho}_{m}+3H(\rho_{m}+p_{m})=0\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \Rightarrow\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \sum_{i=1}^{n}q_{i}=0.$ (31)

Based on the discussion of Refs. [63, 64, 65], we apply the first law of
thermodynamics to each component

$\mathrm{d}E_{i}=T_{i}\mathrm{d}S_{i}-p_{i}\mathrm{d}V,$ (32)

where $E_{i}=\rho_{i}V$ is the total energy of the $i$-th component within the
horizon, $T_{i}$ is the temperature of the $i$-th component, and $S_{i}$ is
the total entropy of the $i$-th component within the horizon. From the first
law (32), we can derive the change rate of the total entropy of matter

$\dot{S}_{m}=\frac{4\pi}{3H^{3}}\sum_{i=1}^{n}\frac{q_{i}}{T_{i}}-4\pi\frac{\dot{H}+H^{2}}{H^{4}}\sum_{i=1}^{n}\frac{\rho_{i}+p_{i}}{T_{i}},$
(33)

where $S_{m}=\sum S_{i}$ is the total entropy of matter within the horizon. On
the other hand, the change rate of the horizon entropy can be obtained from
Eq. (26), given by

$\dot{S}(A)=-\frac{2\pi\dot{H}}{H^{3}}(f_{Q}+2Qf_{QQ}).$ (34)

Finally, the change rate of the total entropy within the horizon can be found

$\dot{S}_{t}=\frac{4\pi}{3H^{3}}\sum_{i=1}^{n}\frac{q_{i}}{T_{i}}-4\pi\frac{\dot{H}+H^{2}}{H^{4}}\sum_{i=1}^{n}\frac{\rho_{i}+p_{i}}{T_{i}}-\frac{2\pi\dot{H}}{H^{3}}(f_{Q}+2Qf_{QQ}),$
(35)

where $S_{t}=S_{m}+S(A)$ is the total entropy within the horizon. The so-
called generalized second law of thermodynamics requires $\dot{S}_{t}\geq 0$,
which is what we are going to discuss.

The composition of the real universe is very complex, which inevitably makes
the analysis of the second law $\dot{S}_{t}\geq 0$ complicated and bloated. In
order to make the analysis of the second law more concise and intuitive, we
study the simplified toy universe in this subsection. Firstly, we assume that
all components have the same temperature, that is, $T_{i}=T_{m}$ for any
$i=1,2,...,n$. Then, by combining the continuity equation (31) and the
background equations (17) and (18), the generalized second law can be
simplified to

$\dot{S}_{t}=\frac{\dot{H}}{H^{4}T_{m}}(\dot{H}+H^{2}-2\pi
HT_{m})(f_{Q}+2Qf_{QQ})\geq 0.$ (36)

Secondly, we also assume that the temperature of matter is equal to the
temperature of the horizon. This assumption has also been widely adopted in
the literature when studying the general second law [65, 59, 66, 61].
Subsequently, the generalized second law can be further simplified to

$\dot{S}_{t}=\frac{\dot{H}^{2}}{2H^{4}T}(f_{Q}+2Qf_{QQ})\geq 0.$ (37)

For a universe with positive temperature, Eq. (37) means

$f_{Q}+2Qf_{QQ}\geq 0.$ (38)

In such a toy universe, the second law of generalized thermodynamics in the
$f(Q)$ universe is equivalent to a very simple inequality (38).

Finally, let’s look at some simple applications of inequality (38). For the
STEGR model which is equivalent to GR, the inequality (38) can be reduced to
$1>0$. This means that the general second law of thermodynamics is always
satisfied in GR. For the $f(Q)$ model that only retains the lowest-order
correction, that is, $f(Q)=Q+\lambda\,Q^{2}$, the inequality (38) is reduced
to $1-36\lambda H^{2}\geq 0$. If the coefficient $\lambda<0$, the generalized
second law always holds. But if the coefficient $\lambda>0$, then the
generalized second law requires $H^{2}<(36\lambda)^{-1}$, that is, the
expansion of the universe cannot be arbitrarily fast. This example shows that
the generalized second law can impose constraints on model parameters or the
evolution of the universe.

###  III.3 $P-V$ Phase Transition

Phase transitions and critical phenomena are fascinating topics in
thermodynamics. After having a preliminary understanding of the thermodynamic
laws, it is natural to ask whether there is a phase transition in the $f(Q)$
universe. In the thermodynamic description of the $f(Q)$ universe, there are
thermodynamic quantities such as pressure, temperature and volume. Therefore,
there may be a phase transition similar to that of a gas-liquid system. In
this subsection, we will analyze the $P-V$ phase transition of the $f(Q)$
universe and calculate its critical exponents.

From Eqs. (22) and (23), we can find the equation of state $P=P(V,T)$ of the
$f(Q)$ universe as

$P=\frac{f}{16\pi}+\frac{f_{Q}}{2\pi R_{A}^{2}}+\frac{3f_{QQ}}{\pi
R_{A}^{4}}+\left(\frac{f_{Q}}{2R_{A}}-\frac{6f_{QQ}}{R_{A}^{3}}\right)T,$ (39)

where the horizon radius $R_{A}$ should be understood as
$R_{A}(V)=(3V/4\pi)^{1/3}$ in this subsection. For the STEGR model, the
equation of state (39) can be simplified to

$P=\frac{1}{8\pi R_{A}^{2}}+\frac{T}{2R_{A}}.$ (40)

Given any non-negative temperature $T$, the pressure $P$ in Eq. (40) is a
monotonically decreasing function of volume $V$. Therefore, there is no $P-V$
phase transition in GR. This result is consistent with the research in Ref.
[67].

Next, we examine the $f(Q)$ model that only retains the lowest-order
correction to GR, that is, the model $f(Q)=Q+\lambda\,Q^{2}$. In this model,
the equation of state (39) is reduced to

$P=\frac{1}{8{\pi}R_{A}^{2}}+\frac{T}{2R_{A}}+\frac{9\lambda}{4{\pi}R_{A}^{4}}-\frac{18{\lambda}T}{R_{A}^{3}}.$
(41)

The necessary conditions for the $P-V$ phase transition are

$\left(\frac{\partial P}{\partial
V}\right)_{T}=\left(\frac{\partial^{2}P}{\partial V^{2}}\right)_{T}=0.$ (42)

For the case of $\lambda>0$, Eq. (42) has the following physically reasonable
solution

$T_{c}=\frac{\sqrt{3+2\sqrt{3}}}{36\pi\sqrt{\lambda}}\leavevmode\nobreak\
,\leavevmode\nobreak\ R_{c}=6\sqrt{(2\sqrt{3}-3)\lambda}\leavevmode\nobreak\
,\leavevmode\nobreak\ P_{c}=\frac{15+8\sqrt{3}}{5184\pi\lambda},$ (43)

where $T_{c}$ is the critical temperature, $R_{c}$ is the horizon radius at
the critical point, and $P_{c}$ is the pressure at the critical point. The
universal dimensionless quantity for this critical point is

$\frac{P_{c}R_{c}}{T_{c}}=1+\frac{\sqrt{3}}{6},$ (44)

which is a constant and independent of the model parameter $\lambda$. The
above results show that there is a $P-V$ phase transition in the model
$f(Q)=Q+\lambda\,Q^{2}$ with $\lambda>0$. This can be seen more intuitively
from the $P-R_{A}$ phase diagram.

Figure 1:  Isotherms in the $P-R_{A}$ phase diagram

Figure 1 shows the isotherms near the critical temperature, where warmer
colors mean higher temperatures and cooler colors mean lower temperatures. It
can be seen that the shape of the isotherms is somewhat similar to those of
the van der Waals gas-liquid system, some black hole systems [68, 69, 70, 71,
72], the special Horndeski universe [23] and the holographic dark energy model
[73]. However, there is a significant difference from the van der Waals
system, that is, the coexistence phase appears above the critical temperature
$T>T_{c}$ in the $f(Q)$ universe, whereas the coexistence phase appears below
the critical temperature $T<T_{c}$ in the van der Waals system. In addition,
the critical temperature $T_{c}\rightarrow+\infty$ if $\lambda\rightarrow
0^{+}$, which is consistent with the previous conclusion that there is no
phase transition in the STEGR model.

Similar to the gas-liquid phase transition, we can define the critical
exponents $\\{\alpha,\beta,\gamma,\delta\\}$ near the critical point as
follows:

$\displaystyle
C_{V}=T\left(\frac{{\partial}S}{{\partial}T}\right)_{V}\leavevmode\nobreak\ $
$\displaystyle\propto$ $\displaystyle\leavevmode\nobreak\ |t|^{-\alpha},$ (45)
$\displaystyle\eta=\frac{V_{l}-V_{s}}{V_{c}}\leavevmode\nobreak\ $
$\displaystyle\propto$ $\displaystyle\leavevmode\nobreak\ |t|^{\beta},$ (46)
$\displaystyle\kappa_{T}=-\frac{1}{V}\left(\frac{\partial V}{\partial
P}\right)_{T}\leavevmode\nobreak\ $ $\displaystyle\propto$
$\displaystyle\leavevmode\nobreak\ |t|^{-\gamma},$ (47) $\displaystyle
p|_{t=0}\leavevmode\nobreak\ $ $\displaystyle\propto$
$\displaystyle\leavevmode\nobreak\ |v|^{\delta},$ (48)

where

$t=\frac{T-T_{c}}{T_{c}}\leavevmode\nobreak\ ,\leavevmode\nobreak\
p=\frac{P-P_{c}}{P_{c}}\leavevmode\nobreak\ ,\leavevmode\nobreak\
v=\frac{V-V_{c}}{V_{c}},$ (49)

are the reduced temperature, pressure, and volume, respectively, while $V_{l}$
and $V_{s}$ are two different volumes with equal pressure and Gibbs free
energy. In the language of gas-liquid systems, $C_{V}$ is the isovolumetric
heat capacity, $\kappa_{T}$ is the isothermal compressibility, and $\eta$ can
be regarded as the order parameter. In the following, we will calculate the
four critical exponents and examine their scaling laws.

Since the entropy is a function only of volume and has nothing to do with
temperature, the isovolumetric heat capacity is zero, i.e. $C_{V}=0$, which
means the first critical exponent $\alpha=0$. In order to obtain the other
three critical ones, we expand the reduced pressure $p$ near the critical
point

$p=-\frac{8}{11}(2\sqrt{3}-1)t+\frac{8}{33}(3+5\sqrt{3})vt-\frac{4}{297}(5+\sqrt{3})v^{3}+...$
(50)

where ”$...$” means unimportant higher-order terms. Note that after we use the
reduced thermodynamic variables, the equation of state becomes independent of
the model parameter $\lambda$. From Eq. (50), one can easily find the
isothermal compressibility and the pressure near the critical point as

$\kappa_{T}\propto\left(\frac{{\partial}p}{{\partial}v}\right)^{-1}{\propto}\leavevmode\nobreak\
t^{-1}\leavevmode\nobreak\ ,\leavevmode\nobreak\ p|_{t=0}\propto v^{3},$ (51)

which gives the critical exponents $\gamma=1$ and $\delta=3$. To obtain the
last critical exponent, we must solve for $v_{s}$ and $v_{l}$, which are
different reduced volumes with the same pressure and the same Gibbs free
energy $G=U+PV-TS$. Clearly, $v_{l}$ and $v_{s}$ satisfy the following two
equations

$p(v_{s},t)=p(v_{l},t)\leavevmode\nobreak\ \leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \int
vdp=\int_{v_{s}}^{v_{l}}v\left(\frac{\partial p}{\partial
v}\right)\mathrm{d}v=0.$ (52)

From Eqs. (50) and (52), one can get a nontrivial solution

$v_{l}=3\sqrt{2\sqrt{3}t}\leavevmode\nobreak\ ,\leavevmode\nobreak\
v_{s}=-3\sqrt{2\sqrt{3}t}.$ (53)

Therefore, $\eta=v_{l}-v_{s}\propto t^{1/2}$, and the last critical exponent
$\beta=1/2$.

In summary, the four critical exponents are

$\alpha=0\leavevmode\nobreak\ ,\leavevmode\nobreak\
\beta=\frac{1}{2}\leavevmode\nobreak\ ,\leavevmode\nobreak\
\gamma=1\leavevmode\nobreak\ ,\leavevmode\nobreak\ \delta=3.$ (54)

Somewhat surprisingly, these critical exponents are exactly the same as those
in the mean field theory. So of course they obey the usual scaling laws

$\displaystyle\alpha+2\beta+\gamma=2,\leavevmode\nobreak\
\alpha+\beta(1+\delta)=2,$
$\displaystyle\gamma(1+\delta)=(2-\alpha)(\delta-1),\leavevmode\nobreak\
\gamma=\beta(\delta-1).$ (55)

In subsection III.1, we showed that the lowest-order correction to GR by the
$f(Q)$ model leads to a logarithmic correction to the entropy, which is
consistent with the lowest-order correction in loop quantum gravity. In this
subsection, we show that the lowest-order correction to GR by the $f(Q)$ model
can yield the $P-V$ phase transition, and the critical exponents are exactly
the same as those in the mean field theory. These unexpected coincidences seem
to hint that the correction brought by the $f(Q)$ model may emerge from some
quantum mechanism, and its microscopic statistics satisfy the mean field
theory in the lowest-order approximation. Of course, we have no evidence to
support this idea so far, and it could just be a coincidence.

##  IV Thermodynamic implication of nontrivial connections

In section III, we analyzed the thermodynamics of the $f(Q)$ universe with
trivial connection from the perspectives of the first law, generalized second
law and $P-V$ phase transition. As mentioned in section II, there are flat
universes with nontrivial connection in the $f(Q)$ model. In this section, we
will explore the thermodynamic implications of these nontrivial connections in
the $f(Q)$ flat universe.

In fact, it is unprecedented to discuss how different connections affect the
thermodynamics of the universe in the same gravity model. Because in most
gravity models, such as $f(R)$ and $f(T)$ models, the connection that
satisfies cosmological symmetry is unique. Without different connections,
naturally there would be no discussion of the thermodynamic effects of
different connections. As a result, our exploration in this section is
somewhat pioneering.

In the $f(Q)$ flat universe, the connection that obeys cosmological symmetry
has three branch solutions, as shown in section II. For the branch 1 solution,
one can find that the background equations in the branch 1 universe are
exactly the same as those in the universe with trivial connection. Therefore,
the branch 1 universe should also be classified as a universe with trivial
connection. To analyze the thermodynamic effects of nontrivial connections, we
should focus on branch 2 and branch 3 universes.

For the branch 2 universe, following the process in subsection III.1, we
obtain the energy and work density as

$\displaystyle
E=\frac{R_{A}f_{Q}}{2}+\frac{R_{A}^{3}}{12}(f-Qf_{Q})+\frac{R_{A}^{3}}{4}\gamma\dot{f}_{Q},$
(56) $\displaystyle
W=\frac{1}{8\pi}\Big{[}f_{Q}(\dot{H}+3H^{2})+\frac{1}{2}(f-Qf_{Q})+H\dot{f}_{Q}\Big{]}.$
(57)

After identifying the internal energy $U=-E$ and the thermodynamic pressure
$P=W$, the first law of thermodynamic (25) can be established as long as

$\displaystyle\mathrm{d}S=\bigg{\\{}2\pi R_{A}f_{Q}+3\pi
f_{QQ}\Big{[}18\gamma-\frac{8}{R_{A}}-9R_{A}\gamma^{2}+(3R_{A}\gamma-2)(R_{A}\ddot{\gamma}+3\dot{\gamma})\frac{R_{A}}{\dot{R}_{A}}\Big{]}\bigg{\\}}\,\mathrm{d}R_{A},\leavevmode\nobreak\
\leavevmode\nobreak\ $ (58)

where $Q=-6R_{A}^{-2}+9\gamma R_{A}^{-1}+3\dot{\gamma}$ and $\gamma$ satisfies
the following equation due to Eq. (16)

$\gamma(\ddot{f}_{Q}+3H\dot{f}_{Q})=0.$ (59)

At first glance, the entropy in Eq. (58) is not-integrable unless $\gamma=0$
or $f_{Q}=\text{const}$. To see this more clearly, let us again consider the
simple model $f(Q)=Q+\lambda Q^{2}$. In this model, the differential of
entropy (58) can be simplified to

$\mathrm{d}S=\mathrm{d}\left(\frac{A}{4}-36\pi\lambda\ln\frac{A}{A_{0}}\right)+\mathrm{d}\hat{S},$
(60)

where

$\mathrm{d}\hat{S}=6\pi\lambda\left\\{24\gamma\dot{R}_{A}+R_{A}\left[(2\dot{\gamma}-9\gamma^{2})\dot{R}_{A}+(3R_{A}\gamma-2)(R_{A}\ddot{\gamma}+3\dot{\gamma})\right]\right\\}\mathrm{d}t.$
(61)

If $\lambda=0$ (equivalent to GR) or $\gamma=0$ (trivial connection), then
$\mathrm{d}\hat{S}=0$, so the entropy $S$ in Eq. (60) can be integrated as a
function of the horizon area. If $\lambda\neq 0$ and $\gamma\neq 0$, then
$\gamma$ should satisfy the following equation due to Eq. (59)

$\dddot{\gamma}+\frac{6}{R_{A}}\ddot{\gamma}+\frac{3(3-2\dot{R}_{A})}{R_{A}^{2}}\dot{\gamma}-\frac{3(3\dot{R}_{A}-2\dot{R}_{A}^{2}+R_{A}\ddot{R}_{A})}{R_{A}^{3}}\gamma+\frac{4(3\dot{R}_{A}-3\dot{R}_{A}^{2}+R_{A}\ddot{R}_{A})}{R_{A}^{4}}=0.\leavevmode\nobreak\
\leavevmode\nobreak\ $ (62)

Since the dependence of $\gamma$ on $\\{R_{A},\dot{R}_{A},\ddot{R}_{A}\\}$ is
realized through a third-order differential equation, it is impossible to
eliminate $\gamma$ in Eq. (61) by Eq. (62). This means that $\hat{S}$ cannot
be integrated as a function of the horizon radius $R_{A}$. For the same
reason, the entropy in Eq. (60) cannot be integrated into a function only
about $R_{A}$ and $\dot{R}_{A}$, that is, it cannot be expressed as a function
$S=S(V,T)$ only about volume and temperature. A more rigorous proof can be
found in Appendix A. From Eq. (58) we can expect the same thing to happen in
other $f(Q)$ models.

The situation of the branch 3 universe is almost the same as that of the
branch 2 universe. In the branch 3 universe, the energy and work density are

$\displaystyle
E=\frac{R_{A}f_{Q}}{2}+\frac{R_{A}^{3}}{12}(f-Qf_{Q})-\frac{R_{A}^{3}}{4}\gamma\dot{f}_{Q},$
(63) $\displaystyle
W=\frac{1}{8\pi}\Big{[}f_{Q}(\dot{H}+3H^{2})+\frac{1}{2}(f-Qf_{Q})+(H-\gamma)\dot{f}_{Q}\Big{]}.$
(64)

After identifying the internal energy $U=-E$ and the thermodynamic pressure
$P=W$, the first law of thermodynamic (25) gives

$\displaystyle\mathrm{d}S=\bigg{\\{}2\pi R_{A}f_{Q}+3\pi
f_{QQ}\Big{[}-\frac{8}{R_{A}}+2\gamma+3R_{A}\gamma^{2}-(2+R_{A}\gamma)(R_{A}\ddot{\gamma}+3\dot{\gamma})\frac{R_{A}}{\dot{R}_{A}}\Big{]}\bigg{\\}}\,\mathrm{d}R_{A},\leavevmode\nobreak\
\leavevmode\nobreak\ $ (65)

where $Q=-6R_{A}^{-2}+9\gamma R_{A}^{-1}+3\dot{\gamma}$ and
$\gamma\ddot{f}_{Q}+(2\dot{\gamma}+5H\gamma)\dot{f}_{Q}=0$. Once again, the
entropy in Eq. (65) appears to be non-integrable unless $\gamma=0$ or
$f_{Q}=\text{const}$. This can be seen more clearly in a simple example of
$f(Q)=Q+\lambda Q^{2}$. In this model, the differential of entropy (58) can be
simplified to

$\mathrm{d}S=\mathrm{d}\left(\frac{A}{4}-36\pi\lambda\ln\frac{A}{A_{0}}\right)+\mathrm{d}\hat{S},$
(66)

with

$\mathrm{d}\hat{S}=6\pi\lambda\left\\{8\gamma\dot{R}_{A}+R_{A}\left[(2\dot{\gamma}+3\gamma^{2})\dot{R}_{A}-(2+R_{A}\gamma)(R_{A}\ddot{\gamma}+3\dot{\gamma})\right]\right\\}\mathrm{d}t,$
(67)

where $\gamma$ satisfies the equation

$\displaystyle\gamma\dddot{\gamma}+2\dot{\gamma}\ddot{\gamma}+\frac{2}{R_{A}}(4\gamma\ddot{\gamma}+3\dot{\gamma}^{2})+\frac{3(5-4\dot{R}_{A})}{R_{A}^{2}}\gamma\dot{\gamma}$
$\displaystyle\quad\quad-\frac{3(5\dot{R}_{A}-2\dot{R}_{A}^{2}+R_{A}\ddot{R}_{A})}{R_{A}^{3}}\gamma^{2}+\frac{8\dot{R}_{A}}{R_{A}^{3}}\dot{\gamma}+\frac{4(5\dot{R}_{A}-3\dot{R}_{A}^{2}+R_{A}\ddot{R}_{A})}{R_{A}^{4}}\gamma=0.$
(68)

For the same reasons as for the branch 2 universe, the entropy $S$ in Eq. (66)
cannot be integrated as a function of $R_{A}$ and $\dot{R}_{A}$ unless
$\gamma=0$ or $\lambda=0$.

In summary, in the $f(Q)$ flat universe with nontrivial connection, the first
law of thermodynamic gives

$\displaystyle\mathrm{d}S=\mathrm{d}S(A)+\mathrm{d}\hat{S}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{with}\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{d}\hat{S}=\mathcal{F}(\gamma,\dot{\gamma},\ddot{\gamma},R_{A},\dot{R}_{A})\mathrm{d}t,$
(69)

where $\gamma$ depends on $\\{R_{A},\dot{R}_{A},\ddot{R}_{A}\\}$ through a
third-order differential equation. Consequently, the entropy cannot be
integrated as a function of the horizon radius and the horizon temperature in
general. The change of entropy depends on more specific process of cosmic
evolutions. Following the discussions in Refs. [22, 13, 14, 74], the new term
$\mathrm{d}\hat{S}$ can be interpreted as the entropy production term. This
reveals that the horizon thermodynamics in the $f(Q)$ universe with nontrivial
connection is non-equilibrium thermodynamics 111 Just as the $f(R)$ and $f(T)$
universes can be described by equilibrium thermodynamic [26, 27, 30], the
$f(Q)$ universe with nontrivial connection can also be described by
equilibrium thermodynamics. However, the realization of the equilibrium
description can be independent of the gravity model and spacetime connection
(see Appendix B), so we do not adopt this perspective in this paper. . This is
a significant feature of nontrivial connections from the perspective of
thermodynamics.

##  V Conclusion

In this paper, we have analyzed the thermodynamics of apparent horizon in the
$f(Q)$ universe with trivial and nontrivial connections. For the case of the
trivial connection, first, we have studied its first law and showed that the
first law in thermal equilibrium holds for any $f(Q)$ models. We have also
obtained the area law of entropy for the most general $f(Q)$ model and showed
that the lowest-order and second-lowest-order corrections to the entropy have
the same form as those in loop quantum gravity. Then, we briefly investigated
the generalized second law, which asserts that the entropy of the horizon plus
the entropy of matter within the horizon does not decrease. After assuming
that all matter is in thermal equilibrium with the horizon, we derived a
simple inequality (38) from the generalized second law that holds for any
$f(Q)$ models. Finally, we have studied the $P-V$ phase transition by analogy
with the gas-liquid system. We have shown that even the $f(Q)$ model that
retains only the lowest-order correction to GR can lead to the $P-V$ phase
transition that does not exist in GR. We have calculated the critical
exponents and found that they are exactly the same as those in mean field
theory. For the case of nontrivial connections, we have examined the
thermodynamic effect of nontrivial connections from the perspective of the
first law. We have proved that the change in entropy cannot be determined
entirely by the change in volume and temperature. This means that the horizon
thermodynamics of the $f(Q)$ universe with nontrivial connection is non-
equilibrium. This may be the first exploration of how different connections
affect the thermodynamics of the universe in the same gravity model.

### Acknowledgements

This work is supported in part by the National Key Research and Development
Program of China under Grant No. 2020YFC2201501 and the National Natural
Science Foundation of China (NSFC) under Grant No. 12347103 and 12205063.

##  Appendix A The condition for $S=S(V,T)$

In this appendix, we will discuss the condition for the integrability of
entropy, and further prove that the entropies in Eqs. (61) and (67) in the
main text are not integrable in the thermodynamic sense.

In the thermodynamics of the horizon, sometimes entropy can be integrated as a
function of the horizon radius. This is a remarkable thermodynamic property
that is not always possible. Once this fails, some might conclude that entropy
is non-integrable. But we think it is too hasty to draw such a conclusion.
Because it is also possible that entropy can be integrated as
$S=S(R_{A},\dot{R}_{A})$, so that entropy is a function of volume
$V=\frac{4}{3}\pi R_{A}^{3}$ and temperature (22). This means that the change
of entropy only depends on the initial and final states of the system, so such
entropy is integrable.

If entropy can be integrated as a function of volume and temperature, that is,
$S=S(R_{A},\dot{R}_{A})$, then the differential of entropy is

$\mathrm{d}S=\frac{\partial S}{\partial R_{A}}\mathrm{d}R_{A}+\frac{\partial
S}{\partial\dot{R}_{A}}\mathrm{d}\dot{R}_{A}=\left(\frac{\partial S}{\partial
R_{A}}\dot{R}_{A}+\frac{\partial
S}{\partial\dot{R}_{A}}\ddot{R}_{A}\right)\mathrm{d}t.$ (70)

As as result, the differential of entropy must have the form

$\mathrm{d}S=\left[f_{1}(R_{A},\dot{R}_{A})+f_{2}(R_{A},\dot{R}_{A})\ddot{R}_{A}\right]\mathrm{d}t.$
(71)

Since the order of partial derivatives of smooth functions is commutative,
that is

$\frac{\partial^{2}S}{\partial
R_{A}\partial\dot{R}_{A}}=\frac{\partial^{2}S}{\partial\dot{R}_{A}\partial{R}_{A}},$
(72)

$f_{1}(R_{A},\dot{R}_{A})$ and $f_{2}(R_{A},\dot{R}_{A})$ must satisfy the
following condition:

$\dot{R}_{A}\frac{\partial
f_{1}}{\partial\dot{R}_{A}}-f_{1}=\dot{R}_{A}^{2}\frac{\partial
f_{2}}{\partial{R}_{A}}.$ (73)

Therefore, if entropy is integrable in the thermodynamic sense, the
differential of entropy should be consistent with the form of Eq. (71) and
satisfy the condition (73).

However, the differential of entropy in Eq. (61) cannot even be consistent
with the form of Eq. (71), because the dependence of $\gamma$ on
$\\{R_{A},\dot{R}_{A},\ddot{R}_{A}\\}$ is implemented through a differential
equation rather than a algebraic equation. Consequently, the differential of
entropy in Eq. (61) is not integrable at least in the thermodynamic sense. For
the same reason, the differential entropy in Eq. (67) is also not integrable
in the thermodynamic sense.

##  Appendix B Thermal equilibrium description of horizon thermodynamics

In this appendix, we will show that the horizon in a flat universe can always
be described by equilibrium thermodynamics.

In a flat universe, the background equations can be written as

$\displaystyle 3H^{2}=8\pi(\rho_{m}+\rho_{MG}),$ (74) $\displaystyle
2\dot{H}+3H^{2}=-8\pi(p_{m}+p_{MG}),$ (75)

where $\rho_{MG}$ and $p_{MG}$ are the effective energy density and effective
pressure resulting from the modified gravity. With the internal energy
$U=-(\rho_{m}+\rho_{MG})V$ and the thermodynamic pressure
$P=\frac{1}{2}(\rho_{m}+\rho_{MG}-p_{m}-p_{MG})$, the first law of
thermodynamics (25) holds, where the entropy $S$ is always equal to $A/4$,
regardless of the specific gravity model. Since the above derivation depends
neither on a specific gravity model nor on spacetime connection, it is
applicable to any flat universes with cosmological symmetry in any gravity
models. The details of the gravity model are contained in $\rho_{MG}$ and
$p_{MG}$.

In this perspective, the characteristics of the $f(Q)$ model and the
distinction between trivial and nontrivial connections are difficult to
capture. This is why we did not adopt the equilibrium description in this
paper.

## References

  * [1] J. M. Bardeen, B. Carter and S. W. Hawking, Commun. Math. Phys. 31, 161-170 (1973) doi:10.1007/BF01645742
  * [2] J. D. Bekenstein, Phys. Rev. D 7, 2333-2346 (1973) doi:10.1103/PhysRevD.7.2333
  * [3] S. W. Hawking, Commun. Math. Phys. 43, 199-220 (1975) [erratum: Commun. Math. Phys. 46, 206 (1976)] doi:10.1007/BF02345020
  * [4] T. Jacobson, Phys. Rev. Lett. 75, 1260-1263 (1995) doi:10.1103/PhysRevLett.75.1260 [arXiv:gr-qc/9504004 [gr-qc]].
  * [5] T. Jacobson, Phys. Rev. Lett. 116, no.20, 201101 (2016) doi:10.1103/PhysRevLett.116.201101 [arXiv:1505.04753 [gr-qc]].
  * [6] T. Padmanabhan, Rept. Prog. Phys. 73, 046901 (2010) doi:10.1088/0034-4885/73/4/046901 [arXiv:0911.5004 [gr-qc]].
  * [7] R. M. Wald, Living Rev. Rel. 4, 6 (2001) doi:10.12942/lrr-2001-6 [arXiv:gr-qc/9912119 [gr-qc]].
  * [8] D. Kubiznak, R. B. Mann and M. Teo, Class. Quant. Grav. 34, no.6, 063001 (2017) doi:10.1088/1361-6382/aa5c69 [arXiv:1608.06147 [hep-th]].
  * [9] S. W. Hawking and D. N. Page, Commun. Math. Phys. 87, 577 (1983) doi:10.1007/BF01208266
  * [10] N. Altamirano, D. Kubiznak, R. B. Mann and Z. Sherkatghanad, Galaxies 2, 89-159 (2014) doi:10.3390/galaxies2010089 [arXiv:1401.2586 [hep-th]].
  * [11] R. G. Cai, Phys. Rev. D 65, 084014 (2002) doi:10.1103/PhysRevD.65.084014 [arXiv:hep-th/0109133 [hep-th]].
  * [12] R. G. Cai, Phys. Lett. B 582, 237-242 (2004) doi:10.1016/j.physletb.2004.01.015 [arXiv:hep-th/0311240 [hep-th]].
  * [13] M. Akbar and R. G. Cai, Phys. Lett. B 648, 243-248 (2007) doi:10.1016/j.physletb.2007.03.005 [arXiv:gr-qc/0612089 [gr-qc]].
  * [14] R. X. Miao, M. Li and Y. G. Miao, JCAP 11, 033 (2011) doi:10.1088/1475-7516/2011/11/033 [arXiv:1107.0515 [hep-th]].
  * [15] T. Padmanabhan, Phys. Rept. 406, 49-125 (2005) doi:10.1016/j.physrep.2004.10.003 [arXiv:gr-qc/0311036 [gr-qc]].
  * [16] R. G. Cai and S. P. Kim, JHEP 02, 050 (2005) doi:10.1088/1126-6708/2005/02/050 [arXiv:hep-th/0501055 [hep-th]].
  * [17] G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2738-2751 (1977) doi:10.1103/PhysRevD.15.2738
  * [18] R. G. Cai, L. M. Cao and Y. P. Hu, Class. Quant. Grav. 26, 155018 (2009) doi:10.1088/0264-9381/26/15/155018 [arXiv:0809.1554 [hep-th]].
  * [19] R. Di Criscienzo, S. A. Hayward, M. Nadalini, L. Vanzo and S. Zerbini, Class. Quant. Grav. 27, 015006 (2010) doi:10.1088/0264-9381/27/1/015006 [arXiv:0906.1725 [gr-qc]].
  * [20] T. Zhu, J. R. Ren and D. Singleton, Int. J. Mod. Phys. D 19, 159-169 (2010) doi:10.1142/S0218271810016336 [arXiv:0902.2542 [hep-th]].
  * [21] M. Akbar and R. G. Cai, Phys. Rev. D 75, 084003 (2007) doi:10.1103/PhysRevD.75.084003 [arXiv:hep-th/0609128 [hep-th]].
  * [22] R. G. Cai and L. M. Cao, Phys. Rev. D 75, 064008 (2007) doi:10.1103/PhysRevD.75.064008 [arXiv:gr-qc/0611071 [gr-qc]].
  * [23] S. B. Kong, H. Abdusattar, Y. Yin, H. Zhang and Y. P. Hu, Eur. Phys. J. C 82, no.11, 1047 (2022) doi:10.1140/epjc/s10052-022-10976-9 [arXiv:2108.09411 [gr-qc]].
  * [24] M. Akbar and R. G. Cai, Phys. Lett. B 635, 7-10 (2006) doi:10.1016/j.physletb.2006.02.035 [arXiv:hep-th/0602156 [hep-th]].
  * [25] Y. Gong and A. Wang, Phys. Rev. Lett. 99, 211301 (2007) doi:10.1103/PhysRevLett.99.211301 [arXiv:0704.0793 [hep-th]].
  * [26] K. Bamba, C. Q. Geng and S. Tsujikawa, Phys. Lett. B 688, 101-109 (2010) doi:10.1016/j.physletb.2010.03.070 [arXiv:0909.2159 [gr-qc]].
  * [27] K. Bamba and C. Q. Geng, JCAP 11, 008 (2011) doi:10.1088/1475-7516/2011/11/008 [arXiv:1109.1694 [gr-qc]].
  * [28] K. Bamba and C. Q. Geng, Phys. Lett. B 679, 282-287 (2009) doi:10.1016/j.physletb.2009.07.039 [arXiv:0901.1509 [hep-th]].
  * [29] K. Bamba, C. Q. Geng, S. Nojiri and S. D. Odintsov, EPL 89, no.5, 50003 (2010) doi:10.1209/0295-5075/89/50003 [arXiv:0909.4397 [hep-th]].
  * [30] K. Bamba, C. Q. Geng and S. Tsujikawa, Int. J. Mod. Phys. D 20, 1363-1371 (2011) doi:10.1142/S0218271811019542 [arXiv:1101.3628 [gr-qc]].
  * [31] K. Bamba and C. Q. Geng, JCAP 06, 014 (2010) doi:10.1088/1475-7516/2010/06/014 [arXiv:1005.5234 [gr-qc]].
  * [32] J. M. Nester and H. J. Yo, Chin. J. Phys. 37, 113 (1999) [arXiv:gr-qc/9809049 [gr-qc]].
  * [33] J. Beltrán Jiménez, L. Heisenberg and T. S. Koivisto, Universe 5, no.7, 173 (2019) doi:10.3390/universe5070173 [arXiv:1903.06830 [hep-th]].
  * [34] S. Capozziello, V. De Falco and C. Ferrara, Eur. Phys. J. C 82, no.10, 865 (2022) doi:10.1140/epjc/s10052-022-10823-x [arXiv:2208.03011 [gr-qc]].
  * [35] J. Beltrán Jiménez, L. Heisenberg and T. Koivisto, Phys. Rev. D 98, no.4, 044048 (2018) doi:10.1103/PhysRevD.98.044048 [arXiv:1710.03116 [gr-qc]].
  * [36] L. Heisenberg, Phys. Rept. 1066, 1-78 (2024) doi:10.1016/j.physrep.2024.02.001 [arXiv:2309.15958 [gr-qc]].
  * [37] J. Beltrán Jiménez, L. Heisenberg, T. S. Koivisto and S. Pekar, Phys. Rev. D 101, no.10, 103507 (2020) doi:10.1103/PhysRevD.101.103507 [arXiv:1906.10027 [gr-qc]].
  * [38] L. Atayde and N. Frusciante, Phys. Rev. D 104, no.6, 064052 (2021) doi:10.1103/PhysRevD.104.064052 [arXiv:2108.10832 [astro-ph.CO]].
  * [39] A. Oliveros and M. A. Acero, Int. J. Mod. Phys. D 33, no.01, 2450004 (2024) doi:10.1142/S0218271824500044 [arXiv:2311.01857 [astro-ph.CO]].
  * [40] W. Khyllep, J. Dutta, E. N. Saridakis and K. Yesmakhanova, Phys. Rev. D 107, no.4, 044022 (2023) doi:10.1103/PhysRevD.107.044022 [arXiv:2207.02610 [gr-qc]].
  * [41] D. Zhao, Eur. Phys. J. C 82, no.4, 303 (2022) doi:10.1140/epjc/s10052-022-10266-4 [arXiv:2104.02483 [gr-qc]].
  * [42] K. Hu, T. Paul and T. Qiu, Sci. China Phys. Mech. Astron. 67, no.2, 220413 (2024) doi:10.1007/s11433-023-2275-0 [arXiv:2308.00647 [hep-th]].
  * [43] M. Hohmann, Phys. Rev. D 104, no.12, 124077 (2021) doi:10.1103/PhysRevD.104.124077 [arXiv:2109.01525 [gr-qc]].
  * [44] L. Heisenberg, M. Hohmann and S. Kuhn, Eur. Phys. J. C 83, no.4, 315 (2023) doi:10.1140/epjc/s10052-023-11462-6 [arXiv:2212.14324 [gr-qc]].
  * [45] M. Hohmann, Int. J. Geom. Meth. Mod. Phys. 18, no.supp01, 2140005 (2021) doi:10.1142/S0219887821400053 [arXiv:2008.12186 [gr-qc]].
  * [46] D. A. Gomes, J. Beltrán Jiménez and T. S. Koivisto, JCAP 12, 010 (2023) doi:10.1088/1475-7516/2023/12/010 [arXiv:2309.08554 [gr-qc]].
  * [47] N. Dimakis, A. Paliathanasis, M. Roumeliotis and T. Christodoulakis, Phys. Rev. D 106, no.4, 043509 (2022) doi:10.1103/PhysRevD.106.043509 [arXiv:2205.04680 [gr-qc]].
  * [48] J. Shi, Eur. Phys. J. C 83, no.10, 951 (2023) doi:10.1140/epjc/s10052-023-12139-w [arXiv:2307.08103 [gr-qc]].
  * [49] D. A. Gomes, J. Beltrán Jiménez, A. J. Cano and T. S. Koivisto, Phys. Rev. Lett. 132, no.14, 141401 (2024) doi:10.1103/PhysRevLett.132.141401 [arXiv:2311.04201 [gr-qc]].
  * [50] H. Rao, C. Liu and C. Q. Geng, Phys. Lett. B 850, 138497 (2024) doi:10.1016/j.physletb.2024.138497 [arXiv:2311.06600 [gr-qc]].
  * [51] D. Bak and S. J. Rey, Class. Quant. Grav. 17, L83 (2000) doi:10.1088/0264-9381/17/15/101 [arXiv:hep-th/9902173 [hep-th]].
  * [52] S. A. Hayward, S. Mukohyama and M. C. Ashworth, Phys. Lett. A 256, 347-350 (1999) doi:10.1016/S0375-9601(99)00225-X [arXiv:gr-qc/9810006 [gr-qc]].
  * [53] S. A. Hayward, Class. Quant. Grav. 15, 3147-3162 (1998) doi:10.1088/0264-9381/15/10/017 [arXiv:gr-qc/9710089 [gr-qc]].
  * [54] K. A. Meissner, Class. Quant. Grav. 21, 5245-5252 (2004) doi:10.1088/0264-9381/21/22/015 [arXiv:gr-qc/0407052 [gr-qc]].
  * [55] A. Ghosh and P. Mitra, Phys. Rev. D 71, 027502 (2005) doi:10.1103/PhysRevD.71.027502 [arXiv:gr-qc/0401070 [gr-qc]].
  * [56] A. Chatterjee and P. Majumdar, Phys. Rev. Lett. 92, 141301 (2004) doi:10.1103/PhysRevLett.92.141301 [arXiv:gr-qc/0309026 [gr-qc]].
  * [57] R. Banerjee and S. K. Modak, JHEP 05, 063 (2009) doi:10.1088/1126-6708/2009/05/063 [arXiv:0903.3321 [hep-th]].
  * [58] S. K. Modak, Phys. Lett. B 671, 167-173 (2009) doi:10.1016/j.physletb.2008.11.043 [arXiv:0807.0959 [hep-th]].
  * [59] H. Mohseni Sadjadi and M. Jamil, EPL 92, no.6, 69001 (2010) doi:10.1209/0295-5075/92/69001 [arXiv:1002.3588 [gr-qc]].
  * [60] S. Das, S. Shankaranarayanan and S. Sur, Phys. Rev. D 77, 064013 (2008) doi:10.1103/PhysRevD.77.064013 [arXiv:0705.2070 [gr-qc]].
  * [61] N. Radicella and D. Pavon, Phys. Lett. B 691, 121-126 (2010) doi:10.1016/j.physletb.2010.06.019 [arXiv:1006.3745 [gr-qc]].
  * [62] A. Sheykhi and M. Jamil, Gen. Rel. Grav. 43, 2661-2672 (2011) doi:10.1007/s10714-011-1190-x [arXiv:1011.0134 [physics.gen-ph]].
  * [63] P. C. W. Davies, Class. Quant. Grav. 4, L225 (1987) doi:10.1088/0264-9381/4/6/006
  * [64] T. Clifton and J. D. Barrow, Phys. Rev. D 75, 043515 (2007) doi:10.1103/PhysRevD.75.043515 [arXiv:gr-qc/0701070 [gr-qc]].
  * [65] H. Mohseni Sadjadi, Phys. Rev. D 73, 063525 (2006) doi:10.1103/PhysRevD.73.063525 [arXiv:gr-qc/0512140 [gr-qc]].
  * [66] U. Debnath, S. Chattopadhyay, I. Hussain, M. Jamil and R. Myrzakulov, Eur. Phys. J. C 72, 1875 (2012) doi:10.1140/epjc/s10052-012-1875-7 [arXiv:1111.3858 [gr-qc]].
  * [67] H. Abdusattar, S. B. Kong, W. L. You, H. Zhang and Y. P. Hu, JHEP 12, 168 (2022) doi:10.1007/JHEP12(2022)168 [arXiv:2108.09407 [gr-qc]].
  * [68] C. H. Nam, Eur. Phys. J. C 78, no.12, 1016 (2018) doi:10.1140/epjc/s10052-018-6498-1
  * [69] J. Xu, L. M. Cao and Y. P. Hu, Phys. Rev. D 91, no.12, 124033 (2015) doi:10.1103/PhysRevD.91.124033 [arXiv:1506.03578 [gr-qc]].
  * [70] S. H. Hendi, R. B. Mann, S. Panahiyan and B. Eslam Panah, Phys. Rev. D 95, no.2, 021501 (2017) doi:10.1103/PhysRevD.95.021501 [arXiv:1702.00432 [gr-qc]].
  * [71] Y. P. Hu, H. A. Zeng, Z. M. Jiang and H. Zhang, Phys. Rev. D 100, no.8, 084004 (2019) doi:10.1103/PhysRevD.100.084004 [arXiv:1812.09938 [gr-qc]].
  * [72] S. W. Wei and Y. X. Liu, Phys. Rev. D 101, no.10, 104018 (2020) doi:10.1103/PhysRevD.101.104018 [arXiv:2003.14275 [gr-qc]].
  * [73] M. Cruz, S. Lepe and J. Saavedra, [arXiv:2312.14257 [gr-qc]].
  * [74] C. Eling, R. Guedens and T. Jacobson, Phys. Rev. Lett. 96, 121301 (2006) doi:10.1103/PhysRevLett.96.121301 [arXiv:gr-qc/0602001 [gr-qc]].