QBS can be written as: $$\begin{aligned}
    S[\omega_1,\omega_2] =
            S_{\rm th}[\omega_1,\omega_2] + 
            S_{\rm dr}[\omega_1,\omega_2],
            \label{eq:QBSDecomposition}\end{aligned}$$ where the first term is completely independent of the drive $f$, and the second term is proportional to $|f|^2$.

[*Quantum bispectrum of drive-independent photon number fluctuations–* ]{} The $f$-independent QBS $S_{\rm th}[\omega_1,\omega_2]$ can be calculated by solving Eq.  with $f=0$. In this case, the cavity relaxes to a thermal steady state with no coherence between different Fock states. Its fluctuations can thus be mapped to a classical Markovian master equation. For such a classical and thermal Markov process, the bispectrum $S_{\mathrm{th} }   [ {\omega}_{1},{\omega}_{2}  ] $ [*must*]{} always be real [@Semerjian2004; @Sinitsyn2016]. In the SM, we also show that in our case, $S_{\mathrm{th} }   [ {\omega}_{1},{\omega}_{2}  ] $ must also be positive semidefinite [@SI]. Letting $\omega_{3} \equiv -\omega_1 -\omega_2$ in all equations that follow, our full calculation for the Keldysh-ordered QBS yields as expected a real, positive function: $$\label{Eq:bispmth}
     S_{\mathrm{th} }   [ {\omega}_{1},{\omega}_{2}  ]  =
        \mathcal{C}_{{\bar n}_\mathrm{th} }  \gamma^2   
        \left( 6 \gamma^2+   \sum\limits_{j=1}^3 { \omega}^2 _j  \right) \Big /  \prod\limits_{j=1}^3 (  \gamma^2  + { \omega}^2 _j   ),$$ with $\mathcal{C}_{{\bar n}_\mathrm{th} }= {\bar n}_\mathrm{th} ( {\bar n}_\mathrm{th} +1  )  ( 2{\bar n}_\mathrm{th} +1) $. We see that 