# Low-lying baryon masses using twisted mass fermions ensembles at the
physical pion mass

Constantia Alexandrou  Department of Physics, University of Cyprus, P.O. Box
20537, 1678 Nicosia, Cyprus  Computation-based Science and Technology Research
Center, The Cyprus Institute, Nicosia, Cyprus  Simone Bacchio  Computation-
based Science and Technology Research Center, The Cyprus Institute, Nicosia,
Cyprus  Georgios Christou  Department of Physics, University of Cyprus, P.O.
Box 20537, 1678 Nicosia, Cyprus  Jacob Finkenrath  Computation-based Science
and Technology Research Center, The Cyprus Institute, Nicosia, Cyprus
University of Wuppertal, Wuppertal, Germany

###### Abstract

We investigate the low-lying baryon spectrum using three $N_{f}=2+1+1$
ensembles simulated with physical values of the quark masses and lattice
spacings of $0.080,\;0.069\text{ and }0.057$ fm. The ensembles are generated
using twisted mass clover-improved fermions and the Iwasaki gauge action. The
spatial length is kept approximately the same at about 5.1 fm to 5.5 fm
fulfilling the condition $m_{\pi}L>3.6$. We investigate isospin splitting
within isospin multiples and verify that for most cases the isospin splitting
for these lattice spacing is consistent with zero. In the couple of cases, for
which there is a non-zero value, in the continuum limit, the mass splitting
goes to zero. The baryon masses are extrapolated to the continuum limit using
the three $N_{f}=2+1+1$ ensembles and are compared to other recent lattice QCD
results. For the strange and charm quark masses we find, respectively,
$m_{s}(2~{}\rm{GeV})=99.2(2.7)$ MeV and $m_{c}(3~{}{\rm GeV)}=1.015(39)$ GeV.
The values predicted for the masses of the doubly charmed $\Xi_{cc}^{\star}$,
$\Omega_{cc}$ and $\Omega_{cc}^{\star}$ baryons are 3.676(55) GeV, 3.703(51)
GeV and 3.803(50) GeV, respectively, and for the triply charmed $\Omega_{ccc}$
baryon is 4.785(71) GeV.

##  I Introduction

Baryon spectroscopy is an active field of research. Experimentally, new
baryons are being studied at LHCb and CERN, by the SELEX and FOCUS experiments
at Fermilab, the BaBaR experiment at SLAC, and the Belle-II at KEK. The
discovery of new baryons like the doubly charmed $\Xi$ resonance by SELEX and
LHCb triggered a revival of the interest in charmed baryon spectroscopy. The
Beijing Spectrometer BES-III has an extensive program in spectroscopy and so
does PANDA planned at GSI [1].

There are also many theoretical studies in baryon spectroscopy using QCD sum
rules [2], as well as relativistic [3, 4] and non-relativistic quark models
[5]. Recent progress in algorithms and access to larger computational
resources have led to tremendous progress in simulating lattice QCD at
physical values of the quark masses. Therefore, lattice QCD is in a good
position to investigate the masses of doubly and triply charmed baryons using
simulations with physical values of the quark masses providing valuable input
directly from first principles. While baryon masses of low-lying hyperons and
singly charmed baryons are well determined experimentally and thus serve as
benchmark quantities, the doubly and triply charmed sector still remain mostly
unexplored. A number of lattice QCD groups have studied the ground states of
spin-1/2 and spin-3/2 charmed baryons using a variety of lattice schemes, with
the most recent ones using dynamical simulations [6, 7, 8, 9, 10, 11, 12, 13,
14]. However, these calculations still performed a chiral extrapolation since
ensembles with multiple lattice spacings at the physical point are limited.

In this work, we use three ensembles generated by the Extended Twisted Mass
Collaboration (ETMC) with two degenerate twisted mass clover-improved light
quarks with mass tuned to reproduce the physical pion mass, plus a strange and
charm quark with masses fixed to approximately their physical ones. We will
refer to these $N_{f}=2+1+1$ ensembles as physical point ensembles. These
physical point ensembles have lattice spacings $a=0.08$ fm, 0.069 fm, and
0.057 fm and approximately equal volume. This enables us, for the first time,
to take the continuum limit directly at the physical point, thus eliminating a
systematic uncertainty arising from the chiral extrapolation that in the
baryon sector can introduce an uncontrolled systematic error. Including a
clover-term helps in the stabilization of the simulations, while it still
preserves automatic $\mathcal{O}(a)$-improvement of the twisted mass action
and reduces the $\mathcal{O}(a^{2})$ lattice artifacts related to the breaking
of the isospin symmetry. This study extends our previous computations on the
low-lying baryon spectrum using $N_{f}=2$ with one physical point ensemble
[15].

Ensemble | $V$ | $\beta$ | $\mu_{\ell}$ | $\mu_{\sigma}$ | $\mu_{\delta}$ | $\kappa$ | $c_{SW}$ | $L\cdot m_{\pi}$ | $m_{\pi}[$\mathrm{M}\mathrm{e}\mathrm{V}$]$
---|---|---|---|---|---|---|---|---|---
cB211.072.64 | $128\times 64^{3}$ | 1.778 | 0.00072 | 0.1246826 | 0.1315052 | 0.1394265 | 1.69 | 3.62 | 140.1 (0.2)
cC211.060.80 | $160\times 80^{3}$ | 1.836 | 0.00060 | 0.106586 | 0.107146 | 0.13875285 | 1.6452 | 3.78 | 136.7 (0.2)
cD211.054.96 | $192\times 96^{3}$ | 1.900 | 0.00054 | 0.087911 | 0.086224 | 0.137972174 | 1.6112 | 3.90 | 140.8 (0.2)
Table 1:  Simulation parameters used by ETMC for generating each ensemble [16]. The first column gives the ensemble name, the second the lattice volume, the third the $\beta$-value, and the fourth, fifth, and sixth the quark mass parameters. In the seventh and eighth columns, we give the critical value of $\kappa$ for which $m_{l}\sim 0$ and the value of $c_{SW}$. In the last but one column we give $Lm_{\pi}$ and in the last column the pion mass. Ensemble | $a_{G}$ | $n_{G}$ | $a_{\text{APE}}$ | $n_{\text{APE}}$
---|---|---|---|---
cB211.072.64 | 1.0 | 95 | 0.5 | 50
cC211.060.80 | 1.0 | 140 | 0.5 | 60
cD211.054.96 | 1.0 | 200 | 0.5 | 60
Table 2: Smearing parameters used for all ensembles.

We follow our previous studies and use Osterwalder-Seiler (OS) valence strange
and charm quarks and choose the physical masses of the $\Omega^{-}$ and the
$\Lambda_{c}^{+}$ baryons to tune the OS valence strange and charm quark
masses, respectively. We also use the nucleon mass to fix the lattice spacing
in order to convert our lattice values to physical units. Isospin breaking in
the twisted mass formulation is a lattice artifact of order ${a^{2}}$. It has
been shown that adding the clover term reduces isospin splitting in the
$\Delta$-multiplet [17] as compared to the $N_{f}=2+1+1$ twisted mass
simulations at a similar lattice spacing. Here we study the effects of isospin
breaking effects to higher accuracy for the $\Delta$ multiplet and in the
strange and charm sectors. We also extract the renormalized strange and charm
quark masses and find agreement with our previous determination [18]. We
compare our final results on the masses of the forty baryons studied in this
work with those of other recent lattice calculations, using a variety of
discretization schemes, as well as with experiments. We find agreement with
experimental results. Our predictions for the doubly charmed
$\Xi_{cc}^{\star}$, $\Omega_{cc}$ and $\Omega_{cc}^{\star}$ baryons are
$M_{\Xi_{cc}^{\star}}=3.676(55)\text{GeV}$,
$M_{\Omega_{cc}}=3.703(51)\text{GeV}$ and
$M_{\Omega_{cc}^{\star}}=3.803(50)\text{GeV}$, respectively, and for the
triply charmed $\Omega_{ccc}$ baryon is
$M_{\Omega_{ccc}}=4.785(71)\text{GeV}$.

The paper is organized as follows: the lattice action employed in this work,
as well as the details of the calculations, including the interpolating fields
and the effective mass analysis procedure, are given in Section II. The
determination of the lattice spacing and the tuning of the strange and charm
quark masses are given in Sections III and IV respectively, and we present our
lattice results in Section V. In Section , we study the isospin splitting of
the spin-1/2 and spin-3/2 baryons. In Section VII we present our results for
the baryon spectrum at the continuum limit and in Section VIII a comparison
between our results and other lattice QCD results is presented. Lastly, we
present our conclusions in Section IX.

##  II Lattice setup

###  II.1 The lattice action

We analyze three gauge ensembles produced by ETMC [19, 16, 20] within the
twisted mass fermion formulation [21], whose parameters are given in Table 1.
The $N_{f}$ = 2 + 1 + 1 ensembles are generated using the Iwasaki gauge
action, the $N_{f}$= 2 mass-degenerate twisted mass fermion action for the
light doublet with a clover term [22] and the $N_{f}$ = 1 + 1 non-degenerate
twisted mass fermion action for the strange and charm quarks. Including a
clover term, reduces lattice artifacts and decreases the mass gap between the
charged and neutral pion. We use 1-loop tadpole boosted perturbation theory to
fix the value of the clover parameter $c_{\rm SW}$. The light quark mass is
tuned to reproduce the isosymmetric pion mass $m_{\pi}=0.135$ MeV, the strange
quark mass is tuned to reproduce the physical value of the ratio between
$D_{s}$ meson mass and decay constant, and the charm quark is fixed by setting
the ratio of charm to strange quark mass $\mu_{c}/\mu_{s}=11.8$ [16].

Ensemble | $\mu^{\prime}_{\ell}$
---|---
cB211.072.64 | 0.0006675
cD211.054.96 | 0.0004964
Table 3: Light quark twisted masses $\mu_{\ell}$ used for correcting the
effective mass [23].

For completeness, we provide the action used in the simulations. The Iwasaki
gauge action is given by

$\begin{split}S_{g}&=\frac{\beta}{3}\sum_{x}\left(b_{0}\sum_{\mu<\nu}\left[1-\real\tr(U_{\mu\nu}^{1\times
1}(x))\right]\right.\\\
&\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left.+b_{1}\sum_{\mu<\nu}\left[1-\real\tr(U_{\mu\nu}^{1\times
2}(x))\right]\right),\end{split}$ (1)

where $\beta=6/g_{0}^{2}$ is the inverse bare coupling constant and we set the
parameters $b_{0}=1-8b_{1}$ and $b_{1}=-0.331$ such that the Iwasaki improved
gauge action [24] is reproduced.

Figure 1: The effective mass for the nucleon using the cD211.054.96 ensemble. On the top panel, we present one- (grey band), two- (red band) and three-state (green) fits for which the extracted nucleon mass is consistent. On the bottom left panel, we show the value extracted from the plateau fit (green points) as a function of the value of $t_{\rm low}$ used in the fit, while, on the bottom right panel, we show the energy difference $\Delta E$ between the first excited state and the ground state extracted from the two- (red) and three-state (black) fits. On the bottom panels, the filled symbols give the selected values for the ground state (bottom left) and $\Delta E$ (bottom right), and the error band is centered around the result of the most probable model. All values are given in lattice units. Ensemble | $t_{\text{max}}/a$ | one-state | two-state | three-state
---|---|---|---|---
$t_{\text{low}}/a$ | $aM_{B}$ | $t_{\text{low}}/a$ | $aM_{B}$ | $t_{\text{low}}/a$ | $aM_{B}$
cB211.072.64 | 34 | 18 | 0.3819(18) | 8 | 0.3796(24) | 3 | 0.3798(34)
cC211.060.80 | 38 | 18 | 0.32756(78) | 8 | 0.3258(12) | 3 | 0.3257(17)
cD211.054.96 | 44 | 23 | 0.2724(11) | 9 | 0.2723(11) | 2 | 0.27330(84)
Table 4: Nucleon effective mass results from the most probable one-, two- and
three-state fits. For each fit, we give in lattice units the maximum and
minimum times $t_{\rm maax}$ and $t_{\rm low}$ used in the fit, as well as the
baryon mass extracted from each fit.

For the fermionic sector, we distinguish between light- $(u,d)$ and heavy-
flavored $(c,s)$ quarks. We construct the light and heavy quark doublets
$\chi_{\ell}=(u,d)^{T}$ and $\chi_{h}=(s,c)^{T}$ in the twisted basis,
respectively. The light quark action reads

$\begin{split}S_{tm}^{\ell}&=\sum_{x}\bar{\chi}_{\ell}(x)\left[D_{W}[U]+\frac{i}{4}c_{SW}\sigma_{\mu\nu}\mathcal{F}^{\mu\nu}+m_{\ell}\right.\\\
&\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\left.\frac{}{}i\mu_{\ell}\gamma_{5}\tau^{3}\right]\chi_{\ell}(x),\end{split}$
(2)

and the heavy quark part of the action reads

$\begin{split}S_{tm}^{h}&=\sum_{x}\bar{\chi}_{h}(x)\left[D_{W}[U]+\frac{i}{4}c_{SW}\sigma_{\mu\nu}\mathcal{F}^{\mu\nu}+m_{h}\right.\\\
&\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-\mu_{\delta}\tau^{1}+\left.\frac{}{}i\mu_{\sigma}\gamma_{5}\tau^{3}\right]\chi_{h}(x).\end{split}$
(3)

The third Pauli matrix $\tau^{3}$ acts in flavor space, $D_{W}$ is the
massless Wilson-Dirac operator and $\mu_{l}$ is the light quark twisted mass.
Automatic $\mathcal{O}(a^{2})$ improvement is achieved by tuning $m_{l}$ to
zero using the partially conserved axial current (PCAC) relation. In the heavy
quark action, the additional term $\mu_{\delta}\tau^{1}$ is acting on flavor
space, thus, the mass degeneracy, in contrast to the light quark sector, is
lifted. The massless Wilson-Dirac operator is given by

$D_{W}[U]=\frac{1}{2}\gamma_{\mu}\left(\nabla_{\mu}+\nabla_{\mu}^{\star}\right)-\frac{ar}{2}\nabla_{\mu}\nabla_{\mu}^{\star},$
(4)

where the forward derivatives read

$\nabla_{\mu}\psi(x)=\frac{1}{a}\left[U_{\mu}^{\dagger}(x)\psi(x+a\hat{\mu}-\psi(x)\right]$
(5)

and

$\nabla_{\mu}^{\star}\psi(x)=-\frac{1}{a}\left[U_{\mu}(x-a\hat{\mu})\psi(x-a\hat{\mu}-\psi(x)\right].$
(6)

At or close to maximal twist, we transform from twisted to physical basis
using $\psi(x)=\frac{1}{\sqrt{2}}\left[1+i\tau^{3}\gamma_{5}\right]\chi(x)$
and
$\bar{\psi}(x)=\bar{\chi}(x)\frac{1}{\sqrt{2}}\left[1+i\tau^{3}\gamma_{5}\right]$.
Throughout this paper, unless otherwise stated, we use quark fields in the
physical basis, in particular when we define baryon interpolators. For the
valence strange and charm quarks, we utilize Osterwalder Seiler (OS) fermions
[25]. We discuss the tuning of the OS strange and charm quark masses below.

###  II.2 Two-point functions and effective masses

In this work, we consider positive parity baryon correlators of two-point
functions of the form

$C_{B}^{\pm}(t,\vec{0})=\sum_{\vec{x}_{f}}\expectationvalue{\frac{1}{4}\tr(1\pm\gamma_{0})\mathcal{J}_{B}(\vec{x}_{f},t_{f})\bar{\mathcal{J}}_{B}(\vec{x}_{i},t_{i})}\,,$
(7)

where $\bar{\mathcal{J}}_{B}$ is the interpolating field of the baryon acting
at the source $(\vec{x}_{i},t_{i})$ and ${\mathcal{J}}_{B}$ at the sink
$(\vec{x}_{f},t_{f})$. We consider three-quark baryon interpolating fields of
the form $\varepsilon_{abc}\left(q_{1,a}^{T}\Gamma q_{2,b}\right)q_{3,c}$. The
$\Gamma$ structures we use are $C\gamma_{5}$ and $C\gamma_{i},i=1,\ldots,3$
for spin-$1/2$ and $3/2$, respectively.

We exploit the symmetries of the action and the anti-periodic boundary
conditions in the temporal direction to write

$C^{+}_{B}(t,\vec{0})=-C^{-}_{B}(T-t,\vec{0}).$ (8)

Thus, we average the correlators in the forward and backward direction to
improve the signal by taking

$C_{B}(t,\vec{0})=C^{+}_{B}(t,\vec{0})-C^{-}_{B}(T-t,\vec{0}).$ (9)

For the spin-3/2 baryons, we project to spin $3/2$ by performing the following
transformations at zero momentum [26]

$C_{\frac{3}{2}}(t)=\frac{1}{3}\tr(C(t))+\frac{1}{6}\sum_{i\neq
j}^{3}\gamma_{i}\gamma_{j}C_{ij}(t),$ (10)

where $t=t_{f}-t_{i}$ and from now on we drop the momentum dependence in the
argument of the two-point function correlator since we only consider
$\vec{p}=\vec{0}$.

Since, in this work, we are interested in the low-lying spectrum and gauge
noise grows with the time $t$, we need to use techniques that suppress excited
state contamination. We apply Gaussian smearing to the quark fields at the
source and sink [27, 28] using

$q_{\text{smear}}(n)=\sum_{m\in\Lambda}F(n,m,U)q(m),$ (11)

where the gauge invariant smearing function reads

$F(n,m,U)=(1+\alpha_{G}H)^{n_{G}}(n,m,U)$ (12)

and $H$ is the hopping function

$H(n,m,U)=\sum_{j=1}^{3}\left(U_{j}(n)\delta_{n+a\hat{j},m}+U^{\dagger}_{j}(n-a\hat{j})\delta_{n-a\hat{j},m}\right).$
(13)

The parameters $a_{G}$ and $n_{G}$ for the three ensembles are given in Table
2. We also apply APE smearing to the links that enter the hopping function
$H(n,m,U)$ with APE smearing parameter $a_{\text{APE}}$ and number of
iterations $n_{\text{APE}}$ for all ensembles given in Table 2.

To extract the low-lying masses, we construct the effective mass using
correlators at zero momentum

$am_{\text{eff}}^{B}(t)=\log(\frac{C(t)}{C(t+1)})\xrightarrow[t\to\infty]{}aM_{B}.$
(14)

We adopt a multi-state fit approach to identify the ground state keeping up to
the second excited state. Namely, we fit the effective mass to a constant, and
to the following expressions taking into account the first and second excited
states, respectively

$am_{\text{eff}}^{B}(t)\approx aM_{B}+\log(\frac{1+\alpha_{2}e^{-\Delta
E_{1}t}}{1+\alpha_{2}e^{-\Delta E_{1}(t+1)}})$ (15)

and

$am_{\text{eff}}^{B}(t)\approx aM_{B}+\log(\frac{1+\alpha_{2}e^{-\Delta
E_{1}t}+\alpha_{3}e^{-\Delta E_{2}t}}{1+\alpha_{2}e^{-\Delta
E_{1}(t+1)}+\alpha_{3}e^{-\Delta E_{2}(t+1)}}).$ (16)

We fix the maximum time used in fit to given value $t_{\rm max}$ and
subsequently, we iteratively change the time from where we start the fit
$t_{\rm low}$ for all three types of fits. We refer to the different fit
ranges as models.

##  III Lattice spacing determination

|  $a$ [fm] |  $a$ [fm]
---|---|---
Ensembles | AIC Selection | Model Averaging
cB211.072.64 | 0.07982(51) | 0.07998(52)
cC211.060.80 | 0.06852(24) | 0.06853(27)
cD211.054.96 | 0.05747(18) | 0.05719(43)
Table 5: The values of the lattice spacing for the three ensembles using the maximum probability fit (second column) and averaging over all models (third column). Ensemble | $\mu_{s}$ | $\mu_{c}$ | $Z_{p}$
---|---|---|---
cB211.072.64 | 0.017 | 0.20 | 0.4788(34)(42)
0.019 | 0.22
cC211.060.80 | 0.015 | 0.18 | 0.4871(26)(42)
0.017 | 0.24
cD211.054.96 | 0.0125 | 0.146 | 0.4967(7)(36)
0.0152 | 0.18
Table 6: The values of the two bare quark masses used for the strange (second
column) and charm (third column) quark flavors for the interpolation to the
tuned value. In the fourth column, we give the renormalization constants [18]
$Z_{P}$ in the $\overline{\text{MS}}$ scheme at
$2\,$\mathrm{G}\mathrm{e}\mathrm{V}$$ for each ensemble, with the value in the
first and second parenthesis being the statistical and systematic error,
respectively.

For the tuning of the lattice spacing, we opt to use the physical nucleon mass
as input. From the value of $m_{\pi}$ given in Table 1, it can be seen that
two ensembles have slightly higher pion mass as compared to the one target by
our simulations, which is $135$ MeV. We correct for this small deviation by
generating at lower statistics two-point functions for lighter twisted masses
with values shown in Table 3 and applying a correction to the effective mass
as follows

$\begin{split}a\delta
m_{\text{eff}}^{B}(t)&=\log(\frac{C(t)}{C(t+1)}\frac{C^{\prime}(t+1)}{C^{\prime}(t)})\xrightarrow[t\to\infty]{}a\Delta
M_{B},\\\ \end{split}$ (17)

where $C^{\prime}(t)$ denotes the correlator evaluated at the lower value of
the twisted mass parameter.

The corrected nucleon effective mass reads

$\begin{split}am_{\text{eff}}^{\prime B}(t)&=am_{\text{eff}}^{B}(t)+a\delta
m_{\text{eff}}^{B}(t)\\\
&=\log(\frac{C^{\prime}(t)}{C^{\prime}(t+1)})+\log(\frac{C(t)}{C(t+1)}\frac{C^{\prime}(t+1)}{C^{\prime}(t)}).\end{split}$
(18)

This correction changes the nucleon mass on a percent level and ensures that
all baryon masses computed with these three ensembles correspond to the same
physical pion mass.

Figure 2: We show the effective mass analysis for the $\Lambda_{c}^{+}$ (top)
and $\Omega^{-}$ (bottom) baryon for the cD211.054.96 and cC211.060.80
ensembles respectively. On the top left we depict the effective mass data for
the two values of $\mu_{c}$ as well as the one-, two-, and three-state fits.
On the top left we show the linear interpolation. On the middle and bottom
panels, we show the results for the mass for the two $\mu_{c}$ values and the
difference between the ground and first excited state extracted from each fit.
The green and red bands represent the most probable model result. The rest of
the notation is the same as that of Fig. 1. Tuning of the strange quark mass
---
Ensemble | $a\mu_{s}$ | $t_{\text{max}}/a$ | one-state | two-state | three-state | $A_{\Omega^{-}}$\mathrm{G}\mathrm{e}\mathrm{V}$$ | $B_{\Omega^{-}}$
$t_{\text{low}}/a$ | $aM_{B}$ | $t_{\text{low}}/a$ | $aM_{B}$ | $t_{\text{low}}/a$ | $aM_{B}$
cB211.072.64 | 0.017 | 35 | 14 | 0.6768(12) | 5 | 0.6741(16) | 3 | 0.6729(22) | 1.6948(38) | 3.957(48)
0.019 | 35 | 14 | 0.6932(11) | 5 | 0.6906(15) | 3 | 0.6899(19)
cC211.060.80 | 0.015 | 43 | 20 | 0.5729(23) | 5 | 0.5750(18) | 2 | 0.5726(41) | 1.6803(51) | 3.866(50)
0.017 | 43 | 20 | 0.5893(21) | 5 | 0.5909(17) | 2 | 0.5888(37)
cD211.054.96 | 0.0125 | 52 | 23 | 0.4784(14) | 7 | 0.4774(14) | 2 | 0.4774(17) | 1.6772(46) | 3.894(48)
0.0152 | 52 | 23 | 0.4997(12) | 7 | 0.4985(13) | 2 | 0.4984(16)
Charm quark mass evaluation
Ensemble | $a\mu_{c}$ | $t_{\text{max}}/a$ | one-state | two-state | three-state | $A_{\Lambda_{c}^{+}}$\mathrm{G}\mathrm{e}\mathrm{V}$$ | $B_{\Lambda_{c}^{+}}$
$t_{\text{low}}/a$ | $aM_{B}$ | $t_{\text{low}}/a$ | $aM_{B}$ | $t_{\text{low}}/a$ | $aM_{B}$
cB211.072.64 | 0.20 | 30 | 12 | 0.8663(21) | 5 | 0.8631(35) | 1 | 0.8598(98) | 2.389(11) | 0.9463(99)
0.22 | 30 | 12 | 0.9020(22) | 5 | 0.8987(38) | 1 | 0.895(12)
cC211.060.80 | 0.18 | 31 | 12 | 0.7584(20) | 3 | 0.7579(22) | 3 | 0.757(14) | 2.4104(74) | 0.9456(59)
0.24 | 31 | 12 | 0.8633(24) | 3 | 0.8629(27) | 3 | 0.862(15)
cD211.054.96 | 0.146 | 40 | 13 | 0.6289(13) | 7 | 0.6232(38) | 2 | 0.6247(46) | 2.445(15) | 1.034(12)
0.18 | 40 | 13 | 0.6916(14) | 7 | 0.6870(42) | 2 | 0.6890(54)
Table 7: We give the $\Omega^{-}$ (top) and $\Lambda_{c}^{+}$(bottom) baryon mass results from the most probable one-, two- and three-state fits. For each fit we give $t_{\rm min}$ and $t_{\rm max}$ and the interpolation parameters from Eq. (IV) as determined from the most probable model. Ensemble | Strange quark | Charm quark
---|---|---
AIC Selection | Model Averaging | AIC Selection | Model Averaging
cB211.072.64 | 0.0894(10) | 0.0894(11) | 1.092(10) | 1.091(19)
cC211.060.80 | 0.0930(13) | 0.0931(14) | 1.0689(72) | 1.0708(93)
cD211.054.96 | 0.0938(12) | 0.0939(13) | 1.047(13) | 1.050(17)
Table 8: Tuned renormalized quark masses for each ensemble in
$\mathrm{G}\mathrm{e}\mathrm{V}$ in the $\overline{\text{MS}}$ scheme. The
strange quark masses are given at $2\,$\mathrm{G}\mathrm{e}\mathrm{V}$$ and
the charm at $2\,$\mathrm{G}\mathrm{e}\mathrm{V}$$.

For the determination of the lattice spacing we implement the following
analysis method: We use the Akaike Information Criterion (AIC) [29, 30] and
calculate the model probability for each fit function and for each fit range

$\displaystyle P(M_{i}|\\{y\\})$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\exp(-\frac{1}{2}\chi_{i}^{2}+N_{\text{d.o.f},i}),$
$\displaystyle Z$ $\displaystyle=$ $\displaystyle\sum_{i}P(M_{i}|\\{y\\}),$
(19)

where $\chi_{i}^{2}$ is the chi-squared and $N_{\text{d.o.f},i}$ is the number
of degrees of freedom used for the specific model. Subsequently, we select the
model with the highest probability and fix the nucleon mass to its physical
value. Another approach is to average over all models and fix the nucleon mass
using the model average (MA) value

$\displaystyle\expectationvalue{\alpha}$ $\displaystyle=$
$\displaystyle\sum_{i}\alpha_{i}P(M_{i}|\\{y\\}),$
$\displaystyle\sigma_{\alpha}^{2}$ $\displaystyle=$
$\displaystyle\sum_{i}(\alpha_{i}^{2}+\sigma_{\alpha,i}^{2})P(M_{i}|\\{y\\})-\expectationvalue{\alpha}^{2}.$
(20)

We show the nucleon effective mass and results of various fits in Fig. 1. In
Table 4, we give the values of the nucleon mass extracted from the plateau,
two- and three-state fits for the case of the fits for the largest
probability.

The values of the lattice spacings that we extract using the nucleon mass as
input are given in Table 5. Since both approaches employed are consistent with
each other we select the MA result as our final values of the lattice spacing
since they take into account systematic errors in the variation of the fit
ranges.

##  IV Determination of the strange and charm quark masses

For tuning the valence strange and charm quark masses we use the $\Omega^{-}$
and $\Lambda_{c}^{+}$ baryon masses respectively. In order to tune the quark
masses we generate two sets of correlators for two different values of the
non-renormalized strange and charm quark masses. The values used are listed in
Tab. 6. Since the data are correlated two sets suffice per ensemble to perform
the linear interpolation. In Fig. 2, we illustrate the analysis of the
effective mass for $\Lambda_{c}^{+}$ and $\Omega^{-}$. The procedure is
analogous to the extraction of the nucleon mass. The most probable one-, two-
and three-state fit results for $\Omega-$ and $\Lambda_{c}^{+}$ baryons are
given in Tab 7.

We implement three different methods to perform the quark mass tuning. In
method-I, we first renormalize the quark masses and interpolate to the
physical baryon masses for each ensemble using

$\displaystyle M_{\Omega}$
$\displaystyle=A_{\Omega}+B_{\Omega}(m_{s}-\tilde{m}_{s}),$ $\displaystyle
M_{\Lambda_{c}}$
$\displaystyle=A_{\Lambda_{c}}+B_{\Lambda_{c}}(m_{c}-\tilde{m}_{c}),$ (21)

where we take $\tilde{m}_{s}=0.095\,$\mathrm{G}\mathrm{e}\mathrm{V}$$ and
$\tilde{m}_{c}=1.2\,$\mathrm{G}\mathrm{e}\mathrm{V}$$. This is done for all
combinations of fits ranges, i.e. for all models. Subsequently, using Eq.
(III) we select the model with the highest probability as our final result.
Alternatively, we can average over all models with the corresponding
probabilities and calculate a model average of the tuned quark masses. The
strange and charm quark masses calculated using method-I are given in Tab. 8.

Figure 3: Continuum extrapolation of the interpolation parameters as described
in method-II for $\Omega^{-}$.

Figure 4: Continuum extrapolation of the interpolation parameters as described
in method-II for $\Lambda_{c}^{+}$.

Figure 5: The continuum limit extrapolation for the mass of the $\Omega^{-}$ (top) and the $\Lambda_{c}^{+}$ (bottom) used for method-III after tuning the mass of the strange and charm quarks. The cross represents the physical value of their mass. | $m_{s}(2~{}{\rm GeV})$ | $m_{c}(3~{}(\rm{GeV})$
---|---|---
Method-I (AIC) | 0.0991(26) | 1.006(27)
Method-I (MA) | 0.0992(27) | 1.015(39)
Method-II | 0.0994(26) | 1.012(24)
Method-III | 0.0994(26) | 1.011(24)
Table 9: We give the values of the renormalized strange and charm quark masses
at the continuum limit in $\mathrm{G}\mathrm{e}\mathrm{V}$, using all three
methods. Both values are given in the $\overline{\rm MS}$ scheme at a scale of
2 GeV and 3 GeV, respectively for the strange and charm quarks.

In method-II, we perform the interpolation for the results of the most
probable model using Eq. (IV) and then extrapolate the interpolation
parameters to the continuum limit using

$\displaystyle
A_{\Omega^{-},\Lambda_{c}^{+}}(a)=c_{\Omega^{-},\Lambda_{c}^{+}}+d_{\Omega^{-},\Lambda_{c}^{+}}a^{2},$
$\displaystyle
B_{\Omega^{-},\Lambda_{c}^{+}}(a)=e_{\Omega^{-},\Lambda_{c}^{+}}+f_{\Omega^{-},\Lambda_{c}^{+}}a^{2}.$
(22)

The continuum extrapolation of the parameters $A_{\Omega^{-},\Lambda_{c}^{+}}$
and $B_{\Omega^{-},\Lambda_{c}^{+}}$ are shown in Fig. 4 and 4. Subsequently,
we calculate the strange quark mass quark mass using

$m_{s,c}=\tilde{m}+\frac{M_{\Omega^{-},\Lambda_{c}^{+}}^{\text{phys}}-A_{s,c}(0)}{B_{s,c}(0)},$
(23)

where $M^{\text{phys}}_{\Omega}$ is the physical value of $\Omega^{-}$. We
repeat the procedure for the tuning of the charm quark using the mass of
$\Lambda_{c}$.

Figure 6: $\Delta$ baryon effective mass analysis. We use the same notation as
Fig. 1.

Figure 7: Results on the effective mass of non-charmed baryon containing a
strange quark. From top to bottom for the $\Sigma$, and $\Sigma^{\star}$
baryons computed using the cC211.060.80 ensemble. The notation is the same as
that in Fig. 2.

Figure 8: Results on the effective mass of charmed baryon containing a strange
quark. From top to bottom, we show results for the $\Sigma_{c}$, and
$\Sigma^{\star}_{c}$ baryons computed using the cC211.060.80 ensemble. The
notation is the same as that in Fig. 2.

In method-III, we adopt an iterative strategy. Namely, we consider a given
value of the renormalized strange and charm quark mass for all ensembles and
use Eq. (IV) to find the $\Omega^{-}$ and $\Lambda_{c}^{+}$ masses that
correspond to this given quark masses. We perform the continuum extrapolation
of the $\Omega^{-}$ and $\Lambda_{c}^{+}$ masses linearly in $a^{2}$ and
iteratively change the values of the renormalized strange and charm quark
masses until the continuum limit of the $\Omega^{-}$ and $\Lambda_{c}^{+}$
masses match their physical values. We present the results of method-III in
Fig. 5.

We present the results on the renormalized $m_{s}$ and $m_{c}$ in the
continuum limit extracted from all methods in Tab 9. As can be seen, all
procedures yield the same values. We take as our final values the ones from
method-I with the model average, namely

$\displaystyle m_{s}(2~{}{\rm GeV})$ $\displaystyle=99.2(2.7)~{}{\rm
MeV}\quad{\rm and}$ (24) $\displaystyle m_{c}(3~{}{\rm GeV})$
$\displaystyle=1.015(39)~{}{\rm GeV},$

in the $\overline{\text{MS}}$ scheme at the continuum limit.

##  V Extraction of baryon masses

Having determined the lattice spacing and tuned the valence strange and charm
quark masses, we discuss in this section the analysis to extract the low-lying
baryon spectrum. For the $\Delta$, where we have only $u$ and $d$ valence
quarks, only method-I is relevant since the sea and valence quarks are the
same and there is no need to perform any interpolation. However, as in the
case of the nucleon mass, we do correct the $\Delta$ effective mass to take
into account the slightly larger pion mass for the cB211.072.64 and
cC211.054.96 ensembles. This correction is noticeable enough once we take the
continuum limit and yield a result consistent with experiments. The procedure
followed is the same as that used for correcting the nucleon effective mass.
In Fig. 6, we show the effective mass for $\Delta$ and the corresponding one-,
two- and three-states fits. For baryons that contain strange and charm quarks,
we employ all three methods described in the previous section. We find the
baryon masses that correspond to the tuned masses of the strange and charm
quarks using Eq. (IV) for the ones that contain only strange or only charm
quarks. For baryons that contain both quark flavors, we use

$M_{B}=A+B(m_{s}-\tilde{m}_{s})+C(m_{c}-\tilde{m}_{c}),$ (25)

with values of $\tilde{m}_{s}$ and $\tilde{m}_{c}$ the same as those used in
Eq. (IV) employing all three methods. Our analysis of the effective masses for
selective cases is presented in Figs. 6, 7, and 8.

##  VI Isospin splitting

Figure 9: Isospin symmetry splitting for strange spin-1/2 baryons (left) and
spin-3/2 baryons (right).

Figure 10: Isospin symmetry splitting for charm spin-1/2 baryons (left) and
spin-3/2 baryons (right).

Figure 11: Continuum extrapolation of the doubly charmed baryons
$\Xi_{cc}^{\star}$, $\Omega_{cc}$ and $\Omega_{cc}^{\star}$ using the results
of method-I with model average.

Figure 12: Continuum extrapolation of the triply charmed baryon $\Omega_{ccc}$
using the results of method-I with model average.

The twisted mass fermion formulation breaks isospin symmetry at finite lattice
spacing. This is an ${\cal O}(a^{2})$ effect and it should vanish in the
continuum limit. We examine the splitting for the strange and charmed baryons
as a function of the lattice spacing and in the continuum limit. In Figs. 9
and 10, we show the results on the mass spitting for the strange and charmed
baryons, respectively, as a function of $a^{2}$. As we can see, for all
spin-3/2 non-charmed baryons the mass splittings are consistent with zero at
all three lattice spacings. For the Delta-baryon, we average the masses of the
$\Delta^{++}$ and $\Delta^{-}$ as well the masses of the $\Delta^{-}$ and
$\Delta^{0}$, since these are degenerate by symmetry. For the spin-1/2
baryons, the only mass splitting that is not consistent with zero is
$\Xi^{0}-\Xi^{-}$, which becomes zero only at the finest lattice spacing. If
we extrapolate linearly to the continuum limit we see that indeed the mass
splitting also in this case becomes consistent with zero. The mass splitting
of the $\Sigma^{+}-\Sigma^{-}$ and $\Sigma^{0}-\Sigma^{+}$ is consistent with
zero for both cB211.072.64 and cC211.060.80 ensembles. The small deviation
from zero for the cD211.054.96 ensemble of less than a standard deviation can,
thus, only be a statistical fluctuation since this ensemble has even smaller
lattice spacing as compared to the other two. Similarly, the mass splitting
between the $\Sigma^{0}$ and the charged $\Sigma^{-}$ is consistent with zero
for all ensembles within one standard deviation. For the charmed baryons, we
also observe that all mass splittings are consistent with zero for all lattice
spacings or it vanishes for the finest lattice spacing i.e. the cD211.054.96
ensemble. For the mass splitting of the $\Xi^{+}_{c}-\Xi^{0}_{c}$, that case
is similar to the one of $\Xi^{0}-\Xi^{-}$ i.e. $M_{\Xi^{+}_{c}-\Xi^{0}_{c}}$
is consistent with zero after taking the continuum limit. We include the
linear extrapolation to the continuum limit in Fig. 10 to highlight the fact
that also this splitting becomes zero at the continuum limit.

Having shown that in the continuum limit isospin splitting vanishes we present
in what follows the results on the baryons averaging among the isospin
mutliplets. For the $\Xi$ and $\Xi_{c}$ multiples, we average only in the
continuum limit.

##  VII Continuum limit

We obtain our final values on the baryon masses by extrapolating linearly in
$a^{2}$ to the continuum limit. We use all three methods except the case of
the $\Delta$ where only method-I is relevant. In Fig. 12 and 12 we present the
continuum extrapolation for the doubly and triply charmed baryons which we
predict their masses.

A comparison among the continuum extracted results using the three methods is
shown in Fig. 13 for the non-charmed baryons, in Fig. 14 for the spin-1/2
charmed baryons, and in Fig. 15 for the spin-3/2 charmed baryons. We observe
very good agreement among values for all three methods. We thus opt to use as
our final values the results obtained from method I with the model averaging,
since within this method, one takes into account all ranges for the fits
(models) eliminating a source of bias that may arise in choosing one
particular fit.

Figure 13: We show the masses of the $\Delta$ extracted with method-I and of
the strange baryons calculated using all three methods. The red circles and
blue squares denote results obtained with method-I i.e. from the fit having
the largest probability and with averaging over all models, respectively.
Black stars denote results obtained using method-II and green diamonds denote
results using method-III. The empty symbols correspond to the baryons we used
as input. Figure 14: We show the masses for the spin-$1/2$ charmed baryons
calculated using the three methods methods. The notation is the same as that
of Fig. 13. Figure 15: We show the masses for the spin-$3/2$ charmed baryons
calculated using the three methods methods. The notation is the same as that
of Fig. 13. Figure 16:  We compare our results for the masses of non-charmed
baryons (top panel) (blue squares) with the results from RQCD [31] (purple
circles), and from BMW [32](yellow stars). The horizontal bands denote
experimental results from PDG [33]. For the resonances $\Delta$,
$\Sigma^{\star}$ and $\Xi^{\star}$ the width of the band is the resonance
width.

Figure 17: We compare our results for the spin-1/2 (top panel) and spin-3/2
(bottom panel) charmed baryons (blue squares) with the results from Ref. [8]
(green diamonds) and Ref. [32] (black crosses). We also compare with the same
experimental results as in Fig. 16. For baryons that are resonances
($\Sigma_{c}^{\star}$, $\Xi_{c}^{\star}$, and $\Omega_{c}^{\star}$) we use the
resonance width as the error.

The values for the baryons masses determined using method I and the model
averaging are given in Table 10.

##  VIII Comparison with other lattice QCD results

Several collaborations have computed baryon masses, using ensembles that
include larger than physical pion masses. This work is the first one to
evaluate the low-lying baryon mass spectrum at the continuum limit using only
physical point ensembles. We compare our results to other recent lattice QCD
results. The BMW collaboration [32] performed the first groundbreaking
computation of baryon masses using unquenched gauge ensembles. They used
several $N_{f}=2+1$ ensembles of clover-improved Wilson fermions simulated
with pion masses down to $190\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ and with
lattice spacings ranging from $0.125\,$\mathrm{f}\mathrm{m}$$ down to
$0.065\,$\mathrm{f}\mathrm{m}$$. They performed simultaneous extrapolations to
the physical pion mass and to the continuum limit and since they had several
volumes they included finite-volume corrections in their analysis. Another
collaboration [8] used a relativistic heavy-quark action for charm quarks and
the clover-improved Wilson fermion action for the light and strange quarks in
order to calculate the charmed baryon spectrum. They used five ensembles with
pion masses down to $220\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ with lattice
spacings in the range $a\approx 0.12-0.06\,$\mathrm{f}\mathrm{m}$$ to
extrapolate to the continuum limit and they performed simultaneous
extrapolations to the physical light and strange quark masses. In Ref. [13], a
relativistic heavy-quark action is employed for the charm sector, while for
the light and strange quarks, they used domain-wall fermions. They used eight
ensembles having two different lattice spacings of $(a\approx 0.11\text{ and
}0.085\,$\mathrm{f}\mathrm{m}$)$ and of pion masses down to
$227\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ to extrapolate to the physical point
and to the continuum limit.

Strange baryons
---
$\Lambda(1.116)$ | $\Sigma(1.193)$ | $\Xi(1.314)$ | $\Delta(1.232)$ | $\Sigma^{\star}(1.385)$ | $\Xi^{\star}(1.530)$
1.1079(84) | 1.187(12) | 1.317(18) | 1.218(61) | 1.360(38) | 1.550(16)
Spin-$1/2$ Charmed baryons
$\Sigma_{c}(2.455)$ | $\Xi_{c}(2.470)$ | $\Xi_{c}^{\prime}(2.578)$ | $\Omega_{c}^{0}(2.695)$ | $\Xi_{cc}(3.622)$ | $\Omega_{cc}$
2.443(51) | 2.460(35) | 2.584(36) | 2.675(32) | 3.634(51) | 3.703(51)
Spin-$3/2$ Charmed baryons
$\Sigma_{c}^{\star}(2.520)$ | $\Xi_{c}^{\star}(2.645)$ | $\Omega_{c}^{\star}(2.770)$ | $\Xi_{cc}^{\star}$ | $\Omega_{cc}^{\star}$ | $\Omega_{ccc}$
2.509(46) | 2.645(55) | 2.747(31) | 3.676(55) | 3.803(50) | 4.785(71)
Table 10: Our final values for the low-lying baryon spectrum determined using
method I with model averaging. The experimental mass from the PDG [33], if
known, is given in the parenthesis next to the baryon symbol in the first row.

Very recently, the RQCD collaboration [31] used 58 $N_{f}=2+1$ ensembles of
clover-improved Wilson fermions with six different lattice spacings starting
from $a=0.10\,$\mathrm{f}\mathrm{m}$$ down to $a\sim
0.04\,$\mathrm{f}\mathrm{m}$$. These ensembles yielded pion masses ranging
from $429\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ down to
$127\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ and using them they performed a chiral
extrapolation to the physical point, infinite volume, and continuum
extrapolations to obtain the strange baryon spectrum.

In Figs. 16 and 17, we compare our results with those of the aforementioned
collaborations as well as with experimental results [33].

We observe excellent agreement between our results and previous lattice QCD
results. They are also in agreement with the experimental values when
available. We find that lattice QCD results on the mass of $\Xi^{*}$ favor a
larger central value by about 20 MeV.

##  IX Conclusions

Figure 18: Results on the strange (left) and charm (right) quark masses (blue
circles) compared to the FLAG21 average (black squares) [34] and our previous
results (red triangles) [18].

Using three $N_{f}=2+1+1$ ensembles simulated with physical values of the
quark masses, we perform for the first time directly at the physical pion mass
the continuum limit for all low-lying baryon masses. The volume for these
three ensembles is approximately the same with a spatial extent of 5.1 fm for
the cB211.072.64 ensemble and 5.5 fm for the other two ensembles. The lattice
spacings are determined using the nucleon mass after applying a correction for
the small miss-tuning in the pion mass. We use the twisted mass fermion
formulation with clover improvement and Osterwalder-Seiler valence strange and
charm quarks. We tune the valence quark masses using the physical masses of
the $\Omega^{-}$ and the $\Lambda_{c}^{+}$ baryons, respectively. The
renormalized strange and charm quark masses that we find are
$m_{s}(2~{}\text{GeV})=99.2(2.7)\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ and
$m_{c}(3~{}\text{GeV})=1.015(39)\,$\mathrm{G}\mathrm{e}\mathrm{V}$$,
respectively in the $\overline{\text{MS}}$ scheme at the continuum limit. We
compare these values in Fig. 18. There is agreement with the values of our
previous analysis [18] that did not include the cD211.054.96 ensemble but
included ensembles simulated at larger than physical pion mass that enabled us
to include ensembles with $a\sim 0.093$ fm. While there is agreement with the
FLAG average [34] for the charm quark mass, the one standard deviation tension
for the strange quark mass persists.

Our values for the baryon masses are in agreement with experimental results.
For the $\Xi^{\star}$ baryon, we find
$M_{\Xi^{\star}}=1.550(16)\,$\mathrm{G}\mathrm{e}\mathrm{V}$$, as compared to
the experimental value of 1.530 GeV and consistent with other two lattice QCD
determinations[32, 31]. Although within about a standard deviation the lattice
QCD results agree with the experimental value, all three results yield a
larger central value as compared to what PDG quotes. This overall agreement
with the experimental results allows us to predict the unmeasured baryon
masses of doubly and triply charmed baryons. The values predicted for the
masses of the doubly charmed $\Xi_{cc}^{\star}$, $\Omega_{cc}$ and
$\Omega_{cc}^{\star}$ baryons are 3.676(55)GeV, 3.703(51) GeV and
3.803(50)GeV, respectively, and for the triply charmed $\Omega_{ccc}$ baryon
4.785(71)GeV.

## Acknowledgments

C.A. acknowledges support by the project 3D-nucleon, id number
EXCELLENCE/0421/0043, co-financed by the European Regional Development Fund
and the Republic of Cyprus through the Research and Innovation Foundation.
G.C. was funded by the projects NextQCD, id EXCELLENCE/0918/0129, and
NiceQuarks id EXCELLENCE/0421/0195, co-financed by the European Regional
Development Fund and the Republic of Cyprus through the Research and
Innovation Foundation. S.B. is funded by the project QC4LGT, id number
EXCELLENCE/0421/0019, co-financed by the European Regional Development Fund
and the Republic of Cyprus through the Research and Innovation Foundation.
S.B. also acknowledges funding by the EuroCC project (grant agreement No.
951740). This project used computer time under the project pr74yo on the
JUWELS Booster system at the Jülich Supercomputing Centre and on the Cyclone
supercomputer at The Cyprus Institute. The authors gratefully acknowledge the
Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this
project by providing computing time through the John von Neumann Institute for
Computing (NIC) on the GCS Supercomputer JUWELS-Booster at Jülich
Supercomputing Centre (JSC).

## References

  * Lutz  _et al._ [2009] M. F. M. Lutz _et al._ (PANDA), (2009), arXiv:0903.3905 [hep-ex] .
  * Wang [2010] Z.-G. Wang, Eur. Phys. J. A 45, 267 (2010), arXiv:1001.4693 [hep-ph] .
  * Martynenko [2008] A. P. Martynenko, Phys. Lett. B 663, 317 (2008), arXiv:0708.2033 [hep-ph] .
  * Ebert  _et al._ [2002] D. Ebert, R. N. Faustov, V. O. Galkin, and A. P. Martynenko, Phys. Rev. D 66, 014008 (2002), arXiv:hep-ph/0201217 .
  * Roberts and Pervin [2008] W. Roberts and M. Pervin, Int. J. Mod. Phys. A 23, 2817 (2008), arXiv:0711.2492 [nucl-th] .
  * Na and Gottlieb [2007] H. Na and S. A. Gottlieb, PoS LATTICE2007, 124 (2007), arXiv:0710.1422 [hep-lat] .
  * Na and Gottlieb [2008] H. Na and S. Gottlieb, PoS LATTICE2008, 119 (2008), arXiv:0812.1235 [hep-lat] .
  * Briceno  _et al._ [2012] R. A. Briceno, H.-W. Lin, and D. R. Bolton, Phys. Rev. D 86, 094504 (2012), arXiv:1207.3536 [hep-lat] .
  * Liu  _et al._ [2010] L. Liu, H.-W. Lin, K. Orginos, and A. Walker-Loud, Phys. Rev. D 81, 094505 (2010), arXiv:0909.3294 [hep-lat] .
  * Basak  _et al._ [2012] S. Basak, S. Datta, M. Padmanath, P. Majumdar, and N. Mathur, PoS LATTICE2012, 141 (2012), arXiv:1211.6277 [hep-lat] .
  * Durr  _et al._ [2012] S. Durr, G. Koutsou, and T. Lippert, Phys. Rev. D 86, 114514 (2012), arXiv:1208.6270 [hep-lat] .
  * Namekawa  _et al._ [2013] Y. Namekawa _et al._ (PACS-CS), Phys. Rev. D 87, 094512 (2013), arXiv:1301.4743 [hep-lat] .
  * Brown  _et al._ [2014] Z. S. Brown, W. Detmold, S. Meinel, and K. Orginos, Phys. Rev. D 90, 094507 (2014), arXiv:1409.0497 [hep-lat] .
  * Pérez-Rubio  _et al._ [2015] P. Pérez-Rubio, S. Collins, and G. S. Bali, Phys. Rev. D 92, 034504 (2015), arXiv:1503.08440 [hep-lat] .
  * Alexandrou and Kallidonis [2017] C. Alexandrou and C. Kallidonis, Phys. Rev. D 96, 034511 (2017), arXiv:1704.02647 [hep-lat] .
  * Finkenrath  _et al._ [2022] J. Finkenrath _et al._ , PoS LATTICE2021, 284 (2022), arXiv:2201.02551 [hep-lat] .
  * Abdel-Rehim  _et al._ [2017] A. Abdel-Rehim _et al._ (ETM), Phys. Rev. D 95, 094515 (2017), arXiv:1507.05068 [hep-lat] .
  * Alexandrou  _et al._ [2021] C. Alexandrou _et al._ (Extended Twisted Mass), Phys. Rev. D 104, 074515 (2021), arXiv:2104.13408 [hep-lat] .
  * Alexandrou  _et al._ [2018] C. Alexandrou _et al._ , Phys. Rev. D 98, 054518 (2018), arXiv:1807.00495 [hep-lat] .
  * Finkenrath [2023] J. Finkenrath, PoS LATTICE2022, 227 (2023).
  * Frezzotti  _et al._ [2000] R. Frezzotti, P. A. Grassi, S. Sint, and P. Weisz, Nucl. Phys. B Proc. Suppl. 83, 941 (2000), arXiv:hep-lat/9909003 .
  * Sheikholeslami and Wohlert [1985] B. Sheikholeslami and R. Wohlert, Nucl. Phys. B 259, 572 (1985).
  * Alexandrou  _et al._ [2023] C. Alexandrou _et al._ (Extended Twisted Mass), Phys. Rev. D 107, 074506 (2023), arXiv:2206.15084 [hep-lat] .
  * Iwasaki [1983] Y. Iwasaki, (1983), arXiv:1111.7054 [hep-lat] .
  * Osterwalder and Seiler [1978] K. Osterwalder and E. Seiler, Annals Phys. 110, 440 (1978).
  * Benmerrouche  _et al._ [1989] M. Benmerrouche, R. M. Davidson, and N. C. Mukhopadhyay, Phys. Rev. C 39, 2339 (1989).
  * Gusken [1990] S. Gusken, Nucl. Phys. B Proc. Suppl. 17, 361 (1990).
  * Alexandrou  _et al._ [1994] C. Alexandrou, S. Gusken, F. Jegerlehner, K. Schilling, and R. Sommer, Nucl. Phys. B 414, 815 (1994), arXiv:hep-lat/9211042 .
  * Jay and Neil [2021] W. I. Jay and E. T. Neil, Phys. Rev. D 103, 114502 (2021), arXiv:2008.01069 [stat.ME] .
  * Neil and Sitison [2022] E. T. Neil and J. W. Sitison, (2022), arXiv:2208.14983 [stat.ME] .
  * Bali  _et al._ [2022] G. S. Bali, S. Collins, P. Georg, D. Jenkins, P. Korcyl, A. Schäfer, E. E. Scholz, J. Simeth, W. Söldner, and S. Weishäupl (RQCD), (2022), arXiv:2211.03744 [hep-lat] .
  * Durr  _et al._ [2008] S. Durr _et al._ , Science 322, 1224 (2008), arXiv:0906.3599 [hep-lat] .
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Aoki  _et al._ [2022] Y. Aoki _et al._ (Flavour Lattice Averaging Group (FLAG)), Eur. Phys. J. C 82, 869 (2022), arXiv:2111.09849 [hep-lat] .


∎

e1e-mail<EMAIL_ADDRESS>e2e-mail: :
<EMAIL_ADDRESS>
11institutetext: Al Ain University 22institutetext: Physics Dept., Damascus
University, Damascus, Syria

# $F(R,w)Gravity$

A new gravity framework

Mahmoud AlHallak e1,e2,addr1,addr2

(Received: date / Accepted: date)

###### Abstract

In this work we present a new framework of the gravity sector by considering
the extension $F(R,w)$, in which $R$ is the Ricci scalar and $w$ is the
equation of state. Three different choices of function $F(R,w)$ are
investigated under the Palatini formalism. The models appear equivalent to
$F(R)$ models of gravity with effective momentum-energy tensors. For linear
dependence of Ricci scalar in which $F(R,w)=k(w)R$, the model appears
equivalent to Einstein-Hilbert action with effective momentum-energy tensor.
Recovering the minimal coupling case of the last choice does not face Jordan-
Einstein frame ambiguities and exhibits natural alignments with general
relativity results in the matter/ radiation dominated eras. We discuss some
astrophysical implications of the model by considering scalar fields as
dominant matter forms. We show that the Higgs inflation could be saved within
the $F(R,w)$ model. We suggest some future investigations exemplified by
constant-roll inflation and universe evolution for $F(R)=f(R)k(w)$ where
$f(R)$ represents the Starobinsky gravitational form. Using the model and
comparing it with pure $F(R)$ gravity, we provide preliminary indications of
$F(R,w)$’s impact. As a final note, we suggest using the Polytropic equation
of state in future works to investigate $F(R,w)$.

###### Keywords:

Modified gravity Palatini formalism Cosmological inflation

##  1 Introduction:

The theory of general relativity GR Einstein:1916vd  has made significant
contributions to cosmology, which cannot be ignored. The theory of general
relativity alters our understanding of space-time and its relationship to
energy in some way that we cannot avoid in any future developments of our
knowledge of the universe.

In spite of the previous notes, GR, under its original minimal metric
formalism, failed to explain some recent cosmological observations, such as,
the late expansion of the universe measured from the light-curves of several
hundred type Ia supernovae SupernovaSearchTeam:1998fmf ;
SupernovaCosmologyProject:1998vns .
Several interesting attempts are made in the scientific community to modify GR
in order to solve and explain these cosmological puzzles.
$F(R)$ gravity Carroll:2003wy ; Hu:2007nk  is the most common modification, in
which the gravitational sector is represented by higher-order curvature. Many
cosmological features are examined within F(R) gravity, such as in Liu:2023ypa
; Odintsov:2023nqe ; Odintsov:2023cli ; AlHallak:2022gbv ; AlHallak:2021hwb ;
AlHallak:2023wsx .
One of the interesting extensions of $F(R)$ models is assuming a non-minimal
coupling between the curvature represented by Ricci scalar and matter
represented by the momentum-energy tensor $T$. This class of models is known
as $F(R,T)$ gravity Harko:2011kv .
Models of this category are investigated both from a metric perspective
Harko:2011kv  as well as a Palatini perspectiveWu:2018idg .
A wide range of cosmological aspects have been studied in $F(R,T)$ gravity. In
this section, we will mention some of these studies that can serve as
motivation for performing similar studies in the new model of $F(R,w)$
gravity.

A study of wormholes that involves electric charge effects is reported in
Azmat:2023ygn . Authors use $F(R,T)$ gravity to investigate the existence and
viability of wormholes. A wormhole shape function is obtained, the fundamental
wormhole conditions are analyzed, and the validity of the null energy
condition is examined. In addition, the authors analyze the gravitational
active mass of wormhole solutions based on their electric charge dynamics.
These results demonstrate that the obtained wormhole solutions exhibit valid
behavior and that the values of the complexity factor are always within
reasonable limits.
Other studies pertaining to wormholes in $F(R,T)$ can be found in
Moraes:2019pao ; Moraes:2017rrv ; Moraes:2016akv ; Sahoo:2017ual ;
Bhar:2021lat .

A study of cosmological inflation is carried out in F(R,T) gravity. Ref
Gamonal:2020itt  investigates slow-roll inflation using an explicit form of
$F(R,T)=R+2\kappa\alpha T$. Authors analyze observable quantities represented
by the spectral index $n_{s}$ and tensor-to-scalar ratio $r$. They computed
these quantities for a general power-law potential, natural and quasi-quartic
Hilltop inflation, and a Starobinsky model.
The ref Deb:2022hna  investigates the Double-Well potential using the same
function of $F(R,T)$ and analyzes $n_{s}$ and $r$ as the primary observables.
In Bhattacharjee:2022lcs  authors employ mimetic $F(R,T)$ gravity coupled with
Lagrange multiplier and mimetic potential to yield viable inflationary
cosmological solutions consistent with latest Planck and BICEP2/Keck Array
data.

As in Singh:2020gxd  and Sun:2015yga , gravity extension $F(R,T)$ has been
explored as a possible explanation for late accelerating expansion of the
universe.
The dark energy models are also investigated in $F(R,T)$ such as in
Maurya:2022pzw  where a perfect dark fluid is examined in Logarithmic $F(R,T)$
Cosmology with $F(R,T)=R-\kappa\alpha\ln(T)$.
Holographic dark energy has been the subject of several studies, as can be
seen in Sharma:2022dzc ; Shaikh:2022ynt ; Varshney:2021xvg ; Shaikh:2019ppk .

A second interesting study in this field is the $F(R,L_{m})$ gravity
Harko:2010mv . In a similar manner to $F(R,T)$ gravity, various cosmological
phenomena are investigated in $F(R,L)$. Studies on wormholes can be found in
Naseer:2023rab ; Kavya:2023tjf ; Kavya:2023mwv . Investigations into dark
energy in the context of$F(R,L_{m})$ are discussed in Maurya:2023lxt ;
Pradhan:2022msm .

A recent proposal by Dioguardi et al. Dioguardi:2023jwa  introduces a model of
$F(R,X)$ to examine the inflationary expansion of the universe, where $X$
represents the kinetic energy of the inflaton field.

In the current work, we introduce a novel gravity framework by considering an
extended form of gravity denoted $F(R,w)$, where $w$ represents the equation
of state. The model exhibits a notable characteristics: The model coincides
naturally with General Relativity during epochs in which matter and radiation
are dominant. A closer look at specific forms of $F(R,w)$ shows changes in the
scales of some cosmological quantities such as the Hubble parameter and the
cosmological acceleration of the universe, but the functions of these
quantities remain unchanged. This statement, however, needs further
investigation before it can be generalized. The model also has a significant
impact when scalar fields become dominant during early stages of the evolution
of the universe (inflaton fields) or during the later expansion phase
(quintessence). Unlike common models of non-minimal coupling, for
$F(R,w)=R+\alpha Rw$ under Palatini formalism, the model appears to be free of
Jordan-Einstein frame ambiguity since it is reduced to the minimal case
without the need for conformal transformation.

The paper is organized as follows. In section. 2, $F(R,w)$ is presented in its
most general form, analyzing it both in the metric and Palatini formalisms.
Section. 3 discusses specific cases of model $F(R,w)$ and derivations of field
equations and Friedmann equations based on FRW flat metric. Further, we
discuss the reduction of the considered models to the minimal coupling case.
Finally, In section. 4, the astronomical implications of the model are
illustrated through an examination of the inflaton field as dominating matter
during early universe evolution. We consider the Higgs inflationary potential,
and we show how it can be saved within $F(R,w)$ gravity. As a follow-up, we
provide some comments for future studies exemplified by constant-roll
inflation, the Starobinsky model, and polytropic equations of state, and we
provide some insights into how $F(R,w)$ impacts them. We end up with a summary
and conclusion in section. 5.

##  2 The model:

As a first step, we will discuss the most general case of $F(R,w)$ gravity, in
which the action is expressed as follows:

$S=\frac{1}{16\pi G}\int
F(\mathcal{R},w)\sqrt{-g}d^{4}x+\int\mathcal{L}_{M}\sqrt{-g}d^{4}x$ (1)

Throughout this section, we discuss the the field equations in both cases of
metric and Palatini formalism, where in the case of metric formalism, the
connection is constructed from the metric tensor through the Levi-Civita
connection $\Gamma^{\lambda}_{\mu\nu}=\\{^{\lambda}\>_{\mu\nu}\\}$, while in
the Palatini formalism, the connection and the metric are considered
independent.

According to the metric formalism, by taking the variation of the action (1)
with respect to the metric tensor, the modified field equations take the
following form:

$\begin{split}F_{\mathcal{R}}(\mathcal{R},w)\mathcal{R}_{\mu\nu}-\frac{1}{2}F(\mathcal{R},w)g_{\mu\nu}+(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu})\\\
F_{\mathcal{R}}(\mathcal{R},w)=8\pi GT_{\mu\nu}\end{split}$ (2)

By contracting equation (2), the relation between the Ricci scalar
$\mathcal{R}$ and the trace $T$ of the stress-energy tensor is as follows:

$F_{\mathcal{R}}(\mathcal{R},w)\mathcal{R}+3\Box
F_{\mathcal{R}}(\mathcal{R},w)-2F(\mathcal{R},w)=8\pi GT$ (3)

By eliminating the term $\Box F_{\mathcal{R}}(\mathcal{R},w)$between equations
(2) and (3), the gravitational field equations can be written in the form:

$\begin{split}F_{\mathcal{R}}(\mathcal{R},w)(\mathcal{R}_{\mu\nu}-\frac{1}{3}\mathcal{R}g_{\mu\nu})+\frac{1}{6}F(\mathcal{R},w)g_{\mu\nu}=\\\
8\pi G(T_{\mu\nu}-\frac{1}{3}Tg_{\mu\nu})\end{split}$ (4)

In the Palatini formalism, the gravitational action (1) is now varied with
respect to the metric tensor $g_{\mu\nu}$ under the assumption $\delta
R_{\mu\nu}=0$.
As a result we obtain the field equations:

$F_{R}(R,w)R_{\mu\nu}-\frac{1}{2}F(R,w)g_{\mu\nu}=8\pi GT_{\mu\nu}$ (5)

The expression for the Ricci scalar is obtained by contracting the equation
(5) with $g_{\mu\nu}$ as,

$RF_{R}(R,w)+2F(R,w)=8\pi GT$ (6)

In terms of the Einstein tensor one can write:

$\begin{split}G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\frac{1}{F_{R}(R,w)}\\\
\bigg{\\{}8\pi G\
(T_{\mu\nu}-\frac{g_{\mu\nu}}{2}T)-\frac{g_{\mu\nu}}{2}F(R,w)\bigg{\\}}\end{split}$
(7)

This paper focuses on palatini formalism due to its interesting
characteristics in $F(R,w)$ models, whereas future studies will examine the
model in metric formalism.

##  3 Particular Cases:

This section discusses some possible choices regarding the function $F(R,w)$.
There are many more forms that could be assumed and studied, but we only
provide three examples here.

###  3.1 $F(R,w)=f(R)k(w)$

Let’s consider the following particular function of $F(R,w)$ gravity as:

$F(R,w)=f(R)k(w)$ (8)

Taking the variation of action (1) with respect to the metric $g_{\mu\nu}$, we
can express the field equation as follows:

$f_{R}R_{\mu\nu}-\frac{1}{2}f(R)g_{\mu\nu}=8\pi GT^{eff}_{\mu\nu}$ (9)

Where:

$T^{eff}_{\mu\nu}=\frac{T_{\mu\nu}}{k(w)}$ (10)

Interestingly, the model can be reduced to an $f(R)$ model with an effective
momentum-energy tensor.
Taking the variation of the action with respect to connection we can find:

$\nabla_{\alpha}(\sqrt{-g}f_{R}g^{\mu\nu})=0$ (11)

In order to solve equation (11), it is convinent to introduce a new metric:

$\sqrt{h}h_{\mu\nu}=\sqrt{-g}f_{R}g_{\mu\nu}$ (12)

The connection for this metric is a Levi-Civita connection. so both metrics
are related by the conformal factor.

$h_{\mu\nu}=f_{R}g_{\mu\nu}$ (13)

For this reason one should assume that the conformal factor $f_{R}\neq 0$.
Considering perfect fluid $T_{\mu}^{\nu}=diag(-\rho,p,p,p)$ we can find the
effective equation of state equivalent to the standard one
$w(\text{eff})=\frac{p_{\text{eff}}}{\rho_{\text{eff}}}=w=\frac{p}{\rho}$
where, $\rho_{\text{eff}}=\rho/K(w)$ and $p_{\text{eff}}=p/K(w)$.

As long as we are in the Jordan frame, we can continue working on the field
equations (9). Alternatively, equations can also be displayed in the Einstein
frame, where gravitational action is expressed formally in Hilbert-Einstein
form.
In the Jordan frame, action (1) is dynamically equivalent to the first order
Palatini gravitational action:

$\begin{split}S=\frac{1}{16\pi G}\int
d^{4}x\sqrt{-g}\big{[}\frac{1}{2}(f(\zeta)+f^{\prime}(\zeta)(R-\zeta))\big{]}+\\\
\int d^{4}x\sqrt{-g}L_{M}^{eff}\end{split}$ (14)

with $f^{\prime\prime}(R)\neq 0$.
The action (14) takes the Brans-Dicke model form by introducing an auxiliary
scalar field $\Phi=f^{\prime}(\zeta)$ and taking into consideration the
constraint $\zeta=R$.

$S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}\big{[}\frac{1}{2}(\Phi
R-U(\Phi)\big{]}+\int d^{4}x\sqrt{-g}\mathcal{L}_{M}^{eff}$ (15)

where $U(\Phi)=\zeta(\Phi)\Phi-f(\zeta(\Phi))$.
The action (15) can be reduced to the minimal coupling case by applying a
conformal transformation $g_{\mu\nu}^{E}=\Phi g_{\mu\nu}^{J}$:

$S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g_{E}}(R_{E}-U_{E}(\Phi))+\int
d^{4}x\sqrt{-g_{E}}\mathcal{L}_{M(E)}^{eff}$ (16)

Where $U_{E}=U/\Phi^{2}$.
In FRW metric, Friedmann equations take the form:

$H_{E}^{2}=\frac{8\pi G}{3}\bigg{(}\rho^{E}_{\Phi}+\rho^{E}_{M}\bigg{)}$ (17)
$(\dot{H}+H^{2})_{E}=-\frac{4\pi
G}{3}\big{(}2\rho^{E}_{\Phi}-\rho^{E}_{M}+3P^{E}_{M}\big{)}$ (18)

Subscript E denotes the evaluation of all quantities in the Einstein frame,
with $\rho^{E}_{M}=\rho_{M}/\Phi k(w)$ and $P^{E}_{M}=P_{M}/\Phi k(w)$.

###  3.2 $F(R,w)=R+\alpha wR$

In a specific case of $f(R)$, linear dependence of the Ricci scalar would be
considered as:

$F(R,w)=R+\alpha wR$ (19)

The choice (19) appears to be a small perturbation to the theory for
$\alpha\ll 1$.
Interesting features can be found in this choice. First of all, The model
coincides naturally with general relativity during the matter and radiation
era, in which $w=0$ and $w=1/3$ respectively. Second, even with the non-
minimal coupling with matter through the equation of state, the model appears
as an effective energy-momentum theory with minimal coupling. Thus, the model
is free from Einstein-Jordan frame ambiguity.
The field equations are given as:

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=8\pi GT^{eff}_{\mu\nu}$ (20)

with

$T^{eff}_{\mu\nu}=\frac{T_{\mu\nu}}{(1+\alpha w)}$ (21)

Friedmann equations can be shown as:

$H^{2}=\frac{8\pi G}{3}\bigg{(}\frac{\rho^{2}}{\rho+\alpha p}\bigg{)}$ (22)
$\dot{H}+H^{2}=-\frac{4\pi G}{3}\rho\bigg{(}\frac{\rho+3p}{\rho+\alpha
p}\bigg{)}$ (23)

This class of models has great influence when considering scalar fields as the
dominant matter in the universe to explain early or late inflationary
expansions.

###  3.3 $F(R,w)=f(R)+k(w)$

As a third case of generalized $F(R,w)$ gravity models, we consider that the
action is given by,

$F(R,w)=f(R)+k(w)$ (24)

For an arbitrary matter source the gravitational field equations are given by:

$f_{R}R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}f(R)=8\pi G(T_{\mu\nu}+K_{\mu\nu})$ (25)

where:

$K_{\mu\nu}=\frac{k(w)}{8\pi G}g_{\mu\nu}$ (26)

Once again, we can recover the minimal coupling case by applying the conformal
transformation $g_{\mu\nu}^{E}=\Phi g_{\mu\nu}^{J}$.
The Friedman equations take the form,

$H_{E}^{2}=\frac{8\pi G}{3}(\rho_{\Phi}+\rho_{M}+\rho_{w})_{E}$ (27)
$(\dot{H}+H^{2})_{E}=-\frac{4\pi
G}{3}(2\rho_{\Phi}-\rho_{M}+3(p_{M}+p_{w}))_{E}$ (28)

where $\rho_{w}=\frac{k(w)}{8\pi G}$, and $p_{w}=\frac{k(w)}{8\pi G}$

##  4 Astrophysical implications:

In order to understand the astrophysical implications of this gravitational
extension, the function $F(R,w)$ needs to be defined specifically. In this
section, we examine the slow-roll inflationary scenario in the simplest case
to demonstrate how $F(R,w)$ gravity impacts the main observable quantities
represented by $n_{s}$ and $r$. Following that, we present three possible
future research directions. Firstly, we provide a brief look at constant-roll
inflation in which the $F(R,w)$ have a greater impact on the model. Next, we
inspect the influence of $F(R,w)$ on Starobinsky’s model, and finally we
discuss the generalized polytropic equation of state.

###  4.1 Slow-roll Inflation:

Using scalar fields as dominated matter, this section examines the early
accelerated expansion of the universe.
Let’s consider a canonical scalar field as follows:

$L_{\phi}=-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-V(\phi)$ (29)

with signature $(-,+,+,+)$.
Following the results of section (3.2) and the function (19) , we find the
following:

$S=\int\sqrt{-g}d^{4}x\bigg{[}\frac{1}{2}R+-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\varphi\nabla_{\nu}\varphi-U(\phi(\varphi))\bigg{]}$
(30)

Where the re-scaled canonical field $\varphi$ is given by:

$\frac{d\chi}{d\phi}=\frac{1}{\sqrt{(1+\alpha w)}}$ (31)

The potential $U$ is given by:

$U(\phi(\varphi))=\frac{V(\phi)}{(1+\alpha w)}$ (32)

As a result we can find Friedmann equation in Jordan frame as:

$H^{2}=\frac{1}{2}\dot{\chi}^{2}+U(\chi)$ (33)
$\dot{H}=-\frac{1}{2}\dot{\chi}$ (34)

The equation of motion of the scalar field:

$\ddot{\chi}+3H\dot{\chi}+U_{\chi}=0$ (35)

Let us now consider the most extreme case where the potential energy dominates
the kinetic energy $U\gg X$, so $w\sim-1$. then the slow-roll approximation
for constant $\alpha$ is given as:

$\varepsilon=\frac{(1-\alpha)}{2}\bigg{(}\frac{V_{\phi}}{V}\bigg{)}^{2}$ (36)
$\eta=(1-\alpha)\bigg{(}\frac{V_{\phi\phi}}{V}\bigg{)}$ (37)

Thus, the spectral index $n_{s}$ and tensor to scalar ratio $r$ are
represented as:

$n_{s}-1=(1-\alpha)(2\eta_{\phi}-6\varepsilon_{\phi})$ (38)
$r=16(1-\alpha)\varepsilon_{\phi}$ (39)

where $\varepsilon_{\phi}=\frac{1}{2}\bigg{(}\frac{V_{\phi}}{V}\bigg{)}^{2}$,
and $\eta_{\phi}=\bigg{(}\frac{V_{\phi\phi}}{V}\bigg{)}$.
Consider Higgs inflation with the potential given by:

$V(\phi)=\lambda(\phi^{2}-\nu)^{2}$ (40)

The potential takes the form of a quartic potential for $\nu\ll\phi$.
Modified gravity models such as Pozdeeva:2020apf ; Tahmasebzadeh:2016zhd ;
AlHallak:2023tca  take into account potentials of this type. In this
section,However we consider it in $F(R,w)$.
The observables $n_{s}$ and $r$ take the following forms:

$n_{s}=(1-\alpha)\bigg{(}1-\frac{24}{\phi^{2}}\bigg{)}$ (41)
$r=(1-\alpha)\bigg{(}\frac{128}{\phi^{2}}\bigg{)}$ (42)

Although the value of $n_{s}$ and $r$ will be affected by the factor
$(1-\alpha)$, the slope of the line representing the relationship between them
remains unaffected ($\simeq 5.33$), so observational results are still exempt
from the Higgs inflation, as we can see in Figure (1).
In this extremely simplified case of $w$, we can still see the benefits of the
model by considering the parameter $\alpha(\phi)$. It is easy to see the
effect by examining the term $(1-\alpha(\phi))$ as follows:

$(1-\alpha(\phi))=(1-\xi\phi^{2})^{2}$ (43)

As can be seen in this case, the model is equivalent numerically to the non-
minimal coupling case in Einstein frame Linde:2011nh , but again, the
treatment is performed in Jordan frame.
The non-minimal coupling to gravity has been extensively studied within
different gravity formalisms, such as Linde:2011nh ; Bauer:2008zj ;
AlHallak:2021sic ; AlHallak:2016wsa ; Karamitsos:2021mtb ; Hertzberg:2010dc ;
Qiu:2010dk .
Figure (1) showcases the theoretical results of the quartic potential within
the context of the minimal coupling case in standard GR in contrast to the
model within $F(R,w)$ gravity. The black dashed line represents the
theoretical results obtained in the minimal coupling case. It is evident that
the model is excluded since there is no intersection between the theoretical
results and the observational contours. The black dot corresponds to the
number of e-foldings $N=60$. The green dashed line corresponds to the results
obtained in the $F(R,w)$ framework, taking into account the equation (43). The
line corresponds to the value of N=60, with $\xi$ increasing from zero at the
top to infinity at the bottom of the graph. The results of the inflationary
model exhibit noteworthy agreement with the observational outcomes at the
$95\%$and $68\%$ confidence levels of the combined results of (Planck TT, TE,
EE + lowE + lensing), represented as ”Green Contours,” along with the results
of (Planck TT, TE, EE + lowE + lensing + BK15), denoted as ”Blue Contours.”
Keeping in mind that $\alpha$ can take different forms, one can expect and
assess how various forms will impact inflationary results. In addition, we can
consider the more natural case where $w=w(\phi)$, as in some inflationary
scenarios, such as warm inflation, in which the inflaton field continuously
dissipates into radiation during the inflationary period, or constant-roll
inflation, where the kinetic term take place during the inflationary period.

Figure 1: A plot of $n_{s}$ and $r$ in the case of the quartic potential . The
black dashed line depicts the inflationary model within the framework of
standard General Relativity (GR), while the green dashed line represents the
model in F(R,w) gravity. The experimental contours are generated by combining
the results of (Planck TT, TE, EE + lowE + lensing) denoted as ”Green
Contours,” along with the results of (Planck TT, TE, EE + lowE + lensing +
BK15) indicated as ”Blue Contours” Planck:2018jri  .

###  4.2 Future Directions and Expected Implications

In this section we provide some suggestions for future works to investigate
important cosmological phenomena in the field of $F(R,w)$. We present here,
some comments about this suggestions leaving the detailed investigations for
future works.

####  4.2.1 Constant-roll inflation:

The slow roll inflationary scenario was discussed in the previous section
quickly to illustrate how $F(R,w)$ models of gravity impact cosmological
inflation.
However, in the slow-roll inflationary scenario the kinetic energy term is
negligible in front of the potential energy and the equation of state is
approximately $-1$. As discussed earlier, $F(R,w)$ affects the inflationary
model under the slow-roll approximation because of the coupling coefficient
$\alpha$, which may be constant or depend on the scalar field.

In this section, we discuss constant-roll inflation, in which kinetic energy
has a significant role to play.
Specifically, we consider $F(R,w)=k(w)R=(R+\alpha Rw)$.
Accordingly, Friedmann equations will be expressed in terms of $\phi$ as
follows:

$3H^{2}=\frac{1}{k(w)}[\frac{1}{2}\dot{\phi}^{2}+V(\phi)]$ (44)
$-2\dot{H}=\frac{1}{k(w)}\dot{\phi}^{2}$ (45)

The scalar field’s equation of motion is:

$\ddot{\phi}+3H\dot{\phi}+k(w)U_{\phi}=0$ (46)

Where $U_{\phi}$ represents the derivative of the potential $U$ with respect
to scalar field $\phi$. A minimum ground state $\dot{w}=0$ is assumed in the
previous equations.
As part of the original work of the constant-roll inflation Motohashi:2014ppa
, authors used the Hamiltonian-Jacobi formalism to extract the Hubble
parameter from the inflaton field $H=H(\phi)$.
In this paper we follow same argument which is applicable as long as
$t=t(\phi)$ is a single-valued function, i.e $\dot{\phi}\neq 0$.
The constant rate of rolling is parameterized by the equation:

$\ddot{\phi}=-(3+\beta)H\dot{\phi}$ (47)

As shown in ref Motohashi:2014ppa , $H=H(\phi)$ is independent of the
inflationary model represented by $V(\phi)$, and its general solution can be
formulated as follows:

$H(\phi)=C_{1}\exp\bigg{(}\sqrt{\frac{3+\beta}{2}}\frac{\phi}{M_{pl}}\bigg{)}+C_{2}\exp\bigg{(}-\sqrt{\frac{3+\beta}{2}}\frac{\phi}{M_{pl}}\bigg{)}$
(48)

Conversely, we can find in the $F(R,w)$ from equation (45) that:

$-2M_{pl}^{2}H_{\phi}=\frac{1}{k(w)}\dot{\phi}$ (49)

We can determine $\ddot{\phi}$ from the time derivatives of equation (49).
The differential equation for the Hubble parameter $H$ can be found by
replacing $\ddot{\phi}$ in equation (47) as,

$H_{\phi\phi}=\frac{(3+\beta)}{2M_{pl^{2}}k(w)}H$ (50)

Solving equation (44) , one can find $\dot{\phi}(\phi)$ as,

$\dot{\phi}^{2}=3\alpha M_{pl}^{2}H^{2}+\sqrt{3\alpha M_{pl}^{2}H^{2}(3\alpha
M_{pl}^{2}H^{2}-8V)-2V}$ (51)

In the case of quartic potential $V=\lambda\phi^{4}$, we solved the
differential equation (50) numerically.
In figure (2), we plot $H$ as a function of $\phi$ in both standard GR
Motohashi:2014ppa  and $F(R,w)$ for comparison.
As shown in figure (2), the Hubble parameter’s function in GR and $F(R,w)$ are
exactly the same for $\alpha=0$, as it expected theoretically.
Here is a numerical example where each parameter value of the model is fixed
to unity, so that $H$ is represented in figure (2). Given $\phi=4$, we found
$\frac{H_{F(R,w)}}{H_{\text{GR}}}=0.2$ for $\alpha=1$. The ratio
$\frac{H_{F(R,w)}}{H_{\text{GR}}}$ depends on $\alpha$, where for $\alpha=2$
results in a ratio $0.1$, while $\alpha=3$ results in a ratio $0.07$.
As a result, we can conclude that the model $F(R,w)$ has an impact on the
evolution of the universe.
On the other hand, when checking the effect of the energy scale $\lambda$ of
the potential $V$, we found that there is no effect at all. Accordingly, the
question arises as to whether the constant-roll inflation in $F(R,w)$ gravity
is independent of the potential, similar to general Relativity, or not. This
conclusion, however, requires more investigation, which can be carried out in
future studies. All we are looking at here is the estimated inflationary
impact of the $F(R,w)$ model.

Figure 2: The Hubble parameter is plotted as a function of the scalar field
for different $\alpha$ values in $F(R,w)$ gravity with $F(R,w)=R+\alpha Rw$
compared to the GR model. A flatter Hubble parameter is seen as the parameter
$\alpha$ increases.

####  4.2.2 Starobinsky cosmological model:

Assuming the function $F(R,T)=f(R)k(w)=(R+\gamma R^{n})(1+\alpha w)$, then the
trace of equation (9) is shown as:

$R=\gamma(n-2)R^{n}=-T^{\text{eff}}$ (52)

For $\gamma\gg 1$ and $n\neq 1$ we can write (52) as:

$R=\bigg{(}-\frac{T^{\text{eff}}}{\gamma(n-2)}\bigg{)}^{\frac{1}{n}}$ (53)

In the Starobinsky model, where $f(R)=R+\gamma R^{2}$, the Ricci scalar is:

$R=-T^{\text{eff}}=\frac{(3P-\rho)\rho}{\alpha P}$ (54)

To simplify the analytical equation, we considered $\alpha w\gg 1$.
Considering $\Lambda$CDM model with flat FRW, we can derive Fiendmann equation
as follows:

$\begin{split}H^{2}=\big{(}\frac{b^{2}}{b+d/2}\big{)}^{2}\frac{3(\rho_{m}+\rho_{r}+\rho_{\Lambda})}{\alpha(\rho_{r}-3\rho_{\Lambda})}\\\
\bigg{[}\frac{3(\rho_{m}+\rho_{r}+\rho_{\Lambda})}{\alpha(\rho_{r}-3+\rho_{\Lambda})}(\rho_{m}+\rho_{\Lambda})^{2}\\\
\frac{(K-3)(K+1)}{2b}+(\rho_{m}+\rho_{\Lambda})+\rho_{r}/b\bigg{]}\end{split}$
(55)

where:

$K=\frac{3\rho_{\Lambda}}{(\rho_{m}+\rho_{\Lambda})}$ (56)
$b=1+\frac{6\gamma}{\alpha}(\rho_{m}+4\rho_{\Lambda})\frac{(\rho_{m}+\rho_{r}+\rho_{\Lambda})}{(\rho_{r}-3\rho_{\Lambda})}$
(57) $d=\frac{\dot{b}}{H}$ (58)

A numerical solution of equation (55) is used to determine the factors of
universe expansion represented by scale factor and Hubble parameter. These
results are presented along with the results within the Starobinsky model in
figures (3) and (4).
The numerical check proceeds by fixing unit values for all parameters except
the parameter $\beta=4$, and assuming current value of scale factor equals to
one. Figure (3) shows two solutions for the scale factor in $F(R,w)$ ”dashed
lines”, of which the ”Brown curve” corresponds to a negative Hubble’s
parameter solution, as seen in figure (4). Starobinsky’s pure $F(R)$ model
also has similar solutions, although with different amplitudes.
A significant impact is expected on the grand scheme of evolution of the
universe as a whole as a result of $F(R,w)$. Nevertheless, a detailed study
should be carried out in the future.

Figure 3: The scale factor for $F(R,w)=f(R)(1+\alpha w)$ is represented by
dashed lines compared with the pure Starobinsky model, which represents
$f(R)=R+\gamma R^{2}$ by solid lines. Figure 4: The Hubble parameter for
$F(R,w)=f(R)(1+\alpha w)$ is represented by dashed lines in contrast to the
pure Starobinsky model, which represents $f(R)=R+\gamma R^{2}$ by solid lines.

####  4.2.3 Polytropic equation of state:

All of the last discussions focused on the perfect fluid’s equation of state
as $P=\kappa\rho$, where $\kappa=0$ represents dust matters, $\kappa=1/3$
corresponds to radiation, and $\kappa=0$ represents vacuum energy. In order to
understand different cosmological aspects, different equations of state may be
required, such as a polytropic equation of state in which
$P=\xi\rho^{\sigma}$. where $\xi$ and $\sigma$ are constants.$\sigma$ usually
is written as, $\sigma=1-1/n$, where $n$ is known as polytropic index.
Reference Mardan:2020noh  shows how with a generalized polytropic equation of
state we can obtain new classes of polytropic models from the solution of
Einstein-Maxwell field equations for charged anisotropic fluid configurations.
Various values of polytropic index $n=1,2$ are used to develop the models.
Models that have been developed have helped regain masses and radii of eight
different stars. Analyses of the models indicate that they are well-behaved
and physically feasible.

In Chavanis:2021jwr  authors develop a cosmological model involving a scalar
field with a power-law potential associated with polytropic equations of
state. In order to simplify the equations of the problem, they consider a fast
oscillation regime called ”spintessence”. Various values of polytropic
constants and polytropic indexes are studied. The LambdaCDM model, the
Chaplygin gas model, and the Bose-Einstein condensate model are recovered.
Considering such an equation of state and examining models in F(R,w) are
worthy of investigation, and a significant impact can be expected.

##  5 Summary and Conclusion:

A modified gravitational model is proposed in this paper. The modification is
represented by the function $F(R,w)$, where $w$ is the equation of state.
Three different cases are discussed with functions
$F(R,w)=f(R)k(w)$,$F(R,w)=f(R)+k(w)$ and $F(R,w)=(R+\alpha Rw)$. The models
are studied under the Palatini formalism.
The first and second choices appear as $F(R)$ models with an effective
momentum-energy tensor. In contrast, the third choice takes the form of a
Hilbert-Einstein action with an effective momentum-energy tensor. The third
case shows no Jordan-Einstein frame ambiguity.
The slow-roll inflation model is examined to illustrate how $F(R,w)$ may
appear to mimic the standard non-minimal coupling to gravity.
We consider the Higgs inflationary potential as a case study and demonstrate
how the $F(R,w)$ model effectively addresses the challenges posed by Higgs
inflation.
Finally, we present three future directions for work and provide some initial
indications of how $F(R,w)$ models could impact cosmological results.
A constant-roll inflation model with $w=w(\phi)$ is discussed, and the Hubble
parameter is numerically impacted by this model.
In the $\Lambda$ CDM model, specifically under the $F(R,w)=F(R)(1+\alpha w)$
framework with $F(R)=R+\gamma R^{2}$, we examine the behavior of the scale
factor and Hubble parameter using the Starobinsky function. We compare the
results to a pure Starobinsky model to clarify the expected effect of
$F(R,w)$.
As a final point, we introduce the Polytropic equation of state as a field
that will be of interest to future works.

## References

  * (1) A. Einstein, Annalen Phys. 49, no.7, 769-822 (1916) doi:10.1002/andp.19163540702
  * (2) A. G. Riess et al. [Supernova Search Team], Astron. J. 116, 1009-1038 (1998) doi:10.1086/300499 [arXiv:astro-ph/9805201 [astro-ph]].
  * (3) S. Perlmutter et al. [Supernova Cosmology Project], Astrophys. J. 517, 565-586 (1999) doi:10.1086/307221 [arXiv:astro-ph/9812133 [astro-ph]].
  * (4) S. M. Carroll, V. Duvvuri, M. Trodden and M. S. Turner, Phys. Rev. D 70, 043528 (2004) doi:10.1103/PhysRevD.70.043528 [arXiv:astro-ph/0306438 [astro-ph]].
  * (5) W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007) doi:10.1103/PhysRevD.76.064004 [arXiv:0705.1158 [astro-ph]].
  * (6) C. Liu and H. Zhang, Eur. Phys. J. C 83, no.12, 1096 (2023) doi:10.1140/epjc/s10052-023-12261-9 [arXiv:2312.12702 [gr-qc]].
  * (7) S. D. Odintsov, D. Sáez-Chillón Gómez and G. S. Sharov, Phys. Dark Univ. 42, 101369 (2023) doi:10.1016/j.dark.2023.101369 [arXiv:2310.20302 [gr-qc]].
  * (8) S. D. Odintsov, V. K. Oikonomou and G. S. Sharov, Phys. Lett. B 843, 137988 (2023) doi:10.1016/j.physletb.2023.137988 [arXiv:2305.17513 [gr-qc]].
  * (9) M. AlHallak, N. Chamoun and M. S. Eldaher, JCAP 10, 001 (2022) doi:10.1088/1475-7516/2022/10/001 [arXiv:2202.01002 [astro-ph.CO]].
  * (10) M. AlHallak, A. AlRakik, N. Chamoun and M. S. El-Daher, Universe 8, no.2, 126 (2022) doi:10.3390/universe8020126 [arXiv:2111.05075 [astro-ph.CO]].
  * (11) M. AlHallak, Phys. Sci. Forum 7, no.1, 35 (2023) doi:10.3390/ECU2023-14048 [arXiv:2302.06524 [astro-ph.CO]].
  * (12) T. Harko, F. S. N. Lobo, S. Nojiri and S. D. Odintsov, Phys. Rev. D 84, 024020 (2011) doi:10.1103/PhysRevD.84.024020 [arXiv:1104.2669 [gr-qc]].
  * (13) J. Wu, G. Li, T. Harko and S. D. Liang, Eur. Phys. J. C 78, no.5, 430 (2018) doi:10.1140/epjc/s10052-018-5923-9 [arXiv:1805.07419 [gr-qc]].
  * (14) H. Azmat, Q. Muneer, M. Zubair, E. Gudekli, I. Ahmad and S. Waheed, Nucl. Phys. B 998, 116396 (2024) doi:10.1016/j.nuclphysb.2023.116396
  * (15) P. H. R. S. Moraes and P. K. Sahoo, Eur. Phys. J. C 79, no.8, 677 (2019) doi:10.1140/epjc/s10052-019-7206-5 [arXiv:1903.03421 [gr-qc]].
  * (16) P. H. R. S. Moraes, W. de Paula and R. A. C. Correa, Int. J. Mod. Phys. D 28, no.08, 1950098 (2019) doi:10.1142/S0218271819500986 [arXiv:1710.07680 [gr-qc]].
  * (17) P. H. R. S. Moraes, R. A. C. Correa and R. V. Lobato, JCAP 07, 029 (2017) doi:10.1088/1475-7516/2017/07/029 [arXiv:1701.01028 [gr-qc]].
  * (18) P. K. Sahoo, P. H. R. S. Moraes and P. Sahoo, Eur. Phys. J. C 78, no.1, 46 (2018) doi:10.1140/epjc/s10052-018-5538-1 [arXiv:1709.07774 [gr-qc]].
  * (19) P. Bhar, P. Rej and P. K. Sahoo, Int. J. Mod. Phys. D 31, no.03, 2250016 (2022) doi:10.1142/S021827182250016X [arXiv:2112.02623 [gr-qc]].
  * (20) M. Gamonal, Phys. Dark Univ. 31, 100768 (2021) doi:10.1016/j.dark.2020.100768 [arXiv:2010.03861 [gr-qc]].
  * (21) B. Deb and A. Deshamukhya, Int. J. Mod. Phys. A 37, no.18, 2250127 (2022) doi:10.1142/S0217751X22501275 [arXiv:2201.04378 [gr-qc]].
  * (22) S. Bhattacharjee, New Astron. 90, 101657 (2022) doi:10.1016/j.newast.2021.101657
  * (23) P. S. Singh and K. Priyokumar Singh, New Astron. 84, 101542 (2021) doi:10.1016/j.newast.2020.101542 [arXiv:2010.13784 [gr-qc]].
  * (24) G. Sun and Y. C. Huang, Int. J. Mod. Phys. D 25, no.03, 1650038 (2016) doi:10.1142/S0218271816500383 [arXiv:1510.01061 [gr-qc]].
  * (25) D. C. Maurya, J. Singh and L. K. Gaur, Int. J. Geom. Meth. Mod. Phys. 20, no.11, 2350192 (2023) doi:10.1142/S021988782350192X [arXiv:2212.05605 [gr-qc]].
  * (26) U. K. Sharma, M. Kumar and G. Varshney, Universe 8, no.12, 642 (2022) doi:10.3390/universe8120642
  * (27) A. Y. Shaikh, New Astron. 91, 101676 (2022) doi:10.1016/j.newast.2021.101676
  * (28) G. Varshney, U. K. Sharma, A. Pradhan and N. Kumar, Chin. J. Phys. 73, 56-73 (2021) doi:10.1016/j.cjph.2021.04.014
  * (29) A. Y. Shaikh and K. S. Wankhade, Found. Phys. 51, no.3, 58 (2021) doi:10.1007/s10701-021-00463-8 [arXiv:1912.08044 [physics.gen-ph]].
  * (30) T. Harko and F. S. N. Lobo, Eur. Phys. J. C 70, 373-379 (2010) doi:10.1140/epjc/s10052-010-1467-3 [arXiv:1008.4193 [gr-qc]].
  * (31) T. Naseer, M. Sharif, A. Fatima and S. Manzoor, Chin. J. Phys. 86, 350-360 (2023) doi:10.1016/j.cjph.2023.10.032
  * (32) N. S. Kavya, V. Venkatesha, G. Mustafa, P. K. Sahoo and S. V. D. Rashmi, Chin. J. Phys. 84, 1-11 (2023) doi:10.1016/j.cjph.2023.05.002 [arXiv:2305.01469 [gr-qc]].
  * (33) N. S. Kavya, V. Venkatesha, G. Mustafa and P. K. Sahoo, Annals Phys. 455, 169383 (2023) doi:10.1016/j.aop.2023.169383 [arXiv:2307.02498 [gr-qc]].
  * (34) D. C. Maurya, Phys. Dark Univ. 42, 101373 (2023) doi:10.1016/j.dark.2023.101373
  * (35) A. Pradhan, D. C. Maurya, G. K. Goswami and A. Beesham, Int. J. Geom. Meth. Mod. Phys. 20, no.06, 2350105 (2023) doi:10.1142/S0219887823501050 [arXiv:2209.14269 [gr-qc]].
  * (36) C. Dioguardi and A. Racioppi, [arXiv:2307.02963 [gr-qc]].
  * (37) H. Motohashi, A. A. Starobinsky and J. Yokoyama, JCAP 09, 018 (2015) doi:10.1088/1475-7516/2015/09/018 [arXiv:1411.5021 [astro-ph.CO]].
  * (38) S. A. Mardan, M. Rehman, I. Noureen and R. N. Jamil, Eur. Phys. J. C 80, no.2, 119 (2020) doi:10.1140/epjc/s10052-020-7647-x
  * (39) P. H. Chavanis, Phys. Rev. D 106, no.4, 043502 (2022) doi:10.1103/PhysRevD.106.043502 [arXiv:2111.01828 [gr-qc]].
  * (40) E. O. Pozdeeva, M. R. Gangopadhyay, M. Sami, A. V. Toporensky and S. Y. Vernov, Phys. Rev. D 102, no.4, 043525 (2020) doi:10.1103/PhysRevD.102.043525 [arXiv:2006.08027 [gr-qc]].
  * (41) B. Tahmasebzadeh and K. Karami, Nucl. Phys. B 918, 1-10 (2017) doi:10.1016/j.nuclphysb.2017.02.018 [arXiv:1608.06543 [gr-qc]].
  * (42) M. AlHallak, [arXiv:2311.06933 [astro-ph.CO]].
  * (43) A. Linde, M. Noorbala and A. Westphal, JCAP 03, 013 (2011) doi:10.1088/1475-7516/2011/03/013 [arXiv:1101.2652 [hep-th]].
  * (44) F. Bauer and D. A. Demir, Phys. Lett. B 665, 222-226 (2008) doi:10.1016/j.physletb.2008.06.014 [arXiv:0803.2664 [hep-ph]].
  * (45) M. AlHallak, A. Al Rakik, S. Bitar, N. Chamoun and M. S. Eldaher, Int. J. Mod. Phys. A 36, no.30, 2150226 (2021) doi:10.1142/S0217751X21502262 [arXiv:2105.00848 [astro-ph.CO]].
  * (46) M. AlHallak and N. Chamoun, JCAP 09, 006 (2016) doi:10.1088/1475-7516/2016/09/006 [arXiv:1604.03966 [hep-ph]].
  * (47) S. Karamitsos and A. Strumia, JHEP 05, 016 (2022) doi:10.1007/JHEP05(2022)016 [arXiv:2109.10367 [hep-th]].
  * (48) M. P. Hertzberg, JHEP 11, 023 (2010) doi:10.1007/JHEP11(2010)023 [arXiv:1002.2995 [hep-ph]].
  * (49) T. Qiu and K. C. Yang, Phys. Rev. D 83, 084022 (2011) doi:10.1103/PhysRevD.83.084022 [arXiv:1012.1697 [hep-th]].
  * (50) Y. Akrami et al. [Planck], Astron. Astrophys. 641, A10 (2020) doi:10.1051/0004-6361/201833887 [arXiv:1807.06211 [astro-ph.CO]].