  The Gaussian kernel plays a central role in machine learning, uncertainty
quantification and scattered data approximation, but has received relatively
little attention from a numerical analysis standpoint. The basic problem of
finding an algorithm for efficient numerical integration of functions
reproduced by Gaussian kernels has not been fully solved. In this article we
construct two classes of algorithms that use $N$ evaluations to integrate
$d$-variate functions reproduced by Gaussian kernels and prove the exponential
or super-algebraic decay of their worst-case errors. In contrast to earlier
work, no constraints are placed on the length-scale parameter of the Gaussian
kernel. The first class of algorithms is obtained via an appropriate scaling of
the classical Gauss-Hermite rules. For these algorithms we derive lower and
upper bounds on the worst-case error of the forms $\exp(-c_1 N^{1/d})
N^{1/(4d)}$ and $\exp(-c_2 N^{1/d}) N^{-1/(4d)}$, respectively, for positive
constants $c_1 > c_2$. The second class of algorithms we construct is more
flexible and uses worst-case optimal weights for points that may be taken as a
nested sequence. For these algorithms we derive upper bounds of the form
$\exp(-c_3 N^{1/(2d)})$ for a positive constant $c_3$.
