I want to numerically calculate the low temperature limit of the following Matsubara sum
$$ S(\Omega) = \pi T \sum_{\omega_n} \frac{4\Delta^2+\Omega^2}{s_1 s_2 (s_1 + s_2)},$$ with $\omega_n = \pi T (2n+1)$, $s_1(\omega_n) = \sqrt{\omega_n^2 + \Delta^2}$ and $s_2(\omega_n) = \sqrt{(\omega_n+\Omega)^2 + \Delta^2}$. In fact, in the limit $T\rightarrow0$, I get
$$ S(\Omega) = \frac{2\sqrt{4\Delta^2 + \Omega^2}}{\Omega} \tanh^{-1}{\frac{\Omega}{\sqrt{4\Delta^2 + \Omega^2}}}.$$
When I calculate the sum above directly I get the same result. This is nice but I am really interested in the real frequency analytical continuation of $S(\Omega)$. To get it I do the usual change $\Omega \rightarrow i \Omega$. This simply gives me
$$ S_a(\Omega) = \frac{2\sqrt{4\Delta^2 - \Omega^2}}{\Omega} \tan^{-1}{\frac{\Omega}{\sqrt{4\Delta^2 - \Omega^2}}}.$$ This is the correct result. Nonetheless, when I try to calculate numerically the Matsubara sum above with the change $\Omega \rightarrow i \Omega$, I only get the correct result for $\Omega < \Delta$. I suspect that it comes from the fact that the real part inside the square root of $s_2(\omega_n)$ changes sign at $\Omega = \Delta$. One trick I used was to do the change of variable $\omega_n \rightarrow \omega_n - \Omega/2$. In this case the change of sign of the real part inside $s_2(\omega_n)$ changes at $\Omega = 2\Delta$. And indeed the analytical result $S_a$ and the direct sum calculation give the same results for $\Omega <2\Delta$. But the results remain different for $\Omega > 2 \Delta$. My question is: what should I do to get perfect agreement between the sum and the analytical solution $S_a$ for all $\Omega$ ?
To be clear, after the analytical continuation $\Omega \rightarrow i \Omega$, I use the positive real part value for $s_2(\omega_n)$ for all $\Omega$.
Edit: I use a python script to calculate the sum where the square root has branch cut in the negatif axis $] - \infty, 0]$. I can provide my script on request.