F^{(R)}_{k,q}.
\label{eq:HE_R}\end{aligned}$$ Here $F^{(M,R)}_{k,q}$ captures the back-reaction of the various separable contributions from the fourth-order correlation functions $$\begin{aligned}
\nonumber
F^{(M)}_{k,q}&=& \frac{2U\psi_0}{\sqrt{L}} \Big(c^\ast_{k-q}n_q + n_{k-q}c^\ast_q - n_k(c^\ast_q+c^\ast_{k-q})\Big)\\\nonumber
&&+ \frac{2U\psi_0^\ast}{\sqrt{L}}\Big(n_{k-q}n_q-n_k(1+n_q +n_{k-q})\\
&& \;\;\;- c_k(c^\ast_q+c^\ast_{k-q})\Big)  \label{eq:F_M}\\\nonumber
F^{(R)}_{k,q}&=& \frac{2U\psi_0}{\sqrt{L}}\Big(c_k + c_q + c_{k+q} + n_{k+q}c_q + c_{k+q}n_q\\\nonumber
&& \;\;\; + n_k c_q + n_k c_{k+q} + c_k n_q + c_k n_{k+q} \Big)\\
&&+\frac{2U\psi_0^\ast}{\sqrt{L}}\Big(c_k c_q + c_k c_{k+q} + c_q c_{k+q} \Big)
\label{eq:F_R}\end{aligned}$$ In principle equations (\[eq:HE\_MF\]-\[eq:HE\_R\]) provide a solution to the full time-dependent problem when appropriate initial conditions are inserted. The focus of the present work is, however, on the steady-state solution. To obtain this, we in practice initialize the system with the mean-field condensate amplitude $\psi_0$ and the Bogoliubov solution (\[eq:momentum\_distr\]) for $n_k$ and $c_k$, and we initially set $M_{k,q}$ and $R_{k,q}$ to zero. We then let the system evolve until it spontaneously reaches its steady-state.

To follow the time-evolution, we have implemented a Runge-Kutta-based routine with adaptive timestep to integrate equations (\[eq:HE\_MF\]-\[eq:HE\_R\]) in time. By plotting a quantity such as $ \delta(t) = 1/(L\Delta t)\sum_k{ \big| n_k^{t + \Delta t} - n_k^t\big| / n_k^t} $ as a function of $t$ for fixed $\Delta t$, we can monitor the convergence. As a criterion 