# Hybrid Parametric Classes of Isotropic Covariance Functions for Spatial
Random Fields

Alfredo Alegría111Corresponding author. Email<EMAIL_ADDRESS>Departamento de Matemática, Universidad Técnica Federico Santa María, Chile
Fabián Ramírez  Departamento de Matemática, Universidad Técnica Federico Santa
María, Chile  Emilio Porcu  Department of Mathematics, Khalifa University, The
Arab Emirates

###### Abstract

Covariance functions are the core of spatial statistics, stochastic processes,
machine learning as well as many other theoretical and applied disciplines.
The properties of the covariance function at small and large distances
determine the geometric attributes of the associated Gaussian random field.
Having covariance functions that allow to specify both local and global
properties is certainly on demand. This paper provides a method to find new
classes of covariance functions having such properties. We term these models
_hybrid_ as they are obtained as scale mixtures of piecewise covariance
kernels against measures that are also defined as piecewise linear combination
of parametric families of measures. In order to illustrate our methodology, we
provide new families of covariance functions that are proved to be richer with
respect to other well known families that have been proposed by earlier
literature. More precisely, we derive a hybrid Cauchy-Matérn model, which
allows us to index both long memory and mean square differentiability of the
random field, and a hybrid Hole-Effect-Matérn model, which is capable of
attaining negative values (hole effect), while preserving the local attributes
of the traditional Matérn model. Our findings are illustrated through
numerical studies with both simulated and real data.

_Keywords_ : Cauchy model; Gaussian scale mixtures; Hole effect; Long memory;
Matérn model; Mean square differentiability.

##  1 Introduction

Covariance functions are central to many disciplines such as spatial
statistics (Cressie,, 1993; Chilès and Delfiner,, 2012; Hristopulos,, 2020),
stochastic processes (Porcu et al., 2018a, ; Porcu et al., 2018b, ), machine
learning (Schaback and Wendland,, 2006; James et al.,, 2013; Barp et al.,,
2022), numerical analysis (Pazouki and Schaback,, 2011; Cockayne et al.,,
2019) and stochastic mechanics (Ostoja-Starzewski,, 2006, with the references
therein). Recent applications in climatology (Guinness and Hammerling,, 2018;
Edwards et al.,, 2019), oceanography (Furrer et al.,, 2007; Di Lorenzo et
al.,, 2014), environmental sciences (Cressie and Kornak,, 2003; Stein,, 2007)
and natural resources engineering (Chen et al.,, 2018; Emery and Séguret,,
2020) witness on the importance of covariance functions.

It is very customary to assume the covariance function to depend on the
distance between any pair of random variables located at two different points
at the input space. Such an assumption is termed isotropy in spatial
statistics and machine learning, and it is termed radial symmetry in other
areas of applied mathematics. The behaviour of the covariance function at
short or long distances (we call this local and global properties,
respectively) is crucial to understand the properties of random processes with
a given covariance function. Specifically, the local properties are related to
the fractal dimension as well as the geometric properties (e.g., mean square
differentiability) of the associated random process, as well as to its sample
paths. On the other hand, the global behaviour of the covariance function
allows to characterize persistency or antipersistency (i.e., the long term
behaviour) of the associated process. Another global behaviour of great
interest is the so-called hole effect, which means that the covariance
function could take negative values in a certain interval.

Finding parametric families of isotropic covariance functions that allow to
index both local and global behaviour is a major challenge that has been
tackled to a very limited extent. The Matérn family has been the cornerstone
in spatial statistics for over half a century now (Stein,, 1999). Its
popularity is due to a parameter that controls the degree of mean square
differentiability and fractal dimension of the corresponding random field
(Stein,, 1999). Recently, Bevilacqua et al., (2022) have shown that the Matérn
class is a special case of a richer class of models that, additionally to
indexing local properties, allow to switch between compact or global supports.
In turn, compactly supported models lead to sparse covariance matrices (Furrer
et al.,, 2006; Kaufman et al.,, 2008) and this implies considerable
computational gains in both estimation and prediction. Unfortunately, the
Matérn class does not allow to index global behaviour of the associated random
process. The Generalized Cauchy family (Gneiting and Schlather,, 2004) allows
to index the fractal dimension and the long memory behaviour. Notably, it does
not allow to index mean square differentiability, as the model is either non
differentiable or infinitely differentiable at the origin. The same properties
are shared by the Dagum model (Berg et al.,, 2008), which does not allow to
index mean square differentiability either. None of the aforementioned models
allow to attain negative spatial dependencies.

Spectral approaches can be a promising avenue to find flexible families of
covariance functions. Laga and Kleiber, (2017) proposed a modified version of
the spectral density associated with the Matérn family. The new class has two
additional parameters that can be loosely interpreted as a continuous version
of a moving average process. More recently, Ma and Bhadra, (2022) have proved
that a two-fold application of Gaussian scale mixtures can provide models with
polynomial decays while preserving the local properties of the candidate
covariance function. Other non conventional properties of covariance functions
have been studied by Alegría, (2020) and Alegría et al., (2021), who proposed
some modified scale mixtures representations to obtain classes of cross-
covariance functions with non-monotonic behaviours (the so-called cross-dimple
effect) for vector-valued random fields. In Schlather and Moreva, (2017),
models that allow for a smooth transition between stationary and intrinsically
stationary Gaussian random fields are derived.

All the previously mentioned parametric classes of covariance functions admit
a scale mixture representation of a Gaussian kernel against a continuous,
positive and bounded measure. Our paper starts from the Schoenberg integral
representation of isotropic covariance functions on $\mathbb{R}^{d}$
(Schoenberg,, 1938), for all natural numbers, $d$. We specifically assume the
Schoenberg measures to be parametric families of measures that are defined
piecewise. Such a strategy is then shown to provide hybrid classes that
generalize classes proposed in earlier literature. We illustrate this
methodology by constructing a model that combines the global attributes of the
Cauchy class and the local properties of the Matérn class. We show that the
proposed model admits a closed form expression and examine its theoretical
properties. Additionally, we study a more flexible formulation, where the
Gaussian kernel involved in the scale mixture is replaced with a covariance
kernel that is also defined piecewise. Following this approach, we derive a
hybrid model with local behaviour of Matérn type, and global behaviour that
allows for covariance functions with negative values. We conduct numerical
experiments with both simulated and real data in order to assess the
statistical performance of the proposed models.

The article is organized as follows. Section 2 contains a concise review of
random fields and covariance functions coming from scale mixtures. Section 3
presents general methodologies to build hybrid covariance models. Then, we
derive the hybrid Cauchy-Matérn and the hybrid Hole-Effect-Matérn classes.
Section 4 guides the reader through some numerical studies. We finally provide
a critical discussion in Section 5, including a description of technical
extensions of the present work such as the multivariate case, where covariance
functions are matrix-valued, and the case of spherically indexed fields, where
isotropy is defined in terms of the geodesic distance.

##  2 Background

Let $\\{Z({\bf s}):{\bf s}\in\mathbb{R}^{d}\\}$ be a (centered) second-order
stationary Gaussian random field on $\mathbb{R}^{d}$. Such a field is
completely characterized by its covariance function (or kernel). The isotropy
of the covariance function is defined through a mapping
$\varphi:[0,\infty)\rightarrow\mathbb{R}$ such that $\text{cov}[Z({\bf
s}),Z({\bf s}^{\prime})]=\varphi(h),$ for every ${\bf s},{\bf
s}^{\prime}\in\mathbb{R}^{d}$, where $h=\|{\bf s}-{\bf s}^{\prime}\|$. The
covariance function must satisfy the positive (semi) definiteness condition:
for any $k\in\mathbb{N}$, $\\{a_{1},\ldots,a_{k}\\}\subset\mathbb{R}$ and
$\\{{\bf s}_{1},\ldots,{\bf s}_{k}\\}\subset\mathbb{R}^{d}$,

$\sum_{i,j=1}^{k}a_{i}a_{j}\varphi(\|{\bf s}_{i}-{\bf s}_{j}\|)\geq 0.$

We use the notation $\varphi(\cdot;\bm{\lambda})$ for a parametric family of
continuous covariance functions, where $\bm{\lambda}\in\mathbb{R}^{p}$ is a
vector of parameters. Further, we make use of the celebrated Schoenberg’
theorem (Schoenberg,, 1938): the functions $\varphi$ that are valid in any
dimension $d\in\mathbb{N}$ are uniquely written as Gaussian scale mixtures of
positive and bounded measures, that is

$\varphi(h;\bm{\lambda})=\int_{0}^{\infty}\exp(-uh^{2})G({\rm
d}u;\bm{\lambda}),\qquad h\geq 0,$

where $\\{G({\rm d}\cdot;\bm{\lambda}),\;\bm{\lambda}\in\mathbb{R}^{p}\\}$ is
a parametric family of measures, that are termed Schoenberg measures in Daley
and Porcu, (2014). Most of the covariance classes listed in the introduction
admit such a representation against a measure that is absolutely continuous
with respect to the Lebesgue measure, that is

$\varphi(h;\bm{\lambda})=\int_{0}^{\infty}\exp(-uh^{2})g(u;\bm{\lambda})\text{d}u,\qquad
h\geq 0,$ (2.1)

for $\\{g(\cdot;\bm{\lambda}),\;\bm{\lambda}\in\mathbb{R}^{p}\\}$ a parametric
family of nonnegative functions. Throughout, we call $g$ the mixing function.

We now describe examples of some parametric classes of functions $\varphi$
that are determined according to (2.1). Special attention is devoted to the
Matérn, Cauchy and Generalized Cauchy models. Other examples, including the
stable and generalized hyperbolic models, can be found in Yaglom, (1987),
Barndorff-Nielsen, (1978), Schlather, (2010) and Porcu et al., 2018b .

######  Example 2.1 (Matérn).

This class of covariance functions is defined as (Matérn,, 1986)

$\varphi_{\mathcal{M}}\left(h;\bm{\lambda}\right)=\frac{2^{1-\nu}}{\Gamma(\nu)}(h/\alpha)^{\nu}K_{\nu}(h/\alpha),\qquad
h\geq 0,$ (2.2)

where $\Gamma$ is the gamma function and $K_{\nu}$ is the modified Bessel
function of the second kind (Abramowitz and Stegun,, 1972). Here,
$\bm{\lambda}=[\alpha,\nu]^{\top}$, with $\alpha$ and $\nu$ being positive
parameters that control the scale (the rate of decay of the covariance in
terms of $h$) and shape of (2.2), respectively. More precisely, $\nu$
regulates the degree of mean square differentiability of the random field
(large values of $\nu$ are associated with smoother sample paths) (Stein,,
1999). When $\bm{\lambda}=[\alpha,1/2]^{\top}$, (2.2) simplifies into the
exponential model, $\exp(-h/\alpha)$. On the other hand, as
$\nu\rightarrow\infty$, a reparameterization of (2.2) tends to the Gaussian
covariance function, defined as $\exp(-h^{2}/\alpha)$.

######  Example 2.2 (Cauchy).

This class of covariance functions is given by (Chilès and Delfiner,, 2012)

$\varphi_{\mathcal{C}}\left(h;\bm{\lambda}\right)=\left(1+h^{2}/\alpha\right)^{-\nu/2},\qquad
h\geq 0,$ (2.3)

with $\bm{\lambda}=[\alpha,\nu]^{\top}$. As in the Matérn model, $\alpha>0$ is
a scale parameter. However, unlike the Matérn model which decays exponentially
with distance (Stein,, 1999), (2.3) has a polynomial decay regulated by
$\nu>0$. When $\nu\in(0,2)$, such a polynomial decay is connected with the
Hurst parameter, a measure of long term memory, given by $H=1-\nu/2$.

######  Example 2.3 (Generalized Cauchy).

This class of covariance functions is defined as (Gneiting and Schlather,,
2004 and the references therein)

$\varphi_{\mathcal{GC}}\left(h;\bm{\lambda}\right)=\left(1+h^{\delta}/\alpha\right)^{-\nu/\delta},\qquad
h\geq 0,$ (2.4)

with $\bm{\lambda}=[\alpha,\nu,\delta]^{\top}$, where $\delta\in(0,2]$,
$\alpha>0$ and $\nu>0$. This generalized class preserves the polynomial decay
of (2.3), but it is more flexible in the sense that the fractal dimension can
be arbitrarily regulated through $\delta$ (see Gneiting and Schlather,, 2004
for details). Maybe surprisingly, this model does not allow to control the
mean square differentiability of the respective random field, as the model is
either non differentiable or infinitely differentiable at the origin.

Additional classes of covariance functions can be obtained from the more
general mixture

$\varphi(h;\bm{\lambda},\bm{\vartheta})=\int_{0}^{\infty}\phi(h;u,\bm{\vartheta})g(u;\bm{\lambda})\text{d}u,\qquad
h\geq 0,$ (2.5)

where $\phi(\cdot;u,\bm{\vartheta})$ is an arbitrary covariance kernel, for
every $u>0$, and $\bm{\vartheta}$ is a vector of parameters. Since the class
of positive definite functions is a convex cone that is closed under the
topology of pointwise convergence, if $\phi$ is valid (positive definite) in
$\mathbb{R}^{d}$ for $d\leq d^{\prime}$, for some $d^{\prime}\in\mathbb{N}$,
then $\varphi$ is valid in $\mathbb{R}^{d}$ for $d\leq d^{\prime}$ as well. We
refer the reader to Emery and Lantuéjoul, (2006) for several explicit
examples.

##  3 Hybrid Classes of Covariance Functions

###  3.1 General Construction

In this study, we propose new parametric classes of isotropic covariance
functions, $\widetilde{\varphi}(\cdot;\bm{\lambda},\bm{\omega},\bm{\xi})$,
determined according to

$\widetilde{\varphi}(h;\bm{\lambda},\bm{\omega},\bm{\xi})=\omega_{1}\int_{0}^{\xi_{1}}\exp(-uh^{2}){g}_{1}(u;\bm{\lambda}_{1})\text{d}u+\omega_{2}\int_{\xi_{2}}^{\infty}\exp(-uh^{2}){g}_{2}(u;\bm{\lambda}_{2})\text{d}u,$
(3.1)

where $g_{1}$ and $g_{2}$ are nonnegative functions on $[0,\xi_{1})$ and
$[\xi_{2},\infty)$, respectively, and
$\bm{\lambda}=[\bm{\lambda}_{1}^{\top},\bm{\lambda}_{2}^{\top}]^{\top}$,
$\bm{\omega}=[\omega_{1},\omega_{2}]^{\top}$ and
$\bm{\xi}=[\xi_{1},\xi_{2}]^{\top}$ are vectors of parameters, with
$\omega_{i},\xi_{i}>0$, for $i=1,2$. In other words, we replace the mixing
function, $g$, in Equation (2.1) with a function $\widetilde{g}$ that is
defined piecewise, i.e.,

$\widetilde{g}(u;\bm{\lambda},\bm{\omega},\bm{\xi})=\omega_{1}\,g_{1}(u;\bm{\lambda}_{1})1_{[0,\xi_{1})}(u)+\omega_{2}\,g_{2}(u;\bm{\lambda}_{2})1_{[\xi_{2},\infty)}(u),\qquad
u\geq 0,$ (3.2)

with $1_{A}(\cdot)$ standing for the indicator function of a set $A$. Note
that $\widetilde{g}$ may have discontinuities as it is built by gluing two
individual pieces. If the functions $g_{i}$ are continuous and bounded on
their domains, a direct application of the dominated convergence theorem
implies that the proposed covariance function (3.1) is continuous on
$[0,\infty)$. Throughout this manuscript, each function $g_{i}$ is positively
proportional to a continuous probability density function. Hence, the
parametric family proposed in Equation (3.1) belongs to the Schoenberg class
as defined through Equation (2.1).

A more general construction considers different kernels in each segment of the
mixture, i.e.,

$\widetilde{\varphi}(h;\bm{\lambda},\bm{\omega},\bm{\xi},\bm{\vartheta})=\omega_{1}\int_{0}^{\xi_{1}}\phi_{1}(h;u,\bm{\vartheta}_{1}){g}_{1}(u;\bm{\lambda}_{1})\text{d}u+\omega_{2}\int_{\xi_{2}}^{\infty}\phi_{2}(h;u,\bm{\vartheta}_{2}){g}_{2}(u;\bm{\lambda}_{2})\text{d}u,$
(3.3)

where
$\bm{\vartheta}=[\bm{\vartheta}^{\top}_{1},\bm{\vartheta}^{\top}_{2}]^{\top}$.
If $\phi_{i}$ is a valid covariance function in $\mathbb{R}^{d}$ for $d\leq
d^{\prime}_{i}$, for some $d_{i}^{\prime}\in\mathbb{N}$, $i=1,2$, then (3.3)
is a valid model in $\mathbb{R}^{d}$ if and only if
$d\leq\text{min}(d^{\prime}_{1},d^{\prime}_{2})$. The continuity of (3.3) can
be justified by following the same arguments used for the continuity of (3.1).

######  Remark 3.1.

Let us point out some additional remarks on this methodology.

  1. 1.

When $\xi_{1}=\xi_{2}=\xi$, this parameter produces a continuous bridge
between two apparently disunited _marginal_ models. More precisely, as it
increases from $0$ to $\infty$, we gradually go from
$\omega_{2}\int_{0}^{\infty}\phi_{2}(h;u,\bm{\vartheta}_{2}){g}_{2}(u;\bm{\lambda}_{2})\text{d}u$
to
$\omega_{1}\int_{0}^{\infty}\phi_{1}(h;u,\bm{\vartheta}_{1}){g}_{1}(u;\bm{\lambda}_{1})\text{d}u$.

  2. 2.

When $\xi_{1}>\xi_{2}$, instead, there is a superposition of the marginal
structures in the interval $[\xi_{2},\xi_{1})$. As $\xi_{2}\rightarrow 0$ and
$\xi_{1}\rightarrow\infty$, we obtain the greatest possible superposition,
which corresponds to a linear combination of the marginal models,
$\omega_{1}\int_{0}^{\infty}\phi_{1}(h;u,\bm{\vartheta}_{1}){g}_{1}(u;\bm{\lambda}_{1})\text{d}u+\omega_{2}\int_{0}^{\infty}\phi_{2}(h;u,\bm{\vartheta}_{2}){g}_{2}(u;\bm{\lambda}_{2})\text{d}u.$

The apparent flexibility of the proposed mixtures is justified by classical
theory on local and global behaviour of covariance functions. In particular, a
direct application of Tauberian theorems (Stein,, 1999) proves that mean
square differentiability of $\widetilde{\varphi}$ will be determined by
$g_{2}$. On the other hand, direct inspection in concert with Equation (4) in
Gneiting and Schlather, (2004) shows that the long term behaviour of
$\widetilde{\varphi}$ is decided by $g_{1}$. The next sections show that it is
possible to provide examples in algebraically closed form that allow to attain
the desired flexibility.

###  3.2 A Hybrid Cauchy-Matérn Class

We present a hybrid Cauchy-Matérn model, for which the acronym ${\cal CM}$ is
used. This model is a special case of (3.1). Let us first introduce the
generalized incomplete gamma function (Chaudhry and Zubair,, 1994),

$\Gamma(a;b;c)=\int_{b}^{\infty}t^{a-1}\exp(-t-ct^{-1})\,\text{d}t,$

and the lower incomplete gamma function,
$\gamma(a,b)=\Gamma(a;0;0)-\Gamma(a;b;0)$.

######  Proposition 3.1.

Let $\bm{\lambda}=[\bm{\lambda}_{1}^{\top},\bm{\lambda}_{2}^{\top}]^{\top}$,
with $\bm{\lambda}_{i}=[\alpha_{i},\nu_{i}]^{\top}$,
$\bm{\omega}=[\omega_{1},\omega_{2}]^{\top}$ and
$\bm{\xi}=[{\xi}_{1},\xi_{2}]^{\top}$ be vectors having positive elements. Let

$\widetilde{\varphi}_{\mathcal{CM}}(h;\bm{\lambda},\bm{\omega},\bm{\xi})=\omega_{1}\,\widetilde{\varphi}^{\,(1)}_{\mathcal{C}}(h;\bm{\lambda}_{1},\xi_{1})+\omega_{2}\,\widetilde{\varphi}^{\,(2)}_{\mathcal{M}}(h;\bm{\lambda}_{2},\xi_{2}),\qquad
h\geq 0,$ (3.1)

where

$\widetilde{\varphi}^{\,(1)}_{\mathcal{C}}(h;\bm{\lambda}_{1},\xi_{1})=\frac{\gamma(\nu_{1}/2,(h^{2}+\alpha_{1})\xi_{1})}{\Gamma(\nu_{1}/2)}\varphi_{\mathcal{C}}(h;\bm{\lambda}_{1})$
(3.2)

and

$\widetilde{\varphi}^{\,(2)}_{\mathcal{M}}(h;\bm{\lambda}_{2},{\xi}_{2})=\varphi_{\mathcal{M}}(h;\bm{\lambda}_{2})-\frac{1}{\Gamma(\nu_{2})}\Gamma\left(\nu_{2};\frac{1}{4\xi_{2}\alpha_{2}^{2}};\frac{h^{2}}{4\alpha_{2}^{2}}\right),$
(3.3)

with $\varphi_{{\cal M}}$ and $\varphi_{{\cal C}}$ being, respectively, the
Matérn and the Cauchy models defined at (2.2) and (2.3). Then,
$\widetilde{\varphi}_{\mathcal{CM}}$ is positive definite in $\mathbb{R}^{d}$
for all $d\in\mathbb{N}$.

######  Proof 3.1.

We provide a proof of the constructive type, by showing that
$\widetilde{\varphi}_{{\cal CM}}$ admits the representation (3.1), with
$g_{1}(u;\bm{\lambda}_{1})=g_{{\cal C}}(u;\bm{\lambda}_{1})$ and
$g_{2}(u;\bm{\lambda}_{2})=g_{{\cal M}}(u;\bm{\lambda}_{2})$, with $g_{{\cal
C}}$ and $g_{{\cal M}}$ that are respectively defined as

$g_{\mathcal{C}}(u;\bm{\lambda}_{1})=\frac{\alpha_{1}^{\nu_{1}/2}}{\Gamma(\nu_{1}/2)}u^{\nu_{1}/2-1}\exp\left(-\alpha_{1}u\right),$
(3.4)

and

$g_{\mathcal{M}}(u;\bm{\lambda}_{2})=\frac{1}{\Gamma(\nu_{2})}\left(\frac{1}{2\alpha_{2}}\right)^{2\nu_{2}}u^{-\nu_{2}-1}\exp\left(-\frac{1}{4u\alpha_{2}^{2}}\right),$
(3.5)

where for both cases all the parameters are positive. To attain the analytical
expression of $\widetilde{\varphi}_{{\cal C}}^{(1)}$, we notice that

$\displaystyle\int_{0}^{\xi_{1}}\exp(-uh^{2}){g}_{\mathcal{C}}(u;\bm{\lambda}_{1})\text{d}u$
$\displaystyle=$
$\displaystyle\varphi_{\mathcal{C}}(h;\bm{\lambda}_{1})\int_{0}^{\xi_{1}}\frac{(h^{2}+\alpha_{1})^{\nu_{1}/2}}{\Gamma(\nu_{1}/2)}u^{\nu_{1}/2-1}\exp(-(h^{2}+\alpha_{1})u)\text{d}u$
$\displaystyle=$
$\displaystyle\varphi_{\mathcal{C}}(h;\bm{\lambda}_{1})\frac{\gamma(\nu_{1}/2,(h^{2}+\alpha_{1})\xi_{1})}{\Gamma(\nu_{1}/2)},$

where the second equality is due to the fact that the integral on the right
hand side of the first line amounts to the cumulative distribution function of
a gamma random variable with parameters $h^{2}+\alpha_{1}$ and $\nu_{1}/2$.

To attain the expression of $\widetilde{\varphi}_{{\cal M}}^{(2)}$, we invoke
Equation (10) in Alegría et al., (2021), so that

$\int_{0}^{\xi_{2}}\exp(-uh^{2}){g}_{\mathcal{M}}(u;\bm{\lambda}_{2})\text{d}u=\frac{1}{\Gamma(\nu_{2})}\Gamma\left(\nu_{2};\frac{1}{4\xi_{2}\alpha_{2}^{2}};\frac{h^{2}}{4\alpha_{2}^{2}}\right).$
(3.6)

The function $\widetilde{\varphi}_{{\cal M}}^{(2)}$ is thus attained by
invoking formula 3.471.9 in Gradshteyn and Ryzhik, (2007), for which we have
$\int_{0}^{\infty}\exp(-uh^{2}){g}_{\mathcal{M}}(u;\bm{\lambda}_{2})\text{d}u=\varphi_{\mathcal{M}}(h;\bm{\lambda}_{2}).$
∎

When $\nu_{2}=n+1/2$, for some $n\in\mathbb{N}$, (3.3) can be expressed in
terms of complementary error functions and modified Bessel functions of first
and second kinds. We refer the reader to Alegría et al., (2021) for a more
detailed study of these special cases.

The flexibility of the proposed structure is now illustrated through the
following result, where we use the notation $f_{1}(h)\sim f_{2}(h)$,
$h\rightarrow\infty$, to represent that, for some positive constant $c_{0}$,
the asymptotic relationship $\lim_{h\rightarrow\infty}f_{1}(h)/f_{2}(h)=c_{0}$
holds.

######  Proposition 3.2.

Let $Z$ be a Gaussian random field with covariance function of the form (3.1).
Then, $Z$ is $\kappa$-times mean square differentiable if and only if
$\nu_{2}>\kappa\geq 0$. Moreover, it is true that
$\widetilde{\varphi}_{\mathcal{CM}}(h;\bm{\lambda},\bm{\omega},\bm{\xi})\sim
h^{-\nu_{1}}$, $h\rightarrow\infty.$ Hence, the Hurst parameter associated
with $Z$ is solely indexed by the parameter $\nu_{1}$.

######  Proof 3.2.

Arguments in Chapter 2 of Stein, (1999) show that an isotropic random field
with covariance function $\varphi$ is $\kappa$-times mean square
differentiable if and only if $\varphi^{(2\kappa)}(0;\bm{\lambda})$ exists and
is finite. Direct inspection in concert with dominated convergence on the
invoked Schoenberg representation (2.1) show that this happens if and only if
the mixing function $g$ satisfies

$\int_{0}^{\infty}u^{\kappa}g(u;\bm{\lambda})\text{d}u<\infty.$ (3.7)

We use the latter argument for the special case of the function
$\widetilde{\varphi}_{{\cal CM}}$, for which the tale of the resulting mixing
function is uniquely determined by the mixing function associated with
$\varphi_{{\cal M}}^{(2)}$ as in Proposition 3.1. Direct inspection shows that
(3.7) is true if and only if $\nu_{2}>\kappa$. The first part of the
proposition is established.

For the second part, note that (3.2) behaves as $h^{-\nu_{1}}$, as
$h\rightarrow\infty$, because the lower incomplete gamma function involved in
such an equation tends to $\Gamma(\nu_{1}/2)$, and the Cauchy class with
parameter $\nu_{1}$ decays as $h^{-\nu_{1}}$. The result follows by noting
that (3.3) is dominated by the traditional Matérn model, which decays
exponentially. ∎

To wrap up, the hybrid Cauchy-Matérn model allows to index both mean square
differentiability and long term behaviour of the associated Gaussian random
field. We also note that these properties are independently addressed by the
two parameters $\nu_{1}$ and $\nu_{2}$, and hence those parameters are
statistically identifiable and allow to decouple local and global properties.

From a statistical viewpoint, a parsimonious choice may be considered by
setting $\omega_{1}=\omega_{2}=\omega$, $\alpha_{1}=\alpha_{2}=\alpha$ and
$\xi_{1}=\xi_{2}=\xi$. Thus, we obtain that Proposition 3.1 provides a five
parameter family where $\omega$ indexes the variance, $\alpha$ the scale,
$\nu_{2}$ the mean square differentiability, and $\nu_{1}$ the Hurst effect,
whereas $\xi$ is a parameter that balances the shapes of the marginal
structures involved in this model. Hence, (3.1) generalizes the Matérn model
in that it allows for polynomial decay while indexing continuously mean square
differentiability.

Figure 3.1 shows the parsimonious hybrid Cauchy-Matérn model for different
values of $\xi$. The traditional Matérn and traditional Cauchy, as well as
their average, which are also special cases of the hybrid construction, are
reported for comparison purposes. Note that the curves have a linear or
parabolic decay near the origin according to $\nu_{2}=1/2$ or $\nu_{2}=3/2$,
respectively, and then the decay is more gradual (polynomial rate) for large
distances according to $\nu_{1}$, which is consistent with the local and
global patterns that are coexisting. We observe that $\xi$ has a manifest
impact on the shape of the covariance function, as it produces some
interesting forms (apparent changes of concavity) that could be useful in
practice.

Figure 3.1: Parsimonious hybrid Cauchy-Matérn model for $\omega=1/2$,
$\alpha=1/8$, $\nu_{1}=3/4$ and different values of $\xi$. (Left)
$\nu_{2}=1/2$ and (Right) $\nu_{2}=3/2$. The dashed lines represent the purely
Cauchy, purely Matérn, and their average. All the models have been
appropriately rescaled in order to obtain correlation functions.

###  3.3 A Hybrid Hole-Effect-Matérn Class

We now present a hybrid class of covariance functions, with local attributes
of Matérn type, attaining negative values at large distances. We use the
acronym $\mathcal{HM}$ for this model, termed hybrid Hole-Effect-Matérn. The
proposed class comes from the mixture (3.3), where $\phi_{1}$ is chosen in
such a way that the resulting model can take negative values.

######  Proposition 3.1.

Let $\bm{\lambda}=[\bm{\lambda}_{1}^{\top},\bm{\lambda}_{2}^{\top}]^{\top}$,
with $\bm{\lambda}_{i}=[\alpha_{i},\nu_{i}]^{\top}$,
$\bm{\omega}=[\omega_{1},\omega_{2}]^{\top}$ and
$\bm{\xi}=[{\xi}_{1},\xi_{2}]^{\top}$ be vectors having positive elements, and
$\bm{\vartheta}=[\tau,\eta]^{\top}$ be a vector of additional parameters. Let

$\widetilde{\varphi}_{\mathcal{HM}}(h;\bm{\lambda},\bm{\omega},\bm{\xi},\bm{\vartheta})=\omega_{1}\,\widetilde{\varphi}_{\mathcal{H}}^{\,(1)}(h;\bm{\lambda}_{1},\xi_{1},\bm{\vartheta})+\omega_{2}\,\widetilde{\varphi}_{\mathcal{M}}^{\,(2)}(h;\bm{\lambda}_{2},\xi_{2}),\qquad
h\geq 0,$ (3.1)

where

$\widetilde{\varphi}^{\,(1)}_{\mathcal{H}}(h;\bm{\lambda}_{1},\xi_{1},\bm{\vartheta})=\frac{\tau}{\Gamma(\nu_{1})}\Gamma\left(\nu_{1};\frac{1}{4\xi_{1}\alpha_{1}^{2}};\frac{\eta
h^{2}}{4\alpha_{1}^{2}}\right)-\frac{1}{\Gamma(\nu_{1})}\Gamma\left(\nu_{1};\frac{1}{4\xi_{1}\alpha_{1}^{2}};\frac{h^{2}}{4\alpha_{1}^{2}}\right),$
(3.2)

and $\widetilde{\varphi}^{\,(2)}_{\mathcal{M}}$ as in (3.3). Then,
$\widetilde{\varphi}_{\mathcal{HM}}$ is positive definite in $\mathbb{R}^{d}$
if and only if $1<\eta<\tau^{2/d}$.

######  Proof 3.3.

We consider the construction (3.3), with both $g_{1}$ and $g_{2}$ of the form
(3.5), and $\phi_{2}$ of Gaussian type. Thus, the derivation of
$\widetilde{\varphi}^{\,(2)}_{\mathcal{M}}$ follows the same arguments
employed in the proof of Proposition 3.1.

Before deriving (3.2), let us introduce the following lemma, which is a
combination of Corollaries 4, 8 and 11 in Posa, (2022).

######  Lemma 3.1.

The mapping $h\mapsto A\exp(-ah^{2})-B\exp(-bh^{2})$ is positive definite in
$\mathbb{R}^{d}$ if and only if

$1<\frac{a}{b}<\left(\frac{A}{B}\right)^{2/d}.$ (3.3)

Although Posa, (2022) focused on dimensions $d\leq 3$, the same proof can be
used in arbitrary dimensions. To obtain the expression (3.2), we take the
following covariance kernel in the first segment of the scale mixture

$\phi_{1}(h;u,\bm{\vartheta})=\tau\exp(-u\eta h^{2})-\exp(-uh^{2}),\qquad
h\geq 0.$ (3.4)

Lemma 3.1 ensures that (3.4) is positive definite in $\mathbb{R}^{d}$,
provided that $u>0$ and $1<\eta<\tau^{2/d}$. Thus,

$\widetilde{\varphi}^{\,(1)}_{\mathcal{H}}(h;\bm{\lambda}_{1},\xi_{1},\bm{\vartheta})=\tau\int_{0}^{\xi_{1}}\exp(-u\eta
h^{2})g_{\mathcal{M}}(u;\bm{\lambda}_{1})\text{d}u-\int_{0}^{\xi_{1}}\exp(-uh^{2})g_{\mathcal{M}}(u;\bm{\lambda}_{1})\text{d}u.$
(3.5)

Finally, we invoke the identity (3.6), and we apply it to each integral
involved in the right hand side of Equation (3.5). ∎

The covariance function (3.4) always takes negative values (Posa,, 2022), so
it is a natural building block to achieve hybrid models with hole effect. The
parameters in $\bm{\vartheta}$ are responsible for the sharpness of the hole
effect. More precisely, as $\eta$ approaches $\tau^{2/d}$, the hole effect is
more pronounced because the positive term in the right hand size of (3.4) has
less dominance. Moreover, when $d=1$, we have the least restrictive condition
on $\eta$, and the resulting hole effect is more marked. It is well known that
the possibility of significant negative correlations vanishes as the dimension
increases (see, e.g., page 45 in Stein,, 1999).

The next proposition characterizes the local attributes of (3.1) and provides
a lower bound for this model.

######  Proposition 3.2.

Let $Z$ be a Gaussian random field with covariance function of the form (3.1).
Then, $Z$ is $\kappa$-times mean square differentiable if and only if
$\nu_{2}>\kappa\geq 0$. Moreover, we have the lower bound

$\widetilde{\varphi}_{\mathcal{HM}}(h;\bm{\lambda},\bm{\omega},\bm{\xi},\bm{\vartheta})\geq\omega_{1}(\tau\eta)^{-1/(\eta-1)}\left(\frac{1-\eta}{\eta}\right)\left(1-\frac{\gamma(\nu_{1};\alpha_{1}/\xi_{1})}{\Gamma(\nu_{1})}\right),\qquad
h\geq 0.$ (3.6)

######  Proof 3.4.

The fact that $\nu_{2}$ controls the mean square differentiability is a direct
consequence of the arguments used in the proof of Proposition 3.2. On the
other hand, to find a lower bound, we note that

$\displaystyle\widetilde{\varphi}_{\mathcal{HM}}(h;\bm{\lambda},\bm{\omega},\bm{\xi},\bm{\vartheta})$
$\displaystyle\geq$ $\displaystyle\omega_{1}\,\inf_{h\geq
0}\widetilde{\varphi}_{\mathcal{H}}^{\,(1)}(h;\bm{\lambda}_{1},\xi_{1},\bm{\vartheta})+\omega_{2}\,\inf_{h\geq
0}\widetilde{\varphi}_{\mathcal{M}}^{\,(2)}(h;\bm{\lambda}_{2},\xi_{2})$
$\displaystyle=$ $\displaystyle\omega_{1}\int_{0}^{\xi_{1}}\inf_{h\geq
0}\phi_{1}(h;u,\bm{\vartheta})g_{1}(u;\bm{\lambda}_{1})\text{d}u.$

A straightforward calculation shows that $\phi_{1}$ attains its minimum value
at $h^{*}=\sqrt{\frac{\log(\tau\eta)}{u(\eta-1)}}$. Thus,

$\phi_{1}(h;u,\bm{\vartheta})\geq\phi_{1}(h^{*};u,\bm{\vartheta})=\tau\exp\left(-\frac{\eta\log(\tau\eta)}{\eta-1}\right)-\exp\left(-\frac{\log(\tau\eta)}{\eta-1}\right)=(\tau\eta)^{-1/(\eta-1)}\left(\frac{1-\eta}{\eta}\right).$

Since $g_{1}$ is given by (3.5), we invoke the formula of the cumulative
distribution function of an inverse gamma random variable to establish that

$\int_{0}^{\xi_{1}}g_{1}(u;\bm{\lambda}_{1})\text{d}u=1-\frac{\gamma(\nu_{1};\alpha_{1}/\xi_{1})}{\Gamma(\nu_{1})}.$

The proof is completed. ∎

Note that as $\xi_{1}\rightarrow\infty$ (i.e., as the hole effect
predominates), the lower bound in Equation (3.6) decreases to
$(\tau\eta)^{-1/(\eta-1)}(1-\eta)/\eta$. On the contrary, as
$\xi_{1}\rightarrow 0$, such a bound increases to zero, i.e., the hole effect
becomes negligible, which is not surprising, because in such a case the Matérn
class is predominant. A similar conclusion can be obtained in the limit case
$\eta\rightarrow 1$.

A parsimonious variant of this model consists of taking
$\omega_{1}=\omega_{2}=\omega$ (variance parameter),
$\alpha_{1}=\alpha_{2}=\alpha$ (scale parameter) and $\nu_{1}=\nu_{2}=\nu$
(smoothness parameter), whereas $\bm{\vartheta}$ regulates the hole effect (as
discussed above) and $\xi_{1}=\xi_{2}=\xi$ has a similar interpretation as in
the hybrid Cauchy-Matérn model.

Figure 3.2 shows the parsimonious hybrid Hole-Effect-Matérn model for
different values of $\xi$. The limit cases described in Remark 3.1 are also
reported, in a similar fashion to Figure 3.1. It can be seen that negative
values coexist with different levels of smoothness at the origin, as expected.

Figure 3.2: Parsimonious hybrid Hole-Effect-Matérn model in dimension one, for
$\omega=1/2$, $\alpha=1/8$, $\tau=2$, $\eta=7/2$ and different values of
$\xi$. (Left) $\nu=1/2$ and (Right) $\nu=3/2$. The dashed lines represent the
limit cases reported in Remark 3.1. All the models have been appropriately
rescaled in order to obtain correlation functions.

##  4 Numerical Experiments

###  4.1 Simulated Data

We conduct simulation studies to assess the performance of maximum likelihood
inference when a hybrid covariance structure is present. We focus on the
parsimonious hybrid Cauchy-Matérn dependence structure, as it will be applied
to real data in the next section. We consider $\omega=1$, $\alpha=1/8$,
$\nu_{1}=3/4$ and the following scenarios for $[\nu_{2},\xi]$: (a) $[1/2,40]$,
(b) $[1/2,120]$, (c) $[3/2,40]$ and (d) $[3/2,120]$. For each scenario, we
simulate 200 independent realizations of a Gaussian random field on 100
uniformly sampled points in the square $[0,3]^{2}$ and estimate the parameters
through maximum likelihood. We then repeat the experiment with 256 spatial
locations. We only estimate $\omega$ , $\alpha$ and $\xi$, whereas $\nu_{1}$
and $\nu_{2}$ are fixed, which is a common practice in geostatistics. Instead
of directly estimating $\xi$, we consider the following alternative
parameterization: $\widetilde{\xi}=\sqrt{\xi}\alpha$, which seems to be a
natural choice according to Equations (3.2) and (3.3).

Figure 4.1 displays the results. The estimates are approximately unbiased and
the variance decreases as the sample size increases from 100 to 256, which is
an expected behaviour. The variability of the estimates substantially
decreases in scenarios (c) and (d), i.e., when the random field is smoother,
which is a typical attribute of likelihood-based estimates in this context
(Bevilacqua and Gaetan,, 2015). On the contrary, such a variability
deteriorates as $\xi$ increases from $40$ to $120$. Figure 4.2 shows the log-
likelihood in terms of $\xi$ and $\alpha$, with fixed $\omega$, for a single
realization of the random field, under scenario (b). Although the surface has
a clear maximum value, the objective function is apparently more flat in the
direction of $\xi$. This could explain the increased variability in scenarios
(b) and (d), with respect to (a) and (c). Despite the previous remarks, in
general, the estimates appear to be reasonable in each scenario and no
identifiability issues are observed.

Figure 4.1: Centered boxplots of the maximum likelihood estimates for the
parsimonious hybrid Cauchy-Matérn model in scenarios (a)-(d).

Figure 4.2: Log-likelihood function, with respect to $\alpha$ and $\xi$, for
scenario (b). Left and right panels correspond to the same plot from different
viewpoints.

We now explore the predictive performance of the proposed class through a
cross validation analysis. We simulate 200 independent realizations on 100
uniformly sampled locations in $[0,3]^{2}$ according to the scenarios (a)-(d)
described above. We assess the accuracy through a leave-one-out prediction
strategy in terms of the mean squared error (MSE), mean absolute error (MAE),
log-score (LSCORE) and continuous ranked probability score (CRPS) (see Zhang
and Wang,, 2010). Small values of these indicators suggest superior
predictions. We evaluate the performance of the hybrid Cauchy-Matérn model,
using the Generalized Cauchy class as benchmark. Thus, for each realization,
we estimate the parameters with both models and proceed to make the
predictions through a simple kriging approach. The Generalized Cauchy model
(2.4) has been augmented with a multiplicative parameter $\omega$, namely
$h\mapsto\omega(1+h^{\delta}/\alpha)^{-\nu/\delta}$, so it is parameterized by
$\omega$ and ${\alpha}$, and $\nu=3/4$ and $\delta=1,2$ are fixed.

Table 4.1 shows that, in each scenario, the proposed hybrid model outperforms
its competitor. All the cross-validation scores substantially decrease in
scenarios (c) and (d). From this brief study, we observe that when the true
underlying covariance has a hybrid structure, an incorrect specification of
the spatial association has a negative impact on the posterior predictions.
Since the behaviour of an isotropic covariance function near the origin has a
strong impact on the quality of predictions (Stein,, 1999), our simulation
experiment suggests that in some circumstances the local shape of the proposed
model cannot be replicated by other appealing existing structures.

Table 4.1: Cross-validation scores for the parsimonious hybrid Cauchy-Matérn and Generalized Cauchy (with $\delta=1,2$) models in scenarios (a)-(d). Scenario | Model | MSE | MAE | LSCORE | CRPS
---|---|---|---|---|---
(a) | Hybrid Cauchy-Matérn | 0.706 | 0.668 | 1.231 | 1.773
| Generalized Cauchy ($\delta=1$) | 0.714 | 0.672 | 1.238 | 1.787
| Generalized Cauchy ($\delta=2$) | 0.718 | 0.674 | 1.241 | 1.798
(b) | Hybrid Cauchy-Matérn | 0.480 | 0.549 | 1.034 | 1.462
| Generalized Cauchy ($\delta=1$) | 0.489 | 0.555 | 1.046 | 1.480
| Generalized Cauchy ($\delta=2$) | 0.497 | 0.559 | 1.055 | 1.506
(c) | Hybrid Cauchy-Matérn | 0.172 | 0.316 | 0.398 | 0.840
| Generalized Cauchy ($\delta=1$) | 0.176 | 0.319 | 0.446 | 0.851
| Generalized Cauchy ($\delta=2$) | 0.176 | 0.319 | 0.414 | 0.862
(d) | Hybrid Cauchy-Matérn | 0.075 | 0.210 | 0.024 | 0.561
| Generalized Cauchy ($\delta=1$) | 0.077 | 0.213 | 0.052 | 0.575
| Generalized Cauchy ($\delta=2$) | 0.082 | 0.219 | 0.103 | 0.612

###  4.2 A Real Data Illustration

The estimation of recoverable resources is a task of fundamental importance in
modern mining processes. A sound evaluation of such resources is crucial from
an economic viewpoint and is critical for assessing the long-term availability
of mineral resources and its impact on society. We consider a data set from a
lateritic nickel deposit mined by open pit in Colombia, which contains
measurements of the grades of nickel, iron, chrome, alumina, magnesia and
silica.

This study focuses on nickel concentrations that are placed at an elevation of
about $120$ meters, where $199$ irregularly spaced observations are available.
We apply a log-transformation to reduce the skewness, and then the sample mean
is subtracted. The resulting values are approximately Gaussian. The left panel
of Figure 4.4 shows the transformed data set. We fit two covariance models:
the former is the parsimonious hybrid Cauchy-Matérn, parameterized by
$\omega$, ${\alpha}$ and $\widetilde{\xi}$, with fixed $\nu_{1}=1/4$ and
$\nu_{2}=1/2$, and the latter is the Generalized Cauchy, parameterized as in
Section 4.1, with fixed $\nu=1/4$ and $\delta=0.95$. The values of the fixed
parameters have been selected after some experimental trials, taking into
account the local behavior of the sample covariance (see Figure 4.3).

Table 4.2 reports the likelihood estimates, with the corresponding standard
errors, and the Akaike information criterion (AIC). We observe that the hybrid
Cauchy-Matérn model outperforms its competitor in terms of AIC. Figure 4.3
shows that the fitted covariance models seem to be reasonably close to the
sample covariance. The fitted models differ substantially near the origin
(distances less than 3 meters), since the hybrid model decays faster. On the
contrary, for larger distances the hybrid model decays slower, although the
difference between the curves becomes slight for distances greater than $15$
meters.

Table 4.2: Parameter estimates and Akaike Information Criterion (AIC) of fitted covariance models. Standard errors are reported in parentheses. Model | $\omega$ | $\alpha$ | $\widetilde{\xi}$ | AIC
---|---|---|---|---
Hybrid Cauchy-Matérn | $1.055$ | $12.31$ | $0.063$ | $-34.03$
| $(0.2516)$ | $(2.801)$ | $(0.042)$ |
Generalized Cauchy | $0.164$ | $1.959$ | $-$ | $-31.93$
| $(0.040)$ | $(0.849)$ | $-$ |
Figure 4.3: Sample (circles) and modeled (solid lines) covariances of log-
nickel concentrations.

In order to compare the models in terms of predictive performance, we conduct
a cross-validation study, in a similar fashion to the experiments performed
with simulated data. Table 4.3 shows evidence, based on a leave-one-out cross-
validation scheme, that the hybrid model has a better performance for this
specific data set. In percentage terms, the MSE shows an improvement of
approximately $3.4\%$. The largest difference occurs when we compare the
LSCORE’s (about $9\%$ improvement).

We conclude this section with an illustration of a downscaled map of log-
nickel concentrations (see Figure 4.4), using the hybrid Cauchy-Matérn model.
The interpolated spatial map, which is obtained through simple kriging, is
exhibited on a spatial grid of approximately $1$ meter ($7500$ locations).
This kriged surface could be useful in small-scale mining processes, as it is
a crucial step for industrial exploration and to quantify mineral reserves.

Table 4.3: Scores for the leave-one-out cross-validation study of log-nickel concentrations. Model | MSE | MAE | LSCORE | CRPS
---|---|---|---|---
Hybrid Cauchy-Matérn | $0.0428$ | $0.1431$ | $-0.1840$ | $0.4113$
Generalized Cauchy | $0.0443$ | $0.1462$ | $-0.1677$ | $0.4159$

Figure 4.4: Log-nickel concentrations (left), with the kriged surface (middle)
and the corresponding variance (right).

##  5 Conclusions and Perspectives

We introduced a simple formalism to build sophisticated parametric families of
covariance functions. We focused on a combination between the Matérn and
Cauchy models, where local (mean square differentiability) and global (long
memory) properties coexist in a single family. We have also illustrated the
use of our methodology by constructing a model that behaves as the Matérn
class at short distances and attains negative values at large distances.
Simulation studies show that a parsimonious hybrid Cauchy-Matérn model has
statistically identifiable parameters. Also, this model provides improvements
in terms of predictive performance in comparison to existing models, when a
hybrid inherent dependence structure is present. We reach similar conclusions
when we apply this methodology to a mining dataset. While similar numerical
studies could be performed for the hybrid Hole-Effect-Matérn model, we avoid
them for the sake of simplicity and brevity. Additional interesting extensions
of this work can be tackled in future investigations. We now provide two
concrete research lines that could emerge from this work.

### Multivariate Hybrid Covariance Models

In many practical situations, two or more variables are simultaneously
recorded. Thus, our findings can be generalized to the case of multivariate
fields $\\{{\bf Z}({\bf s})=(Z_{1}({\bf s}),\ldots,Z_{p}({\bf
s}))^{\top},\;{\bf s}\in\mathbb{R}^{d}\\}$, having an isotropic matrix-valued
covariance function $\Phi:[0,\infty)\to\mathbb{R}^{p\times p}$, that is, ${\rm
cov}[Z_{i}({\bf s}),Z_{j}({\bf s}^{\prime})]=\Phi_{ij}(h)$, $h\geq 0$, where
$h=\|{\bf s}-{\bf s}^{\prime}\|$ and $i,j=1,\ldots,p$. We propose the hybrid
model

$\widetilde{\Phi}(h;\bm{\lambda},\bm{\omega},\bm{\xi})=\omega_{1}\int_{0}^{\xi_{1}}\exp(-uh^{2}){\bf
G}_{1}(u;\bm{\lambda}_{1})\text{d}u+\omega_{2}\int_{\xi_{2}}^{\infty}\exp(-uh^{2}){\bf
G}_{2}(u;\bm{\lambda}_{2})\text{d}u,$

that generalizes (3.1), where the vectors of parameters $\bm{\lambda}_{i}$
must be chosen in such a way that the $p\times p$ matrices ${\bf
G}_{i}(u;\bm{\lambda}_{i})$ are positive semi-definite for every fixed $u\geq
0$. Hence, a straight application of Proposition 4 in Porcu and Zastavnyi,
(2011) would ensure $\widetilde{\Phi}$ to be positive semi-definite. A
multivariate version of the hybrid Cauchy-Matérn covariance function is a
natural candidate. The works of Gneiting et al., (2010) and Moreva and
Schlather, (2022) are relevant to tackle this challenge. A multivariate
version of the formulation (3.3) could be deduced similarly.

### Hybrid Covariance Models on Spheres

For random fields that are indexed by the $d$-dimensional unit sphere,
$\mathbb{S}^{d}$, which is a useful framework when analyzing global data
($\mathbb{S}^{2}$ is used as an approximation of the Earth), the isotropy
assumption is given by ${\rm cov}[Z({\bf s}),Z({\bf
s}^{\prime})]=\psi(\theta)$, ${\bf s},{\bf s}^{\prime}\in\mathbb{S}^{d}$,
where $\psi:[0,\pi]\rightarrow\mathbb{R}$ is a continuous mapping and
$\theta=\arccos({\bf s}^{\top}{\bf s}^{\prime})\in[0,\pi]$ is the geodesic
distance. Schoenberg’s characterization (Schoenberg,, 1942) establishes that a
parametric isotropic covariance function $\psi(;\bm{\lambda})$ is valid in any
dimension $d$, if and only if, it can be written as
${\psi}(\theta;\bm{\lambda})=\sum_{\ell=0}^{\infty}\beta_{\ell}(\bm{\lambda})(\cos\theta)^{\ell}$,
$\theta\in[0,\pi]$, for some nonnegative and summable parametric sequence
$\\{\beta_{\ell}(\bm{\lambda})\\}_{\ell=0}^{\infty}$. Thus, the hybrid models
can be adapted to the spherical context by considering a modified sequence of
the form

$\widetilde{\beta}_{\ell}(\bm{\lambda},\bm{\omega},\bm{\xi})=\omega_{1}\,\beta_{\ell}^{(1)}(\bm{\lambda}_{1})1_{[0,\lfloor\xi_{1}\rfloor)}(\ell)+\omega_{2}\,\beta_{\ell}^{(2)}(\bm{\lambda}_{2})1_{[\lfloor\xi_{2}\rfloor,\infty)}(\ell),\qquad\ell=0,1,\ldots,$

where $\lfloor\xi_{i}\rfloor\geq 0$, for $i=1,2$, with $\lfloor\cdot\rfloor$
standing for the floor function, and $\beta_{\ell}^{(i)}$ being a nonnegative
and summable sequence. The local properties of spherically indexed random
fields, and their connections with the covariance function, have been studied
in past literature (Bingham,, 1973; Guinness and Fuentes,, 2016). However,
global properties such as long memory are less intuitive in this scenario as
the spatial domain is a compact set. Covariance functions with hole effect,
for low-dimensional spheres, could be obtained by adapting formulation (3.3).

## Acknowledgements

Alfredo Alegría was partially supported by the National Agency for Research
and Development of Chile, through grant ANID/FONDECYT/INICIACIÓN/No. 11190686.
Fabián Ramirez was partially supported by the Dirección de Postgrados y
Programas (DPP) of the Universidad Técnica Federico Santa María. Emilio Porcu
is supported by the Khalifa University of Science and Technology under Award
No. FSU-2021-016.

## References

  * Abramowitz and Stegun, (1972) Abramowitz, M. and Stegun, I. A. (1972).  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Dover Publications.
  * Alegría, (2020) Alegría, A. (2020).  Cross-dimple in the cross-covariance functions of bivariate isotropic random fields on spheres.  Stat, 9(1):e301.
  * Alegría et al., (2021) Alegría, A., Emery, X., and Porcu, E. (2021).  Bivariate Matérn covariances with cross-dimple for modeling coregionalized variables.  Spatial Statistics, 41:100491.
  * Barndorff-Nielsen, (1978) Barndorff-Nielsen, O. (1978).  Hyperbolic distributions and distributions on hyperbolae.  Scandinavian Journal of Statistics, 5(3):151–157.
  * Barp et al., (2022) Barp, A., Oates, C. J., Porcu, E., and Girolami, M. (2022).  A Riemann–Stein kernel method.  Bernoulli, 28(4):2181–2208.
  * Berg et al., (2008) Berg, C., Mateu, J., and Porcu, E. (2008).  The Dagum family of isotropic correlation functions.  Bernoulli, 14(4):1134–1149.
  * Bevilacqua et al., (2022) Bevilacqua, M., Caamaño-Carrillo, C., and Porcu, E. (2022).  Unifying compactly supported and Matérn covariance functions in spatial statistics.  Journal of Multivariate Analysis, 189:104949.
  * Bevilacqua and Gaetan, (2015) Bevilacqua, M. and Gaetan, C. (2015).  Comparing composite likelihood methods based on pairs for spatial Gaussian random fields.  Statistics and Computing, 25(5):877–892.
  * Bingham, (1973) Bingham, N. H. (1973).  Positive Definite Functions on Spheres.  Proc. Cambridge Phil. Soc., 73:145–156.
  * Chaudhry and Zubair, (1994) Chaudhry, M. A. and Zubair, S. M. (1994).  Generalized incomplete gamma functions with applications.  Journal of Computational and Applied Mathematics, 55(1):99–123.
  * Chen et al., (2018) Chen, W., Castruccio, S., Genton, M. G., and Crippa, P. (2018).  Current and future estimates of wind energy potential over Saudi Arabia.  Journal of Geophysical Research: Atmospheres, 123(12):6443–6459.
  * Chilès and Delfiner, (2012) Chilès, J.-P. and Delfiner, P. (2012).  Geostatistics: Modeling Spatial Uncertainty, Second Edition.  John Wiley & Sons.
  * Cockayne et al., (2019) Cockayne, J., Oates, C. J., Sullivan, T. J., and Girolami, M. (2019).  Bayesian probabilistic numerical methods.  SIAM review, 61(4):756–789.
  * Cressie, (1993) Cressie, N. (1993).  Statistics for Spatial Data.  Wiley, New York, revised edition.
  * Cressie and Kornak, (2003) Cressie, N. and Kornak, J. (2003).  Spatial Statistics in the Presence of Location Error with an Application to Remote Sensing of the Environment.  Statistical Science, 18(4):436–456.
  * Daley and Porcu, (2014) Daley, D. and Porcu, E. (2014).  Dimension walks and Schoenberg spectral measures.  Proceedings of the American Mathematical Society, 142(5):1813–1824.
  * Di Lorenzo et al., (2014) Di Lorenzo, E., Combes, V., Keister, J., Strub, P., Andrew, T., Peter, F., Marck, O., Furtado, J. ., Bracco, A., Bograd, S., Peterson, W., Schwing, F., Taguchi, B., Hormázabal, S., and Parada, C. (2014).  Synthesis of Pacific Ocean Climate and ecosystem dynamics.  Oceanography, 26(4):68–81.
  * Edwards et al., (2019) Edwards, M., Castruccio, S., and Hammerling, D. (2019).  A multivariate Global Spatio-Temporal Stochastic Generator for Climate Ensembles.  Journal of Agricultural, Biological and Environmental Sciences.
  * Emery and Lantuéjoul, (2006) Emery, X. and Lantuéjoul, C. (2006).  Tbsim: A computer program for conditional simulation of three-dimensional gaussian random fields via the turning bands method.  Computers & Geosciences, 32(10):1615–1628.
  * Emery and Séguret, (2020) Emery, X. and Séguret, S. A. (2020).  Geostatistics for the Mining Industry: Applications to Porphyry Copper Deposits.  CRC Press.
  * Furrer et al., (2006) Furrer, R., Genton, M. G., and Nychka, D. (2006).  Covariance tapering for interpolation of large spatial datasets.  Journal of Computational and Graphical Statistics, 15(3):502–523.
  * Furrer et al., (2007) Furrer, R., Sain, S. R., Nychka, D., and Meehl, G. A. (2007).  Multivariate bayesian analysis of atmosphere–ocean general circulation models.  Environmental and Ecological Statistics, 14(3):249–266.
  * Gneiting et al., (2010) Gneiting, T., Kleiber, W., and Schlather, M. (2010).  Matérn cross-covariance functions for multivariate random fields.  Journal of the American Statistical Association, 105:1167–1177.
  * Gneiting and Schlather, (2004) Gneiting, T. and Schlather, M. (2004).  Stochastic models that separate fractal dimension and the Hurst effect.  SIAM review, 46(2):269–282.
  * Gradshteyn and Ryzhik, (2007) Gradshteyn, I. and Ryzhik, I. (2007).  Table of Integrals, Series, and Products.  Academic Press, Amsterdam, 7th. edition.
  * Guinness and Fuentes, (2016) Guinness, J. and Fuentes, M. (2016).  Isotropic covariance functions on spheres: Some properties and modeling considerations.  Journal of Multivariate Analysis, 143:143–152.
  * Guinness and Hammerling, (2018) Guinness, J. and Hammerling, D. (2018).  Compression and conditional emulation of climate model output.  Journal of the American Statistical Association, 113(521):56–67.
  * Hristopulos, (2020) Hristopulos, D. T. (2020).  Random fields for spatial data modeling.  Springer.
  * James et al., (2013) James, G., Witten, D., Hastie, T., and Tibshirani, R. (2013).  An introduction to statistical learning, volume 112.  Springer.
  * Kaufman et al., (2008) Kaufman, C. G., Schervish, M. J., and Nychka, D. W. (2008).  Covariance tapering for likelihood-based estimation in large spatial data sets.  Journal of the American Statistical Association, 103(484):1545–1555.
  * Laga and Kleiber, (2017) Laga, I. and Kleiber, W. (2017).  The modified Matérn process.  Stat, 6(1):241–247.
  * Ma and Bhadra, (2022) Ma, P. and Bhadra, A. (2022).  Beyond Matérn: On a Class of Interpretable Confluent Hypergeometric Covariance Functions.  Journal of the American Statistical Association, (in press).
  * Matérn, (1986) Matérn, B. (1986).  Spatial Variation — Stochastic Models and Their Application to Some Problems in Forest Surveys and Other Sampling Investigations.  Springer.
  * Moreva and Schlather, (2022) Moreva, O. and Schlather, M. (2022).  Bivariate covariance functions of Pólya type.  Journal of Multivariate Analysis, in press.
  * Ostoja-Starzewski, (2006) Ostoja-Starzewski, M. (2006).  Material spatial randomness: From statistical to representative volume element.  Probabilistic engineering mechanics, 21(2):112–132.
  * Pazouki and Schaback, (2011) Pazouki, M. and Schaback, R. (2011).  Bases for kernel-based spaces.  Journal of Computational and Applied Mathematics, 236(4):575–588.
  * (37) Porcu, E., Alegria, A., and Furrer, R. (2018a).  Modeling temporally evolving and spatially globally dependent data.  International Statistical Review, 86(2):344–377.
  * (38) Porcu, E., Bevilacqua, M., and Hering, A. S. (2018b).  The Shkarofsky-Gneiting class of covariance models for bivariate Gaussian random fields.  Stat, 7(1):e207.
  * Porcu and Zastavnyi, (2011) Porcu, E. and Zastavnyi, V. (2011).  Characterization theorems for some classes of covariance functions associated to vector valued random fields.  Journal of Multivariate Analysis, 102(9):1293–1301.
  * Posa, (2022) Posa, D. (2022).  Special classes of isotropic covariance functions.  Preprint, Research Square.
  * Schaback and Wendland, (2006) Schaback, R. and Wendland, H. (2006).  Kernel techniques: from machine learning to meshless methods.  Acta numerica, 15:543–639.
  * Schlather, (2010) Schlather, M. (2010).  Some covariance models based on normal scale mixtures.  Bernoulli, 16(3):780–797.
  * Schlather and Moreva, (2017) Schlather, M. and Moreva, O. (2017).  A parametric model bridging between bounded and unbounded variograms.  Stat, 6(1):47–52.
  * Schoenberg, (1938) Schoenberg, I. J. (1938).  Metric Spaces and Completely Monotone Functions.  Annals of Mathematics, 39(4):811–841.
  * Schoenberg, (1942) Schoenberg, I. J. (1942).  Positive definite functions on spheres.  Duke Mathematical Journal, 9(1):96–108.
  * Stein, (1999) Stein, M. L. (1999).  Statistical Interpolation of Spatial Data: Some Theory for Kriging.  Springer, New York.
  * Stein, (2007) Stein, M. L. (2007).  Spatial variation of total column ozone on a global scale.  The Annals of Applied Statistics, 1(1):191–210.
  * Yaglom, (1987) Yaglom, A. M. (1987).  Correlation Theory of Stationary and Related Random Functions, Volume I: Basic Results, volume 131.  Springer.
  * Zhang and Wang, (2010) Zhang, H. and Wang, Y. (2010).  Kriging and cross-validation for massive spatial data.  Environmetrics, 21(3-4):290–304.


#  The constrained cosmological model in Lyra geometry

J. K. Singh<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas
University of Technology, New Delhi-110078, India  Shaily
<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas
University of Technology, New Delhi-110078, India  Shri Ram
<EMAIL_ADDRESS>Department of Mathematical Sciences, Indian
Institute of Technology (BHU), Varanasi 221 005, India  Joao R. L. Santos
<EMAIL_ADDRESS>UFCG-Universidade Federal de Campina Grande-Unidade
Academica de Fisica, 58429-900 Campina Grande, PB, Brazil  Jéferson A. S.
Fortunato<EMAIL_ADDRESS>PPGCosmo, CCE, Universidade Federal
do Espírito Santo (UFES), 29.075-910, Vitória, ES, Brazil

###### Abstract

In this article, we study a flat homogeneous FLRW model in Lyra geometry which
is described by a time-dependent displacement vector. We consider an
appropriate parametrization of the energy density of scalar field
$\rho_{\phi}$ in terms of the cosmic scale factor. The result shows two
transitions from deceleration to acceleration. Furthermore, we constrain the
model parameter $\alpha$ and the displacement field vector $\beta$ using the
recent supernovae data, Hubble data set of 77 points and their joint data
which predicts the accelerated expanding phase of the universe in late times.
The effective equation of state parameter $\omega_{eff}$ speculate
$\Lambda$CDM in late times. Finally, we use the statefinder diagnostic to
differentiate our model from the various dark energy models.

PACS numbers: 04.20.-q, 04.50.Kd, 98.80.Es
Keywords: Lyra geometry, Observational constraints, Dark energy, FLRW space-
time

##  I Introduction

From various observations, the accelerating behavior of the Universe is
confirmed. However, figuring out the reason for acceleration is still a big
challenge in cosmology. In standard cosmology, dark energy models are
supported as a resolution to this fundamental question. Here, the cosmological
constant $\Lambda$, which can be associated with the vacuum quantum energy, is
most prominent to describe the dark energy [1, 2]. Another reassuring way to
understand the recent probes for cosmic expansion is Einstein theory of
general relativity [3, 4]. The Einstein’s GR models seize up at large cosmic
scales and a more generic action characterizes the gravitational field. To
generalize the Einstein-Hilbert action of GR, there are various models where
standard action is replaced by the generic function $f(R)$, where $R$ is the
Ricci scalar [5, 6, 7, 8, 9, 10, 11, 12]. Using these $f(R)$ models, the early
inflation of the universe and the concept of dark energy are discussed. In the
extension of $f(R)$ gravity theory, an explicit coupling of matter Lagrangian
density $L_{m}$ is included with $f(R)$ [13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23]. For a proper description of the early universe, homogeneous and
isotropic cosmological models play a significant role.

It is noteworthy that among all of the alternative theories of gravitation,
Lyra geometry is one of the important theory, which is proposed by Lyra [24].
Lyra geometry comes with a modification of Riemannian geometry and may be
considered as a modification of Weyl’s geometry. In Weyl’s geometry, the
change in the length of vectors under parallel transport is non-zero and
depends on a new vector quantity $\zeta_{i}$, which is used for
electromagnetic potential. In 1951, Lyra proposed a different way to maintain
the vector length integrability by adopting a gauge function $\beta(t)$ as an
intrinsic part of the manifold’s geometric structure [25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36]. Further, Sen [37] formulate a static cosmological
model where he obtained the EFE through the variational principle. Singh and
Shri Ram [38] discussed the spatially homogeneous Bianchi type-I metric in the
normal gauge for Lyra’s geometry.

In the literature, a reconstructing technique or model-independent technique
has been employed. Nowadays, this model-independent technique is of great
interest to configure of some DE candidates. Initially, this approach is
discussed by Starobinsky [16]. Many authors proposed this parametrization
scheme on density, pressure, deceleration parameter, Hubble parameter, and
scale factor. This method has two categories: (1) non-parametric and (2)
parametric. The non-parametric method has certain restrictions because it
could not be able to explain all the incidents of the present Universe while
the parametric method includes a specific parametrization of the distinct
types of cosmological parameters which are used to explore the cosmological
models [39].

The present work is organized as follows: In Sec. II, we present the model
field equations using Lyra’s geometry frame transformation for flat FLRW
space-time metric. We parametrize the energy density for the scalar field and
obtain the analytical solutions of the Hubble parameter and other cosmological
parameters. In Sec. III, we discuss a brief introduction of Type Ia supernova
data and observational Hubble data to derive the constraints on the model
parameters. In Sec. IV, we study the evolution of various quantities of the
universe. Finally, we discuss and conclude our findings in Sec. V.

##  II Field Equations in Lyra geometry

The gravitational action proposed by Sen [37] is defined as

$A=\int(\beta^{4}R+L_{\phi}+L_{m})\sqrt{-g}d^{4}x,$ (1)

Using Lyra’s reference frame transformation the field equations obtained as:

$R_{\gamma\delta}-\frac{1}{2}Rg_{\gamma\delta}+\frac{3}{2}\zeta_{\gamma}\zeta_{\delta}-\frac{3}{4}g_{\gamma\delta}\zeta_{\lambda}\zeta^{\lambda}=\kappa
T_{\gamma\delta},$ (2)

where $\kappa=8\pi G$ which is normalized to $1$ for further calculation,
$R_{\gamma\delta}$ is Ricci tensor, $R$ is scalar curvature, and
$\zeta_{\gamma}$ is the time-like displacement vector, defined as
$(\beta(t),0,0,0)$.

The flat FLRW space-time metric is given by

$ds^{2}=-a^{2}(t)(dx^{2}+dy^{2}+dz^{2})+dt^{2},$ (3)

where $a(t)$ is the scale factor of the universe. From Eq. (2), we get

$3H^{2}-\frac{3\beta^{2}}{4}=\rho_{eff}=\rho_{\phi}+\rho_{m},$ (4)
$2\dot{H}+3H^{2}+\frac{3\beta^{2}}{4}=-p_{eff}=p_{\phi}.$ (5)

Overhead dot is used for derivative w.r.t. time. $\rho_{eff}$ and $p_{eff}$
signifies the total amount of energy density and pressure contained in the
universe. Also $\rho_{\phi}$, $\rho_{m}$, $p_{\phi}$ depict the energy density
of the scalar field, energy density of matter, and pressure of scalar field
respectively. According to the present Universe is filled with two types of
fluid, one is corresponding to the scalar field and another is pressureless
cold matter.

The action of the scalar field is defined as:

$L_{\phi}=\int[\frac{1}{2}\partial_{\nu}{\phi}\partial^{\nu}{\phi}-V(\phi)]\sqrt{-g}d^{4}x.$
(6)

In the study regarding an isotropy and homogeneity of the model, we require a
dominant unit of the scalar field $\phi$. As dependency of scalar field is on
time therefore this may assume as perfect fluid with $\rho_{\phi}$ and
${\phi}$. And for FLRW cosmology energy density and pressure for scalar field
can be obtained as

$\rho_{\phi}=\frac{\dot{\phi}^{2}}{2}+V(\phi),$ (7)
$p_{\phi}=\frac{\dot{\phi}^{2}}{2}-V(\phi),$ (8)

where $V(\phi)$ and $\frac{1}{2}\dot{\phi}^{2}$ are potential energy and
kinetic energy depending on scalar field $\phi$. Now the conservation
equations for matter field and scalar field read as

$\dot{\rho_{m}}+3H\rho_{m}=0,$ (9)
$\dot{\rho_{\phi}}+3H(\rho_{\phi}+p_{\phi})=0.$ (10)

As EoS parameter $\omega_{\phi}=\frac{p_{\phi}}{\rho_{\phi}}$, thus Eq. (10)
leads

$\dot{\rho_{\phi}}+3H(1+\omega_{\phi})\rho_{\phi}=0,$ (11)

or

$\omega_{\phi}=-(1+a\frac{1}{3\rho_{\phi}}\frac{d\rho_{\phi}}{da}),$ (12)

Since we have four unknowns in three independent Eqs. (9), (10), (12), thus in
order to solve the system of independent equations we need one more constraint
equation. Therefore, we consider an appropriate parametrization of
$\rho_{\phi}$ as

$\rho_{\phi}(a)=e^{-\alpha a}tan^{-1}(a^{-\alpha}).$ (13)

This kind of approach is also been used by Singh et al. [39]. In this
expression $\alpha\in(0,0.2)$ is the model parameter and will be constrained
from observational datasets. The value of $\rho_{\phi}$ can be written in
terms of redshift by using the relation $a=\frac{a_{0}}{(1+z)}$, where
$a_{0}=1$ is the present value of scale factor, as

$\rho_{\phi}(z)=e^{\frac{-\alpha}{1+z}}tan^{-1}(1+z)^{\alpha},$ (14)

and the present value of energy density i.e. at z=0,

$\rho_{\phi_{0}}=\frac{\pi}{4}e^{-\alpha}.$ (15)

Eq. (14) and (15) yield

$\rho_{\phi}(z)=\frac{4}{\pi}\rho_{\phi_{0}}e^{\frac{\alpha
z}{1+z}}tan^{-1}(1+z)^{\alpha},$ (16)

Also using Eq. (9), energy density of matter field $\rho_{m}$ can be
calculated in terms of redshift $z$ as

$\rho_{m}(z)=\rho_{m_{0}}(1+z)^{3}$ (17)

Together Eqs. (4), (16) and (17), we have

$3H^{2}=\frac{3\beta^{2}}{4}+\rho_{m_{0}}(1+z)^{3}+\frac{4}{\pi}\rho_{\phi_{0}}e^{\frac{\alpha
z}{1+z}}tan^{-1}(1+z)^{\alpha},$ (18)

Now, by taking the gauge function $\beta^{2}=3\beta_{0}H_{0}^{2}a^{-2\alpha}$
and considering the density parameter $\Omega=\frac{\rho}{\rho_{c}}$ which
plays the key role to explain the whole content of the Universe. Here
$\rho_{c}=\frac{3H^{2}}{8\pi G}$ is the critical density of the Universe.

Eq. (18) can be written in terms of density parameter of scalar field and
matter,

$H=\sqrt{\frac{1}{4}\beta_{0}(1+z)^{2\alpha}+\frac{4}{3\pi}\Omega_{\phi_{0}}e^{\frac{\alpha
z}{z+1}}\tan^{-1}(z+1)^{\alpha}+\frac{1}{3}{\Omega_{m_{0}}}(z+1)^{3}},$ (19)

where $H_{0}$, $\Omega_{m_{0}}$ and $\Omega_{\phi_{0}}$ are the current values
of Hubble parameter, density parameters of matter and scalar field
respectively. We consider the Hubble parameter $H$ as dimensionless quantity
by dividing Eq. (18) by $\sqrt{3}H_{0}$, which yields Eq. (19).

The expression for deceleration parameter $q$ given by

$q=-\frac{a\ddot{a}}{\dot{a}^{2}}=-1+\frac{1+z}{H}\frac{dH}{dz}.$ (20)

Thus, the deceleration parameter $q$ and EoS parameter for scalar field
$\omega_{\phi}$ can be evaluated as

$q(z)=-\frac{-\frac{1}{2}\alpha\beta_{0}(\frac{1}{z+1})^{-2\alpha-1}-\frac{4\alpha\Omega_{\phi_{0}}e^{\alpha(1-\frac{1}{z+1})}(\frac{1}{z+1})^{-\alpha-1}}{3\pi((\frac{1}{z+1})^{-2\alpha}+1)}-\frac{4\alpha\Omega_{\phi_{0}}e^{\alpha(1-\frac{1}{z+1})}\tan^{-1}((\frac{1}{z+1})^{-\alpha})}{3\pi}-\Omega_{m_{0}}(z+1)^{4}}{2(z+1)(\frac{1}{4}\beta_{0}(\frac{1}{z+1})^{-2\alpha}+\frac{4\Omega_{\phi_{0}}e^{\alpha(1-\frac{1}{z+1})}\tan^{-1}((\frac{1}{z+1})^{-\alpha})}{3\pi}+\frac{1}{3}\Omega_{m_{0}}(z+1)^{3})}-1,$
(21)

$\displaystyle\omega_{\phi}(z)=-\frac{\pi
e^{\alpha(\frac{1}{z+1}-1)}(-\frac{4\alpha\Omega\phi
e^{\alpha(1-\frac{1}{z+1})}(\frac{1}{z+1})^{-\alpha-1}}{\pi((\frac{1}{z+1})^{-2\alpha}+1)}-\frac{4\alpha\Omega_{\phi_{0}}e^{\alpha(1-\frac{1}{z+1})}\tan^{-1}((\frac{1}{z+1})^{-\alpha})}{\pi})}{12\Omega_{\phi_{0}}(z+1)^{2}\tan^{-1}((\frac{1}{z+1})^{-\alpha})}$
$\displaystyle\times\frac{\sqrt{\frac{1}{4}\beta_{0}(\frac{1}{z+1}t)^{-2\alpha}+\frac{4\Omega\phi
e^{\alpha(1-\frac{1}{z+1})}\tan^{-1}((\frac{1}{z+1})^{-\alpha})}{3\pi}+\frac{1}{3}\Omega_{m_{0}}(z+1)^{3}}}{12\Omega_{\phi_{0}}(z+1)^{2}\tan^{-1}((\frac{1}{z+1})^{-\alpha})}-1.$
(22)

In the next section, we constrain the model parameter $\alpha$ and $\beta_{0}$
using different observational datasets and obtain their best fit value.
Thereafter use these values to examine the behaviours of various physical
parameters.

##  III Statistical observation for model parameters

Observational cosmology plays a major role to explore the evolution of the
Universe. Therefore, in the current section, we are using the expression for
the Hubble parameter from Eq. (19) to perform statistical analysis with three
observational datasets namely the Hubble dataset (77 points) [40], Type Ia
Supernovae dataset [41] and their joint datasets and to acquire the present
best fit value of the model parameters $\alpha$ and $\beta_{0}$. Here, for the
likelihoods minimization of $\chi^{2}$ for the goodness of fit of the model,
we use the Markov chain Monte Carlo (MCMC) method to investigate the parameter
space from the python package emcee [42], which is extensively used in
observational cosmology.

In this process, the present value of Hubble parameter, density parameters of
matter, and scalar field are taken as $67.4KmSec^{-1}Mpc^{-1}$, 0.315, 0.685
respectively, which are given in recent Planck data [43].

(a)

(b)

(c)

Figure 1: Figs. (a), (b) and (c) illustrate the constraints on the model for
the data $OHD$, $SNeIa$, and joint dataset $OHD+SNeIa$ respectively
corresponding to $1\sigma$ and $2\sigma$ confidence regions.

###  III.1 Supernovae type Ia Data

We use $SNeIa$ Union $2.1$ compilation dataset and constrained both the model
parameters. For flat universe Chi-square function
$\chi_{SN}^{2}(\alpha,\beta_{0},H_{0})$ can be stated as

$\chi_{SN}^{2}(\alpha,\beta_{0},H_{0})=\sum\limits_{i=1}^{580}\left[\frac{\mu_{th}(\alpha,\beta_{0},H_{0},z_{i})-\mu_{obs}(z_{i})}{\sigma_{\mu(z_{i})}}\right]^{2}.$
(23)

where $SN$ is stand for the observational $SNeIa$ dataset. $\mu_{th}$ and
$\mu_{obs}$ represents theoretical and observed distance modulus of the model.
The standard error in the observed value is denoted by $\sigma_{\mu(z_{i})}$.
Also, the theoretical distance modulus $\mu(z)$ is defined by

$\mu(z)=m-M=5LogD_{l}(z)+\mu_{0},$ (24)

where $m$ and $M$ are used for the apparent and absolute magnitudes of a
standard candle respectively. The luminosity distance $D_{l}(z)$ for flat
Universe and the nuisance parameter $\mu_{0}$ are given by

$D_{l}(z)=(1+z)c\int_{0}^{z}\frac{1}{H(z^{*})}dz^{*},$ (25)

and

$\mu_{0}=5Log\Big{(}\frac{H_{0}^{-1}}{1Mpc}\Big{)}+25,$ (26)

respectively. Here, we perform a global fitting to dictate the model
parameters using the Markov chain Monte Carlo method and Python implementation
of the ensemble sampler for the MCMC method with the EMCEE library, introduced
by Foreman-Mackey et al. to show the best-fit value of model parameters
written in Table-I.

###  III.2 Hubble Dataset

Since Hubble data is directly based on differential ages of the galaxies and
related to the expansion history of the universe thus this is very useful to
understand the dark section i.e. dark energy, dark matter, and dark ages of
the universe. The value $H$ can be written in terms of $z$ as

$H(z)=-(1+z)^{-1}\frac{dz}{dt}.$ (27)

The best-fit values of the model parameters are determined by minimizing the
Chi-square value

$\chi_{HD}^{2}(\alpha,\beta_{0})=\sum\limits_{i=1}^{77}\frac{[H(\alpha,\beta_{0},z_{i})-H_{obs}(z_{i})]^{2}}{\sigma_{z_{i}}^{2}},$
(28)

where $H(\alpha,\beta_{0},z_{i})$ and $H_{obs}$ represent the theoretical and
observed values respectively and $\sigma_{z_{i}}$ indicates the standard
deviation for every $H(z_{i})$. For $H(z)$ 77 data points, the best fit values
of $\alpha$ and $\beta_{0}$ are estimated in Table-I.

Table 1: Best fit values of model parameters. Dataset |  $\alpha$ |  $\beta_{0}$ |
---|---|---|---
$SNeIa$ |  $0.635^{+0.060}_{-0.060}$ |  $0.1935^{+0.0083}_{-0.0083}$ |
$H(z)$ |  $0.687^{+0.013}_{-0.0025}$ |  $0.157^{+0.044}_{-0.0098}$ |
$H(z)$ \+ $SNeIa$ |  $0.686^{+0.014}_{-0.0027}$ |  $0.1868^{+0.0025}_{-0.0031}$ |

The Gelman-Rubin convergence test [44] is a widely used statistical tool in
Bayesian inference that allows to access the convergence of the chains. The
test is based on the idea that multiple MCMC chains with different starting
points should converge to the same posterior distribution if they have been
run for long enough, that is, after a number of steps. Basically, this test
measures, for each parameter of the considered model, the quantity called
potential scale reduction $\hat{R}$, which is the ratio between the variance
$\rm W$ within a chain and the variance $\rm Var(\theta)$ among the chains:

$\hat{R}=\sqrt{\frac{\rm\hat{V}ar(\theta)}{W}}.$ (29)

The value of $\hat{R}$ for a perfectly converged chain should be $1$, so
$\hat{R}\approx 1.1$ corresponds to the maximum allowed value for this
parameter for which the convergence of the chain is achieved.

We then apply the Gelman-Rubin convergence test to chains of the fits of
$\alpha$ and $\beta_{0}$ from Hubble data alone. The test yielded
$\hat{R}=1.03008001$, for $\alpha$ and $\hat{R}=1.02498708$ for $\beta_{0}$,
respectively. As $1000$ steps with $100$ steps of burn-in were used to run the
MCMC chains, in Fig. 2 we verify that convergence, for both parameters, is
reached just after approximately $250$ steps when crossing the horizontal line
that represents the convergence criterion.

Figure 2: The evolution of the potential scale reduction $\hat{R}$ with the
chain length for $\alpha$ and $\beta$.

###  III.3 Joint Datasets ($H(z)+SNeIa$)

The $\chi^{2}$ function for joint analysis is given by

$\chi_{HS}^{2}=\chi_{HD}^{2}+\chi_{SN}^{2}.$ (30)

Using this joint statistical analysis, the stronger constraints of the model
parameter can be obtained. For joint datasets, the best fit values of $\alpha$
and $\beta_{0}$, can be seen in Fig. 1(c) and Table-I. Now to compare our
model with $\Lambda CDM$, we use error bar plots for Hubble datasets and
Supernovae datasets. The values of model parameters are constrained from
observational datasets. In Fig. 3, the error bar plots of the Hubble and Type
Ia Supernovae dataset show the best fit plots.

(a)  (b)

Figure 3: The error bar plots in respect of $H(z)$ and Type Ia Supernovae
datasets show the similarity between our model and $\Lambda$CDM.

##  IV Evolution of the Universe

In this section, we study the evolution of various cosmological quantities
using the best-fit values of the model parameters obtained from different
observational datasets. Energy density for the scalar field is decreasing with
respect to redshift $z$. In Fig. 4(a), It is clearly visible that
$\rho_{\phi}$ monotonically decreases and approaches zero in late times, which
indicates that amount of dark energy density of the scalar field diminished in
the future. In Fig. 4(b), the isotropic pressure of the scalar field is
negative for all values of $z$, which indicates the presence of dark energy up
to late times.

(a)  (b)

Figure 4: The graphs of $\rho_{\phi}$ and $p_{\phi}$ w.r.t. $z$.

Also the plots of EoS parameter $\omega_{\phi}$ and
$\omega_{eff}=\frac{p_{\phi}}{\rho_{\phi}+\rho_{m}}$ show the quintessence
model in late times. Finally, our model shows the quintessence model in late
times as it transits from the quintessence region to a perfect fluid state and
again converges to the quintessence region in the future (see Fig. 5).

Table 2:  The present values of the cosmological parameters. Dataset |  $H$ |  $q$ |  $\omega_{eff}$
---|---|---|---
$SNeIa$ |  $70.802$ |  $-0.1961$ |  $-0.5283$
$H(z)$ |  $70.2202$ |  $-0.1604$ |  $-0.5384$
$H(z)$ \+ $SNeIa$ |  $70.6396$ |  $-0.1640$ |  $-0.5270$

(a)  (b)

Figure 5: The graphs of $\omega_{\phi}$ and $\omega_{eff}$ vs. $z$ for
distinct observational data.

Fig. 6(a) shows the evolution of the deceleration parameter for different
observational datasets. In this model, we observe that the Universe has three
transit points during the evolution of each observational data. In early time
universe transited from deceleration to acceleration, at present, it shows the
acceleration phase. In late time again the universe is having a transition
from acceleration to deceleration and afterwords deceleration to acceleration.
These transition points are shown in Table II. In our model, we observe that
the present value of $q$ is not very near to 0.5 and it is relevant to point
out that $q\approx-0.5$ stands for $\Lambda$CDM model and constraints over
this parameter using the CMB data. Therefore, new theories of gravity may
admit different values for the deceleration parameter, as we can see in the
following works [45, 46]. Moreover, the crises over $H_{0}$ in different
regimes of redshift may also open the possibility of different values for
$q_{0}$ and $\omega_{0}$, as we can see in [47, 48].

Table 3:  Transition points of deceleration parameter for different datasets. Dataset |  $z_{tr_{1}}$ |  $z_{tr_{2}}$ |  $z_{tr_{3}}$
---|---|---|---
$SNeIa$ |  $0.7578$ |  $-0.5815$ |  $-0.8824$
$H(z)$ |  $0.7093$ |  $-0.4942$ |  $-0.8824$
$H(z)$ \+ $SNeIa$ |  $0.7189$ |  $-0.5136$ |  $-0.8824$

(a)  (b)

Figure 6: The plot of $q$ and jerk parameter w.r.t. $z$.

To understand the evolution of the expanding universe in a more significant
way, we explore some cosmographic parameters which contain the higher-order
derivatives of scale factor $a$. Here, we intend to discuss a kinematic
quantity that contains the third-order derivative and is known as the jerk
parameter. The Jerk parameter $j$ is defined as

$j=\frac{\dddot{a}}{aH^{3}}=-q+2q(1+q)+(1+z)\frac{dq}{dz}.$ (31)

The standard value of the Jerk parameter $j$ for $\Lambda$CDM model is one,
thus digression from $j=1$ explores the evolution of different kinds of dark
energy models. Fig. 6(b) highlights that at the early time our model is
similar to $\Lambda$CDM model and afterward this deviates from $\Lambda$CDM
till the late times.

To differentiate our model with $\Lambda$CDM, we also use statefinder
diagnostic technique, which is established by Sahni et al. [49, 50]. The two
geometrical diagnostic parameters $(s,r)$, known as statefinder parameters are
defined as

$r=\frac{\dddot{a}}{aH^{3}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}s=\frac{r-1}{3(q-\frac{1}{2})},~{}q\neq\frac{1}{2}$
(32)

Using statefinder diagnostic parameters, we can compare the goodness of
various dark energy models with $\Lambda$CDM. In $s-r$ plot, the points
$(0,1)$ and $(1,1)$ represent $\Lambda$CDM model and $SCDM$ respectively. From
Eqs. (21), (31) and (32), $r$ and $s$ can be calculated, and by using the
best-fit values of model parameters we draw these trajectories (see Fig.
7(a)). The arrows on the curves show the direction of evolution. In our model
$r<1$ and $s>0$ at an early time, which corresponds to the quintessence DE
model. The present values of ${(s,r)}$ are evaluated
${(0.5214,-0.0330),(0.4729,-0.0125),(0.5190,-0.0338)}$ from $OHD$, $SNeIa$ and
their joint datasets respectively, which show the deviations from $\Lambda$CDM
at present.

In Fig. 7(b), it is highlighted that the $q-r$ curves start from $SCDM$
$(1,\frac{1}{2})$ and enter into the acceleration zone for all datasets in the
future. In this plot, the dotted horizontal line $r=1$ denotes $\Lambda$CDM.
Here, we observe that the trajectories are crossing the transition line thrice
from early evolution up to late times. As $q$ is negative in late times,
therefore this model shows an acceleration phase of the universe in the
future.

(a)  (b)

Figure 7: The Statefinder plots $s-r$ and $q-r$

##  V Conclusion

In this model, we have explored the late-time cosmic expansion of the universe
in the framework of a flat FLRW space-time metric with time-dependent
displacement vectors based on Lyra’s geometry. We have derived EEF from the
principal action (1) and find solutions. Using the recent observational
datasets $OHD$, $SNeIa$ and $OHD+SNeIa$, we have obtained the best-fit values
of model parameters $\alpha$, $\beta$ by applying MCMC method with emcee
package in Python. Furthermore, we investigated the behavior of cosmological
parameters for the constrained values of model parameters, which are evaluated
in Tables I, and II.

The energy density of scalar field $\rho_{\phi}$ decreases monotonically from
high redshift to low redshift and $\rho_{\phi}\rightarrow 0$ as
$z\rightarrow-1$, $\forall$ observations. The pressure $p_{\phi}$ is having
negative finite value throughout the range, which indicates the accelerating
behavior of the Universe (see Fig. 4). The EoS parameter of the scalar field
$\omega_{\phi}$ depicts that our model is a quintessence dark energy model as
$\omega_{\phi}$ is having a value between -1 to 0 for the large range of $z$.
Its value becomes positive for a while, and at the end, it converges to -1 for
all datasets. In addition, it clearly highlighted the evolution profile of the
deceleration parameter involving three transitions for each dataset which are
shown in Table III. The jerk parameter shows that our model is similar to
$\Lambda$CDM in the early universe and it deviates in the future (see Fig.
5,6). The statefinder diagnostic graphs show that for all datasets, $s-r$
trajectories are looking in the form of loops in our model and passing through
the $\Lambda$CDM model. Also, three $q-r$ trajectories, which are also loops
for the present model, starts from SCDM and show acceleration at its final
phase (see Fig. 7). Thus, in this paper, we conclude that this model is an
accelerated expanding model having the characteristics of the quintessence
model studied in Lyra geometry, which is different from $\Lambda$CDM model. It
is relevant to point that the approach here presented could be applied in
other theories of gravity, such as $f(R)$ [51], $f(R,T)$ [52], $f(Q)$ [53,
54], and $f(Q,T)$ [55, 56]. Such applications could give us extra parameters
to test the viability of Lyra geometry with respect to experimental data. We
hope to report these extra contributions in near future.

###### Acknowledgements.

JRLS would like to thank CNPq (Grant no. 309494/2021-4), and
PRONEX/CNPq/FAPESQ-PB (Grant nos. 165/2018, and 0015/2019) for financial
support. JASF thanks FAPES for financial support.

## References

  * [1] S. Weinberg, Rev. Mod. Phys. 61, 1-23 (1989).
  * [2] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753-1936 (2006).
  * [3] H. A. Buchdahl, Mon. Not. Roy. Astron. Soc. 150, 1 (1970)
  * [4] R. Kerner, Gen. Rel. Grav. 14, 453-469 (1982).
  * [5] H. Kleinert and H. J. Schmidt, Gen. Rel. Grav. 34, 1295-1318 (2002).
  * [6] S. M. Carroll, V. Duvvuri, M. Trodden and M. S. Turner, Phys. Rev. D 70, 043528 (2004).
  * [7] S. Capozziello, S. Nojiri, S. D. Odintsov and A. Troisi, Phys. Lett. B 639, 135-143 (2006).
  * [8] L. Amendola, D. Polarski and S. Tsujikawa, Phys. Rev. Lett. 98, 131302 (2007).
  * [9] S. Nojiri and S. D. Odintsov, Phys. Rev. D 68, 123512 (2003).
  * [10] V. Faraoni, Phys. Rev. D 74, 023529 (2006).
  * [11] P. J. Zhang, Phys. Rev. D 76, 024007 (2007).
  * [12] L. Amendola and S. Tsujikawa, Phys. Lett. B 660, 125-132 (2008).
  * [13] S. Tsujikawa, Phys. Rev. D 77, 023507 (2008).
  * [14] T. Liu, X. Zhang and W. Zhao, Phys. Lett. B 777, 286-293 (2018).
  * [15] S. Capozziello and S. Tsujikawa, Phys. Rev. D 77, 107501 (2008).
  * [16] A. A. Starobinsky, JETP Lett. 86, 157-163 (2007).
  * [17] S. Nojiri and S. D. Odintsov, Phys. Lett. B 657, 238-245 (2007).
  * [18] S. Nojiri and S. D. Odintsov, Phys. Rev. D 77, 026007 (2008).
  * [19] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani and S. Zerbini, Phys. Rev. D 77, 046009 (2008).
  * [20] J. Santos, J. S. Alcaniz, M. J. Reboucas and F. C. Carvalho, Phys. Rev. D 76, 083513 (2007).
  * [21] S. Capozziello, V. F. Cardone and V. Salzano, Phys. Rev. D 78, 063504 (2008).
  * [22] O. Bertolami, C. G. Boehmer, T. Harko and F. S. N. Lobo, Phys. Rev. D 75, 104016 (2007).
  * [23] T. Harko, Phys. Lett. B 669, 376-379 (2008).
  * [24] G. Lyra, Math. Z. 54, 52 (1951).
  * [25] R. R. Cuzinatto, E. M. De Morais and B. M. Pimentel, Phys. Rev. D 103 (2021) no.12, 124002.
  * [26] A. A. Penzias and R. W. Wilson, Astrophys. J. 142, 419-421 (1965).
  * [27] K. Bhamra, Australian Journal of Physics 27, 541 (1974).
  * [28] A. Beesham, Astrophysics and Space Science 127, 355 (1986).
  * [29] W. D. Halford, Australian Journal of Physics 23, 863 (1970).
  * [30] D. R. Reddy and R. Venkateswarlu, Astrophysics and Space Science 149, 287 (1988).
  * [31] T. Singh and G. P. Singh, J. Math. Phys. 32, 2456-2458 (1991).
  * [32] T. Singh and G. P. Singh, Nuovo Cimento, B106, 617 (1991).
  * [33] T. Singh and G. P. Singh, Int. J. Theor. Phys. 31, 1433-1446 (1992).
  * [34] D. Sen, Canadian Mathematical Bulletin 3, 255 (1960).
  * [35] D. K. Sen and K. A. Dunn, J. Math. Phys. 12, no.4, 578-586 (1971).
  * [36] G. P. Singh and D. Kalyani, Pramana 49, 205-212 (1997).
  * [37] D. K. Sen, Zeitschrift f¨ur Physik 149, 311 (1957).
  * [38] J. K. Singh, Int. J. Theor. Phys. 48, 905-912 (2009).
  * [39] J. K. Singh and R. Nagpal, Eur. Phys. J. C 80, no.4, 295 (2020).
  * [40] Shaily, M. Zeyauddin and J. K. Singh, [arXiv:2207.05076 [gr-qc]].
  * [41] D. M. Scolnic et al. [Pan-STARRS1], Astrophys. J. 859 (2018) no.2, 101.
  * [42] D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, Publ. Astron. Soc. Pac. 125 (2013), 306-312.
  * [43] N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4].
  * [44] A. Gelman and D. B. Rubin, Statist. Sci. 7 (1992), 457-472
  * [45] C. P. Singh and M. Srivastava, Eur. Phys. J. C 78, no.3, 190 (2018)
  * [46] S. K. Sahu, S. K. Tripathy, P. K. Sahoo and A. Nath, Chin. J. Phys. 55, 862-869 (2017)
  * [47] P. D. Alvarez, B. Koch, C. Laporte, F. Canales and A. Rincon, JCAP 10, 071 (2022)
  * [48] D. Camarena and V. Marra, Phys. Rev. Res. 2, no.1, 013028 (2020)
  * [49] U. Alam, V. Sahni, T. D. Saini and A. A. Starobinsky, Mon. Not. Roy. Astron. Soc. 344 (2003), 1057.
  * [50] V. Sahni, T. D. Saini, A. A. Starobinsky, U. Alam, JETP Lett. 77 201 (2003).
  * [51] M. Rinaldi, G. Cognola, L. Vanzo and S. Zerbini, JCAP 08, 015 (2014).
  * [52] P. H. R. S. Moraes and J. R. L. Santos, Eur. Phys. J. C 76, 60 (2016).
  * [53] J. Beltrán Jiménez, L. Heisenberg and T. Koivisto, Phys. Rev. D 98, no.4, 044048 (2018).
  * [54] S. Mandal, P. K. Sahoo and J. R. L. Santos, Phys. Rev. D 102, no.2, 024057 (2020).
  * [55] Y. Xu, G. Li, T. Harko and S. D. Liang, Eur. Phys. J. C 79, no.8, 708 (2019).
  * [56] S. Arora, J. R. L. Santos and P. K. Sahoo, Phys. Dark Univ. 31, 100790 (2021).