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

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.