Proving PSD-ness of matrices
To get the obvious out of the way first, to show that the $n\times n$ matrix $M$ is PSD, it is enough to show that it is the Gram matrix of $n$ vectors, i.e. $m_{ij} = \langle u_i, u_j\rangle$ where $u_1, \ldots, u_n \in \mathbb{R}^n$. (This is a "sum of squares" proof of positive semi-definiteness, because it shows that the bi-linear form $x^T M x$ can be written as a sum of squares of linear functions of $x$.  I am wondering if you meant something more sophisticated by "sum of squares".)

You get a useful class of PSD matrices from positive definite functions. A function $f: \mathbb{R}^d \to \mathbb{R}$ is positive definite if for all sequences $u_1, \ldots, u_n \in \mathbb{R}^d$, the $n\times n$ matrix $M$ defined by $m_{ij} = f(u_i - u_j)$ is PSD. It is not very hard to see that if $f$ is the Fourier transform of a finite measure, then $f$ is positive definite. I.e. it is sufficient that there is a finite Borel measure $\mu$  on $\mathbb{R}^d$ such that for all $x$
$$
f(x) = \int{e^{i\langle x, y\rangle} d\mu(y)}.
$$
A classical theorem of Bochner shows that this condition is necessary as well.

A special class of positive definite functions are the functions $f(x) = g(\|x\|_2)$, for which the matrix $M$ becomes $m_{ij} = g(\|u_i - u_j\|_2)$. Another classical result of Schoenberg shows that such an $f$ is positive definite if and only if there exists a finite Borel measure $\mu$ on $[0, \infty)$ so that for every $t \geq 0$
$$
g(t) = \int_{0}^\infty{e^{-t^2s}d\mu(s)}.
$$

The original motivation for Schoenberg's theorem was to determine what functions $g$ have the property that for any $u_1, \ldots, u_n \in \mathbb{R}^d$, the metric space with metric $d(u_i, u_j) = g(\|u_i - u_j\|_2)$ embeds isometrically into Euclidean space.

This blog post has a simple proof of Schoenberg's theorem and many links.