When using Gaussian processes, the covariance matrix $\mathbf{\Sigma}$ is often defined via a covariance function $K$ as follows $$ \mathbf{\Sigma}_{ij} = K(\underline{x}_i, \underline{x}_j) $$ where $\underline{x}_i, \underline{x}_j$ are coordinates of two points in some space of interest, and belong to some finite set of $n$ such points, resulting in a $n \times n$ covariance matrix.
Is it known whether the Gaussian (or 'squared exponential') covariance function $$ K(\underline{x}_i, \underline{x}_j) = a^2 \exp{\left(-\frac{(\underline{x}_i - \underline{x}_j)^{\top}(\underline{x}_i - \underline{x}_j)}{2l^2}\right)} $$ will always produce a positive-definite $\mathbf{\Sigma}$ for any chosen set of points and values of $a, l$?
Thanks!