Integral involving product of arctangent and Gaussian

Question

During some research, it became desirable to compute

$$ \int_{\mathbb{R}^2} \arctan{ (y / x) } \exp{ [ - c ( (x - a)^2 + (y - b)^2 ) ]} d x dy $$

I am aware of the solution in the case when $a = b = 0$ (e.g., this question Distribution of $\arctan(X/Y)$, and, of course, the answer is zero since $\arctan$ is an odd function). However, it would be really fantastic if there ended up being a closed form in terms of arbitrary $a ,b , c \in \mathbb{R}$. I've tried most of the tricks I know (Feynman's trick, various substitutions, etc.) and, when I tried running the computation in Mathematica, it ends up timing out. Would anyone know of a potential strategy or solution?

Update: One approach that seemed promising is as follows. Note that the previous integral may be rewritten as

$$ I(a,b) = \int_{\mathbb{R}^2} \arctan{(y+b)/(x+a)} \exp{- (x^2 + y^2) } dx dy . $$

Observe that $I(a,b)$ satisfies Laplace's equation in the variables $a$ and $b$: $(\partial_a^2 + \partial_b^2) I(a,b) = 0$. We can try to separate variables:

$$ I(a,b) = ( \alpha \cos{(\lambda a)} + \beta \sin{(\lambda a)} ) ( \gamma \cosh{( \lambda b)} + \delta \sinh{( \lambda b)} ) , $$

for some $\alpha , \beta , \gamma, \delta, \lambda \in \mathbb{R}$ that are to be determined. Note $$ I( a , 0) = 0, ~~ \forall a \in \mathbb{R} \implies \gamma ( \alpha \cos{(\lambda a)} + \beta \sin{(\lambda a)} ) = 0 , ~~ \forall a \in \mathbb{R} \implies \gamma = 0 , $$ $$ I(0 , b) = 0, ~~ \forall b \in \mathbb{R} \implies \alpha ( \gamma \cosh{( \lambda b)} + \delta \sinh{( \lambda b)} ) = 0 , \forall b \in \mathbb{R} \implies \alpha = 0 . $$ Hence, \begin{align*} I (a,b) = A \sin{(\lambda a)} \sinh{(\lambda b)} , \end{align*} for some $A$ and $\lambda$ to be determined. However, I am struggling to see what other computable conditions on $I$ are possible to determine $A$ and $\lambda$. There are also some implicit assumptions in this approach regarding interchanging some limits and assuming the solution was separable. Does anyone have any ideas on how to determine $A$ and $\lambda$?

Update 2: Another fact that seems to be useful is

$$ \lim_{a \rightarrow \infty} \arctan{ \frac{y + a}{x + a} } = \frac{\pi}{4} , $$

and the numerical integral of $I(a,a)$ for large choices of $a$ seems to be converging to $\pi^2 / 4$. However, this immediately tells us that the formula $I(a,b) = A \sin{(\lambda a)} \sinh{(\lambda b)}$ isn't correct since $\sinh{(x)}$ diverges for $x \rightarrow \infty$.

Doing some more numerical integration, it seems $I(a,b)$ are odd functions of $a$ and $b$ separately, and an even function of $(a,b)$. Moreover, the value $I(1,1)$ is approximately $1.75221$.

Update 3: I've managed to get far in the calculation of the integral. In fact, I've managed to prove that, to leading order the integral $I(a,b)$ goes like $\pi \arctan{(b/a)}$. Here is how the argument goes.

Start with the formula $$ \int_0^{\infty} \mathrm{erf} ((y+b)s) e^{-(x+a)^2 s^2} d s = \frac{1}{(x+a) \sqrt{\pi}} \arctan{ \frac{y+b}{x+a} } , $$ where $$ \mathrm{erf} (z) = \frac{2}{\sqrt{\pi}} \int_0^z e^{- s^2} d s $$ is the standard error function. This ​can be found in http://keithbriggs.info/documents/erf-integrals.pdf (and also https://nvlpubs.nist.gov/nistpubs/jres/75B/jresv75Bn3-4p149_A1b.pdf). Using this, we write (and compute)

$$ I(a,b) = \sqrt{\pi} \int_0^{\infty} \int_{\mathbb{R}^2} (x+a) \mathrm{erf} ((y+b) s) e^{-(x+a)^2 s^2} e^{-(x^2 + y^2)} dxdyds $$ $$ = \sqrt{\pi} \int_0^{\infty} \left( \int_{\mathbb{R}} \mathrm{erf} ((y+b) s) e^{-y^2} d y \right) \left( \int_{\mathbb{R}} (x+a) e^{-(x+a)^2 s^2} e^{-x^2} d x \right) d s $$ $$ = a \pi \int_0^{\infty} \left( \int_{\mathbb{R}} \mathrm{erf} ((y+b) s) e^{-y^2} d y \right) \left( \frac{1}{(1+s^2)^{3/2}} \exp{ \left( \frac{-a^2 s^2}{1+s^2} \right) } \right) d s . $$

To proceed, we need another formula from https://nvlpubs.nist.gov/nistpubs/jres/75B/jresv75Bn3-4p149_A1b.pdf:

$$ \int_{\mathbb{R}} \mathrm{erf} ((y+b) s) e^{-y^2} d y = \frac{1}{s} \int_{\mathbb{R}} \mathrm{erf}{(y)} \exp{ \left[ - \left( \frac{y}{s} - b \right)^2 \right] } d y = - \sqrt{\pi} \mathrm{erf}{ \left( - \frac{bs}{\sqrt{1 + s^2}} \right) } . $$

Therefore,

$$ I(a,b) = a \pi^{3/2} \int_0^{\infty} \frac{1}{(1+s^2)^{3/2}} \mathrm{erf}{ \left( \frac{bs}{\sqrt{1 + s^2}} \right) } \exp{ \left( \frac{-a^2 s^2}{1+s^2} \right) } d s $$

Differentiating the previous expression with respect to $b$:

$$ \frac{d}{db} I(a,b) = 2 \pi a \int_0^{\infty} \frac{s}{(1+s^2)^2} \exp{ \left( \frac{- (a^2 + b^2) s^2}{1+s^2} \right) } d s = a \pi \frac{1 - e^{- (a^2 + b^2)} }{ (a^2 + b^2)} . $$

Finally, integrating with respect to $b$ yields:

$$ I(a,b) = \pi a \int_0^b \frac{1 - e^{- (a^2 + s^2)} }{a^2 + s^2} d s = \pi \arctan{\frac{b}{a}} - \pi a e^{-a^2} \int_0^b \frac{e^{- s^2} }{a^2 + s^2} d s . $$

I think this is excellent progress and am fairly satisfied. The question now is:

$$ \int_0^b \frac{e^{- s^2} }{a^2 + s^2} d s = ??? $$

Answer 1

Together with the answer to your other question in terms of Owen's T-function we can write the final result as $$I(a,b) = \pi \left[\arctan \left(\frac{b}{a}\right) - 2 \pi \operatorname{T} \left(\sqrt{2} a, \frac{b}{a}\right)\right] \, .$$

In particular, $I(a,b) \sim \frac{\pi b}{a}$ as $a \to \infty$, $$ I(a,\infty) = \frac{\pi^2}{2} \operatorname{erf}(a)$$ and $$ I(a,a) = \frac{\pi^2}{4} \operatorname{erf}^2(a) \, .$$

Answer 2

Recall the formula:

$$ \int_0^{\infty} \frac{1}{s} e^{- (x+a) s} \sin{( (y+b) s )} d s = \arctan{ \frac{y+b}{x+a} } . $$

Then,

$$ I (a,b) = \int_0^{\infty} \frac{1}{s} \int_{\mathbb{R}^2} e^{- (x+a) s} \sin{( (y+b) s )} \exp{ \left[ - (x^2 + y^2) \right] } d x d y d s \\ = \int_0^{\infty} \frac{1}{s} \left( \int_{\mathbb{R}} \sin{( (y+b) s )} e^{-y^2} d x \right) \left( \int_{\mathbb{R}} e^{- (x+a) s} e^{-x^2} d x \right) d s \\ = \pi \int_0^{\infty} \frac{1}{s} e^{-as} \sin{(bs)} d s = \pi \arctan{(b/a)} . $$