# Generalizability of predictive models for intensive care unit patients

Alistair E. W. Johnson, Tom J. Pollard
Institute of Medical Engineering & Science
Massachussetts Institute of Technology
Cambridge, MA, USA
{aewj<EMAIL_ADDRESS>
Tristan Naumann
Microsoft Research
Redmond, WA, USA
<EMAIL_ADDRESS>

###### Abstract

A large volume of research has considered the creation of predictive models
for clinical data; however, much existing literature reports results using
only a single source of data. In this work, we evaluate the performance of
models trained on the publicly-available eICU Collaborative Research Database.
We show that cross-validation using many distinct centers provides a
reasonable estimate of model performance in new centers. We further show that
a single model trained across centers transfers well to distinct hospitals,
even compared to a model retrained using hospital-specific data. Our results
motivate the use of multi-center datasets for model development and highlight
the need for data sharing among hospitals to maximize model performance.

##  1 Introduction

The intensive care unit (ICU) admits severely ill patients in order to provide
radical, life-saving treatment. ICUs frequently have high staff to patient
ratio in order to facilitate continuous monitoring of all patients and to
ensure that deterioration in patient condition is detected, and corrected,
before it becomes fatal. This approach is demonstrated to improve patient
outcomes [12]. As a result, the ICU is a data-rich environment.

Major efforts have been placed in utilizing this data to quantify patient
health and to predict future outcomes. One of the most immediately relevant
outcomes to the ICU is patient mortality. The APACHE system, first published
in [15], provides predictions for patient mortality based upon data collected
in the ICU. While the initial system was based off expert rules, later updates
used data driven methods [30]. Numerous other prediction systems have also
been developed [14, 17, 16, 28, 18, 7]; for a review of severity illness
scores in the ICU, see [26] and [13].

With recent advances in both machine learning and hardware for data archiving,
research has begun to return to building better prediction models using more
detailed granular data. The use of mortality prediction models to evaluate
ICUs as a whole has found great success, both for capturing patient severity
of illness and for risk-adjustment of patient populations [6, 25, 9, 2, 27,
21, 22, 19, 3, 4, 1, 20].

We hypothesize that a model trained on ICU data from one institution will
perform considerably worse when assessed at other ICUs, even when those
institutions use the same vendor for collection of data. This expectation of
model generalizability reflects known distributional differences in data that
result from differing care and data-recording practices among ICUs. Such
distributional differences complicate model generalizability in the absence of
specialized approaches [29, 5]. We further investigate a learning approach for
deploying a prediction model, comparing a transferred model to one that is
developed using only local data. We hypothesize that the local model will
initially have poor performance, but with sufficient sample size will
outperform the transferred model.

##  2 Data

We utilize the eICU Collaborative Research Database (eICU-CRD) v2.0 [11], a
publicly available multi-center database sourced from the Philips eICU
programme, a tele-medicine initiative where healthcare workers remotely
monitor acutely ill patients. All data in eICU-CRD has been de-identified. The
eICU-CRD contains 200,859 distinct unit stays, from which we exclude non-ICU
stays and ICU stays for which an APACHE IVa score is not available. The last
exclusion criteria implicitly applies APACHE IVa exclusion criteria in order
to mitigate issues around repeated sampling, ensure sufficient support, and
remove admissions which appear in the database for administrative purproses
[30]. After applying these exclusions, we removed hospitals which had less
than 500 stays, as we aimed to evaluate performance hospital-wise. The final
cohort contained 50,106 ICU stays across 46 hospitals. Details of this cohort
can be found in Appendix A, which contains demographic data generated using
tableone [24].

### Feature Extraction

For each ICU stay, we extracted data from a fixed window of length $W=24$
(hours) starting on ICU admission. Features extracted from the window of fixed
size $W$ are detailed in Appendix B. The features were extracted from a number
of physiologic and laboratory measurements. Notably, no explicit data
regarding treatment was extracted (e.g. use of vasopressors, mechanical
ventilation, dialysis, etc). In addition to those listed in Appendix 2, we
extracted gender, age, race, and whether the hospital admission was for an
elective surgery (binary covariate). A total of 82 features were extracted.

##  3 Methods

We evaluated logistic regression (LR) models trained using L2 regularization,
created using scikit-learn v0.20.0 [23]. The outcome evaluated was in-hospital
mortality. The area under the receiver operator characteristic curve (AUROC)
was used to evaluate discrimination of models, while the Standardized
Mortality Ratio (SMR)111Standardized Mortality Ratio (SMR) is the observed
proportion of patients who died divided by the expected proportion of patient
death (given probabilistic risks of patient mortality). was used to evaluate
the calibration of models.

### Experiments

We split our cohort evenly into training and test sets hospital-wise. We train
a single LR model on the training data and evaluate it on each hospital
separately via the AUROC and the SMR. We experimented with domain transfer but
it did not noticably impact model performance (Appendix C).

For our second experiment, we select a single hospital and set aside 500
patients. We then iteratively sample sets of 200 patients without replacement
as our recalibration dataset, i.e. we increase the size of this recalibration
set by 200 patients each iteration. At all iterations, we (i)
recalibrate222Recalibration involves training a univariable logistic
regression with hospital mortality as the outcome, effectively shifting and
scaling model predictions. the previous model developed using all hospitals in
the training set (“transferred“ model), and (ii) develop a new model using
only data available at the current iteration (“local” model). We evaluate
models on the 500 patients initially held out. To ensure stability of
estimates we repeat this procedure 50 times after randomly reshuffling data
and subsequently plot average performance for the transferred and local
models.

We recognize that it can be difficult to fully specify cohort definitions and
analyses approaches in critical care [10], and as such we have made our code
available for closer study [8].333https://github.com/alistairewj/icu-model-
transfer

##  4 Results

The AUROC across 5-fold CV on eICU-CRD was 0.854 and the SMR was 0.998. Figure
1 shows AUROC and SMR on held-out hospitals were distributed around the cross-
validation performance. The AUROC ranged from 0.688 to 0.933, and the SMR
ranged from 0.560 to 1.43.

Figure 1: Five fold cross validation (CV) performance of a LR model on the
training set (red hatched circle) and for a single model learned from the
training set and evaluated on held out hospitals (blue circles). Dashed lines
indicate the average of held out performance, and circle size corresponds to
hospital size (red hatched circle size is decimated for clarity).

Figure 2 shows the performance of the recalibration experiment on hospital 73,
the largest hospital in the test set. SMRs for both models are reasonable,
with the local model initially overpredicting mortality (SMR $>$ 1.0), and
eventually converging to the ideal value (1.0). The AUROC of the transferred
model is constant at 0.849, which is expected as the recalibration procedure
is a monotonic transform which cannot improve pairwise measures such as AUROC.
The AUROC of the local model is initially poor (0.614), and after 2,400
patients has significantly improved to 0.796, but is still below the
transferred model’s AUROC. Other hospitals assessed in this way exhibited
similar patterns (two other hospitals are shown in Appendix D).

Figure 2: Performance curves when (i) recalibrating a model developed using
other hospitals in eICU-CRD (“transferred model”), and (ii) when training a
model using only local hospital data (“local model”). Data shown from top left
to bottom left: SMRs for the local model (blue triangles), SMRs for the
transferred model (red circles), AUROC for the transferred model (red line),
and AUROC for the local model (blue line).

##  5 Discussion and Related Work

We have evaluated the performance of a model trained using data from many
centers and applied to held out hospitals from the continental United States.
We found that the cross-validation estimate of performance consistently over-
estimated performance on the external hospitals but not by an impractical
amount. Notable was the large variance in hospital performance, and there was
no consistent pattern of degradation in performance.

Our recalibration experiment assessed how well models transfer to held out
hospitals as compared to locally developed models. Local models tended to be
well calibrated after reasonable numbers of patients. In terms of
discrimination, the local model would likely continue to improve if more data
were available in eICU-CRD, but it is noteworthy that even after large sample
sizes these models still had inferior discrimination to the transferred
models. An extremely common approach in healthcare is for hospitals to
quarantine their data, often citing reasonable privacy concerns, and develop
models internally. While the amount of data archived from healthcare
institutions has dramatically increased in the past decade, many prediction
problems of interest require specialized cohorts which results in low sample
sizes. Isolating model development to individual hospitals appears to be
detrimental to overall model performance, and significant gains are available
when data from multiple hospitals is pooled for model development. Our results
indicate that research into safe and secure ways of sharing patient data may
be a good avenue for improving performance of predictive models.

There are limitations to our study. We only evaluated mortality prediction
models, and did not assess other types of predictive models. However, we
believe our results will hold for other predictive models in the ICU. Our
feature set was also a small subset of what is available in eICU-CRD.
Specifically, we did not utilize any data regarding patient treatment, high
granularity vital signs, or patient intake. The incorporation of information
from these sources may improve model performance and generalizability.
Finally, our results only assess logistic regression as a modelling
methodology, though we have found similar results for non-linear approaches
such as gradient boosting.

##  6 Conclusion

We have shown that classic in-hospital mortality models trained on publicly-
accessible data exhibit varying degrees of decreased performance when applied
to external hospital data, both with and without hospital-specific
normalization. The results agree with common intuition about relative model
performance decreasing on external datasets, but provide empirical assessment
of the magnitude of this effect. We further showed that large sample sizes are
needed from individual hospitals before performance is close to that achieved
by a transfered model from a pooled collection of distinct training data. This
result highlights the importance of data sharing in the development of high
performance predictive models.

## Acknowledgements

This research was funded in part by grants from the National Institutes of
Health (NIH): NIH-R01-EB017205, NIH-R01-EB001659, and NIH-R01-GM104987.

## References

  * [1] Karla L Caballero Barajas and Ram Akella.  Dynamically modeling patient’s health state from electronic medical records: A time series approach.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 69–78. ACM, 2015.
  * [2] Luca Citi and Riccardo Barbieri.  Physionet 2012 challenge: Predicting mortality of icu patients using a cascaded svm-glm paradigm.  In Computing in Cardiology (CinC), 2012, pages 257–260. IEEE, 2012\.
  * [3] Marzyeh Ghassemi, Tristan Naumann, Finale Doshi-Velez, Nicole Brimmer, Rohit Joshi, Anna Rumshisky, and Peter Szolovits.  Unfolding physiological state: Mortality modelling in intensive care units.  In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 75–84. ACM, 2014.
  * [4] Marzyeh Ghassemi, Marco A F Pimentel, Tristan Naumann, Thomas Brennan, David A Clifton, Peter Szolovits, and Mengling Feng.  A multivariate timeseries modeling approach to severity of illness assessment and forecasting in icu with sparse, heterogeneous clinical data.  In Proceedings of the… AAAI Conference on Artificial Intelligence. AAAI Conference on Artificial Intelligence, volume 2015, page 446\. NIH Public Access, 2015.
  * [5] Jen J Gong, Thoralf M Sundt, James D Rawn, and John V Guttag.  Instance weighting for patient-specific risk stratification models.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 369–378. ACM, 2015.
  * [6] Caleb W Hug and Peter Szolovits.  Icu acuity: Real-time models versus daily models.  In AMIA, 2009.
  * [7] A. E. W. Johnson, Andrew A Kramer, and Gari D Clifford.  A new severity of illness scale using a subset of acute physiology and chronic health evaluation data elements shows comparable predictive accuracy.  Critical Care Medicine, 41(7):1711–1718, 2013.
  * [8] Alistair Johnson.  alistairewj/icu-model-transfer: ICU Model Transfer v1.0, December 2018\.
  * [9] Alistair E W Johnson, Nic Dunkley, Louis Mayaud, Athanasios Tsanas, Andrew A Kramer, and Gari D Clifford.  Patient specific predictions in the intensive care unit using a bayesian ensemble.  In Computing in Cardiology (CinC), 2012, pages 249–252. IEEE, 2012\.
  * [10] Alistair E. W. Johnson, Tom J. Pollard, and Roger G. Mark.  Reproducibility in critical care: a mortality prediction case study.  In Finale Doshi-Velez, Jim Fackler, David Kale, Rajesh Ranganath, Byron Wallace, and Jenna Wiens, editors, Proceedings of the 2nd Machine Learning for Healthcare Conference, volume 68 of Proceedings of Machine Learning Research, pages 361–376, Boston, Massachusetts, 18–19 Aug 2017. PMLR.
  * [11] Alistair EW Johnson, Tom J Pollard, and Roger G Mark.  eicu collaborative research database, 2017.
  * [12] R L Kane, T A Shamliyan, C Mueller, S Duval, and T J Wilt.  The association of registered nurse staffing levels and patient outcomes: Systematic review and meta-analysis.  Medical Care, 45(12):1195–1204, December 2007.
  * [13] Mark T Keegan, Ognjen Gajic, and Bekele Afessa.  Severity of illness scoring systems in the intensive care unit.  Critical care medicine, 39(1):163–169, 2011.
  * [14] W A Knaus, D P Wagner, E A Draper, J E Zimmerman, M Bergner, C A Bastos, P G snd Sirio, D J Murphy, T Lotring, A Damiano, and F E Harrell Jr.  The APACHE III prognostic system: Risk prediction of hospital mortality for critically ill hospitalized adults.  Chest, 100(6):1619–1636, 1991.
  * [15] W A Knaus, J E Zimmerman, D P Wagner, E A Draper, and D E Lawrence.  APACHE-acute physiology and chronic health evaluation: A physiologically based classification system.  Critical Care Medicine, 9:591–597, 1981.
  * [16] J R Le Gall, S Lemeshow, and F Saulnier.  A new Simplified Acute Physiology Score (SAPS-II) based on a European/North American multicenter study.  JAMA, 270(24):2957–63, 1993.
  * [17] J R Le Gall, P Loirat, F Nicolas, C Granthil, F Wattel, R Thomas, P Glaser, P Mercier, J Latournerie, P Candau, et al.  Use of a severity index in 8 multidisciplinary resuscitation centers.  Presse medicale (Paris, France: 1983), 12(28):1757–1761, 1983.
  * [18] Jean Roger Le Gall, Janelle Klar, Stanley Lemeshow, Fabienne Saulnier, Corinne Alberti, Antonio Artigas, and Daniel Teres.  The logistic organ dysfunction system: a new way to assess organ dysfunction in the intensive care unit.  Jama, 276(10):802–810, 1996.
  * [19] Li-Wei H Lehman, Mohammed Saeed, William J Long, Joon Lee, and Roger G Mark.  Risk stratification of icu patients using topic models inferred from unstructured progress notes.  In AMIA. Citeseer, 2012.
  * [20] Yuan Luo, Yu Xin, Rohit Joshi, Leo Celi, and Peter Szolovits.  Predicting icu mortality risk by grouping temporal trends from a multivariate panel of physiologic measurements.  In AAAI, pages 42–50, 2016.
  * [21] Martin Macaš, Jakub Kuzilek, Tadeáš Odstrčilík, and Michal Huptych.  Linear bayes classification for mortality prediction.  In Computing in Cardiology (CinC), 2012, pages 473–476. IEEE, 2012\.
  * [22] Sean McMillan, Chih-Chun Chia, Alexander Van Esbroeck, Ilan Rubinfeld, and Zeeshan Syed.  Icu mortality prediction using time series motifs.  In Computing in Cardiology (CinC), 2012, pages 265–268. IEEE, 2012\.
  * [23] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  Journal of Machine Learning Research, 12(Oct):2825–2830, 2011.
  * [24] Tom J Pollard, Alistair EW Johnson, Jesse D Raffa, and Roger G Mark.  tableone: An open source python package for producing summary statistics for research papers.  JAMIA Open, 2018.
  * [25] Ikaro Silva, George Moody, Daniel J Scott, Leo A Celi, and Roger G Mark.  Predicting in-hospital mortality of icu patients: The physionet/computing in cardiology challenge 2012.  In Computing in Cardiology (CinC), 2012, pages 245–248. IEEE, 2012\.
  * [26] K Strand and H Flaatten.  Severity scoring in the icu: A review.  Acta Anaesthesiologica Scandinavica, 52(4):467–478, 2008.
  * [27] Srinivasan Vairavan, Larry Eshelman, Syed Haider, Abigail Flower, and Adam Seiver.  Prediction of mortality in an intensive care unit using logistic regression and a hidden markov model.  In Computing in Cardiology (CinC), 2012, pages 393–396. IEEE, 2012\.
  * [28] J-L Vincent, R Moreno, J Takala, S Willats, A De Mendoca, H Bruining, C K Reinhart, P M Suter, and L G Thijs.  The SOFA (Sepsis-related Organ Failure Assessment) score to describe organ dysfunction/failure.  Intensive Care Medicine, 22:707–710, 1996.
  * [29] Jenna Wiens, John Guttag, and Eric Horvitz.  A study in transfer learning: Leveraging data from multiple hospitals to enhance hospital-specific predictions.  Journal of the American Medical Informatics Association, 21(4):699–706, 2014.
  * [30] J E Zimmerman and A A Kramer.  Outcome prediction in critical care: the Acute Physiology and Chronic Health Evaluation models.  Current Opinion in Critical Care, 14:491–497, 2008.

##  Appendix A Demographics for eICU-CRD

Demographics for the cohort are shown in Table 1.

Demographic | Value
---|---
Number of ICU stays | 50106
Hospital LOS (days) | 5.53 [3.27,9.13]
Male | 27833 (55.55)
Female | 22273 (44.45)
Age (years) | 63.00 [52.00,74.00]
Elective Surgery | 11348 (22.65)
Race |
White | 40998 (81.82)
African Americans | 6720 (13.41)
Asian | 631 (1.26)
Hispanic | 1757 (3.51)
Table 1: Demographics of the eICU-CRD cohort. Binary variables are represented
as a count with percentage, normally distributed continuous variables are
represented as means with standard deviations, and non-normally distributed
variables are represented as medians with 25th and 75th percentiles.

##  Appendix B Full list of features and window sizes

Time window |  Feature extracted |  Variables
---|---|---
$[t_{i,w}-W,t_{i,w}]$ |  First, Last |  Heart rate, Systolic/Diastolic/Mean blood pressure, Respiratory rate, Temperature, Oxygen Saturation, Glasgow coma scale
$[t_{i,w}-W-24,t_{i,w}]$ |   First, last |  Base excess, Calcium, Partial pressure of oxygen in arterial blood, Partial pressure of carbon dioxide in arterial blood, pH, Ratio of partial pressure of oxygen to fraction of oxygen inspired, Total carbon dioxide concentration
$[t_{i,w}-W-24,t_{i,w}]$ |   First, last |  Albumin, Blood urea nitrogen, Immature band forms, Bicarbonate, Bilirubin, Creatinine, Hematocrit, Hemoglobin, Lactate, Platelet, Potassium, International Normalized Ratio, Sodium, White blood cell count
$[t_{i,w}-W,t_{i,w}]$ |  Sum |  Urine output
Table 2: Features extracted during the first day of the patient’s ICU stay.
Blood gases and laboratory measurements have the same feature extraction
(first, last), but are separated for clarity.  Note that some features have
had their window extended backward an extra 24 hours due to infrequent
collection.

Features extracted used consistent functional forms: the first, last, or sum
(in the case of urine output) value was extracted from all measurements of the
variable made within first 24 hours. As laboratory values are less frequently
sampled, this window was extended 24 hours backward for these measurements.

##  Appendix C Domain transfer for improving model performance

For each hospital in the test set of eICU-CRD, we rescaled the data using
means and standard deviations as calculated on the test set hospital. We
compare this to our standard approach which uses the mean and standard
deviation of the entire training set to rescale the data. Figure 3 shows the
performance across test hospitals as measured by the AUROC and SMR. We found
no noticably difference between the two methods.

Figure 3: Evaluation of domain transfer on eICU-CRD. Five fold cross
validation (CV) performance of a LR model on the training set (red hatched
circle) and for a single model learned from the training set and evaluated on
held out hospitals (blue circles). Held out hospital data was rescaled using
their respective means and standard deviations. Dashed lines indicate the
average of held out performance, and circle size corresponds to hospital size
(red hatched circle size is decimated for clarity).

##  Appendix D Transfer performance in other hospitals

Figure 4 shows AUROC performance of a model developed on the eICU-CRD training
set and recalibrated iteratively to a held-out hospital dataset as described
in the paper. Figure 5 shows the same results for a distinct hospital.

Figure 4: Performance curve when (i) recalibrating a model developed using
other hospitals in eICU-CRD (“transferred model”), and (ii) when training a
model using only local hospital data (“local model”). The top red line shows
the AUROC for the transferred model, and the lower blue line shows the AUROC
for the local model. Figure 5: Performance curve when (i) recalibrating a
model developed using other hospitals in eICU-CRD (“transferred model”), and
(ii) when training a model using only local hospital data (“local model”). The
top red line shows the AUROC for the transferred model, and the lower blue
line shows the AUROC for the local model.


# Concurrence for well-formed CAFs:
Naive Semantics

Rafael Kiesel, Anna Rapberger

(TU Wien, Institute of Logic and Computation, Austria
<EMAIL_ADDRESS>
###### Abstract

In the area of claim-based reasoning in abstract argumentation, a claim-based
semantics is said to be concurrent in a given framework if all its variants
yield the same extensions. In this note, we show that the concurrence problem
with respect to naive semantics is $\coNP$-hard for well-formed CAFs. This
solves a problem that has been left open in [2].

##  1 Introduction

_Claim-augmented argumentation frameworks (CAFs)_ [4] extend _abstract
argumentation frameworks (AFs)_ [1] by a function that assigns a claim to each
argument. Formally, CAFs correspond to directed labeled graphs where the nodes
correspond to arguments, the labels correspond to the arguments claims, and
the arcs indicate (unidirectional) conflicts between them. As with arguments
in AFs, the acceptance status of claims is decided via _argumentation
semantics_. There are often several possibilities to lift AF semantics to
claim-level [3]; considering _naive semantics_ which is based on maximization
of conflict-free sets gives rise to two claim-based variants: The so-called
_inherited_ variant, which performs maximization on argument-level, and the
_claim-level semantics_ , which performs maximization over claim-sets.
Deciding whether two different variants of a semantics yield the same claim-
based outcome for a given CAF (the so-called _concurrence problem_) can be
computationally challenging: As shown in [2], the concurrence problem is
$\coNP$-complete with respect to naive semantics.

In this note, we show that deciding concurrence with respect to naive
semantics is $\coNP$-hard even for _well-formed CAFs_ , an important sub-class
of CAFs that impose restrictions on the attack relation. Well-formed CAFs
satisfy a quite natural behavior of conflicts: A CAF is well-formed if
arguments with the same claims attack the same arguments. This closes the
complexity gap in [2] for the concurrence problem.

##  2 Preliminaries

We introduce abstract argumentation and claim-based reasoning following [1]
and [3].

### Abstract Argumentation.

We fix a non-finite background set $U$. An argumentation framework (AF) [1] is
a directed graph $F=(A,R)$ where $A\subseteq U$ represents a set of arguments
and $R\subseteq A\times A$ models attacks between them. For two arguments
$a,b\in A$, if $(a,b)\in R$ we say that $a$ attacks $b$ as well as $a$ attacks
(the set) $E$ given that $b\in E\subseteq A$. A set $E\subseteq A$ is
_conflict-free_ in $F$ ($E\in\mathit{cf}(F)$) iff for no $a,b\in E$, $(a,b)\in
R$. A _semantics_ is a function $\sigma$ with $F\mapsto\sigma(F)\subseteq
2^{A}$. In this note we focus on _naive semantics_.

######  Definition 1.

For an AF $F=(A,R)$, a set $E\subseteq A$ is _naive_
($E\in\mathit{na}(\mathcal{F})$) iff $E$ is $\subseteq$-maximal in
$\mathit{cf}(F)$.

### Reasoning about claims.

A _claim-augmented argumentation framework (CAF)_ [4] is a triple
$\mathcal{F}=(A,R,\mathit{cl})$ where $F=(A,R)$ is an AF and $\mathit{cl}$ is
a function which assigns a claim to each argument in $A$. The claim-function
is extended to sets in the natural way, i.e., for a set $E\subseteq A$, we let
$\mathit{cl}(E)=\\{\mathit{cl}(a)\mid a\in E\\}$. A CAF is _well-formed_ iff
$a^{+}_{F}=b^{+}_{F}$ for all $a,b\in A$ with $\mathit{cl}(a)=\mathit{cl}(b)$,
i.e., arguments with the same claim attack the same arguments. The literature
offers several ways to extend semantics for AFs that involve claims
respectively arguments to a different extent. We introduce _inherited_ and
_claim-level semantics_ for naive semantics that perform maximization in
different stages of the evaluation.

######  Definition 2.

For a CAF $\mathcal{F}=(A,R,\mathit{cl})$, $F=(A,R)$, and a semantics
$\sigma$, we let $\sigma\\!_{c}(\mathcal{F})=\\{\mathit{cl}(E)\mid
E\in\sigma(F)\\}$. A set $S\subseteq\mathit{cl}(A)$ is

  * •

_i-naive_ ($S\in\mathit{na}\\!_{c}(\mathcal{F})$) iff
$S\in\mathit{na}\\!_{c}(\mathcal{F})$, i.e., there is $E\subseteq A$ with
$S=\mathit{cl}(E)$ and $E$ is $\subseteq$-maximal in $\mathit{cf}(F)$;

  * •

_cl-naive_ ($S\in\operatorname{\mathit{cl-naive}}(\mathcal{F})$) $S$ is
$\subseteq$-maximal in $\mathit{cf}\\!_{c}(\mathcal{F})$.

As shown in [3], $\operatorname{\mathit{cl-
naive}}(\mathcal{F})\subseteq\mathit{na}\\!_{c}(\mathcal{F})$.

##  3 The Concurrence Problem for Naive Semantics

An interesting problem that arises when considering different variants of
semantics for CAFs is the so-called _concurrence problem_.

$\textit{Con}_{\sigma}^{\Delta}$, $\Delta\in\\{\mathit{CAF},\mathit{wf}\\}$
---
Input: |  A CAF $\mathcal{F}$ (if $\Delta=\mathit{CAF}$)/a well-formed CAF $\mathcal{F}$ (if $\Delta=\mathit{wf}$)
Output: |  true iff $\sigma\\!_{c}(\mathcal{F})=\operatorname{\mathit{cl-\sigma}}(\mathcal{F})$

For naive semantics, the problem can be formulated as follows: Given a CAF
$\mathcal{F}$, is it the case that maximization on argument-level yields the
same accepted sets as maximization on claim-level?

In [2] it has been shown that $\textit{Con}_{\sigma}^{\mathit{CAF}}$ is
$\coNP$-complete. Since concurrence for well-formed CAFs with respect to naive
semantics is a special case of CAFs we obtain upper bounds for
$\textit{Con}_{\sigma}^{\mathit{wf}}$, that is,
$\textit{Con}_{\sigma}^{\mathit{wf}}$ is in $\coNP$.

We restate the following proposition [2].

######  Proposition 1.

For a CAF $\mathcal{F}=(A,R,\mathit{cl})$,
$\mathit{na}\\!_{c}(\mathcal{F})=\operatorname{\mathit{cl-
naive}}(\mathcal{F})$ if and only if $naive\\!_{c}(\mathcal{F})$ is
incomparable.

Thus it suffices to verify incomparability of
$\mathit{na}\\!_{c}(\mathcal{F})$. An $\NP$ procedure for the complementary
problem is by a standard guess and check procedure: Guess $E,G\subseteq A$ and
check (i) $E,G\in\sigma((A,R))$ and (ii)
$\mathit{cl}(E)\subset\mathit{cl}(G)$. The former can be checked in time
polynomial in the number of arguments in $\mathcal{F}$.

Next we show that verifying incomparability of
$\mathit{na}\\!_{c}(\mathcal{F})$ is $\NP$-hard even if $\mathcal{F}$ is well-
formed. We will first define the base reduction, which we slightly extend to
obtain the the reduction for the well-formed case.

######  Reduction 1.

Let $\varphi$ be given by a set of clauses $C=\\{cl_{1},\dots,cl_{n}\\}$ over
atoms in $X$. We construct $(A,R,\mathit{cl})$ with

$\displaystyle A=$ $\displaystyle X\cup\bar{X}\cup
C\cup\\{\varphi\\}\cup\\{a_{1},a_{2}\\},\text{ with }\bar{X}=\\{\bar{x}\mid
x\in X\\},$ $\displaystyle R=$ $\displaystyle\\{(x,cl)\mid cl\in C,x\in
cl\\}\cup\\{(\bar{x},cl)\mid cl\in C,\neg x\in cl\\}$
$\displaystyle\cup\\{(x,\bar{x}),(\bar{x},x)\mid x\in
X\\}\cup\\{(cl_{i},\varphi)\mid i\leq n\\}\cup\\{(\varphi,a_{2})\\},$

and $\mathit{cl}(x)=x$,
$\mathit{cl}(\bar{x})=\bar{x},\mathit{cl}(cl_{i})=cl_{i},\mathit{cl}(\varphi)=\varphi$
and $\mathit{cl}(a_{i})=a$.

An example of this reduction is given in Figure 1.

$\varphi$$a_{2}$a$a_{1}$a$cl_{1}$$cl_{2}$$cl_{3}$$x_{1}$$\bar{x}_{1}$$x_{2}$$\bar{x}_{2}$$x_{3}$$\bar{x}_{3}$$x_{4}$$\bar{x}_{4}$
Figure 1:  Reduction 1 for a formula $\varphi$ which is given by the clauses
$\\{\\{x_{1},x_{3},x_{4}\\},\\{\bar{x}_{3},\bar{x}_{4},\bar{x}_{2})\\},\\{\bar{x}_{1},\bar{x}_{3},x_{2}\\}\\}$.

The intuition behind this reduction is the following: We take the smallest
sub-CAF that is not concurrent. This corresponds to
$(\\{a_{1},a_{2},\varphi\\},\\{(\varphi,a_{2})\\},\mathit{cl})$ with
$\mathit{cl}(a_{i})=a$ and $\mathit{cl}(\varphi)=\varphi$, with naive
extensions $\\{a_{1},\varphi\\}$ and $\\{a_{1},a_{2}\\}$. Observe that the
second extension requires that $\varphi$ is not taken. We then add the
additional part corresponding to a propositional formula in CNF that only
allows $\varphi$ to not be included without additional changes, when the
formula is satisfiable. Thus allowing us to reduce the concurrence problem to
UNSAT.

In the example in Figure 1 we see that
$\\{a_{1},\varphi,x_{1},x_{2},x_{3},\bar{x_{4}}\\}$ and
$\\{a_{1},a_{2},x_{1},x_{2},x_{3},\bar{x_{4}}\\}$ are naive extensions because
$\\{x_{1},x_{2},x_{3},\bar{x_{4}}\\}$ is a satisfying assignment of the
formula.

Generally, the proof proceeds as follows.

######  Proposition 2.

$\textit{Con}_{\mathit{na}}^{\mathit{wf}}$ is $\coNP$-hard.

###### Proof.

For hardness, we present a reduction from UNSAT: Let $\varphi$ be given by a
set of clauses $C=\\{cl_{1},\dots,cl_{n}\\}$ over literals in $X$. W.l.o.g. we
can assume that $\varphi$ does not contain tautological clauses, i.e., there
is no $cl_{i}$, $i\leq n$ with $x,\bar{x}\in cl_{i}$ for any $x\in X$. Let
$(A,R,\mathit{cl})$ be defined as in Reduction 1. We will show $\varphi$ is
unsatisfiable iff $\mathit{na}\\!_{c}(\mathcal{F})$ is incomparable.

First asssume $\varphi$ is satisfiable and consider a model $M$ of $\varphi$.
Let $E=M\cup\\{\bar{x}\mid x\notin M\\}\cup\\{\varphi,a_{1}\\}$. Clearly, $E$
is conflict-free; moreover, as $M$ satisfies each clause $cl_{i}$ there is
either $x\in cl_{i}$ with $x\in M$ or $\bar{x}\in cl_{i}$ with $x\notin M$,
thus $E$ attacks each $cl_{i}$. Since, also $\varphi$ and $a_{1}$ are in $E$
we the only argument left to consider is $a_{2}$, which is however attacked by
$\varphi$ and can therefore not be included, while preserving conflict-
freeness. We can conclude that $E$ is a subset-maximal conflict-free set.
Moreover, $E^{\prime}=M\cup\\{\bar{x}\mid x\notin M\\}\cup\\{a_{1},a_{2}\\}$
is also a subset-maximal conflict-free set, since still every $cl_{i}$ is
attacked and $\varphi$ attacks $a_{2}$. It follows that
$\mathit{na}\\!_{c}(\mathcal{F})$ is not incomparable since
$\mathit{cl}(E)=M\cup\\{\bar{x}\mid x\notin M\\}\cup\\{\varphi,a\\}$ is a
strict superset of $M\cup\\{\bar{x}\mid x\notin
M\\}\cup\\{a\\}=\mathit{cl}(E^{\prime})$ and both are contained in
$\mathit{na}\\!_{c}(\mathcal{F})$.

Now assume that $\varphi$ is unsatisfiable. Let $E$ be a subset-maximal
conflict-free set. If $\varphi\in E$ it follows that none of the $cl_{i}$ are
in $E$. Therefore, it holds that for each $x$ exactly one of $x$ and $\bar{x}$
is in $E$. Furthermore, it always holds that $a_{1}$ is in $E$. This means
that $\mathit{cl}(E)=\\{\varphi,a\\}\cup\\{x\mid x\in
X^{\prime}\\}\cup\\{\bar{x}\mid x\not\in X^{\prime}\\}$ for some subset
$X^{\prime}$ of $X$.

Case I: Assume there is a naive extension $E^{\prime}$ such that
$\mathit{cl}(E)\subsetneq\mathit{cl}(E^{\prime})$. Then $E^{\prime}$ must
contain $\mathit{cl}_{i}$ for some $i$, $x$ for some $x\not\in X^{\prime}$ or
$\bar{x}$ for some $x\in X^{\prime}$. This is not possible, since this would
imply $E\subsetneq E^{\prime}$, which is a contradiction to the assumption
that $E$ and $E^{\prime}$ are naive extensions.

Case II: Assume there is a naive extension $E^{\prime}$ that contains
$\varphi$ such that $\mathit{cl}(E^{\prime})\subsetneq\mathit{cl}(E)$. Then
there must be some $x$ with $x\in X^{\prime}$ or $\bar{x}$ with $x\not\in
X^{\prime}$ that is not contained $E^{\prime}$. This is not possible, since
this would imply $E^{\prime}\subsetneq E$, which is a contradiction to the
assumption that $E$ and $E^{\prime}$ are naive extensions.

Case III: Assume there is a naive extension $E^{\prime}$ that does not contain
$\varphi$ such that $\mathit{cl}(E^{\prime})\subsetneq\mathit{cl}(E)$. This is
impossible, as $\varphi$ is unsatisfiable, which entails that there is a
clause $cl_{i}$ that is not satisfied by $X^{\prime}$. Therefore, $cl_{i}$ is
not attacked by $E\setminus\varphi$ and also not $E^{\prime}$. Thus
$E^{\prime}$ contains at least one $cl_{i}$, which means that also
$\mathit{cl}(E^{\prime})$ contains $cl_{i}$.

We see that any naive extension $E$ that contains $\varphi$ is incomparable to
any other naive extension on claim level.

Next we show that the same holds for any naive extension $E$ that does not
contain $\varphi$. It follows that $E$ contains $a_{1},a_{2}$. Since
$\mathcal{F}^{\prime}=(A\setminus\\{a_{1},a_{2},\varphi\\},R\setminus\\{(\varphi,a_{2})\\}\cup\\{(cl_{i},\varphi)\mid
i=1,\dots,n\\},\mathit{cl})$, the rest of the CAF, is such that all arguments
have distinct claims, there is a one to one correspondence between
$\mathit{na}\\!_{c}(\mathcal{F}^{\prime})$ and
$\mathit{na}(\mathcal{F}^{\prime})$, which implies that they are all
incomparable with one another. ∎

##  4 Conclusion

We see that even for well formed CAFs the concurrence problem is
$\coNP$-complete, as it is the case for general CAFs. This shows that while
well formedness is an interesting property that only allows a fragment of the
CAFs that might be deemed more reasonable, it does not lead to CAFs that are
simpler, at least with respect to the concurrence problem.

## Acknowledgements

This work has been supported by the Austrian Science Fund (FWF) Grant
W1255-N23.

## References

  * [1] Phan Minh Dung.  On the acceptability of arguments and its fundamental role in nonmonotonic reasoning, logic programming and n-person games.  Artificial Intelligence, 77(2):321–358, 1995.
  * [2] Wolfgang Dvorák, Alexander Greßler, Anna Rapberger, and Stefan Woltran.  The complexity landscape of claim-augmented argumentation frameworks.  In AAAI 2021, Proc., pages 6296–6303. AAAI Press, 2021.
  * [3] Wolfgang Dvořák, Anna Rapberger, and Stefan Woltran.  Argumentation semantics under a claim-centric view: Properties, expressiveness and relation to SETAFs.  In KR 2020, Proc., pages 341–350. IJCAI.org, 2020.
  * [4] Wolfgang Dvořák and Stefan Woltran.  Complexity of abstract argumentation under a claim-centric view.  Artificial Intelligence, 285:103290, 2020.


# Cosmological constraints on the $R^{2}$-corrected Appleby-Battye model

B. W. Ribeiro  A. Bernui  and M. Campista

###### Abstract

Nowadays, efforts are being devoted to the study of alternative cosmological
scenarios, in which, modifications of General Relativity (GR) theory have been
proposed to explain the late cosmic acceleration, without assuming the
existence of the dark energy (DE) component. We investigate the
$R^{2}$-corrected Appleby-Battye model, or $R^{2}$-AB model, which consists of
an $f(R)$ model with only one extra free parameter $b$, besides the
cosmological parameters of the flat-$\Lambda$CDM model: $H_{0}$ and
$\Omega_{m,0}$. Regarding this model, it was already shown that a positive
value for $b$ is required for the model to be consistent with Solar System
tests, moreover, the condition for the existence of a de Sitter state requires
$b\geq 1.6$. To impose observational constraints on the $R^{2}$-AB model, we
consider in our analyses two data sets: cosmic chronometer $H(z)$ data for the
background level, and $[f\sigma_{8}](z)$ data, for the perturbative level. The
first one provides $b=1.6^{+3.1}_{-0.0}$ and the cosmological parameters
$\\{H_{0},\Omega_{m,0}\\}$ in agreement to Planck values, while the second
one, indicates $b=1.76^{+2.91}_{-0.15}$ and the parameters
$\\{\Omega_{m,0},\sigma_{8,0}\\}$ also in agreement to Planck values; in the
last case the data was marginalized over the parameter $H_{0}$. Additionally,
we perform illustrative analyses that compare this $f(R)$ model with the
flat-$\Lambda$CDM model, considering several values of the parameter $b$, for
diverse cosmological functions like the Hubble function $H(z)$, the equation
of state $w_{eff}(z)$, the parametrized growth rate of cosmic structures
$[f\sigma_{8}](z)$, and $\sigma_{8}(z)$. The overall conclusion is that the
$R^{2}$-AB model is a promising $f(R)$ model that deserves to continue being
tested with diverse cosmological data.

##  1 Introduction

The recognition of the flat-$\Lambda$CDM model as the best model to describe
the Universe is almost a consensus, and this is undoubted because, besides
explaining the recent phase of accelerated cosmic expansion discovered by
Riess et al. [1] and Perlmutter et al. [2], this model is in agreement with a
plethora of observations, e.g., CMB, BAO, SNe, gravitational lensing, etc.,
motivating to call it the concordance cosmological model. However, this does
not mean that it will remain the concordance model. The large volume of
available data and the evolution of statistical procedures have transported
cosmology from the era of precision to the era of accuracy. Recent results
have revealed several tensions at large and small scales, hitherto hidden [3,
4]. These tensions, added to the unknown nature of the dark energy (DE)
component, have motivated studies of alternative scenarios (see, e.g., [7, 5,
6]).

The $\Lambda$CDM-type cosmological models, which include the flat-$\Lambda$CDM
case, assume the General Relativity (GR) theory as the metric theory. On the
other hand, alternative cosmological models propose to give up this
hypothesis, by considering other metric theories as extra dimensions, extra
fields, and higher orders corrections, termed modified gravity (MG) theories.
In some cases, the cosmic accelerated expansion comes from the gravity model
assumed, and the DE is no longer needed.

Generically, an MG theory assumes that the metric theory of the $n+1$
dimensional space-time describing the universe phenomena at quantum and large
scales comes from a suitable modification of the GR. This is because at Solar
System scale, and also in the distant past $z\gg 1$, GR has passed with honors
a set of astrophysical observations confirming its validity. Thus, the first
MG theories were proposed to renormalize the GR theory and to obtain a
classical theory of gravitation as the low-energy limit of the quantum gravity
[8]. Recently, some MG theories have been proposed to explain the cosmic
acceleration of the universe [9]. However, this does not necessarily solve the
issue because a MG model can satisfactorily explain the dynamics of the
universe at the background level, but can not correctly explain phenomena at
the perturbative level as, e.g., in the case of negative power index law
$f(R)$ models, that, due to the absence of a regular phase of matter, are
incompatible with the formation of structures [10], and also in the case of
logarithmic $f(R)$ model [11]. Other models, instead, are more interesting
because they adequately explain cosmological data, both in the background and
perturbative levels [12, 13, 14, 15, 16].

One of the simplest modifications of the GR theory is obtained by replacing
the term $R-2\Lambda$ in the Einstein-Hilbert (EH) Lagrangian by an arbitrary
function of the Ricci scalar $R$, $f(R)$. These theories, best known as $f(R)$
theories, are conformally equivalent to Einstein’s theory, with the addition
of an extra degree of freedom in the gravitational sector, the scalaron, a
canonical scalar field whose potential is uniquely determined by the scalar
curvature $R$ [17, 18].

The first successful $f(R)$ model was proposed by Starobinsky in order to
explain the primordial inflationary era, with the generic form $R+\alpha
R^{2}$, $\alpha>0$ [19], and from then on, several $f(R)$ models have been
proposed, considering from simple polynomial laws to more complicated
functions of the Ricci scalar (see, e.g., Amendola et al. [10], Starobinsky
[20], Hu-Sawicki [21], Appleby-Battye [22], Li-Barrow [23], Amendola-Tsujikawa
[24], Tsujikawa [25], Cognola et al. [26], Linder [27], Elizalde et al. [28],
Xu-Chen [29], Nautiyal et al. [30], Gogoi-Goswami [31] and Oikonomou [32,
33]).

Notice, however, that there is a degeneracy because diverse $f(R)$ models
correctly describe the accelerated cosmic expansion as the flat-$\Lambda$CDM
does. In these cases, one has to go to the perturbative level to decide which
model reproduces better the matter clustering in the observed universe [35,
34, 36].

It is important to mention that the alternative models have a larger number of
parameters when compared to the flat–$\Lambda$CDM, with only one free
parameter, and this is the reason why these models undergo a statistically
less efficient process of best-fitting cosmological data. However, this is not
always true when Bayesian statistical analysis is considered since in Bayesian
approach the comparison can vary significantly with the prior choices.

In this work we shall study the $R^{2}$-corrected Appleby-Battye (AB) model,
proposed in [37], that for the sake of simplicity, it will be denoted
throughout the text by $R^{2}$-AB model. This model results from the
improvement of the original AB model [22], where a term proportional to
$R^{2}$, with a sufficiently small coefficient to ensure the existence of the
primordial inflation, was added to solve the weak curvature singularity
problem [37, 38, 39], present in a number of $f(R)$ models. The reasons for
analyzing this model are diverse; firstly, the fact that it passed many
important tests (e.g., classical and semi-classical stability, solar system
constraints, correct primordial nucleosynthesis of light elements, and has
radiation, matter, and DE epochs), makes this model a viable alternative to
explain the current accelerated epoch. Besides, there are no analyses of this
model involving cosmological data to investigate its model parameters. We
emphasize that, if this model shows a good agreement with the observational
data, it will provide a geometrical explanation for the accelerated expansion,
not being necessary to assume a (non-physical) dark energy component in the
universe.

This work is organized as follows. In section 2 we introduce the formal basis
for $f(R)$ Gravity within the framework of the metric formalism. In section 3
a brief description of the $f(R)$ Cosmology is presented. Throughout these
sections, we highlight important sets of constraints that the $f(R)$ models
must obey. Next, in section 4, we present the main aspects of the $R^{2}$-AB
model and compare it theoretically with the flat-$\Lambda$CDM model. In
section 5 we provide two datasets and details on their compilations,
measurements, surveys, and cosmological tracers. Finally, we describe the
statistical methodology of the analysis performed and show our results in
section 6, as well as address our conclusions in section 7. Additionally, we
provide some plots in appendix A emphasizing the role of the free parameter of
the $R^{2}$-AB model versus the flat-$\Lambda$CDM model.

Throughout this work, we assume natural units, in which $\hbar=c=1$. Greek
symbols range from $0$ to $i$, with $0$ being the cosmic time and $i=1,2,3$,
the three-dimensional space. The covariant derivative is designated by
$\nabla_{\mu}$ and $\square\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the
d’Alembertian operator.

##  2 $f(R)$ gravity in brief: the state of the art

The modified Einstein-Hilbert (EH) action for $f(R)$ gravity [40, 41, 8, 9,
17, 18] is

$S=\int d^{4}x\sqrt{-g}\left[\frac{M^{2}_{\mbox{\footnotesize
P}}}{2}f(R)+\mathcal{L}_{M}\right],$ (2.1)

where $M_{\mbox{\footnotesize P}}^{2}\equiv(8\pi G)^{-1}$ is the Planck mass
and $\mathcal{L}_{M}$ is the matter Lagrangian describing the material
content. Varying the modified EH action with respect to the metric variable,
we get the following field equations

$f^{\prime}(R)R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}f(R)-\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\right)f^{\prime}(R)=\frac{T_{\mu\nu}}{M^{2}_{\mbox{\footnotesize
P}}}\,,$ (2.2)

and the trace equation

$Rf^{\prime}(R)-2f(R)+3\,\square
f^{\prime}(R)=\frac{T}{M^{2}_{\mbox{\footnotesize P}}}\,,$ (2.3)

where $f^{\prime}(R)\equiv df(R)/dR$ and $T\equiv T^{\mu}_{\mu}$ is the trace
of the stress-energy tensor. In the absence of matter, the exact solution of
eq. (2.2) is

$Rf^{\prime}(R)-2f(R)=0\,,$ (2.4)

where the positive real roots of this equation yield the well-known de Sitter
vacuum solutions. These solutions form the basis for the description of the
early and late acceleration phases of the universe. However, as pointed out in
[42], there are acceptable solutions if

$\left.\frac{f^{\prime}(R)}{f^{\prime\prime}(R)}\right|_{R=R_{*}}>R_{*}\,,$
(2.5)

where $f^{\prime\prime}(R)\equiv d^{2}f(R)/{dR}^{2}$, and $R_{*}$ is a
positive real root of eq. (2.4). In fact, the $f(R)$ function characterizing
the MG model is not arbitrary. Instead, it must satisfy a set of rules to
ensure both theoretical consistency and phenomenological viability of the
model.

  * (i)

be stable in the interval of $R$ of cosmological interest, i.e.,

$f^{\prime}(R)>0\,,\ f^{\prime\prime}(R)>0;$ (2.6)

  * (ii)

have a stable Newtonian limit, i.e.,

$\left|f(R)-R\right|\ll R\,,\ \left|f^{\prime}(R)-1\right|\ll 1\,,\
Rf^{\prime\prime}(R)\ll 1\,,$ (2.7)

for $R\gg R_{0}$, where $R_{0}$ is the curvature scalar today;

  * (iii)

be indistinguishable from GR at the current level of accuracy of laboratory
experiments and tests of the Solar System phenomena involving gravity;

  * (vi)

recover the GR in the high-curvature regime, i.e.,

$\lim_{R\gg R_{0}}f(R)=R-2\Lambda\,.$ (2.8)

This constraint also means that $f^{\prime}(\infty)=1$, implying that
$0<f^{\prime}(R)<1$ for all $R$; and

  * (v)

have a stable (or metastable) asymptotic de Sitter future.

In the first condition (i), $f^{\prime}(R)>0$, ensures that gravity be
attractive (i.e., the gravitons are not ghosts) [43, 44, 45, 46], while the
second, $f^{\prime\prime}(R)>0$, avoids the scalaron from becoming a tachyon,
in the high-curvature regime [47, 48, 49]. Next, given the success of the
Newtonian theory in explaining the observed non-homogeneities at small scale
and compact objects, it is necessary to impose that $f(R)$ model recovers this
theory when $R\gg R_{0}$ at these scales, as indicated at (ii). The GR theory
is very well-tested in the laboratory and in the Solar System, where no
significant deviations from the theory have been observed to date. Besides,
observational data from the cosmic microwave background (CMB) regarding
processes of the early universe strongly agree with the robust predictions of
the concordance flat-$\Lambda$CDM model, such as the big bang nucleosynthesis
(BBN), a primitive radiation-dominated age, and another middle age matter-
dominated. Furthermore, we need the $f(R)$ theory to recover the GR in the
weak-field regime and in the distant past, as determined at (iii) and (iv).
Finally, for a description of the current DE-dominated era, the $f(R)$ model
needs to have a stable (or metastable) de Sitter phase. As we shall see in the
next section, new constraints on the $f(R)$ will be placed exploring the
cosmological viability of the model.

##  3 $f(R)$ cosmology

In this section, we briefly discuss the cosmic dynamics in the $f(R)$
cosmology, finding the evolution equations in the background and perturbative
levels for a generic $f(R)$ model.

###  3.1 Cosmological background

In order to derive the dynamics of the cosmological background in the $f(R)$,
let’s first consider the FLRW metric, describing a statistically homogeneous
and isotropic universe, given by

${ds}^{2}=-{dt}^{2}+a^{2}(t)\left(\frac{{dr}^{2}}{1-Kr^{2}}+r^{2}{d\Omega}^{2}\right)\,,$
(3.1)

where $K=0,\pm 1$ and $a(t)$ is the scale factor. Assuming a flat spatial
section ($K=0$) and replacing eq. (3.1) into eq. (2.2), we get the Friedmann
equations

$3H^{2}=\frac{\rho_{M}+\rho_{DE}}{M_{P}^{2}},$ (3.2)
$2\dot{H}+3H^{2}=-\frac{P_{M}+P_{DE}}{M_{P}^{2}}\,,$ (3.3)

where the overdot denotes differentiation with respect to cosmic time, $t$,
and

$\frac{\rho_{DE}}{M^{2}_{\mbox{\footnotesize
P}}}\equiv\frac{1}{2}\left[Rf^{\prime}(R)-f(R)\right]-3H\dot{f}^{\prime}(R)+3H^{2}\left[1-f^{\prime}(R)\right],$
(3.4) $\frac{P_{DE}}{M^{2}_{\mbox{\footnotesize
P}}}\equiv-\frac{1}{2}\left[Rf^{\prime}(R)-f(R)\right]+\ddot{f}^{\prime}(R)+2H\dot{f}^{\prime}(R)-(2\dot{H}+3H^{2})\left[1-f^{\prime}(R)\right].$
(3.5)

Dividing eq. (3.5) by eq. (3.4) we get the equation of state for the DE
gravitational component in the $f(R)$ theory, $w_{DE}=P_{DE}/\rho_{DE}$, where
$w_{DE}$ is its state parameter; this DE component obeys the continuity
equation

$\dot{\rho}_{DE}+3H\left(\rho_{DE}+P_{DE}\right)=0\,.$ (3.6)

In particular, in the GR case, where $f(R)\equiv R-2\Lambda$, the eqs. (3.4)
and (3.5) reduce to the perfect fluid equations, with $P_{DE}=-\rho_{DE}$,
hence $w_{DE}\equiv w_{\Lambda}=-1$. Finally, the deceleration parameter is
given by

$q\left(t\right)\equiv-\frac{\ddot{a}}{aH^{2}}=-1-\frac{\dot{H}}{H^{2}}.$
(3.7)

The $f(R)$ is expected to describe the current accelerated epoch,
$q(t_{0})<0$, as well as reduce to GR at the high-curvature regime (e.g., as
in the solar system neighborhood), and for high-redshift data $z\gg 1$.

The effect of a general function $f(R)$ on cosmological dynamics can be
analyzed from a geometric point of view. It is useful to define a new set of
functions $(m,s)$ and to study the solutions of $m=m(s)$ on this plane. Using
the eqs. (3.2), (3.3), (3.4), and (3.5) we define

$m\equiv\frac{Rf^{\prime\prime}(R)}{f^{\prime}(R)}\,,\
s\equiv-\frac{Rf^{\prime}(R)}{f(R)}\,.$ (3.8)

According to [10], some rules111since the newly established rules are
constraints imposed on $f(R)$, we adopt the numbering sequence, as presented
in section 2. are established for a viable cosmic dynamics, namely:

  * (vi)

A matter-dominated middle epoch, necessary for structure formation, is
achieved if

$m(s)\approx+0\,,\ \frac{dm(s)}{ds}>-1\,,$ (3.9)

at $s=-1$; and

  * (vii)

The matter-dominated phase, will be followed by a late accelerated phase, only
if

$0\leq m(s)\leq 1\,,$ (3.10)

at $s=-2$, or

$m(s)=-s-1\,,\ \frac{\sqrt{3}-1}{2}<m(s)\leq 1\,,\ \frac{dm(s)}{ds}<-1\,.$
(3.11)

The $\Lambda$CDM model correspond to the straight line $m=0,\,\forall s$, on
the plane $(m,s)$. Since the viable $f(R)$ models are those that approach the
$\Lambda$CDM at high-curvature (i.e., when $R\gg R_{0}$), $m(s)$ needs to be
close to zero during the matter-driven epoch, as indicated in the first of the
conditions in the eq. (3.9). The second condition, in turn, is required to
connect the matter-dominated epoch and the current DE-dominated epoch. On the
other hand, the constraints in the eq. (3.10) indicate a purely de Sitter
vacuum, $w_{DE}=-1$, while the eqs. (3.11) corresponds to a non-phantom
attractor, $w_{DE}>-1$. Regarding the conditions (3.10) and (3.11), the $f(R)$
model needs to satisfy only one of them.

###  3.2 Cosmological perturbations

Perturbations around the FLRW background produced in the primordial universe
explain the relevant part of the CMB power spectrum [50, 51, 52]. When
crossing the horizon at the inflationary period, such fluctuations were
frozen, becoming the primordial seeds to the growth of the structures on large
scales after their re-enter. At the linear level, and using the Newtonian
gauge, the line element can be written as [53, 54]

${ds}^{2}=-\left(1+2\Psi\right){dt}^{2}+a^{2}\left(1-2\Phi\right)\delta_{ij}{dx}^{i}{dx}^{j},$
(3.12)

where

$\delta_{ij}=\left\\{\begin{array}[]{ll}1,&\text{for}\ i=j,\\\ 0,&\text{for}\
i\neq j,\end{array}\right.$ (3.13)

is the Kronecker delta function, and $|\Psi|,|\Phi|\ll 1$ are small functions
of the three-space, $x^{i}$, and cosmic time, $t$, termed the Bardeen
potentials [55].

The space-time metric fluctuations give rise to fluctuations in the material
content that evolve via gravitational instability [54, 55, 56, 57], and can be
described by the perturbed components of the stress-energy tensor,

$\displaystyle T^{0}_{0}$ $\displaystyle=$
$\displaystyle-\left(1+\delta\right)\,,$ (3.14) $\displaystyle T^{0}_{j}$
$\displaystyle=$ $\displaystyle-\left(\rho+P\right)v^{i}\,,$ (3.15)
$\displaystyle T^{i}_{0}$ $\displaystyle=$
$\displaystyle\left(\rho+P\right)v_{j}\,,$ (3.16) $\displaystyle T^{i}_{j}$
$\displaystyle=$ $\displaystyle\left(P+\delta
P\right){\delta}^{i}_{j}+P{\Pi}^{i}_{j}\,,$ (3.17)

where $\delta\equiv\delta\rho/\rho$ is the matter density contrast, that, like
the four-velocities ($v^{i}$), pressure fluctuations ($\delta P$) and the
anisotropic stresses, depend on the space-time coordinates. In the absence of
anisotropic stresses, ${\Pi}^{i}_{j}=0$, one has $\Psi=\Phi$ [54].

For small scale fluctuations ($k^{2}\gg a^{2}H^{2}$) of a matter-dominated
fluid one arrives to the equation [58]

$\ddot{\delta}_{m}+2H\dot{\delta}_{m}-4\pi G_{eff}\rho_{m}\delta_{m}=0\,,$
(3.18)

where

$G_{eff}\left(a,k\right)\equiv\frac{G}{f^{\prime}(R)}\frac{1+4\left(\frac{k^{2}}{a^{2}R}\right)m}{1+3\left(\frac{k^{2}}{a^{2}R}\right)m}\,,$
(3.19)

is called the effective gravitational constant, that depends on the scale
factor $a$ and on the scale $k$.

An observable quantity useful to compare a model with current observations is
the growth rate of cosmic structures (or simply growth function), $f_{g}(a)$,
defined by [59]

$f_{g}(a)\equiv\frac{d\ln{\delta_{m}\left(a\right)}}{d\ln{a}}\,.$ (3.20)

A good approximation for the growth function is given by
$f_{g}(z)\simeq\Omega_{m}^{\gamma}(z)$, where $\gamma(z)$ is the growth index
[60]. In DE models based on GR theory $\gamma(z)$ is approximated by the
constant value $\gamma\simeq 3\,(1-w_{DE})/(5-6w_{DE})$ [61]. Thus, the
$\Lambda$CDM model corresponds to $\gamma\simeq 6/11$.

The matter fluctuations amplitude, $\sigma_{\mathcal{R}}(a)$, in turn, are
directly observed from the CMB power spectrum and is given by [62]

$\sigma_{\mathcal{R}}\left(a\right)=\sigma_{\mathcal{R},0}\left[\frac{\delta_{m}\left(a\right)}{\bar{\delta}_{m}\left(1\right)}\right]\,,$
(3.21)

where $\sigma_{\mathcal{R},0}$ corresponds to the amplitude at $z=0$ where
$\mathcal{R}$ represents a physical scale. The constant $\bar{\delta}_{m}(1)$
is a Planck normalization factor so that at $z=0$, in the $\Lambda$CDM model,
$\sigma_{\mathcal{R}}(0)\equiv\sigma_{\mathcal{R},0}$. It is common to perform
the measurements of several cosmological tracers at the physical scale222the
Hubble constant $H_{0}$ is usually written as $H_{0}\equiv 100\,h$ km/s/Mpc,
where $h$ is a real adimensional number $\mathcal{R}=8$ Mpc/h. Thus the
product between $f_{g}(a)$ and $\sigma_{8}(a)$ results in

$[f{\sigma_{8}}](a)=\frac{\sigma_{8,0}}{\bar{\delta}_{m}(1)}\left[\frac{d\,\delta_{m}(a)}{d\ln{a}}\right]\,,$
(3.22)

which measures the matter density perturbation rate at the physical scale of
$8$ Mpc$/h$. This combination is used more often than simply $f_{g}(a)$ to
derive constraints for theoretical model parameters because of data
availability.

##  4 $R^{2}$-corrected AB model

The first $f(R)$ model proposed by Appleby and Battye in [22] is a two-
parameter model:

$\displaystyle
f^{\,I}_{AB}(R)\equiv\frac{R}{2}+\frac{\epsilon_{AB}}{2}\ln{\left[\frac{\cosh{\left(\frac{R}{\epsilon_{AB}}-b\right)}}{\cosh
b}\right]}\,,$ (4.1)

that at large $R$, mimic the GR theory with a non-true cosmological constant.
The free model parameters $b$ and $\epsilon_{AB}$, are related by

$\displaystyle\epsilon_{AB}\equiv\frac{R_{vac}}{b+\ln\,(\,2\cosh b\,)}\,,$
(4.2)

where $R_{vac}\equiv 12H_{0}^{2}>\epsilon_{AB}$ is the current vacuum scalar
curvature (roughly speaking, $R_{vac}\sim R_{0}$). The eq. (4.2) relating $b$
and $\epsilon_{AB}$ tells us that this model has only one free parameter $b$.

It is important to mention that some $f(R)$ models well established in the
literature, like Hu-Sawicki [21] and Starobinsky [20] models, present three
parameters.

It has been shown in [37, 38] that a sudden weak curvature singularity forms
generically when the $f^{\prime\prime}(R)$ becomes zero for some finite value
of $R$ so that the condition $f^{\prime\prime}(R)>0$ is marginally violated,
leading to two more elementary problems: both an unbounded growth of the
scalaron mass and an undesired overabundance of this particle at high-
curvatures. In other words, it means that the scalaron can behave like a
tachyon and that the amplitude of the oscillations in $R$ grows indefinitely
when $R\gg R_{vac}$. For the AB model in eq. (4.1) and many others discussed
in the literature, this occurs for $R$ within the range of cosmological
interest, pointing out an incompleteness of these models. However, it was
observed that the simple addition of a term proportional to $R^{2}$ in the
$f(R)$, with a sufficiently small coefficient to ensure the existence of the
primordial inflation, solves this type of singularity [37, 39].

In this way it was proposed the $R^{2}$-corrected Appleby-Battye, or
$R^{2}$-AB, model [37]

$f_{AB}(R)\equiv f^{\,I}_{AB}(R)+\frac{R^{2}}{6M^{2}}\,,$ (4.3)

where $b\geq 1.6$ is required for the $R^{2}$-AB model to reproduce the recent
accelerated phase, and $M$ characterizes a scale mass coinciding with the
scalaron rest-mass whenever low curvature modifications to GR can be
neglected.

The scalaron mass, $M_{s}(R)$, is given by

$M_{s}^{2}\left(R\right)\equiv\frac{f^{\prime}\left(R\right)-Rf^{\prime\prime}\left(R\right)}{3f^{\prime\prime}\left(R\right)},$
(4.4)

so that, at high-curvature regime ($R\gg R_{0}$), we must have
$M_{s}=1/\sqrt{3f^{\prime\prime}(R)}$ [17]. Then, one can verify that the
upper limit for the scalaron mass, $M_{s}\leq M$, was made possible by adding
the Starobinsky-like term $R^{2}/6M^{2}$.

Assuming that the universe is filled with a fluid, consisting of ordinary
matter and radiation, represented by subscript index $m$ and $r$,
respectively, and a dark energy component, by subscript index DE, the
background dynamical equations are given by:

$3H^{2}=\frac{\rho_{m}+\rho_{r}+\rho_{DE}}{M_{P}^{2}}\,,$ (4.5)
$2\dot{H}+3H^{2}=-\frac{P_{m}+P_{r}+P_{DE}}{M_{P}^{2}}\,.$ (4.6)

Previously, it was shown by eqs. (3.4) and (3.5), that the dark energy
component has a purely geometrical origin. For the model (4.3):

$\displaystyle\frac{\rho_{DE}}{M_{P}^{2}}=-\left(\frac{1}{3M^{2}}+\frac{\text{sech}^{2}{x}}{2\epsilon_{AB}}\right)3H\dot{R}$
$\displaystyle-\frac{\epsilon_{AB}}{4}\ln{\left(\frac{\cosh{x}}{\cosh{b}}\right)}+\frac{R}{4}-\frac{R^{2}}{12M^{2}}$
(4.7)
$\displaystyle+(\dot{H}+H^{2})\left[\frac{R}{M^{2}}+\frac{3(\tanh{x}-1)}{2}\right]\,,$
$\displaystyle\frac{P_{DE}}{M_{P}^{2}}=\left(\frac{1}{3M^{2}}+\frac{\text{sech}^{2}{x}}{2\epsilon_{AB}}\right)\xi_{R}$
$\displaystyle+\frac{\epsilon_{AB}}{4}\ln{\left(\frac{\cosh{x}}{\cosh{b}}\right)}-\frac{R}{4}+\frac{R^{2}}{12M^{2}}$
(4.8)
$\displaystyle-\left(\frac{R}{3M^{2}}+\frac{\tanh{x}-1}{2}\right)\xi_{H}-\frac{\dot{R}^{2}\tanh{x}\,\text{sech}^{2}{x}}{\epsilon_{AB}^{2}}\,,$

where we define the auxiliary variables $x\equiv R/\epsilon_{AB}-b$,
$\xi_{R}\equiv\ddot{R}+2H\dot{R}$, and $\xi_{H}\equiv\dot{H}+3H^{2}$.
Throughout the domain eras of radiation and early matter, when $R_{vac}\ll
R\ll M^{2}$, the expressions in (4.7) and (4.8) can be approximated,
respectively, by

$\frac{\rho_{DE}}{M_{P}^{2}}\simeq\frac{R_{vac}}{4}+\frac{1}{M^{2}}\left[(\dot{H}+H^{2})R-H\dot{R}-\frac{R}{12}\right]\\\
-\frac{1}{e^{2x}}\left[\frac{6H\dot{R}}{\epsilon_{AB}}+3(\dot{H}+H^{2})+\frac{\epsilon_{AB}}{4}\right],$
(4.9)
$\frac{P_{DE}}{M_{P}^{2}}\simeq-\frac{R_{vac}}{4}+\frac{1}{M^{2}}\left[\frac{(\xi_{R}-\xi_{H}R)}{3}+\frac{R^{2}}{12}\right]-\frac{1}{e^{2x}}\left[\frac{2\xi_{R}}{\epsilon_{AB}}+\xi_{H}-\left(\frac{2\dot{R}}{\epsilon_{AB}}\right)^{2}+\frac{\epsilon_{AB}}{4}\right].$
(4.10)

Since $R$ and $M$ are both very large in these epochs, we can expect a strong
suppression of the last two terms in eqs. (4.9) and (4.10), and, consequently,
to obtain the state equation $w_{DE}\simeq-1$. It means that $R^{2}$-AB model
recovers GR at high-$z$. On the other hand, at later times, when $R\sim
R_{vac}$, the scalaron mass becomes small and significant deviations from
$w_{DE}=-1$ must be expected.

Figure 1: Plot of $H(z)$ functions for the $R^{2}$-AB and flat-$\Lambda$CDM
models. We have used the Planck Collaboration best-fit data: $H_{0}=67.4$
km/s/Mpc and $\Omega_{m,0}=0.315$ [63]; in addition we set $b=2$ and a scale
mass such that $\mu=10^{-7}$. Figure 2: Comparison plot of $w_{eff}(z)$ from
the $R^{2}$-AB model and $w_{\Lambda}=-1$ from fiducial $\Lambda$CDM. We have
used $H_{0}=67.4$ km/s/Mpc, $\Omega_{m,0}=0.315$, $b=2$, and $\mu=10^{-7}$.

We now discuss the initial conditions and the most convenient mass scale for
the curvature range. Firstly, we fix $z_{i}=4$ as the initial redshift in our
analyses, because we can safely assume that the $\Lambda$CDM describes well
the observed universe at this epoch. In this way we will have the three
necessary initial values $H_{i}\equiv H_{GR}(t_{i})$,
$\dot{H}_{i}\equiv\dot{H}_{GR}(t_{i})$ and
$\ddot{H}_{i}\equiv\ddot{H}_{GR}(t_{i})$, where $t_{i}\equiv t(z_{i})$. Then,
for growing modes in the perturbation level, we assume initial conditions such
that $\delta_{m}(a_{i})=a_{i}$ and $d\,\delta_{m}(a_{i})/da=1$ for $a_{i}\ll
1$ [62]. As suggested in [37], we can set a mass scale such that
$\mu\equiv\epsilon_{AB}/M^{2}=10^{-7}$.

Figure 3: Evolution of the matter density fluctuations, $\delta_{m}(z)$, for
both $R^{2}$-AB and flat-$\Lambda$CDM models at the scale $k=0.125$ Mpc-1. We
set $\Omega_{m,0}=0.315$, $b=2$, and $\mu=10^{-7}$. It is worth mentioning
that $\delta_{m}(z)$ is independent of $H_{0}$. Figure 4: Evolution of
$[f\sigma_{8}](z)$ for $R^{2}$-AB and flat-$\Lambda$CDM models at the scale
$k=0.125$ Mpc-1. In this plot we have used $\Omega_{m,0}=0.315$, $b=2$,
$\sigma_{8,0}=0.811$, and $\mu=10^{-7}$. Like $\delta_{m}(z)$,
$[f\sigma_{8}](z)$ is also independent of $H_{0}$.

In figure 1, we present the $H(z)$ function , given by eqs. (4.5) and (4.6),
for the $R^{2}$-AB model. Additionally, for comparison purposes, we include
the fiducial $H(z)$ from the $\Lambda$CDM. Afterward, in figure 2, we show the
behavior of the state parameter, $w_{DE}(z)$, where, it is possible to note
that the $R^{2}$-AB model does not reproduces the standard $\Lambda$CDM
history around the interval $z<2$. This situation changes very little as we
vary the parameters involved $H_{0}$, $\Omega_{m,0}$, and $b$. We can also
observe that as we increase the value of $b$, the solution from $f(R)$
approaches that one of $\Lambda$CDM model (see appendix A). In order to
understand this behavior, let us consider the limit case $b\gg 1$, or
equivalently, $x\gg 1$. In this case, $R^{2}$-AB model defined in eq. (4.3)
can be expanded as

$f(R)\simeq R-\frac{R_{vac}}{2}+\Gamma(R)\,,$ (4.11)

where $\Gamma(R)\equiv(\epsilon_{AB}/2)e^{-2x}+R^{2}/6M^{2}$ and $\Gamma(R)\ll
R$. Then, for $b$ sufficiently large, the eq. (4.11) can be written as
$f(R)\simeq R-2\Lambda$, where $\Lambda\equiv R_{vac}/4$. In this context,
$\Lambda$ is said to be a non-true cosmological constant, as it assumes a
geometric (i.e., gravitational) role rather than be part of the material
content. Another possibility to recover RG is reached by setting $R\gg
R_{vac}$, as mentioned in eq. (2.8).

Due to degeneracy at the background level, we must look for new cosmological
tracers to obtain the best fit for the model parameters and their
uncertainties. Our first choice is to consider the cosmological perturbations
through the matter density contrast, $\delta_{m}(z)$, and the parametrized
growth rate of cosmic structures, $[f\sigma_{8}](z)$. Following this goal, we
obtain the eq. (3.18) by solving eq. (3.22) for both $f(R)$ and
flat-$\Lambda$CDM (for which $G_{eff}\equiv G$) models. It is worth mentioning
that in MG the structure formation depends on scale $k$ through the effective
gravitational constant, $G_{eff}(z,k)$ (see eq. (3.19)). The plots of
$\delta_{m}(z)$ and $[f\sigma_{8}](z)$, shown in figures 3 and 4, were
obtained by assuming $k=0.125$ Mpc-1 and the Planck Collaboration best-fit
$\sigma_{8,0}=0.811$ [63].

As in the cosmological background, increasing the model parameter $b$ causes
the blue curves to overlap with the red ones, pointing out that this parameter
provides a measurement of the similarity (or difference) between GR and
$R^{2}$-AB model. Thus, the $f(R)$ model (4.3) recovers the GR whenever $R\gg
R_{vac}$ or $b\gg 1$.

##  5 Cosmological datasets

One powerful technique to measure $H(z)$, independently of the assumption of a
cosmological model, is the cosmic chronometers (CC) method. This approach is
based on the relationship

$H\left(z\right)=-\frac{1}{\left(1+z\right)}\frac{dz}{dt}\simeq-\frac{1}{\left(1+z\right)}\frac{\Delta
z}{\Delta t}\,,$ (5.1)

obtained from the definition $H\equiv\dot{a}/a$, where the derivative term,
$dz/dt\simeq\Delta z/\Delta t$, can be determined from two passively-evolving
galaxies, i.e., galaxies with old stellar populations and low star formation
rates, whose redshifts are slightly different and whose ages are well-known.
Furthermore, the chosen galaxies must have an age difference much smaller than
their passively-evolving time [71]. Of course, to estimate the age of galaxies
it is necessary to assume a stellar population synthesis (SPS) model. In table
1 we list $31$ measurements on $H(z)$ obtained with the CC methodology, where
the age of galaxies was obtained assuming BC$03$ SPS model [72]. Thus, these
measurements contain systematic uncertainties related only to SPS model and to
possible contamination due to the presence of young stars in quiescent
galaxies [73, 74].

Table 1: The $31$ measurements on $H(z)$ obtained from the CC approach using the BC$03$ SPS model. $z$ | $H(z)$ | Ref. | $z$ | $H(z)$ | Ref.
---|---|---|---|---|---
$0.07$ | $69.0\pm 19.6$ | [64] | $0.4783$ | $80.9\pm 9.0$ | [67]
$0.09$ | $69.0\pm 12.0$ | [65] | $0.48$ | $97.0\pm 62.0$ | [69]
$0.12$ | $68.6\pm 26.2$ | [64] | $0.593$ | $104.0\pm 13.0$ | [66]
$0.17$ | $83.0\pm 8.0$ | [65] | $0.68$ | $92.0\pm 8.0$ | [66]
$0.179$ | $75.0\pm 4.0$ | [66] | $0.781$ | $105.0\pm 12.0$ | [66]
$0.199$ | $75.0\pm 5.0$ | [66] | $0.875$ | $125.0\pm 17.0$ | [66]
$0.2$ | $72.9\pm 29.6$ | [64] | $0.88$ | $90.0\pm 40.0$ | [69]
$0.27$ | $77.0\pm 14.0$ | [65] | $0.9$ | $117.0\pm 23.0$ | [65]
$0.28$ | $88.8\pm 36.6$ | [64] | $1.037$ | $154.0\pm 20.0$ | [66]
$0.352$ | $83.0\pm 14.0$ | [66] | $1.3$ | $168.0\pm 17.0$ | [65]
$0.3802$ | $83.0\pm 13.5$ | [67] | $1.363$ | $160.0\pm 33.6$ | [70]
$0.4$ | $95.0\pm 17.0$ | [65] | $1.43$ | $177.0\pm 18.0$ | [65]
$0.4004$ | $77.0\pm 10.2$ | [67] | $1.53$ | $140.0\pm 14.0$ | [65]
$0.4247$ | $87.1\pm 11.2$ | [67] | $1.75$ | $202.0\pm 40.0$ | [65]
$0.4497$ | $92.8\pm 12.9$ | [67] | $1.965$ | $186.5\pm 50.4$ | [70]
$0.47$ | $89.0\pm 49.6$ | [68] |  |  |

The most common approach to study the clustering evolution of cosmic
structures is through the normalized growth rate, $f\sigma_{8}$, given in eq.
(3.22). A possible way to obtain this one is by first measuring the velocity
scale parameter, given by $v_{\beta}=f_{g}/\beta$, where $\beta$ is the bias
factor, and then use it in the relation

$[f\sigma_{8}](z)=v_{\beta}(z)\,\sigma_{8}^{tr}(z)\,,$ (5.2)

where $\sigma_{8}^{tr}(z)$ is the matter fluctuation amplitude of the
cosmological tracer, e.g., HI line extra-galactic sources (EGS), luminous red
galaxies (LRG), quasars (QSOs), Supernovae Ia (SNe-Ia) and emission-line
galaxies (ELG).

The parametrized growth rate data, $[f\sigma_{8}](z)$, are most often obtained
using the redshift-space distortions (RSD) effect observed in galaxy surveys
[76, 77, 78, 79, 80, 81, 89, 90, 35]. In table 2 we can see the data
compilation from [75] considering measurements of $[f\sigma_{8}](z)$. This
compilation follows a methodology, in which double counting is eliminated and
possible biases are reduced, thus ensuring the reliability of the dataset (see
section $3$ of [75]).

Table 2: The $20$ measurements on $[f\sigma_{8}](z)$, obtained from several surveys and cosmological tracers, compiled by [75]. $z$ | $[f\sigma_{8}](z)$ | Ref. | Survey | Cosmological Tracer
---|---|---|---|---
$0.02$ | $0.398\pm 0.065$ | [76] | SNe-Ia $+$ IRAS | SNe-Ia $+$ galaxies
$0.025$ | $0.39\pm 0.11$ | [77] |  $6$dFGS | Voids
$0.067$ | $0.423\pm 0.055$ | [78] |  $6$dFGS | Galaxies
$0.10$ | $0.37\pm 0.13$ | [79] | SDSS-veloc | DR7 Galaxies
$0.15$ | $0.53\pm 0.16$ | [80] | SDSS-IV | eBOSS DR16 MGS
$0.32$ | $0.384\pm 0.095$ | [81] | BOSS-LOWZ | DR10 and DR11
$0.38$ | $0.497\pm 0.045$ | [80] | SDSS-IV | eBOSS DR16 galaxies
$0.44$ | $0.413\pm 0.08$ | [82] | WiggleZ | Bright ELG
$0.57$ | $0.453\pm 0.022$ | [83] | CMASS-BOSS | DR12 voids $+$ galaxies
$0.59$ | $0.488\pm 0.060$ | [84] | SDSS-CMASS | DR12
$0.70$ | $0.473\pm 0.041$ | [80] | SDSS-IV | eBOSS DR16 LRG
$0.73$ | $0.437\pm 0.072$ | [82] | WiggleZ | Bright ELG
$0.74$ | $0.50\pm 0.11$ | [85] | SDSS-IV | eBOSS DR16 voids
$0.76$ | $0.440\pm 0.040$ | [86] | VIPERS v$7$ | Galaxies
$0.85$ | $0.52\pm 0.10$ | [85] | SDSS-IV | eBOSS DR16 voids
$0.978$ | $0.379\pm 0.176$ | [87] | SDSS-IV | eBOSS DR14 QSOs
$1.05$ | $0.280\pm 0.080$ | [86] | VIPERS v$7$ | Galaxies
$1.40$ | $0.482\pm 0.116$ | [88] | FastSound | ELG
$1.48$ | $0.30\pm 0.13$ | [85] | SDSS-IV | eBOSS DR16 voids
$1.944$ | $0.364\pm 0.106$ | [87] | SDSS-IV | eBOSS DR14 QSOs

##  6 Analyses and Results

A very powerful technique for determining the most likely parameter set of a
model is the Markov-Chain Monte Carlo (MCMC) method. The general goal of the
MCMC method is to sample the posterior probability distribution function (PDF)
around their most likely values, which allows us to obtain the best-fit model
parameters with robust uncertainties.

Figure 5: MCMC analyses for the $R^{2}$-AB model considering CC $H(z)$ data,
with the priors $H_{0}=[50,\,80]$, $\Omega_{m,0}=[0.1,\,0.6]$, and
$b=[1.6,\,6]$.

The MCMC is a Bayesian statistical procedure based on the Bayes theorem [91]

$P(\vartheta_{i}|D_{i},\alpha)=\frac{P(D_{i}|\vartheta_{i},\alpha)P(\vartheta_{i}|\alpha)}{P(D_{i}|\alpha)}\,,$
(6.1)

where $P(\vartheta_{i}|D_{i},\alpha)$ is the posterior PDF,
$P(D_{i}|\vartheta_{i},\alpha)$ is the likelihood, $P(\vartheta_{i}|\alpha)$
is called prior, and $P(D_{i}|\alpha)$ is the evidence; $\vartheta_{i}$
represent the model parameters set, $D_{i}$ is the dataset and $\alpha$
denotes the prior information (model). Since the evidence, $P(D_{i}|\alpha)$,
is independent on model, we can ignore it as a normalizing constant.

Assuming that all observations $D_{i}$ are independent, the form which best
describes the likelihood is the Gaussian distribution

$\mathcal{L}(\vartheta_{j})=\exp{\left\\{-\frac{1}{2}\sum_{i=1}^{n}\left[\frac{E(\vartheta_{j}|\alpha)-D_{i}}{\sigma_{i}}\right]^{2}\right\\}}\,,$
(6.2)

where $E(\vartheta_{j}|\alpha)$ is the expected value, based on a model, and
$\sigma_{i}$ is the set of errors associated with $D_{i}$. Here we disregard
the constant $1/\sqrt{2\pi}\sigma_{i}$, as it is noted that there is no
significant difference in the result when we consider it. Next, the summation
in the exponential argument is defined as chi-square ($\chi^{2}$), simplifying
the likelihood as $\mathcal{L}(\vartheta_{j})\equiv
P(D_{i}|\vartheta_{j},\alpha)=\exp{(-\chi^{2}/2)}$. Additionally, it is usual
to consider that the prior sets have the same probability of occurrence, so
that $P(\vartheta_{i}|\alpha)$ has the form of the Dirac delta distribution,

$P(\vartheta_{i}|\alpha)=\left\\{\begin{array}[]{ll}1&\text{if}\quad\vartheta_{i}^{(0)}<\vartheta_{i}<\vartheta_{i}^{(1)},\\\
0&\text{otherwise},\end{array}\right.$ (6.3)

where $\vartheta_{i}^{(0)}$ and $\vartheta_{i}^{(1)}$ are the flat prior
intervals. Hence the posterior PDF is given by
$P(\vartheta_{j}|D_{i},\alpha)=\mathcal{L}(\vartheta_{j})$, i.e, only by the
likelihood.

There exist two ways in order to estimate the optimum value for the model
parameters: maximum likelihood and least chi-square. The latter is more common
and can be done using a Python library (e.g., scipy.optimize). Next, we shall
explore the parametric space of the model parameters, sampling the posterior
distribution around that value, following the MCMC method and the Metropolis-
Hastings algorithm. The confidence regions has been drawn assuming
$\mathcal{L}=\mathcal{L}_{max}+{\Delta\mathcal{L}}_{0}$, where the constant
${\Delta\mathcal{L}}_{0}$ is determined by the cumulative probability density.
To implement the MCMC routine we use Python as well.

By defining the Hubble parameter as $H(a)\equiv H_{0}\,\varepsilon(a)$ it can
be shown that $f\sigma_{8}$ is independent on $H_{0}$, since the term
$G_{eff}/G$ is $H_{0}$ independent as well. In our analyses we set
$\vartheta_{i}=[H_{0},\,\Omega_{m,0},\,b]$, for the Hubble function, and
$\vartheta_{i}=[\Omega_{m,0},\,\sigma_{8,0},\,b]$, for the growth rate, where
we have used the priors $H_{0}=[50,\,80]$, $\Omega_{m,0}=[0.1,\,0.6]$,
$\sigma_{8,0}=[0.7,\,1]$ and $b=[1.6,\,6]$.

In figures 5 and 6 we show the results of our MCMC analyses for $R^{2}$-AB
model, considering $H(z)$ and $[f\sigma_{8}](z)$ datasets, respectively, and
we summarize them in table 3. As seen in this table, the most likely values
obtained for the model parameters are: $H_{0}=66.6^{+3.3}_{-4.4}$ kms-1Mpc-1,
$\Omega_{m,0}=0.324^{+0.098}_{-0.052}$ and $b=1.6^{+3.1}_{-0.0}$, from the
$H(z)$ analysis, and $\Omega_{m,0}=0.304^{+0.080}_{-0.067}$,
$\sigma_{8,0}=0.791\pm 0.040$ and $b=1.76^{+2.91}_{-0.15}$, from the
$[f\sigma_{8}](z)$ analysis.

Figure 6: MCMC analyses for the $R^{2}$-AB model considering $[f\sigma_{8}](z)$ data from several cosmic tracers and the priors $\Omega_{m,0}=[0.1,\,0.6]$, $\sigma_{8,0}=[0.7,\,1]$, and $b=[1.6,\,6]$. Table 3: Best-fit values from our MCMC likelihood analyses. Dataset | $H_{0}$ | $\Omega_{m,0}$ | $\sigma_{8,0}$ | $b$
---|---|---|---|---
$H(z)$ | $66.6^{+3.3}_{-4.4}$ | $0.324^{+0.098}_{-0.052}$ | — | $1.60^{+3.1}_{-0.0}$
$[f\sigma_{8}](z)$ | — | $0.304^{+0.080}_{-0.067}$ | $0.791\pm 0.040$ | $1.76^{+2.91}_{-0.15}$

##  7 Conclusions

The flat-$\Lambda$CDM, with its mysterious dark energy component in the form
of a cosmological constant, is not the final model of cosmology. Efforts are
being devoted to study alternative cosmological scenarios where GR theory is
modified in a way to explain the observed universe, both at the background and
perturbative levels, but having GR as a suitable limit.

One class of promising candidates to explain the accelerated expansion of the
universe is based on an interpretation of gravity different from that provided
by GR, a theoretical approach known by the generic name of modified gravity
theory (MG theory). This new geometric scenario for the space-time has to
satisfy phenomenological rigorous criteria [40, 41, 8, 9]:

  * •

because GR is a well-established theory for strong gravitational fields and
small scales, any attempt to modify it should contain GR as a limiting theory
at suitable scales and strong gravitational fields;

  * •

in the distant past, $z\gg 1$, the MG theory should have a behavior concordant
with a matter dominated era;

  * •

at large scales and from a recent past, $z<z_{t}$ (= transition redshift), the
behavior expected is such that the MG theory explains the accelerated
expansion phase of the universe, a feature well established by different
cosmological tracers (background level);

  * •

at the perturbative level, the MG theory should satisfactorily explain the
growth rate of cosmic structures data.

This cosmic phenomenology, expected to be satisfied by any contender of the
concordance model of cosmology, the $\Lambda$CDM model, makes non-trivial the
search for good $f(R)$ candidates. One of such good candidates could be the
$R^{2}$-AB model, here investigated using cosmological data to constrain its
model parameters [22, 37, 38].

In general, an $f(R)$ –or other alternative– model undergo the problem that
having a large number of parameters makes the statistical best-fit process be
less efficient than that one made with the flat-$\Lambda$CDM with just one
free parameter. For this our interest here in studying the corrected Appleby-
Battye model, the $R^{2}$-AB model, a model with 2 parameters and 1
constraining relation, which determines a model with one free parameter: $b$.
In this way, we performed MCMC analyses of the $R^{2}$-AB model and observed
that, for both $H(z)$ and $[f\sigma_{8}](z)$ studies, the values obtained for
$H_{0}$, $\Omega_{m,0}$, and $\sigma_{8,0}$ were fully concordant with the
flat-$\Lambda$CDM cosmological parameters obtained by the Planck Collaboration
[63].

Moreover, our MCMC statistical analyses of the section 6 show that the
$R^{2}$-AB model is reasonably well constrained by the cosmological data
applied: (i) at the background level the $H(z)$ data returns a best-fit value
for the model parameter $b=1.6^{+3.1}_{-0.0}$, and (ii) at the perturbative
level the $[f\sigma_{8}](z)$ measurements returns a best-fit value
$b=1.76^{+2.91}_{-0.15}$, in agreement with the value obtained at the
background level, with both values of the $b$ parameter within the interval
where the $R^{2}$-AB model satisfies all the phenomenological criteria
mentioned above [10, 22, 37].

Our analyses reveal that the $R^{2}$-AB model is in good agreement with
observational data, both at the background and perturbative levels, for this
we conclude that it is an alternative cosmological model that provides a
geometrical explanation for the universe accelerated expansion, not being
necessary to assume a dark energy component in the universe.

Figure 7: Behavior of the Hubble function, $H(z)$, from the $R^{2}$-AB model
for different values of the parameter $b$, keeping fixed $H_{0}=67.4$ km/s/Mpc
and $\Omega_{m,0}=0.315$. The dashed red curve corresponds to $H(z)$ from the
flat-$\Lambda$CDM model. Figure 8: Equation of state from the $R^{2}$-AB model
for different values of the parameter $b$, keeping fixed $H_{0}=67.4$ km/s/Mpc
and $\Omega_{m,0}=0.315$, for comparison we plot the case $w_{\Lambda}=-1$ for
the flat-$\Lambda$CDM model. Figure 9: $\sigma_{8}(z)$ function, at the scale
$k=0.125$ Mpc-1, from the $R^{2}$-AB model considering several values of $b$,
and for comparison also plot the flat-$\Lambda$CDM model case, where we have
used $\Omega_{m,0}=0.315$ and $\sigma_{8,0}=0.811$. Figure 10:
$[f\sigma_{8}](z)$ function, at scale of $k=0.125$ Mpc-1, from the $R^{2}$-AB
model considering illustrative values of $b$, and for comparison also plot the
standard flat-$\Lambda$CDM case, where we have used $\Omega_{m,0}=0.315$ and
$\sigma_{8,0}=0.811$.

##  Appendix A  Effect of different $b$ values on cosmic observables

We find interesting to compare the evolution of some observables, like
$H,\,w_{eff},\,\sigma_{8}$, and $[f\sigma_{8}]$, for different values of the
model parameter $b$. For this, in figures 7, 8, 9 and 10 we show these plots
to observe the effect of $b$ on some cosmological observables obtained in the
$R^{2}$-AB model. These figures illustrate this comparative analysis, where we
are keeping the values of the other parameters, $H_{0}$, $\Omega_{m,0}$,
$\sigma_{8,0}$, fixed; moreover, we keep the scale mass as $\mu=10^{-7}$.

## Acknowledgments

BR and AB would like to thank the Brazilian Agencies CAPES and CNPq for their
respective fellowships. MC would like to acknowledges the Observatório
Nacional for the hospitality.

## References

  * [1] A. G. Riess et. al., _Observational evidence from supernovae for an accelerating universe and a cosmological constant_ , _AJ_ 116 (1998) 1009 [astro-ph/9805201].
  * [2] S. Perlmutter et. al., _Measurements of $\Omega$ and $\Lambda$ from 42 high-redshift supernovae_, _AJ_ 517 (1999) 565 [astro-ph/9812133].
  * [3] A. G. Riess, L. Breuval, W. Yuan, S. Casertano, L. M. Macri, D. Scolnic, T. Cantat-Gaudin, R I. Anderson and M. Cruz Reyes, _Cluster cepheids with high precision Gaia parallaxes, low zero-point uncertainties, and Hubble Space Telescope photometry_ , _AJ_ 938 (2022) 36 [arXiv:2208.01045].
  * [4] L. A. Anchordoqui, E. Di Valentino, S. Pan and W. Yang, _Dissecting the $H_{0}$ and $S_{8}$ tensions with Planck + BAO + supernova type Ia in multi-parameter cosmologies_, _JHEAp_ 32 (2021) 28 [arXiv:2107.13932].
  * [5] W. Yang, E. Di Valentino, S. Pan, A. Shafieloo and X. Li, _Generalized emergent dark energy model and the Hubble constant tension_ , _Phys. Rev. D_ 104 (2021) 063521 [arXiv:2103.03815].
  * [6] S. D. Odintsov, D. Sáez-Chillón Gómez and G. S. Sharov, _Analyzing the $H_{0}$ tension in $F(R)$ gravity models_, _Nucl. Phys. B_ 966 (2021) 115377 [arXiv:2011.03957].
  * [7] A. Bernui, E. Di Valentino, W. Giarè, S. Kumar and R. C. Nunes, _Solution of $H_{0}$ tension with evidence of dark sector interaction from 2D BAO measurements_, arXiv:2301.06097.
  * [8] S. Capozziello and M. De Laurentis, _Extended theories of gravity_ , _Phys. Rep._ 509 (2011) 167 [arXiv:1108.6266].
  * [9] T. Clifton and P. G. Ferreira, A. Padilla and C. Skordis, _Modified gravity and cosmology_ , _Phys. Rep._ 513 (2012) 1 [arXiv:1106.2476].
  * [10] L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, _Conditions for the cosmological viability of $f(R)$ dark energy models_ _Phys. Rev. D_ bf 75 (2007) 083504 [gr-qc/0612180].
  * [11] A. de la Cruz-Dombriz, A. Dobado and A. L. Maroto, _Comment on “Viable Singularity-Free $f(R)$ Gravity without a Cosmological Constant”_, _Phys. Rev. Lett._ 103 (2009) 1 [arXiv:0910.1441].
  * [12] P. Bessa, M. Campista and A. Bernui, _Observational constraints on Starobinsky $f(R)$ cosmology from cosmic expansion and structure growth data_, _EPJC_ 82 (2022) 506 [arXiv:2112.00822].
  * [13] L. Perenon, M. Martinelli, R. Maartens, S. Camera and C. Clarkson, _Measuring dark energy with expansion and growth_ , _Phys. Dark Universe_ 37 (2022) 101119 [arXiv:2206.12375].
  * [14] G. Alestas, L. Kazantzidis and S. Nesseris, _Machine learning constraints on deviations from General Relativity from the large scale structure of the universe_ , _Phys. Rev. D_ 106 (2022) 103519 [arXiv:2209.12799].
  * [15] D. Staicova, _DE models with combined $H0\cdot r_{d}$ from BAO and CMB Dataset and Friends_, _Universe_ 8 (2022) 631 [arXiv:2211.08139].
  * [16] D. Staicova and D. Benisty, _Constraining the dark energy models using Baryon Acoustic Oscillations: an approach independent of $H_{0}\cdot r_{d}$_, _A $\&$A_ 668 (2022) A135 [arXiv:2107.14129].
  * [17] A. De Felice and S. Tsujikawa, _$f(R)$ theories_, _Living Rev. Relativ._ 13 (2010) 1 [arXiv:1002.4928].
  * [18] T. P. Sotiriou and V. Faraoni, _$f(R)$ theories of gravity_, _RMP_ 82 (2008) 451 [arXiv:0805.1726].
  * [19] A. Starobinsky, _A new type of isotropic cosmological models without singularity_ , _Phys. Lett. B_ 91 (1980) 99.
  * [20] A. Starobinsky, _Disappearing cosmological constant in $f(R)$ gravity_, _JETP Lett._ 86 (2007) 157 [arXiv:0706.2041].
  * [21] W. Hu and I. Sawicki, _Models of $f(R)$ cosmic acceleration that evade solar system tests_, _Phys. Rev. D_ 76 (2007) 064004 [arXiv:0705.1158].
  * [22] S. Appleby and R. Battye, _Do consistent $F(R)$ models mimic General Relativity Plus $\Lambda$?_, _Phys. Lett. B_ 654 (2007) 7 [arXiv:0705.3199].
  * [23] B. Li and J. Barrow, _Cosmology of $f(R)$ gravity in the metric variational approach_, _Phys. Rev. D_ 75 (2007) 084010 [gr-qc/0701111].
  * [24] L. Amendola and S. Tsujikawa, _Phantom crossing, equation-of-state singularities, and local gravity constraints in $f(R)$ models_, _Phys. Lett. B_ 660 (2008) 125 [arXiv:0705.0396].
  * [25] S. Tsujikawa, _Observational signatures of $f(R)$ dark energy models that satisfy cosmological and local gravity constraints_, _Phys. Rev. D_ 77 (2008) 023507 [arXiv:0709.1391].
  * [26] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani and S. Zerbini, _Class of viable modified $f(R)$ gravities describing inflation and the onset of accelerated expansion_, _Phys. Rev. D_ 77 (2008) 046009 [arXiv:0712.4017].
  * [27] E. V. Linder, _Exponential gravity_ , _Phys. Rev. D_ 80 (2009) 123528 [arXiv:0905.2962].
  * [28] E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani and S. Zerbini, _Nonsingular exponential gravity: a simple theory for early and late-time accelerated expansion_ , _Phys. Rev. D_ 83 (2011) 086006 [arXiv:1012.2280].
  * [29] Q. Xu and B. Chen, _A new exponential gravity_ , _Commun. Theor. Phys_ 61 (2014) 141 [arXiv:1203.6706].
  * [30] A. Nautiyal, S. Panda and A. Patel, _A new $f(R)$ model in the light of local gravity test and late-time cosmology_, _Int. J. Mod. Phys. D_ 27 (2018) 1750185 [arXiv:1607.05583].
  * [31] D. Gogoi and U. Dev Goswami, _A new $f(R)$ gravity model and properties of gravitational waves in it_, _Eur. Phys. J. C_ 80 (2020) 1 [arXiv:2006.04011].
  * [32] V. K. Oikonomou, _An exponential $F(R)$ dark energy model_, _Gen. Relativ. Gravit._ 45 (2013) 2467 [arXiv:1304.4089].
  * [33] V. K. Oikonomou, _Unifying inflation with early and late dark energy epochs in axion $F(R)$ gravity_, _Phys. Rev. D_ 103 (2021) 044036 [arXiv:2012.00586].
  * [34] F. Avila, C. Novaes, A. Bernui, E. de Carvalho and J. P. Nogueira-Cavalcante, _The angular scale of homogeneity in the local Universe with the SDSS blue galaxies_ , _MNRAS._ 488 (2019) 1481 [arXiv:1906.10744].
  * [35] G. A. Marques and A. Bernui, _Tomographic analyses of the CMB lensing and galaxy clustering to probe the linear structure growth_ , _JCAP._ 2020 (2020) 052 [arXiv:1908.04854].
  * [36] E. de Carvalho, A. Bernui, F. Avila, C. P. Novaes and J. P. Nogueira-Cavalcante, _BAO angular scale at z eff = 0.11 with the SDSS blue galaxies_ _A $\&$A._ bf 649 (2021) A20 [arXiv:2103.14121].
  * [37] S. Appleby, R. Battye and A. Starobinsky, _Curing singularities in cosmological evolution of $F(R)$ gravity_ _JCAP_ bf 2010 (2010) 005 [arXiv:0909.1737].
  * [38] S. Appleby and R. Battye, _Aspects of cosmological expansion in $F(R)$ gravity models_, _JCAP_ 2008 (2008) 019 [arXiv:0803.1081].
  * [39] A. Dev, D. Jain, S. Jhingan, S. Nojiri, M. Sami and I. Thongkool, _Delicate $f(R)$ gravity models with a disappearing cosmological constant and observational constraints on the model parameters_, _Phys. Rev. D_ 78 (2008) 083515 [arXiv:0807.3445].
  * [40] V. Faraoni and S. Capozziello, _Beyond Einstein Gravity: A Survey of Gravitational Theories for Cosmology and Astrophysics_ , Springer. (2011).
  * [41] E. Papantonopoulos, _Modifications of Einstein’s Theory of Gravity at Large Distances_ , Springer. (2014).
  * [42] V. Muller, H. J. Schmidt and A. Starobinsky, _The stability of the de Sitter space-time in fourth order gravity_ , _Phys. Lett. B_ 202 (1988).
  * [43] A. Nunez and S. Solganik, _The content of $f(R)$ gravity_ (2004), hep-th/0403159.
  * [44] A. Krause and S. P. Ng, _Ghost cosmology: exact solutions, transitions between standard cosmologies and ghost evolution_ , _Int. J. Mod. Phys. A_ 21 (2006) 1091 [hep-th/0409241].
  * [45] B. Himmetoglu, C. R. Contaldi and M. Peloso, _Ghost instabilities of cosmological models with vector fields nonminimally coupled to the curvature_ , _Phys. Rev. D_ 80 (2009) 123530 [arXiv:0909.3524].
  * [46] N. Deruelle, M. Sasaki, Y. Sendouda and A. Youssef, _Inflation with a Weyl term, or ghosts at work_ , _JCAP_ 2011 (2011) 040 [arXiv:1012.5202].
  * [47] A. D. Dolgov and M. Kawasaki, _Can modified gravity explain accelerated cosmic expansion?_ , _Phys. Lett. B_ 573 (2003) 1 [astro-ph/0307285].
  * [48] G. J. Olmo, _Post-Newtonian constraints on cosmologies in metric and Palatini formalism_ , _Phys. Rev. D_ 72 (2005) 083505 [gr-qc/0505135].
  * [49] V. Faraoni, _Solar system experiments do not yet veto modified gravity models_ , _Phys. Rev. D_ 74 (2006) 023529 [gr-qc/0607016].
  * [50] W. Hu, _Wandering in the background: a Cosmic Microwave Background explorer”, howpublished_ , Ph.D. thesis, University of California, USA (1995), [astro-ph/9508126].
  * [51] W. Hu and M. White, _Acoustic signatures in the Cosmic Microwave Background_ , _ApJ_ 471 (1996) 30 [astro-ph/9602019].
  * [52] W. Hu and S. Dodelson, _Cosmic Microwave Background anisotropies_ , _Annu. Rev. Astron. Astrophys._ 40 (2002) 171 [astro-ph/0110414].
  * [53] H. Kodama and M. Sasaki, _Cosmological perturbation theory_ , _PTP Supplement_ 78 (1984) 1.
  * [54] V. F. Mukhanov, H. A. Feldman and R. H. Brandenberger, _Theory of cosmological perturbations_ , _Phys. Rep._ 215 (1992) 203.
  * [55] J. M. Bardeen, _Gauge-invariant cosmological perturbations_ , _Phys. Rev. D_ 22 (1980) 1882.
  * [56] P. J. E. Peebles, _The gravitational instability of the universe_ , _ApJ_ 147 (1967) 859.
  * [57] Y. B. Zel’Dovich, _Gravitational instability: an approximate theory for large density perturbations_ , _A $\&$A_ 5 (1970) 84.
  * [58] S. Tsujikawa, _Matter density perturbations and effective gravitational constant in modified gravity models of dark energy_ , _Phys. Rev. D_ 76 (2007) 023514 [arXiv:0705.1032].
  * [59] M. A. Strauss and J. A. Willick, _The density and peculiar velocity fields of nearby galaxies_ , _Phys. Rep._ 261 (1995) 271 [arXiv:9502079].
  * [60] L. Wang and P. J. Steinhardt, _Cluster abundance constraints for cosmological models with a time-varying, spatially inhomogeneous energy component with negative pressure_ , _ApJ_ 508 (1998) 483 [arXiv:9804015].
  * [61] E. V. Linder and R. N. Cahn, _Parameterized beyond-Einstein growth_ , _Astropart. Phys._ 28 (2007) 4 [astro-ph/0701317].
  * [62] S. Nesseris, G. Pantazis and L. Perivolaropoulos, _Tension and constraints on modified gravity parametrizations of $G_{eff}(z)$ from growth rate and Planck data_, _Phys. Rev. D_ 96 (2017) 023542 [arXiv:1703.10538].
  * [63] N. Aghanim el al., _Planck $2018$ results. VI. Cosmological parameters_, _A $\&$A_ 641 (2020) [arXiv:1807.06209].
  * [64] C. Zhang, H. Zhang, S. Yuan, S. Liu, T. Zhang and Y. Sun, _Four new observational $H(z)$ data from luminous red galaxies in the Sloan Digital Sky Survey data release seven_, _Res. Astron. Astrophys._ 14 (2014) 1221 [arXiv:1207.4541].
  * [65] J. Simon, L. Verde and R. Jimenez, _Constraints on the redshift dependence of the dark energy potential_ , _Phys. Rev. D_ 71 (2005) 123001 [astro-ph/0412269].
  * [66] M. Moresco et al., _Improved constraints on the expansion rate of the universe up to $z\sim 1.1$ from the spectroscopic evolution of cosmic chronometers_, _JCAP_ 2012 (2012) 006 [arXiv:1201.3609].
  * [67] M. Moresco, L. Pozzetti, A. Cimatti, R. Jimenez, C. Maraston, L. Verde, D. Thomas, A. Citro, R. Tojeiro and D. Wilkinson, _A $6\%$ measurement of the Hubble parameter at $z~{}0.45$: direct evidence of the epoch of cosmic re-acceleration_, _JCAP_ 2016 (2016) 014 [arXiv:1601.01701].
  * [68] A. Ratsimbazafy, S. I. Loubser, S. M. Crawford, C. M. Cress, B. A. Bassett, R. C. Nichol and P. Väisänen, _Age-dating luminous red galaxies observed with the Southern African large telescope_ , _MNRAS_ 467 (2017) 3239 [arXiv:1702.00418].
  * [69] D. Stern, R. Jimenez, L. Verde, M. Kamionkowski and S. A. Stanford, _Cosmic chronometers: constraining the equation of state of dark energy. I: $H(z)$ measurements_, _JCAP_ 2010 (2010) 008 [arXiv:0907.3149].
  * [70] M. Moresco, _Raising the bar: new constraints on the Hubble parameter with cosmic chronometers at $z~{}2$_, _MNRAS Lett._ 450 (2015) L16 [arXiv:1503.01116].
  * [71] R. Jimenez and A. Loeb, _Constraining cosmological parameters based on relative galaxy ages_ , _ApJ_ 573 (2002) 37 [astro-ph/0106145].
  * [72] G. Bruzual and S. Charlot, _Stellar population synthesis at the resolution of 2003_ , _MNRAS_ 344 (2003) 1000 [astro-ph/0309134].
  * [73] A. Gómez-Valent, _Quantifying the evidence for the current speed-up of the universe with low and intermediate-redshift data. A more model-independent approach_ , _JCAP_ 2019 (2019) 026 [arXiv:1810.02278].
  * [74] Y. Yang and Y. Gong, _The evidence of cosmic acceleration and observational constraints_ , _JCAP_ 2020 (2020) 059 [arXiv:1912.07375].
  * [75] F. Avila, A. Bernui, A. Bonilla and R. C. Nunes, _Inferring $S_{8}(z)$ and $\gamma(z)$ with cosmic growth rate measurements using machine learning_, _EPJC_ 82 (2022) [arXiv:2201.07829].
  * [76] S. J. Turnbull, M. J. Hudson, H. A. Feldman, M. Hicken, R. P. Kirshner and R. Watkins, _Cosmic flows in the nearby universe from type Ia supernovae_ , _MNRAS_ 420 (2012) 447 [arXiv:1111.0631].
  * [77] I. Achitouv, C. Blake, P. Carter, J. Koda and F. Beutler, _Consistency of the growth rate in different environments with the $6$dF Galaxy Survey: measurement of the void-galaxy and galaxy-galaxy correlation functions_, _Phys. Rev. D_ 95 (2016) 083502 [arXiv:1606.03092].
  * [78] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, G. B. Poole, L. Campbell, Q. Parker, W. Saunders and F. Watson, _The $6$dF Galaxy Survey: $z\approx 0$ measurements of the growth rate and $\sigma_{8}$_, _MNRAS_ 423 (2012) 3430 [arXiv:1204.4725].
  * [79] M. Feix, A. Nusser and E. Branchini, _Growth rate of cosmological perturbations at $z\sim 0.1$ from a new observational test_, _Phys. Rev. Lett._ 115 (2015) 011301 [arXiv:1503.05945].
  * [80] S. Alam et al., _The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR $12$ galaxy sample_, _MNRAS_ 470 (2017) 2617 [arXiv:1607.03155].
  * [81] A. G. Sánchez et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological implications of the full shape of the clustering wedges in the data release 10 and 11 galaxy samples_ , _MNRAS_ 440 (2014) 2692 [arXiv:1312.4854].
  * [82] C. Blake et al., _The WiggleZ dark energy survey: joint measurements of the expansion and growth history at $z<1$_, _MNRAS_ 425 (2012) 405 [arXiv:1204.3674].
  * [83] S. Nadathur, P. M. Carter, W. J. Percival, H. A. Winther and J. Bautista, _Beyond BAO: improving cosmological constraints from BOSS data with measurement of the void-galaxy cross-correlation_ , _Phys. Rev. D_ 100 (2019) 023504 [arXiv:1904.01030].
  * [84] C.-H. Chuang et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: single-probe measurements from CMASS anisotropic galaxy clustering_ , _MNRAS_ 461 (2016) 3781 [arXiv:1312.4889].
  * [85] M. Aubert et al., _The completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: growth rate of structure measurement from Cosmic voids_ , _MNRAS_ 513 (2022) 186 [arXiv:1312.4889].
  * [86] M. J. Wilson, _Geometric and growth rate tests of General Relativity with recovered linear cosmological perturbations_ , Ph.D. thesis, Edinburgh University, USA (2017), [arXiv:1610.08362].
  * [87] G.-B. Zhao et al., _The clustering of the SDSS-IV extended Baryon Oscillation Spectroscopic Survey DR14 quasar sample: a tomographic measurement of cosmic structure growth and expansion rate based on optimal redshift weights_ , _MNRAS_ 482 (2018) 3497 [arXiv:1801.03043].
  * [88] T. Okumura et al., _The subaru FMOS Galaxy Redshift Survey (fastsound). IV. New constraint on gravity theory from redshift space distortions at $z\sim 1.4$_, _PASJ_ 68 (2016) 38 [arXiv:1511.08083].
  * [89] F. Avila, A. Bernui, E. de Carvalho and C. P. Novaes, _The growth rate of cosmic structures in the local universe with the ALFALFA survey_ , _MNRAS._ 505 (2021) 3404 [arXiv:2105.10583].
  * [90] F. Avila, A. Bernui, R. C. Nunes, E. de Carvalho and C. P. Novaes, _The homogeneity scale and the growth rate of cosmic structures_ , _MNRAS._ 509 (2022) 2994 [arXiv:2111.08541].
  * [91] C. Robert and G. Casella, _Monte Carlo Statistical Methods_ , Springer New York. (2005).