What happens if you square an RBF kernel function? Let's say we use a kernel regularization algorithm such as ridge regression to minimize some loss in an RBF kernel:
$$\min_{h \in H} \frac{1}{n} \sum_i (h(x_i) - y(x_i))^2 + ||h||^2_K$$
We get some $h$ in the RBF kernel space. What now happens if we square h? For example, if we look at $h^2$, is it still in an RBF kernel space or not?
Due to the representer theorem we know that $h$ is a sum of kernel products.
A guess is to say that $h(x) = e^{-||x - z||^2/\sigma^2}$, and if we look at $h^2$ we get that $h^2(x) = e^{-2||x - z||^2/\sigma^2}$. 
But if $h$ contains two terms: $h(x) = e^{-||x - z||^2/\sigma^2} + e^{-||x - y||^2/\sigma^2}$, and we look at $h^2(x)$ we find that:
$$h^2(x) = e^{-2||x - z||^2/\sigma^2} + e^{-2||x - y||^2/\sigma^2} - 2e^{(-||x-y||^2-||x-z||^2)/\sigma^2}$$
And I can't see now we have two terms, how to rewrite it to a sum of RBF kernel products. So does it stay in the same RBF space? And how about the norm of h?
 A: To start with, you're slightly off on the representer theorem; it means that $h$ is a linear combination of kernel functions. That is, if your input data points are $\{z_i\}_{i=1}^n$, then
$$
h(x) = \sum_{i=1}^n \alpha_i \, k(x, z_i)
$$
for some real coefficients $\alpha_i$.
Thus
$$
h^2(x) = \left( \sum_i \alpha_i \, k(x, z_i) \right)^2
= \sum_i \sum_j \alpha_i \,\alpha_j\, k(x, z_i) \, k(x, z_j)
,$$
which will be in a given RKHS as long as $k(x, z_i) \, k(x, z_j)$ is in the RKHS.
Now specialize to the class of RBF kernels
$$
k_\sigma(x, z) = \exp\left( - \frac{1}{2 \sigma^2} \lVert x - z \rVert^2 \right)
,$$
and define the corresponding RKHS as $\mathcal H_\sigma$.
(Note that your definition didn't have the 2 in the denominator, which is nonstandard and somewhat misleading; I'll use this one.)
Now, it turns out that:
$$
\begin{align}
k_\sigma(x, z_i) \, k_\sigma(x, z_j)
&= k_{\sqrt{2}\sigma}(z_i, z_j) \, k_{\sigma/\sqrt{2}}\left(x, \frac{z_i + z_j}{2} \right)
\end{align}
.$$
Thus
\begin{align}
h^2(x)
&= \sum_{ij} \underbrace{\alpha_i \alpha_j k_{\sqrt{2}\sigma}(z_i, z_j)}_{\beta_{ij}} \, k_{\sigma/\sqrt{2}}\Big(x, \underbrace{\frac{z_i + z_j}{2}}_{\zeta_{ij}} \Big)
\\&= \sum_{ij} \beta_{ij} \, k_{\sigma/\sqrt{2}}(x, \zeta_{ij})
,\end{align}
and $h^2 \in \mathcal H_{\sigma/\sqrt{2}}$.
In the interest of brevity, define $\tau := \sigma / \sqrt{2}$.
Now, turn to the following paper:

Tanaka, Imai, Kudo, and Miyakoshi, Theoretical analyses on a class of nested RKHS's, ICASSP 2011 (doi, author's pdf) 

They show a few interesting facts in their Section 5:
First, $H_\sigma \subset H_\tau$ for all $\tau < \sigma$ (Theorem 5), and moreover $\lVert f \rVert_{\mathcal H_\sigma} \ge \lVert f \rVert_{\mathcal H_\tau}$ for every $f \in \mathcal H_{\sigma}$ (consequence of Theorem 4).
This doesn't tell us anything in particular about whether $h^2 \in H_\sigma$, but it does let us use their Proposition 1 as follows:
Let $\nu := \sqrt{ \sigma^2 + \sigma^2 - (\sigma / \sqrt{2})^2} = \sqrt{3/2} \sigma$. Then
$$\langle k_{\sigma}(\cdot, x), k_{\sigma}(\cdot, y) \rangle_{\mathcal H_{\tau}} = k_\nu(x, y).$$
Thus
\begin{align}
\lVert h \rVert_{\mathcal H_\tau}^2
&= \lVert \sum_i \alpha_i k_\sigma(\cdot, z_i) \rVert_{\mathcal H_\tau}^2
\\&= \sum_{ij} \alpha_i \alpha_j \langle k_\sigma(\cdot, z_i), k_\sigma(\cdot, z_j)\rangle_{\mathcal H_\tau}
\\&= \sum_{ij} \alpha_i \alpha_j k_\nu(z_i, z_j)
\end{align}
whereas
\begin{align}
\lVert h^2 \rVert_{\mathcal H_\tau}^2
&= \lVert \sum_{ij} \beta_{ij} k_\tau(\cdot, \zeta_{ij}) \rVert_{\mathcal H_\tau}^2
\\&= \sum_{abij} \beta_{ab} \beta_{ij} \langle k_\tau(\cdot, \zeta_{ab}), k_\tau(\cdot, \zeta_{ij})\rangle_{\mathcal H_\tau}
\\&= \sum_{abij} \beta_{ab} \beta_{ij} k_\tau(\zeta_{ab}, \zeta_{ij})
\\&= \sum_{abij} \alpha_a \alpha_b \alpha_i \alpha_j k_{\sqrt{2}\sigma}(z_a, z_b) k_{\sqrt{2}\sigma}(z_i, z_j) k_\tau(\zeta_{ab}, \zeta_{ij})
.\end{align}
The diagonal terms ($a = i$, $b = j$) of the latter sum look kind of like the first sum, but the off-diagonals can be either positive or negative, so that doesn't lead to an immediate bound. At least in this form, though, the norms are numerically comparable.
