Suppose we have the Gaussian kernel, $k(x, y) = \exp\left(-\frac{||x-y||^2}{2\sigma^2}\right), x, y \in \mathbb{R}^d$. Also, assume that $k(.,.)$ has a Gaussian distribution, i.e., $k \sim \mathcal{N}(\mu, \lambda^2)$. I was trying to find the distributions from which sampling $x \mbox{ and } y$ would give such a kernel.
Let's try to solve the simpler case where $d = 1$. Assume that $X, Y \sim g$. After solving for a bit and simplifying a lot of equations, I get the following equality,
$$\int_{-\infty}^{\infty}g(y)\left[g(y + x) + g(y - x)\right]dy = \frac{x\exp\left({-\frac{x^2}{2\sigma^2} -\frac{(l- \mu)^2}{2\lambda^2}}\right)}{\sqrt{2\pi\lambda^2}\sigma^2}$$ where $l = \exp\left({-\frac{x^2}{2\sigma^2}}\right)$.
I spent a lot of time trying to solve this, but my familiarity with solving such equations is very limited and can't figure out if one can even get a closed form for $g$. Any help would be great!
P.S.: If we can solve this, is it possible to extend it to any general $d$?