I will try to give you some guidelines that will (hopefully) help you understand some key differences between Weinberg's work and the authors' in the article you cite (despite not having a complete understanding of the latter). Sorry for not having a definite answer, but I think I do not have the understanding (or time to think) for giving you one. But, I do believe (and hope) that the following will prove to be usefull.
First, Weinberg performs the $k^0$ integration in the soft factor (that is related to the virtual soft divergences) first by closing the contour in 2 ways depending on whether or not the particles $i$ and $j$ belong both in the same state (i.e. incoming/outgoing) or not. In the latter case there is no imaginary part whatsoever (see Weinberg's book on QFT, Chapter 13.2 for the same derivation of the virtual infrared divergences in the case of QED!). Despite this not being the case in the paper you cite, the authors do not perform the $k^0$ integral first, but rather, they use the trick
$$\frac{1}{ABC}=\int_0^{\infty}dx_1\int_0^{\infty}dx_2
\frac{2}{(A+Bx_1+Cx_2)^3}$$
for $A=(k^2+i\varepsilon),\ B=2k\cdot p_i+i\varepsilon$ and $C=-2k\cdot p_j+i\varepsilon$. Then, they complete the square and shift the momentum with respect to which they integrate, such that the first line of Eq. (4.3) is formed. So far so good, with the minor difference that I got $i\varepsilon(1+x_1+x_2)$ instead of their $i\varepsilon$ they have as the last term of the denominator, but since $\varepsilon$ is infinitesimally small, both are equivalent! The sign is what matters. However, they seem to stop keeping track of the $i\varepsilon$ factors upon performing the integral with respect to the shifted momentum. As mentioned earlier, the integral with respect to the shifted momentum is not performed first by integrating the $k^0$ component and then by performing the integration with respect to the spatial components of the momentum.
Instead, what they are doing is a Wick rotation followed by expressing all the components of the (Euclidean version of the) momentum in spherical coordinates. The result is expressed in terms of some Gamma functions (I copy from P&S)
$$\int\frac{d^dk}{(2\pi)^d}\frac{1}{(k^2-\Delta)^3}=
\frac{-i}{(4\pi)^{d/2}}\frac{\Gamma(1+\epsilon)}{\Gamma(3)}
\frac{1}{\Delta^{1+\epsilon}}$$
where $\Delta=(x_1p_i-x_2p_j)^2-i\varepsilon$
if one does the calculations carefully.
Having said that, and keeping in mind that the imaginary part should exist regardless the way we perform the integral, I think that all the information about the imaginary part is encoded in that last term of the denominator $\Delta$! The authors seemingly do not take it into account, as they focus on the real part of the divergences associated to the virtual soft insertions (and I think this is evident from the last line of Eq. (4.3), in which they neglect the imaginary part of the denominator).
This is my guess and hopefully I haven't said anything wrong! You can check yourself the calculations. I hope this helps
P.S. #1: If you study Weinberg's book, there is a sign convention difference between the latter and the paper you study!
P.S. #2: P&S is the Peskin's and Schroeder's book on QFT
