I want to decompose the Laplacian kernel $K(x,y)=e^{-|x-y|_2}$, with $x,y \in \mathbb{R^d}$, over the Gegenbauer polynomials of degree k: $Q_k^{(d)}(\vec{x}\cdot \vec{y})$, which are linked with the spherical harmonics $Y_{km}$ as follows: $$Q_k^{(d)}(\vec{x}\cdot \vec{y})=\frac{\alpha}{\alpha+k}\sum_{m=1}^{N(d,k)}Y_{km}(\vec{x})Y_{km}(\vec{y}),$$ with $\alpha=(d-2)/2$ and $N(d,k)$ the number of degree $k$ harmonics. The notation $|.|_2$ stands for the Euclidean norm.
In other words, I want to compute the following $g_k$: $$ e^{-|x-y|_2} =\sum_k Q_k^{(d)}(\vec{x}\cdot \vec{y}) g_k(|\vec{x}|_2,|\vec{y}|_2). $$ Does anyone have some suggestions? Even an ansatz for the form of $g_k$.