Following @Abdelmalek's great advice in the comments above it is enough to compute the following triple integral: \begin{align*} I=\frac{\pi^s}{\Gamma(s)}\int_{0}^{\infty}\left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} (x_1+x_2 i)^{\ell}\cdot e^{-\pi t(x_1^2+x_2^2+a^2)} e^{2\pi i (x_1\zeta_1+x_2\zeta_2)} dx_1 dx_2 \right) t^s\frac{dt}{t}. \end{align*} So we have \begin{align*} I=&=\frac{\pi^s}{\Gamma(s)}\int_{0}^{\infty}t^{-1}\cdot\left(\frac{i}{t}(\zeta_1+i\zeta_2)\right)^{\ell}\cdot e^{-\frac{\pi}{t}(\zeta_1^2+\zeta_2^2)}\cdot e^{-\pi ta^2}t^{s}\frac{dt}{t}\\ &=\frac{(i)^{\ell}\pi^s}{\Gamma(s)}(\zeta_1+i\zeta_2)^{\ell}\int_{0}^{\infty} e^{-\frac{\pi}{t}(\zeta_1^2+\zeta_2^2)-\pi ta^2}\cdot t^{s-\ell-1} \frac{dt}{t} \\ &=\frac{(i)^{\ell}\pi^s}{\Gamma(s)}(\zeta_1+i\zeta_2)^{\ell}\cdot 2\left(\frac{|\zeta|}{ta}\right)^{s-\ell-1} K_{s-\ell-1}\left(2\pi|\zeta|a\right) \end{align*} where $\zeta=\zeta_1+i\zeta_2$.
So modulo some minor mistakes made above, out of excitement, it means that the initial messy formulas that I had obtained could be vastly simplified, as what I was hoping for. In any case this new approach is just better and simpler!
I can't wait for somebody to publish a book with tables of double integrals (or even multiple integrals)!