Consider a related integral $\displaystyle\;J(a) = \int_0^\infty e^{-ax}\tan^{-1}(x) dx\;$.
We can integrate it by part to get:
$$J(a) = \frac{1}{a}\int_0^\infty e^{-ax}\frac{dx}{1+x^2}$$
Using this, we can rewrite our integral as
$$I(a)
= \int_0^\infty \left(\sum_{n=1}^\infty e^{-nax}\right)\tan^{-1}(x)dx
= \sum_{n=1}^\infty J(na)
= \sum_{n=1}^\infty \frac{1}{na}\int_0^\infty e^{-nax}\frac{dx}{1+x^2}$$
Replace $x$ by $x/n$ in each terms, this can be transformed as
$$I(a) = \frac{1}{a}\int_0^\infty e^{-at}\left(\sum_{n=1}^\infty \frac{1}{t^2+n^2}\right) dt
$$
Start with a infinite product expansion of $\displaystyle\;\frac{\sinh(\pi x)}{\pi x}\;$. By taking logarithm and differentiate,
$$\frac{\sinh(\pi x)}{\pi x} = \prod_{n=1}^\infty \left(1 + \frac{x^2}{n^2}\right)
\quad\implies\quad
\sum_{n=1}^\infty \frac{1}{x^2 + n ^2 } = \frac{1}{2x}\left( \pi\coth(\pi x) - \frac{1}{x}\right)
$$
we can rewrite the sum inside above integrand of $I(a)$ and obtain
$$I(a) = \frac{\pi}{2a}F\left(\frac{a}{\pi}\right)
\quad\text{ where }\quad
F(x) = \int_0^\infty e^{-xt} \left(\coth(t) - \frac{1}{t}\right) \frac{dt}{t}
$$
Notice
$$-\frac{dF(x)}{dx} = \int_0^\infty e^{-xt}\left(\coth(t) - \frac{1}{t}\right)dt$$
and compare this with a integral representation of digamma function
$\displaystyle\;\psi(x) = \frac{d\log\Gamma(x)}{dx}\;$,
$$\psi(x) = \int_0^\infty\left(\frac{e^{-t}}{t} - \frac{e^{-xt}}{1-e^{-t}}\right)dt
\quad\implies\quad
\psi\left(\frac{x}{2}\right) =
\int_0^\infty\left(\frac{e^{-2t}}{t} - \frac{2e^{-xt}}{1-e^{-2t}}\right)dt
$$
We get
$$\begin{array}{rrl}
&-\frac{dF(x)}{dx}
&= \int_0^\infty \frac{e^{-2t} - e^{-xt}}{t}dt - \frac12
\left[\psi\left(\frac{x}{2}\right) + \psi\left(\frac{x}{2}+1\right)\right]\\
&&= \log\frac{x}{2} - \frac12\left[\psi\left(\frac{x}{2}\right) + \psi\left(\frac{x}{2}+1\right)\right]\\
\implies &
F(x) &= \text{const.} + \log\Gamma\left(\frac{x}{2}\right) + \log\Gamma\left(\frac{x}{2}+1 \right) - x\log\left(\frac{x}{2e}\right)
\end{array}$$
Since $\lim\limits_{x\to\infty}F(x) = 0$, we can use Stirling's approximation to fix the constant in $F(x)$ to $-\log(2\pi)$. As a result, We find
$$I(a) = \frac{\pi}{2a}\left[ 2\log\Gamma\left(\frac{a}{2\pi}+1\right)
- \log(a)\right] -\frac12\log\left(\frac{a}{2\pi e}\right)
$$
For $\alpha, \beta \geq 0$ and $\gamma > 0$ let
$$ f_{\alpha,\beta,\gamma} (z) = \frac{\exp \left(\mathrm{i} z \frac{z^2 - \alpha^2}{z^2 - \beta^2}\right)}{z^2 + \gamma^2} \, , \, z \in \mathbb{C} \setminus \{\pm\beta, \pm \mathrm{i} \gamma\} \, . $$
By symmetry we have
$$ I (\alpha,\beta,\gamma) = \frac{1}{2} \int \limits_{-\infty}^\infty f_{\alpha,\beta,\gamma} (x) \, \mathrm{d} x \, .$$
Clearly, the integral exists for every combination of parameters and is bounded by $\frac{\pi}{2 \gamma}$ .
We can compute its value using the residue theorem. First note that we have
$$ \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) = \frac{1}{2 \mathrm{i} \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) $$
and (by Jordan's lemma)
$$ \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z = 0 \, ,$$
where $\Gamma_R$ is a semi-circle of radius $R > 0$ in the upper half-plane. A naive application of the residue theorem therefore yields
$$ I (\alpha,\beta,\gamma) = \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z \right]= \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) \, .$$
While the result is indeed correct, this calculation is only valid for $\alpha = \beta$ . In this case we obtain (as already mentioned in the question) $I(\alpha,\alpha,\gamma) = \frac{\pi}{2 \gamma}\mathrm{e}^{- \gamma}$ .
If $\alpha \neq \beta$ , $f_{\alpha,\beta,\gamma}$ has essential singularities on the real axis, namely at $\pm \beta$ . Thus we need to deform our contour using small semi-circles $\gamma_\varepsilon (\pm \beta)$ of radius $\varepsilon > 0$ centred at these points and show that their contribution to the integral vanishes as $\varepsilon \to 0$ .
For simplicity, I will only discuss the case $\beta = 0$ in detail. We want to find
$$ \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z = \varepsilon \int \limits_0^\pi \frac{\mathrm{e}^{\mathrm{i} \left[\phi + \varepsilon \mathrm{e}^{\mathrm{i}\phi} - \alpha^2 \varepsilon^{-1} \mathrm{e}^{-\mathrm{i} \phi}\right]}}{\gamma^2+\varepsilon^2 \mathrm{e}^{2 \mathrm{i} \phi}} \, \mathrm{d} \phi = \frac{\varepsilon}{\gamma^2} \int \limits_0^\pi \mathrm{e}^{\mathrm{i} \left[\phi - \alpha^2 \varepsilon^{-1} \mathrm{e}^{-\mathrm{i} \phi}\right]} \left[1 + \mathcal{O} (\varepsilon) \right] \, \mathrm{d} \phi \, .$$
The leading-order term can actually be calculated analytically: for $z \in \mathbb{C}$ we have
$$ \int \limits_0^\pi \mathrm{e}^{\mathrm{i} \left[\phi - z \mathrm{e}^{-\mathrm{i} \phi}\right]} \, \mathrm{d} \phi = \mathrm{i} \left[(2 \operatorname{Si}(z) - \pi) z + 2 \cos(z)\right] \, . $$
The sine integral $\operatorname{Si}$ satisfies $\lim_{x \to \infty} \operatorname{Si}(x) = \frac{\pi}{2}$, so we obtain
$$ \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z \sim \frac{\mathrm{i} \varepsilon}{\gamma^2} \left[(2 \operatorname{Si}(\alpha^2 \varepsilon^{-1}) - \pi) \alpha^2 \varepsilon^{-1} + 2 \cos(\alpha^2 \varepsilon^{-1})\right] \stackrel{\varepsilon \to 0}{\longrightarrow} 0 $$
as desired. Now we are allowed to apply the residue theorem , which yields
\begin{align}
I(\alpha,0,\gamma) &= \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,0,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z \right] \\
&= \frac{\pi}{2 \gamma} \exp \left(- \frac{\alpha^2 + \gamma^2}{\gamma}\right) \, .
\end{align}
Almost the same calculation yields
\begin{align}
I(\alpha,\beta,\gamma) &= \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (\beta)} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (-\beta)} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z \right] \\
&= \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right)
\end{align}
for $\beta > 0$ .
Combining these results we conclude that the integral is given by
$$ I(\alpha,\beta,\gamma) = \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) $$
for every $\alpha,\beta \geq 0$ and $\gamma > 0$ .
Best Answer
Yes, you can apply Differentiation Under the Integral Sign as I show you bellow. Just a quick note, I'm applying a slightly modified version of Leibniz Rule, but in essence it's pretty much the same.
$$f(\alpha,\beta)=\int_{0}^{\infty}\frac{\arctan(\alpha x)\arctan(\beta x)}{x^2}dx$$
$$=\int_{0}^{\infty}\int_{0}^{\alpha}\int_{0}^{\beta}\frac{dydzdx}{(1+y^2x^2)(1+z^2x^2)}=\int_{0}^{\alpha}\int_{0}^{\beta}dydz\int_{0}^{\infty}\frac{dx}{(1+y^2x^2)(1+z^2x^2)}$$
$$=\int_{0}^{\alpha}\int_{0}^{\beta}\frac{dydz}{y^2-z^2}\int_{0}^{\infty}\left(\frac{y^2}{1+y^2x^2}-\frac{z^2}{1+z^2x^2}\right)dx$$
$$=\int_{0}^{\alpha}\int_{0}^{\beta}dydz\left[\frac{y\arctan(yx)-z\arctan(zx)}{(y-z)(y+z)}\right]^{\infty}_{0}=\frac{\pi}{2}\int_{0}^{\alpha}\int_{0}^{\beta}\frac{dydz}{y+z}$$
$$=\frac{\pi}{2}\int_{0}^{\alpha}dz[log(z+\beta)-log(z)]=\frac{\pi}{2}[(\beta+z)log(\beta+z)-z-zlog(z)+z]^{\alpha}_{0}$$
$$\boxed{f(\alpha,\beta)=\int_{0}^{\infty}\frac{\arctan(\alpha x)\arctan(\beta x)}{x^2}dx=\frac{\pi}{2}log\left(\frac{\left(\alpha+\beta\right)^{\alpha+\beta}}{\alpha^\alpha\beta^\beta}\right)}$$