Can this integral be expanded as a power series? I am looking at this integral:
$$I(x_1^2, x_2^2) = \frac{1}{4}x_1^2 x_2^2 \int\limits_{-\infty}^{\infty} d\tau_3 \int\limits_{-\infty}^{\tau_3} d\tau_4 \int\limits_{-\infty}^{\tau_4} d\tau_5 \int\limits_{-\infty}^{\tau_5} d\tau_6\ (I_{13}I_{25} + I_{15}I_{23})(I_{14}I_{26} + I_{16}I_{24}), \tag{1}$$
with:
$$I_{ij} := \frac{1}{(2\pi)^2} \frac{1}{x_i^2 + \tau_j^2}. \tag{2}$$
I would like to expand $I(1,x_2^2)$ as a power series centered at $x_2 \sim 0^+$, i.e.:
$$I(1,x_2^2) = c_0 + c_1 x_2 + c_2 x_2^2 + ... \tag{3}$$
with $x_2 > 0$. I know from numerical integration that:
$$I(1,0) = \frac{1}{2^{12} \pi^4} = c_0. \tag{4}$$
I would like to find as many coefficients as possible in eq. (3). My naive attempt was to power expand the integrand of eq. (1), which gives:
$$\begin{align}
I(1,x_2^2) = \frac{x_2^2}{4(2\pi)^8} \sum_{k=0}^\infty \sum_{l=0}^\infty (-1)^{k+l} x_2^{2(k+l)} \int\limits_{-\infty}^{\infty} d\tau_3 \int\limits_{-\infty}^{\tau_3} d\tau_4 \int\limits_{-\infty}^{\tau_4} d\tau_5\\ \int\limits_{-\infty}^{\tau_5} d\tau_6\ \left( \frac{\tau_5^{-2(k+1)}}{1+\tau_3^2} + \frac{\tau_3^{-2(k+1)}}{1+\tau_5^2} \right) \left( \frac{\tau_6^{-2(l+1)}}{1+\tau_4^2} + \frac{\tau_4^{-2(l+1)}}{1+\tau_6^2} \right).
\end{align} \tag{5}$$
The idea was to numerically integrate the integrals for given $k$ and $l$. Unfortunately this cannot be right, since there is no $c_0$ term, and the integrals seem to diverge for all $k,l$.
So can such an expansion as a power series be done, and if yes how?
EDIT:
So I have managed to reduce $(1)$ to a one-dimensional integral analytically (recall that I assume $x_1,x_2>0$):
$$\begin{align}
I(x_1^2,x_2^2) = \frac{1}{256\pi^6} \int\limits_{-\infty}^\infty d\tau_3\ \Biggl\lbrace x_1 I_{13} \left(\tan^{-1} \frac{\tau_3}{x_2} \right)^2 \left( 2\tan^{-1} \frac{\tau_3}{x_1} + \pi \right)\\
 + x_2 I_{23} \left(\tan^{-1} \frac{\tau_3}{x_1} \right)^2 \left( 2\tan^{-1} \frac{\tau_3}{x_2} + \pi \right) \Biggr\rbrace. \end{align}\tag{6}$$
From there I should be able to compute the coefficients. For example the coefficient $c_1$ should obey:
$$c_1 = \left. \frac{\partial}{\partial x_2} I(1,x_2^2) \right|_{x_2 = 0}. \tag{7}$$
However when I differentiate $(6)$ and then integrate for decreasing values of $x_2$, the result does not seem to converge. I am moderately confident that this is not a numerical artefact. I observe the same thing for the second derivative.
Any idea why this is happening, or what it could mean?
 A: For $x_1, x_2 > 0$, we have
$$ I(x_1^2,x_2^2) = \frac{F(x_1/x_2) + F(x_2/x_1)}{2^{10}\pi^7}, $$
where
$$ F(x) := \int_{-\infty}^{\infty} \mathrm{d}t \, \frac{\arctan^2(xt)}{t^2+1}. $$
1. The behavior of $F(x)$ as $x \to 0^+$ is easier to study. Indeed, $F(0) = 0$ is obvious, and
$$ \frac{F(x)}{x} \stackrel{(u=xt)}= \int_{-\infty}^{\infty} \mathrm{d}u \, \frac{\arctan^2(u)}{u^2+x^2} \xrightarrow{x \to 0^+} \int_{-\infty}^{\infty} \mathrm{d}u \, \frac{\arctan^2(u)}{u^2} = 2\pi \log 2 $$
shows that
$$F(x) = (2\pi \log 2)x + o(x) \quad \text{as} \quad x \to 0^+.$$
In fact, it can be shown that $F(x)$ extends to an analytic function on all of $\mathbb{C}\setminus(-\infty, 1]$ that agrees with the integral definition along the positive line $(0, \infty)$.
2. The behavior of $F(1/x)$ as $x \to 0^+$ is much trickier. First, we have $F(+\infty) = \frac{\pi^3}{4}$. Then
\begin{align*}
F(1/x) - F(+\infty)
&= 2 \int_{0}^{\infty} \mathrm{d}t \, \frac{\arctan^2(t/x) - (\pi/2)^2}{t^2+1} \\
&= -4x \int_{0}^{\infty} \mathrm{d}t \, \frac{\arctan(t/x)\arctan(t)}{t^2+x^2} \tag{IbP} \\
&= -4x \int_{0}^{\infty} \mathrm{d}t \, \frac{\arctan(t/x)(\arctan(t) - t \mathbf{1}_{[0,1]}(t))}{t^2+x^2} \\
&\quad - 4x \int_{0}^{1} \mathrm{d}t \, \frac{\arctan(t/x)t}{t^2+x^2} .
\end{align*}
In the last line, the first integral without the prefactor $4x$ converges as $x \to 0^+$ by the Dominated Convergence Theorem. Next,
$$ \int_{0}^{1} \mathrm{d}t \, \frac{\arctan(t/x)t}{t^2+x^2}
= \int_{0}^{1/x} \mathrm{d}u \, \frac{\arctan(u)u}{u^2+1}
\sim \frac{\pi}{2} \log(1/x) $$
as $x \to 0^+$, where the asymptotic equivalence in the last step follows from the L'Hopital's Rule. Combining altogether,
$$ F(1/x) = \frac{\pi^3}{4} + (2\pi + o(1)) x \log x \quad\text{as}\quad x \to 0^+. $$
Conclusion. Using the above estimates, we get
$$ I(1,x_2) = \frac{1}{2^{12}\pi^4} + \frac{1 + o(1)}{2^9\pi^6} x_2 \log x_2 \quad\text{as}\quad x_2 \to 0^+. $$
