eventually relax to a thermalized one satisfying $n_k^\text{class} \sim k_BT/\epsilon_k$. Apart from the fact that this Rayleigh-Jeans-like law does not match the expected Bose statistics, the TWA can also lead to negative values for the physical occupation of high-$k$ modes after subtraction of $1/2$ vacuum noise (\[eq:TWA\_dens\]). Therefore one can reliably use the TWA to compute time evolutions only up to a limited time, such that no thermalization sets in for the high momentum modes. Also on the calculation of phonon damping rates this problem has a direct impact, as the unphysical late-time thermalization of the classical field dramatically affects the Landau processes [@castin_TWA].

Hierarchy of correlation functions {#sec:HOC}
==================================

Given the dramatic failure of the TWA classical field approach discussed in the previous section, we need to develop a more sophisticated method, apt to capture the quantum nature of the field operators more accurately. The idea is to go back to the quantum equation of motion (\[eq:q\_fluct\]) and to expand up to higher orders in the quantum fluctuations.

The method
----------

While a linearized form of eq. (\[eq:q\_fluct\]) was sufficient to reproduce the quadratic correlation functions (\[eq:bog\_distr\]), which describe the effect of a nonzero quasiparticle occupation, one can expect 