We implement an extension of the pseudofermion functional renormalization group method for quantum spin systems that takes into account two-loop diagrammatic contributions. An efficient numerical treatment of the additional terms is achieved within a nested graph construction which recombines different one-loop interaction channels. In order to be fully self-consistent with respect to self-energy corrections, we also include certain three-loop terms of Katanin type. We first apply this formalism to the antiferromagnetic J1−J2 Heisenberg model on the square lattice and benchmark our results against the previous one-loop plus Katanin approach. Even though the renormalization group (RG) equations undergo significant modifications when including the two-loop terms, the magnetic phase diagram, comprising Néel ordered and collinear ordered phases separated by a magnetically disordered regime, remains remarkably unchanged. Only the boundary position between the disordered and the collinear phases is found to be moderately affected by two-loop terms. On the other hand, critical RG scales, which we associate with critical temperatures Tc, are reduced by a factor of ∼2 indicating that the two-loop diagrams play a significant role in enforcing the Mermin-Wagner theorem. Improved estimates for critical temperatures are also obtained for the Heisenberg ferromagnet on the three-dimensional simple cubic lattice where errors in Tc are reduced by ∼34%. These findings have important implications for the quantum phase diagrams calculated within the previous one-loop plus Katanin approach which turn out to be already well converged.
