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?