The abstract theory of RKHSs has been developed over a number of years outside the domain of physics (e.g., for the study of conformal mappings [11], integral equations [12], and partial differential equations [13]). Loéve has demonstrated that the RKHS theory provides a unified framework for stochastic processes and signal processing [14]. Duchon [15] and Meinguet [16] have formulated generalized smooth surface spline functions using the RKHS technique in Sobolev space [17]. Based on these rigorous foundations, the RKHS method has been used for a variety of applications, especially in data interpolation and smoothing [18, 19, 20, 21, 22, 23, 24]. The RKHS method provides a rigorous and effective framework for smooth multivariate interpolation of arbitrarily scattered data and for accurate approximation of general multidimensional functions. These features can be attributed to the reproducing property
of the reproducing kernel (rk) 
associated with an RKHS 
[4-7], for
every function 
belonging to .
The subscript 
by the inner product 
indicates that the inner product operation applies to
functions of . The relation in Eq. (1) may be
viewed as a continuum limit of interpolation where the rk
acts to weight the values of 
such that is reliably (exactly) generated.
The reproducing property, Eq. (1), imparts a
rich physically based structure in the associated Hilbert space that possesses
many important properties (e.g., the uniqueness and positive definiteness
of the reproducing kernel which are important for its practical utility).
It is readily seen that the rk
plays a role similar to that of the Dirac delta function
in a typical 
Hilbert space [25].
However, in contrast to the delta function, the reproducing kernels are
continuous bounded functions and can be tailored to carry physical
information on the particular system or problem [26].
As a result, the projection of an arbitrary function into an RKHS renders
the development of approximation or interpolation
schemes tailored to the inherent properties of the function.
Reproducing kernels can in general be constructed from any complete
system of linearly independent or compact
functions (e.g., orthogonal polynomials [5]
and wavelets [27]). In particular, for
accurate representations of the smoothly varying PES
, reproducing kernels of appropriate smoothness
can be constructed in the context of
Sobolev spaces (i.e., RKHSs endowed with inner products
involving derivatives of continuously differentiable functions [17]).
For example, we have successfully constructed
one-dimensional reproducing kernels (Taylor splines):
for distance-like variables and
for angle-like variables [8].
Here, the integers 
and 
characterize the asymptotic behavior of the kernels, 
and 
are,
respectively, the larger and smaller of x and 
, B(a,b) is the
beta function, and 
is the Gauss hypergeometric
function. These two rks have been successfully adopted,
in tensor-product forms, for
constructing accurate and fast interpolated multi-dimensional
polyatomic potential energy surfaces
for various molecular systems from high-level scattered or
gridded ab initio
points [2, 3, 8, 9, 10].
Smooth global multi-dimensional reproducing kernels
have been successfully utilized in other contexts for multivariate
interpolation (e.g., in computer aided geometric
design [28, 29])
and to solve differential equations by collocation [30].
These rks usually are of simple and easily computable closed
forms [31, 32, 33, 34], e.g.
(1) Sobolev spline: 
a modified
Bessel function of the third kind and (2) generalized inverse multi-quadric:
an odd positive
integer and c a real number. Both the Sobolev spline [35] and
generalized inverse multi-quadric [36] have only recently been
adopted for successfully solving
bound-state Schrödinger equations.
Numerical experiments have shown that Sobolev spline
and generalized inverse multi-quadric are numerically stable, due to
their mild decaying nature, with respect to
the number and distribution of grid
points [34, 35, 36].
In order for the RKHS potential energy surfaces (PESs) to be useful in dynamic and spectroscopic studies, these PESs (and their derivatives) must possess at least the following four attributes: (i) conformity - correct symmetry properties and asymptotic behavior; (ii) accuracy - absolute accuracy in the important regions; (iii) smoothness - continuity and differentiability as a function of molecular coordinates; and (iv) efficiency - fast to evaluate at any molecular configuration [39]. The last property plays a crucial role in any molecular simulation study. Specifically, the goal here is to fully explore and apply the powerful RKHS method for accurately constructing smooth, and rapidly interpolated potential energy surfaces using high-level ab initio data.