The following integral appears naturally within the computation of the Fourier series coefficients of a real analytic $\mathrm{SL}_2(\mathbb{C})$-Eisenstein series: \begin{align*} \int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty} \frac{(x_1+ i x_2)^{\ell}}{(x_1^2+x_2^2+a^2)^s} e^{ i x_1\zeta_1}dx_1 \right) e^{ i x_2\zeta_2} dx_2, \end{align*} where $i=\sqrt{-1}$, $\zeta_1,\zeta_2,a\in\mathbb{R}$, $\ell\in\mathbb{Z}_{\geq 0}$ and $s\in\mathbb{C}$ with $\Re(2s-\ell-1)>1$. Using tables of one dimensional integrals, it is possible to write down this double integral as a "linear combination" of $K$-Bessel functions $K_{\nu}\left(a\sqrt{\zeta_1^2+\zeta_2^2}\right)$ for various order $\nu$ which is very similar to some of the formulas which appear in Friedberg's paper [1], see below.
Here is ma question:
Q: Is it possible to express this double integral in a very concise way using an appropriate hypergeometric function (e.g. a Meijer G-function with the appropriate parameters or something alike) ?
[1] S. Friedberg, On Maass wave forms and the imaginary quadratic Doi-Naganuma lifting. Math. Ann. 263 (1983), no. 4, 483–508