to diagonalize this propagator, at least for $R_k=0$. To that end we note that the inverse propagator has the position space representation $$\tilde G_{k=0}^{-1}(p) = p_0-\frac{1}{2(M_1+M_2)}\vec p^2-\frac{1}{2\mu}(-\vec \nabla_y^2)-V(\vec y)$$ which is diagonalized, of course, by the solutions of the stationary Schrödinger equation $$\tilde G_{k=0}^{-1}(p) g_n(\vec y)=(p_0-\frac{1}{2(M_1+M_2)}\vec p^2-E_n) g_n(\vec y).$$ Here, $n$ is a combined index that labels all quantum numbers. For a spherical symmetric potential this includes the radial quantum number as well as those for angular momentum. In the basis $$\Phi(p,\vec y) = \sum_n \phi_n(p) g_n(\vec y)$$ we can write $$(\tilde G_{k=0}(p))_{n m} = (p_0-\frac{1}{2(M_1+M_2)}\vec p^2-E_n) \delta_{n m}.$$ For completeness we note the explicit form of the Yukawa coupling term in the effective action in this basis $$\Gamma_k^{\Phi\psi\psi} = \sum_n \int_{p,\vec y} \left\{\phi^*_n(p) g^*_n(\vec y) V(\vec y) (\psi_1\psi_2)(p,\vec y)+c.c.\right\}.$$ For a constant cutoff function $R_k=k^2$ the basis of functions $g_n$ diagonalizes also $\tilde G_k$ for nonzero $k$. However, the situation becomes more complicated for other choices of $R_k$ where the basis would depend on the center of mass momentum $p$ and the flow parameter $k$.

At this point we could in principle undo the Hubbard-Stratonovich transformation and thus go back to our original formulation of the theory in terms 