  A Gaussian quantum theory of bosonic modes has been widely used to describe
quantum optical systems, including coherent Ising machines (CIMs) that consist
of $\chi^{(2)}$ degenerate optical parametric oscillators (DOPOs) as nonlinear
elements. However, Gaussian models have been thought to be invalid in the
extremely strong-gain-saturation limit. Here, we develop an extended Gaussian
model including two third-order fluctuation products, $\langle \delta
\hat{X}^3\rangle$ and $\langle \delta \hat{X}\delta \hat{P}^2\rangle$, which we
call self-skewness and cross-skewness, respectively. This new model which we
call skew-Gaussian model more precisely replicates the success probability
predicted by the quantum master equation (QME), relative to Gaussian models. We
also discuss the impact of skew variables on the performance of CIMs.
