the spontaneous occurrence of Beliaev-Landau scattering processes by identifying their footprint on steady-state properties of the non-equilibrium photon gas, such as, most notably, the momentum distribution and the quantum statistics.

Semiconductor microcavities and superconducting circuits are among the most promising platforms for realizing scalable arrays of quantum-optical building blocks, suitable for large-scale quantum simulations (see [@schmidt2013; @Iacopo_QFL; @hartmann2016; @angelakis2016] for recent reviews). Motivated by recent experimental advances, much theoretical effort has been devoted to developing numerical techniques to simulate large-scale driven-dissipative quantum systems. Several reasons cooperate to make these systems computationally much more challenging than the corresponding equilibrium ones. First of all the total photon number is not conserved, resulting in a much larger effective Hilbert space, and secondly the steady state is a mixed state, thus requiring the evaluation of a full density matrix rather than a single wavefunction. This led to the development of new numerical tools such as, among others, variational approaches based on matrix product states in 1D [@verstraete_diss; @orus_infinite; @cui_variational; @mascarenhas_matrix], resummation techniques [@li_2016_resummation], self-consistent projection operator theory [@degenfeld_2014_self], extensions of the variational principle [@weimer_variational] and the corner-space renormalization method for 2D lattices [@finazzi_corner].

In addition to the above-cited exact methods, approximated techniques based on 