# Microscopic description of $\alpha$, $2\alpha$, and cluster decays of
216-220Rn and 220-224Ra

J. Zhao  Center for Circuits and Systems, Peng Cheng Laboratory, Shenzhen
518055, China  J.-P. Ebran  Université Paris-Saclay, CEA, Laboratoire Matière
en Conditions Extrêmes, 91680, Bruyères-le-Châtel, France  L. Heitz  IJCLab,
Université Paris-Saclay, CNRS/IN2P3, 91405 Orsay Cedex, France  E. Khan
IJCLab, Université Paris-Saclay, CNRS/IN2P3, 91405 Orsay Cedex, France  F.
Mercier  IJCLab, Université Paris-Saclay, CNRS/IN2P3, 91405 Orsay Cedex,
France  T. Nikšić  Physics Department, Faculty of Science, University of
Zagreb, 10000 Zagreb, Croatia  D. Vretenar  Physics Department, Faculty of
Science, University of Zagreb, 10000 Zagreb, Croatia

(August 27, 2024)

###### Abstract

$\alpha$ and cluster decays are analyzed for heavy nuclei located above 208Pb
on the chart of nuclides: 216-220Rn and 220-224Ra, that are also candidates
for observing the $2\alpha$ decay mode. A microscopic theoretical approach
based on relativistic Energy Density Functionals (EDF), is used to compute
axially-symmetric deformation energy surfaces as functions of quadrupole,
octupole and hexadecupole collective coordinates. Dynamical least-action paths
for specific decay modes are calculated on the corresponding potential energy
surfaces. The effective collective inertia is determined using the
perturbative cranking approximation, and zero-point and rotational energy
corrections are included in the model. The predicted half-lives for
$\alpha$-decay are within one order of magnitude of the experimental values.
In the case of single $\alpha$ emission, the nuclei considered in the present
study exhibit least-action paths that differ significantly up to the scission
point. The differences in $\alpha$-decay lifetimes are not only driven by Q
values, but also by variances of the least-action paths prior to scission. In
contrast, the $2\alpha$ decay mode presents very similar paths from
equilibrium to scission, and the differences in lifetimes are mainly driven by
the corresponding Q values. The predicted 14C cluster decay half-lives are
within three orders of magnitudes of the empirical values, and point to a much
more complex pattern compared to the $\alpha$-decay mode.

##  I  Introduction

In a series of recent studies, a microscopic approach based on relativistic
Energy Density Functionals (EDF) has successfully been applied to the
phenomenon of $\alpha$ radioactivity. In this framework, self-consistent
calculations of potential energy surfaces (PES), as functions of various
deformation degrees of freedom, are used to predict the least-action dynamical
path leading to a particular decay. Using the Wentzel-Kramers-Brillouin (WKB)
approximation [1], $\alpha$-decay half-lives have been analyzed in 104Te,
108Xn, 212Po and 224Ra nuclei [2, 3]. The relativistic EDF-based microscopic
framework used to model $\alpha$ decay, has also been applied to the
description of spontaneous and induced fission [4, 5, 6, 7, 8, 9, 10, 11, 12].

This approach has recently been used to predict a new mode of decay, where two
$\alpha$-particles are simultaneously emitted back-to-back [3]. The new
$2\alpha$ decay mode should be verified experimentally, which raises the
question of the optimal candidate nucleus that must have: i) the shortest
$2\alpha$ decay lifetime, ii) the largest branching ratio, and iii) realistic
production rates at present-day exotic beam facilities. The most obvious
candidate is 216Rn, which corresponds to 2 $\alpha$ particles added to the
doubly magic 208Pb. Neighboring nuclei, such as 218Rn, and their
$\alpha$-decay ‘parent nuclei’, for instance 220,222Ra, could also be
considered.

From the theoretical point of view, it is also relevant to investigate whether
the present approach is able to describe cluster decay, which is known to
occur in this area of the nuclear chart. Since the relativistic EDF-based
method can be used to compute both fission and $\alpha$-decay observables, the
intermediate case, cluster emission, should be investigated as well. This will
also provide a link with previous calculations of cluster decay based on the
Gogny EDF, performed in a similar framework [13]. It should be noted, however,
that the relativistic EDF-based model appears to be the only fully microscopic
approach able to describe $\alpha$-decay in a quantitative agreement with
available data. This feature could be related to the ability of relativistic
EDFs to describe the formation of $\alpha$ cluster states in nuclei [14, 15,
16, 17, 18, 19, 20, 21]. Especially, the relation between the formation of
$\alpha$ clusters states and $\alpha$ decay has been explored in Refs. [22,
23].

The aim of the present work is threefold: i) to assess the quality of the
microscopic relativistic EDF-based calculations of $\alpha$ emission
lifetimes, by comparing with empirical values for the nuclei of interest,
namely 216-220Rn and 220-224Ra, as well as 212Po; ii) to explore whether the
relativistic EDF approach could also predict cluster decay, and the level of
agreement with data, as well as possible open questions raised by a universal
microscopic description of fission, cluster, $2\alpha$ and $\alpha$ decays;
and iii) to analyze possible $2\alpha$ decays of nuclei in the region
${\rm{}^{208}Pb}+2\alpha={\rm{}^{216}Rn}$, which should be favored as a 2
$\alpha$ emitter candidate. On this last point, 216Rn, but also its
$\alpha$-parent nucleus, as well as slightly more neutron-rich nuclei, will be
analyzed. The aim is to identify optimal nuclei that also have large
production rates, usually from a primary reaction on a Uranium target.

The present study is based on calculations performed in the Relativistic
Hartree-Bogoliubov (RHB) framework. For a detailed review, we refer the reader
to Refs. [24, 25, 26]. Self-consistent calculations of deformation energy
surfaces are performed using the DD-PC1 [27] relativistic functional and a
separable finite-range pairing interaction. The potential energy surfaces
(PES) are calculated using quadrupole, octupole and hexadecupole deformations
as collective degrees of freedom. Section II describes the formalism and the
details of the calculation. Predictions of $\alpha$-decay half-lives for
216-220Rn and 220-224Ra are discussed in section III. Section IV explores the
symmetric $2\alpha$ decay mode, and in section V we analyze 14C cluster
emission from 222Ra and 224Ra. A brief summary of the principal results is
included in section VI.

##  II  Relativistic Hartree-Bogoliubov calculations of decay rates

The process of decay by an $\alpha$-particle or cluster emission is modelled
along a dynamical path $L$, determined by minimizing the action integral [28,
29]:

$S(L)=\int_{s_{\rm in}}^{s_{\rm out}}{1\over\hbar}\sqrt{2\mathcal{M}_{\rm
eff}(s)\left[V_{\rm eff}(s)-E_{0}\right]}ds,$ (1)

where $\mathcal{M}_{\rm eff}(s)$ and $V_{\rm eff}(s)$ are the effective
collective inertia and potential, respectively. $E_{0}$ is the ground-state
energy of the nucleus, and the integration limits correspond to the classical
inner ($s_{\rm in}$) and outer turning points ($s_{\rm out})$, defined by
$V_{\rm eff}(s)=E_{0}$.

The effective inertia is computed from the multidimensional inertia tensor
$\mathcal{M}$ [28, 30, 31, 32, 33]

$\mathcal{M}_{\rm eff}(s)=\sum_{ij}\mathcal{M}_{ij}{dq_{i}\over
ds}{dq_{j}\over ds}\;,$ (2)

where $q_{i}(s)$ denotes the collective coordinate as a function of the path’s
length. The collective inertia tensor is calculated from the self-consistent
and deformation-constrained RHB solutions for the quasiparticle wave functions
and energies, using the adiabatic time-dependent Hartree-Fock-Bogoliubov
(ATDHFB) method [34]. In the perturbative cranking approximation, the
collective inertia reads [35]

$\mathcal{M}=\hbar^{2}{\it M}_{(1)}^{-1}{\it M}_{(3)}{\it M}_{(1)}^{-1},$ (3)

where

$\left[{\it M}_{(k)}\right]_{ij}=\sum_{\mu\nu}{\left\langle
0\left|\hat{Q}_{i}\right|\mu\nu\right\rangle\left\langle\mu\nu\left|\hat{Q}_{j}\right|0\right\rangle\over(E_{\mu}+E_{\nu})^{k}}\;.$
(4)

$|\mu\nu\rangle$ are two-quasiparticle states, and $E_{\mu}$ and $E_{\nu}$ the
corresponding quasiparticle energies. $\hat{Q}_{i}$ denotes the multipole
operators that represent collective degrees of freedom. The effective
collective potential $V_{\rm eff}$ is obtained by subtracting the vibrational
zero-point energy (ZPE) from the total RHB deformation energy. Following the
prescription of Refs. [36, 37, 32, 33], the ZPE is computed using the Gaussian
overlap approximation,

$E_{\rm ZPE}={1\over 4}{\rm Tr}\left[{\it M}_{(2)}^{-1}{\it M}_{(1)}\right].$
(5)

In addition, the rotational energy correction (REC) is also subtracted

$E_{\rm REC}=\frac{<J^{2}>}{2\mathcal{I}_{Y}}\;,$ (6)

where the expectation values of the angular momentum operator are evaluated
with respect to the self-consistent and deformation-constrained RHB states,
and $\mathcal{I}_{Y}$ is the Yoccoz moment of inertia [38] for a given
multipole operator $\hat{Q}_{i}$:

$\mathcal{I}_{Y}(\hat{Q}_{i})=\sum_{\mu\nu}{2\left(\left|\left\langle
0\left|\hat{Q}_{i}\right|\mu\nu\right\rangle\right|^{2}\right)^{2}\over\left|\left\langle
0\left|\hat{Q}_{i}\right|\mu\nu\right\rangle\right|^{2}(E_{\mu}+E_{\nu})}\;.$
(7)

The REC approximation of Eq. (6) is known to be valid for large deformations
[13], and replaces full angular momentum projection which, in the case of
heavy nuclei, is generally computationally prohibitive.

It should be noted that the REC was not taken into account in our previous
calculations of $\alpha$ decay based on relativistic EDFs [2, 3]. Therefore,
the present results should be, in principle at least, more accurate. A typical
effect of this correction is to decrease the lifetimes of $\alpha$-decay by
about one order of magnitude. For instance, in the case of 224Ra, the
inclusion of the REC leads to a decrease of the $\alpha$ emission lifetime
from 9.5 d to 5.7 d. We have verified that the rotational energy correction
globally improves the results in comparison to data. Therefore, in the present
work, all calculations are performed by including the REC.

The microscopic self-consistent solutions of the deformation-constrained RHB
equations, that is, the single-quasiparticle energies and wave functions on
the entire energy surface as functions of the quadrupole, octupole, and/or
hexadecapole deformations, provide the microscopic input for computing the
collective inertia, zero-point energy and rotational energy corrections.

In practice, the least-action path is built from the inner turning point to
the scission point, whose position is determined by monitoring the integrated
density distribution, that is, the density distribution of the emerging
fragment. As scission point, we select the point with emerging fragment mass
equal to the mass of the $\alpha$ particle, $2\alpha$ particles, or 14C
cluster. For the single $\alpha$ or cluster emission, beyond scission the
configuration with two well separated fragments becomes the lowest energy
solution, and the energy curve up to the outer turning point $s_{\rm out}$ can
be approximated by the classical expression for two uniformly charged spheres:

$V_{\rm eff}(\beta_{3})=e^{2}{Z_{1}Z_{2}\over R}-Q,$ (8)

where $R$ represents the distance between the centers of mass of the
fragments, and the second term is the experimental $Q$ value. We use Eqs. (9)
and (10) of Ref. [13] to approximate the relation between $R$ and the octuple
moment $Q_{30}$,

$Q_{30}=f_{3}R^{3},$ (9)

with

$f_{3}={A_{1}A_{2}\over A}{(A_{1}-A_{2})\over A},$ (10)

and $\beta_{30}=4\pi Q_{30}/3AR^{3}$. The corresponding effective collective
mass reads

$\mathcal{M}_{\rm eff}={\mu\over 9Q^{4/3}_{30}f^{2/3}_{3}},$ (11)

where $\mu=m_{n}A_{1}A_{2}/(A_{1}+A_{2})$ is the reduced mass of the two
fragments, and $m_{n}$ denotes the nucleon mass [13].

In the case of emission of two $\alpha$ particles, to calculate the
contribution of the action from the scission point to the outer turning point
$s_{\rm out}$, we consider the superposition of each $\alpha$-plus-nucleus
Coulomb interaction, namely:

$V_{\rm eff}(\beta_{2})=2e^{2}{Z_{1}Z_{2}\over R}-Q_{2\alpha},$ (12)

where $R$ represents the distance between the centers of mass of the fragments
(the index 1 refers to the heavy fragment, and 2 to the $\alpha$-particle).
The approximate relation between R and the quadrupole moment is

$Q_{20}=2A_{2}R^{2}\;.$ (13)

The corresponding effective collective mass reads

$\mathcal{M}_{\rm eff}={\mu\over 8A_{2}Q_{20}},$ (14)

where $\mu=m_{n}A_{2}/2$ is the reduced mass of the (2$\alpha$\+ heavy)
fragments.

Therefore, the path involved in the action integral of Eq. (1) consists of the
least-action path from $s_{\rm in}$ to scission, and the energy is
approximated by the Coulomb potential from scission to $s_{\rm out}$ [13]. The
decay half-life is calculated as $T_{1/2}=\ln 2/(nP)$, where $n$ is the number
of assaults on the potential barrier per unit time [30, 31, 32, 33], and $P$
is the barrier penetration probability in the WKB approximation

$P={1\over 1+\exp[2S(L)]}.$ (15)

We choose $E_{0}=1$ MeV in Eq. (1) for the value of the collective ground
state energy. For the vibrational frequency $\hbar\omega=1$ MeV, the
corresponding value of $n$ is $10^{20.38}$ s-1.

In the present study, for the single $\alpha$ emission the collective space is
3-dimensional $(\beta_{20},\beta_{30},\beta_{40})$, while for the symmetric
2$\alpha$ emission, because of reflection symmetry, the collective space can
be built from the coordinates $\beta_{20}$ and $\beta_{40}$. In the case of
14C emission, two models will be analyzed: i) a 2D $(\beta_{20},\beta_{30})$
calculation, where the scission point is determined by a discontinuity in
$\beta_{40}$, to be compared with previous Gogny EDF-based results [13] and,
in the case of 222Ra ii) a full 3D $(\beta_{20},\beta_{30},\beta_{40})$
calculation, in which scission is determined as a point where the integrated
density of the emitted fragment corresponds to 14 nucleons.

We note that the present deformation-constrained calculations are performed
using a method with linear constraints that has successfully been applied to
fission (see Ref. [35] for details). The Dirac-Hartree-Bogoliubov equations
are solved by expanding the nucleon spinors in the basis of a 3D harmonic
oscillator, with 20 oscillators shells. In addition to the particle-hole
channel determined by the choice of the EDF, a separable pairing interaction
of finite-range [40, 41] is used that reproduces the pairing gap in nuclear
matter as calculated with the D1S parametrization of the Gogny force [42, 41].
As in Ref. [3], we have fine-tuned the neutron and proton pairing strengths to
reproduce the empirical pairing gaps of the isotope 224Ra. Compared to the
original values, this corresponds to an increase of the neutron and proton
pairing strengths by 9% and 12%, respectively. This modification is consistent
with the conclusions of the recent global study of the separable pairing
interaction, when used with relativistic energy density functionals [43, 44].
We have used the modified pairing strength parameters for all the Ra and Rn
nuclei considered in the present work.

##  III  $\alpha$-decay

In the case of $\alpha$ decay, the importance of including at least three
deformation degrees of freedom has been shown in our previous calculations [2,
3]. The quadrupole degree of freedom corresponds to the elongation of the
nucleus, while octupole deformations are essential to model the formation of
an $\alpha$ particle on the surface of a nucleus. The formation of a neck
between the $\alpha$ particle and the daughter nucleus before emission, gives
rise to a more pronounced hexadecapole deformation. For the $\alpha$ decay of
220Ra, figure 1 displays several projections of the PES in the axially
symmetric quadrupole-octupole plane, for different values of the hexadecapole
deformation parameter $\beta_{40}$. The least-action path remains close to the
energy minimum of the deformation surfaces, with some deviations due to the
effect of the collective inertia. Along the path from the inner ($s_{\rm in}$)
turning point to the scission point, the value of $\beta_{20}$ increases
slightly from 0.12 to 0.16. The increase of $\beta_{30}$ is more pronounced,
from 0.12 to 0.32, and even more so for $\beta_{40}$ from 0.2 to 0.7. Hence,
in the case of $\alpha$ decay from 220Ra, the major deformation degrees of
freedom at work are the octupole and hexadecapole ones. The scission point on
the ($\beta_{20}$, $\beta_{30}$, $\beta_{40}$) surface is found at
($0.16,0.31,0.68$), and the least-action integral up to this point has a value
of 8.10, whereas the one from scission to the outer turning point ($s_{\rm
out})$ is 14.18. The predicted half-life is 60 ms, to be compared with the
experimental value of 18 ms.

Figure 1: Deformation energy surface of 220Ra in the axially-symmetric
quadrupole-octupole plane, for selected values of the hexadecapole deformation
$\beta_{40}$. Calculations are performed using the RHB model based on the DD-
PC1 functional and a separable pairing interaction. Contours join points on
the surface with the same energy, and red circles indicate the points on the
dynamical (least-action) path for $\alpha$ emission. The insets display the
intrinsic nucleon densities at selected points on the dynamical path.

It is interesting to compare the predictions for the half-lives of
$\alpha$-decay of 220Ra and 212Po. In the latter case, the scission point is
very close to the one of 220Ra: (0.16,0.33,0.70). The least-action integral
calculated form the inner turning point to scission is 6.58, and 9.46 from
scission to the outer turning point. We note that the action up to the
scission point is smaller than in the case of 220Ra. This is due to the almost
spherical shape of the equilibrium minimum of 212Po: the least-action path
starts from a small $\beta_{20}$ value, close to zero, because of the vicinity
of the doubly magic nucleus 208Pb. The difference between the calculated
least-action integrals is then magnified by the well-known Q-value effect
after the scission point, making the total action much smaller in the case of
212Po. About 25% of the difference in the total action between 220Ra and 212Po
is due to nuclear dynamics before the scission point, rather than to the
Q-value effect. The correct description of the dynamical path before the
scission point is crucial, as demonstrated by the predicted half-life: 0.2
$\mu$s for 212Po, to be compared to the experimental value of 0.3 $\mu$s.

Figure 2 illustrates the calculation of the $\alpha$-decay of 222Ra. The
evolution of the path is similar to the one of 220Ra. The scission point is
obtained at ($0.15,0.31,0.68$) on the ($\beta_{20}$, $\beta_{30}$,
$\beta_{40}$) surface, which is almost identical to the case of 220Ra. The
least-action integral up to the scission point is 9.03, and the one in the
outer region starting from this point is 17.02. The corresponding half-life is
122 s, in comparison to the experimental value 38 s. Again, a correct
description of the least-action path prior to the scission point is crucial to
predict a realistic half-life within one order of magnitude of the empirical
value.

Figure 2: Same as in the caption to Fig. 1 but for 222Ra

For the $\alpha$-decay of 216Rn and 218Rn, in figures 3 and 4, respectively,
we display the corresponding projections of the 3D deformation energy
surfaces. The case of 216Rn is similar to the one of 212Po, with a starting
point of the least-action path close to sphericity. The scission point is
found at ($0.16,0.32,0.70$) on the ($\beta_{20}$, $\beta_{30}$, $\beta_{40}$)
surface, the action integral up to the scission point is 8.12, and the one in
the outer region 11.83. The predicted half-life is 185 $\mu$s, and the
experimental one is 45 $\mu$s. For 218Rn, the PES is comparable to those of
the Ra nuclei discussed above. The scission point is at ($0.16,0.32,0.70$),
the action up to the scission point is 8.58, and the one from this point is
14.32. Here, despite the difference in the dynamical paths before the scission
point (cf. Figs. 3 and 4), the calculated actions before the scission point
are close, though a bit larger in the case of 218Rn. After this point, because
of the smaller Q-value, the action is also larger in the case of 218Rn. The
calculated half-life 35ms is in agreement with the experimental value.

Figure 3: Same as in the caption to Fig. 1 but for 216Rn

Figure 4: Same as in the caption to Fig. 1 but for 218Rn

The case of 220Rn is close to the 218Rn one, with similar PES and small
variations on the dynamical path. The scission point is at ($0.12,0.32,0.70$),
the action up to the scission point is 8.03, and the one from this point is
17.48. This larger value, compared to the previous 216,218Rn cases, is due to
a smaller Q-value. The predicted half-life is 40.8 s, and the experimental one
is 55.6 s.

Table 1 summarizes the $\alpha$ decay half-lives calculated using the present
method. All the theoretical values are of the same order of magnitude as the
experimental ones. This result demonstrates the validity of the microscopic
approach, which explicitly takes into account various important effects in a
consistent framework: pairing, collective inertia, energy corrections, etc. It
should also be noted that there are no free parameters specifically adjusted
to reproduce the empirical half-lives. The only input is a universal energy
density functional, and the pairing interaction.

Table 1:

Experimental and theoretical values for the $\alpha$-emission half-lives,
predicted 2$\alpha$ symmetric emission half-lives, and the corresponding
branching ratios. BR denotes the log10 of the branching ratio.  Texp Tα log10
(T2α[s]) BR 212Po 0.3 $\mu$s 0.2 $\mu$s 17.23 -23.8 216Rn 45 $\mu$s 185 $\mu$s
2.67 -6.4 218Rn 35 ms 35 ms 5.06 -6.5 220Rn 55.6 s 40.8 s 8.04 -6.4 220Ra 18
ms 60 ms 6.1 -7.3 222Ra 38 s 122 s 9.23 -7.1 224Ra 3.6 d 5.7 d 13.03 -7.3

##  IV  Symmetric two-$\alpha$ decay

The simultaneous symmetric emission of two $\alpha$ particles is a new mode of
decay which has been predicted recently [3]. In this process, the two $\alpha$
particles are emitted back to back, and the least-action path is calculated on
the surface of quadrupole and hexadecapole deformations. The octupole degree
of freedom plays no role here. Considering the condition of positive Q2α
value, almost all of the one-$\alpha$ emitter nuclei are predicted to be also
2-$\alpha$ emitters, with a few exceptions [3]. However, the branching ratio
for this mode is generally very small, about 10-8 in 224Ra. This is
nevertheless of the same order of magnitude as the branching ratios for
cluster radioactivity, which has already been detected. Therefore, it is
necessary to focus on the best 2 $\alpha$-decay candidates for a possible
measurement. Figure 5 displays the Q2α values on the nuclear chart, for nuclei
exhibiting positive Q2α values, where Q2α is defined by:

$Q_{2\alpha}=Q_{\alpha 1}+Q_{\alpha 2}+\Delta E\;.$ (16)

Qα1(2) are the Q-values of the daughter and grand-daughter nuclei for single
$\alpha$ decay, determined from experimental masses [47]. $\Delta$E is the
difference between the sum of the excitation energies of the daughter and
grand-daughter nuclei in the sequence of single $\alpha$ decays, and the
excitation energy of the grand-daughter nucleus in the $2\alpha$ decay. We
consider, for simplicity, $\alpha$ or $2\alpha$ transitions involving only the
ground states of the daughter or grand-daughter nuclei (that is, $\Delta
E$=0). Considering the empirical Q2α values in Fig. 5, it appears that the
most favored region for $2\alpha$ decay is located in a narrow strip north-
west of 208Pb. The optimal experimental candidate should be as close as
possible to this region, with an intensity of production large enough to
detect a few possible $2\alpha$ decays. Unstable nuclei around this region are
best produced from a Uranium primary beam. Taking into account these
constraints, some of the best candidates could be 216Rn, which corresponds to
208Pb+2$\alpha$, followed by 218Rn, 220Ra, 222Ra, and 224Ra. We note that
nuclei around 218Th or 222U are even more favorable when considering Q2α
values. However, their current production rate is too low to perform dedicated
experiments on the $2\alpha$ decay mode.

Figure 5: $N>50$ and $Z>50$ nuclides with positive Q2α values.

In Figs. 6-9 we show the RHB deformation energy surfaces in the quadrupole-
hexadecapole plane for the nuclei 220,222Ra and 216,218Rn. In these four
cases, the least-action paths for the symmetric $2\alpha$ emission are very
similar: a slight initial decrease of $\beta_{20}$, followed by a linear
increase of both $\beta_{20}$ and $\beta_{40}$, leading to very large values
of the hexadecapole deformation parameter at scission. The scission point is
located at $\beta_{20}=0.30$ and $\beta_{40}=1.34$ in the case of 220Ra, and
almost the same for the other three nuclei (less than 2% change in the
deformation parameters). Correspondingly, the least-action integrals up to the
scission point are also rather similar: 13.93, 14.36, 12.70, and 12.61 for
220,222Ra and 216,218Rn, respectively. However, it is the Q2α values that play
the main role in differentiating the $2\alpha$ decay rates in these nuclei.
The least-action integrals from the scission point to the outer turning points
are: 16.76, 19.92, 14.02, and 16.87 for 220,222Ra and 216,218Rn, respectively.
The corresponding half-lives, listed in Table 1, can directly be correlated to
these values of the action integral: the shortest half-life is calculated for
216Rn, whereas the longest one is for 222Ra. The above conclusions also apply
to the case of 220Rn with an action before the scission point of 13.01, the
one after the scission point being 19.90.

Figure 6: Reflection symmetric deformation energy surface of 220Ra in the
quadrupole-hexadecapole axially-symmetric plane. Calculations have been
performed using the RHB model based on the DD-PC1 functional, and a separable
pairing interaction. The grey curve denotes the dynamical (least-action) path
for 2$\alpha$ emission from the equilibrium deformation to scission. The
insets display the intrinsic nucleon densities at selected points on the
dynamical path.

Figure 7: Same as in the caption to Fig. 6 but for 222Ra.

Figure 8: Same as in the caption to Fig. 6 but for 216Rn.

Figure 9: Same as in the caption to Fig. 6 but for 218Rn.

The predicted $2\alpha$ decay branching ratios (BR) listed in Table 1, show
that the most unfavorable case for this mode appears to be 212Po. This could
be understood because of the 208Pb+1$\alpha$ configuration of this nucleus. In
contrast, 216Rn, which corresponds to 208Pb+2$\alpha$, is predicted to be the
optimal case for simultaneous symmetric $2\alpha$ decay. However, neighboring
nuclei such as 218,220Rn or 220-224Ra also exhibit BRs that do not differ by
more than an order of magnitude. This is due to the simultaneous increase of
the action integrals in both the $\alpha$ and $2\alpha$ channels, when
compared to 216Rn. In summary, 216-220Rn and 220-224Ra should be considered as
candidates for an experimental search of the $2\alpha$ decay mode.

It would also be interesting to consider additional geometric degrees of
freedom, such as triaxiality, but this is beyond the scope of this work, and
hardly feasible from a computational point of view. However, the inclusion of
additional degrees of freedom would generally decrease the action calculated
along the decay path and, hence, the predicted half-lives for the two-$\alpha$
decay. The half-lives calculated in the present study should be considered as
upper limits for the two-$\alpha$ decay mode.

Three-body effects are also relevant for two-particle, e.g., two-proton decay.
However, it is difficult to draw a direct analogy between two-proton and
two-$\alpha$ decays, because preformation of $\alpha$ clusters has a large
impact on the decay probability which, of course, is not the case for two-
proton decay. In the present calculation, all possible interactions are
considered before the scission point, because the least action path is
calculated starting from nucleonic degree of freedom, displaying the emergence
of the two $\alpha$ along this path. After the scission point, each
$\alpha$-nucleus interaction is taken into account, while the interaction
between the two $\alpha$ is neglected. This corresponds to a simplified case
of a full three-body treatment. Translated in the description of two-proton
decay, this would rather correspond to the so-called direct emission model
[45], which yields a correct order of magnitude for the decay half-lives [46],
although only the full three-body treatment provides the correct momentum
density distribution of the emitted protons.

##  V  Cluster decay

As an example of cluster decay, we consider the emission of 14C from 222Ra and
224Ra. Cluster radioactivity was previously explored in a study based on
least-action calculations with the Gogny energy density functional [13]. Fig.
10 displays the reflection asymmetric deformation energy surface of 224Ra in
the quadrupole-octupole axially-symmetric plane. The grey curve denotes the
dynamical least-action path for 14C emission from equilibrium deformation to
scission. The scission point is located at
($\beta_{20}=0.45,\beta_{30}=0.82$). The corresponding contribution to the
dimensionless action integral (1) is 13.97, and from the scission to the outer
turning point the action is 27.95. The corresponding cluster half-live log10
$\left(T_{{}^{14}C}\left[s\right]\right)=15.87$, is in excellent agreement
with the experimental value log10 $\left(T_{exp}\left[s\right]\right)=15.86$.
In a similar calculation with the Gogny functional, the value log10
$\left(T_{{}^{14}C}\left[s\right]\right)=15.06$ (cf. Fig. 13 of [13]) was
obtained for the 14C decay half-life of 224Ra.

Figure 10: Reflection symmetric deformation energy surface of 224Ra in the
quadrupole-octupole axially-symmetric plane. The grey curve denotes the
dynamical (least-action) path for 14C emission from equilibrium deformation to
scission.

In Fig. 11 we show the corresponding deformation energy surface of 222Ra, and
the least action path followed till the scission point, located at
($\beta_{20}=0.48,\beta_{30}=0.79$), close to the value obtained for 224Ra.
The corresponding contribution to the dimensionless action (1) is 15.28 up to
the scission point, and 24.04 from the scission point to the outer turning
point. The predicted cluster half-live is log10
$\left(T_{{}^{14}C}\left[s\right]\right)=13.61$, to be compared with the
experimental value log10 $\left(T_{exp}\left[s\right]\right)=11.01$. This
appears to be on the same level of agreement with experiment as the result
obtained using the Gogny functional: log10
$\left(T_{{}^{14}C}\left[s\right]\right)=8.9$ [13]. The energy surfaces and
the least-action paths for 14C emission are very similar for 222Ra and 224Ra.
A Q-value that is 3 MeV larger in the case of 222Ra, explains its smaller
action from the scission point to emission, and correspondingly its shorter
cluster emission half-life.

Figure 11:  Same as in the caption to Fig. 10 but for 222Ra.

Just as in the case of $\alpha$-decay, the hexadecupole degree of freedom
could also play an important role for heavier cluster emission, such as 14C.
We have thus carried out 3D hexadecapole-octupole-quadrupole calculations of
14C emission from 222Ra. Fig. 12 displays the projections of the deformation
energy surface onto the quadrupole-octupole plane, for several values of the
hexadecapole parameter $\beta_{40}$. Compared to $\alpha$ emission, the
deformation parameters reach much larger values: the scission point on the
surface ($\beta_{20}$, $\beta_{30}$, $\beta_{40}$) is at ($0.49,0.87,1.66$).
The least-action integral up to the scission point is 18.24, and the one form
the scission point to the outer turning point is 22.48. The predicted half-
life is log10 $\left(T_{{}^{14}C}\left[s\right]\right)=14.82$. Hence, the
inclusion of the hexadecapole degree of freedom produces a limited effect,
with an increase of about one order of magnitude for the half-life. This
increase is mainly due to the different locations of the scission point in 2D
and 3D calculations.

Figure 12: Deformation energy surface of 222Ra in the quadrupole-octupole
axially-symmetric plane, for selected values of the hexadecapole deformation
$\beta_{40}$. Contours join points on the surface with the same energy and red
circles indicate the points on the dynamical (least-action) path for 14C
emission. The insets display the intrinsic nucleon densities at selected
points on the dynamical path.

In the case of heavier cluster emission, a better agreement with data could
possibly be obtained by improving the treatment of collective inertia. It
would be interesting to explore the effect of the perturbative and
nonperturbative cranking approximations for the multidimensional inertia
tensor, on the calculated half-lives for cluster decay. However, such an
analysis is numerically rather involved and beyond the scope of the present
study. Table 2 summarizes the predicted 14C half-lives, and Fig 13 presents an
overview of the theoretical results for the half-lives of 216-220Rn, 220-224Ra
and 212Po, in comparison with available data. The overall good agreement with
the data, especially in the $\alpha$-decay case, shows that the present
approach provides a consistent microscopic framework for calculating $\alpha$
and cluster emission lifetimes, taking explicitly into account the dynamical
path on a multidimensional deformation energy surface, the effective
collective mass, and the pairing interaction.

Table 2:

Experimental, 2D and 3D calculated log10 T[s] of the half-lives for 14C
cluster emission in the case 222Ra and 224Ra.
Texp T${}^{2D}_{C}$ T${}^{3D}_{C}$ 222Ra 11.01 13.61 14.82 224Ra 15.86 15.87

Figure 13: Theoretical results for the $\alpha$, $2\alpha$ and cluster decay
half-lives of 216-220Rn, 220-224Ra and 212Po, in comparison with available
experimental values.

##  VI  Summary

$\alpha$, 2-$\alpha$ and 14C spontaneous emissions from heavy nuclei have been
explored using a microscopic framework, based on the RHB model with the DD-PC1
functional and a separable pairing force. The corresponding 3D deformation
energy surfaces (quadrupole, octupole, and hexadecapole collective degrees of
freedom) enable the calculation of least-action integrals along dynamical
paths. The predicted half-lives for $\alpha$-decay are within one order of
magnitude of the experimental values. The predictions for heavier cluster
emission deviate from the experimental values by 1 to 4 orders of magnitude.
This may point to a possible improvement in the calculation of the effective
collective inertia in the latter case.

Two main effects impact the decay half-life: the dynamical path from the inner
turning point up to the scission point and, for the part of the path beyond
the scission point, the Q-value of the transition. The relative importance of
these two effects has been analyzed. $2\alpha$ symmetric decay has been
predicted within the same framework. All nuclei considered in the present
study are also predicted to be $2\alpha$ emitters. The corresponding branching
ratio is more favorable for the Rn isotopes, but the Ra nuclei appear also as
good candidates for the experimental observation of this rare decay mode. This
has been confirmed in a recent phenomenological study of the double
$\alpha$-decay half-life [48].

###### Acknowledgements.

This work has been supported in part by the QuantiXLie Centre of Excellence, a
project co-financed by the Croatian Government and European Union through the
European Regional Development Fund - the Competitiveness and Cohesion
Operational Programme (KK.01.1.1.01), and by the Institut Universitaire de
France. J. Z. acknowledges support by the National Natural Science Foundation
of China under Grant No. 12005107 and No. 11790325.

## References

  * [1] D. S. Delion, Lect. Notes in Physics 819, 1 (2010).
  * [2] F. Mercier, J. Zhao, R.-D. Lasseri, J.P. Ebran, E. Khan, T. Nikšić and D. Vretenar, Phys. Rev. C 102, 011301(R) (2020).
  * [3] F. Mercier, J. Zhao, J.P. Ebran, E. Khan, T. Nikšić and D. Vretenar, Phys. Rev. Lett. 127, 012501 (2021).
  * [4] J. Zhao, B.-N. Lu, T. Nikšić, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 93, 044315 (2016).
  * [5] H. Tao, J. Zhao, Z. P. Li, T. Nikšić, and D. Vretenar, Phys. Rev. C 96, 024319 (2017) .
  * [6] J. Zhao, T. Nikšić, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 99, 014618 (2019).
  * [7] J. Zhao, J. Xiang, Z.-P. Li, T. Nikšić, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 99, 054613 (2019).
  * [8] J. Zhao, T. Nikšić, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 101, 064605 (2020).
  * [9] J. Zhao, T. Nikšić, and D. Vretenar, Phys. Rev. C 104, 044612 (2021).
  * [10] J. Zhao, T. Nikšić, and D. Vretenar, Phys. Rev. C 105, 054604 (2022).
  * [11] Z. X. Ren, J. Zhao, D. Vretenar, T. Nikšić, P. W. Zhao, and J. Meng, Phys. Rev. C 105, 044313 (2022).
  * [12] Z. X. Ren, D. Vretenar, T. Nikšić, P. W. Zhao, J. Zhao, and J. Meng, Phys. Rev. Lett. 128, 172501 (2022).
  * [13] M. Warda, and L. M. Robledo, Phys. Rev. C 84, 044608 (2011).
  * [14] J.P. Ebran, E. Khan, T. Nikšić and D.Vretenar, Nature 487, 341 (2012).
  * [15] J.P. Ebran, E. Khan, T. Nikšić and D. Vretenar, Phys. Rev. C 87, 044307 (2013).
  * [16] J.P. Ebran, E. Khan, T. Nikšić and D. Vretenar, Phys. Rev. C 90, 054329 (2014).
  * [17] J.P. Ebran, E. Khan, T. Nikšić and D. Vretenar, Phys. Rev. C 89, 031303(R) (2014).
  * [18] P.W. Zhao, N. Itagaki, J. Meng, Phys. Rev. Lett 115, 022501 (2015).
  * [19] E.F. Zhou, J.M. Yao, Z.P. Li et al., Phys. Lett. B 753, 227 (2016).
  * [20] P.Marević, J.P. Ebran, E. Khan, T. Nikšić and D. Vretenar, Phys. Rev. C 97, 024334 (2018).
  * [21] P.Marević, J.P. Ebran, E. Khan, T. Nikšić and D. Vretenar, Phys. Rev. C 99, 034317 (2019).
  * [22] J.P. Ebran, E. Khan, R.-D. Lasseri, and D. Vretenar, Phys. Rev. C 97, 061301(R) (2018).
  * [23] J.P. Ebran, E. Khan, and R.-D. Lasseri, J. Phys. G 48, 025106 (2021).
  * [24] D. Vretenar, A.V. Afanasjev, G. A. Lalazissis, P. Ring, Phys. Rep. 409, 101 (2005).
  * [25] J. Meng, Relativistic Density Functional for Nuclear Structure, World Scientific (2016).
  * [26] J.P. Ebran, in Energy Density Functional Methods for Atomic Nuclei, IOP Publishing (2019).
  * [27] T. Nikšić, D. Vretenar, and P. Ring, Phys. Rev. C 78, 034318 (2008).
  * [28] M. Brack, J. Damgaard, A. S. Jensen, H. C. Pauli, V. M. Strutinsky, and C. Y. Wong, Rev. Mod. Phys. 44, 320 (1972).
  * [29] T. Ledergerber and H.-C. Pauli, Nucl. Phys. A 207, 1 (1973).
  * [30] A. Baran, Phys. Lett. B 76, 8 (1978).
  * [31] A. Baran, K. Pomorski, A. Lukasiak, and A. Sobiczewski, Nucl. Phys. A 361, 83 (1981).
  * [32] J. Sadhukhan, K. Mazurek, A. Baran, J. Dobaczewski, W. Nazarewicz, and J. A. Sheikh, Phys. Rev. C 88, 064314 (2013).
  * [33] J. Sadhukhan, J. Dobaczewski, W. Nazarewicz, J. A. Sheikh, and A. Baran, Phys. Rev. C 90, 061304(R) (2014).
  * [34] A. Baran, J. A. Sheikh, J. Dobaczewski, W. Nazarewicz, and A. Staszczak, Phys. Rev. C 84, 054321 (2011).
  * [35] J. Zhao, B.-N. Lu, T. Nikšić, and D. Vretenar, Phys. Rev. C 92, 064315 (2015).
  * [36] A. Staszczak, A. Baran, and W. Nazarewicz, Phys. Rev. C 87, 024320 (2013).
  * [37] A. Baran, A. Staszczak, J. Dobaczewski, and W. Nazarewicz,Int. J. Mod. Phys. E 16, 443 (2007).
  * [38] J. L. Egido and L. M. Robledo, Lect. Notes Phys. 641, 269 (2004).
  * [39] P. Ring and P. Schuck, The nuclear many-body problem, Springer-Verlag (1980).
  * [40] T. Duguet, Phys. Rev. C69, 054317 (2004)
  * [41] Y. Tian, Z. Y. Ma, and P. Ring, Phys. Lett. B 676, 44 (2009).
  * [42] J.-F. Berger, M. Girod, and D. Gogny, Comput. Phys. Commun. 63, 365 (1991).
  * [43] L. J. Wang, B. Y. Sun, J. M. Dong, and W. H. Long, Phys. Rev. C 87, 054331 (2013).
  * [44] A. V. Afanasjev and O. Abdurazakov, Phys. Rev. C 88, 014320 (2013).
  * [45] M. Pfützner, M. Karny, L.V. Grigorenko, K. Riisager, Rev. Mod. Phys. 84, 567 (2012).
  * [46] M. Pfützner, Phys. Scr. T152, 014014 (2013).
  * [47] A.H. Wapstra, G. Audi, and C. Thibault. Nuclear Physics A729, 129 (2003).
  * [48] V.Yu. Denisov, Phys. Lett. B 835, 137569 (2022).


#  Matrix valued positive definite kernels related to the generalized Aitken’s
integral for Gaussians

V. A. Menegatto & C. P. Oliveira

We introduce a method to construct general multivariate positive definite
kernels on a nonempty set $X$ that employs a prescribed bounded completely
monotone function and special multivariate functions on $X$. The method is
consistent with a generalized version of Aitken’s integral formula for
Gaussians. In the case where $X$ is a cartesian product, the method produces
nonseparable positive definite kernels that may be useful in multivariate
interpolation. In addition, it can be interpreted as an abstract multivariate
generalization of the well-established Gneiting’s model for constructing
space-time covariances commonly cited in the literature. Many parametric
models discussed in statistics can be interpreted as particular cases of the
method.

Keywords: multivariate positive definite kernels; conditionally negative
definite functions; Aitken’s formula; Schur exponential; Oppenheim’s
inequality; Gneiting’s model.
2010 MSC: 42A82, 47A56

##  1 Introduction

Let $X$ be a nonempty set and write $M_{q}(\mathbb{C})$ to denote the set of
all $q\times q$ matrices with complex entries. A kernel
$K=[K_{m,n}]_{m,n=1}^{q}:X\times X\to M_{q}(\mathbb{C})$ is positive definite
if for every positive integer $N$ at most the cardinality of $X$ and distinct
points $x_{1},\ldots,x_{N}$ in $X$, the block matrix
$[[K_{m,n}(x_{\mu},x_{\nu})]_{\mu,\nu=1}^{N}]_{m,n=1}^{q}$ of order $Nq$ is
positive semi-definite, that is,

$\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}K(x_{\mu},x_{\nu})c_{\nu}=\sum_{m,n=1}^{q}\sum_{\mu,\nu=1}^{N}\overline{c_{\mu}^{m}}c_{\nu}^{n}K_{m,n}(x_{\mu},x_{\nu})\geq
0,$ (1.1)

whenever $c_{1},\ldots,c_{N}$ are column vectors in $\mathbb{C}^{q}$ and
$c_{\mu}=[c_{\mu}^{1}\,\ldots\,c_{\mu}^{q}]^{\intercal}$. The star notation
refers to conjugate transposition of column vectors in $\mathbb{C}^{q}$. If
the matrices $[[K_{m,n}(x_{\mu},x_{\nu})]_{\mu,\nu=1}^{N}]_{m,n=1}^{q}$ are
all positive definite, that is, the inequalities in (1.1) are strict when at
least one of the vectors $c_{\mu}$ is nonzero, then the positive definite
kernel $K$ is termed strictly positive definite on $X$. The two classes of
kernels introduced above will be denoted by $PD_{q}(X)$ and $SPD_{q}(X)$,
respectively. Kernels in these classes correspond to the standard positive
definite kernels studied in [4] when we set $q=1$ and identify
$M_{q}(\mathbb{C})$ with $\mathbb{C}$. The importance of matrix valued
positive definite kernels in their various formats may be ratified in the
references [2, 16, 17, 23]. Examples of kernels in $PD_{q}(X)$ and
$SPD_{q}(X)$ can be easily constructed. If $A$ is a positive semi-definite
matrix in $M_{q}(\mathbb{C})$, then the constant kernel

$K(x,x^{\prime})=A,\quad x,x^{\prime}\in X,$

belongs to $PD_{q}(X)$. If $f_{1},\ldots,f_{q}$ are kernels in $PD_{1}(X)$,
then the kernel $K$ given by the formula

$K(x,x^{\prime})=\mbox{Diag}(f_{1}(x,x^{\prime}),\ldots,f_{q}(x,x^{\prime})),\quad
x,x^{\prime}\in X,$

belongs to $PD_{q}(X)$. Further, if all the $f_{m}$ belong to $SPD_{1}(X)$,
then $K$ belongs to $SPD_{q}(X)$. Moving the other way around, if $K$ is a
kernel in $PD_{q}(X)$ and $c\in\mathbb{C}^{q}$, then
$f(x,x^{\prime})=c^{*}K(x,x^{\prime})c$, $x,x^{\prime}\in X$, defines a
function in $PD_{1}(X)$. If $c\neq 0$ and $K$ belongs to $SPD_{q}(X)$, then
$f$ actually belongs to $SPD_{1}(X)$.

The purpose of this paper is to introduce methods to construct abstract
matrix-valued mappings with the additional requirement of positive
definiteness and strict positive definiteness. In many cases, the method
yields flexible models, once it encompasses models found in geophysical
sciences, including probabilistic weather forecasting, data assimilation,
statistical analysis of climate model output, etc, when one makes the right
choice for $X$ and set a metric structure in it.

The method itself will be based on bounded completely monotone functions and
special matrix valued functions attached to the notion of conditional negative
definiteness. Recall that the complete monotonicity of a function
$f:(0,\infty)\to\mathbb{R}$ is characterized by two properties: $f$ is
$C^{\infty}$ and $(-1)^{n}f^{(n)}(t)\geq 0$ for $n=0,1,\ldots$ and
$t\in(0,\infty)$. Throughout the paper, we will not distinguish between a
bounded completely monotone function and its unique continuous extension to
$[0,\infty)$. A kernel $K=[K_{m,n}]_{m,n=1}^{q}:X\times X\to
M_{q}(\mathbb{C})$ is conditionally negative definite if it is Hermitian and
the matrices $[[K_{m,n}(x_{\mu},x_{\nu})]_{\mu,\nu=1}^{N}]_{m,n=1}^{q}$ are of
negative type, that is, the quadratic forms (1.1) are nonpositive whenever the
vectors $c_{\mu}$ satisfy $\sum_{\mu=1}^{N}c_{\mu}=0$. The conditionally
negative definite kernel $K$ is strictly conditionally negative definite if
the matrices $[[K_{m,n}(x_{\mu},x_{\nu})]_{\mu,\nu=1}^{N}]_{m,n=1}^{q}$ are of
strict negative type for $N\geq 2$, that is, the quadratic forms are negative
whenever $N\geq 2$ and at least one $c_{\mu}$ is nonzero. These two classes of
kernels will be denoted by $CND_{q}(X)$ and $SCND_{q}(X)$, respectively.
Examples of kernels in $CND_{1}(X)$ and $SCND_{1}(X)$ can be found in [4]
while connections between the classes $PD_{1}(X)$ and $CND_{1}(X)$ are
described in [3, 4, 10]. As for examples in the classes $CND_{q}(X)$ and
$SCND_{q}(X)$ one may imitate the procedures adopted for producing kernels in
$PD_{q}(X)$ and $SPD_{q}(X)$ previously mentioned.

All the major results we intend to prove here will be based on a
generalization of Aitken’s integral formula for computing Gaussians: if $A$ is
a positive definite matrix in $M_{q}(\mathbb{R})$ and $b$ is a vector in
$\mathbb{R}^{q}$, then

$\int_{\mathbb{R}^{q}}e^{\displaystyle{-u^{\intercal}Au+i\,b^{\intercal}u}}du=\frac{\pi^{q/2}}{\sqrt{\det\,A}}e^{\displaystyle{-b^{\intercal}(4A)^{-1}b}}.$

Aitken’s integral itself corresponds to the formula above in the case $b=0$. A
proof for the generalized Aitken’s integral formula can be reached by
mimicking the proof of Aitken’s integral in [22, p. 340] but an independent
proof is available in [14]. This reference also contains univariate versions
of some of the results to be described here.

Before we proceed to the outline of the paper, it is worth mentioning that if
$X$ is actually a cartesian product of sets, the method to be presented here
lead to nonseparable kernels, a desirable property in applications. Meanwhile,
in some specific cases, the method will become a generalization of the well
established Gneiting’s contribution in [7] on the construction of kernels in
$PD_{1}(\mathbb{R}^{q}\times\mathbb{R}^{d})$. Gneiting’s classical result is
as follows: for a bounded completely monotone function
$\phi:(0,\infty)\to\mathbb{R}$ and a positive valued function $f$ with a
completely monotone derivative, it asserts that the formula

$G_{r}((x,y),(x^{\prime},y^{\prime}))=\frac{1}{f(\|y-y^{\prime}\|^{2})^{r}}\phi\left(\frac{\|x-x^{\prime}\|^{2}}{f(\|y-y^{\prime}\|^{2})}\right),\quad
x,x^{\prime}\in\mathbb{R}^{q};\,y,y^{\prime}\in\mathbb{R}^{d},$ (1.2)

defines a kernel $G_{r}$ in $PD_{1}(\mathbb{R}^{q}\times\mathbb{R}^{d})$,
whenever $r\geq d/2$ and $\|\cdot\|$ denotes the usual norms in both
$\mathbb{R}^{q}$ and $\mathbb{R}^{d}$. The boundedness of $\phi$ is required
in order to make $\phi(0^{+})<\infty$. The references [12, 18] include some
extensions and generalizations of this important result along with additional
references on the topic.

The paper proceeds as follows. Section 2 begins with the description of two
additional notions to be employed in the paper, one for families of vector
functions and another for families of matrix functions, along with examples.
The first major result of the paper is Theorem 2.4: it describes a method to
construct kernels in $PD_{p}(Y)$ from bounded completely monotone functions,
special families of vector functions on $Y$ and special families of matrix
functions on $Y$. Further, it provides a sufficient condition in order that
the resulting kernel be in $SPD_{p}(Y)$. At the end of the section we discuss
some examples and detach a relevant consequence of Theorem 2.4. The main
result in Section 3 expands Theorem 2.4 via integration with respect to a
convenient measure. We separate a special simpler version of the theorem in
Corollary 3.3. Section 4 describes extensions of Theorems 2.4 and 3.1 that
lead to kernels in $PD_{p}(X\times Y)$. Applications and a multivariate
abstract extension of the classical Gneiting’s result are described.

##  2 The main result for positive definiteness on a single set

This section contains the first main contribution in the paper to be made
explicit in Theorem 2.4. It provides a method to construct functions in
$PD_{q}(Y)$ using completely monotonic functions via Aitken’s formula. A
sufficient condition for strict positive definiteness is included. The
contribution itself demands two notions for families of functions with domain
$Y$ which we now discuss.

For a matrix function $G$ in $CND_{q}(Y)$ and a vector $u$ from
$\mathbb{C}^{q}$, the kernel

$(y,y^{\prime})\in Y\times Y\mapsto u^{*}G(y,y^{\prime})u$

belongs to $CND_{1}(Y)$. Further, the kernel belongs to $SCND_{1}(Y)$ whenever
$G$ belongs to $SCND_{q}(Y)$ and $u$ is nonzero. Theorem 2.4 will demand a
family $\\{G_{m,n}:m,n=1,\ldots,p\\}$ for which all the matrix kernels

$(y,y^{\prime})\in Y\times
Y\mapsto\left[u^{\intercal}G_{m,n}(y,y^{\prime})u\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $CND_{p}(Y)$. Since this is not easily achievable, the following
example is apposite.

######  Example 2.1.

Define

$G_{m,n}(y,y^{\prime})=g_{m}(y)+g_{n}(y^{\prime}),\quad y,y^{\prime}\in X,$

where the $g_{m}:Y\to M_{q}(\mathbb{C})$ are functions subject to our choice.
If $y_{1},\ldots,y_{N}$ are distinct points in $Y$, $c_{1},\ldots,c_{N}$ are
vectors in $\mathbb{C}^{p}$ such that $\sum_{\mu=1}^{N}c_{\mu}=0$, and
$u\in\mathbb{C}^{q}$, then

$\displaystyle\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[u^{\intercal}G_{m,n}(y_{\mu},y_{\nu})u\right]_{m,n=1}^{p}c_{\nu}$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{p}\sum_{\nu=1}^{N}c_{\nu}^{n}\sum_{\mu=1}^{N}\sum_{m=1}^{p}\overline{c_{\mu}^{m}}u^{\intercal}g_{m}(y_{\mu})u$
$\displaystyle\hskip
8.53581pt+\sum_{m=1}^{p}\sum_{\mu=1}^{N}\overline{c_{\mu}^{m}}\sum_{\nu=1}^{N}\sum_{n=1}^{p}c_{\nu}^{n}u^{\intercal}g_{n}(y_{\nu})u=0,$

that is, the matrix function

$(y,y^{\prime})\in X\times
X\mapsto\left[u^{\intercal}G_{m,n}(y,y^{\prime})u\right]_{m,n=1}^{p},$

belongs to $CND_{p}(Y)$.

######  Example 2.2.

Set $G_{m,n}=0$ when $m\neq n$ and pick each $G_{m,m}$ in the class
$CND_{q}(Y)$. Keeping the $c_{\mu}$ and the $y_{\mu}$ as in Example 2.1, it is
easily seen that

$\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[u^{\intercal}G_{m,n}(y_{\mu},y_{\nu})u\right]_{m,n=1}^{p}c_{\nu}=\sum_{m=1}^{p}\sum_{\mu,\nu=1}^{N}\overline{c_{\mu}^{m}}c_{\nu}^{m}u^{\intercal}G_{m,m}(y_{\mu},y_{\nu})u\leq
0.$

Thus, the matrix function

$(y,y^{\prime})\in X\times
X\mapsto\left[u^{\intercal}G_{m,n}(y,y^{\prime})u\right]_{m,n=1}^{p},$

belongs to $CND_{p}(Y)$.

Theorem 2.4 will also need special families $\\{H_{m,n}:m,n=1,\ldots,p\\}$ of
vector functions $H_{m,n}:Y\times Y\to\mathbb{C}^{q}$. As matter of fact, we
will require families for which all the matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,H_{m,n}(y,y^{\prime})^{*}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $PD_{p}(Y)$. Again, this is not easy to achieve, reason why a simple
example is handy.

######  Example 2.3.

Let us set

$H_{m,n}(y,y^{\prime})=h_{m}(y)-h_{n}(y^{\prime}),\quad y,y^{\prime}\in Y,$

where $h_{m}:Y\to\mathbb{R}^{q}$, $m=1,\ldots,p$. If $y_{1},\ldots,y_{N}$ are
distinct points in $Y$ and $c_{1},\ldots,c_{N}$ are vectors in
$\mathbb{C}^{p}$, then

$\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[e^{\displaystyle{i\,H_{m,n}(y_{\mu},y_{\nu})^{\intercal}u}}\right]_{m,n=1}^{p}c_{\nu}=\left|\sum_{\mu=1}^{N}\sum_{m=1}^{p}\overline{c_{\mu}^{m}}\,e^{\displaystyle{i\,h_{m}(y_{\mu})^{\intercal}u}}\right|^{2}\geq
0,\quad u\in\mathbb{R}^{q},$

that is, the kernels

$(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,H_{m,n}(y,y^{\prime})^{*}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $PD_{p}(Y)$.

We observe that if the matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,H_{m,n}(y,y^{\prime})^{*}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $PD_{p}(Y)$, then each $H_{m,n}$ must be anti-symmetric in the sense
that

$\mbox{Re\,}H_{m,n}(y,y^{\prime})=-\mbox{Re\,}H_{m,n}(y^{\prime},y),\quad
y,y^{\prime}\in Y.$

In particular,

$\mbox{Re\,}H_{m,n}(y,y)=0,\quad m,n=1,\ldots,p;\,y\in Y.$

Some specific properties of Hadamard exponentials will be needed. We recall
that if $A$ is a matrix in $M_{q}(\mathbb{C})$, then its Hadamard exponential
is the matrix

$e^{\circ A}:=[e^{A_{\mu\nu}}]_{\mu,\nu=1}^{q}.$

Let $M_{q}(\mathbb{R})$ denote subset of $M_{q}(\mathbb{C})$ composed of real
matrices only. If $A\in M_{q}(\mathbb{R})$ is symmetric and of negative type,
then the Hadamard exponential of $-A$ is positive semi-definite. It is
positive definite if, and only if,

$A_{\mu\mu}+A_{\nu\nu}<2A_{\mu\nu},\quad\mu\neq\nu.$

These facts are proved in Lemma 2.5 in [19] albeit [13] analyzed similar
properties earlier. As an obvious consequence, we have that if $A\in
M_{q}(\mathbb{R})$ is of strict negative type, then the Hadamard exponential
of $-A$ is positive definite. Finally, if a real symmetric block matrix
$A=[[A_{mn}(\mu\nu)]_{\mu,\nu=1}^{N}]_{m,n=1}^{q}$ is of negative type, then
the Hadamard exponential of $-A$ is positive definite if, and only if,

$A_{mm}(\mu\mu)+A_{nn}(\nu\nu)<2A_{mn}(\mu\nu),\quad|m-n|+|\mu-\nu|>0.$ (2.3)

Below, we will use the symbol $\bullet$ to denote the Schur product of two
matrices of same size.

######  Theorem 2.4.

Let $\phi$ be a bounded and completely monotone function. For each $m,n$ in
$\\{1,\ldots,p\\}$, let $G_{m,n}:Y\times Y\to M_{q}(\mathbb{R})$ be a matrix
function with range containing positive definite matrices only and
$H_{m,n}:Y\times Y\to\mathbb{R}^{q}$ a vector function. Assume the matrix
functions

$(y,y^{\prime})\in Y\times
Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $CND_{p}(Y)$ and that

$(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,H_{m,n}(y,y^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $PD_{p}(Y)$. The following assertions hold for the kernel $K:Y\times
Y\to M_{p}(\mathbb{R})$ given by the formula

$K(y,y^{\prime})=\left[\frac{\phi\left(H_{m,n}(y,y^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(y,y^{\prime})\right)}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}\right]_{m,n=1}^{p},\quad y,y^{\prime}\in Y.$

  * $(i)$

$K$ belongs to $PD_{p}(Y)$.

  * $(ii)$

If $\phi$ is not identically 0 and there exists an open subset $U$ of
$\mathbb{R}^{q}\setminus\\{0\\}$ so that

$u^{\intercal}[G_{m,m}(y,y)+G_{n,n}(y^{\prime},y^{\prime})-2G_{m,n}(y,y^{\prime})]u<0,\quad(m,y)\neq(n,y^{\prime});\,u\in
U,$

then $K$ belongs to $SPD_{p}(Y)$.

###### Proof.

We begin proving Assertion $(i)$ in the case where $\phi$ is a constant
function, that is, the case in which

$K(y,y^{\prime})=\left[\frac{\phi(0)}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}\right]_{m,n=1}^{p},\quad y,y^{\prime}\in Y.$

Since each matrix $G_{m,n}(y,y^{\prime})$ is positive definite, we may apply
Aitken’s integral formula to obtain

$K(y,y^{\prime})=\frac{\phi(0)}{\pi^{q/2}}\left[\int_{\mathbb{R}^{q}}e^{\displaystyle{-u^{\intercal}G_{m,n}(y,y^{\prime})u}}du\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y.$ (2.4)

If $y_{1},\ldots,y_{N}$ are distinct points in $Y$ and $c_{1},\ldots,c_{N}$
are vectors in $\mathbb{R}^{p}$, then

$\displaystyle\sum_{\mu,\nu=1}^{N}c_{\mu}^{\intercal}K(y_{\mu},y_{\nu})c_{\nu}$
$\displaystyle=$
$\displaystyle\frac{\phi(0)}{\pi^{q/2}}\sum_{\mu,\nu=1}^{N}\sum_{m,n=1}^{p}c_{\mu}^{m}c_{\nu}^{n}\int_{\mathbb{R}^{q}}e^{\displaystyle{-u^{\intercal}G_{m,n}(y_{\mu},y_{\nu})u}}du$
$\displaystyle=$
$\displaystyle\frac{\phi(0)}{\pi^{q/2}}\int_{\mathbb{R}^{q}}\sum_{\mu,\nu=1}^{N}c_{\mu}^{\intercal}E_{u}(y_{\mu},y_{\nu})c_{\nu}du,$

where

$E_{u}(y,y^{\prime})=\left[e^{\displaystyle{-u^{\intercal}G_{m,n}(y,y^{\prime})u}}\right]_{m,n=1}^{p}=e^{\displaystyle{\circ[-u^{\intercal}G_{m,n}(y,y^{\prime})u}]_{m,n=1}^{p}},\quad
u\in\mathbb{R}^{q}.$

One of the assumptions on the $G_{m,n}$ now yields that

$\sum_{\mu,\nu=1}^{N}c_{\mu}^{\intercal}E_{u}(y_{\mu},y_{\nu})c_{\nu}\geq
0,\quad u\in\mathbb{R}^{q},$

and Assertion $(i)$ follows in this case. In the general case, the Bernstein-
Widder Theorem ([21, p. 3]) implies that

$K(y,y^{\prime})=\left[\frac{1}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}\int_{[0,\infty)}e^{\displaystyle{-H_{m,n}(y,y^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(y,y^{\prime})\,s}}d\sigma(s)\right]_{m,n=1}^{p}$

for some finite and positive measure $\sigma$ on $[0,\infty)$. On the other
hand, the generalized Aitken’s integral formula provides the alternative
representation

$K(y,y^{\prime})=\left[\frac{1}{\pi^{q/2}}\int_{[0,\infty)}\left(\int_{\mathbb{R}^{q}}e^{\displaystyle{-u^{\intercal}G_{m,n}(y,y^{\prime})u}}e^{\displaystyle{2i\,\sqrt{s}H_{m,n}(y,y^{\prime})^{\intercal}u}}du\right)d\sigma(s)\right]_{m,n=1}^{p}.$

If the $y_{\mu}$ are as before and the $c_{\mu}$ are now complex vectors, the
quadratic form

$Q:=\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}K(y_{\mu},y_{\nu})c_{\nu}$

becomes

$\displaystyle Q$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{q/2}}\sum_{\mu,\nu=1}^{N}\sum_{m,n=1}^{p}\overline{c_{\mu}^{m}}c_{\nu}^{n}\int_{[0,\infty)}\int_{\mathbb{R}^{q}}e^{-\displaystyle{u^{\intercal}G_{m,n}(y_{\mu},y_{\nu})u}}e^{\displaystyle{i\,2\sqrt{s}H_{m,n}(y_{\mu},y_{\nu})^{\intercal}u}}dud\sigma(s)$
$\displaystyle=$
$\displaystyle\frac{1}{\pi^{q/2}}\int_{[0,\infty)}\int_{\mathbb{R}^{q}}\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[E_{u}(y_{\mu},y_{\nu})\bullet
E_{u}^{s}(y_{\mu},y_{\nu})\right]c_{\nu}\,dud\sigma(s),$

where

$E_{u}^{s}(y,y^{\prime})=\left[e^{\displaystyle{i\,2\sqrt{s}H_{m,n}(y,y^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y;s\geq 0.$

The assumption on the $H_{m,n}$ settles the positive semi-definiteness of each
matrix $E_{u}^{s}(y_{\mu},y_{\nu})$ while the Schur Product Theorem ratifies
the positive semi-definiteness of each Schur product
$E_{u}(y_{\mu},y_{\nu})\bullet E_{u}^{s}(y_{\mu},y_{\nu})$. These arguments
validate the inequality $Q\geq 0$.
Let us keep the notation used above to prove Assertion $(ii)$. Assume further
that the $c_{\mu}$ are not all zero vectors. If there exists an open subset
$U$ of $\mathbb{R}^{q}\setminus\\{0\\}$ so that

$u^{\intercal}[G_{m,m}(y,y)+G_{n,n}(y^{\prime},y^{\prime})-2G_{m,n}(y,y^{\prime})]u<0,\quad(m,y)\neq(n,y^{\prime});\,u\in
U,$

we can infer via (2.3) that the block matrix

$E_{u}(y_{\mu},y_{\nu})=\left[\left[e^{\displaystyle{-u^{\intercal}G_{m,n}(y_{\mu},y_{\nu})u}}\right]_{\mu,\nu=1}^{N}\right]_{m,n=1}^{p}$

is positive definite whenever $u\in U$. Thus, if $\phi$ is constant and not
identically 0, then $Q>0$ by Formula (2.4). If $\phi$ is nonconstant, first we
invoke our assumption on the $H_{m,n}$ in order to see that the diagonal
entries in each block matrix

$E_{u}^{s}(y_{\mu},y_{\nu})=\left[\left[e^{\displaystyle{i\,2\sqrt{s}H_{m,n}(y_{\mu},y_{\nu})^{\intercal}u}}\right]_{\mu,\nu=1}^{N}\right]_{m,n=1}^{p}$

are all equal to 1. An application of Oppenheim’s inequality ([9, p. 509])
shows that the Schur product $E_{u}(y_{\mu},y_{\nu})\bullet
E_{u}^{s}(y_{\mu},y_{\nu})$ is positive definite for $u\in U$ and $s\geq 0$.
In particular,

$\int_{\mathbb{R}^{q}}\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[E_{u}(y_{\mu},y_{\nu})\bullet
E_{u}^{s}(y_{\mu},y_{\nu})\right]c_{\nu}\,du>0,\quad s\geq 0.$

Since $\sigma$ is not the zero measure we may go one step further and infer
that $Q>0$. ∎

######  Remark 2.5.

Theorem 17 in [20] is a very special case of Theorem 2.4-$(i)$.

Next, we present some examples that illustrate our findings.

######  Example 2.6.

For $m=1,\ldots,p$, let $g_{m}:Y\to M_{q}(\mathbb{R})$ be a function with
range containing positive definite matrices only and
$h_{m}:Y\to\mathbb{R}^{q}$ an arbitrary function. Setting
$G_{m,n}(y,y^{\prime})=g_{m}(y)+g_{n}(y^{\prime})$, $y,y^{\prime}\in Y$, and
$H_{m,n}(y,y^{\prime})=h_{m}(y)-h_{n}(y^{\prime})$, $y,y^{\prime}\in Y$, the
assumptions in Theorem 2.4 are satisfied. Thus, the formula

$\left[\frac{\phi\left((h_{m}(y)-h_{n}(y^{\prime}))^{\intercal}(g_{m}(y)+g_{n}(y^{\prime}))^{-1}(h_{m}(y)-h_{n}(y^{\prime}))\right)}{\sqrt{\det[(g_{m}(y)+g_{n}(y^{\prime})]}}\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y,$

defines a kernel in $PD_{p}(Y)$ whenever $\phi$ is bounded completely monotone
function. The inequalities in Theorem 2.4-$(ii)$ cannot be matched in this
abstract example.

######  Example 2.7.

For $m,n=1,\ldots,p$, let us set

$G_{m,n}(y,y^{\prime})=g_{m,n}(y,y^{\prime})I_{q},\quad y,y^{\prime}\in Y,$

where each $g_{m,n}$ is a positive valued kernel on $Y$ and $(y,y^{\prime})\in
Y\times Y\mapsto[g_{m,n}(y,y^{\prime})]_{m,n=1}^{p}$ belongs to $CND_{p}(Y)$.
Observe that for each $m$ and $n$,

$u^{\intercal}G_{m,n}(y,y^{\prime})u=\|u\|^{2}g_{m,n}(y,y^{\prime}),\quad
u\in\mathbb{R}^{q};y,y^{\prime}\in Y.$

On the other hand, if $c_{1},\ldots,c_{N}$ are column vectors satisfying
$\sum_{\mu=1}^{n}c_{\mu}=0$ and $y_{1},\ldots,y_{n}$ belong to $Y$, then

$\sum_{m,n=1}^{p}\sum_{\mu,\nu=1}^{N}\overline{c_{\mu}^{m}}c_{\nu}^{n}u^{\intercal}G_{m,n}(y_{\mu},y_{\nu})u=\|u\|^{2}\sum_{m,n=1}^{p}\sum_{\mu,\nu=1}^{N}\overline{c_{\mu}^{m}}c_{\nu}^{n}g_{m,n}(y_{\mu},y_{\nu})\leq
0.$

that is, each kernel

$(y,y^{\prime})\in Y\times
Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},$

belongs to $CND_{p}(Y)$. If the $H_{m,n}$ satisfy the assumptions of Theorem
2.4, then it is promptly seen that the formula

$K(y,y^{\prime})=\left[\frac{1}{g_{m,n}(y,y^{\prime})^{q/2}}\phi\left(\frac{\|H_{m,n}(y,y^{\prime})\|^{2}}{g_{m,n}(y,y^{\prime})}\right)\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y,$

defines a matrix kernel in $PD_{p}(Y)$ whenever $\phi$ is a bounded completely
monotone function.

######  Example 2.8.

If we take $H_{m,n}$ as in Example 2.6, then the kernel $K$ in Example 2.7
takes the form

$K(y,y^{\prime})=\left[\frac{1}{g_{m,n}(y,y^{\prime})^{q/2}}\phi\left(\frac{\|h_{m}(y)-h_{n}(y^{\prime})\|^{2}}{g_{m,n}(y,y^{\prime})}\right)\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y.$

This example has an structure that resembles that of Gneiting’s model in [7]
for the construction of space-time covariances. We can get even closer by
setting $g_{m,n}:=g$ for all $m$ and $n$, where $g:Y\to(0,\infty)$ belongs to
$CND_{1}(Y)$, a choice that leads to

$K(y,y^{\prime})=\frac{1}{g(y,y^{\prime})^{q/2}}\left[\phi\left(\frac{\|h_{m}(y)-h_{n}(y^{\prime})\|^{2}}{g(y,y^{\prime})}\right)\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y.$

The setting adopted in both Examples 2.7 and 2.8 is a particular case of that
detached in Theorem 2.9 below. Needless to say that the theorem can be
interpreted as a multivariate version of the Gneiting’s criterion in [7].

######  Theorem 2.9.

Let $\phi$ be a bounded and completely monotone function. Let $g$ be a
positive valued kernel in $CND_{1}(Y)$ and for each $m,n$ in
$\\{1,\ldots,p\\}$, define

$G_{m,n}(y,y^{\prime})=g(y,y^{\prime})I_{q},\quad y,y^{\prime}\in Y.$

If $H_{m,n}:Y\times Y\to\mathbb{R}^{q}$ is a vector function such that the
matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,H_{m,n}(y,y^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $PD_{p}(Y)$, then the following assertions hold for the kernel
$K:Y\times Y\to M_{p}(\mathbb{R})$ given by the formula

$K(y,y^{\prime})=\frac{1}{g(y,y^{\prime})^{q/2}}\left[\phi\left(\frac{\|H_{m,n}(y,y^{\prime})\|^{2}}{g(y,y^{\prime})}\right)\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y.$

  * $(i)$

$K$ belongs to $PD_{p}(Y)$.

  * $(ii)$

If $\phi$ is not identically 0 and
$g(y,y)+g(y^{\prime},y^{\prime})-2g(y,y^{\prime})<0$ for $y\neq y^{\prime}$,
then $K$ belongs to $SPD_{p}(Y)$.

##  3 An extension of the main result via integration

Here, we extend the results proved in Section 2 by introducing a scale mixture
in the formula that defines the positive definite kernels.

Our first contribution here is as follows.

######  Theorem 3.1.

Let $\rho$ be a nonzero positive measure on $(0,\infty)$ and $\phi$ a bounded
and completely monotone function. For each $m,n$ in $\\{1,\ldots,p\\}$, let
$G_{m,n}:Y\times Y\to M_{q}(\mathbb{R})$ be a matrix function with range
containing positive definite matrices only, $H_{m,n}:Y\times
Y\to\mathbb{R}^{q}$ a vector function and $\\{P_{m,n}^{s}\\}_{s>0}$ a family
of kernels on $Y$ such that each function $s\in(0,\infty)\mapsto
P_{m,n}^{s}(y,y^{\prime})$ is $\rho$-integrable. If the matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$ $(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,H_{m,n}(y,y^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

and

$(y,y^{\prime})\in Y\times
Y\mapsto[P_{m,n}^{s}(y,y^{\prime})]_{m,n=1}^{p},\quad s>0,$

belong to $CND_{p}(Y)$, $PD_{p}(Y)$, and $PD_{p}(Y)$, respectively, then the
kernel $K:Y\times Y\to M_{p}(\mathbb{R})$ given by the formula

$\displaystyle K(y,y^{\prime})$ $\displaystyle=$
$\displaystyle\left[\frac{1}{\sqrt{\det G_{m,n}(y,y^{\prime})}}\right.$
$\displaystyle\times$
$\displaystyle\left.\int_{(0,\infty)}\phi\left(H_{m,n}(y,y^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(y,y^{\prime})\,s\right)P_{m,n}^{s}(y,y^{\prime})d\rho(s)\right]_{m,n=1}^{p}$

belongs to $PD_{p}(Y)$.

###### Proof.

Let $y_{1},\ldots,y_{N}$ be distinct points in $Y$, $c_{1},\ldots,c_{N}$
vectors in $\mathbb{C}^{p}$ and set
$Q:=\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}K(y_{\mu},y_{\nu})c_{\nu}$. Direct
calculation shows that

$Q=\int_{(0,\infty)}\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[I^{s}(y_{\mu},y_{\nu})\bullet
P^{s}(y_{\mu},y_{\nu})\right]c_{\nu}d\rho(s),$

where

$I^{s}(y,y^{\prime})=\left[\frac{\phi\left(\sqrt{s}H_{m,n}(y,y^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}\sqrt{s}H_{m,n}(y,y^{\prime})\right)}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}\right]_{m,n=1}^{p}$

and

$P^{s}(y,y^{\prime})=\left[P_{m,n}^{s}(y,y^{\prime})\right]_{m,n=1}^{p},\quad
y,y^{\prime}\in Y;s>0.$

As in the proof of Theorem 2.4, the matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,\sqrt{s}H_{m,n}(y,y^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q};s>0,$

belong to $PD_{p}(Y)$. However, since the assumptions on the $G_{m,n}$ are the
same as those in Theorem 2.4, we can apply Theorem 2.4-$(i)$ in order to see
that each matrix $I^{s}(y_{\mu},y_{\nu})$ is positive semi-definite. As for
$P^{s}(y_{\mu},y_{\nu})$, $s>0$, they are positive semi-definite as well by
our assumption on the family $\\{P_{m,n}^{s}\\}_{s>0}$. Thus, the Schur
Product Theorem implies that

$\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[I^{s}(y_{\mu},y_{\nu})\bullet
P^{s}(y_{\mu},y_{\nu})\right]c_{\nu}\geq 0,\quad s>0.$

Therefore, $Q\geq 0$. ∎

As for strict positive definiteness, the following consequence of Theorem 3.1
holds.

######  Theorem 3.2.

If $\phi$ is not identically zero, then the following additional assertions
hold for the kernel $K$ defined in Theorem 3.1:

  * $(i)$

If there exists an open subset $A$ of $\mathbb{R}^{q}\setminus\\{0\\}$ so that

$u^{\intercal}[G_{m,m}(y,y)+G_{n,n}(y^{\prime},y^{\prime})-2G_{m,n}(y,y^{\prime})]u<0,\quad(m,y)\neq(n,y^{\prime});\,u\in
U,$

and a $\rho$-measurable subset $A$ of $(0,\infty)$ so that $\rho(A)>0$ and

$P_{m,m}^{s}(y,y)>0,\quad m\in\\{1,\ldots,p\\};\,y\in Y;\,s\in A,$

then $K$ belongs to $SPD_{p}(Y)$.

  * $(ii)$

If there exists a $\rho$-measurable subset $A$ of $(0,\infty)$ so that
$\rho(A)>0$ and

$(y,y^{\prime})\in Y\times Y\mapsto[P_{m,n}^{s}(y,y^{\prime})]_{m,n=1}^{p}\in
SPD_{p}(Y),\quad s\in A,$

then $K$ belongs to $SPD_{p}(Y)$.

###### Proof.

Let the $x_{\mu}$ and the $c_{\mu}$ be as in the proof of Theorem 3.1.
Further, assume at least one $c_{\mu}$ is nonzero. If the assumptions in $(i)$
hold, then Theorem 2.4-$(ii)$ implies that each matrix
$I^{s}(y_{\mu},y_{\nu})$ is positive definite while the diagonal entries in
$P^{s}(y_{\mu},y_{\nu})$ are all positive for $s\in A$. Therefore, by
Oppenheim’s inequality, we can assert that

$\sum_{\mu,\nu=1}^{N}c_{\mu}^{*}\left[I^{s}(y_{\mu},y_{\nu})\bullet
P^{s}(y_{\mu},y_{\nu})\right]c_{\nu}>0,\quad s\in A.$

Since the measure $\rho$ is nonzero, $Q>0$. If the assumptions in $(ii)$ hold,
we may reach the very same conclusion once the diagonal elements in
$I^{s}(y_{\mu},y_{\nu})$, $s>0$, are given by

$\frac{\phi(0)}{\sqrt{\det G_{m,m}(y_{\mu},y_{\mu})}}>0,\quad
m=1,\ldots,p;\mu=1,\dots,q.$

Indeed, Oppenheim’s inequality once again would imply that $Q>0$. ∎

A specially chosen family $\\{G_{m,n}:m,n=1,\ldots,p\\}$ in Theorem 3.1 leads
to the following improved abstract multivariate version of Gneiting’s
criterion in [7].

######  Corollary 3.3.

Let $\phi:(0,\infty)\to\mathbb{R}$ be a bounded and completely monotone
function. For $m,n=1,2,\ldots,p$, set
$G_{m,n}(y,y^{\prime})=g_{m,n}(y,y^{\prime})I_{q}$, $y,y^{\prime}\in Y$, where
each $g_{m,n}$ is a positive valued kernel in $CND_{1}(Y)$, let
$H_{m,n}:Y\times Y\to\mathbb{R}^{q}$ be a vector function and
$\\{P_{m,n}^{s}\\}_{s>0}$ a family of kernels on $Y$ such that each function
$s\in(0,\infty)\mapsto P_{m,n}^{s}(y,y^{\prime})$ is $\rho$-integrable. If the
matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$ $(y,y^{\prime})\in Y\times
Y\mapsto\left[e^{\displaystyle{i\,H_{m,n}(y,y^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

and

$(y,y^{\prime})\in Y\times
Y\mapsto[P_{m,n}^{s}(y,y^{\prime})]_{m,n=1}^{p},\quad s>0,$

belong to $CND_{p}(Y)$, $PD_{p}(Y)$, and $PD_{p}(Y)$, respectively, then the
kernel $K:Y\times Y\to M_{p}(\mathbb{R})$ given by the formula

$K(y,y^{\prime})=\left[\frac{1}{g_{m,m}(y,y^{\prime})^{q/2}}\int_{0}^{\infty}\phi\left(\frac{\|H_{m,n}(y,y^{\prime})\|^{2}s}{g_{m,n}(y,y^{\prime})}\right)P_{m,n}^{s}(y,y^{\prime})d\rho(s)\right]_{m,n=1}^{p}$

belongs to $PD_{p}(Y)$. Further, if $\phi$ is not identically 0, the following
two additional assertions hold:

  * $(i)$

If $g_{m,m}(y,y)+g_{n,n}(y^{\prime},y^{\prime})-2g_{m,n}(y,y^{\prime})<0$ when
$(m,y)\neq(n,y^{\prime})$, and there exists a $\rho$-measurable subset $A$ of
$(0,\infty)$ so that $\rho(A)>0$ and

$P_{m,m}^{s}(y,y)>0,\quad m\in\\{1,\ldots,p\\};\,y\in Y;\,s\in A,$

then $K$ belongs to $SPD_{p}(Y)$.

  * $(ii)$

If there exists a $\rho$-measurable subset $A$ of $(0,\infty)$ so that
$\rho(A)>0$ and

$(y,y^{\prime})\in Y\times Y\mapsto[P_{m,n}^{s}(y,y^{\prime})]_{m,n=1}^{p},$

belongs to $SPD_{p}(Y)$ for $s\in A$, then $K$ belongs to $SPD_{p}(Y)$.

##  4 The main results in the case of a product of sets

An easy way to construct kernels in $PD_{p}(X\times Y)$ is given by the
product of a kernel in $PD_{p}(X)$ with another one in $PD_{p}(Y)$, a fact
that can be ratified via the Schur Product Theorem. The separable kernels
produced by this method may be not suitable if one needs strong interactions
between $X$ and $Y$. The main result in this section will provide a version of
Theorem 2.4 that leads to kernels in $PD_{p}(X\times Y)$ and, except for very
particular cases, the kernels produced by this version will be nonseparable.
In particular, the aforementioned interactions are possible. The result
explains, from a mathematical point of view, some important practical models
adopted in the statistical literature. The proofs will be omitted once they
are very similar to those of the theorems proved in Sections 2 and 3.

######  Theorem 4.1.

Let $\phi:(0,\infty)\to\mathbb{R}$ be a bounded and completely monotone
function. For each $m,n$ in $\\{1,\ldots,p\\}$, let $G_{m,n}:Y\times Y\to
M_{q}(\mathbb{R})$ be a matrix function with range containing positive
definite matrices only and $H_{m,n}:X\times X\to\mathbb{R}^{q}$ a vector
function. If the matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $CND_{p}(Y)$ and

$(x,x^{\prime})\in X\times
X\mapsto\left[e^{\displaystyle{i\,H_{m,n}(x,x^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $PD_{p}(X)$, then the kernel $K:(X\times Y)^{2}\to
M_{p}(\mathbb{R})$ given by

$K((x,y),(x^{\prime},y^{\prime}))=\left[\frac{\phi\left(H_{m,n}(x,x^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(x,x^{\prime})\right)}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}\right]_{m,n=1}^{p}$

belongs to $PD_{p}(X\times Y)$.

In the Example LABEL:exa below, we illustrate Theorem 4.1 in the case
$X=\mathbb{R}$ and $Y=S^{d}$, the unit sphere in $\mathbb{R}^{d+1}$.

######  Example 4.2.

Define $H_{m,n}(x,x^{\prime})=h_{m}(x)-h_{n}(x^{\prime})$,
$x,x^{\prime}\in\mathbb{R}$, where each $h_{m}:\mathbb{R}\to\mathbb{R}^{q}$ is
an arbitrary function. If $\delta$ denotes the geodesic distance in $S^{d}$,
set

$G_{m,n}(y,y^{\prime})=[m+n+\delta(y,y^{\prime})]I_{q},\quad y,y^{\prime}\in
S^{d}.$

It is well known that $(y,y^{\prime})\in S^{d}\times
S^{d}\mapsto\delta(y,y^{\prime})$ belongs to $CND_{1}(S^{d})$ (see Section 4
in [1]). Hence, each $G_{m,n}$ has range containing positive definite matrices
only. On the other hand, according to Examples 2.7 and 2.8, each kernel

$(y,y^{\prime})\in Y\times
S^{d}\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belongs to $CND_{p}(S^{d})$. It follows that

$K((x,y),(x^{\prime},y^{\prime}))=\left[\frac{1}{[m+n+\delta(y,y^{\prime})]^{q/2}}\phi\left(\frac{\|h_{m}(x)-h_{n}(x^{\prime})\|^{2}}{m+n+\delta(y,y^{\prime})}\right)\right]_{m,n=1}^{p}$

belongs to $PD_{p}(\mathbb{R}\times S^{d})$. The choice

$h_{m}(x)=(x,0,\ldots,0)^{\intercal},\quad x\in\mathbb{R};\,m=1,\ldots,p,$

leads to the simpler example

$K((x,y),(x^{\prime},y^{\prime}))=\left[\frac{1}{[m+n+\delta(y,y^{\prime})]^{q/2}}\phi\left(\frac{(x-x^{\prime})^{2}}{m+n+\delta(y,y^{\prime})}\right)\right]_{m,n=1}^{p}$

in $PD_{p}(\mathbb{R}\times S^{d})$.

A version of Theorem 3.1 for kernels acting on the product $X\times Y$ is as
follows.

######  Theorem 4.3.

Let $\rho$ be a nonzero positive measure on $(0,\infty)$ and $\phi$ a bounded
and completely monotone function. For each $m,n$ in $\\{1,\ldots,p\\}$, let
$G_{m,n}:Y\times Y\to M_{q}(\mathbb{R})$ be a matrix function with range
containing positive definite matrices only, $H_{m,n}:X\times
X\to\mathbb{R}^{q}$ vector functions and $\\{P_{m,n}^{s}\\}_{s>0}$ a family of
kernels on $X\times Y$ such that each function $s\in(0,\infty)\mapsto
P_{m,n}^{s}((x,y),(x^{\prime},y^{\prime}))$ is $\rho$-integrable. If the
matrix functions

$(y,y^{\prime})\in Y\times
Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$ $(x,x^{\prime})\in X\times
X\mapsto\left[e^{\displaystyle{i\,H_{m,n}(x,x^{\prime})^{\intercal}u}}\right]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

and

$((x,y),(x^{\prime},y^{\prime}))\in Y\times
Y\mapsto[P_{m,n}^{s}((x,y),(x^{\prime},y^{\prime}))]_{m,n=1}^{p},\quad s>0,$

belong to $CND_{p}(Y)$, $PD_{p}(Y)$, and $PD_{p}(X\times Y)$, respectively,
then the matrix kernel $K=[K_{m,n}]_{m,n=1}^{p}:(X\times Y)^{2}\to
M_{p}(\mathbb{R})$ given by the formula

$\displaystyle K_{m,n}((x,y)$ ,
$\displaystyle(x^{\prime},y^{\prime}))=\frac{1}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}$
$\displaystyle\times\int_{(0,\infty)}\phi\left(H_{m,n}(x,x^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(x,x^{\prime})\,s\right)P_{m,n}^{s}((x,y),(x^{\prime},y^{\prime}))d\rho(s)$

belongs to $PD_{p}(X\times Y)$.

We now move to some specific applications of Theorem 4.3.

######  Example 4.4.

Here we will employ the formula deduced in Theorem 1.1 in [6]:

$\mathcal{M}_{\nu}(r\sqrt{u})=\frac{r^{2\nu}}{2^{2\nu}\Gamma(\nu)}\int_{0}^{\infty}e^{\displaystyle{-s\;u}}e^{\displaystyle{-r^{2}/4s}}s^{-\nu-1}ds,\quad
r,u>0.$

that defines the so-called Matérn function. This function is studied in
details in [6]. We may apply Theorem 3.1 with $\phi(u)=\exp(-u)$, $u>0$ and
$d\rho(s)=e^{-r^{2}/4s}s^{-1}ds$. If for $x,x^{\prime}\in X$ and
$y,y^{\prime}\in Y$ we set

$2v_{m,n}((x,y),(x^{\prime},y^{\prime})):=v_{m}(x,y)+v_{n}(x^{\prime},y^{\prime}),$

where $v_{m}:X\times Y\to(0,\infty)$, for all $m$, and

$P^{s}_{m,n}((x,y),(x^{\prime},y^{\prime})):=\frac{r^{\displaystyle{2v_{m,n}((x,y),(x^{\prime},y^{\prime}))}}s^{\displaystyle{-v_{m,n}((x,y),(x^{\prime},y^{\prime}))}}}{2^{\displaystyle{2v_{m,n}((x,y),(x^{\prime},y^{\prime}))}}},$

for $x,x^{\prime}\in X$ and $y,y^{\prime}\in Y$, it is easily seen that the
kernels

$((x,y),(x^{\prime},y^{\prime}))\in(X\times
Y)^{2}\mapsto\left[P_{m,n}^{s}((x,x^{\prime}),(y,y^{\prime}))\right]_{m,n=1}^{p},\quad
s>0,$

belong to $PD_{p}(X\times Y)$. If each $s\in(0,\infty)\mapsto
s^{-v_{m,n}((x,y),(x^{\prime},y^{\prime}))/2}$ is $\rho$-integrable, Theorem
4.3 implies that the formula

$\displaystyle K_{m,n}((x,y)$ ,
$\displaystyle(x^{\prime},y^{\prime}))=\frac{\Gamma(v_{m,n}((x,y),(x^{\prime},y^{\prime}))}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}$
$\displaystyle\times\mathcal{M}_{v_{m,n}((x,y),(x^{\prime},y^{\prime}))}(r(H_{m,n}(x,x^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(x,x^{\prime}))^{1/2})$

defines a kernel
$K((x,y),(x^{\prime},y^{\prime}))=[K_{m,n}((x,y),(x^{\prime},y^{\prime}))]_{m,n=1}^{p}$
that belongs to $PD_{p}(X\times Y)$, as long as the $G_{m,n}$ and the
$H_{m,n}$ satisfy the assumptions of the theorem. We could also modify the
$P_{m,n}^{s}$ by introducing a matrix $[r_{m,n}]_{m,n=1}^{p}$ with positive
entries, by setting

$P^{s}_{m,n}((x,y),(x^{\prime},y^{\prime})):=\frac{r_{m,n}^{\displaystyle{2v_{m,n}((x,y),(x^{\prime},y^{\prime}))}}s^{\displaystyle{-v_{m,n}((x,y),(x^{\prime},y^{\prime}))}}}{2^{\displaystyle{2v_{m,n}((x,y),(x^{\prime},y^{\prime}))}}},$

for $x,x^{\prime}\in X$ and $y,y^{\prime}\in Y$, as long as the kernels

$((x,y),(x^{\prime},y^{\prime})\in(X\times
Y)^{2}\mapsto\left[P_{m,n}^{s}((x,y),(x^{\prime},y^{\prime}))\right]_{m,n=1}^{p},\quad
s>0,$

stay in $PD_{p}(X\times Y)$. In this case, the outcome of Theorem 3.1 would be
that the formula

$\displaystyle K_{m,n}((x,y)$ ,
$\displaystyle(x^{\prime},y^{\prime}))=\frac{\Gamma(v_{m,n}((x,y),(x^{\prime},y^{\prime}))}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}$
$\displaystyle\times\mathcal{M}_{v_{m,n}((x,y),(x^{\prime},y^{\prime}))}(r_{m,n}(H_{m,n}(x,x^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(x,x^{\prime}))^{1/2})$

defines a kernel
$K((x,y),(x^{\prime},y^{\prime}))=[K_{m,n}((x,y),(x^{\prime},y^{\prime}))]_{m,n=1}^{p}$
in $PD_{p}(X\times Y)$, if we keep the assumptions on the $G_{m,n}$ and the
$H_{m,n}$ required in the theorem. An specific and simple example in the
space-time setting can be produced in analogy with Theorem 1 in [5]: set
$Y=\mathbb{R}^{d}$, $X=\mathbb{R}$,

$G_{m,n}(y,y^{\prime})=g(\|y-y^{\prime}\|^{2})I_{q},\quad
y,y^{\prime}\in\mathbb{R}^{d};m,n=1,\ldots,p,$

where $g:(0,\infty)\to(0,\infty)$ has a completely monotone derivative and

$H_{m,n}(x,x^{\prime})=x-x^{\prime},\quad
x,x^{\prime}\in\mathbb{R};m,n=1,\ldots,p.$

Since $(y,y^{\prime})\in\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto
g(\|y-y^{\prime}\|^{2})$ belongs to $CND_{1}(\mathbb{R}^{d})$ by a result of
Micchelli ([15]), it follows that the matrix kernels $(y,y^{\prime})\in
Y\times Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p}$,
$u\in\mathbb{R}^{q}$, belong to $CND_{p}(\mathbb{R}^{d})$. If we put

$v_{m,n}((x,y),(x^{\prime},y^{\prime}))=\frac{v_{m}+v_{n}}{2},\quad
x,x^{\prime}\in\mathbb{R};y,y^{\prime}\in Y;m,n=1,\ldots,p,$

in which each $v_{m}$ is a positive constant and properly specify
$[r_{m,n}]_{m,n=1}^{p}$, then for $x,x^{\prime}\in\mathbb{R}$ and
$y,y^{\prime}\in\mathbb{R}^{d}$ the formula

$P^{s}_{m,n}((x,y),(x^{\prime},y^{\prime})):=\frac{r_{m,n}^{\displaystyle{v_{m}+v_{n}}}s^{\displaystyle{-(v_{m}+v_{n})/2}}}{2^{\displaystyle{v_{m}+v_{n}}}},\
\ m,n=1,\ldots,p,$

defines kernels

$((x,y),(x^{\prime},y^{\prime}))\in(\mathbb{R}\times\mathbb{R}^{d})^{2}\mapsto[P^{s}_{m,n}((x,y),(x^{\prime},y^{\prime}))]_{m,n=1}^{p},\quad
s>0$

in $PD_{p}(\mathbb{R}\times\mathbb{R}^{d})$. An application of Theorem 4.3
would lead to

$K_{m,n}((x,y),(x^{\prime},y^{\prime}))=\frac{\Gamma((v_{m}+v_{n})/2)}{g(\|y-y^{\prime}\|^{p/2}}\mathcal{M}_{(v_{m}+v_{n})/2}\left(r_{mn}\frac{\|x-x^{\prime}\|^{2}}{g(\|y-y^{\prime}\|^{2})}\right)$

with
$K((x,y),(x^{\prime},y^{\prime}))=[K_{m,n}((x,y),(x^{\prime},y^{\prime}))]_{m,n=1}^{p}$
in $PD_{p}(\mathbb{R}\times\mathbb{R}^{d})$. We observe that the factor
$\Gamma((v_{m}+v_{n})/2)$ can be eliminated as long as we can specify
$[r_{mn}]_{m,n=1}^{p}$ in such a way that
$[r_{m,n}^{v_{m}+v_{n}}/\Gamma((v_{m}+v_{n})/2)]_{m,n=1}^{p}$ is a positive
definite matrix. Theorem 1 in [11] is another construction that fits into
Theorem 4.3. Details on that will be left to the readers.

######  Example 4.5.

The so-called generalized Cauchy function ([8, p.337]) is given by

$\frac{1}{(1+cu^{\gamma})^{\nu}}=\frac{c^{-\nu}}{\Gamma(\nu)}\int_{0}^{\infty}e^{-s\,u^{\gamma}}s^{\nu}d\rho(s),\quad
u\geq 0,$

where $c>0$, $\nu>1$, $\gamma\in(0,1]$, and $d\rho(s)=s^{-1}\exp(-\,s/c)$. In
order to apply Theorem 4.3 we now set $\phi(u)=e^{-u^{\gamma}}$, $u>0$ and

$P_{m,n}^{s}((x,y),(x^{\prime},y^{\prime})=\left(\frac{s}{c}\right)^{v_{m}(x,y)+v_{n}(x^{\prime},y^{\prime})},\quad
s>0;\,y,y^{\prime}\in Y,$

where $v_{m}:X\times Y\to(0,\infty)$ is chosen in such a way that each
$s\in(0,\infty)\mapsto s^{v_{m}(x,y)/2}$ is $\rho$-integrable. The outcome is
that

$\displaystyle K_{m,n}((x,y),(x^{\prime},y^{\prime}))$ $\displaystyle=$
$\displaystyle\frac{\Gamma(v_{m}(x,y)+v_{n}(x^{\prime},y^{\prime}))}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}$
$\displaystyle\times\frac{1}{\left(1+c(H_{m,n}(x,x^{\prime})^{\intercal}G_{m,n}(y,y^{\prime})^{-1}H_{m,n}(x,x^{\prime})^{\gamma}\right)^{v_{m}(x,y)+v_{n}(x^{\prime},y^{\prime})}},$

defines a kernel
$K((x,y),(x^{\prime},y^{\prime}))=[K_{m,n}((x,y),(x^{\prime},y^{\prime}))]_{m,n=1}^{p}$
in $PD_{p}(X\times Y)$, if we keep the assumptions on the $G_{m,n}$ and the
$H_{m,n}$ required in the theorem. Arguments similar to those developed in the
second half of Example 4.4 leads to an example aligned with Theorem 2 in [5].

##  5 A further extension

As a final remark let us point an improvement that one can make in all the
theorems proved in this paper. If for each $m$ and $n$ in $\\{1,\ldots,p\\}$,
$G_{m,n}:Y\times Y\to M_{q}(\mathbb{R})$ is a matrix function with range
containing positive definite matrices only, Theorem 2.4 justifies the
following fact: if the matrix kernels

$(y,y^{\prime})\in Y\times
Y\mapsto[u^{\intercal}G_{m,n}(y,y^{\prime})u]_{m,n=1}^{p},\quad
u\in\mathbb{R}^{q},$

belong to $CND_{p}(Y)$, then the kernel $K$ given by

$K(y,y^{\prime})=\left[\frac{1}{\sqrt{\det
G_{m,n}(y,y^{\prime})}}\right]_{m,n=1}^{p},\quad y,y^{\prime}\in Y,$

belongs to $PD_{p}(Y)$. Under the same setting, it follows from the Schur
Product Theorem that

$K_{l}(y,y^{\prime})=\left[\frac{1}{[\det
G_{m,n}(y,y^{\prime})]^{l/2}}\right]_{m,n=1}^{p},\quad y,y^{\prime}\in Y.$

belongs to $PD_{p}(Y)$ whenever $l\in\\{1,2,\ldots\\}$. In particular, we can
introduce the same power $l/2$ in the assertions of all the theorems proved in
the paper.

## References

  * [1] Alexander, R.; Stolarsky, K. B., Extremal problems of distance geometry related to energy integrals. Trans. Amer. Math. Soc. 193 (1974), 1-31.
  * [2] Alfonsi, A; Klock, F.; Schied, A., Multivariate transient price impact and matrix-valued positive definite functions. Math. Oper. Res. 41 (2016), no. 3, 914-934.
  * [3] Belton, A.; Guillot, D.; Khare, A.; Putinar, M., A Panorama of Positivity. I: Dimension Free. In: Aleman A., Hedenmalm H., Khavinson D., Putinar M. (eds). Analysis of Operators on Function Spaces. Trends in Mathematics. Birkháuser, Cham, 117-165, 2019.
  * [4] Berg, C.; Christensen, J. P. R.; Ressel, P., Harmonic analysis on semigroups. Theory of positive definite and related functions. Graduate Texts in Mathematics, 100. Springer-Verlag, New York, 1984.
  * [5] Bourotte, M.; Allard, D.; Porcu, E., A flexible class of non-separable cross-covariance functions for multivariate space-time data. Spat. Stat. 18 (2016), Part A, 125-146.
  * [6] Cho, Yong-Kum; Kim, Dohie; Park, Kyungwon; Yun, Hera, Schoenberg representations and Gramian matrices of Matérn functions, arXiv:1702.05894v1 [math.CA] (2017).
  * [7] Gneiting, T., Nonseparable, stationary covariance functions for space-time data. J. Amer. Stat. Assoc. 97 (2002), no. 458, 590-600.
  * [8] Gradshteyn, I. S.; Ryzhik, I. M., Table of integrals, series, and products. Translated from the Russian. Seventh edition. Elsevier/Academic Press, Amsterdam, 2007.
  * [9] Horn, R. A.; Johnson, C. R., Matrix analysis. Second edition. Cambridge University Press, Cambridge, 2013.
  * [10] Kapil, Y., Pal, R., Aggarwal, A., Singh, M., Conditionally negative definite functions. Mediterr. J. Math. 15 (2018), no. 5, Paper No. 199, 12 pp.
  * [11] Kleiber W.; Nychka, D., Nonstationary modeling for multivariate spatial processes. J. Multivariate Anal. 112 (2012), 76-91.
  * [12] Menegatto, V. A., Positive definite functions on products of metric spaces via generalized Stieltjes functions. Proc. Amer. Math. Soc. 148 (2020), no. 11, 4781-4795.
  * [13] Menegatto, V. A., Strictly positive definite kernels on the Hilbert sphere. Appl. Anal. 55 (1994), no. 1-2, 91-101.
  * [14] Menegatto, V. A.; Oliveira, C. P., An extension of Aitken’s integral for Gaussians and positive definiteness. Methods Appl. Anal. (2021), to appear.
  * [15] Micchelli, C. A., Interpolation of scattered data: distance matrices and conditionally positive definite functions, Constr. Approx. 2 (1986), 11-22.
  * [16] Micheli, M.; Glaunés, M. A., Matrix-valued kernels for shape deformation analysis. Geometry, Imaging, and Computing 1 (2014), no. 1, 57-139.
  * [17] Minh, H. Q.; Bazzani, L.; Murino, V., A unifying framework in vector-valued reproducing kernel Hilbert spaces for manifold regularization and co-regularized multi-view learning, J. Mach. Learn. Res. 17 (2016), 1-72.
  * [18] Porcu, E.; Furrer, R.; Nychka, D., 30 years of space-time covariance functions. Wiley Interdisciplinary Reviews: Computational Statistics 13 (2020), no. 3, e1512.
  * [19] Reams, R., Hadamard inverses, square roots and products of almost semidefinite matrices. Linear Algebra Appl. 288 (1999), no. 1-3, 35-43 .
  * [20] Schlather, M., Some covariance models based on normal scale mixtures. Bernoulli 16 (2010), no. 3, 780-797.
  * [21] Schilling, R. L.; Song, R.; Vondracek, Z., Bernstein functions. Theory and applications. Second edition. De Gruyter Studies in Mathematics, 37. Walter de Gruyter $\&$ Co., Berlin, 2012.
  * [22] Searle, S. R., Matrix algebra useful for statistics. Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. John Wiley $\&$ Sons, Ltd., Chichester, 1982.
  * [23] D. Wittwar, G. Santin, and B. Haasdonk, Interpolation with uncoupled separable matrix-valued kernels. Dolomites research notes on approximation. Special issue of the “Seminari Padovani di Analisi Numerica 2018” (SPAN2018), Volume 11, 23-39, 2018.

V. A. Menegatto
Departamento de Matemática - ICMC-USP - São Carlos
Caixa Postal 668
13560-970, São Carlos SP, Brazil
E-mail<EMAIL_ADDRESS>

C. P. Oliveira
Instituto de Matemática e Computação - UNIFEI
Av. BPS, 1303, Pinheirinho
37500-903, Itajubá MG, Brazil.
E-mail<EMAIL_ADDRESS>