Self-interaction

There are two types of self-interaction terms:

$$\Delta \Omega^{(4)}_A = \frac{v_{q_0}^4 \rho_{q} \rho_{q}\rho_{-q}\rho_{-q}}{4 \beta} \sum_{\omega_n, k} \frac{1}{(i\omega_n - \epsilon_{k})(i\omega_n - \epsilon_{k+q})^2 (i\omega_n - \epsilon_{k+2q})},$$

$$\Delta \Omega^{(4)}_B = \frac{v_{q_0}^4 \rho_{q} \rho_{-q}\rho_{q}\rho_{-q}}{4 \beta} \sum_{\omega_n, k} \frac{1}{(i\omega_n - \epsilon_k)^2 (i\omega_n - \epsilon_{k+q})^2},$$

They can be converted into contour integrals by introducing the Fermi function,

$$\Delta \Omega^{(4)}_A \sim \sum_k \int_C dz \frac{n_F(z)}{(z - \epsilon_{k})(z - \epsilon_{k+q})^2 (z - \epsilon_{k+2q})}\\\ = \sum_k \frac{n_F(\epsilon_{k})}{(\epsilon_{k} - \epsilon_{k+q})^2 (\epsilon_{k} - \epsilon_{k+2q})} + \frac{n_F(\epsilon_{k+2q})}{(\epsilon_{k+2q} - \epsilon_{k})(\epsilon_{k+2q} - \epsilon_{k+q})^2}\\\ - \frac{n_F(\epsilon_{k+q})}{(\epsilon_{k+q} - \epsilon_{k})^2 (\epsilon_{k+q} - \epsilon_{k+2q})} - \frac{n_F(\epsilon_{k+q})}{(\epsilon_{k+q} - \epsilon_{k}) (\epsilon_{k+q} - \epsilon_{k+2q})^2} + \frac{n'_F(\epsilon_{k+q})}{(\epsilon_{k+q} - \epsilon_{k}) (\epsilon_{k+q} - \epsilon_{k+2q})} ,$$
and
$$\Delta \Omega^{(4)}_B \sim \sum_k \int_C dz \frac{n_F(z)}{(z - \epsilon_{k})^2(z - \epsilon_{k+q})^2 }\\\ = \sum_k -\frac{n_F(\epsilon_{k})}{(\epsilon_{k} - \epsilon_{k+q})^3} - \frac{n_F(\epsilon_{k+q})}{(\epsilon_{k+q} - \epsilon_{k})^3}+ \frac{n'_F(\epsilon_{k})+n'_F(\epsilon_{k+q})}{(\epsilon_{k} - \epsilon_{k+q})^2} .$$

The origin of different terms can be understood from a simple 2-level system analysis (see attached note).

It is interesting to note, that at finite temperature, even though the denominators are singular, the functions under the sums are in fact finite (singularities of the same type as $[n_F(E_1) - n_F(E_2)]/(E_1 - E_2)\approx n_F'(E_1)$ at $E_2 \to E_1$ that occurs in susceptibility). Even though using symmetries of the electron dispersion it is possible to rewrite these more compactly, for the purposes of numerics the function should be kept as is, to be able to integrate them with minimal amount trouble. The form of the functions can be played with using the attached Mathematica script.

The summations over momenta, which become integrals in the thermodynamic limit are singular. $\Delta \Omega^{(4)}_B$ can be computed analytically in the limit of $T=0$ for $q> 2k_F$,

$$ \Delta \Omega^{(4)}_B = \frac{v^4 \| \rho_q \|^4}{2} (B_1 + B_2) \\\ = \frac{v^4 \| \rho_q \|^4}{2} \left( \frac{4}{\pi^{} q^3} \frac{k_{F^{}}^2}{[q^2 - 4 k_F^2]^{3 / 2}} - \frac{2}{\pi q^{}} \frac{1}{[q^2 - 4 k_F^2]^{3 / 2}} \right) \\\ = \frac{v^4 \| \rho_q \|^4}{4 (2 \pi)^2} \times \frac{16 \pi}{q^3} \frac{2 k_{F^{}}^2 - q^2}{[q^2 - 4 k_F^2]^{3 / 2}} $$
The details are in the attached TeXmacs note.

Numerically, momentum integrals have to be regularized. One possible way is
$\frac{1}{(E_1 - E_2)^n} \to Re \frac{1}{(E_1 - E_2 + i\Gamma)^n}$. $\Gamma$ should be kept ridiculously small to allow for the natural regularization by temperature. Then, the analytical and numerical results indeed agree. With this regularization, and "good" choice of integrands, the mutual interaction matches the self-interaction terms in the appropriate limit.