Double integral of two Gaussians and few complex poles

Question

Recently encountered an integral: $$ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \dfrac{ e^{-i(x_1+x_2)k} \exp\left(-\frac{(x_1-x_0)^2}{2\sigma^2} -\frac{(x_2-x_0)^2}{2\sigma^2}\right) }{(x_1+x_2-\alpha-\beta)(x_1-\gamma)(x_2-\delta)} d x_1 d x_2,$$ where $k, \sigma$ are Real, $\alpha, \beta, \gamma, \delta$ are complex, but their imaginary parts are negative.

I have tried to change variables to $x_+=x_1+x_2, x_-=x_1-x_2$, it results: $$ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \dfrac{ e^{-ix_+k} \exp\left( - \frac{x_-^2}{4 \sigma^2} - \frac{(x_+-2x_0)^2}{4 \sigma^2} \right) }{(x_+-\alpha-\beta)(x_++x_--\gamma)(x_+-x_--\delta)} 2 d x_+ d x_-$$.

But did not help much. It seems impossible to factorize the integrand into $f_1(x_1) f_2(x_2)$ because of $x_1+x_2$.

What I know is that the integral of a Gaussian with a complex pole can be reduced to complex valued error function. For exmaple, $\int_{-\infty}^{\infty} \frac{e^{-x^2/\sigma^2}}{x-\alpha}dx = - i \pi w \left( -\frac{\alpha}{\sigma} \right)$, where $w(z)$ is Faddeeva/Kramp function, here $\text{Im}(\alpha)<0$, and $\sigma$ is real, as before.

Even if I somehow can get rid of one of the integral, the second one will contain Gaussian, error function, and pole(s). What are the chances to have an answer for this integral in terms of special functions?

UPDATE: After expanding into the minimal denominators with $x_1, x_2$, I realized that the only part I can not compute has the following form:

$$ \int_{-\infty}^{\infty} \dfrac{w\left( \eta + \xi x_2 \right) e^{-\frac{(x_2 - \nu )^2}{2 \sigma^2}} }{(x_2 - \mu)} d x_2 $$, where $w(z)$ is Faddeeva function, $\nu, \eta, \nu, \mu$ are complex, while $\sigma$ is real. If there is a special function that represents this integral, than the original integral is doable.

Answer 1

$\newcommand{\al}{\alpha}\newcommand{\be}{\beta}\newcommand\la\lambda\newcommand\R{\mathbb R}$By shifting and rescaling, without loss of generality $x_0=0$ and $\sigma=1$. By partial fraction decomposition (and, concerning the denominator $x_1+x_2-\al-\be$, the change of variables from $x_1,x_2$ to, say, $x_1+x_2,x_1-x_2$), we reduce the problem to that of finding the integral \begin{equation} I(t):=I_a(t):=\int_\R dx\,\frac{e^{-itx-x^2/2}}{x-a} \end{equation} for $t\in\R$ and complex $a$ with $\Im a\ne0$.

So, it is enough to find
\begin{equation} J(t):=e^{ita}I(t)=\int_\R dx\,\frac{e^{-it(x-a)-x^2/2}}{x-a} \end{equation} for such $t,a$. Note that
\begin{equation} J(t)\to0 \end{equation} as $t\to-\infty$ (by the Riemann--Lebesgue lemma) and \begin{equation} J'(t)=-i\int_\R dx\,e^{-it(x-a)-x^2/2}=-i\sqrt{2\pi}\,e^{ita-t^2/2}. \end{equation} Thus, \begin{equation} J(t)=\int_{-\infty}^t ds\,J'(s) =i\pi e^{-a^2/2} \Big(\text{erf}\Big(\frac{ia-t}{\sqrt{2}}\Big)-1\Big). \end{equation}
(I have checked these calculations numerically, to make sure I did not miss a constant factor or something like that.)