This is to address integrability issues not addressed in @PaulEnta answer to this problem.
As observed, after a change of variables using polar coordinates , the desired integral $I$ may be expressed as the convolution at $(c,0)$ of the functions
\begin{align}
f_a(x,y)&=\frac{a}{2\pi}\frac{1}{(x^2+y^2+a^2)^{3/2}}\\
g_b(x,y)&= \frac{1}{\sqrt{x^2+y^2+b^2}}
\end{align}
that is, $I=(f_a*g_b)(c,0)=\int_{\mathbb{R}^2}f_a(x,y)g_b(c-x,-y)\,dxdy$
Connections to well known Fourier transforms:
- It is known from Fourier analysis (Stein and Weiss for instance) that
$f_a$ is related to Poisson kernel in 2D, and that it is the Fourier transform of a $L_1$ function:
Let $\Phi_a(\boldsymbol{x})=e^{-2\pi a|\boldsymbol{x}|}$, which is in $L_1(\mathbb{R}^2)$.
$$
\widehat{\Phi_1}(\boldsymbol{y})=\int_{\mathbb{R}^2} e^{-2\pi i \boldsymbol{x}\cdot\boldsymbol{t}} e^{-2\pi|\boldsymbol{x}|},d\boldsymbol{x} =\frac{1}{2\pi}\frac{1}{(1+|\boldsymbol{y}|^2)^{3/2}}
$$
and so,
$$
\widehat{\Phi_a}(\boldsymbol{y})=a^{-2}\widehat{\Phi_1}(a^{-1}\boldsymbol{y})=\frac{a}{2\pi(a^2+|\boldsymbol{y}|^2)^{3/2}}=f_a(\boldsymbol{y})$$
Observed that $f_a\in L_1(\mathbb{R}^2)$ as a simple calculation using polar coordinates shows.
- The function $g_b$ is the Fourier transform of an $L_1$ function:
Let $\Psi_b(\boldsymbol{x})=\frac{e^{-2\pi b|\boldsymbol{x}|}}{2\pi b|x|}$ for $\boldsymbol{x}\neq\boldsymbol{0}$, which is in $L_1(\mathbb{R}^2)$.
$$
\widehat{\Psi_1}(\boldsymbol{y})=\frac{1}{2\pi \sqrt{1+|\boldsymbol{y}|^2}}
$$
This particular Fourier transform is computed here. Hence
$$
\widehat{\Psi_b}(\boldsymbol{y})=b^{-2}\widehat{\Phi_1}(b^{-1}\boldsymbol{y}) =\frac{1}{2\pi b \sqrt{b^2+|\boldsymbol{y}|^2}}=\frac{1}{2\pi b}g_b(\boldsymbol{y})
$$
Notice however that $g_b\not\in L_1(\mathbb{R}^2)$.
- By Young's convolution theorem $f_a*g_b\in\mathcal{C}_\infty(\mathbb{R}^2)$. In order to invoke the Fourier inversion theorem, some integrability considerations are in order:
Lemma: Assume $\Phi,\Psi\in L_1(\mathbb{R}^n)$ and $\widehat{\Phi}\in L_1(\mathbb{R}^n)$. Then $\Phi\cdot\Psi\in L_1(\mathbb{R}^n)$, and
\begin{align}
\widehat{\Phi\cdot\Psi}=\widehat{\Phi}*\widehat{\Psi}\tag{1}\label{one}
\end{align}
Proof: The Fourier inversion theorem for $L_1$ implies that $\Phi$ is almost surely equal to function (which we still denote by $\Phi$) in $\mathcal{C}_\infty(\mathbb{R}^n)\cap L_1(\mathbb{R}^n)$ and
$$
\Phi(\boldsymbol{x})=\int_{\mathbb{R}^n}\widehat{\Phi}(\boldsymbol{y})e^{2\pi i \boldsymbol{x}\cdot\boldsymbol{y}}\,d\boldsymbol{y}
$$
It follows that $\Phi\cdot\Psi\in L_1(\mathbb{R}^n)$ and so,
\begin{align}
\widehat{\Phi\cdot\Psi}(\boldsymbol{t})=\int_{\mathbb{R}^n}\Phi(\boldsymbol{x})\Psi(\boldsymbol{x}) e^{-2\pi\boldsymbol{x}\cdot\boldsymbol{t}}\,d\boldsymbol{x}=\int_{\mathbb{R}^n}\Psi(\boldsymbol{x}) e^{-2\pi\boldsymbol{x}\cdot\boldsymbol{t}}\Big(\int_{\mathbb{R}^n}\widehat{\Phi}(\boldsymbol{y})e^{2\pi i \boldsymbol{x}\cdot\boldsymbol{y}}\,d\boldsymbol{y}\Big)
\,d\boldsymbol{x}
\end{align}
Since $(\boldsymbol{x},\boldsymbol{y})\mapsto \Psi(\boldsymbol{x})\widehat{\Phi}(\boldsymbol{y})$ is in $L_1(\mathbb{R}^n\times\mathbb{R}^n)$, we can apply Fubini's theorem to obtain
\begin{align}
\widehat{\Phi\cdot\Psi}(\boldsymbol{t})= \int_{\mathbb{R}^n}\widehat{\Phi}(\boldsymbol{y}) \Big(\int_{\mathbb{R}^n}\Psi(\boldsymbol{x})e^{-2\pi i \boldsymbol{x}\cdot(\boldsymbol{t}-\boldsymbol{y})}\,d\boldsymbol{x}\Big)
\,d\boldsymbol{y}=\int_{\mathbb{R}^n}\widehat{\Phi}(\boldsymbol{y}) \widehat{\Psi}(\boldsymbol{t}-\boldsymbol{y})\,d\boldsymbol{y}=\widehat{\Phi}*\widehat{\Psi}(\boldsymbol{t}).
\end{align}
Remark: The equality in \eqref{one} equality among continuous functions (pointwise).
- Now that we have taken all integrability considerations out of the way, we can proceed as in @PaulEnta's (our lemma outlines the strategy)
\begin{align}
f_a*g_b(t_1,t_2)&=2\pi b\widehat{\Phi_a\cdot\Psi_b}(t_1,t_2)=2\pi b\int_{\mathbb{R}^2}\frac{e^{-2\pi|\boldsymbol{x}|(a+b)}}{2\pi b|\boldsymbol{x}|}e^{-2\pi i\boldsymbol{x}\cdot(t_1,t_2)}\,d\boldsymbol{x}\\
&=2\pi(a+b)\widehat{\Psi_{a+b}}(t_1,t_2)=\frac{1}{\sqrt{(a+b)^2+t^2_1 + t^2_2 }}
\end{align}