Truncated correlation hierarchy schemes for driven-dissipative multimode quantum systems

We present a method to describe driven-dissipative multi-mode systems by considering a truncated hierarchy of equations for the correlation functions. We consider two hierarchy truncation schemes with a global cutoff on the correlation order, which is the sum of the exponents of the operators involved in the correlation functions: a 'hard' cutoff corresponding to an expansion around the vacuum, which applies to a regime where the number of excitations per site is small; a 'soft' cutoff which corresponds to an expansion around coherent states, which can be applied for large excitation numbers per site. This approach is applied to describe the bunching-antibunching transition in the driven-dissipative Bose-Hubbard model for photonic systems. The results have been successfully benchmarked by comparison with calculations based on the corner-space renormalization method in 1D and 2D systems. The regime of validity and strengths of the present truncation methods are critically discussed.


I. INTRODUCTION
Recently there has been a strong interest in driven-dissipative manybody systems, particularly nonequilibrium photonic systems. These systems can be realized with various platforms, including semiconductor microcavities and superconducting circuits (see for example Refs. [1][2][3] for recent reviews). Typically, the steady-state of such systems can be far from equilibrium due to the competition between driving and dissipation, which can not be neglected. In many regimes, they can be described by a Lindblad master equation for the manybody density-matrix. The possibility of coupling many of these modes (e.g., in lattices of coupled nonlinear resonators) is currently receiving a lot of attention for the realization of quantum many-body states of light.
Due to the inherent out-of-equilibrium nature of these systems many of the numerical approaches developed for systems at equilibrium are not necessarily applicable. These systems are typically computationally more challenging since the total number of photons is typically not conserved resulting in a much larger effective Hilbert space. Moreover, they are in general in a mixed state described by a density matrix while a closed system at low temperature is described by the groundstate wavefunction. As a result, exact diagonalisation approaches are only applicable for relatively small systems. Examples of numerical tools for larger systems are the TEBD algorithm [4,5] and a DMRG-type variational approach for the steady-state [6,7]. The underlying structure of both these techniques is based on the so-called matrix product states: such approaches are efficient in 1D, but encounter difficulties for 2D systems [8,9]. For the description of 2D lattices, the corner-space renormalization method has been recently introduced [10].
In this paper we investigate truncation schemes for the infinite hierarchy of equations coupling correlations functions. Both schemes are based on a global cutoff on the correlation order, that is the sum of the exponents of the operators involved in the correlation functions. In particular, one truncation scheme corresponds to an expansion around the vacuum (hard cutoff) and is most appropriate when the number of excitations is small. The second truncation scheme is based on a self-consistent 'soft' cutoff, which corresponds to an expansion around coherent states and applies when the number of excitations per site is large. The numerical implementation of this scheme is straightforward with respect to the more advanced numerical approaches discussed above. Moreover, it allows to explore regimes that are very hard to reach with techniques that are based on a cutoff in the number of photons, such as exact diagonalization. The presented numerical scheme is very versatile and allows to study various theoretical models with a wide range of applications. Recently a scheme with a hard cutoff has been applied in Ref. [11] to obtain a perturbative expansion in the hopping for an infinite 2D driven-dissipative Bose-Hubbard lattice. Note that the truncation schemes presented in the present paper are different from the approach in Ref. [11] where a local cutoff for the exponents of the operators is enforced for each site separately.
As an example we apply the method to the drivendissipative Bose-Hubbard model.
In particular the physics of the anti-bunching to bunching transition is examined: such an effect was investigated for or a single resonator (e.g. in Ref. [12]) and explored more recently in Ref. [13] for a lattice of dissipative two-level systems. Both a 1D chain as well as a 2D square lattice of resonators are investigated and the results are compared with the results obtained from the corner-space renormalization method. Two sets of system parameters are considered, one resulting in a relatively small photon density, ideally suited for the scheme with a hard cutoff, and one with a relatively high density which is well captured by the soft cutoff truncation. In all cases, the presence of additional modes is shown to reduce both the bunching and the anti-bunching and for finite systems an oscillation of the second order correlation functions as a function of the cavity-drive frequency detuning is shown. As the size of the system is increased these oscillations become less pronounced and eventually vanish in the thermodynamic limit, resulting in an overall bunching to anti-bunching transition at zero detuning, corresponding to the middle of the single-particle energy band.
This paper is structured as follows: in Section II we introduce the theoretical approach in a general way; in Section III, we specialize the theory to the driven-dissipative Bose-Hubbard model; in Section IV the numerical results are presented and critically discussed. Finally, in Section V, the conclusions are presented together with an outlook regarding the perspectives of the method.

II. THEORETICAL FRAMEWORK AND DESCRIPTION OF THE METHOD
A. The Lindblad master equation If a system is weakly interacting with a large bath, such that retardation effects can be neglected and the Markov approximation is valid, the degrees of freedom of the environment can be traced out to get an effective Lindblad master equation for the reduced density operatorρ of the system (we seth = 1): whereĤ is the system Hamiltonian andĈ j are the jump operators that describe the interaction with the environment. From the master equation (1) an expression can be obtained for the time derivative of the expectation value of an arbitrary operatorÔ that acts on the system: Although this is a general picture, from now on we will focus on the class of systems consisting of interacting bosons hopping on a lattice. Within the formalism of second quantization, we introduce the bosonic annihilation operators {â i } and creation operators â † i , where i indicates the mode index. Such operators satisfy the bosonic commutation rules [â i ,â j ] = â † i ,â † j = 0; and â i ,â † j = δ i,j with δ i,j the Kronecker delta. We now consider a HamiltonianĤ that consists of a number conserving partĤ Sys and a drive termĤ Drive . The particle conserving part typically contains one-and two-particle terms, namely: where the sums run over all the modes of the system. The matrix h (1) i,j can be used to describe the kinetic energy and external potentials while the matrix h (2) i,j introduces interactions.
For the drive term a coherent pump field is considered: with F i the drive amplitude for mode i. Note that this is written in the frame rotating at the drive frequency (in the rotating rame, the system Hamiltonian conserves the general form in Eq. (3)).
For the environment we assume that the bosonic excitations can escape to a bath at temperature zero. This is described by jump operators that are proportional to the annihilation operatorsĈ j = √ γ jâj , with γ j the corresponding dissipation rate. It is straightforward to also introduce jump operators that populate the modes incoherently (proportional to the creation operators), which can describe a non-resonant pump or the thermal excitations of the environment. Also other types of dissipation processes can be considered, such as pure dephasing [14].

B. Equations of motion for the correlation functions
All the observables that are local in time can be written in terms of the correlation functions Applying expression (2), we get the equations of motion for the correlation functions. For the one-particle term of the number conserving part (3), we find: The two-particle part gives: The part describing the coherent drive results in: Finally, for the term corresponding to the dissipation in expression (2), we need: Summing all the terms, we find for the time derivative of an arbitrary correlation function: where the following notation has been introduced: To make the notation more compact and valid for all the integer indexes, note that formally there are a few terms with negative exponents, but they are of course equal to zero because the pre-factors multiplying them vanish as it can be readily verified. Eq. (10) corresponds to a linear set of equations for the correlation functions. Due to the two-particle partĤ (2) Sys in the Hamiltonian each correlation function is a function of higher order correlation functions, thus resulting in a hierarchy with an infinite set of coupled equations. If one is interested in the properties of the steady-state the time derivatives can be set to zero.
In order to find a numerical solution for the correlation functions, a truncation scheme of the infinite hierarchy of equations must be introduced. We will introduce a global cutoff on the correlation order N C , i.e., we will consider only the equations of motion for correlation functions with i n i ≤ N C and i m i ≤ N C . Note that this is different from a local cutoff, as for example considered in Ref. [11], for which all the correlation functions with n i ≤ N C and m i ≤ N C are retained, for all i. In the case of multiple modes the global cutoff allows to capture higher order correlation functions with a smaller set of equations.
It is clear that if we consider the correlation functions up to the order N C , in their equations a dependence on higher order correlation functions remains. In the following, we propose two complementary schemes to deal with them.

C. Global 'hard' cutoff
The first rather straightforward possibility corresponds to setting these correlation functions equal to zero, which we denote as a hard cutoff. Namely, we impose This cutoff is expected to be adequate in the limit of a small number of excitations.
This scheme result in a total number of correlation functions N (HC) for which the equations of motion have to be solved, given by: Hence, for a given cutoff order N C the number of equations grows according to a power law as a function of the size of the system. This is a consequence of the global cutoff: for a local cutoff, instead, the number of equations grows exponentially with the system size. This shows that truncation schemes based on a global cutoff can be very efficient for problems that involve a low order of correlation.

D. Global self-consistent 'soft' cutoff
When the steady-state is in a coherent state, the correlation functions decouple: Hence, in a regime where the number of excitations is large and the system is not so far from a coherent state, we are led to consider the following 'soft' cutoff condition: for i n i > N C or i m i > N C . Note that this equation gives an expression for or i m i > N C in terms of lower order correlation functions and thus closes the system of equations. This scheme can be thought of as an expansion in correlations around a coherent state, thus allowing the study of correlations in regimes with a high photon density that are unattainable with numerical approaches based on a cutoff in number of photons. Note that if the cutoff is set to N c = 1 this reduces to the non-equilibrium Gross-Pitaevskii approximation [1]. A consequence of the soft cutoff is that the set of equations of motion becomes nonlinear and has to be solved self-consistently for { â k }. We note that this approach with the soft cutoff is related to a cumulant expansion, as for example pursued in Ref. [15]. A consequence of the soft cutoff is that the set of equations of motion becomes nonlinear and has to be solved self-consistently for { â k }. In practice this is obtained by using the Gross-Pitaevskii prediction (corresponding to N C = 1) for the { â k } as an initial guess. These values are updated by solving the system of equations with the 'old' values for { â k } in (14) which results in a new set { â k }. This procedure is then iterated until convergence is reached. For the regimes considered in this paper this is obtained already after a few iterations.
In this case, the total number of correlation functions N (SC) that have to be solved self-consistently is: This can be understood by following the same reasoning as presented in the previous subsection for a hard cutoff and adding the number of Eqs. (14) that are needed to close the system of equations. Similar to Eq. (13) for a hard cutoff, the dominant contribution to Eq. (15) for a large system size N S is: Since this approach results in a larger set of equations with respect to the hard cutoff scheme it is limited to smaller systems and/or a smaller cutoff N C .

III. DRIVEN-DISSIPATIVE BOSE-HUBBARD MODEL
As an example we apply the presented method to the driven-dissipative Bose-Hubbard model. The corresponding system Hamiltonian is (in the frame rotating at the pump frequency): The parameters are the detuning ∆ = ω p − ω c between the pump and the cavity frequency, the on-site interaction strength U and the hopping parameter J. The sum i,j is restricted to the nearest neighbours. Note that this Hamiltonian indeed has the form considered in Eq. (3). Furthermore, in the rotating frame the coherent drive is described by Eq. (4).
We will consider one-dimensional chains and twodimensional arrays, both with periodic boundary conditions. The length of a chain is denoted as L = N a, with N the number of sites and a the lattice constant. The linear (U = 0) system Hamiltonian is diagonal in the reciprocal space, with k n = n2π/L the allowed wave numbers with n ∈ {0, 1, ..., N − 1} and contains resonances at ∆ = −2J cos(ak n ). For the two dimensional clusters, N x sites in the x-direction and N y sites in the y-direction are considered (denoted as a N x × N y cluster). The allowed k-vectors are: (k nx , k ny ) = 2π/a × (n x /N x , n y /N y ), with n x ∈ {0, 1, ..., N x − 1} and n y ∈ {0, 1, ..., N y − 1}. In this case the linear resonances are at ∆ = −2J cos (ak nx ) − 2J cos ak ny . For a homogeneous system with uniform drive only the modes corresponding to n = 0 (1D) or (n x , n y ) = (0, 0) (2D) are driven, the other modes can only get populated through nonlinear interactions.
In the following, we will calculate the local density n i and the normalised second-order correlation function g (2) i,j , defined as: Since we will consider translationally invariant systems the value of i is not important and will be omitted for the density, i.e. n i = n.

IV. NUMERICAL RESULTS
In this section we present results for the steady-state of the driven-dissipative Bose-Hubbard model with a homogeneous drive (F i = F ) in 1 and 2 dimensions with periodic boundary conditions. Note that in 1D the periodic boundary conditions can be realised experimentally by arranging the resonators in a ring geometry [16,17]. We will consider two complementary sets of system parameters that are ideally suited for the two presented truncation schemes.
The convergence of the results as a function of N C has been verified numerically for the smaller system sizes. In-deed the considered regimes are well described by the presented truncation schemes already with a cutoff N C = 2. For the largest considered system sizes, the accuracy of the results have been successfully benchmarked by comparison to calculations based on the corner-space renormalization method [10] (convergence of the results has been carefully checked by increasing the dimension M of the corner space, with relative errors below 1%. ) i,i (b) and the nearest neighbour g (2) i,i+1 as a function of the detuning ∆ (in units of γ) for the driven dissipative Bose-Hubbard 1D chain with periodic boundary conditions. For the density there is no significant difference for different lattice sizes and only a very small difference with the linear spectrum (shown as a dashed curve but very hard to distinguish). For the correlation functions the length of the chain is denoted in the figures. The system parameters are: J = γ, U = 1γ, F = 0.1γ. The black dots are results obtained with the corner-space renormalization method [10] for the largest considered system of 16 coupled cavities (with a corner space dimension M = 1000; relative errors below 1%).

Hard Cutoff
In Fig. 1 the results are presented for 1D chains with up to 16 sites as a function of the detuning ∆, for the system parameters J = γ, U = γ and F = 0.1γ. Due to the small density, the hard cutoff truncation scheme has been considered. For the local density in Fig. 1 (a) only the result for a single cavity is presented since it is practically independent of the chain length and furthermore deviates only slightly from the linear spectrum (also presented in Fig. 1 (a)).
The on-site normalized second order correlation function g (2) i,i is presented in Fig. 1 (b) . In contrast to the behavior of the density in Fig. 1 (a), the second-order correlation function dramatically depends on the chain length. For sufficiently negative detuning anti-bunching is observed and for sufficiently positive detuning we find bunching. This behavior is already well-known for a single cavity with a anti-bunching to bunching transition around the resonance [12]. As the length is increased an overall attenuation of the bunching and anti-bunching is found. Moreover, an oscillating behavior appears at the detunings corresponding to the linear resonances at ∆ = 2J cos(ak n ). Note that since the coherent drive is homogeneous only the n = 0 mode is pumped and all other modes are only populated through scattering and are completely absent for the linear system. The amplitude of these oscillations decreases as the number of modes increases, indicating a bunching to anti-bunching transition around ∆ = 0 in the thermodynamic limit N → ∞. Note that the results are in excellent agreement with the predictions of the corner-space renormalization method (black dots).
In Fig. 1 (c) the nearest neighbour normalised second order correlation function g (2) i,i+1 is presented as a function i,i (b) and the nearest neighbour g (2) i,i+1 (c) second order correlation functions as a function of the detuning ∆ for the driven dissipative Bose-Hubbard lattice in 1D with periodic boundary conditions. For the density there is practically no difference for different lattice sizes but the nonlinear shift with respect to the linear spectrum (dashed curve) can be clearly seen. For the correlation functions the size of the lattice is denoted in the figures. The system parameters are: J = γ, U = 0.05γ, F = 1.2γ. The black dots are results obtained with the corner-space renormalization method for the largest considered system of 9 coupled cavities (with a corner-space dimension M = 500; relative error below 1%).
of the detuning. This reveals that g (2) i,i+1 < 1 for large detuning |∆| > 2J, with a minimum just outside the edges of the linear energy band at ∆ = ±2J. For small detunings we find g (2) i,i+1 > 1. Similar to the results for the local g (2) i,i in Fig. 1 (b), the overall values of g (2) i,i+1 decrease as the chain length is increased. Also here, an oscillating behavior is observed for the finite systems at the detunings corresponding to the linear eigenmodes, with a decreasing amplitude as the length is increased.
In Fig. 2 (a) the normalised second order correlation function g (2) i,i+j is presented as a function of the distance j for three values of the detuning (∆/γ = −2.2, 0 and 2.2). The results show an overall decay to g (2) i,i+j → 1 for j → ∞, with superimposed oscillations. Similar oscillations in the non-local g (2) i,i+j have been predicted for various one dimensional driven-dissipative nonlinear photonic system [18][19][20] and they typically indicate a regular ordering of the photons.

Soft Cutoff
We now perform a similar analysis as presented in the previous subsection, but for the system parameters J = γ, U = 0.05γ, F = 1.2γ. This corresponds to a relatively large photon density and is ideally suited for the soft cutoff truncation scheme. As discussed before, due to the soft cutoff we are limited to smaller chain sizes as compared with the hard cutoff and we consider chains with up to 9 sites. In Fig. 3 (a) the density is presented as a function of the detuning for which the dependence on the chain size is again negligible. Note that the resulting density is much larger than the one found in Fig.  1 (a) and the nonlinear shift of the resonance peak with respect to the linear system is now clearly visible. In Fig. 3 (b) and (c) the local and nearest neighbour normalised second order correlation functions are presented. Due to the larger photon density the system is closer to a coherent state and the normalised correlation functions are much closer to 1 with respect to the results in Fig.  1. However, the same qualitative behavior is found as in Fig. 1. In Fig. 2 (b) the dependence of the g (2) i,i+j on the distance j is presented which is also qualitatively similar to the result in Fig. 2 (a). Note that also here an excellent agreement is found with the results obtained with the corner-space renormalization method for the largest system size (black dots).

Hard Cutoff
Here, we again consider the system parameters J = γ, U = 0.1γ, F = γ and the hard cutoff truncation scheme, as in Sec. IV A 1, but now for 2D clusters with up to 16 sites arranged in a square geometry (N x = N y ). In Fig. 4 (a) the density is presented as a function of the detuning ∆ which is again practically independent of the size of the cluster and indistinguishable from the linear result. In Fig. 4 (b) and (c) the local and nearest neighbour second order correlation functions are presented as a function of the detuning ∆. This again exhibits the same qualitative behavior as found for 1D chains with an anti-bunching to bunching transition that becomes less pronounced as the cluster size is increased and superimposed oscillations at the linear resonances. These oscillations again diminish for increasing size indicating a anti-bunching to bunching transition around ∆ = 0 in the thermodynamic limit N x , N y → ∞. For the nonlocal g (2) i,j ( i, j denotes nearest neighbours) we again find g (2) i,j < 1 far detuned from the linear resonances (|∆| > 4J), with a minimum just outside the edge of the linear energy band ∆ = ±4J; instead, g (2) i,j > 1 for detuning within the linear energy band. i,i (b) and the nearest neighbour g (2) <i,j> (c) second order correlation functions as a function of the detuning ∆ for the driven dissipative Bose-Hubbard lattice in 2D with periodic boundary conditions. The geometry of the considered clusters are characterized as Nx * Ny. For the density there is practically no difference for different lattice sizes and only a small difference with the linear spectrum (shown as a dashed curve, practically indistinguishable). The system parameters are: J = γ, U = 1γ, F = 0.1γ. The black dots are results obtained with the corner-space renormalization method for the largest considered system of 16 coupled cavities (with a corner space dimension M = 1500; relative error below 1%) .

Soft Cutoff
As before in Sec. IV A 2, the system parameters J = γ, U = 0.05γ, F = 1.2γ are considered and the system is studied with the soft cutoff truncation scheme for clusters with up to 9 sites arranged in a 3 × 3 geometry. In Fig.  4 (a) the density is presented which is again practically independent of the size of the cluster and much larger with respect to the previous subsection with a clear nonlinear shift of the resonance peak. Also the local and nearest neighbour second order correlation functions in Fig. 4 (b) and (c) exhibit the same qualitative behavior as discussed before. Note that the results presented in both Figs. 3 and 4 for the largest considered system sizes are in excellent agreement with the predictions of the corner-space renormalization method. Notice the reported results correspond to a mean-number of photons per site exceeding 5. Such larger number of bosons per site cannot be described by brute-force integration of the master equation in the Fock basis, because the Hilbert size dimension becomes prohibitively large. i,i (b) and the nearest neighbour g (2) <i,j> (c) second order correlation functions as a function of the detuning ∆ for the driven dissipative Bose-Hubbard lattice in 2D with periodic boundary conditions. The geometry of the considered clusters are characterized as Nx * Ny. For the density there is practically no difference for different lattice sizes and a clear nonlinear shift with respect to the linear spectrum (shown as a dashed curve). The system parameters are: J = γ, U = 0.05γ, F = 1.2γ. The black dots are results obtained with the corner-space renormalization method for the largest considered system of 9 coupled cavities (with a corner-space dimension M = 800; relative error below 1%).

V. CONCLUSIONS AND OUTLOOK
We have presented a numerical approach suitable for the description of few-body correlations in systems out of equilibrium, based on the equations of motion for the correlation functions. Two possible complementary truncation schemes based on a global cutoff were discussed, one ideally suited for small densities and the other for large densities. As an example we have applied this scheme to the driven-dissipative Bose-Hubbard model in both oneand two-dimensions with periodic boundary conditions and examined the fate of the bunching to anti-bunching transition in extended lattices. An oscillating behavior is observed for the local second-order correlation function g (2) i,i as a function of the detuning for intermediate system sizes. The oscillations tend to disappear as the size of the system increases. An extrapolation of this behav-ior leads to an overall bunching to anti-bunching transition around ∆ = 0 in the thermodynamic limit. For the nearest neighbor correlations we found g (2) <i,j> < 1 for detuning outside of the single-particle energy band and g (2) <i,j> > 1 for detunings within the band. For 1D chains also the dependence of the correlations on the distance was studied, revealing an oscillating behavior typically observed for these systems. The presented results were successfully benchmarked with the corner-space renormalization method [10] by comparing a representative set of points.
The presented results reveal that there are physical problems for which the considered truncation schemes can be more efficient with respect to other numerical approaches. An important advantage with respect to other truncation schemes is that for a given correlation order the number of equations increases as a power law versus the number of sites in contrast to the typical exponential dependence for truncation approaches with a local cutoff. Furthermore, the numerical implementation is rather simple and straightforward with respect to more advanced numerical schemes such as the ones based on matrix product operators (MPO) in 1D and the cornerspace renormalization method (1D and 2D). The results presented in this paper for the largest considered systems are obtained in a computational time of the order of minutes on a desktop computer. Although with MPO approaches and the corner-space renormalization method it is possible to tackle larger lattices and larger correlations, the present truncation schemes with global cutoff can be much faster (a couple of orders of magnitude for the lattices of moderate size considered in this paper).
The presented method can be generalized to various other models and physical situations. One can consider multiple modes per site, allowing for example the study of the driven-dissipative Jaynes Cummings model [21,22], an opto-mechanical coupling [3,23] or the photon polarization [17]. Also dissipative spin systems such as for example the dissipative XYZ model [24,25] or collective phases with hard-core bosons [26] are a possible extension. Other possibilities are the consideration of spatially inhomogeneous systems with site dependent parameters (for example to study transport properties [27]), other lattice geometries (for example to study the effect of geometric frustration [19,28]), other boundary conditions (for example closed or twisted boundary conditions [16] or a (cluster) Gutzwiller mean field approach [25,29,30]) or an incoherent population of the modes (for example thermal excitations or an incoherent drive). Also regimes with higher order correlations can be examined which however limits the size of the considered lattices. If one is interested in the time-dependent correlation functions the coupled differential equations of motion (10) can be solved numerically with the presented truncation schemes. This allows to examine dynamical phenomena such as the dynamical hysteresis typically encountered in the optical bistability regime [31].