# The Low Energy Spectrum of Trapped Bosons
in the Gross-Pitaevskii Regime

Christian Brennecke

Institute for Applied Mathematics, University of Bonn,
Endenicher Allee 60, 53115 Bonn, Germany

###### Abstract

Bogoliubov Theory [9] provides important predictions for the low energy
properties of the weakly interacting Bose gas. Recently, Bogoliubov’s
predictions could be justified rigorously in [7] for translation invariant
systems in the Gross-Pitaveskii regime, where $N$ bosons in
$\Lambda=[0;1]^{3}\subset\mathbb{R}^{3}$ interact through a potential whose
scattering length is of size $N^{-1}$. In this note, we review recent results
from [14], a joint work with B. Schlein and S. Schraven, which extends the
analysis of [7] to systems of bosons in $\mathbb{R}^{3}$ that are trapped by
an external potential.

##  1 Introduction and Main Results

The rigorous understanding of low energy properties of dilute Bose gases has
been an active field of research in the last decades after the first
experimental realization of Bose-Einstein condensates in trapped atomic gases
[2, 15]. As a mathematical model for such experimental setups, we consider
systems of $N$ trapped bosons moving in $\mathbb{R}^{3}$ in the Gross-
Pitaevskii regime. The energy of the system is described by

$H_{N}=\sum_{j=1}^{N}\left(-\Delta_{x_{j}}+V_{\text{ext}}(x_{j})\right)+\sum_{1\leq
i<j\leq N}N^{2}V(N(x_{i}-x_{j})),$ (1.1)

acting in a dense subspace of $L^{2}_{s}(\mathbb{R}^{3N})$, the subspace of
$L^{2}(\mathbb{R}^{3N})$-functions which are invariant under permutations of
the $N$ particle coordinates. We assume that the trapping potential
$V_{\text{ext}}\in L^{\infty}_{loc}(\mathbb{R}^{3})$ is such that
$V_{\text{ext}}(x)\to\infty$ for $|x|\to\infty$ and that $V\in
L^{3}(\mathbb{R}^{3})$ is pointwise non-negative, radial and of compact
support.

We are interested in the computation of the low energy spectrum of $H_{N}$ in
the limit $N\to\infty$. Notice that, due to the presence of a trapping
potential, the particles typically move in a region of volume
$\mathcal{O}(1)$. By a simple scaling argument, the system is equivalent to a
system of $N$ particles interacting through the unscaled interaction $V$ and
trapped in a region of volume $\mathcal{O}(N^{3})$. That is, the Gross-
Pitaevskii limit can be thought of as a joint thermodynamic and low density
limit (with the particle density of size $N^{-2}$). In particular, in view of
Bogoliubov theory [9], one may expect that the low energy spectrum depends on
the interaction to leading order only through its scattering length.

Recall that the scattering length $\mathfrak{a}_{0}$ of $V$ is defined through
the solution $f$ of the zero energy scattering equation

$\Big{(}-\Delta+\frac{1}{2}V\Big{)}f=0\hskip 14.22636pt\text{with}\hskip
14.22636pt\lim_{|x|\to\infty}f(x)=1.$ (1.2)

For $x\in\mathbb{R}^{3}$ outside the support of $V$, $f$ is explicitly given
by

$f(x)=1-\frac{\mathfrak{a}_{0}}{|x|}$

for some constant $\mathfrak{a}_{0}$, the _scattering length of $V$_. A simple
computation shows that

$\displaystyle\mathfrak{a}_{0}=\frac{1}{8\pi}\int_{\mathbb{R}^{3}}dx\,V(x)f(x)$

and that the solution of (1.2) with the rescaled potential $N^{2}V(N.)$ is
equal to $f(N.)$, i.e.

$\Big{(}-\Delta+\frac{1}{2}N^{2}V(N.)\Big{)}f(N.)=0.$

This implies that the scattering length of the interaction $N^{2}V(N.)$ is
equal to $\mathfrak{a}_{0}N^{-1}$. The size $N^{-1}$ of the scattering length
characterizes the Gross-Pitaeavskii regime.

The leading order contribution to the ground state energy $E_{N}$ of $H_{N}$
is well-known and has been derived in [29, 28, 31], showing that

$\lim_{N\to\infty}\frac{E_{N}}{N}=\inf_{\psi\in
L^{2}(\mathbb{R}^{3}):\|\psi\|_{2}=1}\mathcal{E}_{GP}(\psi).$

Here, $\mathcal{E}_{GP}$ denotes the Gross-Pitaevskii functional which is
defined by

$\mathcal{E}_{GP}(\psi)=\int_{\mathbb{R}^{3}}\left(|\nabla\psi(x)|^{2}+V_{\text{ext}}(x)|\psi(x)|^{2}+4\pi\mathfrak{a}_{0}|\psi(x)|^{4}\right)dx.$

By standard arguments, $\mathcal{E}_{GP}$ admits a unique normalized, strictly
positive minimizer which we denote from now on by $\varphi_{0}\in
L^{2}(\mathbb{R}^{3})$. It solves the Euler-Lagrange equation

$\displaystyle-\Delta\varphi_{0}+V_{\text{ext}}\varphi_{0}+8\pi\mathfrak{a}_{0}|\varphi_{0}|^{2}\varphi_{0}=\varepsilon_{GP}\varphi_{0},\hskip
14.22636pt\text{where}\hskip
14.22636pt\varepsilon_{GP}=\mathcal{E}_{GP}(\varphi_{0})+4\pi\mathfrak{a}_{0}\|\varphi_{0}\|_{4}^{4}.$

In order to determine the next to leading order contribution to the ground
state energy as well as the excitation energies of $H_{N}$, we need to use the
fact that approximate ground states exhibit complete Bose-Einstein
condensation (BEC). Indeed, it follows from [27, 28, 31] that every sequence
of normalized wave functions $\psi_{N}\in L^{2}_{s}(\mathbb{R}^{3N})$ with

$\lim_{N\to\infty}\frac{1}{N}\langle\psi_{N},H_{N}\psi_{N}\rangle=\mathcal{E}_{\text{GP}}(\varphi_{0})\,,$

that is, every sequence of approximate ground states, exhibits complete Bose-
Einstein condensation into the minimizer $\varphi_{0}$ of $\mathcal{E}_{GP}$.
Mathematically, this means that the one-particle reduced density
$\gamma_{N}^{(1)}=\operatorname{tr}_{2,\dots,N}|\psi_{N}\rangle\langle\psi_{N}|$
associated with $\psi_{N}$ satisfies

$\lim_{N\to\infty}\langle\varphi_{0},\gamma_{N}^{(1)}\varphi_{0}\rangle=1.$
(1.3)

Notice that this is equivalent to
$\lim_{N\to\infty}\gamma_{N}^{(1)}=|\varphi_{0}\rangle\langle\varphi_{0}|$ in
the trace class topology.

In the sequel, we make use of a quantitative version of (1.3) with optimal
rate of convergence that follows from [13]. More precisely, the main result of
[13] is that (under the assumptions $(1)$ and $(2)$ in (1.6) stated below)
there exists some $C>0$ such that

$H_{N}\geq
N\mathcal{E}_{\text{GP}}(\varphi_{0})+C^{-1}\sum_{j=1}^{N}\big{(}1-|\varphi_{0}\rangle\langle\varphi_{0}|\big{)}_{x_{j}}-C.$
(1.4)

In particular, if $\psi_{N}\in L^{2}_{s}(\mathbb{R}^{3N})$ is a normalized
sequence of approximate ground states such that
$\langle\psi_{N},H_{N}\psi_{N}\rangle\leq
N\mathcal{E}_{\text{GP}}(\varphi_{0})+\zeta$ for some $\zeta>0$, then
$\gamma_{N}^{(1)}$ satisfies

$1-\langle\varphi_{0},\gamma_{N}^{(1)}\varphi_{0}\rangle\leq
C\frac{(C+\zeta)}{N}.$

The lower bound (1.4) generalizes the main results of [5, 8], dealing with
translation invariant systems, to systems of trapped particles in
$\mathbb{R}^{3}$. For sufficiently small $\mathfrak{a}_{0}$, (1.4) has
previously been proved in [30] (adapting a complete-the-square argument of
[11] to systems of trapped bosons) and in [22] (for translation invariant
systems, simplifying the arguments of [5]). The bound (1.4) also follows from
[32] (with no additional smallness restriction on the size of
$\mathfrak{a}_{0}$) on which we comment briefly below. For further recent
results on BEC, including regimes beyond the Gross-Pitaevskii scaling, see [1,
18].

To state our main results, we denote by $H_{GP}$ the self-adjoint operator

$H_{\text{GP}}=-\Delta+V_{\text{ext}}+8\pi\mathfrak{a}_{0}\varphi_{0}^{2}-\varepsilon_{GP}$
(1.5)

and we make our assumptions on $V$ and $V_{ext}$ more precise:

$\begin{split}(1)\;&V\in L^{3}(\mathbb{R}^{3}),V(x)\geq 0\text{ for
}a.e.\;x\in\mathbb{R}^{3},V(x)=V(y)\text{ if }|x|=|y|,\\\ &V\text{ compactly
supported},\\\ (2)\;&V_{\text{ext}}\in
C^{1}\big{(}\mathbb{R}^{3};[0;\infty)\big{)},V_{\text{ext}}(x)\to\infty\text{
as }|x|\to\infty,\\\
&\exists\;C>0\;\forall\;x,y\in\mathbb{R}^{3}:V_{\text{ext}}(x+y)\leq
C(V_{\text{ext}}(x)+C)(V_{\text{ext}}(y)+C),\\\ &\nabla V_{\text{ext}}\text{
has at most exponential growth as }|x|\to\infty,\\\
(3)\;&(H_{\text{GP}}+1)^{-3/4-\varepsilon}e^{-\alpha|x|}\text{ is a Hilbert-
Schmidt operator, for all $\varepsilon>0$}\\\ &\text{and some
$\alpha>0$}.\end{split}$ (1.6)

######  Theorem 1.1.

_([14, Theorem 1.1])_. Let $V$, $V_{\text{ext}}$ and $H_{GP}$ satisfy (1.6).
Then, there exists $\rho>0$ so that, in the limit $N\to\infty$, the ground
state energy $E_{N}$ of $H_{N}$ is given by

$\begin{split}E_{N}&=N\mathcal{E}_{GP}(\varphi_{0})-4\pi\mathfrak{a}_{0}\|\varphi_{0}\|_{4}^{4}+E_{\text{Bog}}+\mathcal{O}\big{(}N^{-\rho}\big{)},\end{split}$
(1.7)

where $E_{\text{Bog}}$ is the finite constant (independent of $N$)

$\begin{split}E_{\text{Bog}}=\frac{1}{2}\lim_{\delta\to
0}&\left\\{\operatorname{tr}_{\perp\varphi_{0}}\Big{[}\Big{(}H_{\text{GP}}^{1/2}(H_{\text{GP}}+16\pi\mathfrak{a}_{0}1\\!\\!1_{\delta}\varphi_{0}^{2})H_{\text{GP}}^{1/2}\Big{)}^{1/2}-H_{\text{GP}}-8\pi\mathfrak{a}_{0}1\\!\\!1_{\delta}\varphi_{0}^{2}\Big{]}\right.\\\
&\hskip
156.49014pt\left.+\frac{(8\pi\mathfrak{a}_{0})^{2}}{2}\operatorname{tr}\Big{[}\varphi_{0}^{2}1\\!\\!1_{\delta}\frac{1}{-\Delta}1\\!\\!1_{\delta}\varphi_{0}^{2}\Big{]}\right\\}\end{split}$
(1.8)

with $1\\!\\!1_{\delta}$ denoting the approximation of the identity with
$1\\!\\!1_{\delta}(x;y)=(2\pi\delta)^{-3/2}e^{-(x-y)^{2}/2\delta^{2}}$ and
where $\operatorname{tr}_{\bot\varphi_{0}}$ denotes the trace over the
orthogonal complement $\\{\varphi_{0}\\}^{\bot}$.

Moreover, the spectrum of $H_{N}-E_{N}$ below a threshold $\zeta>0$ (assuming
$\zeta\leq CN^{\rho/5-\varepsilon}$ for some sufficiently small
$\varepsilon>0$) consists of finite sums of the form

$\sum_{i=1}^{\infty}n_{i}e_{i}+\mathcal{O}\big{(}N^{-\rho}(1+\zeta^{5}+N^{-5\rho}\zeta^{30})\big{)},$
(1.9)

where $(e_{j})_{j\in\mathbb{N}}$ denote the eigenvalues of the operator $E$
which is defined by

$E=\Big{(}H_{\text{GP}}^{1/2}\big{(}H_{\text{GP}}+16\pi\mathfrak{a}_{0}\varphi_{0}^{2}\big{)}H_{\text{GP}}^{1/2}\Big{)}^{1/2}$
(1.10)

and where $n_{i}\in\mathbb{N}$ with $n_{i}\neq 0$ for finitely many
$i\in\mathbb{N}$.

Remarks:

  1. 1)

The form (1.9) of the excitation spectrum of $H_{N}$ was conjectured in [21],
providing results comparable to Theorem 1.1 in the mean field limit (see also
[37, 16, 26] for previous and related results on the derivation of the
excitation spectrum for mean field systems of bosons). In the mean field
regime, more precise expansions in powers of $N^{-1}$ for the ground state
energy and low energy excitation energies have been obtained in [35, 10]. For
further background and results on mean field systems, see for instance the
recent survey article [36].

  2. 2)

Theorem 1.1 generalizes the main results of [7] (dealing with translation
invariant systems) to systems of trapped bosons. On a formal level, one can
recover the results of [7] by setting $V_{ext}=0$ in $\Lambda=[0;1]^{3}$,
$V_{ext}=\infty$ otherwise and by imposing periodic boundary conditions in
$\Lambda$. In this setting, the condensate wave function is described by
$\varphi=1_{|\Lambda}\in L^{2}(\Lambda)$ and $H_{GP}=-\Delta$. In particular,
$H_{GP}$ commutes with the multiplication operator
$16\pi\mathfrak{a}_{0}\varphi_{0}^{2}=16\pi\mathfrak{a}_{0}$ and $E$, defined
in (1.10), is the self-adjoint operator that multiplies in Fourier space by
$\sqrt{|p|^{4}+16\pi\mathfrak{a}_{0}|p|^{2}}$, for $p\in 2\pi\mathbb{Z}^{3}$.
For trapped systems as in [14], on the other hand, $H_{GP}$ and
$16\pi\mathfrak{a}_{0}\varphi_{0}^{2}$ do not commute anymore which makes the
analysis technically more involved, compared to [7].

  3. 3)

In (1.8), we cannot take directly the limit $\delta\to 0$, replacing
$1\\!\\!1_{\delta}$ with the identity, because the resulting operator is not a
trace class operator. However, using the assumptions (1.6), we can compute the
limit defining $E_{\text{Bog}}$ explicitly and find that

$\begin{split}E_{\emph{Bog}}&=\frac{\kappa^{2}}{4}(8\pi\mathfrak{a}_{0})^{2}\operatorname{tr}\left[\varphi_{0}^{2}(-\Delta+\kappa^{2})^{-1}(-\Delta)^{-1}\varphi_{0}^{2}\right]\\\
&\hskip
14.22636pt+\frac{1}{4}(8\pi\mathfrak{a}_{0})^{2}\left\|\varphi_{0}(-\Delta+\kappa^{2})^{-1}[\varphi_{0}^{2},-\Delta](-\Delta+\kappa^{2})^{-1/2}\right\|_{\text{HS}}^{2}\\\
&\hskip
14.22636pt+\frac{1}{4}(8\pi\mathfrak{a}_{0})^{2}\left\|\varphi_{0}(-\Delta+\kappa^{2})^{-1}\nabla\varphi_{0}\right\|_{\text{HS}}^{2}\\\
&\hskip
14.22636pt+\frac{(8\pi\mathfrak{a}_{0})^{2}}{4}\operatorname{tr}\,\Big{[}\varphi_{0}^{2}\Big{(}\frac{1}{-\Delta+\kappa^{2}}-\frac{1}{H_{\text{GP}}+\kappa^{2}}\Big{)}\varphi_{0}^{2}\Big{]}\\\
&\hskip
14.22636pt+\frac{(8\pi\mathfrak{a}_{0})^{2}}{4\kappa^{2}}\|\varphi_{0}^{3}\|^{2}+\frac{(8\pi\mathfrak{a}_{0})^{2}}{4}\|(H_{\text{GP}}+\kappa^{2})^{-1/2}Q\varphi_{0}^{3}\|^{2}\\\
&\hskip
14.22636pt-\frac{(8\pi\mathfrak{a}_{0})^{2}\kappa^{2}}{4}\operatorname{tr}_{\perp\varphi_{0}}\Big{[}\varphi_{0}^{2}Q\frac{1}{H_{\text{GP}}(H_{\text{GP}}+\kappa^{2})}Q\varphi_{0}^{2}\Big{]}\\\
&\hskip
14.22636pt-\frac{(8\pi\mathfrak{a}_{0})^{2}}{\pi}\int_{0}^{\infty}ds\,\sqrt{s}\,\operatorname{tr}_{\perp\varphi_{0}}\frac{H_{\text{GP}}^{1/2}}{s+H_{\text{GP}}^{2}}\Big{[}\Big{[}\varphi_{0}^{2},\frac{H_{\text{GP}}}{s+H_{\text{GP}}^{2}}\Big{]},\varphi_{0}^{2}\Big{]}\frac{H_{\text{GP}}^{1/2}}{s+H_{\text{GP}}^{2}}\\\
&\hskip
14.22636pt-\frac{4(8\pi\mathfrak{a}_{0})^{3}}{\pi}\int_{0}^{\infty}ds\,\sqrt{s}\operatorname{tr}_{\perp\varphi_{0}}\Big{[}\frac{1}{s+H_{\text{GP}}^{2}}H_{\text{GP}}^{1/2}\varphi_{0}^{2}H_{\text{GP}}^{1/2}\Big{]}^{3}\\\
&\hskip
170.71652pt\times\frac{1}{s+H_{\text{GP}}^{1/2}(H_{\text{GP}}+16\pi\mathfrak{a}_{0}\varphi_{0}^{2})H_{\text{GP}}^{1/2}}\end{split}$
(1.11)

for any $\kappa>0$ large enough. In particular, all contributions on the right
hand side of (1.11) are finite.

  4. 4)

The form of the excitation spectrum (1.9) has been derived independently in
[32], valid under slightly more general assumptions on $V$ and $V_{ext}$,
compared to (1.6). We refer to [32, 14] for detailed comparisons of the
assumptions and the proofs.

  5. 5)

The formulas (1.7) and (1.11) show that $E_{N}$ depends on the interaction $V$
only through its scattering length $\mathfrak{a}_{0}$, up to errors vanishing
in the limit $N\to\infty$. In particular, (1.7) is the analogue of the Lee-
Huang-Yang formula [23, 24, 25]

$E_{N,L}=4\pi\rho\mathfrak{a}_{0}\Big{(}1+\frac{128}{15\sqrt{\pi}}(\rho\mathfrak{a}_{0}^{3})^{1/2}+o\big{(}(\rho\mathfrak{a}_{0}^{3})^{1/2}\big{)}\Big{)}$
(1.12)

for the ground state energy $E_{N,L}$ of $N$ bosons trapped in
$\Lambda_{L}=[0;L]^{3}$, in the thermodynamic limit where $N,L\to\infty$ with
the density $\rho=NL^{-3}$ fixed. (1.12) has been rigorously established in
[38, 3] (upper bound) and [19, 20] (lower bound).

In the following sections, we outline the main steps that lead to the proof of
Theorem 1.1. Our proof is based on a rigorous implementation of Bogoliubov’s
method [9],previously established in [7] for translation invariant systems in
the Gross-Pitaevskii regime. For the details of the implementation, we refer
to [14].

##  2 Excitation Hamiltonians

For our analysis, it is convenient to switch to a Fock space setting,
factoring out the Bose-Einstein condensate. This allows us to focus on the
analysis of the orthogonal excitations around the condensate. We follow here
[26] and consider the unitary map
$U_{N}:L^{2}_{s}(\mathbb{R}^{3N})\to\mathcal{F}_{\perp\varphi_{0}}^{\leq N}$
defined by $U_{N}\psi_{N}=(\alpha_{0},\alpha_{1},\dots,\alpha_{N})$ if

$\psi_{N}=\alpha_{0}\varphi_{0}^{\otimes
N}+\alpha_{1}\otimes_{s}\varphi_{0}^{\otimes(N-1)}+\ldots+\alpha_{N}$

with $\alpha_{j}\in L^{2}_{\perp\varphi_{0}}(\mathbb{R}^{3})^{\otimes j}$, for
$j=0,1,\dots,N$. Here, $\otimes_{s}$ denotes the symmetric tensor product and
$\mathcal{F}_{\perp\varphi_{0}}^{\leq N}$ denotes the truncated excitation
Fock space

$\mathcal{F}_{\perp\varphi_{0}}^{\leq
N}=\bigoplus_{n=0}^{N}L^{2}_{\perp\varphi_{0}}(\mathbb{R}^{3})^{\otimes_{s}n}$

that is built over the orthogonal complement
$L^{2}_{\perp\varphi_{0}}(\mathbb{R}^{3})=\\{\varphi_{0}\\}^{\bot}$. The
property (1.3) of BEC translates in the Fock space setting to the statement
that

$\lim_{N\to\infty}\Big{(}1-\langle\varphi_{0},\gamma_{N}^{(1)}\varphi_{0}\rangle\Big{)}=\lim_{N\to\infty}N^{-1}\langle
U_{N}\psi_{N},\mathcal{N}U_{N}\psi_{N}\rangle=0,$

where $\mathcal{N}$ denotes the number of particles operator in
$\mathcal{F}_{\perp\varphi_{0}}^{\leq N}$, defined by
$(\mathcal{N}\xi)^{(n)}=n\xi^{(n)}$ for every
$\xi=(\xi^{(1)},\dots,\xi^{(N)})\in\mathcal{F}_{\perp\varphi_{0}}^{\leq N}$.
In other words, condensation into $\varphi_{0}$ means that the average number
of excitations is small compared to $N$, the total number of particles.

In $\mathcal{F}_{\perp\varphi_{0}}^{\leq N}$, we then consider
$\mathcal{L}_{N}=U_{N}H_{N}U_{N}^{*}$ which reads in second quantized form

$\begin{split}\mathcal{L}_{N}=\;&\big{\langle}\varphi_{0},\big{(}-\Delta+V_{\text{ext}}+\frac{1}{2}\big{(}N^{3}V(N\cdot)*|\varphi_{0}|^{2}\big{)}\big{)}\varphi_{0}\big{\rangle}(N-\mathcal{N})\\\
&\qquad-\frac{1}{2}\big{\langle}\varphi_{0},\big{(}N^{3}V(N\cdot)*|\varphi_{0}|^{2}\big{)}\varphi_{0}\big{\rangle}(\mathcal{N}+1)(1-\mathcal{N}/N)\\\
&+\bigg{(}\sqrt{N}b\left(\left(N^{3}V(N\cdot)*|\varphi_{0}|^{2}-8\pi\mathfrak{a}_{0}|\varphi_{0}|^{2}\right)\varphi_{0}\right)\\\
&\hskip
14.22636pt-\frac{\mathcal{N}+1}{\sqrt{N}}b\left(\left(N^{3}V(N\cdot)*|\varphi_{0}|^{2}\right)\varphi_{0}\right)+\text{h.c.}\bigg{)}\\\
&+\int dx\;a_{x}^{*}(-\Delta_{x})a_{x}+\int
dx\;V_{\text{ext}}(x)a_{x}^{*}a_{x}\\\ &\hskip 14.22636pt+\int
dxdy\;N^{3}V(N(x-y))\varphi_{0}^{2}(y)\Big{(}b_{x}^{*}b_{x}-\frac{1}{N}a_{x}^{*}a_{x}\Big{)}\\\
&\hskip 14.22636pt+\int
dxdy\;N^{3}V(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\Big{(}b_{x}^{*}b_{y}-\frac{1}{N}a_{x}^{*}a_{y}\Big{)}\\\
&\hskip 14.22636pt+\frac{1}{2}\int
dxdy\;N^{3}V(N(x-y))\varphi_{0}(y)\varphi_{0}(x)\Big{(}b_{x}^{*}b_{y}^{*}+b_{x}b_{y}\Big{)}\\\
&+\int
dxdy\;N^{5/2}V(N(x-y))\varphi_{0}(y)\big{(}b_{x}^{*}a^{*}_{y}a_{x}+a^{*}_{x}a_{y}b_{x}\big{)}\\\
&+\frac{1}{2}\int
dxdy\;N^{2}V(N(x-y))a_{x}^{*}a_{y}^{*}a_{y}a_{x}.\end{split}$ (2.1)

Here, $a_{x},a_{y}^{*}$, for $x,y\in\mathbb{R}^{3}$, denote the usual creation
and annihilation operators, satisfying the canonical commutation relations
$[a_{x},a_{y}^{*}]=\delta(x-y)$ and $[a_{x},a_{y}]=[a_{x}^{*},a_{y}^{*}]=0$,
and the operators $b_{x},b^{*}_{y}$ denote modified creation and annihilation
operators, defined by

$b_{x}=\sqrt{1-\mathcal{N}/N}\,a_{x},\qquad
b^{*}_{x}=a^{*}_{x}\,\sqrt{1-\mathcal{N}/N}.$

The operators $b_{x},b^{*}_{y}$ preserve the number of particles truncation in
$\mathcal{F}_{\perp\varphi_{0}}^{\leq N}$ and we have

$[b_{x},b_{y}^{*}]=(1-N^{-1}\mathcal{N})\delta(x-y)-N^{-1}a_{y}^{*}a_{x},\hskip
14.22636pt[b_{x},b_{y}]=[b_{x}^{*},b_{y}^{*}]=0.$

Hence, on states with a low number of excitations $\mathcal{N}\ll N$, the
modified field operators satisfy the usual commutation relations up to errors
that vanish in the limit $N\to\infty$.

Extracting the low energy spectrum of $H_{N}$ by analyzing $\mathcal{L}_{N}$
directly seems difficult, because we face a number of serious problems. First,
the constant contribution in (2.1),

$N\big{\langle}\varphi_{0},\big{(}-\Delta+V_{\text{ext}}+\frac{1}{2}\big{(}N^{3}V(N\cdot)*|\varphi_{0}|^{2}\big{)}\big{)}\varphi_{0}\big{\rangle}=N\mathcal{E}_{GP}(\varphi_{0})+\mathcal{O}(N),$

is off by a quantity of order $\mathcal{O}(N)$ from the correct leading order
energy $N\mathcal{E}_{GP}(\varphi_{0})$. Second, the terms that are linear in
the creation and annihilation operators,

$\sqrt{N}b\left(\left(N^{3}V(N\cdot)*|\varphi_{0}|^{2}-8\pi\mathfrak{a}_{0}|\varphi_{0}|^{2}\right)\varphi_{0}\right)+\text{h.c.},$

are a priori of size $\mathcal{O}(N^{1/2})$ and, similarly, the simple form
bounds

$\begin{split}&\pm\frac{1}{2}\int
dxdy\;N^{3}V(N(x-y))\varphi_{0}(y)\varphi_{0}(x)\big{(}b_{x}^{*}b_{y}^{*}+b_{x}b_{y}\big{)}\\\
&\hskip 5.69046pt\leq\frac{1}{2}\int
dxdy\;N^{2}V(N(x-y))a_{x}^{*}a_{y}^{*}a_{y}a_{x}+\frac{N}{2}\|V\|_{1}\|\varphi_{0}\|_{\infty}^{2}\|\varphi_{0}\|_{2}^{2}\end{split}$

indicate that the off-diagonal pairing terms contribute to order
$\mathcal{O}(N)$. The reason for these difficulties is that through the map
$U_{N}$, we expand $H_{N}$ around the energy of the state
$\varphi_{0}^{\otimes N}$, which neglects the correlations among the
particles. The correlations, however, are crucial even for extracting the
leading order term $N\mathcal{E}_{GP}(\varphi_{0})$.

In order to resolve these difficulties, the first step of our analysis
consists of renormalizing $\mathcal{L}_{N}$ through conjugation by a suitable
unitary map that extracts the missing correlation energies of order
$\mathcal{O}(N)$. Here, we make use of ideas that were first implemented in
the dynamical context [4, 12]. To motivate the approach heuristically, let’s
consider the right hand side in (2.1) for the moment as an operator in
$\mathcal{F}^{\leq N}=\bigoplus_{n=0}^{N}L^{2}_{s}(\mathbb{R}^{3n})$, ignoring
the orthogonality constraints in $\mathcal{F}_{\bot\varphi_{0}}^{\leq N}$.
Setting

$\mathcal{H}_{N}=\int dx\;a_{x}^{*}(-\Delta_{x})a_{x}+\int
dx\;V_{\text{ext}}(x)a_{x}^{*}a_{x}+\frac{1}{2}\int
dxdy\;N^{2}V(N(x-y))a_{x}^{*}a_{y}^{*}a_{y}a_{x}$

and $\widetilde{\eta}(x;y)=(-N)(1-f)(N(x-y))\varphi_{0}(x)\varphi_{0}(y)$,
recalling that $f$ denotes the solution of the zero energy scattering equation
(1.2), as well as

$B(\widetilde{\eta})=\frac{1}{2}\int
dxdy\;\widetilde{\eta}(x;y)b^{*}_{x}b^{*}_{y}-\text{h.c.},$

we are going to conjugate the terms on the right hand side in (2.1) by the
unitary operator exponential $e^{B(\widetilde{\eta})}:\mathcal{F}^{\leq
N}\to\mathcal{F}^{\leq N}$, using the second order Taylor expansion

$e^{-B(\widetilde{\eta})}O\,e^{B(\widetilde{\eta})}\approx
O+[O,B(\widetilde{\eta})]+\frac{1}{2}\big{[}[O,B(\widetilde{\eta})],B(\widetilde{\eta})\big{]}$

for observables $O$ in $\mathcal{F}^{\leq N}$. The main observation that
follows from [4, 12] is that conjugation by
$e^{-B(\widetilde{\eta})}(\cdot)e^{B(\widetilde{\eta})}$ leads to the
renormalizations

$\begin{split}&e^{-B(\widetilde{\eta})}\Big{(}\sqrt{N}b\left(\left(N^{3}V(N\cdot)*|\varphi_{0}|^{2}-8\pi\mathfrak{a}_{0}|\varphi_{0}|^{2}\right)\varphi_{0}\right)+\text{h.c.}\Big{)}e^{B(\widetilde{\eta})}\\\
&\hskip 14.22636pt+\int
dxdy\;N^{5/2}V(N(x-y))\varphi_{0}(y)e^{-B(\widetilde{\eta})}\big{(}b_{x}^{*}a^{*}_{y}a_{x}+a^{*}_{x}a_{y}b_{x}\big{)}e^{B(\widetilde{\eta})}\\\
&\approx\int
dxdy\;N^{5/2}V(N(x-y))\varphi_{0}(y)\big{(}b_{x}^{*}a^{*}_{y}a_{x}+a^{*}_{x}a_{y}b_{x}\big{)}\end{split}$
(2.2)

and

$\begin{split}&N\big{\langle}\varphi_{0},\big{(}-\Delta+V_{\text{ext}}+\frac{1}{2}\big{(}N^{3}V(N\cdot)*|\varphi_{0}|^{2}\big{)}\big{)}\varphi_{0}\big{\rangle}+e^{-B(\widetilde{\eta})}\mathcal{H}_{N}e^{B(\widetilde{\eta})}\\\
&\hskip 14.22636pt+\frac{1}{2}\int
dxdy\;N^{3}V(N(x-y))\varphi_{0}(y)\varphi_{0}(x)e^{-B(\widetilde{\eta})}\Big{(}b_{x}^{*}b_{y}^{*}+b_{x}b_{y}\Big{)}e^{B(\widetilde{\eta})}\\\
&\approx N\mathcal{E}_{GP}(\varphi_{0})+\mathcal{H}_{N},\end{split}$ (2.3)

while the remaining parts on the r.h.s. in (2.1) are essentially left
unchanged, up to errors of order $\mathcal{O}(1)$ which are well under
control. In other words, conjugation by
$e^{-B(\widetilde{\eta})}(\cdot)e^{B(\widetilde{\eta})}$ leads to the
cancellation of the large linear and quadratic pairing terms in (2.1), and it
renormalizes the constant $\mathcal{O}(N)$ contribution to the correct energy
$N\mathcal{E}_{GP}(\varphi_{0})$.

To make the above heuristics rigorous, we need to make a few technical
adjustments. First of all, instead of conjugating with the exponential of
$B(\widetilde{\eta})$, we choose a related operator $B=B(\eta)$ with kernel
$\eta\in(Q\otimes Q)L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})$, for
$Q=1-|\varphi_{0}\rangle\langle\varphi_{0}|$. This ensures that $e^{B}$ is a
unitary map from $\mathcal{F}_{\bot\varphi_{0}}^{\leq N}$ to itself. Moreover,
if $\|\eta\|_{2}$ is sufficiently small, we know from [12, 7, 8, 13] that one
has the exact expansion

$e^{-B}b(g)e^{B}=b(\cosh_{\eta}(g))+b^{*}(\sinh_{\eta}(\overline{g}))+d_{\eta}(g),$
(2.4)

for every $g\in L^{2}_{\perp\varphi_{0}}(\mathbb{R}^{3})$, where
$b(g)=\int_{\mathbb{R}^{3}}dx\,\overline{g}(x)b_{x}$ and where

$\begin{split}\cosh_{\eta}=\sum_{n=0}^{\infty}\frac{\eta^{(2n)}}{(2n)!},\hskip
14.22636pt\sinh_{\eta}=\sum_{n=0}^{\infty}\frac{\eta^{(2n+1)}}{(2n+1)!}.\end{split}$
(2.5)

Moreover, the remainder $d_{\eta}(g)$ in (2.4) is small on states with a small
number $\mathcal{N}\ll N$ of excitations, in the sense that
$\|d_{\eta}(g)\xi\|\leq CN^{-1}\|g\|\|(\mathcal{N}+1)^{3/2}\xi\|$ for all
$\xi\in\mathcal{F}_{\bot\varphi_{0}}^{\leq N}$.

To choose a kernel $\eta\in(Q\otimes
Q)L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})$ with small
$L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})$ norm that still leads to the
cancellations (2.2) and (2.3) (up to lower order corrections), one could
directly modify $\widetilde{\eta}$ from above, for example by imposing a
suitable cutoff in position space (such an approach is carried out in [32]).
We use instead previously established tools from the works [12, 7, 8, 13] and
consider the Neumann ground state $f_{\ell}$ solving

$\Big{(}-\Delta+\frac{1}{2}V\Big{)}f_{\ell}=\lambda_{\ell}f_{\ell}$

in the ball $B_{N\ell}$ of radius $N\ell$, fixing $0<\ell<1$, and normalized
such that $f_{\ell}(x)=1$ if $|x|=N\ell$ (notice that $f_{\ell}$ is radially
symmetric). Note that, by scaling, $f_{\ell}(N.)$ solves

$\Big{(}-\Delta+\frac{1}{2}N^{2}V(N.)\Big{)}f_{\ell}(N.)=N^{2}\lambda_{\ell}f_{\ell}(N.)$

in the ball $B_{\ell}$ of radius $\ell$. We then set $w_{\ell}=1-f_{\ell}$
which satisfies the pointwise bounds

$0\leq w_{\ell}(x)\leq\frac{C}{|x|+1}\quad\text{ and }\quad|\nabla
w_{\ell}(x)|\leq\frac{C}{|x|^{2}+1}.$

Defining $\eta=(Q\otimes Q)k\in(Q\otimes
Q)L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})$ for
$k(x;y)=-Nw_{\ell}(N(x-y))\varphi_{0}(x)\varphi_{0}(y)$, it is then simple to
check that $\|\eta\|_{2}\leq C\ell^{1/2}$ so that for $\ell>0$ sufficiently
small (independently of $N$), we can make use of the expansions (2.4).

Using the unitary map $e^{B}$ with $B=B(\eta)$ defined above, we can now
analyze the excitation Hamiltonian
$\mathcal{G}_{N}=e^{-B}\mathcal{L}_{N}e^{B}$. Before stating its main
properties, let us define

$\begin{split}\mathcal{K}\\!=\\!\int
dx\;a^{*}_{x}(-\Delta_{x})a_{x},\,\mathcal{V}_{\text{ext}}\\!=\\!\int
dx\,V_{\text{ext}}(x)a^{*}_{x}a_{x},\,\mathcal{V}_{N}\\!=\\!\frac{1}{2}\int
dxdy\;N^{2}V(N(x-y))a^{*}_{x}a^{*}_{y}a_{y}a_{x}\end{split}$

so that
$\mathcal{H}_{N}=\mathcal{K}+\mathcal{V}_{\text{ext}}+\mathcal{V}_{N}$. The
following proposition was proved in [14, Prop 2.5].

######  Proposition 2.1.

Assume (1.6), let $\mathcal{G}_{N}=e^{-B}\mathcal{L}_{N}e^{B}$ and let $E_{N}$
denote the ground state energy of $H_{N}$, defined in (1.1). Then, the
following holds true:

  * a)

We have that $|E_{N}-N\mathcal{E}_{GP}(\varphi_{0})|\leq C$ and

$\mathcal{G}_{N}-E_{N}=\mathcal{H}_{N}+\Delta_{N},$

where the error term $\Delta_{N}$ is such that for every $\delta>0$, there
exists $C>0$ with

$\pm\Delta_{N}\leq\delta\mathcal{H}_{N}+C(\mathcal{N}+1).$

Furthermore, for every $k\in\mathbb{N}$ there exists a $C>0$ such that

$\pm\text{ad}^{\,(k)}_{\,i\mathcal{N}}(\mathcal{G}_{N})=\pm\text{ad}^{\,(k)}_{\,i\mathcal{N}}(\Delta_{N})=\pm\big{[}i\mathcal{N},\dots\big{[}i\mathcal{N},\Delta_{N}\big{]}\dots\big{]}\leq
C(\mathcal{H}_{N}+1).$ (2.6)

  * b)

Let $\sigma=\sinh_{\eta}$ and $\gamma=\cosh_{\eta}$ be defined as in (2.5) and
let $\kappa_{\mathcal{G}_{N}}$ denote the constant

$\begin{split}\kappa_{\mathcal{G}_{N}}&=N\big{\langle}\varphi_{0},(-\Delta+V_{\text{ext}}+\widehat{V}(0)\varphi_{0}^{2}/2)\varphi_{0}\big{\rangle}-4\pi\mathfrak{a}_{0}\|\varphi_{0}\|_{4}^{4}\\\
&\hskip
14.22636pt+\operatorname{tr}\big{(}\sigma(-\Delta+V_{\text{ext}}-\varepsilon_{GP})\sigma\big{)}\\\
&\hskip
14.22636pt+\operatorname{tr}\big{(}\gamma\big{[}N^{3}V(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{]}\sigma)\\\
&\hskip
14.22636pt+\operatorname{tr}\big{(}\sigma\big{[}N^{3}V(N.)\ast\varphi_{0}^{2}+N^{3}V(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{]}\sigma\big{)}\\\
&\hskip 14.22636pt+\frac{1}{2}\int
dxdy\,N^{2}V(N(x-y))|\langle\sigma_{x},\gamma_{y}\rangle|^{2}.\end{split}$
(2.7)

Here, $(N^{3}V(N.)\ast\varphi_{0}^{2})$ acts as multiplication operator and we
identify kernels like $N^{3}V(N(x-y))\varphi_{0}(x)\varphi_{0}(y)$ with their
associated Hilbert Schmidt operators. Let

$\begin{split}\Phi&=\gamma\big{(}-\Delta+V_{\text{ext}}-\varepsilon_{GP}\big{)}\gamma+\sigma\big{(}-\Delta+V_{\text{ext}}-\varepsilon_{GP}\big{)}\sigma\\\
&\hskip
11.38092pt+\gamma\big{(}N^{3}V(N.)\ast\varphi_{0}^{2}+N^{3}V(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\gamma\\\
&\hskip
11.38092pt+\sigma\big{(}N^{3}V(N.)\ast\varphi_{0}^{2}+N^{3}V(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\sigma\\\
&\hskip
11.38092pt+\Big{(}\gamma\big{(}N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\sigma+\emph{h.c.}\Big{)},\end{split}$

and

$\begin{split}\Gamma&=\gamma\big{(}N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\gamma\\\
&\hskip
14.22636pt+\sigma\\!\big{(}N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\\!\sigma\\\
&\hskip
14.22636pt+\big{[}\sigma\big{(}-\Delta+V_{\text{ext}}-\varepsilon_{GP}\big{)}\gamma+\emph{h.c.}\big{]}\\\
&\hskip
14.22636pt+\big{[}\sigma\big{(}N^{3}V(N.)\ast\varphi_{0}^{2}+N^{3}V(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\gamma+\emph{h.c.}\big{]}.\\\
\end{split}$ (2.8)

With $\Phi$ and $\Gamma$, we define the quadratic Fock space Hamiltonian

$\begin{split}\mathcal{Q}_{\mathcal{G}_{N}}=\int
dxdy\,\Phi(x;y)b^{*}_{x}b_{y}+\frac{1}{2}\int
dxdy\,\Gamma(x;y)\big{(}b^{*}_{x}b^{*}_{y}+b_{x}b_{y}\big{)},\end{split}$
(2.9)

and we denote by $\mathcal{C}_{\mathcal{G}_{N}}$ the cubic Fock space operator

$\mathcal{C}_{\mathcal{G}_{N}}=\int
dxdy\,N^{5/2}V(N(x-y))\varphi_{0}(y)b^{*}_{x}b^{*}_{y}\big{(}b(\gamma_{x})+b^{*}(\sigma_{x})\big{)}+\emph{h.c.},$
(2.10)

where $\gamma_{x}(\cdot)=\gamma(\cdot;x)$ and
$\sigma_{x}(\cdot)=\sigma(\cdot;x)$ for $x\in\mathbb{R}^{3}$. Then,
$\mathcal{G}_{N}$ is equal to

$\begin{split}\mathcal{G}_{N}&=\kappa_{\mathcal{G}_{N}}+\mathcal{Q}_{\mathcal{G}_{N}}+\mathcal{C}_{\mathcal{G}_{N}}+\mathcal{V}_{N}+\mathcal{E}_{\mathcal{G}_{N}}\end{split}$
(2.11)

for some self-adjoint operator $\mathcal{E}_{\mathcal{G}_{N}}$ which is
bounded in $\mathcal{F}_{\bot\varphi_{0}}^{\leq N}$ by

$\pm\mathcal{E}_{\mathcal{G}_{N}}\leq
CN^{-1/2}(\mathcal{H}_{N}+\mathcal{N}^{2}+1)(\mathcal{N}+1).$ (2.12)

Proposition (2.1) implements the heuristic renormalizations in (2.2) and (2.3)
rigorously, noticing that
$\kappa_{\mathcal{G}_{N}}=N\mathcal{E}_{GP}(\varphi_{0})+\mathcal{O}(1)$ and
$\|\Gamma\|_{2}\leq C$, uniformly in $N\in\mathbb{N}$, for
$\kappa_{\mathcal{G}_{N}}$ and $\Gamma$ defined in (2.7) and (2.8),
respectively (see [14] for the details). Moreover, the error
$\mathcal{E}_{\mathcal{G}_{N}}$ in (2.11) satisfying (2.12), is negligible on
low energy eigenstates. Indeed, this is a direct consequence of the following
proposition which was proved in [14, Theorem 2.6] by combining the main result
(1.4) of [13] with the commutator estimates (2.6).

######  Proposition 2.2.

Assume (1.6) and let $E_{N}$ denote the ground state energy of $H_{N}$,
defined in (1.1). For $\zeta>0$, let $\psi_{N}\in L^{2}_{s}(\mathbb{R}^{3N})$
with $\|\psi_{N}\|=1$ be such that

$\psi_{N}={\bf 1}_{(-\infty;E_{N}+\zeta]}(H_{N})\psi_{N}$

and let $\xi_{N}=e^{-B}U_{N}\psi_{N}\in\mathcal{F}_{\bot\varphi_{0}}^{\leq N}$
denote the renormalized excitation vector related to $\psi_{N}$. Then, for
every $j\in\mathbb{N}$ there is a constant $C>0$ such that

$\displaystyle\langle\xi_{N},(\mathcal{H}_{N}+1)(\mathcal{N}+1)^{j}\xi_{N}\rangle\leq
C(1+\zeta^{j+1}).$

With Bogoliubov’s arguments [9] in mind, we may be tempted to neglect the
cubic and quartic contributions $\mathcal{C}_{\mathcal{G}_{N}}$ and
$\mathcal{V}_{N}$ on the r.h.s. in (2.10) and diagonalize the remaining
quadratic operator $\mathcal{Q}_{\mathcal{G}_{N}}$ in (2.9) to determine the
spectrum of $\mathcal{G}_{N}$ up to order $\mathcal{O}(1)$. Proceeding this
way does, however, not yield the right spectrum, because the terms
$\mathcal{C}_{\mathcal{G}_{N}}$ and $\mathcal{V}_{N}$ still contain important
energy contributions of order $\mathcal{O}(1)$. This has been a key
observation for the analysis in [7] (see also [38] and [33, 34] related to
this point) so that, in order to conclude Theorem 1.1, we first need to
extract the missing $\mathcal{O}(1)$ energies.

To this end, we proceed similarly as in the first step, now having the goal to
renormalize the cubic contribution $\mathcal{C}_{\mathcal{G}_{N}}$ in (2.10).
To understand the main idea on a heuristic level, consider for simplicity
first the cubic operator

$\int dxdy\,N^{5/2}V(N(x-y))\varphi_{0}(y)b^{*}_{x}b^{*}_{y}b_{x}+\text{h.c.}$
(2.13)

that is obtained from $\mathcal{C}_{\mathcal{G}_{N}}$ by replacing
$\gamma_{x}\approx\delta_{x}$ and $\sigma_{x}\approx 0$ (to leading order in
$\eta$). By analogy to the quadratic renormalization through $e^{B}$, a
straightforward computation involving the scattering equation (1.2) implies
that

$\begin{split}&\int
dxdy\,(-N(1-f)(N(x-y))\varphi_{0}(y)[\mathcal{K}_{N}+\mathcal{V}_{N},b^{*}_{x}b^{*}_{y}b_{x}-\text{h.c.}]\\\
&\hskip 14.22636pt\approx-\int
dxdy\,N^{5/2}V(N(x-y))\varphi_{0}(y)b^{*}_{x}b^{*}_{y}b_{x}+\text{h.c.},\end{split}$

up to errors of order $\mathcal{O}(1)$. What this indicates is that we can
cancel the cubic term (2.13) by conjugating $\mathcal{G}_{N}$ through a
unitary operator exponential that is now cubic in the (modified) creation and
annihilation operators, having the form

$\widetilde{A}=\frac{1}{\sqrt{N}}\int
dxdy\,(-N(1-f)(N(x-y))\varphi_{0}(y){b}_{x}^{*}{b}_{y}^{*}b_{x}-\text{h.c.}$
(2.14)

From the analysis of the quadratic renormalization and the fact that our goal
is to cancel the full cubic term $\mathcal{C}_{\mathcal{G}_{N}}$ in
$\mathcal{G}_{N}$ (and not only the term in (2.13)), it is clear that we need
to make a few technical adjustments to $\widetilde{A}$ in (2.14). First of
all, to stay consistent with the first renormalization step, instead of
working with the correlation factor $(-N(1-f))(N.)$, we use the function
$(-Nw_{\ell})(N.)$ on which we impose additionally an $N$-dependent low
momentum cutoff. More precisely, for some small $\varepsilon>0$, we set
$\chi_{H}(p)=\chi(|p|>N^{\varepsilon})$, denote by $\check{\chi}_{H}$ its
inverse Fourier transform and define the kernel

$\widetilde{k}_{H}(x;y)=\left(-Nw_{\ell}(N.)*\check{\chi}_{H}\right)(x-y)\varphi_{0}(y).$
(2.15)

Through the cutoff $\check{\chi}_{H}$, it turns out that it is enough to
compute only a few commutators in the expansion of the cubic conjugation,
because most of the error terms become small (in orders of $N$). Second, to
cancel not only the cubic term in (2.13), but the full term
$\mathcal{C}_{\mathcal{G}_{N}}$ in (2.10), we need to replace the
$b_{x}$-field in (2.14) by $b(\gamma_{x})+b^{*}(\sigma_{x})$. Also here, it
turns out to simplify the analysis if we impose high momentum cutoffs on the
kernels $\gamma_{x}$ and $\sigma_{x}$. To this end, we choose
$0<\tau<\varepsilon$, we set $g_{L}(p)=e^{-p^{2}/N^{2\tau}}$ and we denote by
$\check{g}_{L}$ its inverse Fourier transform. With the notation

$\sigma_{L}=\sigma*_{2}\check{g}_{L},\qquad\gamma_{L}=\gamma*_{2}\check{g}_{L},$

where $*_{2}$ denotes convolution in the second variable, we then define the
cubic operator

$A=\frac{1}{\sqrt{N}}\int
dxdy\,\widetilde{k}_{H}(x;y)\widetilde{b}_{x}^{*}\widetilde{b}_{y}^{*}\left[\widetilde{b}(\gamma_{L,x})+\widetilde{b}^{*}(\sigma_{L,x})\right]-\text{h.c.}$

where $\widetilde{b}_{x}=b(Q_{x})=\int dz\,Q(x;z)b_{z}$, with
$Q=1-|\varphi_{0}\rangle\langle\varphi_{0}|$. Notice that the use of the
operators $\widetilde{b}_{x},\widetilde{b}^{*}_{y}$ ensures that $A$ is a
unitary map from $\mathcal{F}_{\perp\varphi_{0}}^{\leq N}$ to itself.

The following proposition summarizes the main properties of the renormalized
excitation Hamiltonian $\mathcal{J}_{N}=e^{-A}\mathcal{G}_{N}e^{A}$. For its
proof, see [14, Prop. 2.8].

######  Proposition 2.3.

Assume (1.6) and let $0<6\tau\leq\varepsilon\leq\frac{1}{2}$. Then

$\mathcal{J}_{N}=\kappa_{\mathcal{J}_{N}}+\mathcal{Q}_{\mathcal{J}_{N}}+\mathcal{V}_{N}+\mathcal{E}_{\mathcal{J}_{N}},$
(2.16)

where

$\begin{split}\kappa_{\mathcal{J}_{N}}&=\kappa_{\mathcal{G}_{N}}-\emph{tr}(\sigma
N^{3}(Vw_{\ell})(N.)*\varphi_{0}^{2}\sigma)-\emph{tr}(\sigma
N^{3}(Vw_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\sigma)\end{split}$ (2.17)

with $\kappa_{\mathcal{G}_{N}}$ defined in (2.7) and where the quadratic
operator $\mathcal{Q}_{\mathcal{J}_{N}}$ is given by

$\begin{split}\mathcal{Q}_{\mathcal{J}_{N}}=\int
dxdy\,\widetilde{\Phi}(x;y)b^{*}_{x}b_{y}+\frac{1}{2}\int
dxdy\,\widetilde{\Gamma}(x;y)\big{(}b^{*}_{x}b^{*}_{y}+b_{x}b_{y}\big{)}\end{split}$

for

$\begin{split}\widetilde{\Phi}&=\gamma\big{(}-\Delta+V_{\text{ext}}-\varepsilon_{GP}\big{)}\gamma+\sigma\big{(}-\Delta+V_{\text{ext}}-\varepsilon_{GP}\big{)}\sigma\\\
&\hskip
11.38092pt+\gamma\big{(}8\pi\mathfrak{a}_{0}\varphi_{0}^{2}+N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\gamma\\\
&\hskip
11.38092pt+\sigma\big{(}8\pi\mathfrak{a}_{0}\varphi_{0}^{2}+N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\sigma\\\
&\hskip
11.38092pt+\Big{(}\gamma\big{(}N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\sigma+\emph{h.c.}\Big{)},\end{split}$
(2.18)

and

$\begin{split}\widetilde{\Gamma}&=\gamma\big{(}N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\gamma\\\
&\hskip
14.22636pt+\sigma\big{(}N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\sigma\\\
&\hskip
14.22636pt+\big{[}\sigma\big{(}-\Delta+V_{\text{ext}}-\varepsilon_{GP}\big{)}\gamma+\emph{h.c.}\big{]}\\\
&\hskip
14.22636pt+\big{[}\sigma\big{(}8\pi\mathfrak{a}_{0}\varphi_{0}^{2}+N^{3}(Vf_{\ell})(N(x-y))\varphi_{0}(x)\varphi_{0}(y)\big{)}\gamma+\emph{h.c.}\big{]}.\\\
\end{split}$ (2.19)

Moreover, the self-adjoint operator $\mathcal{E}_{\mathcal{J}_{N}}$ is bounded
by

$\pm\mathcal{E}_{\mathcal{J}_{N}}\leq
CN^{-\min\\{\frac{\tau}{2},\frac{1}{2}-\varepsilon\\}}(\mathcal{H}_{N}+1)(\mathcal{N}+1)^{3}.$
(2.20)

Compared to Prop. 2.1, notice that the decomposition (2.16) of
$\mathcal{J}_{N}$ does not contain a cubic contribution anymore. In other
words, conjugating $\mathcal{G}_{N}$ by the cubic exponential
$e^{-A}(\cdot)e^{A}$ establishes rigorously the cubic renormalization outlined
after (2.13). To control the error term $\mathcal{E}_{\mathcal{J}_{N}}$ in
(2.16), we make use of the error bound (2.20) and the following analogue of
Prop. 2.2 (see [14, Prop. 2.9] for the proof).

######  Proposition 2.4.

Assume (1.6) and let $E_{N}$ denote the ground state energy of $H_{N}$. For
some $\zeta>0$, let $\psi_{N}\in L^{2}_{s}(\mathbb{R}^{3N})$ with
$\|\psi_{N}\|=1$ be such that

$\psi_{N}={\bf 1}_{(-\infty;E_{N}+\zeta]}(H_{N})\psi_{N}.$

Let $\xi_{N}=e^{-A}e^{-B}U_{N}\psi_{N}\in\mathcal{F}_{\bot\varphi_{0}}^{\leq
N}$ be the renormalized excitation vector related to $\psi_{N}$. Then, for any
$j\in\mathbb{N}$ there exists some $C>0$ such that

$\langle\xi_{N},(\mathcal{N}+1)^{j}(\mathcal{H}_{N}+1)\xi_{N}\rangle\leq
C(1+\zeta^{j+3}).$

##  3 Diagonalization and Proof of Theorem 1.1

Using the decomposition (2.16) in Prop. 2.3, we can now determine the spectrum
$\sigma(\mathcal{J}_{N})$ of
$\mathcal{J}_{N}=e^{-A}e^{-B}U_{N}H_{N}U_{N}^{*}e^{B}e^{A}$ (and hence of
$H_{N}$, by unitary equivalence), up to errors that vanish as $N\to\infty$. We
choose the momentum cutoff parameters (introduced around (2.15))
$\varepsilon=6/13$ and $\tau=\varepsilon/6$, so that we can apply Propositions
2.1, 2.2, 2.3 and 2.4.

To obtain upper and lower bounds on $\sigma(\mathcal{J}_{N})$, we apply the
min-max principle and compare the eigenvalues of $\mathcal{J}_{N}$ with those
of the quadratic Fock space Hamiltonian

$\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}=\kappa_{\mathcal{J}_{N}}+\int
dxdy\,\widetilde{\Phi}(x,y)a_{x}^{*}a_{y}+\frac{1}{2}\int dxdy\
\widetilde{\Gamma}(x,y)(a_{x}^{*}a_{y}^{*}+\text{h.c.}),$

defined in the Fock space
$\mathcal{F}_{\bot\varphi_{0}}=\bigoplus_{n=0}^{\infty}L^{2}_{\bot\varphi_{0}}(\mathbb{R}^{3})^{\otimes_{s}n}$
(without a restriction on the number of particles). Here,
$\kappa_{\mathcal{J}_{N}}$ denotes the constant from (2.17), and the kernels
of $\widetilde{\Phi}$ and $\widetilde{\Gamma}$ were defined in (2.18) and
(2.19), respectively.

To see that it is enough to analyze
$\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}$, notice first that, to get a lower
bound on the min-max values of $\mathcal{J}_{N}$, we can drop the non-negative
potential energy $\mathcal{V}_{N}$ in (2.16) and we can control the error
$\mathcal{E}_{\mathcal{J}_{N}}$ in (2.16) through Prop. 2.4. Since
$\mathcal{F}_{\bot\varphi_{0}}^{\leq N}\subset\mathcal{F}_{\bot\varphi_{0}}$,
we obtain a lower bound on the spectrum of $\mathcal{J}_{N}$ by computing the
spectrum of $\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}$ in
$\mathcal{F}_{\bot\varphi_{0}}$. To obtain upper bounds on the min-max values
of $\mathcal{J}_{N}$, on the other hand, we construct suitable subspaces in
$\mathcal{F}_{\bot\varphi_{0}}^{\leq N}$, built from the eigenvectors of
$\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}$ truncated to
$\mathcal{F}_{\bot\varphi_{0}}^{\leq N}$. It turns out that the potential
energy $\mathcal{V}_{N}$ is negligible on such eigenspaces and that the upper
and lower bounds on the min-max values of $\mathcal{J}_{N}$ coincide, up to
errors vanishing in the limit $N\to\infty$ (see [14, Section 3] for the
details).

Hence, let us focus on $\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}$ and let us
recall how to determine its spectrum. As a quadratic Fock space Hamiltonian,
$\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}$ is exactly diagonalizable. To
diagonalize it, we follow [21] and recall the definitions of $H_{\text{GP}}$
and $E$ in (1.5) and (1.10), respectively. We then identify the operators
$H_{GP},E,\widetilde{\Phi},\widetilde{\Gamma}$ with operators mapping
$L^{2}_{\bot\varphi_{0}}(\mathbb{R}^{3})$ back into itself (so that, in
particular, $H_{\text{GP}}$ and $E$ are invertible in
$L^{2}_{\bot\varphi_{0}}(\mathbb{R}^{3})$). Moreover, we set

$\begin{split}\widetilde{D}&=\widetilde{\Phi}-\widetilde{\Gamma}\hskip
14.22636pt\text{and}\hskip
14.22636pt\widetilde{E}=\big{(}\widetilde{D}^{1/2}(\widetilde{D}+2\widetilde{\Gamma})\widetilde{D}^{1/2}\big{)}^{1/2}\end{split}$

as well as $A={\widetilde{D}}^{1/2}\widetilde{E}^{-1/2}$ and
$\alpha=\log\big{(}|A^{*}|\big{)}$. Finally, Writing $A=W|A|$ for some partial
isometry in $L_{\perp\varphi_{0}}^{2}(\mathbb{R}^{3})$, by the polar
decomposition, we construct a standard Bogoliubov transformation $\mathcal{U}$
that diagonalizes $\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}$. Indeed,
denoting by $(\varphi_{j})_{j\in\mathbb{N}}$ the eigenbasis of $\widetilde{E}$
and setting $a^{\sharp}_{j}=a^{\sharp}(\varphi_{j})$ for
$\sharp\in\\{\cdot,*\\}$ as well as
$\alpha_{ij}=\langle\varphi_{i},\alpha\,\varphi_{j}\rangle$, we set

$\mathcal{W}=\Gamma(W),\hskip
14.22636ptX=\frac{1}{2}\sum_{i,j=1}^{\infty}\alpha_{ij}a^{*}_{i}a^{*}_{j}-\text{h.c.},\hskip
14.22636pt\mathcal{U}=e^{X}\mathcal{W}.$

Here, $\Gamma(W)$ denotes the second quantization of the partial isometry $W$,
acting in the $n$-particle sector of $\mathcal{F}_{\bot\varphi_{0}}$ as
$W^{\otimes n}$. With $\mathcal{U}$ defined above, one verifies that

$\begin{split}\mathcal{U}^{*}\widetilde{\mathcal{Q}}_{\mathcal{J}_{N}}\mathcal{U}&=\kappa_{\mathcal{J}_{N}}+\frac{1}{2}\operatorname{tr}_{\bot\varphi_{0}}\bigg{(}\frac{1}{2}\big{(}\widetilde{D}^{1/2}\widetilde{E}\widetilde{D}^{-1/2}+\widetilde{D}^{-1/2}\widetilde{E}\widetilde{D}^{1/2}\big{)}-\widetilde{D}-\widetilde{\Gamma}\bigg{)}+d\Gamma(\widetilde{E}),\end{split}$
(3.1)

where $\operatorname{tr}_{\bot\varphi_{0}}$ denotes the trace in
$L^{2}_{\bot\varphi_{0}}(\mathbb{R}^{3})$ and where $d\Gamma(\widetilde{E})$
denotes the second quantization of $\widetilde{E}$, acting as
$\sum_{i=1}^{n}\widetilde{E}_{x_{i}}$ in the $n$-particle sector of
$\mathcal{F}_{\bot\varphi_{0}}$.

Now, using that the eigenvalues of $d\Gamma(\widetilde{E})$ are equal to those
of $d\Gamma(E)$ (with $E$ defined in (1.10)), up to errors that vanish as
$N\to\infty$, we conclude (1.9). Theorem 1.1 then follows by showing that the
constant on the right hand side in (3.1), that is

$\kappa_{\mathcal{J}_{N}}+\frac{1}{2}\operatorname{tr}_{\bot\varphi_{0}}\bigg{(}\frac{1}{2}\big{(}\widetilde{D}^{1/2}\widetilde{E}\widetilde{D}^{-1/2}+\widetilde{D}^{-1/2}\widetilde{E}\widetilde{D}^{1/2}\big{)}-\widetilde{D}-\widetilde{\Gamma}\bigg{)},$

is equal to the constant on the right hand side in (1.7), up to errors that
vanish as $N\to\infty$. For the details of this computation, see [14, Section
3].

## References

  * [1] A. Adhikari, C. Brennecke, B. Schlein. Bose-Einstein Condensation Beyond the Gross-Pitaevskii Regime. _Ann. Henri Poincaré_ 22, 1163-1233 (2021).
  * [2] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, E. A. Cornell. Observation of Bose-Einstein condensation in a dilute atomic vapor. _Science_ 269 (1995), 198–201.
  * [3] G. Basti, S. Cenatiempo, B. Schlein. A new second-order upper bound for the ground state energy of dilute Bose gases. _Forum Math. Sigma_ 9 (2021).
  * [4] N. Benedikter, G. de Oliveira and B. Schlein. Quantitative derivation of the Gross-Pitaevskii equation. _Comm. Pure Appl. Math._ 64, no. 8, 1399-1482 (2014).
  * [5] C. Boccato, C. Brennecke, S. Cenatiempo, B. Schlein. Complete Bose-Einstein condensation in the Gross-Pitaevskii regime. Comm. Math. Phys. 359, no. 3, 975-1026 (2018).
  * [6] C. Boccato, C. Brennecke, S. Cenatiempo, B. Schlein. The excitation spectrum of Bose gases interacting through singular potentials. _J. Eur. Math. Soc._ 22 (2020), no. 7, DOI: 10.4171/JEMS/966.
  * [7] C. Boccato, C. Brennecke, S. Cenatiempo, B. Schlein. Bogoliubov Theory in the Gross-Pitaevskii limit. _Acta Mathematica_ 222, no. 2, 219-335 (2019).
  * [8] C. Boccato, C. Brennecke, S. Cenatiempo, B. Schlein. Optimal Rate for Bose-Einstein Condensation in the Gross-Pitaevskii Regime. _Comm. Math. Phys._ 376, 1311-1395 (2020).
  * [9] N. N. Bogoliubov. On the theory of superfluidity. Izv. Akad. Nauk. USSR 11 (1947), 77. Engl. Transl. J. Phys. (USSR) 11, 23 (1947).
  * [10] L. Boßmann, S. Petrat, R. Seiringer. Asymptotic expansion of low-energy excitations for weakly interacting bosons. _Forum Math. Sigma_ 9 (2021).
  * [11] B. Brietzke, S. Fournais, J.P. Solovej. A Simple 2nd Order Lower Bound to the Energy of Dilute Bose Gases. _Comm. Math. Phys._ 376, 323-351 (2020).
  * [12] C. Brennecke, B. Schlein. Gross-Pitaevskii dynamics for Bose-Einstein condensates. _Analysis & PDE_ 12, no. 6, 1513-1596 (2019).
  * [13] C. Brennecke, B. Schlein, S. Schraven. Bose-Einstein Condensation with Optimal Rate for Trapped Bosons in the Gross-Pitaevskii Regime. Preprint: arXiv:2102.11052.
  * [14] C. Brennecke, B. Schlein, S. Schraven. Bogoliubov Theory for Trapped Bosons in the Gross-Pitaevskii Regime. Preprint: arXiv:2108.11129.
  * [15] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn and W. Ketterle. Bose-Einstein condensation in a Gas of Sodium Atoms. _Phys. Rev. Lett._ 75 (1995), no. 22, 3969–3973.
  * [16] J. Dereziński, M. Napiórkowski. Excitation Spectrum of Interacting Bosons in the Mean-Field Infinite-Volume Limit. Annales Henri Poincaré 15, 2409-2439 (2014).
  * [17] L. Erdős, B. Schlein, H.-T. Yau. Ground-state energy of a low-density Bose gas: a second order upper bound. Phys. Rev. A 78, 053627 (2008).
  * [18] S. Fournais. Length scales for BEC in the dilute Bose gas. Preprint: arXiv:2011.00309.
  * [19] S. Fournais, J.P. Solovej. The energy of dilute Bose gases. _Ann. of Math._ 192, No. 3, 893-976 (2020).
  * [20] S. Fournais, J.P. Solovej. The energy of dilute Bose gases II: The general case. Preprint: arXiv:2108.12022.
  * [21] P. Grech, R. Seiringer. The excitation spectrum for weakly interacting bosons in a trap. Comm. Math. Phys. 322, no. 2, 559-591 (2013).
  * [22] C. Hainzl. Another proof of BEC in the GP-limit. _J. Math. Phys._ 62, 051901 (2021).
  * [23] K. Huang, C. N. Yang. Quantum-Mechanical Many-Body Problem with Hard-Sphere Interaction. _Phys. Rev._ 105, no. 3, 767-757 (1957).
  * [24] T. D. Lee, K. Huang, C. N. Yang. Eigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature Properties. _Phys. Rev._ 106, no. 6, 1135-1145 (1957).
  * [25] T. D. Lee, C. N. Yang. Many-Body Problem in Quantum Mechanics and Quantum Statistical Mechanics. _Phys. Rev._ 105, 1119 - 1120 (1957).
  * [26] M. Lewin, P. T. Nam, S. Serfaty, J.P. Solovej. Bogoliubov spectrum of interacting Bose gases.  _Comm. Pure Appl. Math._ 68, 3, 413 - 471 (2014).
  * [27] E. H. Lieb and R. Seiringer.  Proof of Bose-Einstein condensation for dilute trapped gases.  _Phys. Rev. Lett._ 88, 170409 (2002).
  * [28] E. H. Lieb and R. Seiringer.  Derivation of the Gross-Pitaevskii equation for rotating Bose gases.  _Comm. Math. Phys._ 264:2, 505-537 (2006).
  * [29] E. H. Lieb, R. Seiringer, and J. Yngvason.  Bosons in a trap: A rigorous derivation of the Gross-Pitaevskii energy functional.  _Phys. Rev. A_ 61, 043602 (2000).
  * [30] P. T. Nam, M. Napiórkowski, J. Ricaud, A. Triay. Optimal rate of condensation for trapped bosons in the Gross-Pitaevskii regime. Preprint arXiv:2001.04364.
  * [31] P. T. Nam, N. Rougerie, R. Seiringer. Ground states of large bosonic systems: The Gross-Pitaevskii limit revisited. _Analysis and PDE._ 9, no. 2, 459-485 (2016).
  * [32] P. T. Nam, A. Triay. Bogoliubov excitation spectrum of trapped Bose gases in the Gross-Pitaevskii regime. Preprint arXiv:2106.11949.
  * [33] M. Napiórkowski, R. Reuvers, J. P. Solovej. The Bogoliubov free energy functional I. Existence of minimizers and phase diagrams. _Arch. Rational Mech Anal_ 229, 1037-1090 (2018).
  * [34] M. Napiórkowski, R. Reuvers, J. P. Solovej. The Bogoliubov free energy functional II. The dilute limit. _Comm. Math. Phys._ 360, 347–403 (2018).
  * [35] A. Pizzo. Bose particles in a box III. A convergent expansion of the ground state of the Hamiltonian in the mean field limiting regime. Preprint arxiv:1511.07026.
  * [36] N. Rougerie. Scaling limits of bosonic ground states, from many-body to nonlinear Schrödinger. _EMS Surv. Math. Sci._ 7 (2020), no. 2, pp. 253–408.
  * [37] R. Seiringer. The Excitation Spectrum for Weakly Interacting Bosons. Comm. Math. Phys. 306 , 565-578 (2011).
  * [38] H.-T. Yau, J. Yin. The second order upper bound for the ground state energy of a Bose gas. J. Stat. Phys. 136, no. 3, 453-503 (2009).


# Constraining $f(Q,T)$ gravity from energy conditions

Simran Arora 0000-0003-0326-8945<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus,Hyderabad-500078, India J.R.L. Santos 0000-0002-9688-938X
<EMAIL_ADDRESS>UFCG - Universidade Federal de Campina Grande -
Unidade Acadêmica de Física, 58429-900 Campina Grande, PB, Brazil. P.K. Sahoo
0000-0003-2130-8832<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus,Hyderabad-500078, India

###### Abstract

We are living in a golden age for experimental cosmology. New experiments with
high accuracy precision are been used to constrain proposals of several
theories of gravity, as it has been never done before. However, important
roles to constrain new theories of gravity in a theoretical perspective are
the energy conditions. Throughout this work, we carefully constrained some
free parameters of two different families of $f(Q,T)$ gravity using different
energy conditions. This theory of gravity combines the gravitation effects
through the non-metricity scalar function $Q$, and manifestations from the
quantum era of the Universe in the classical theory (due to the presence of
the trace of the energy-momentum tensor $T$). Our investigation unveils the
viability of $f(Q,T)$ gravity to describe the accelerated expansion our
Universe passes through. Besides, one of our models naturally provides a
phantom regime for dark energy and satisfies the dominant energy condition.
The results here derived strength the viability of $f(Q,T)$ as a promising
complete theory of gravity, lighting a new path towards the description of the
dark sector of the Universe.

$f(Q,T)$ gravity; Energy Conditions; Equation of state; Phantom regime.

##  I Introduction

Since the remarkable measurements from Supernova Cosmology Project riess_98
and High Redshift Supernova Team perl_99 , we have consciously known that our
Universe passes through an accelerated phase of expansion, whose agent is
named dark energy. Dark energy corresponds to approximately $70\,\%$ of the
content of the so-called dark sector of the Universe, and its understanding is
one of the actual biggest problems in science. A simple path to describe the
nature of the dark energy consists in to add a cosmological constant to
Einstein’s General Relativity (GR), yielding to the $\Lambda$CDM model.
However, the cosmological constant brings several other issues related to its
nature. Among them, we highlight the cosmic coincidence problem, and its huge
discrepancy between cosmological observations and quantum field theory
predictions, which is about $120$ orders of magnitude adler_95 .

Apart from these listed problems, GR stills the most well succeed theory to
describe the Universe. It was confirmed by several surveys such as PLANCK
Collaboration Planck/2018 , Dark Energy Survey Des/2018 , besides the recent
beautiful measurements of gravitational waves from LIGO/VIRGO Collaboration
lv_papers , and the first image of a black hole obtained by the Event Horizon
Telescope eht_papers . Although, GR does not yield to a renormalizable quantum
theory for gravity, opening space to several alternative theories desiring to
describe gravity at a quantum level.

A promising theory of gravity was introduced by Jimenez et al. Jimenez/2018 ,
and called symmetric teleparallel gravity or $f(Q)$, where the gravitation
interaction is mediated by the non-metricity term $Q$. Such a theory rapidly
inspired several works and it has been constantly tested. Among such tests, we
highlight the work done by Lazkoz et al Lazkoz/2019 , were several $f(Q)$
models were constrained through redshift comparison with data from the
expansion rate, Type Ia Supernovae, Quasars, Gamma-Ray Bursts, Baryon Acoustic
Oscillations data, and Cosmic Microwave Background distance. Besides, the
$f(Q)$ model also unveiled a compatible description of an accelerated phase
when submitted to energy conditions constraints as shown in Mandal/2020 .

As Jimenez noted in his work, $f(Q)$ theories share the background equations
with $f(T)$ theories, which at the perturbative level often lead to strong
coupling problems. It is shown that in $f(Q)$ models on a general FLRW
background, the strong coupling issues experienced in $f(T)$ theories are
absent. However, on maximally symmetric backgrounds such as Minkowski and de-
Sitter, they do appear. It can be ensured that the same coupling issues can
also be found in $f(Q,T)$ on a maximally symmetric background. Similar to
$f(Q)$, two scalar degrees of freedom that are absent in $f(T)$ are also
propagated on large scales by the $f(Q,T)$ model. These two degrees of freedom
are the ones that vanish around maximally symmetric backgrounds, and thus
cause the problem of strong-coupling. Second, a gauge symmetry given by a
restricted diffeomorphism is retained by the maximally symmetric backgrounds.
As a result of a residual gauge symmetry that then roots the strong coupling
issues, these findings allow us to provide a better understanding of the
disappearance of degrees of freedom around these backgrounds.

In the search for a complete theory of gravity emerged the $f(Q,T)$, recently
presented by Yixin Xu et al. Yixin/2019 . Such a theory couples the
gravitation effects through the non-metricity function $Q$, and manifestations
from the quantum domain in the classical theory (due to the presence of the
trace of the energy-momentum tensor $T$). This new formulated $f(Q,T)$ gravity
is constructed in a manner similar to $f(R,T)$Harko_Lobo/2018 ; Harko/2011 ,
but with the standard Ricci scalar replaced by the non-metricity that defines
the gravity symmetric teleparallel formulation. In $f(Q,T)$ theory, the
coupling between $Q$ and $T$ contributes to the non-conservation of the
energy-momentum tensor, similar to the regular curvature trace of the energy-
momentum tensor couplings. Many studies on couplings between matter and non-
metricity are also conducted. T. Harko et al.T/2018 studied coupling matter in
modified $Q$ gravity, novel couplings between nonmetricity and
matterTharko/2019 . Minimal coupling in the presence of non-metricity and
torsion is investigated by A. DelhomDelhom/2020 . Thomas P SotiriouThomas/2009
also worked on $f(R)$ gravity, torsion and non-metricity.
The $f(Q,T)$ gravity has been shown compatible with the accelerating expansion
phase Yixin/2019 , besides it is also in agreement with important different
phases our Universe passes through as the baryogenesis Sahoo/epjc_2020 .
Moreover, such a theory is compatible with measurements of the Hubble
parameter for different redshifts as one can see in Simran/2020 . Beyond these
successful tests, an important role any alternative theory of gravity should
obey is the energy condition constraints Capozziello/2018 . These constraints
are crucial to determine the proper regimes allowed for a new theory of
gravity, to describe its attractive nature, and to assign the causal and the
geodesic structure of space-time. Furthermore, the energy conditions also
allow us to confront a new theory of gravity against $\Lambda$CDM model.

Therefore, in this work, we intend to study carefully all the energy
conditions constraints on different forms of $f(Q,T)$ gravity. Our analyses
were carried using the actual values of the Hubble, and the deceleration
parameters. The energy conditions enable us to impose constraints over our
free parameters, unveiling the viability of the $f(Q,T)$ models. We also
verify the compatibility of our results with $\Lambda$CDM model. The
discussions along this study are organized in the following nutshell: in
section II we introduce generalities about the $f(Q,T)$ gravity. In section
III we use the Raychaudhuri equations to find our energy conditions embedding
the non-metricity and the trace of the energy-momentum tensor contributions.
The constraints on $f(Q,T)$ models are discussed in details in section IV. A
comparison between the $f(Q,T)$ models and the $\Lambda$CDM model is presented
in section V, where we also depicted the equation of state parameters for the
models here studied. Section VI is dedicated to our final remarks and
perspectives.

##  II Overview of $f(Q,T)$ Gravity

The $f(Q,T)$ gravity is described through the following action Yixin/2019 ,

$S=\int\left(\frac{1}{16\pi}f(Q,T)+L_{m}\right)d^{4}x\sqrt{-g}.$ (1)

where $f$ is an arbitrary function of the non-metricity $Q$, and of the trace
of the energy-momentum tensor $T$, besides $L_{m}$ represents the Lagrangian
of a given matter, and $g=det(g_{\alpha\beta})$. As it was discussed by
Jimenez et al. Jimenez/2018 , the non-metricity function is such that

$Q\equiv-g^{\alpha\beta}(L^{\mu}_{\,\,\nu\alpha}L^{\nu}_{\,\,\beta\mu}-L^{\mu}_{\,\,\nu\mu}L^{\nu}_{\,\,\alpha\beta}),$
(2)

where $L^{\mu}_{\,\,\nu\gamma}$ is the disformation tensor whose explicit form
is

$L^{\mu}_{\,\,\nu\gamma}=-\frac{1}{2}g^{\mu\lambda}(\nabla_{\gamma}g_{\nu\lambda}+\nabla_{\nu}g_{\lambda\gamma}-\nabla_{\lambda}g_{\nu\gamma}).$
(3)

Another key ingredient to describe the symmetric teleparallel is the non-
metricity tensor, which is defined as

$Q_{\gamma\mu\nu}=\nabla_{\gamma}g_{\mu\nu}\,,$ (4)

and whose traces are

$Q_{\mu}={{Q_{\mu}}^{\alpha}}_{\alpha}\,;\qquad\widetilde{Q}_{\mu}=Q^{\alpha}_{\,\alpha\,\mu}\,.$
(5)

Here $Q_{\mu}$ is refereed as the Weyl vector, while $\widetilde{Q}_{\mu}$ is
called the second non-metricity vector. We can also define a superpotential
related with the non-metricity tensor as

$4{P^{\mu}}_{\alpha\beta}=-{Q^{\mu}}_{\alpha\beta}+2{Q_{(\alpha}}^{\,\mu}_{\,\,\,\beta)}-Q^{\mu}g_{\alpha\beta}-\widetilde{Q}^{\mu}g_{\alpha\beta}-\delta^{\mu}_{(\alpha}Q_{\beta)}\,,$
(6)

yielding to the quadratic form for the non-metricity function Jimenez/2018

$Q=-Q_{\mu\alpha\beta}P^{\mu\alpha\beta}\,.$ (7)

Moreover, as it is known the energy-momentum tensor can be written as

$T_{\alpha\beta}=-\frac{2}{\sqrt{-g}}\dfrac{\delta(\sqrt{-g}L_{m})}{\delta
g^{\alpha\beta}}$ (8)

and its variation in respect to the metric tensor is such that

$\frac{\delta\,g^{\,\mu\nu}\,T_{\,\mu\nu}}{\delta\,g^{\,\alpha\,\beta}}=T_{\,\alpha\beta}+\Theta_{\,\alpha\,\beta}\,,$
(9)

where

$\Theta_{\alpha\beta}=g^{\mu\nu}\frac{\delta T_{\mu\nu}}{\delta
g^{\alpha\beta}}.$ (10)

Therefore, taking the variation of action (1) with respect to the metric, we
find the field equations

$8\pi
T_{\alpha\beta}=-\frac{2}{\sqrt{-g}}\nabla_{\mu}(f_{Q}\sqrt{-g}P^{\mu}_{\,\,\alpha\beta}-\frac{1}{2}fg_{\alpha\beta}\\\
+f_{T}(T_{\alpha\beta}+\Theta_{\alpha\beta})-f_{Q}(P_{\alpha\mu\nu}Q_{\beta}^{\,\,\,\mu\nu}-2Q^{\mu\nu}_{\,\,\,\alpha}P_{\mu\nu\beta}).$
(11)

where $f_{Q}=\dfrac{df}{dQ}$.

Now, let us assume a Friedmann-Lemaitre-Robertson-Walker metric (FLRW) given
by,

$ds^{2}=-N^{2}(t)dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j},$ (12)

where $a(t)$ is the scale factor of the Universe and $N(t)$ is the lapse
function. Jiménez et al. Jose stated that if we choose the coincident gauge
$\Gamma^{\alpha}_{\,\,\,\mu\nu}=0$, the diffeomorphism gauge’s symmetry gets
compromised in the generic theory. Hence, we cannot arbitrarily choose the
time parameterization as used in General Relativity. The non-metricity scalar
is defined to be $Q=6\frac{H^{2}}{N^{2}}$. As we used diffeomorphism to repair
the coincident gauge, we cannot choose any unique lapse function. However, as
$Q$ maintains a residual time reparametrization invariance stated in
Jimenez/2018 , we are allowed to set $N=1$, as a matter of simplicity.
Therefore, we yield to the metric,

$ds^{2}=-dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j}\,.$ (13)

Furthermore, the matter content of the Universe is assumed as been a perfect
fluid, whose energy-momentum tensor is $T_{\,\alpha\beta}=diag(-\rho,p,p,p)$.

Therefore, substituting Eqs. (7), (8), and (10) into (9), we yield to the
modified Friedmann equations for such a theory, which are explicitly
represented as

$3\,H^{\,2}=\frac{1}{2\,F}\,\left(-8\,\pi\,\rho+\frac{f}{2}-\frac{2\,\overline{G}}{1+\overline{G}}\,\left(\dot{F}\,H+F\,\dot{H}\right)\right)\,,$
(14)
$\dot{H}+3\,H^{\,2}=\frac{8\,\pi\,p}{2\,F}+\frac{f}{4\,F}-\frac{\dot{F}}{F}\,H\,.$
(15)

From the Friedmann equations, one can prove that the density and pressure also
satisfy the generalized energy balance equation for $f(Q,T)$ gravity
Yixin/2019 , whose form is

$\dot{\rho}+3\,H\,(\rho+p)=B_{\mu}\,v^{\mu}\,,$ (16)

where

$\displaystyle B_{\mu}v^{\mu}$ $\displaystyle=$
$\displaystyle\frac{\overline{G}}{16\,\pi\,\left(1+\overline{G}\right)\,\left(1+2\overline{G}\right)}$
$\displaystyle\times\left[\dot{S}-\frac{\left(3\overline{G}+2\right)\dot{\overline{G}}}{\left(1+\overline{G}\right)\overline{G}}\,S+6\,H\,S\right]\,,$

with $S=2\,\left(\dot{F}\,H+\dot{H}\,F\right)$ .

The previous equations can be rewritten in analogy to Einstein’s general
relativity in the following way

$3\,H^{\,2}=-\frac{8\,\pi\,}{2}\,\rho_{eff}\,,$ (18)
$\dot{H}+3\,H^{\,2}=\frac{8\,\pi}{2}\,p_{eff}\,,$ (19)

resulting in

$\rho_{eff}=\frac{\rho}{F}-\frac{f}{16\,\pi\,F}+\frac{\overline{G}}{1+\overline{G}}\,\frac{\dot{F}\,H+F\,\dot{H}}{4\,\pi\,F}\,,$
(20)
$p_{eff}=\frac{p}{F}+\frac{f}{16\,\pi\,F}-\frac{\dot{F}}{4\,\pi\,F}\,H\,,$
(21)

as the effective density, and pressure. Here ($\cdot$) dot represents a
derivative with respect to time, besides $F=f_{Q}$, and
$8\pi\overline{G}=f_{T}$ denote differentiation with respect to $Q$, and $T$,
respectively. Moreover, we are able to observe that the contributions coming
from the $f(Q,T)$ model are embedded into $\rho_{eff}$, and $p_{eff}$. We also
highlight that this same identification procedure was applied to generalized
Gauss-Bonnet gravity by Kaczmarek et. al. Kaczmarek/2020 , and also in $f(Q)$
gravity Mandal/2020 . Another interesting feature related with $\rho_{eff}$,
and $p_{eff}$ is that they obey the conservation equation

$\dot{\rho}_{eff}+3\,H\,\left(\rho_{eff}+p_{eff}\right)=0\,.$ (22)

Such a behavior is going to be used in the next section, to address the
physical interpretations related with the different energy conditions.

##  III Weyl-type $f(Q,T)$ and Raychaudhuri Equation

Energy conditions in modified gravity are the tools which empower the casual
and geodesic structure of space-time. These conditions are formulated with the
help of Raychaudhuri equations which describe the action of congruence and
attractiveness of the gravity for timelike, spacelike, or lightlike curves.

In order to study the implications of the non-metricity, and of the trace of
the energy-momentum tensor for the Raychaudhuri equation, we are going to
derive such an equation in the Weyl framework. The Weyl geometry is based on a
connection where the orientation and the magnitude of a vector can change,
under a parallel transportation. Therefore, such a geometry is a natural
framework to describe the effects of torsion and of the non-metricity on the
particle dynamics.

Recently, Yixin Xu et al. Yixin/2020 extended the formulation of $f(Q,T)$ in
the framework of Weyl geometry. There, the authors were able to derive the
general field equations for gravity in the Weyl geometry coupled to the matter
energy-momentum tensor, as well as, to analyze the cosmological implications
of different families of $f(Q,T)$ gravity. In another recent work Jin-Zhao
Yang et al. Yang/2021 , complemented the formulation raised in Yixin/2020 and
also in the beautiful work of Iosifidis et al. Damianos/2018 , studying the
geodesic deviation equation, the Raychaudhuri equation, and the tidal forces
for the $f(Q,T)$ gravity in Weyl framework. In order to implement our set of
energy conditions for the $f(Q,T)$ models introduced in Yixin/2019 , based on
the constraints over the Raychaudhuri equation, let us briefly reproduce some
definitions and results from Yixin/2020 ; Yang/2021 ; Damianos/2018 .

The Weyl geometry is a generalization of the Riemannian geometry, where the
connection is redefined as

$\widetilde{\Gamma}^{\lambda}_{\,\mu\,\nu}\equiv\Gamma^{\lambda}_{\,\mu\,\nu}+g_{\mu\nu}w^{\lambda}-\delta^{\lambda}_{\mu}w_{\nu}-\delta^{\lambda}_{\nu}w_{\mu}\,,$
(23)

where $w_{\mu}$ is the so-called Weyl vector, and
$\Gamma^{\lambda}_{\,\mu\,\nu}$ is determined through the metric
$g_{\mu\,\nu}$. Taking the Weyl connection, one is able to show that the
covariant derivative applied to the metric tensor is

$\widetilde{\nabla}_{\lambda}\,g_{\,\mu\nu}=2\,w_{\lambda}g_{\,\mu\nu}\,;\qquad\widetilde{\nabla}_{\lambda}\,g^{\,\mu\nu}=-2\,w_{\lambda}\,g^{\,mu\nu}\,.$
(24)

Therefore, in this framework the non-metricity tensors are written as

$\displaystyle
Q_{\,\lambda\mu\nu}=-\widetilde{\nabla}_{\lambda}\,g_{\,\mu\nu}=-2\,w_{\lambda}g_{\,\mu\nu}\,;$
(25) $\displaystyle
Q^{\,\lambda\mu\nu}=\widetilde{\nabla}^{\lambda}\,g^{\,\mu\nu}=-2\,w^{\lambda}g^{\,\mu\nu}\,.$

Such relations also yield to

$Q=-6w^{2}\,,$ (26)

as the non-metricity scalar.

The presence of the non-metricity affects the length of vectors in Weyl
framework, when they are submitted to parallel transport. In such a case, the
four velocity is given by

$u_{\mu}\,u^{\mu}=g_{\mu\nu}u^{\,\mu}u^{\,\nu}=-l^{2}\,,\qquad
l=l(x^{\alpha})\,,$ (27)

where $u=dx^{\mu}/d\lambda$ with $\lambda$ as the affine parameter, and
$l(x^{\alpha})$ is an arbitrary function of space and time coordinates.
Moreover, in the presence of the non-metricity the associated projection
tensor has the form

$h_{\mu\nu}=g_{\mu\nu}+\frac{1}{l^{2}}\,u_{\mu}\,u_{\nu}\,,$ (28)

which results in the following properties:

$\displaystyle h_{\mu\nu}=h_{\nu\mu}\,;$ (29) $\displaystyle
h_{\mu\nu}u^{\mu}=0\,;$ $\displaystyle h_{\mu\nu}h^{\mu\nu}=3$ $\displaystyle
h_{\mu\lambda}h^{\lambda\nu}=h_{\mu}^{\nu}=\delta_{\mu}^{\nu}+\frac{1}{l^{2}}u_{\mu}u^{\nu}\,.$

Besides, the non-metricity also implies in two different acceleration vectors,
since raising and lowering indexes do not commute. Such acceleration vectors
are defined as Damianos/2018

$A^{\mu}=u^{\prime\,\mu}=u^{\lambda}\widetilde{\nabla}_{\lambda}u^{\mu}\,,$
(30)

as the contravariant or 4-acceleration vector, and

$a_{\mu}=u^{\prime}_{\,\mu}=u^{\lambda}\widetilde{\nabla}_{\lambda}u_{\mu}\,,$
(31)

as the covariant or hyper 4-acceleration vector. In these last equations,
prime means a derivative in respect to the affine parameter $\lambda$. These
definitions yield to the following constraint between the acceleration
vectors,

$A^{\mu}=a^{\mu}+Q^{\,\nu\lambda\mu}u_{\nu}u_{\lambda}\,.$ (32)

From Eqs. (30) and (31), we are able to derive that

$A^{\mu}u_{\mu}=-\frac{(l^{2})^{\prime}}{2}+\frac{1}{2}\,Q_{\mu\nu\lambda}u^{\mu}u^{\nu}u^{\lambda}\,,$
(33)

and

$a^{\mu}u_{\mu}=-\frac{(l^{2})^{\prime}}{2}-\frac{1}{2}\,Q_{\mu\nu\lambda}u^{\mu}u^{\nu}u^{\lambda}\,.$
(34)

Therefore, we directly observe that

$\left(A^{\mu}+a^{\mu}\right)u_{\mu}=-(l^{2})^{\prime}\,;\qquad\left(A^{\mu}-a^{\mu}\right)u_{\mu}=Q_{\mu\nu\lambda}u^{\mu}u^{\nu}u^{\lambda}\,.$
(35)

Another key ingredients to build the Raychaudhuri equation in this framework
are the covariant derivatives of $u_{\mu}$ and $a_{\mu}$ vectors, whose
temporal and spacial components are

$\widetilde{\nabla}_{\nu}u_{\mu}={\cal
D}_{\nu}u_{\mu}-\frac{1}{l^{2}}\left(u_{\mu}\xi_{\nu}+a_{\mu}u_{\nu}\right)-\frac{1}{l^{4}}\left(u^{\lambda}a_{\lambda}\right)u_{\mu}u_{\nu}\,,$
(36) $\widetilde{\nabla}^{\mu}a_{\mu}={\cal
D}^{\mu}a_{\mu}+\frac{1}{l^{2}}\,A^{\mu}a_{\mu}-\frac{1}{l^{2}}a_{\mu}u^{\mu}\,,$
(37)

where ${\cal
D}_{\nu}u_{\mu}=h_{\nu}^{\varphi}h_{\mu}^{\lambda}\widetilde{\nabla}_{\varphi}u_{\lambda}$
is the projected covariant derivative, besides
$\xi_{\mu}=u^{\nu}\widetilde{\nabla}_{\mu}u_{\nu}$, and
$\xi_{\mu}u^{\mu}=a_{\mu}u^{\mu}$. The projected covariant derivative of
vector $u_{\mu}$ can be decomposed as

$D_{\nu}u_{\mu}=\frac{1}{3}\,\left(\theta+\frac{1}{l^{2}}a_{\lambda}u^{\lambda}\right)\,h_{\mu\nu}+\sigma_{\mu\nu}+\omega_{\mu\nu}\,,$
(38)

for $\theta$, $\sigma_{\mu\nu}$, and $\omega_{\mu\nu}$ as the expansion,
shear, and rotation, associated to the vector field $u^{\mu}$, respectively.

The previous definitions allow us to write the following geodesic equation

$A^{\mu}=\frac{d^{2}x^{\mu}}{d\lambda^{2}}+\widetilde{\Gamma}^{\mu}_{\nu\lambda}u^{\nu}u^{\lambda}=f^{\mu}\,,$
(39)

where $f^{\mu}$ is interpreted as an extra force induced by the non-minimal
coupling between $Q$ and $T$ at the $f(Q,T)$ gravity.

The Raychaudhuri equation in this framework can be derived from the curvature
tensor

$\left(\widetilde{\nabla}_{\mu}\widetilde{\nabla}_{\nu}-\widetilde{\nabla}_{\nu}\widetilde{\nabla}_{\mu}\right)u_{\lambda}=-\widetilde{R}_{\beta\lambda\nu\mu}u^{\beta}\,,$
(40)

where we are not considering torsion effects. Contracting the last equation
with $g^{\lambda\nu}u_{\mu}$ yields to

$g^{\lambda\nu}u_{\mu}\left(\widetilde{\nabla}_{\mu}\widetilde{\nabla}_{\nu}-\widetilde{\nabla}_{\nu}\widetilde{\nabla}_{\mu}\right)u_{\lambda}=-\widetilde{R}_{\beta\lambda\nu\mu}u^{\beta}u^{\mu}g^{\lambda\nu}\,.$
(41)

With these relations in hand, one can determine that the Raychaudhuri equation
is explicitly written as

$\displaystyle\theta^{\prime}=-\frac{\theta^{2}}{3}-R_{\mu\nu}u^{\mu}u^{\nu}-\sigma_{\mu\nu}\sigma^{\mu\nu}+\omega_{\mu\nu}\omega^{\mu\nu}$
(42)
$\displaystyle+\widetilde{\nabla}^{\mu}a_{\mu}-\frac{2}{3\,l^{2}}\theta\,a_{\mu}u^{\mu}+\frac{2}{3\,l^{4}}\left(a_{\mu}u^{\mu}\right)^{2}+\frac{2}{l^{2}}a_{\mu}\xi^{\mu}-\widetilde{Q}^{\prime}_{\mu}u^{\mu}$
$\displaystyle+\frac{1}{3}\left(\theta+\frac{1}{l^{2}}a_{\nu}u^{\nu}\right)\left(Q_{\mu}-\widetilde{Q}_{\mu}\right)\,u^{\mu}-Q_{\mu\nu\lambda}\left(\sigma^{\mu\nu}+\omega^{\mu\nu}\right)u^{\lambda}$
$\displaystyle-\frac{1}{l^{2}}\,Q_{\mu\nu\lambda}u^{\mu}u^{\nu}\left(a^{\lambda}+\xi^{\lambda}\right)+Q_{\mu\nu\lambda}u^{\mu}\sigma^{\nu\lambda}$
$\displaystyle+\frac{1}{l^{2}}Q_{\mu\nu\lambda}\left(u^{\mu}\xi^{\nu}+a^{\mu}u^{\nu}\right)u^{\lambda}+u^{\mu}u^{\nu}\widetilde{\nabla}^{\lambda}Q_{\mu\nu\lambda}+Q_{\mu}^{\lambda\beta}Q_{\beta\lambda\nu}u^{\mu}u^{\nu}\,.$

The detailed calculation of the previous equation can be found in Refs.
Yang/2021 ; Damianos/2018 . Now, taking the Weyl and the second non-metricity
vectors from Eq. (5), and constraining the Raychaudhuri equation for
autoparallel curves, where particle have zero path acceleration
$\left(A^{\mu}=f^{\mu}=0\right)$ Damianos/2018 , Eq. (42) is reduced to

$\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{\prime}=-\frac{1}{3}\,\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{2}-R_{\mu\nu}u^{\mu}u^{\nu}-2\,\left(\sigma^{2}-\omega^{2}\right)\,.$
(43)

Then, following the procedure introduced by Iosifidis et al. Damianos/2018 ,
we can assume that for an irrotational and shear-free scenario, (43) yields to
the constraint

$\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{\prime}+\frac{1}{3}\,\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{2}\leq
0\,,$ (44)

if

$R_{\mu\nu}u^{\mu}u^{\nu}\geqslant 0\,,$ (45)

resulting in an attractive behavior for gravity. Therefore, this approach
unveils a generalization of the constraints for the energy conditions. Note
that by changing $\theta-2\,\frac{l^{\prime}}{l}\rightarrow\theta$ into (43),
we can recover the standard Raychaudhuri equation for GR. Moreover, the same
procedure here discussed can be applied to the case of a null-vector
$k^{\mu}$, leading to the following reduced form of the Raychaudhuri equation

$\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{\prime}=-\frac{1}{3}\,\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{2}-R_{\mu\nu}k^{\mu}k^{\nu}-2\,\left(\sigma^{2}-\omega^{2}\right)\,.$
(46)

Consequently, we are able to establish the constraint

$\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{\prime}+\frac{1}{3}\,\left(\theta-2\,\frac{l^{\prime}}{l}\right)^{2}\leq
0\,,$ (47)

for

$R_{\mu\nu}k^{\mu}k^{\nu}\geqslant 0\,,$ (48)

generalizing the approach of Santos et al. Santos .

Following the methodology presented in Mandal/2020 , the $f(Q,T)$ models are
going to be restricted to the set of energy conditions bellow

  * •

null energy condition (NEC) $\Leftrightarrow\rho_{eff}+p_{eff}\geq 0$;

  * •

weak energy condition (WEC) $\Leftrightarrow\rho_{eff}+p_{eff}\geq 0$ and
$\rho_{eff}\geq 0$;

  * •

dominant energy condition (DEC) $\Leftrightarrow\rho_{eff}\geq|p_{eff}|$ and
$\rho_{eff}\geq 0$;

  * •

strong energy condition (SEC) $\Leftrightarrow\rho_{eff}+3\,p_{eff}\geq 0$.

Through Eq. (22), we are able to observe an expanding behavior for the
Universe where energy density always decreases, if the null energy condition
is obeyed Rubakov/2014 . Such a behavior is also corroborated by weak energy
condition which also establishes that the local energy density should be
always positive, which means that local observers measure positive mass
densities Tipler/1978 . Beyond the expanding behavior constrained by WEC, the
dominant energy condition also imposes that the pressure cannot exceed energy
density, in order to state that matter flows along timelike or null world
lines Capozziello/2018 . Moreover, when we deal with a perfect fluid, in this
case mapped by $\rho_{eff}$, and $p_{eff}$, the strong energy condition states
that we must have, besides $\rho_{eff}+p_{eff}\geq 0$, the sum
$\rho_{eff}+3\,p_{eff}\geq 0$, which means that SEC includes NEC and excludes
large negative pressures regimes Crislane/2017 . Therefore, in order to
describe a Universe that is dominated by negative pressure, the SEC condition
needs to be violated Capozziello/2018 .

Then, by substituting Eqs. (11), and (13) in the previous relations, we
established the following set of energy conditions

  * •

null energy condition (NEC) $\Leftrightarrow\rho+p\geq 0$;

  * •

weak eneregy condition (WEC) $\Leftrightarrow\rho+p\geq 0$;

  * •

dominant energy condition (DEC) $\Leftrightarrow\rho\geq|p|$.

Moreover, WEC, DEC and SEC energy conditions demand the extra constraints

$\mbox{DEC}\Leftrightarrow F\leq 0\,;$ (49) $\mbox{WEC and
DEC}\Leftrightarrow\rho-\frac{f}{16\,\pi}+\frac{\overline{G}}{1+\overline{G}}\,\frac{\dot{F}\,H+F\,\dot{H}}{4\,\pi}\geq
0\,,$ (50)

and

$\mbox{SEC}\Leftrightarrow\rho+3\,p+\frac{f}{8\,\pi}-\frac{3-\overline{G}}{1+\overline{G}}\,\frac{\dot{F}\,H}{4\,\pi}+\frac{\overline{G}}{1+\overline{G}}\,\frac{F\,\dot{H}}{4\,\pi}\geq
0\,.$ (51)

Through these constraints, we can realize how different $f(Q,T)$ models modify
the standard energy conditions derived from the Raychaudhuri equations.

##  IV Constraining $f(Q,T)$ Gravity Models

In the framework of FRW metric, one can also use the constraints of energy
conditions to restrict certain models in $f(Q,T)$ gravity. A cosmological
quantity which is essential to properly describe the energy conditions in a
phenomenological perspective is the deceleration parameter, whose definition
is Sharif

$q=-\frac{1}{H^{2}}\frac{\ddot{a}}{a}\,.$ (52)

Alternatively, the time derivative of the Hubble parameter can be rewritten as

$\dot{H}=-H^{2}(1+q).$ (53)

Beyond these ingredients, in order to constraint the energy conditions with
phenomenological observations, we are going to consider that
$H=H_{0}=67.9\,km\,s^{\,-1}\,Mpc^{\,-1}$, and $q=q_{0}=-0.503$ as the present
values for the Hubble, and the deceleration parameters, respectively
Planck/2018 ; Capozziello/2019 .

###  IV.1 $f(Q,T)=mQ+bT$

As a first model, we are going to work with $f(Q,T)=mQ+bT$, where $m$, and $b$
are free parameters. This model was introduced by Xu et al. Yixin/2019 , and
it naturally describes an exponential expanding Universe, with $\rho\propto
e^{-H_{0}\,t}$ Yixin/2019 . The present model yields to $F=m$ and
$8\pi\overline{G}=b$, besides, by taking the modified Friedmann equations
(14), and (15) together with Eqs. (52), (53) into the energy conditions, we
find the constraints

$NEC\Leftrightarrow-\frac{H_{0}^{2}m(b(q_{0}+4)+24\pi)}{2\left(b^{2}+12\pi
b+32\pi^{2}\right)}-\frac{H_{0}^{2}m(3bq_{0}+8\pi(2q_{0}-1))}{2\left(b^{2}+12\pi
b+32\pi^{2}\right)}\geq 0\,,$ (54)

$WEC\Leftrightarrow-\frac{H_{0}^{2}m(b(q_{0}+4)+24\pi)}{2\left(b^{2}+12\pi
b+32\pi^{2}\right)}-\frac{H_{0}^{2}m(3bq_{0}+8\pi(2q_{0}-1))}{2\left(b^{2}+12\pi
b+32\pi^{2}\right)}\geq 0\,,\qquad and\qquad-\frac{3mH_{0}^{2}}{4\pi}\geq
0\,,$ (55)

$DEC\Leftrightarrow\frac{H_{0}^{2}m(3bq_{0}+8\pi(2q_{0}-1))}{2\left(b^{2}+12\pi
b+32\pi^{2}\right)}-\frac{H_{0}^{2}m(b(q_{0}+4)+24\pi)}{2\left(b^{2}+12\pi
b+32\pi^{2}\right)}\geq 0\,,\qquad and\qquad-\frac{3mH_{0}^{2}}{4\pi}\geq
0\,,$ (56)

$SEC\Leftrightarrow(8\pi)^{-1}\left(6H_{0}^{2}m-\frac{2bH_{0}^{2}m(q_{0}+1)}{b+8\pi}+\frac{2H_{0}^{2}m(-2bq_{0}+b-24\pi
q_{0})}{b+8\pi}\right)\geq 0.$ (57)

Figure 1: Density parameter, and Energy conditions for $f(Q,T)=mQ+bT$. The
graphics were depicted with the present values of $H_{0}$ and $q_{0}$
parameters.

We can observe the behavior of $\rho$, WEC, DEC, and SEC energy conditions in
the graphics depicted in Fig. 1. There we realize how the density $\rho$
decreases for specific values of $m$, and $b$, corroborating with the
exponential expansion behavior derived by Xu et al. Yixin/2019 . Moreover, the
energy conditions allow us to constrain the free parameters $m$, and $b$.
Through DEC we found that $m$ should be negative. Also (54), (55), and (57)
assure the range of model parameters as $b>-4\,\pi$ and $m\leq 0$, satisfying
NEC, WEC, and DEC. Furthermore, the constrained parameters result in the
violation of SEC, which is compatible with the accelerated phase our Universe
passes through Visser/2000 . Another remarkable feature coming from the energy
conditions, is that the constrained parameters $m$, and $b$ corroborate, and
fine tune the observational bounds for $f(Q,T)$ gravity investigated in
Simran/2020 .

###  IV.2 $f(Q,T)=Q^{n+1}+bT$

As a second model, let us deal with $f(Q,T)=Q^{n+1}+bT$ where $n$ and $b$ are
free parameters. Such a model was proposed by Xu et al. Yixin/2019 , and
considers non-linear contributions due to the torsion in the gravity sector,
moreover, it was constrained through measurements of the Hubble parameter for
different redshifts Simran/2020 . This specific model also can describe an
accelerating Universe with $\rho\propto t^{\,2\,(n-1)}$ as one can see in
Yixin/2019 . The present model yields to $f=f_{Q}=(n+1)Q^{n}$ and
$8\pi\overline{G}=f_{T}=b$, and by working with Eqs. (14), (15), (52), and
(53) we are able to derive the energy conditions below

$NEC\Leftrightarrow-\frac{2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(b(nq_{0}+n+q_{0}+4)+24\pi)}{b^{2}+12\pi
b+32\pi^{2}}-\frac{\displaystyle\begin{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(3b(nq_{0}+n+q_{0})\\\
+8\pi(2n(q_{0}+1)+2q_{0}-1))\end{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(3b(nq_{0}+n+q_{0})\\\
+8\pi(2n(q_{0}+1)+2q_{0}-1))}{b^{2}+12\pi b+32\pi^{2}}\geq 0\,,$ (58)

$WEC\Leftrightarrow-\frac{2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(b(nq_{0}+n+q_{0}+4)+24\pi)}{b^{2}+12\pi
b+32\pi^{2}}-\frac{\displaystyle\begin{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(3b(nq_{0}+n+q_{0})\\\
+8\pi(2n(q_{0}+1)+2q_{0}-1))\end{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(3b(nq_{0}+n+q_{0})\\\
+8\pi(2n(q_{0}+1)+2q_{0}-1))}{b^{2}+12\pi b+32\pi^{2}}\geq 0\,,\\\ \qquad
and\qquad-\frac{2^{n-2}3^{n+1}(n+1)\left(H_{0}^{2}\right)^{n+1}}{\pi}\geq
0\,,$ (59)

$DEC\Leftrightarrow\frac{\displaystyle\begin{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(3b(nq_{0}+n+q_{0})\\\
+8\pi(2n(q_{0}+1)+2q_{0}-1))\end{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}(3b(nq_{0}+n+q_{0})\\\
+8\pi(2n(q_{0}+1)+2q_{0}-1))}{b^{2}+12\pi
b+32\pi^{2}}-\frac{\displaystyle\begin{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}\\\
(b(nq_{0}+n+q_{0}+4)+24\pi)\end{multlined}2^{n-1}3^{n}(2n+1)\left(H_{0}^{2}\right)^{n+1}\\\
(b(nq_{0}+n+q_{0}+4)+24\pi)}{b^{2}+12\pi b+32\pi^{2}}\geq 0\,,\\\ \qquad
and\qquad-\frac{2^{n-2}3^{n+1}(n+1)\left(H_{0}^{2}\right)^{n+1}}{\pi}\geq
0\,,$ (60)

$SEC\Leftrightarrow-\frac{2^{n-2}3^{n+1}(n+1)\left(H_{0}^{2}\right)^{n+1}(b(2n(q_{0}+1)+q_{0}-1)+8\pi(q_{0}-1))}{\pi(b+8\pi)}\geq
0\,.$ (61)

Figure 2: Density parameter, and Energy conditions for $f(Q,T)=Q^{n+1}+bT$.
The graphics were depicted with the present values of $H_{0}$ and $q_{0}$
parameters.

Analogously to our first case, we depicted the density parameter, as well as
the WEC, DEC, and SEC energy conditions, whose features can be appreciated in
Fig. 2. There we observe a slow decreasing of $\rho$ in respect of parameters
$n$, and $b$ corresponding to an expansion regime smoother than our first
case, such a behavior corroborates with features analyzed by Xu et al.
Yixin/2019 . Moreover, Eqs.(58), (59), (60), and (61), unveil that NEC and DEC
are satisfied, while WEC is partially obeyed ($\rho>0$) if $n\leq-1$, and
$b>-4\,\pi\,$. Yet in the energy conditions, we can also see that SEC is again
violated, confirming that our Universe experiences an accelerated phase. If we
consider $b=0$ the model reduces the $f(Q,T)$ gravity to $f(Q)$ which presents
a strong non-minimal coupling problem. S. Mandal Mandal/2020 addressed this
model, explaining how different energy conditions with $n>0$ give rise to
accelerated expansion due to SEC violations. Furthermore, the WEC violation
along with positive density, makes this $f(Q,T)$ gravity naturally behaves
like scalar-tensor gravity models Whinnett/2004 ; Mandal/2020 . Despite the
viability of such a theory in respect to the energy conditions, the
constrained values of parameter $n$ are out of the phenomenological bounds
established in Simran/2020 . There, Arora et al. used Hubble parameter and
Supernovae data sets to constraint $n$ and $b$ as $n\,\in\,(1,4)$ and
$b\,\in(0,2)$. Therefore, the constraints here imposed for $n$ and $b$ create
a tension in use such a $f(Q,T)$ model as a proper description of gravity. We
are going to present some extra comments concerning this tension in the next
section.

##  V Comparison with $\Lambda$CDM model

As a matter of completeness, let us compare our constraints with the
$\Lambda$CDM model. This model is so far the most well succeed to describe the
evolution of the Universe at different phases. A direct way to link an
$f(Q,T)$ gravity with the $\Lambda$CDM model consists to take the special case
$f(Q,T)=f_{\Lambda}(Q)=-Q$ Lazkoz/2019 . Such a regime yields to the following
energy conditions

  * •

NEC: $2\,(1+q)\,H^{2}\geq 0$ ,

  * •

WEC: $3\,H^{2}\geq 0\,,$ and $2\,(1+q)\,H^{2}\geq 0$ ,

  * •

SEC: $6\,q\,H^{2}\geq 0$ ,

  * •

DEC: $2\,(2-q)\,H^{2}\geq 0$ .

One can observe that all energy conditions are satisfied with the present
values of $H$ and $q$ except SEC, corroborating with the description of an
accelerated expansion. This behavior is compatible with the first model here
analyzed for the $f(Q,T)$ gravity.

Another interesting cosmological parameter which is bounded by experiments is
the equation of state parameter $\omega$. Recent observations from Planck
Collaboration inform that $\omega\simeq-1$ Planck/2018 . Therefore, the EoS
parameter is considered a suitable candidate for comparing our models with
$\Lambda$CDM. The EoS parameter ($\omega$) is defined as
$\omega=\frac{p}{\rho}$.
By taking our previous relations for density, and pressure, we are able to
find that the EoS parameter for the model A i.e. $f(Q,T)=mQ+bT$ and model B
i.e. $f(Q,T)=Q^{n+1}+bT$ are respectively written as

$\omega=\frac{4H_{0}^{2}m(3bq_{0}+8\pi(2q_{0}-1))}{b\left(\frac{6(3b+16\pi)H_{0}^{2}m}{b}+4H_{0}^{2}m(q_{0}+1)-6H_{0}^{2}m\right)}.$
(62)

and

$\omega=\frac{3b(nq_{0}+n+q_{0})+8\pi(2n(q_{0}+1)+2q_{0}-1)}{b(nq_{0}+n+q_{0}+4)+24\pi}.$
(63)

Figure 3: EoS parameter for model $f(Q,T)=mQ+bT$ Figure 4: EoS parameter for
$f(Q,T)=Q^{n+1}+bT$

The features of the EoS parameters here derived are presented in Figs 3, and
4. In Fig. 3 we have $\omega\approx-1$ from upper values, showing a negative
pressure phase compatible with the description of $\Lambda$CDM model.
Moreover, in Fig. 4 we observe that $\omega\approx-1$ from lower values,
surprisingly unveiling a behavior compatible with a phantom era for the dark
energy.

As it is known, the $\Lambda$CDM model, where $\Lambda$ represents a strictly
constant vacuum energy, yielding to a negative pressure regime for the
Universe, forbids the existence of a phantom era Cardenas/2020 ; Carroll/2004
; Carroll/2003 . Moreover, a phantom era for the dark energy could lead to a
troublesome description for the Universe when we deal with a standard
description of gravity plus a fundamental phantom field since DEC energy
condition would be violated Carroll/2004 ; Sahni/2000 . Furthermore, in
Einstein’s gravity the equation of state parameter is constrained to be
$\geq-1$, which means that DEC is satisfied Carroll/2003 . However, once the
current experimental bounds for $\omega$ establish that

$\omega=-1.03_{-0.03}^{+0.03}\,,\qquad\mbox{SNe data
\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Planck/2018}{\@@citephrase{(}}{\@@citephrase{)}}}}\,,$
(64)

we still have room for a phantom description of dark energy Wang/2019 . In the
phantom regime the energy density gradually increases, making the Universe
accelerates so fast breaking the particles interactions. Overall, there are
three categories of phantom regime based on the time evolution of the Hubble
parameter, namely Wang/2019 ; Frampton/2012 : Big rip if
$H(t)\rightarrow\infty$ as $t_{rip}$ is constant - in this scenario even the
spacetime rips apart; Little rip if $H(t)\rightarrow\infty$ as
$t_{rip}\rightarrow infty$ \- this type of rip prevent singularities in future
values of time; Pseudo rip $H(t)\rightarrow\mbox{constant}$ as
$t_{rip}\rightarrow\infty$ \- such a rip dissociates bound structures that are
held together by a binding force. By analyzing the cosmological model for
$f(Q,T)=Q^{n+1}+bT$ studied by Yixin et al. Yixin/2019 , we realize that

$H(t)\rightarrow 0\,\qquad\mbox{for}\qquad t_{rip}\rightarrow\infty\,,$ (65)

characterizing a compatible pseudo rip evolution for the constrained values of
parameters $b$, and $n$ here presented.

Withal, despite these issues concerning the phantom related with $\Lambda$CDM
model, here we show that the$f(Q,T)=Q^{n+1}+bT$ gravity satisfies DEC if
$n\leq-1$, and $b>-4\,\pi\,$. Therefore, such a theory of gravity naturally
enables us to describe a phantom era for the dark energy, without the need of
extra dimensions or phantom scalar fields. The phantom phase is also suggested
to reconcile the tension between local and global measurements of the current
Hubble value $H_{0}$Amr/2019 . Hence, we believe that the tension between this
$f(Q,T)$ model and observational data for the Hubble parameter at different
redshifts, lies in the compatibility of this $f(Q,T)$ with a phantom era
description for the dark energy.

##  VI Conclusion

An essential role to establish a consistent theory of gravity is the energy
condition. As new theories of gravity are bubbling in the literature, it is
relevant to put them up to test through constraints over different energy
conditions. In this work, we computed the strong, the weak, the null, and the
dominant energy conditions for two $f(Q,T)$ gravity models. The $f(Q,T)$ is a
promising new theory for gravity based on the combination of the non-metricity
function $Q$ with the trace of the energy-momentum $T$.

The models here considered were proposed by Xu et al. Yixin/2019 , and
constrained by observational data of the Hubble parameter in Simran/2020 .
Firstly we worked with $f(Q,T)=mQ+bT$, where $m$, and $b$ are free parameters.
The energy conditions yield us to constraint these free parameters as
$b>-4\,\pi$, and $m\leq 0$. The previous values result in the violation of
SEC, corroborating with an accelerated phase of expansion for the Universe.
Besides, such a model is suitable to describe the Universe with respect to
energy conditions as $\Lambda$CDM.

As a second case, we worked with $Q^{n+1}+bT$, whose free parameters should be
constrained to $n\leq-1$, and $b>-4\,\pi\,$, to satisfy DEC, and NEC energy
conditions. In this case, WEC energy condition is partially obeyed while SEC
is again violated. The violation of WEC makes this model naturally behaves
like scalar-tensor gravity theories. Moreover, the model is compatible with
the dark energy era once SEC is not satisfied. A surprisingly feature comes
from the equation of state parameter for this model, which describes a phantom
regime for the dark energy, allowing extra acceleration for the expansion of
the Universe without violates DEC.

The results here presented allowed us to verify the viability of different
families of $f(Q,T)$ gravity models, lighting new paths for a complete
description of gravity compatible with the dark energy era, which embeds
effects from the quantum era of the Universe. The constraints for our free
parameters yield to several testable families for $f(Q,T)$ gravity, opening
space even for models compatible with a phantom regime for the dark energy.
Moreover, it would be interesting to investigate carefully the coupling of
$f(Q,T)$ with inflation fields or with dust, looking for possible analytic
models or for cosmological parameters constraints. It would be also
interesting to impose constraints on such theories of gravity with
observational data from low redshifts, such as BAO measurements at $z=0.1-2.5$
which are expected to be performed in the near feature by BINGO Bingo/2019 ,
and CHIME Zhang/2019 telescopes. We hope to report on some of these
investigations in the near future.

###### Acknowledgements.

S. A. acknowledges CSIR, Govt. of India, New Delhi, for awarding Junior
Research Fellowship. JRLS would like to thank CNPq (Grant no. 420479/2018-0),
CAPES, and PRONEX/CNPq/FAPESQ-PB (Grant nos. 165/2018, and 0015/2019) for
financial support. PKS acknowledges CSIR, New Delhi, India for financial
support to carry out the Research project[No. 03(1454)/19/EMR-II
Dt.02/08/2019].

## References

  * (1) A. G. Riess et al., The Astron. J., 116 (1998) 1009.
  * (2) G. Perlmutter et al., ApJ, 517 (1999) 565.
  * (3) Ronald J. Adler, Brendan Casey, and Ovid C. Jacob, Am. J. Phys, 63 (1995) 620.
  * (4) N. Aghanim et al. [Planck Collaboration], Planck 2018 results. VI. Cosmological parameters, astro-ph.CO/1807.06209.
  * (5) E. Baxter et al. [Dark Energy Survey], Dark Energy Survey Year 1 Results, astroph.CO/1802.05257.
  * (6) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett., 116 (2016) 061102; B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett., 119 (2017) 161101; B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett., 123 (2019) 011102.
  * (7) The Event Horizon Telescope Collaboration et al., ApJL, 875 (2019) L1; The Event Horizon Telescope Collaboration et al., ApJL, 875 (2019) L5.
  * (8) J. B. Jiménez, L. Heisenberg, T. Koivisto, Phys. Rev. D, 98 (2018) 044048.
  * (9) R. Lazkoz et al., Phys. Rev. D, 100 (2019) 104027.
  * (10) T. Harko, F.S.N. Lobo, Extensions of $f(R)$ gravity: Curvature- Matter Couplings and Hybrid Metric-Palatini Theory (Cambridge University Press, Cambridge 2018).
  * (11) T. Harko et al., Phys. Rev. D, 84 (2011) 024020.
  * (12) T. Harko et al., Phys. Rev. D, 98 (2018) 084043.
  * (13) T. Harko et al., arXiv:1901.00805.
  * (14) A. Delhom, Eur. Phys. J. C, 80 (2020) 728.
  * (15) Thomas P Sotiriou, Class. Quantum Grav., 26 (2009) 152001.
  * (16) Sanjay Mandal, P.K. Sahoo, and J.R.L. Santos, Phys. Rev. D, 102 (2020) 024057.
  * (17) Y. Xu et al., Eur. Phys. J. C, 79 (2019) 708.
  * (18) Simran Arora, et al., Physics of the Dark Universe, 30 (2020) 100664.
  * (19) Snehasish Bhattacharjee, P. K. Sahoo, Eur. Phys. J. C, 80 (2020) 289.
  * (20) S. Capozziello et al., Phys. Lett. B, 781 (2018) 99.
  * (21) J. B. Jiménez et al., Phys. Rev. D, 101 (2020) 103507.
  * (22) Kaczmarek, A.Z., Szczesniak, Sci Rep, 10 (2020) 18076.
  * (23) Yixin Xu, Tiberiu Harko, Shahab Shahidi, Shi-Dong Liang, Eur. Phys. J. C, 80 (2020) 449.
  * (24) Jin-Zhao Yang, Shahab Shahidi, Tiberiu Harko, Shi-Dong Liang, arXiv:2101.09956 [gr-qc].
  * (25) Damianos Iosifidis, Christos G. Tsagas, Anastasios C. Petkou, Phys. Rev. D, 98 (2018) 104037.
  * (26) J. Santos et al., Phys. Rev. D, 76 (2007) 083513.
  * (27) V. A. Rubakov, Physics-Uspekhi, 57 (2014) 128.
  * (28) Frank J. Tipler, Phys. Rev. D, 17 (1978) 25218.
  * (29) Crislane S. Santos et al., Gen. Relativ. Gravit., 49 (2017) 50.
  * (30) M. Sharif et al., Eur. Phys. J. C, 76 (2016) 640.
  * (31) S. Capozziello et al., Int. J. Mod. Phys. D, 28 (2019) 1930016.
  * (32) M. Visser, C. Barcelo, COSMO-99, (2000) 98.
  * (33) A. W. Whinnett, D. F. Torres, The Astrophys. J, 603 (2004) L133.
  * (34) Víctor H. Cárdenas, et al., Phys. Rev. D, 101, 083530 (2020).
  * (35) Sean M. Carroll, Antonio De Felice, and Mark Trodden, Phys.Rev.D, 71 (2005) 023525.
  * (36) S. M. Carroll, M. Hoffman, and M. Trodden, Phys. Rev. D, 68 (2003) 023509.
  * (37) V. Sahni, A.A. Starobinsky, Int. J. Mod. Phys. D, 9 (2000) 373.
  * (38) Amr El-Zant et al., The Astrophys. J, 871 (2019) 210.
  * (39) Jun-Chao Wang, and Xin-He Meng, Eur. Phys. J. C, 79 (2019) 848.
  * (40) P. H. Frampton, K. J. Ludwick, R. J. Scherrer, Phys. Rev. D, 85 (2012) 083001.
  * (41) M. W. Peel, et al. J. Astron. Instrum., 08 (2019) 1940005.
  * (42) Jing-Fei Zhang, et al. Phys. Lett B, 799 (2019) 135064.