equations for the quadratic correlation functions $n_k = \langle \hat{\phi}_k^\dagger \hat{\phi}_k \rangle$ and $c_k = \langle \hat{\phi}_k \hat{\phi}_{-k} \rangle$, $$\begin{aligned}
\label{eq:2nd_bog}
\partial_t n_k &=& -\gamma n_k + 2{\text{Im}}\big[U\psi_0^2 c^\ast_k\big]\\
i\partial_t c_k &=& (2\epsilon_k +2U|\psi_0|^2 - i\gamma)c_k + U\psi_0^2(2n_k + 1),\end{aligned}$$ whose steady-state solution reads $$\label{eq:momentum_distr}
n_k = \frac{1}{2}\frac{(Un_0)^2}{\omega_k^2 + \gamma^2/4},\;\;\; c_k = -\frac{U\psi_0^2}{2}\frac{\epsilon_k + Un_0 + i\gamma/2}{\omega_k^2 + \gamma^2/4}.$$

Corrections to the Bogoliubov approximation and Beliaev-Landau processes
------------------------------------------------------------------------

The next-order correction to the Bogoliubov approximation is given by Hamiltonian interaction terms that comprise three fluctuation operators and only one condensate mode, which go under the name of Beliaev and Landau scatterings. Beliaev scattering is the collision of a fluctuation with momentum ${\mathbf{k}}$ with a condensate particle into a pair excitations with momenta ${\mathbf{q}}$ and ${\mathbf{k}}-{\mathbf{q}}$, such that total momentum is conserved (see Fig. \[fig:bel\_decay\] (a)). Landau scattering is the opposite process: two fluctuations with momenta ${\mathbf{k-q}}$ and ${\mathbf{q}}$ scatter into a condensate mode and an excitation with momentum ${\mathbf{k}}$.

In a closed system, all scattering processes must occur on-shell, i.e. they conserve both energy and momentum. For the Beliaev-Landau processes, this implies the following relation: $$\label{eq:E_cons}
\omega_{\mathbf{k}}  = \omega_{\mathbf{q}} +  \omega_{{\mathbf{k}} - {\mathbf{q}}},$$ with $\omega_k$ the quasiparticle oscillation frequency. Because 