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 = ??? $$
Best Answer
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) \, .$$