Spontaneous Beliaev-Landau scattering out of equilibrium

TL;DR

This paper studies Beliaev-Landau scattering in a driven-dissipative photon gas, revealing non-equilibrium effects and non-Gaussian features using an advanced theoretical approach beyond standard approximations.

Contribution

It introduces a systematic expansion beyond the Bogoliubov approximation to accurately capture Beliaev-Landau processes in out-of-equilibrium photon systems.

Findings

Detected characteristic peaks and dips in the momentum distribution due to scattering processes.

Identified non-Gaussian features in cavity output field statistics.

Showed limitations of the truncated Wigner approximation for this physics.

Abstract

We investigate Beliaev-Landau scattering in a gas of interacting photons in a coherently driven array of nonlinear dissipative resonators, as described by the 1D driven-dissipative Bose-Hubbard model. Due to the absence of detailed balance in such an out-of-equilibrium setup, steady-state properties can be much more sensitive to the underlying microscopic dynamics. As the popular truncated Wigner approximation dramatically fails in capturing this physics, we present an alternative approach, based on a systematic expansion beyond the Bogoliubov approximation, which includes the third-order correlation functions in the dynamics. As experimentally accessible signatures of Beliaev-Landau processes, we report a small but nonnegligible correction to the Bogoliubov prediction for the steady-state momentum distribution showing a characteristic series of peaks and dips, as well as non-Gaussian features in the statistics of the cavity output field.

Equations

$$\hat{H}$$

$$\mathcal{D}[\hat{\rho }]=\frac{\gamma }{2}\sum _j(2{\hat{a}}_j\hat{\rho }{\hat{a}}_j^{†}-\hat{\rho }{\hat{n}}_j-{\hat{n}}_j\hat{\rho }).$$

$${\partial }_t\hat{\rho }=-i[\hat{H},\hat{\rho }]+\mathcal{D}[\hat{\rho }]$$

$$i\dot{\psi}_{j}=-\Big{(}\delta+i\frac{\gamma}{2}\Big{)}\psi_{j}-J(\psi_{j+1}+\psi_{j-1})+U|\psi_{j}|^{2}\psi_{j}+\Omega_{j}$$

$$n_0((\delta -U{n}_0+2J{)}^2+{\gamma }^2/4)=|\Omega {|}^2,$$

$$\Delta =\delta -U{n}_0+2J$$

$${\hat{a}}_j={\psi }_0+\frac{1}{\sqrt{L}}\sum _k{e}^{ikj}{\hat{\varphi }}_k$$

$$i{\partial }_t{\hat{\varphi }}_k$$

$${ϵ}_k=-\delta +U{n}_0-2J\cos k.$$

$$\displaystyle\big{\langle}\hat{\xi}_{k}(t)\hat{\xi}_{k^{\prime}}(t^{\prime})\big{\rangle}$$

$$\displaystyle\big{\langle}\hat{\xi}_{k}(t)\hat{\xi}_{k^{\prime}}^{\dagger}(t^{\prime})\big{\rangle}$$

$${\hat{\varphi }}_k=u_k{\hat{\chi }}_k+v_k{\hat{\chi }}_{-k}^{†}$$

$$i({\partial }_t+\gamma /2){\hat{\chi }}_k={\omega }_k{\hat{\chi }}_k+u_k{\hat{\xi }}_k-v_k{\hat{\xi }}_{-k}^{†}.$$

$${\omega }_k$$

$$u_k,v_k$$

$$n_k^{(\chi )}=⟨{\hat{\chi }}_k^{†}{\hat{\chi }}_k⟩=v_k^2,c_k^{(\chi )}=⟨{\hat{\chi }}_k{\hat{\chi }}_k⟩=\frac{u_k{v}_k\gamma }{\gamma +2i{\omega }_k}$$

$${\partial }_t{n}_k$$

$$i{\partial }_t{c}_k$$

$$n_k=\frac{1}{2}\frac{(U{n}_0{)}^2}{{\omega }_k^2+{\gamma }^2/4},c_k=-\frac{U{\psi }_0^2}{2}\frac{{ϵ}_k+U{n}_0+i\gamma /2}{{\omega }_k^2+{\gamma }^2/4}.$$

$${\omega }_{\mathbf{k}}={\omega }_{\mathbf{q}}+{\omega }_{\mathbf{k}-\mathbf{q}},$$

$$\begin{split}id\varphi_{j}(t)=\Big{[}-\Big{(}\delta+i\frac{\gamma}{2}\Big{)}\varphi_{j}(t)-J\big{(}\varphi_{j+1}(t)+\varphi_{j-1}(t)\big{)}\\ +U\big{(}|\varphi_{j}(t)|^{2}-1\big{)}\varphi_{j}(t)+\Omega_{j}(t)\Big{]}dt+\sqrt{\frac{\gamma}{2}}dW_{j}(t),\end{split}$$

$$\big{\langle}\varphi_{j}^{\ast}\varphi_{j}\big{\rangle}_{W}=\frac{1}{2}\Big{(}\langle\hat{a}_{j}^{\dagger}\hat{a}_{j}\rangle+\langle\hat{a}_{j}\hat{a}_{j}^{\dagger}\rangle\Big{)}=\langle\hat{n}_{j}\rangle+\frac{1}{2}.$$

$$M_{k,q}=⟨{\hat{\varphi }}_{k-q}^{†}{\hat{\varphi }}_q^{†}{\hat{\varphi }}_k⟩,R_{k,q}=⟨{\hat{\varphi }}_{-k-q}{\hat{\varphi }}_q{\hat{\varphi }}_k⟩.$$

$$i{\partial }_t{\psi }_0$$

$$i{\partial }_t{n}_k$$

$$i{\partial }_t{c}_k$$

$$i{\partial }_t{M}_{k,q}$$

$$i{\partial }_t{R}_{k,q}$$

$$F_{k,q}^{(M)}$$

$$F_{k,q}^{(R)}$$

$$\begin{gathered} ⟨{\hat{\chi }}_{k-q}^{†}{\hat{\chi }}_q^{†}{\hat{\chi }}_k⟩= \\ \frac{2U}{\sqrt{L}}\frac{{\psi }_0{A}_{k,q}+{\psi }_0^{\ast }{B}_{k,q}}{{\omega }_k-{\omega }_q-{\omega }_{k-q}-\frac{3i}{2}\gamma }, \end{gathered}$$

$$\begin{split}\Big{\langle}\big{(}\hat{\phi}_{q}^{\dagger}+e^{-i\theta}\hat{\phi}_{k}^{\dagger}\big{)}\big{(}\hat{\phi}_{q}+e^{i\theta}\hat{\phi}_{k}\big{)}\big{(}e^{i\chi}\Omega\hat{\phi}_{k-q}^{\dagger}+e^{-i\chi}\Omega^{*}\hat{\phi}_{k-q}\big{)}\Big{\rangle}\\ =2\Re\Big{\{}\Omega e^{i(\theta+\chi)}\langle\hat{\phi}_{k-q}^{\dagger}\hat{\phi}_{q}^{\dagger}\hat{\phi}_{k}\rangle\Big{\}}=2\Re\Big{\{}\Omega e^{i(\theta+\chi)}M_{k,q}\Big{\}},\end{split}$$

$$i\frac{\partial C}{\partial t}$$

Full text

I Introduction

The convenient assumption of detailed balance, valid for a quantum many-body system at thermal equilibrium, ensures that any microscopic process is balanced by its reverse process, thus making the specific underlying dynamics irrelevant for the equilibrium ensemble. For systems far from equilibrium no such claims can be made. Although some similarities can exist Sieberer et al. (2016); Foss-Feig et al. (2016); Lebreuilly et al. (2016), driving and dissipation in general prevent the system from approaching a complete thermal equilibrium Carmichael (2009). As a result, the phase-space distribution of out-of-equilibrium systems is typically much more sensitive to the actual microscopic driving, dissipation and equilibration processes that drive the system towards the steady-state.

In this manuscript we study a weakly interacting driven-dissipative Bose-Hubbard model in the superfluid regime. Going beyond the widely applied Bogoliubov approximation, which assumes noninteracting quasiparticles, we show that scattering processes involving three quasiparticles, known as Beliaev-Landau scattering Beliaev (1958), leave a small imprint on the steady-state momentum distribution and are responsible for non-Gaussian features in the photon statistics of the cavity output field.

Beliaev-Landau scattering processes were originally predicted in the many-body theory of quantum fluids Pines and Nozières (1998) to be responsible for the finite lifetime of phonons in systems of either bosonic Pitaevskii and Stringari (1997); Giorgini (1998); Kónya et al. (2014) or fermionic particles Kurkjian et al. (2016) and were experimentally observed with ultracold atomic gasesHodby et al. (2001); Katz et al. (2002). In these works, the damping was typically detected after a population of phonons was introduced in the system by an external perturbation. Here, on the contrary, we aim at capturing 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 Schmidt and Koch (2013); Carusotto and Ciuti (2013); Hartmann (2016); Noh and Angelakis (2016) 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 et al. (2004); Orus and Vidal (2008); Cui et al. (2015); Mascarenhas et al. (2015), resummation techniques Li et al. (2016), self-consistent projection operator theory Degenfeld-Schonburg and Hartmann (2014), extensions of the variational principle Weimer (2015) and the corner-space renormalization method for 2D lattices Finazzi et al. (2015).

In addition to the above-cited exact methods, approximated techniques based on the truncated Wigner approximation are also very popular tools to evaluate corrections beyond the Bogoliubov approximation in both conservative cold-atom Sinatra et al. (2002); Polkovnikov (2010) and lossy optical systems Gardiner and Haken (1991); Carusotto and Ciuti (2013). For the latter, this technique has been applied in various contexts, including the study of condensation and superfluid properties Carusotto and Ciuti (2005); Wouters and Savona (2009); Gladilin et al. (2014); He et al. (2015), dynamical phase transitions Sieberer et al. (2013); Dagvadorj et al. (2015); Foss-Feig et al. (2016), and even genuine quantum effects such as the dynamical Casimir emission Koghee and Wouters (2014) and Hawking radiation Gerace and Carusotto (2012); Grišins et al. (2016). In contrast to these successes, we will show in this work that the truncated Wigner approximation, when naively adopted to study Beliaev-Landau scattering, may dramatically overestimate the corrections to Bogoliubov theory and even lead to unphysical results.

As an alternative approach we discuss how a truncated hierarchy of correlations can serve as a consistent expansion beyond the Bogoliubov approximation del Valle and Hartmann (2013); Pucci et al. (2016); Casteels et al. (2016). In particular, we will illustrate that the truncation of the hierarchy at the third-order correlation functions, i.e. one order beyond the Bogoliubov approximation, is sufficient to incorporate the corrections attributed to Beliaev-Landau processes, provided an adequate truncation scheme is employed.

The structure of our paper is as follows. In Sec. II we present the 1D driven-dissipative Bose-Hubbard model and in Sec. III we illustrate how the non-equilibrium condition allows for on-shell Beliaev-Landau scattering in one-dimension. We next explain how the truncated Wigner approximation fails to describe these processes in Sec. IV. In Sec. V we introduce the third-order correlation functions to incorporate Beliaev and Landau scattering and construct a hierarchy of correlation functions. In Sec. VI we discuss the imprint of Beliaev-Landau scattering on measurable quantities and discuss the expected signal for realistic parameters inspired from state-of-the-art semiconductor devices. Conclusions are finally drawn in Sec. VII. Appendix A reports additional numerical TWA data in the absence of open Beliaev-Landau channels. Appendices B and C summarize technical details on the hierarchy of correlations and on the different truncation schemes. Details on the calculation of the effect of disorder are given in Appendix D.

II The Model

We consider a 1D coupled array of $L$ nonlinear, single-mode photon cavities under a coherent drive with frequency $\omega_{L}$. The resonator frequencies $\omega_{c}$ are assumed to be uniform throughout the chain. After a unitary transformation to remove the time-dependence of the drive, we obtain the driven Bose-Hubbard Hamiltonian (we set $\hbar=1$ throughout the article)

[TABLE]

The operators $\hat{a}^{\dagger}_{j}(\hat{a}_{j})$ create (annihilate) a particle at site $j$ of the chain and $\hat{n}_{j}=\hat{a}^{\dagger}_{j}\hat{a}_{j}$ is the local number operator. Photons in the chain can tunnel to their neighbouring sites with a hopping strength given by $J$. The notation $\langle j,l\rangle$ means that the summation runs over all neighbouring sites. The two-body interaction strength for photons confined inside the same cavity is given by $U$. The amplitude of the driving field at each site is $\Omega_{j}$, while its detuning from the onsite single-photon resonance is given by $\delta=\omega_{L}-\omega_{c}$. For simplicity we impose periodic boundary conditions, such that $\hat{a}_{1}=\hat{a}_{L+1}$.

The dissipative nature of the setup implies that injected photons have a finite lifetime inside the cavity array before they escape. In the Born-Markov approximation, the coupling of the system to its environment at zero temperature is described by the dissipator in the Lindblad form Gardiner and Haken (1991); Breuer and Petruccione (2002)

[TABLE]

The full dynamics of the density matrix $\hat{\rho}$ is then governed by a master equation, which includes both the unitary evolution under $\hat{H}$ and the photonic losses

[TABLE]

A mean-field description of the problem can be derived in terms of coherent fields $\psi_{j}=\langle\hat{a}_{j}\rangle$ by assuming that all normal-ordered operator products factorize Carusotto and Ciuti (2013). This leads to the following equations of motion

[TABLE]

For this work we assume a uniform drive field $\Omega_{j}=\Omega$ at all sites in the chain. In the steady state we thus find one or two stable homogeneous density solutions in the mean-field description, depending on the amplitude and detuning of the pump. They are found as solutions of

[TABLE]

where $n_{0}=|\psi_{0}|^{2}$ in terms of the spatially uniform mean-field steady state $\psi_{j}=\psi_{0}$. The parameter

[TABLE]

is the renormalized laser detuning from the interaction-blueshifted optical resonance. We restrict our analysis to the case $\Delta<0$, such that the system is in the optical limiter regime or in the high-density branch of the hysteresis loop of a bistable regime Carusotto and Ciuti (2013). As we will briefly review later, this restriction asserts a gapped spectrum of excitations Carusotto and Ciuti (2004). In the numerical analysis that follows, the drive amplitude $\Omega$ is always implicitly determined by choosing a value for $n_{0}$.

III Bogoliubov theory and Beliaev-Landau processes

III.1 Bogoliubov dispersion and the non-equilibrium steady-state

The case of a uniform drive field allows for a convenient parametrization of the full quantum field in terms of a homogeneous, coherent field $\psi_{0}$ and quantum fluctuations. Expanding the latter in their $\hat{\phi}_{k}$ spatial Fourier components, one can write

[TABLE]

where the sum over $k$ is restricted to the interval $[-\pi,\pi]$ with a spacing equal to $2\pi/L$ in the case of a finite number $L$ cavities and periodic boundary conditions.

While the time evolution of the mean-field $\psi_{0}$ is governed by the classical equation (4), the dynamics of the quantum fluctuations is governed by a quantum Langevin equation Gardiner and Haken (1991); Carusotto and Ciuti (2013)

[TABLE]

where we have set

[TABLE]

The Markovian losses are responsible for quantum noise with Gaussian statistics, represented by the operators $\hat{\xi}_{k}$, that assume the following zero-temperature statistics

[TABLE]

As usual in Bogoliubov-like approaches, the interaction terms in (8) are ordered in increasing number of fluctuation operators. When the number of photons in the condensate $|\psi_{0}|^{2}$ is much larger than the number of fluctuations, one expects the effect of higher-order terms to be negligible Castin and Dum (1998).

The first-order correction to the mean-field, summarized on the first line of (8), incorporates processes where two condensate particles collide and produce a pair of excitations with counter-propagating wavevectors $k$ and $-k$ and viceversa. In a Hamiltonian formalism, this corresponds to only retaining quadratic terms in the fluctuation operators $\hat{\phi}_{k}$. Restricting to these terms in (8) and dropping the ones on the second and third line, which contain terms with more than one fluctuation operator, results in a set of linear equations for the fluctuation fields $\hat{\phi}_{k}$. As done in Busch et al. (2014), this set of linear stochastic equations is solved by means of a Bogoliubov transform to new operators

[TABLE]

that diagonalize the equations of motion

[TABLE]

Here we have defined the quasiparticle energies $\omega_{k}$ along with the transformation functions $u_{k}$, $v_{k}$ as

[TABLE]

In Fig. 1(a) the Bogoliubov spectrum (14) is shown for different values of the renormalized detuning $\Delta$ (6). In contrast to equilibrium systems, note that a spectral gap is generally present in the Bogoliubov dispersion and only closes for $\Delta\rightarrow 0^{(-)}$, i.e. when the drive is exactly on resonance with the interaction-shifted mode Carusotto and Ciuti (2013).

Due to the noise operators in (13), one has a finite occupation of Bogoliubov modes with non-trivial anomalous correlations in the stationary regime,

[TABLE]

It is important to note that the occupation of the Bogoliubov modes here, in contrast with an equilibrium system, is not at all set by a finite-temperature Boltzmann-Gibbs distribution but by the interplay of interactions, hopping, driving and dissipation Busch et al. (2014). This different origin is apparent in the slow, power-law decay of the occupation of high momentum modes, much slower than the usual exponential \exp{\big{[}-E(k)/k_{B}T\big{]}} of equilibrium systems.

Moving back to the original $\hat{\phi}_{k}$ operators, one can derive a closed system of linear differential 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$,

[TABLE]

whose steady-state solution reads

[TABLE]

III.2 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. 1 (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:

[TABLE]

with $\omega_{k}$ the quasiparticle oscillation frequency. Because of the absence of a spectral gap and the convexity of the Bogoliubov dispersion of conservative continuum systems, criterium (20) in continuum models can only be satisfied in two or more spatial dimensions, while in 1D sytems phonons can only decay through higher-order scattering processes Ristivojevic and Matveev (2016). However, there exist specifically engineered 1D optical lattices with a nonconvex (but gapless) spectrum, such that energy and momentum conservation can be simultaneously satisfied Tsuchiya and Griffin (2004); Arahata and Nikuni (2008).

The situation is different in driven-dissipative systems, where the Bogoliubov spectrum is typically gapped when the drive is below resonance. In a continuous 1D setup, the presence of a finite spectral gap in combination with a convex excitation spectrum always allows for third-order scattering process that satisfy (20). When the spectrum is not convex, as is the case in a lattice model, the situation is somehow more complicated. The subtle interplay between the spectral gap and the degree of nonconvexity determines whether resonant third-order scattering channels are present.

The allowed wavevectors $k$ and $q$ that exactly satisfy the energy and momentum conservation condition (20) in one-dimension are indicated in Fig. 1b). Importantly, maximal and minimal values can be deduced for $k$ and $q$ from the contours, which set limits on allowed in and out states for Beliaev-Landau scattering. Only excitations with a wavevector $k$ for which $k_{\text{min}}<k<k_{\text{max}}$ can scatter resonantly to excitations with wavevectors $q$ and $k-q$ through Beliaev decay. Likewise, only excitations with a wavevector $q$ for which $q_{\text{min}}<q<q_{\text{max}}$ can combine with an excitation at $k-q$ to form one at $k$ through Landau scattering. However, the driven-dissipative nature of our setup allows energy not to be strictly conserved, so that scattering processes are possible within a finite linewidth $\gamma$ around the energy-conservation point.

In Fig. 1c) we show how the extremal input and output momenta shift as a function of the interaction-renormalized detuning, defined in $\Delta$ (6), a parameter that can be tuned in experiment by changing the laser frequency $\omega_{L}$. When the drive is too far below resonance, i.e. when $\Delta<\Delta_{0}$, with $\Delta_{0}<0$ a critical value that can be derived from the dispersion relation, the spectral gap is too large as compared to the bandwidth and no resonant Beliaev-Landau scattering channels exist. In the limit of $\Delta\rightarrow 0$, for which the dispersion relation is linear, we find the contour of an equilibrium condensate from Ref. Tsuchiya and Griffin (2004). We anticipate at this point that the experimental possibility of shifting the limiting scattering momenta in a well-controlled manner provides a genuine signature of Beliaev-Landau scattering.

Before continuing with our analysis, we would like to draw attention to an important consideration. Given a closed quantum system (e.g. a gas of ultracold atoms or a superfluid liquid Helium sample), the presence of detailed balance will unavoidably restrict the effect of Beliaev-Landau scattering to driving the system back into its thermal state once it is kicked out of equilibrium. Therefore most works on this physics are related to phonon-decay experiments, where one studies how externally injected phonons are damped through scattering with the condensate (Beliaev) or with the thermal cloud (Landau) Pitaevskii and Stringari (1997); Giorgini (1998) or how a thermal equilibrium is reached again after a sudden global quench Menegoz and Silva (2015).

Consequently, it is exactly the absence of detailed balance in a driven-dissipative context which motivates us to study the effects of spontaneous Beliaev-Landau processes in the steady-state regime of the cavity array. In this section we have illustrated that two crucial conditions for these scatterings to be possibly relevant are indeed satisfied: 1) there is finite occupation of Bogoliubov modes over the entire Brillouin zone, as given in (16), and 2) there are regions in phase space for which energy and momentum are conserved according to (20), which allows the Bogoliubov modes to scatter and redistribute (quasi)resonantly.

IV The Truncated Wigner Method

A possible approach to compute corrections beyond Bogoliubov and quantify the observable signatures of the Beliaev-Landau processes is the so-called Truncated Wigner approximation (TWA) Gardiner and Haken (1991); Vogel and Risken (1989); Sinatra et al. (2002); Polkovnikov (2010). This approach is based on a one-to-one mapping of the quantum master equation for the density matrix (3) onto a partial differential equation for the corresponding Wigner distribution. The resulting Fokker-Planck equation can be efficiently simulated if the terms with a third-order derivative are neglected. Since these terms are proportional to the single-particle interaction constant $U$, one expects this approximation to be accurate for sufficiently weak values of $U$ Carusotto and Ciuti (2013). This leads to a stochastic differential equation for a classical field $\varphi_{j}(t)$

[TABLE]

where the stochastic Wiener increment $dW_{j}(t)$ is white Gaussian noise with variance \big{\langle}dW^{\ast}_{j}(t)dW_{j^{\prime}}(t)\big{\rangle}=\delta_{j,j^{\prime}}dt and a random phase. Average values of the field $\varphi_{j}$ correspond to expectation values of symmetrically ordered products of quantum operators. In particular, for the number operator we find

[TABLE]

As a consequence, the quantum vacuum is represented by a finite occupation of $1/2$ for the classical field $\varphi_{j}$. As long as $n_{j}\gg 1/2$, one does not expect this to cause problems but, when performing a TWA simulation to estimate the effects of Beliaev-Landau scattering, one finds surprisingly large corrections to the occupation numbers of quantum fluctuations, as can be seen in Fig. 2a). Even worse is that the occupation numbers may become negative at certain values of the momentum, even though the used parameters are well inside the supposed region of validity of TWA. In Fig. 2 we fixed $Un_{0}=10\gamma$ and show the results for two interaction constants $U=0.02\gamma$ and $U=0.1\gamma$, such that the mean-field predictions for the number of particles per site are $n_{0}=500$ and $n_{0}=100$, respectively. The results were obtained by averaging out over a total number of about $10^{6}$ samples, which were collected by integrating (21) in time with small enough time step $\Delta t$, and then taking a statistically independent sample each $\tau_{s}=5\gamma^{-1}$.

As expected, the magnitude of the correction to the Bogoliubov theory is proportional to the single-photon interaction constant $U$ (or, equivalently, to the inverse of $n_{0}$ at a given mean-field energy $Un_{0}$). For both values of $U$, the negative occupation of some high-$k$ modes is a clearly unphysical prediction of the TWA.

To better understand the physical origin of this breakdown, one needs to take a closer look at the nature of the underlying physical processes. Through Beliaev scattering, a quasiparticle at a high momentum $k$ decays into two quasiparticles with smaller momenta $q$ and $k-q$. Since the occupation decreases for larger $k$-modes, one expects the importance of this effect to be suppressed at higher momenta. However, within the TWA the quantum field is represented as a classical field for which the occupation of high-$k$ modes does not decay to [math] but to $1/2$, which represents the quantum vacuum fluctuations (see Eq. (22)). This finite occupation of all modes, even the highest-$k$ ones, results in the possibility of a nonphysical decay of the quantum vacuum through spontaneous Beliaev processes. The final states of these collisions are quasiparticles with smaller momenta, which explains the massive pileup in the momentum distribution around $q_{\text{min}}$, at the cost of a strong negative dip around $k_{\text{max}}$. Of course, the TWA-simulated momentum distribution recovers relatively well to the Bogoliubov result for all $k$ values outside of the region $[q_{\text{min}},k_{\text{max}}]$ for which there are no resonant Beliaev-Landau scatterings possible. In Appendix A we include a simulation of a model without energy-conserving Beliaev-Landau channels and we conclude that in this case the occupation of all modes is positive and much better convergence to the Bogoliubov result is achieved. While this inaccuracy of the TWA is not expected to affect the predictions for dynamical Casimir and Hawking emission Koghee and Wouters (2014); Gerace and Carusotto (2012); Grišins et al. (2016) that are at the level of Bogoliubov theory, special care will be needed in the more advanced study of back-reaction effects in analog models of gravity Balbinot et al. (2005).

As far as we we know, there is no simple solution to this intrinsic problem of the TWA. Note that related problems with the TWA are known also in the conservative case of ultracold atomic Bose gases Sinatra et al. (2002). The equipartition theorem for the classical fields states in fact that the momentum distribution should eventually relax to a thermalized one satisfying $n_{k}^{\text{class}}\sim k_{B}T/\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 (22). 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 Sinatra et al. (2002).

V Hierarchy of correlation functions

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 (8) and to expand up to higher orders in the quantum fluctuations.

V.1 The method

While a linearized form of eq. (8) was sufficient to reproduce the quadratic correlation functions (16), which describe the effect of a nonzero quasiparticle occupation, one can expect that the third-order correlation functions are needed to correctly describe interactions between quasiparticles. In particular, the matrix $M_{k,q}^{(\chi)}=\langle\hat{\chi}_{k-q}^{\dagger}\hat{\chi}_{q}^{\dagger}\hat{\chi}_{k}\rangle$ can be used to represent the scattering of a quasiparticle with momentum $k$ into two quasiparticles with momenta $q$ and $k-q$ and vice versa, i.e. Beliaev and Landau scattering. To facilitate our discussion, from now on we go back from the Bogoliubov basis to the original basis of $\hat{\phi}_{k}$ operators. This requires including two distinct third-order correlators in the dynamics, namely

[TABLE]

Making use of equation (8) for the time evolution of quantum fluctuations, one readily derives differential equations for the correlation functions up to third order (see Appendix B).

If the lowest order $\psi_{0}$ in Eq. (7) is kept fixed to the mean-field value, the inclusion of the third-order correlator leads to a finite value for the first-order correlator $\phi_{0}=\langle\hat{\phi}_{0}\rangle$ as well. Another convenient way of choosing the ansatz (7) is to set $\phi_{0}=0$ by definition, thus capturing the variation of the condensate wavefunction directly in $\psi_{0}$. This goes at the cost of adding back-reaction terms to the generalized Gross-Pitaevskii equation. At the level of approximation considered in this section, both approaches are equivalent, but setting $\phi_{0}=0$ asserts that we are dealing with connected second- and third-order correlation functions, which is a better controlled truncation Köhler and Burnett (2002). We refer the interested readers to Appendix B for more details on the method and to Appendix C for a comparison of different truncation schemes.

Correlation functions of order four, which enter into the equations of motion of the second and third-order correlators, are factorized into different possible products of second-order correlation functions. With this procedure we explicitly neglect the connected part of the fourth-order correlation function, but we keep its main contribution coming from separable correlations. The fifth-order correlator, entering in the equation of motion for the third-order correlation functions, can be instead safely neglected. Already in factorized form it would reduce to various products of second and third order, which constitute negligible corrections to dominant terms in the equations of motion. See Appendix B for more details on the implications of these approximations.

Within this framework, the Gross-Pitaevskii equation for the homogeneous condensate background $\psi_{0}$, extended with the back-reaction terms reads

[TABLE]

The second-order correlation functions (17) are now coupled to the third-order correlation functions (23)

[TABLE]

Note that the factorized contribution of the fourth-order correlator enters here in the equations of motion as a small correction to the couplings $\psi_{0}^{2}$ and $|\psi_{0}|^{2}$. At equilibrium these corrections are well-studied in the Hartree-Fock-Bogoliubov method Andersen (2004).

Finally we also find the equations of motion for the third-order correlation functions, in which the fifth-order back-reaction is neglected

[TABLE]

Here $F^{(M,R)}_{k,q}$ captures the back-reaction of the various separable contributions from the fourth-order correlation functions

[TABLE]

In principle equations (V.1-28) 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 (19) 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 (V.1-28) 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 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}$.

V.2 Results

In Fig. 2 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. 2b) that the overall redistribution of particles is from high $k$ to small $k$, which means that the Beliaev decay of high-momentum quasiparticles is dominant. Around the two extremal momenta of the input states, $q_{\text{min}}$ and $q_{\text{max}}$, this manifests itself as a peak in the momentum distribution, while there is a dip at the extremal momenta $k_{\text{min}}$ and $k_{\text{max}}$ of the output states. This is a consequence of the relatively large density of states for possible output (input) states for scatterings with input (output) momenta around $k_{\text{min}}$ or $k_{\text{max}}$ ($q_{\text{min}}$ or $q_{\text{max}}$), as one can deduce in Fig.1b from the slow bending of the contours at these extremal values. As energy does not need to be exactly conserved in an open system, scattering is also possible slightly outside the interval $[q_{\text{min}},k_{\text{max}}]$, with a width set by the linewidth $\gamma$. This characteristic series of peaks and dips in the steady-state momentum distribution appears to be a promising experimental signature of Beliaev-Landau scattering processes in a novel context of non-equilibrium quantum fluids. On the other hand, in the limit of small and large momenta ($k\rightarrow 0$ and $k\rightarrow\pi$) Beliaev-Landau processes are not allowed, so the Bogoliubov result is accurately recovered.

When trying to gain insight into the nature of out-of-equilibrium Beliaev-Landau scattering, it is worthwhile to take a closer look at the scattering matrix $M_{k,q}$, shown in Fig. 3. In addition to a central peak as a consequence of nonresonant decay, the contour representing energy and momentum conservation (see Fig. 1b)) is clearly manifested as a band of enhanced magnitude of $M_{k,q}$.

To clarify this, we take a step back and go again to the basis of Bogoliubov operators $\hat{\chi}_{k}$. By pursuing transformation (12) consistently, we find a closed set of equations equivalent to (V.1-28), but in terms of the $\hat{\chi}_{k}$. Although the full evaluation is much more cumbersome, as a consequence of the appearance of various products of the $u_{k}$ and $v_{k}$ transformation functions, one easily sees that the third-order correlation function must be of the form

[TABLE]

where the $A_{k,q}$ and $B_{k,q}$ are coefficients of order one which result from the combination of the Bogoliubov $u_{k},v_{k}$ factors corresponding to the different terms originating from the factorization of the fourth-order correlation functions in the Bogoliubov basis.

From the denominator of expression (31) one readily concludes that Beliaev-Landau scatterings are concentrated around the energy-conserving contours from (20). As the hierarchy of correlations (V.1-28) is built in the basis of the $\hat{\phi}_{k}$ operators, scatterings to negative energy states are also possible through the Bogoliubov transformation (12). These contours can be obtained by setting $\omega_{k}\rightarrow-\omega_{-k}$ and/or $\omega_{q}\rightarrow-\omega_{-q}$ in (20) and are visible as less pronounced bands of enhanced matrix elements in Fig. 3.

While the Bogoliubov approximation is a consistent expansion beyond mean-field that captures corrections which scale as $\sim U\psi_{0}^{2}$, we now conclude from (31) that terms scaling as $\sim U\psi_{0}/\sqrt{L}$, the next order in the expansion, are included with the present method. In particular, we have shown that the redistribution of occupation numbers is caused by quasiresonant Beliaev-Landau scattering. In our framework, this effect is accounted for in (25) and (26) by the back-reaction of the third-order correlation matrices, which capture the scatterings, on the momentum distribution and anomalous correlation.

VI Experimental signatures of Beliaev-Landau scatterings

To make the theoretical analysis discussed in the previous sections more concrete, we dedicate this section to providing some clear experimental indications on how the small Beliaev-Landau signal can be extracted in a realistic experimental context. We consider a one-dimensional cavity array consisting of a chain of coupled semiconductor microcavities, such as presented in Ref. Amo and Bloch (2016). The cavities, $L$ in total, are positioned at a distance $\Delta x=1$ $\mu$m from each other and are irradiated by a laser with frequency $\hbar\omega_{L}=1.6$ eV. Photons in the cavities have an average lifetime of $20$ ps, corresponding to a linewidth of $\hbar\gamma=33$ $\mu$eV and display a single-photon nonlinearity of $U=3.3$ $\mu$eV, such that $U\approx 0.1\gamma$. Furthermore we set $J=30\hbar\gamma\approx 1m$eV, $\Delta=-10\hbar\gamma\approx-330\mu$eV and assume an average number of photons per cavity of $n_{0}=100$, such that $Un_{0}\approx 10\gamma\approx 330\mu$eV. This is a case that we have already discussed in the theoretical analysis presented in Sec. V (see Fig. 2b red dashed line and Fig. 3b)).

As we have illustrated previously, the most straightforward approach to observe a signature of Beliaev-Landau scattering in the cavity array is to measure the momentum distribution of the quantum fluctuations and observe the characteristic series of peaks and dips. With the proposed parameters, we predict that a deviation of about $2\%$ from the Bogoliubov result (19) can be observed around the minimal final-state momentum of Beliaev decay $q_{\text{min}}$ (see Fig. 2b, red dashed line). The momentum distribution can be detected through an angle-resolved measurement of the far-field emission, as sketched in Fig. 4a): a photon with (adimensional) in-plane momentum $k$, will fly out of the cavity array at an angle given by $\sin{\theta_{k}}=ck/(\omega_{L}\Delta x)$, with $c$ the speed of light in vacuum Carusotto and Ciuti (2013).

For the proposed setup, we have that all quantum fluctuations can be detected by restricting the field of view to a cone of aperture $\theta_{\text{max}}\approx 23^{\circ}$. Importantly, the dominant signal of the condensate mode at $k=0$ is concentrated about the perpendicular axis and can be filtered out through post-selection.

The momentum-space density of photons escaping from the cavity array is approximately given by $\frac{d\Phi}{dk}=L\,n_{k}\gamma/(2\pi)$ where $2\pi/L$ is the momentum-space separation between adjacent modes for an array of $L$ cavities. For an array of $L=128$ cavities, an angular resolution of $\delta k=0.025\,(2\pi)$ larger than the $k$-space mode separation but well smaller than the width of the Beliaev features, and $n_{k}\approx 0.1$ around $q_{\text{min}}$ (see Fig. 2b) we expect a significant photon flux of about $\Phi=1.5\cdot 10^{10}\,$s*-1*. The number of photon clicks per time unit is then given by $N=\varepsilon_{\text{eff}}\Phi$ with $\varepsilon_{\text{eff}}$ some overall efficiency factor incorporating uncontrolled photon losses and detection efficiency. The signal can be integrated in time until a sufficient amount of photons is collected.

Of course, as seen in Fig. 2b, the experimental signal from Beliaev-Landau processes is enhanced with a larger nonlinearity. Although experimentally challenging, a stronger nonlinearity can in principle be achieved by reducing the size of the microcavities or by increasing the excitonic fraction of polaritons Rodriguez et al. (2017). Another more speculative possibility is to use the platform of superconducting circuits, where high nonlinearities are naturally achieved Fitzpatrick et al. (2017).

A crucial point of concern is that the Beliaev-Landau peaks, being rather small in size, can be washed away by a sufficient amount of disorder. In particular, when a small random potential $V_{i}$ is applied, for instance by variations of the cavity resonance $\omega_{c}$ from site to site, it will perturb the momentum distribution and imprint additional peaks. We can estimate that the disorder amplitude has to satisfy $\sqrt{\langle V_{i}^{2}\rangle}\lesssim\hbar\omega_{k}^{\text{peak}}\sqrt{\delta n_{k}^{\text{peak}}/n_{0}}\approx 3$ $\mu eV$, with $\omega_{k}^{\text{peak}}$ the frequency of the mode at the Beliaev-Landau peak and $\delta n_{k}^{\text{peak}}$ the height of the peak (on the order of $2\cdot 10^{-3}$, see the red line in the inset of Fig. 2b) ). In Appendix D we provide more details about the derivation of this estimation. In a recent experiment with a setup similar to ours, a standard deviation of about $30\mu$eV for the disorder potential was reported Baboux et al. (2017), a factor of about 10 larger than required for our estimations. However, given that the origin of the Beliaev-Landau peaks is different in nature than the disorder background, there are two additional strategies one can employ to isolate them.

First of all, the disorder peaks are different for each realization of a cavity array, while the Beliaev-Landau signal should not depend on this. If one therefore fabricates many copies of the same cavity array on the same sample, averaging over the different copies will cancel out effects from disorder, while the Beliaev-Landau peaks are left in place. If the different copies are positioned adjacent to one another, this amounts to displacing the laser beam from one array to the next.

Secondly, the exact position of the Beliaev-Landau peaks has a well-defined dependence on the mean-field parameters $Un_{0}$, $\Delta$ and $J$, as plotted in Fig. 1c). As a simple example, one could e.g. try to follow the shift of Beliaev-Landau peaks while varying the detuning $\Delta$ by modifying the laser frequency $\omega_{L}$ according to (6).

Finally, an alternative and conceptually more sophisticated strategy to observe Beliaev-Landau scattering processes is proposed in Fig. 4, where we present a sketch of a possible optical setup to directly measure the third-order correlator $M_{k,q}$ (see Fig. 3). The measurement would consist of detecting subtle correlations between the relative phase of the emissions at $k$ and $q$ and a homodyne measurement on the $k-q$ emission mixed with the coherent pump. The detection of simultaneous clicks in detector 1 and the difference signal of detectors 2 and 3 provides a measurement of the quantity

[TABLE]

where the second step is obtained after omitting all correlations that are not momentum-conserving, since they must be zero in a spatially uniform sample. The phases $\theta$ and $\phi$ are introduced by the two phase shifters in the setup and allow for the measurement of different quadratures of $M_{k,q}$. Any deviation from zero of the quantity (32) at non-vanishing angles $k,q,k-q\neq 0$ would be a manifest indication of the non-Gaussian nature of the cavity output field and would provide an indication of quasiresonant Beliaev-Landau scattering.

As high-order interference experiments of this kind go beyond standard quantum optical set-ups, a quantitative study of the expected signal and noise for a realistic experimental setup lies outside the scope of this work and will be the subject of a future study.

VII Conclusions

In this work we have theoretically studied the effect of Beliaev-Landau processes in a coherently driven fluid of light in a one-dimensional array of weakly nonlinear optical or microwave cavities. In contrast to the equilibrium case, where one typically looks at the decay of additional excitations externally generated in the fluid, here characteristic and experimentally accessible signatures of the Beliaev-Landau processes are identified in observable properties of the nonequilibrium steady-state.

Remarkably, the momentum distribution (visible in the angular distribution of the far-field emission pattern) shows a characteristic series of peaks and dips, which we attribute to the absence of detailed balance in an out-of-equilibrium setup. Also the higher-order correlators of the field (visible as non-Gaussian features in the photon statistics of the emitted light) exhibit nontrivial features stemming from quasiparticle interactions. Supported by our estimations, we expect that the predicted signal is within the reach of state-of-the-art experimental setups of coupled-cavity arrays with semiconductor microcavities or superconducting circuits.

From the theoretical point of view, our results pinpoint unexpected limitations to the use of the truncated Wigner method to describe scattering processes between quasi-particles. Given the importance of the TWA method as a tool for numerical studies of quantum fluctuation phenomena, future work will address improved schemes to overcome these difficulties.

The calculations are performed by truncating the hierarchy of correlations of the driven-dissipative Bose-Hubbard model to one order beyond Bogoliubov, i.e. by including the third-order correlation functions and employing a consistent truncation and factorization scheme for the higher-order correlation functions. Future work will extend this technique to spatially inhomogeneous configurations presently of great interest in the context of analog models of gravity.

Acknowledgements.

MVR gratefully acknowledges support in the form of a Ph. D. fellowship of the Research Foundation - Flanders (FWO) and hospitality at the BEC Center in Trento. WC and MW acknowledge financial support from the FWO-Odysseus program. IC was funded by the EU-FET Proactive grant AQuS, Project No. 640800, and by Provincia Autonoma di Trento, partially through the project “On silicon chip quantum optics for quantum computing and secure communications (SiQuro)”.

Appendix A TWA without Beliaev-Landau channels

In the main text we have pointed out how Beliaev-Landau scattering processes lie at the basis of the failure of the truncated Wigner method. To motivate this statement better, we present a TWA simulation of a system which does not contain resonant Beliaev-Landau channels that fulfil condition (20). As illustrated in Fig. 1, the Bogoliubov spectrum (14) determines the contour of resonant third-order scattering processes. By modifying the mean-field parameters $Un_{0}$, $J$ and/or $\delta$ we can enter into a regime where no resonant third-order scattering exists. Here we simply choose to replace the value $J=30\gamma$ that was used throughout the main text with $J=10\gamma$, so to have $\Delta<\Delta_{0}$ ( see Fig. 1c)).

As can be seen in Fig. 5, the absence of on-shell Beliaev-Landau channels leads to a much better agreement with the prediction of the Bogoliubov approximation and does not suffer from unphysical negative densities. The small deviation from the Bogoliubov distribution can probably be attributed to the nonresonant scattering of the $1/2$ vacuum noise.

Appendix B Derivation of the correlation hierarchy

With ansatz (7) from the main text, we find the equation of motion for the quantum fluctuations $\hat{\phi}_{k}$ (8). By repeatedly applying the product rule, one can obtain the equations of motion for the correlation functions of the quantum fluctuations. An alternative, completely equivalent approach would be to evaluate \partial_{t}\langle\hat{O}\rangle=i\big{\langle}\big{[}\hat{H},\hat{O}\big{]}\big{\rangle}+\text{tr}\big{\{}\hat{O}\mathcal{D}[\hat{\rho}]\big{\}}. Due to spatial homogeneity, only momentum-conserving operator products are included in this construction.

For a general correlation function $C=\left\langle\prod_{k}\hat{\phi}_{k}^{\dagger a_{k}}\hat{\phi}_{k}^{b_{k}}\right\rangle$ of order $N$, with N=\sum_{k}\big{(}a_{k}+b_{k}\big{)}, one can derive the following recurrence relation

[TABLE]

Here we adopted the notation, following Ref. Casteels et al. (2016)

[TABLE]

Up to third order, the explicit evaluation of the expression above yields the following equations of motion for the correlators.

• •

*First order: * A finite value for the zero-momentum component is found

[TABLE]

• •

*Second order: * We find for the density of fluctuations

[TABLE]

Likewise, for the anomalous averages

[TABLE]

• •

*Third order: * For the third-order correlation functions we can derive the equations of motion in the same way

[TABLE]

and

[TABLE]

Through its equations of motion, a correlator of order $N$ couples to correlators up to order $N+2$. In principle, this hierarchy of equations continues to infinite order if the number of particles is not conserved. Therefore, to obtain a closed set of equations, the hierarchy must be truncated in some way.

Appendix C Different truncation schemes and comparison

We briefly discuss two different truncation schemes and motivate the choice for the one presented in the main text.

• •

The hard cutoff (HC): The most straightforward approach is to set all correlation functions of order $N$ higher than $N_{c}$ (i.e. orders $N_{c}+1$ and $N_{c}+2$) to zero in the equations of motion. The major benefit of employing this truncation scheme is that, by construction, it produces a linear system of equations, which is relatively easily solved numerically. This allows one to evaluate the result as a function of $N_{c}$ for small enough system sizes. As discussed in Ref. Casteels et al. (2016), this approach is expected to be efficient when the number of excitations is small. More precisely, by pursuing this truncation scheme, one is implicitly assuming a small occupation of all modes $n_{k}<1$, otherwise the factorizable part of the correlation functions quickly grows as a function of $N$ and the calculation may not converge at high $N_{c}$. Nevertheless, even when $n_{k}<1$ is satisfied, no a priori assumptions can be made about the connected part of the correlation functions. Therefore one must always verify convergence of the result by increasing $N_{c}$. Following Ref. Casteels et al. (2016), we refer to this truncation scheme as the ‘hard cutoff’.

• •

The factorized cutoff (FC): The method employed in the manuscript differs in a few aspects from the one described above. First of all, we do not use the linearized equation for the first-order correlator $\langle\hat{\phi}_{0}\rangle$ (33), but integrate in time the full Gross-Pitaevskii equation with back-reaction terms, i.e. Eq. (V.1), along with the different correlation functions. As explained in the main text, this has the advantage of having $\langle\hat{\phi}_{0}\rangle=0$ by definition, making correlators up to order $3$ de facto connected, because their factorizable part vanishes. On the other hand, the field $\psi_{0}$ in equations (25)-(28) is now time dependent, so that the system of equations is no longer linear.

Furthermore, we perform two different approximations to close the system of equations at order $N=3$, by consistently including the back-reaction of the $N=4$ and $N=5$ in the equations for the $N=2$ and $N=3$ correlation functions. While the normal parts of higher order correlators were set to zero in the HC scheme, we start here by including them in factorized form

For a general $N=4$ correlation function, bearing in mind that $\langle\hat{\phi}_{0}\rangle=0$, we find that it can be written as,

[TABLE]

and likewise for other fourth-order correlators. The subscript $c$ (first line) denotes the connected, non-factorizable part of a correlator and is neglected in our truncation scheme. In the equations of motion for the second-order correlation function the factorization of the $N=4$ correlators from (36) and (39) into products of $N=2$ correlators yields the Hartree-Fock-Bogoliubov-like terms in (25)-(26). In turn, for the third-order, the $N=4$ correlators, given in (42)- (45), produce the drive terms $F^{(M,R)}_{k,q}$ from (29)-(30). Note that in the latter we have omitted factorizations of the form $\sim\frac{U\psi_{0}}{\sqrt{L}}c_{q}\sum_{l}n_{l}\delta_{k,0}$, and similar terms with $\sim\delta_{k,q}$ and $\sim\delta_{q,0}$, which drive the diagonal terms of $M(R)_{k,q}$. They are not related to Beliaev-Landau scatterings and we have checked that they merely give a negligible extra shift to $\psi_{0}$ and slightly renormalize the value of $n_{k}$ and $c_{k}$ in $k=0$, while leaving points at $k\neq 0$ essentially unaffected.

Also the fifth-order correlator, entering in the equations of motion for the third-order correlator, can be approximated by its factorizable form, which produces a total of ten different products of 2nd and 3rd order correlators. Two different groups of terms arise with this procedure. For instance, the first of the three terms entering on line (42) is

[TABLE]

The dots indicate more terms of the same kind, with summations over a $N=3$ (first line) or a $N=2$ correlator (second line). Hence we conclude that, after gathering all those terms, they can be captured by making changes of the kind

[TABLE]

The second approximation consists of neglecting the corrections coming from the factorized fifth-order correlator, which is justified by assuming that the condensate density is much larger than the density of fluctuations. We have evaluated all these fifth-order contributions and verified that their influence on the second-order and third-order correlation function is negligible for the parameters that are used.

To check the consistency of the method we used in the main text, we compare it in Fig. 6 and Fig. 7 with the HC truncation scheme at higher orders for a chain of only 10 sites. This allows us to obtain results within a reasonable computation time for truncation orders up to $N_{c}=6$ in the HC scheme. The factorized method from the main text appears to agree very well with the HC scheme for $N_{c}=5$, even though we only included connected correlators up to order 3. From this we conclude that, at least in the parameter regime that we have considered, it is a good approximation to neglect both the connected $N=4$ and the full $N=5$ correlator. On the other hand, we see a large deviation from the HC with $N_{c}=4$, even though both methods include correlators up to the same order $N=4$. This can be attributed to the inaccuracy with which the fourth-order is evaluated in the HC scheme, i.e by bluntly neglecting all higher-orders. Therefore the inclusion of the fourth order in factorized form directly, as was employed in our FC scheme, turns out to be a much better approximation than obtaining it through a HC scheme with $N_{c}=4$.

From Fig. 6 and Fig. 7 we deduce that results obtained with different truncation schemes start to deviate from each other at $U/\gamma=0.5$. This limits the range of parameters in which our approach provides a quantitatively accurate description. While for stronger $U/\gamma$, the error due to truncation of higher orders terms becomes increasingly important, still at $U/\gamma=0.5$ we see that all truncation schemes reproduce at least qualitatively the same result; in particular, the $N_{c}=3$ factorized cutoff still reasonably well agrees with the $N_{c}=5$ HC scheme, and even with $N_{c}=6$. We have therefore chosen to use the $N_{c}=3$ factorized cut-off approximation for the study of Beliaev-Landau processes in larger systems, even though there is a small, but non-negligible quantitative deviation from higher-order HC schemes.

Keeping in mind that higher-order HC schemes are numerically cumbersome and can not be applied to large systems, our main motivation for sticking to the $N_{c}=3$ factorized cut-off approximation is that, to our knowledge, no efficient methods exist to simulate driven-dissipative quantum dynamics in a system with intermediate interactions ($U=0.5\gamma$) and a large particle number (128 cavities, each containing $20$ photons on average). A possible alternative approach would be to develop a variational method with matrix product operators in the spirit of Mascarenhas et al. (2015); Cui et al. (2015), in which the matrix product state of the quantum excitations is determined self-consistently by coupling it back to the coherent condensate. This will hopefully be the subject of future work.

We therefore conclude that, even though we constructed a hierarchy in terms of correlation functions up to order 3 only, the accuracy is comparable to (or even better as) higher-order methods in the HC scheme. Obviously, the reduced numerical complexity in the developed truncation scheme realizes a significant computational speedup as compared to these higher-order methods, thus allowing us to tackle much larger systems and address the physically most relevant questions.

Appendix D The influence of disorder

Given the Fourier transform $V_{k}$ of a random potential, $V_{j}=\frac{1}{\sqrt{L}}\sum_{k}V_{k}e^{ikj}$ that is applied to the cavity array. We find that the mean field follows the equation of motion

[TABLE]

In the linear-response regime, the non-uniform polariton field can be formulated as $\psi_{j}=\psi_{0}+\frac{1}{\sqrt{L}}\sum_{k}\delta\psi_{k}e^{ikj}$. After substitution in (49) and collecting terms up to linear order in $\delta\psi_{k}$ and $V_{k}$, we derive a linear set of equations for each mode

[TABLE]

with the response matrix

[TABLE]

and $\epsilon_{k}$ given in (9). Solving (50) yields the response of the density distribution to the disorder potential in the linear regime

[TABLE]

with $\omega_{k}$ given in (14). Since all energy scales are larger than $\gamma$ and $\omega_{k}\approx\epsilon_{k}$ for the purposes of this qualitative analysis, we can further approximate \delta n_{k}\sim n_{0}\big{(}V_{k}/\omega_{k}\big{)}^{2}. For white uncorrelated noise it therefore follows that roughly $\sqrt{\langle V_{j}^{2}\rangle}\lesssim\omega_{k}^{\text{peak}}\sqrt{\delta n_{k}^{\text{peak}}/n_{0}}$ if we want the peaks of disorder to be smaller than the peaks of Beliaev-Landau scattering.