that the third-order correlation functions are needed to correctly describe interactions between quasiparticles. In particular, the matrix $M_{k,q}^{(\chi)} = \langle \hat{\chi}_{k-q}^\dagger \hat{\chi}_{q}^\dagger \hat{\chi}_{k} \rangle$ can be used to represent the scattering of a quasiparticle with momentum $k$ into two quasiparticles with momenta $q$ and $k-q$ and vice versa, i.e. Beliaev and Landau scattering. To facilitate our discussion, from now on we go back from the Bogoliubov basis to the original basis of $\hat{\phi}_k$ operators. This requires including two distinct third-order correlators in the dynamics, namely $$\label{eq:3rd_corrs}
M_{k,q} = \langle \hat{\phi}_{k-q}^\dagger \hat{\phi}_{q}^\dagger \hat{\phi}_{k} \rangle,\;\;\; R_{k,q} = \langle \hat{\phi}_{-k-q} \hat{\phi}_{q} \hat{\phi}_{k} \rangle.$$ Making use of equation (\[eq:q\_fluct\]) for the time evolution of quantum fluctuations, one readily derives differential equations for the correlation functions up to third order (see Appendix \[app:HOC\]).

If the lowest order $\psi_0$ in Eq. (\[eq:ansatz\]) is kept fixed to the mean-field value, the inclusion of the third-order correlator leads to a finite value for the first-order correlator $\phi_0=\langle \hat{\phi}_0 \rangle$ as well. Another convenient way of choosing the ansatz (\[eq:ansatz\]) is to set $\phi_0 = 0$ by definition, thus capturing the variation of the condensate wavefunction directly in $\psi_0$. This goes at the cost of adding back-reaction terms to 