  Nonlinear quantum optical systems are of paramount relevance for modern
quantum technologies, as well as for the study of dissipative phase
transitions. Their nonlinear nature makes their theoretical study very
challenging and hence they have always served as great motivation to develop
new techniques for the analysis of open quantum systems. In this article we
apply the recently developed self-consistent projection operator theory to the
degenerate optical parametric oscillator to exemplify its general applicability
to quantum optical systems. We show that this theory provides an efficient
method to calculate the full quantum state of each mode with high degree of
accuracy, even at the critical point. It is equally successful in describing
both the stationary limit and the dynamics, including regions of the parameter
space where the numerical integration of the full problem is significantly less
efficient. We further develop a Gaussian approach consistent with our theory,
which yields sensibly better results than the previous Gaussian methods
developed for this system, most notably standard linearization techniques.
