# Solve $I=\int_{0}^{1}\int_{0}^{1}\frac{xy}{\sqrt{x^2+y^2}}\arctan{\sqrt{x^2+y^2}}\arctan{\frac{y^2}{x^2}}dxdy$

I tried to solve this integral by using polar coordinate system, but it becomes harder: $$I=\int_{0}^{1}\int_{0}^{1}\frac{xy}{\sqrt{x^2+y^2}}\arctan{\sqrt{x^2+y^2}}\arctan{\frac{y^2}{x^2}}dxdy$$ After transformation $$I=2\int_{0}^{\frac{\pi}{4}}\arctan{(\tan^2\theta)}\sin\theta \cos\theta d\theta\int_{0}^{\frac{1}{cos\theta}}r^2\arctan{(r)}dr$$ And after IBP, $$\int r^2\arctan{(r)}dr=\frac{1}{6}\left(2r^3\arctan(r)-r^2+\ln(r^2+1)\right)$$ and replace upper and lower limits, it becomes a mess. Did I misdirect? Need some help from everyone, thank you very much.

The integrand function is not symmetric with respect to the line $$y=x$$, so $$I$$ is not two times the integral over the triangle $$\{(x,y):0\leq x\leq y\leq 1\}$$. On the other hand, note that \begin{align*} I&=\int_{0}^{1}\int_{0}^{x}\frac{xy}{\sqrt{x^2+y^2}}\arctan{\sqrt{x^2+y^2}}\arctan\left(\frac{y^2}{x^2}\right)\,dydx\\ &\qquad+\int_{0}^{1}\int_{0}^{y}\frac{xy}{\sqrt{x^2+y^2}}\arctan{\sqrt{x^2+y^2}}\left(\frac{y^2}{x^2}\right)\,dxdy\\ &=\int_{0}^{1}\int_{0}^{x}\frac{xy}{\sqrt{x^2+y^2}}\arctan{\sqrt{x^2+y^2}}\left(\arctan\left(\frac{y^2}{x^2}\right)+\arctan\left(\frac{y^2}{x^2}\right)\right)dydx. \end{align*} Now recall that $$\arctan(t)+\arctan(1/t)=\pi/2$$ for $$t>0$$, therefore \begin{align*} I&=\frac{\pi}{2}\int_{0}^{1}\int_{0}^{x}\frac{xy}{\sqrt{x^2+y^2}}\arctan{\sqrt{x^2+y^2}}dydx\\ &=\frac{\pi}{2}\int_{0}^{\pi/4}\cos(\theta)\sin(\theta)\left(\int_{0}^{1/\cos(\theta)}r^2\arctan(r)\,dr\right)d\theta\\ &=\frac{\pi}{2}\int_{0}^{\pi/4}\cos(\theta)\sin(\theta)\left[\frac{1}{6}\left(2r^3\arctan(r)-r^2+\ln(r^2+1)\right)\right]_{0}^{1/\cos(\theta)}\,d\theta\\ &=\frac{\pi}{12}\int_{1}^{\sqrt{2}}\left(\frac{\ln(1+t^2)}{t^3}+2\arctan(t)-\frac{1}{t}\right)\,dt\\ &=\frac{\pi}{48}\left(8\ln(2)-7\ln(3)+8\sqrt{2}\arctan(\sqrt{2})-2\pi\right) \end{align*} where $$t=1/\cos(\theta)$$.
For the second part, we have also another way which seems to be easier \begin{align*} I&=\frac{\pi}{2}\int_{0}^{1}\int_{0}^{x}\frac{xy}{\sqrt{x^2+y^2}}\arctan{\sqrt{x^2+y^2}}dydx\\ &=\frac{\pi}{2}\int_{0}^{1}r^2\arctan(r)\left(\int_{0}^{\pi/4}\cos(\theta)\sin(\theta)d\theta\right)\,dr\\ &\qquad+\frac{\pi}{2}\int_{1}^{\sqrt{2}}r^2\arctan(r)\left(\int_{\arccos(1/r)}^{\pi/4}\cos(\theta)\sin(\theta)d\theta\right)\,dr\\ &=\frac{\pi}{2}\int_{0}^{1}r^2\arctan(r)\left[-\frac{\cos^2(\theta)}{2}\right]_{0}^{\pi/4}\,dr+\frac{\pi}{2}\int_{1}^{\sqrt{2}}r^2\arctan(r)\left[-\frac{\cos^2(\theta)}{2}\right]_{\arccos(1/r)}^{\pi/4}\,dr\\ &=\frac{\pi}{8}\int_{0}^{1}r^2\arctan(r)\,dr+\frac{\pi}{8}\int_{1}^{\sqrt{2}}\arctan(r)\left(2-r^2\right)\,dr\\ &=\frac{\pi}{48}\left(8\ln(2)-7\ln(3)+8\sqrt{2}\arctan(\sqrt{2})-2\pi\right). \end{align*}