2D: numerics for 3rd order term
We naturally expect singularities when $q_0 = \sqrt{3} k_F$ (equilateral triangle with sides $q_0$ inscribed into the Fermi circle), when $q_0 = 2k_F$, and when $q_0 \to 0$. Indeed that is what we find numerically. Interestingly, we do not find any divergences at finite $q_0$ as $T\to 0$. That is unlike the forth order terms in GL where we find divergence when the resonant condition is satisfied.
File 1 contains the MATLAB script used to generate the $q_0$ dependence of 3rd order term in GL. From Figure 1, one can see that even though there are visible singularities at $q_0 = \sqrt{3} k_F$ and $q_0 = 2k_F$, there are no divergences at these values of $q_0$.