I will use method 1. Check Douglas Zare's answer for a proof using method 2.
I will prove the case when $x,y$ are real numbers, so $k(x,y)=\exp(-(x-y)^2/2\sigma^2)$. The general case follows mutatis mutandis from the same argument, and is worth doing.
Without loss of generality, suppose that $\sigma^2=1$.
Write $k(x,y)=h(x-y)$, where $$h(t)=\exp\left(-\frac{t^2}{2}\right)=\mathrm{E}\left[e^{itZ}\right] $$ is the characteristic function of a random variable $Z$ with $N(0,1)$ distribution.
For real numbers $x_1,\dots,x_n$ and $a_1,\dots,a_n$, we have
$$
\sum_{j,k=1}^n a_j\,a_k\,h(x_j-x_k) = \sum_{j,k=1}^n a_j\,a_k\,\mathrm{E} \left[ e^{i(x_j-x_k)Z}\right] = \mathrm{E} \left[ \sum_{j,k=1}^n a_j\,e^{i x_j Z}\,a_k\,e^{-i x_k Z}\right]
= \mathrm{E}\left[ \left| \sum_{j=1}^n a_j\,e^{i x_j Z}\right|^2\right] \geq 0 \, ,
$$
which entails that $k$ is a positive semidefinite function, aka a kernel.
To understand this result in greater generality, check out Bochner's Theorem: http://en.wikipedia.org/wiki/Positive-definite_function