$\int\limits_{-\infty}^\infty \left(f(\frac{x-\mu}{1+\Psi/2})-f(\frac{x+\mu}{1-\Psi/2})\right)\frac{x\gamma}{(x-x_{0})^{2}+a}dx$

Question

EDIT: I realized from numerical implementation that the step from \begin{align} \mathcal{I}_2=&\frac{\gamma}{2}\int\limits_{-\infty}^\infty \left(f_T(\frac{x-\mu}{1+\Psi/2})-f_T(\frac{x+\mu}{1-\Psi/2})\right)\left[\frac{1}{x-\left(x_{0}+ \mathrm{i}\frac{\gamma}{2}\right)}+\frac{1}{x-\left(x_{0}- \mathrm{i}\frac{\gamma}{2}\right)}\right]dx\nonumber \end{align} is not the same to the 'equivalent expression' \begin{gather} \mathcal{I}_2=\mathcal{I}^A+\mathcal{I}^R\\ \mathcal{I}^{A/R}=\frac{\gamma}{2}\int\limits_{-\infty}^\infty f_T(z)\left[\frac{1}{z-\frac{x_0-\mu\pm\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2}}-\frac{1}{z-\frac{x_0+\mu\pm\mathrm{i}\frac{\gamma}{2}}{1-\Psi/2}}\right]dz \end{gather} (in fact the first expression gives the black line of the plot and the second gives the red line of the plot). Can someone see the difference?

ORIGINAL QUESTION:I want to solve the following integral

\begin{align} W=\int\limits_{-\infty}^\infty \left(f_T(\frac{x-\mu}{1+\Psi/2})-f_T(\frac{x+\mu}{1-\Psi/2})\right)\frac{x\gamma}{(x-x_{0})^{2}+\frac{\gamma^{2}}{4}}dx \end{align}

with $f_{T}(x)=\frac{1}{1+e^{\frac{x}{T}}}$ the Fermi function and $x_{0},\mu,\gamma, T\in\Re_{\ge0}$ and $\Psi\in[-2,2]$.

I have calculated the analytical solution for the integral in terms of the digammas $\psi(z)$ functions (see derivation below), but when comparing the derived expression with the numerical solution using the quadrature method from the scipy routine in python for the set of parameters $x_{0}=\mu=0,\gamma= T=1$ and $\Psi\in[-2,2]$ there is an evident discrepancy between both approaches, plot at the end. I suspect that my analytical derivation is not completely correct, although when I set $\Psi=0$ the agreement between both approaches is perfect. Can anybody help me to understand where I'm wrong?

The required integral is derived after an auxiliary integral $\mathcal{I}_1$ (correctly checked). Problems may start from $\mathcal{I}_2$.

\begin{align} \mathcal{I}_1=\int\limits_{-\infty}^\infty f_T(x-\mu)\frac{\gamma}{(x-x_{0})^{2}+\frac{\gamma^{2}}{4}}dx \end{align}

Using the substitution $x = \mu + T u$ and the abbreviations $u_0 = \frac{x_0 - \mu}{T}$ and $\beta = \frac{\gamma}{2T}$, we can write the integral as

\begin{align} \mathcal{I}_{1} = 2 \beta \int \limits_{-\infty}^\infty \frac{\mathrm{d} u}{(\mathrm{e}^{u} + 1)[(u-u_0)^2 + \beta^2]} \equiv 2 \beta \int \limits_{-\infty}^\infty g_{\beta,u_0}(u) \, \mathrm{d} u \, . \end{align}

$g_{\beta,u_0}$ has simple poles at $u_0 \pm \mathrm{i} \beta$ with residues

\begin{align} \operatorname{Res}(g_{\beta,u_0},u_0 \pm \mathrm{i} \beta) = \pm \frac{1}{2 \mathrm{i} \beta (\mathrm{e}^{u_0 \pm \mathrm{i} \beta} + 1)} \end{align}

and at $\pm (2n+1) \pi \mathrm{i}$ (so only at the odd-integer multiples of $\pi \mathrm{i}$ !) with residues

\begin{align} \operatorname{Res}(g_{\beta,u_0},\pm (2n+1) \pi \mathrm{i}) = - \frac{1}{[(2n+1) \pi \mathrm{i} \mp u_0]^2 + \beta^2} \end{align}

for $n \in \mathbb{N}_0$ . The integrals of $g_{\beta,u_0}$ along semi-circles (avoiding the poles on the imaginary axis) vanish in the limit of large radii, so we can use the residue theorem to evaluate the integral.

Closing the contour in the upper half-plane (the lower half-plane works just as well) yields

\begin{align} \mathcal{I}_{1} =& 2 \beta \, 2 \pi \mathrm{i} \left[\operatorname{Res}(g_{\beta,u_0},u_0 + \mathrm{i} \beta) + \sum \limits_{n=0}^\infty \operatorname{Res}(g_{\beta,u_0},(2n+1) \pi \mathrm{i}) \right] \\ =&4 \pi \beta\mathrm{i} \left[ \frac{1}{2 \mathrm{i} \beta (\mathrm{e}^{u_0 + \mathrm{i} \beta} + 1)}- \sum \limits_{n=0}^\infty \frac{1}{[(2n+1) \pi \mathrm{i} - u_0]^2 + \beta^2} \right]\\ =&4 \pi \beta\mathrm{i} \left[ \frac{1}{2 \mathrm{i} \beta (\mathrm{e}^{u_0 + \mathrm{i} \beta} + 1)}+ \sum \limits_{n=0}^\infty\frac{\frac{1}{4\pi\beta}}{\left(n+\left[\frac{1}{2}+\frac{\mathrm{i}u_0+\beta}{2\pi}\right]\right)}-\sum \limits_{n=0}^\infty\frac{\frac{1}{4\pi\beta}}{\left(n+\left[\frac{1}{2}+\frac{\mathrm{i}u_0-\beta}{2\pi}\right]\right)} \right] \end{align}

Now we use the series formula for the digamma function $\psi$

\begin{align} \psi(z)=\sum \limits_{n=0}^\infty\left(\frac{1}{n+1}-\frac{1}{n+z}\right)-\gamma \end{align}

to find

\begin{align} \mathcal{I} = 2 \pi \left[\frac{1}{\mathrm{e}^{u_0 + \mathrm{i} \beta} + 1} - \frac{1}{2 \pi \mathrm{i}} \left(\psi \left(\frac{1}{2} + \frac{\beta + \mathrm{i} u_0}{2 \pi}\right) - \psi \left(\frac{1}{2} + \frac{-\beta + \mathrm{i} u_0}{2 \pi}\right) \right)\right]\end{align}

Finally, we apply the reflection formula

\begin{align} \psi(1-z)=\psi(z)+\pi\cot{\left(\pi z\right)} \end{align}

to the second digamma function together with the properties

\begin{gather} \operatorname{Im} \left[\psi(x+ \mathrm{i}y)\right]=\frac{ \mathrm{i}}{2}\left(\psi(x- \mathrm{i}y)-\psi(x+ \mathrm{i}y)\right)\\ \cot{\left(\frac{\pi}{2}+z\right)}=-\tan(z)\\ \frac{1}{e^{z}+1}+\frac{\mathrm{i}}{2}\tan\left(\frac{-\mathrm{i}z}{2}\right)=\frac{1}{2} \end{gather}

and simplify the result: \begin{align} \mathcal{I}&= 2 \pi \left[\frac{1}{\mathrm{e}^{u_0 + \mathrm{i} \beta} + 1} - \frac{1}{2 \mathrm{i}} \tan\left(\frac{\beta - \mathrm{i} u_0}{2}\right) -\frac{1}{2 \pi \mathrm{i}} \left(\psi \left(\frac{1}{2} + \frac{\beta + \mathrm{i} u_0}{2 \pi}\right) - \psi \left(\frac{1}{2} + \frac{\beta - \mathrm{i} u_0}{2 \pi}\right) \right)\right] \\ &= 2 \pi \left[\frac{1}{2} - \frac{1}{\pi} \operatorname{Im} \left(\psi \left(\frac{1}{2} + \frac{\beta + \mathrm{i} u_0}{2 \pi}\right) \right) \right] = \pi - 2 \operatorname{Im} \left[\psi \left(\frac{1}{2} + \frac{\beta + \mathrm{i} u_0}{2 \pi}\right) \right] \, . \end{align}

Returning to the original parameters, we end up with

\begin{align} \boxed{\mathcal{I}_{1}(\gamma,x_{0}, \mu,T) =\int \limits_{-\infty}^\infty \frac{1}{\mathrm{e}^{(x - \mu)/T} + 1} \frac{\gamma}{(x-x_0)^2 + \gamma^2/4} \, \mathrm{d} x = \pi - 2 \operatorname{Im} \left[\psi \left(\frac{1}{2} + \frac{\frac{\gamma}{2} + \mathrm{i} (x_0 - \mu)}{2 \pi T}\right) \right] }\, . \end{align}

Note that the dependence on the parameters is explicit and in the special case $x_0 = \mu$ the result is simply $\pi$. This can also be shown using elementary methods, so the implicit assumption $u_0 + \mathrm{i} \beta \not\in (2 \mathbb{Z} + 1) \pi \mathrm{i}$ used in the computation of the residues is justified.

The integral we are interested in is

\begin{align} \mathcal{I}_2=\int\limits_{-\infty}^\infty \left(f_T(\frac{x-\mu}{1+\Psi/2})-f_T(\frac{x+\mu}{1-\Psi/2})\right)\frac{x\gamma}{(x-x_{0})^{2}+\frac{\gamma^{2}}{4}}dx \label{eq:2} \end{align} We can rewrite the integral $I_{2}$ decomposing the second factor as the sum of two complex parts. This is related to the advanced and retarded Green's functions

\begin{align} \frac{x\gamma}{(x-x_{0})^{2}+\frac{\gamma^{2}}{4}} = &\frac{\gamma}{2}\left[\frac{1}{x-\left(x_{0}+ \mathrm{i}\frac{\gamma}{2}\right)}+\frac{1}{x-\left(x_{0}- \mathrm{i}\frac{\gamma}{2}\right)}\right] - \mathrm{i}x_{0}\left[\frac{1}{x-\left(x_{0}+ \mathrm{i}\frac{\gamma}{2}\right)}-\frac{1}{x-\left(x_{0}- \mathrm{i}\frac{\gamma}{2}\right)}\right]\\ =&\frac{\gamma}{2}\left[G^{A}+G^{R}\right]- \mathrm{i}x_{0}\left[G^{A}-G^{R}\right] \end{align}

On the other hand, we can rewrite the Cauchy distribution used in the previous integral as function of the Green's functions too and using that the second contribution is related to the previously calculated integral $\mathcal{I}_{1}$

\begin{align} \frac{\gamma}{(x-x_{0})^{2}+\frac{\gamma^{2}}{4}} =- \mathrm{i}\left[\frac{1}{x-\left(x_{0}+ \mathrm{i}\frac{\gamma}{2}\right)}-\frac{1}{x-\left(x_{0}- \mathrm{i}\frac{\gamma}{2}\right)}\right] =- \mathrm{i}\left[G^{A}-G^{R}\right]. \end{align}

So

\begin{align} \mathcal{I}_2=&\frac{\gamma}{2}\int\limits_{-\infty}^\infty \left(f_T(\frac{x-\mu}{1+\Psi/2})-f_T(\frac{x+\mu}{1-\Psi/2})\right)\left[\frac{1}{x-\left(x_{0}+ \mathrm{i}\frac{\gamma}{2}\right)}+\frac{1}{x-\left(x_{0}- \mathrm{i}\frac{\gamma}{2}\right)}\right]dx\nonumber\\&+x_{0}\left(I_{1}\left(\gamma,x_{0}, \mu,T(1+\Psi/2)\right)-I_{1}(\gamma,x_{0}, -\mu,T(1-\Psi/2))\right)\\ =&\frac{\gamma}{2}\int\limits_{-\infty}^\infty f_T(z)\left[\left(\frac{1}{z-\frac{x_0-\mu+\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2}}+\frac{1}{z-\frac{x_0-\mu-\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2}}\right)-\left(\frac{1}{z-\frac{x_0+\mu+\mathrm{i}\frac{\gamma}{2}}{1-\Psi/2}}+\frac{1}{z-\frac{x_0+\mu-\mathrm{i}\frac{\gamma}{2}}{1-\Psi/2}}\right)\right]dz\nonumber\\&+x_{0}\left(\mathcal{I}_{1}\left(\gamma,x_{0}, \mu,T_L\right)-\mathcal{I}_{1}(\gamma,x_{0}, -\mu,T_R)\right) \end{align}

Note that the integral $\mathcal{I}_2$ is convergent because the difference of the Fermi functions contributes asymptotically as $x^{-1}$ and the Green's function contributes another $x^{-1}$. The integral of a single Green's function with a single Fermi function would be logarithmically divergent. In the second equality the Fermi functions arguments have been shifted to the Green's functions terms we can split the integral in two parts such that

\begin{gather} \mathcal{I}_2=\mathcal{I}^A+\mathcal{I}^R+x_{0}\left(\mathcal{I}_{1}\left(\gamma,x_{0}, \mu,T_L\right)-\mathcal{I}_{1}(\gamma,x_{0}, -\mu,T_R)\right)\\ \mathcal{I}^{A/R}=\frac{\gamma}{2}\int\limits_{-\infty}^\infty f_T(z)\left[\frac{1}{z-\frac{x_0-\mu\pm\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2}}-\frac{1}{z-\frac{x_0+\mu\pm\mathrm{i}\frac{\gamma}{2}}{1-\Psi/2}}\right]dz=\frac{\gamma}{2}\int\limits_{-\infty}^\infty g^{A/R}(z)dz \end{gather}

closing the integration contour in the upper half and using now the series relation for the Fermi function

\begin{align} f_{T}(z)=\frac{1}{1+e^{\frac{z}{T}}} =\frac{1}{2}- \frac{\mathrm{i}}{2\pi}\sum \limits_{n=-\infty}^\infty \frac{1}{n+\frac{1}{2} +\mathrm{i}\frac{z}{2\pi T}} \end{align}

we have

\begin{gather} \mathcal{I}^{A}=\frac{\gamma}{2}2\pi\mathrm{i}\left( \operatorname{Res}(g^{A},\frac{x_0-\mu+\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2})-\operatorname{Res}(g^{A},\frac{x_0+\mu+\mathrm{i}\frac{\gamma}{2}}{1-\Psi/2})+\sum\limits_{n=0}^\infty \operatorname{Res}(g^{A},(2n+1)\pi\mathrm{i}T)\right)\\ =\gamma\pi\mathrm{i}\left(f_{T}(\frac{x_0-\mu+\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2})-f_{T}(\frac{x_0+\mu+\mathrm{i}\frac{\gamma}{2}}{1-\Psi/2})-T\sum\limits_{n=0}^\infty\frac{1}{(2n+1)i\pi T-\frac{x_{0}-\mu+i\gamma/2}{1+\Psi/2}}+T\sum\limits_{n=0}^\infty\frac{1}{(2n+1)i\pi T-\frac{x_{0}+\mu+i\gamma/2}{1-\Psi/2}}\right)\\ =\gamma\pi\mathrm{i}\left(f_{T}(\frac{x_0-\mu+\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2})-f_{T}(\frac{x_0+\mu+\mathrm{i}\frac{\gamma}{2}}{1-\Psi/2})+\frac{i}{2\pi}\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{-\gamma/2+i(x_{0}-\mu)}{2\pi T(1+\Psi/2)}}-\frac{i}{2\pi}\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{-\gamma/2+i(x_{0}+\mu)}{2\pi T(1-\Psi/2)}}\right)\\ =\gamma\pi\mathrm{i}\frac{i}{2\pi}\left(-\sum \limits_{n=-\infty}^\infty \frac{1}{n+\frac{1}{2} +\mathrm{i}\frac{\frac{x_0-\mu+\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2}}{2\pi T}}+\sum \limits_{n=-\infty}^\infty \frac{1}{n+\frac{1}{2} +\mathrm{i}\frac{\frac{x_0+\mu+\mathrm{i}\frac{\gamma}{2}}{1+\Psi/2}}{2\pi T}}+\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{-\gamma/2+i(x_{0}-\mu)}{2\pi T(1+\Psi/2)}}-\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{-\gamma/2+i(x_{0}+\mu)}{2\pi T(1-\Psi/2)}}\right)\\ =\frac{\gamma}{2\pi}\left(\sum \limits_{n=-\infty}^\infty \frac{1}{n+\frac{1}{2} +\frac{-\gamma/2+i(x_{0}-\mu)}{2\pi T_L}}-\sum \limits_{n=-\infty}^\infty \frac{1}{n+\frac{1}{2} +\frac{-\gamma/2+i(x_{0}+\mu)}{2\pi T_R}}-\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{-\gamma/2+i(x_{0}-\mu)}{2\pi T_L}}+\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{-\gamma/2+i(x_{0}+\mu)}{2\pi T_R}}\right)\\ =\frac{\gamma}{2}\left(\sum \limits_{n=-\infty}^{-1} \frac{1}{n+\frac{1}{2} +\frac{-\gamma/2+i(x_{0}-\mu)}{2\pi T_L}}-\sum \limits_{n=-\infty}^{-1} \frac{1}{n+\frac{1}{2} +\frac{-\gamma/2+i(x_{0}+\mu)}{2\pi T_R}}\right)\\ =\frac{\gamma}{2}\left(\sum \limits_{n=1}^{\infty} \frac{-1}{n-\frac{1}{2} +\frac{\gamma/2-i(x_{0}-\mu)}{2\pi T_L}}+\sum \limits_{n=1}^{\infty} \frac{1}{n-\frac{1}{2} +\frac{\gamma/2-i(x_{0}+\mu)}{2\pi T_R}}\right)\\ =\frac{\gamma}{2}\left(\psi(\frac{1}{2} +\frac{\gamma/2-i(x_{0}-\mu)}{2\pi T_L})-\psi(\frac{1}{2} +\frac{\gamma/2-i(x_{0}+\mu)}{2\pi T_R})\right) \end{gather}

with $T_{L/R}=T(1\pm\Psi/2)$. In the last step we have used

\begin{align} \psi(z+1)=-\gamma+\sum\limits_{n=1}^\infty\left(\frac{1}{n}-\frac{1}{n+z}\right) \end{align}

for the second contribution

\begin{gather} \mathcal{I}^{R}=\frac{\gamma}{2}2\pi\mathrm{i}\left(\sum\limits_{n=0}^\infty \operatorname{Res}(g^{R},(2n+1)\pi\mathrm{i}T)\right)\\ =\gamma\pi\mathrm{i}\left(\frac{i}{2\pi}\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{\gamma/2+i(x_{0}-\mu)}{2\pi T(1+\Psi/2)}}-\frac{i}{2\pi}\sum\limits_{n=0}^\infty\frac{1}{n+\frac{1}{2}+\frac{\gamma/2+i(x_{0}+\mu)}{2\pi T)}}\right)\\ =\frac{\gamma}{2}\left(\psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}-\mu)}{2\pi T_L})-\psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}+\mu)}{2\pi T_R})\right) \end{gather}

So

\begin{gather} \mathcal{I}^{A}+\mathcal{I}^{R}=\frac{\gamma}{2}\left(\psi(\frac{1}{2}+\frac{\gamma/2-i(x_{0}-\mu)}{2\pi T_L})-\psi(\frac{1}{2}+\frac{\gamma/2-i(x_{0}+\mu)}{2\pi T_R})+\psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}-\mu)}{2\pi T_L})-\psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}+\mu)}{2\pi T_R})\right)\\ =\gamma\operatorname{Re} \left(\psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}-\mu)}{2\pi T_L})\right)-\gamma\operatorname{Re} \left( \psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}+\mu)}{2\pi T_R}\right) \end{gather}

Finally

\begin{gather} \mathcal{I}_2=\gamma\operatorname{Re} \left(\psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}-\mu)}{2\pi T_L})\right)-\gamma\operatorname{Re} \left( \psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}+\mu)}{2\pi T_{R}})\right)+x_{0}\left(\mathcal{I}_{1}\left(\gamma,x_{0}, \mu,T_L\right)-\mathcal{I}_{1}(\gamma,x_{0}, -\mu,T_R)\right) \end{gather} So

\begin{gather} \boxed{W=\gamma\operatorname{Re} \left(\psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}-\mu)}{2\pi T_L})\right)-\gamma\operatorname{Re} \left( \psi(\frac{1}{2}+\frac{\gamma/2+i(x_{0}+\mu)}{2\pi T_{R}})\right)+x_{0}\left(\mathcal{I}_{1}\left(\gamma,x_{0}, \mu,T_L\right)-\mathcal{I}_{1}(\gamma,x_{0}, -\mu,T_R)\right)} \end{gather}

(The motivation is related to a quantum problem, but people downvoted previous question arguing that this is related to physics, so if somebody wants me to explain the physical interpretation of the parameters please do not hesitate.)

Answer 1

Solution by @ComplexYetTrivial in math.stackexchange: The step in the edit is indeed the only mistake, the rest of your calculation seems to be correct. In order to see what the problem is we write $$ \mathcal{J}^{A/R} = \frac{\gamma}{2} \int \limits_{-\infty}^{\infty} \left[f_T \left(\frac{x - \mu}{1 + \Psi/2}\right) - f_T \left(\frac{x + \mu}{1 - \Psi/2}\right)\right] \, \frac{\mathrm{d} x}{x - x_0 \mp \mathrm{i} \frac{\gamma}{2}} \, .$$ You correctly noted that these integrals are convergent due to the difference of the two Fermi functions. If there was only one Fermi function, they would diverge logarithmically at the lower limit.

However, you then proceed to perform two different changes of variables ($x = \pm \mu + (1 \pm \Psi/2) z$) for the first and the second term. This amounts to splitting the convergent integral into two divergent integrals, doing the substitutions and then combining them again, which is not allowed and changes the result!

There are two simple ways to fix this:

1. We can integrate by parts to obtain an integral with two convergent parts, split it up, do the substitutions and then reverse the integration by parts (the missing term will appear in the last step).
2. We can introduce a variable but finite lower limit (thus making both parts of the integral convergent individually), follow your steps and send the lower limit to infinity at the end.

I will show the second possibility here, since it illustrates nicely how the error arises. We have \begin{align} \mathcal{J}^{A/R} &= \lim_{r \to \infty} \frac{\gamma}{2} \int \limits_{-r}^{\infty} \left[f_T \left(\frac{x - \mu}{1 + \Psi/2}\right) - f_T \left(\frac{x + \mu}{1 - \Psi/2}\right)\right] \, \frac{\mathrm{d} x}{x - x_0 \mp \mathrm{i} \frac{\gamma}{2}} \\ &= \lim_{r \to \infty} \frac{\gamma}{2} \left[~\int \limits_{-r}^{\infty} \, \frac{f_T \left(\frac{x - \mu}{1 + \Psi/2}\right)}{x - x_0 \mp \mathrm{i} \frac{\gamma}{2}} \, \mathrm{d} x - \int \limits_{-r}^{\infty} \, \frac{f_T \left(\frac{x + \mu}{1 - \Psi/2}\right)}{x - x_0 \mp \mathrm{i} \frac{\gamma}{2}} \, \mathrm{d} x \right] \\ &= \lim_{r \to \infty} \frac{\gamma}{2} \left[~\int \limits_{- \frac{r + \mu}{1 + \Psi/2}}^{\infty} \, \frac{f_T \left(z\right)}{z - \frac{x_0 - \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 + \Psi/2}} \, \mathrm{d} z - \int \limits_{- \frac{r - \mu}{1 - \Psi/2}}^{\infty} \, \frac{f_T \left(z\right)}{z - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}} \, \mathrm{d} z\right] \\ &= \lim_{r \to \infty} \frac{\gamma}{2} \left[~\int \limits_{- \frac{r + \mu}{1 + \Psi/2}}^{\infty} f_T (z) \left(\frac{1}{z - \frac{x_0 - \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 + \Psi/2}} - \frac{1}{z - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}}\right) \, \mathrm{d} z + \int \limits_{- \frac{r + \mu}{1 + \Psi/2}}^{- \frac{r - \mu}{1 - \Psi/2}} \frac{f_T \left(z\right)}{z - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}} \, \mathrm{d} z \right]. \end{align} Letting $r \to \infty$ in the first integral we simply obtain your result ($\mathcal{I}^{A/R}$). The second term is the correction we are looking for: \begin{align} \mathcal{J}^{A/R} - \mathcal{I}^{A/R} &= \lim_{r \to \infty} \frac{\gamma}{2} \int \limits_{- \frac{r + \mu}{1 + \Psi/2}}^{- \frac{r - \mu}{1 - \Psi/2}} \frac{f_T \left(z\right)}{z - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}} \, \mathrm{d} z \\ &= \lim_{r \to \infty} \frac{\gamma}{2} \int \limits_{- \frac{r + \mu}{1 + \Psi/2}}^{- \frac{r - \mu}{1 - \Psi/2}} \frac{\mathrm{d} z}{z - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}} - \lim_{r \to \infty} \frac{\gamma}{2} \int \limits_{- \frac{r + \mu}{1 + \Psi/2}}^{- \frac{r - \mu}{1 - \Psi/2}} \frac{1 - f_T \left(z\right)}{z - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}} \, \mathrm{d} z \, . \end{align} Since $1 - f_T(z) \leq \mathrm{e}^{z/T}$, the second limit vanishes and we are left with \begin{align} \mathcal{J}^{A/R} - \mathcal{I}^{A/R} &= \lim_{r \to \infty} \frac{\gamma}{2} \int \limits_{- \frac{r + \mu}{1 + \Psi/2}}^{- \frac{r - \mu}{1 - \Psi/2}} \frac{\mathrm{d} z}{z - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}} = \lim_{r \to \infty} \frac{\gamma}{2} \log \left(\frac{- \frac{r - \mu}{1 - \Psi/2} - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}}{- \frac{r + \mu}{1 + \Psi/2} - \frac{x_0 + \mu \pm \mathrm{i} \frac{\gamma}{2}}{1 - \Psi/2}}\right) \\ &= \frac{\gamma}{2} \log \left(\frac{1 + \Psi/2}{1-\Psi/2}\right) = \gamma \operatorname{artanh} \left(\frac{\Psi}{2}\right) \, . \end{align}

With this addition the correct final result is \begin{align} W &= \int\limits_{-\infty}^\infty \left[f_T\left(\frac{x-\mu}{1+\Psi/2}\right)-f_T\left(\frac{x+\mu}{1-\Psi/2}\right)\right] \frac{\gamma x}{(x-x_0)^2 + \gamma^2/4} \, \mathrm{d} x \\ &=2 \gamma \operatorname{artanh} \left(\frac{\Psi}{2}\right) \!+ \!(\gamma \operatorname{Re} {} \! -\! 2 x_0 \operatorname{Im}) \left[\operatorname{\psi} \left(\frac{1}{2} \! + \! \frac{\gamma/2 \! + \! \mathrm{i} (x_0 \! - \! \mu)}{2 \pi T (1 \! + \! \Psi/2)}\right) \! - \! \operatorname{\psi} \left(\frac{1}{2} \! + \! \frac{\gamma/2 \! + \! \mathrm{i} (x_0 \! + \! \mu)}{2 \pi T (1 \! - \! \Psi/2)} \right) \right] . \end{align}