we set a fixed $\epsilon$ and stop the evolution once $\delta(t) < \epsilon$. Typically $\delta(t) \sim \exp( -\kappa t)$ and therefore convergence is rapidly achieved. For a system with 128 cavities we need about 2 minutes of CPU time on a standard computer, without any optimization, to get an accuracy $\delta < 10^{-6}$.

Results
-------

In Fig. \[fig:2nd\] we present a comparison between the (inaccurate and unreliable) TWA result and the one obtained with the present hierarchy-of-correlations (HOC) approach. First of all, the corrections from the Bogoliubov theory are again as expected proportional to the single-photon interaction constant $U$ for a given mean-field energy $Un_0$. Importantly, the predictions of the HOC do not suffer from unphysical negative occupation numbers and quantitatively the corrections turn out to be much smaller than the ones found in TWA. From a qualitative point of view, we see that they are similar in shape to the TWA ones, but far less pronounced. We therefore conclude that within the TWA, the physical scattering processes between quasiparticles are overwhelmed by the unphysical Beliaev-like decay of the $1/2$ vacuum noise, which is indeed significantly larger than the actual occupation number of the excitations.

Furthermore, we deduce from Fig. 