Spatial correlations in driven-dissipative photonic lattices

We study the nonequilibrium steady-state of interacting photons in cavity arrays as described by the driven-dissipative Bose-Hubbard and spin-$1/2$ XY model. For this purpose, we develop a self-consistent expansion in the inverse coordination number of the array ($\sim 1/z$) to solve the Lindblad master equation of these systems beyond the mean-field approximation. Our formalism is compared and benchmarked with exact numerical methods for small systems based on an exact diagonalization of the Liouvillian and a recently developed corner-space renormalization technique. We then apply this method to obtain insights beyond mean-field in two particular settings: (i) We show that the gas--liquid transition in the driven-dissipative Bose-Hubbard model is characterized by large density fluctuations and bunched photon statistics. (ii) We study the antibunching--bunching transition of the nearest-neighbor correlator in the driven-dissipative spin-$1/2$ XY model and provide a simple explanation of this phenomenon.


Introduction
In recent years interacting photonic lattices have emerged as a versatile platform for the study of many-body phenomena out of equilibrium [1][2][3][4][5]. First prototype quantum simulators have been realized experimentally based on cavity and circuit QED technologies [6][7][8][9][10][11][12]. The increasing experimental interest in assembling cavities to form lattices is also a strong motivation to develop novel theoretical tools. The key object governing the dynamics of such driven-dissipative systems is typically the Liouvillian superoperator [13], which describes the dynamical evolution of the system density matrix ρ through a master equation. Solving the master equation exactly is a formidable numerical task [14]. While exact diagonalization and quantum-trajectory algorithms [15][16][17][18] allow to successfully address this problem for small system sizes, large scale numerical methods based on matrix-product-states (MPS) [19][20][21][22][23][24] are typically limited to one dimension (1D). Recently developed methods such as the corner-space renormalization technique [25] may provide a promising alternative also in two dimensions (2D). On the other hand, decoupling mean-field theory, which is correct in infinite lattice dimensions, is a simple yet valuable tool to gain a first insight into the qualitative physics at work. It has been successfully applied to various lattice models such as the Bose-Hubbard and Jaynes-Cummings-Hubbard model [26][27][28][29][30][31] as well as related spin models [32][33][34]. Recent efforts to improve on the mean-field approximation include perturbative [35,36], projective [37], cluster [38], variational [39] and equations-ofmotion approaches [40].
Here, we develop a systematic expansion around the decoupling mean-field solution of the Lindblad master equation in powers of the inverse dimensionality parameter z 1 (with z being the number of nearest neighbors in a lattice). Such an expansion accounts for quantum fluctuations in a systematic way and provides access to a whole new class of observables, i.e., spatial correlation functions. For systems in (quasi-) equilibrium, which are fully described by the Hamiltonian alone, the z 1 expansion has a diagrammatic interpretation in terms of linked-clusters and was used to calculate the ground-state and elementary excitations of Fermi-Hubbard [41], Bose-Hubbard [42] and Jaynes-Cummings-Hubbard [43] models. In the nonequilibrium context, this Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
technique was employed in [44,45] to calculate quenched dynamics of atoms in optical lattices and in [46] to characterize the transition from low to high density phases in a driven, dissipative Rydberg system.
In this work, we expand on previous efforts by developing a method to solve for the density matrix in a selfconsistent way. We calculate the nonequilibrium steady-state (NESS) of the driven-dissipative Bose-Hubbard model (BHM) up to second order in z 1 and show that the self-consistency condition substantially improves the results by comparing to exact diagonalization (ED) in 1D and the corner-space method in 2D. We then apply our method to two specific problems: (i) we calculate the compressibility of the driven-dissipative BHM and show that the photonic gas-liquid transition is characterized by largely enhanced density fluctuations with bunched photon statistics; (ii) we study the antibunching-bunching transition of the driven-dissipative spin-1/2 XY model in 1D and 2D and provide a simple explanation based on a dimer model.
The remainder of the paper is structured as follows. In section 2, we introduce two models for interacting photons in cavity arrays, the driven dissipative Bose-Hubbard and the spin-1/2 XY model. In section 3, we discuss the self-consistent z 1 expansion and benchmark our method by comparing with numerical results based on ED and the corner-space renormalization technique. In section 4, we address the effects of site-site correlations in the gas-liquid transition of the driven-dissipative BHM. In sections 5, we study the drivendissipative spin-1/2 XY model to discuss the antibunching-bunching transition in one and 2D. In section 6 we summarize the results of the paper and provide an outlook for future work.

Model
We investigate the steady-state of the coherently pumped and dissipative BHM describing photons hopping on a lattice of nonlinear cavities with local coherent pump and decay. The lattice Hamiltonian reads Here, each site i is coherently pumped with strength f as described by the last term in h i , which is expressed in terms of the bosonic operator a i and the associated number operator n a a i i i = † . In a frame rotating with the drive frequency d w the cavity frequency is renormalized to d c w w D = -, while U is the local Kerr nonlinearity. The second term in H describes the hopping to z nearest-neighbor cavities with amplitude J J; ij =the additional factor z 1 in(1) ensures that the bandwidth of the photon dispersion is J 2 , independent of z, and guarantees a regular limit z  ¥. The dissipative dynamics for the density matrix ρ is accounted for via Lindblad's master equation, fluctuations yielding quantitatively correct results in a large parameter range even for small lattice sizes. While we focus here on the BHM and XYM, the technique is rather generic and applicable to a wide range of driven, dissipative lattice models with limited range hopping.
We start by defining the reduced density matrices of one lattice site tr Such a scaling of correlations is known from the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy of statistical mechanics [47], with the difference that it here applies to lattice sites instead of particles. Starting from (3), we obtain the equation of motion for the reduced density matrices up to order z 1 2 [44,45], i.e., Above, we introduced the notation ij ) [ ] as in [44,45]. In the mean-field limit of infinite coordination number (z  ¥) all connected density matrices are zero and one only needs to solve(6a), which is nonlinear and can have multiple solutions. However, in order to account for spatial correlations, one needs to evaluate the density matrix to higher order in z 1 and also solve the equations of motion for the connected density matrices. In a first step, we make use of the scaling hierarchy In the following, we substantially improve this first approximation by keeping explicitly all underlined terms to second order in z 1 in the system of equations above, i.e., by neglecting only the third order term in (6c) ( ijkk c r ¢ ). We then solve the coupled system of equations in a self-consistent way taking the following steps: correspond to the first step, which was explained in the previous paragraph. In order to implement self-consistency we now explain the second step, i.e., (iv) insert ijk c r back in (6b) and obtain an updated ; plug ij c r in (6a) and get a new i r . In (iv)-(v) all the underlined terms are kept. Starting from the updated i r , the procedure (ii)-(v) is iterated till convergence is reached. This yields a solution of the hierarchy equations (6) correct to second (2nd) order z 1 2 , i.e., with an error on the density matrix of order z 3  -( ). Without the steps (iii)-(iv) the solution of the hierarchy equations is correct to first (1st) order z 1 , i.e., with an error on the density matrix of order z 2  -( ).
Step (i) alone is correct to zeroth order and equivalent to a Gutzwiller mean-field (MF) decoupling of the hopping term in the Hamiltonian (1). The sequence of steps performed in this self-consistent scheme is illustrated in figure 1(a).
In figure 1(b) we show how the numerical complexity of the method scales with the number of lattice sites N and compare to an exact numerical solution of the master equation. The size of the full Hamiltonian is given by ) and thus increases exponentially with the number of lattice sites (see blue dots in figure 1(b)). An exact solution of the master equation is then obtained by writing the density matrix as a vector r of density matrix elements with length M 2 such that (3)  . A diagonalization of the (non-hermitian) Liouvillian then yields in general complex eigenvalues and eigenfunctions which fully determine the exact solution.
Let us now estimate the computational complexity of the z 1 expansion assuming a translationally invariant density matrix with periodic boundary conditions. In this case, corresponds to the order of the expansion. Consequently, the computational effort scales only polynomially with the number of lattice sites. In figure 1(b), we compare the dimensions of the Liouvillian operators for n p =1. For example, the calculation of the XYM on a square lattice with 7×7 sites would involve a very large Liouvillian operator with M 10 15 » , which would be far beyond sparse ED methods and even stochastic techniques based on quantum trajectories [15] where M≈10 6 forms an upper limit. On the other hand, the density matrix of such a large system can be easily computed using the coordination number expansion to second order even on a standard laptop computer (e.g., see results in table 1). If translational symmetry is broken, the complexity of the method increases with M N n 1 , as one cannot assume that (i) the matrix i r is site-independent and that (ii) the connected density matrices depend only on the distance between the sites. However, the increase in complexity does not compromise the polynomial scaling of M with the number of lattice sites N.
In table 1, we compare the z 1 expansion with (i) the ED method for a 1D chain and (ii) with numerical data available for a 2D square lattice from the so-called corner-space renormalization method (CM) developed in [25]. The CM is a numerical algorithm which uses the exact solution of the master equation for a small lattice and extrapolates it to larger system sizes. At each extrapolation step of the algorithm, two small lattices are merged to form a larger one, while truncating the basis of the joint Hilbert space to a small number of most probable states (i.e., the corner-space). For better comparison, we chose the same parameters as in [25] for both dimensions. Shown are results for the photon density with the number of lattice sites N, for the z 1 expansion (curves) and exact diagonalization (ED, symbols) for a system with maximally one particle per site n p =1. The scaling for larger n p is similar.
describing instantaneous (zero time delay) correlations between sites i and j. The latter is measurable in a Hanbury Brown-Twiss setup [48,49]. The average in (7) and (8) is taken with respect to the NESS of equation (6) with 0 r = . For the parameters considered in table 1, our self-consistent z 1 expansion improves the mean-field result (MF) substantially and agrees well with the exact numerical findings in both dimensions. In 1D, we find quantitative agreement with the exact result up to the second and the third decimal for weak to moderate hopping rates (J k ). Small discrepancies start to show up for larger hopping rates (J 3k , see rows marked with an asterisk * in table 1) and strong site to site correlations (g 2 01 2( ) ). Such a behavior is expected as the z 1 expansion treats the non-local hopping term perturbatively. In 2D, the comparison with the CM method works similarly well. The convergence of the method after a few iteration steps is demonstrated exemplarily in figure 2. The self-consistency scheme considerably improves the first and second order results of the z 1 expansion and converges rather fast. In the following two sections we apply this technique to study the gas-liquid transition in the BHM and the antibunching-bunching transition in the XYM.

BHM: gas-liquid transition
In this section, we study the gas-liquid transition of photons as described by the driven-dissipative BHM [28,39,50]. The gas (liquid) phase is characterized by low (high) photon densities of the NESS. The transition between the two phases can be driven by the coherent pump parameter f/U at fixed detuning U D . For a single cavity, an exact solution provides a smooth crossover between the two phases when the pump strength is Table 1. Density n, on-site and nearest-neighbor correlators g j 0 2 ( ) with j=0, 1 for the BHM. Shown are results obtained from the z 1 expansion applied to a 1D array with 6 sites and cutoff n p =2 (a) and to a 2D square lattice (b) with 4×4 sites (U 20 k = , n p =3), 3×3 sites (U 10 k = , n p =5), and 7×7 sites (U 1 k = , n p =4), where n p is the local photon cutoff. The z 1 results are compared with ED in (a) and with data from the corner-space method (CM) [25] in (b). We have used the parameters 5 k D = , f 2 k = , and J 1 k = . The rows marked with an asterisk * show results for a large hopping J 3 k = where the agreement is less favorable.  increased [51]. In the lattice case, decoupling mean-field theory predicts that the gas-liquid crossover transforms into a hysteretic transition beyond a critical value of the intercavity hopping J J c = . The phasediagram in the f-J plane (at fixed detuning U D ) is shown in figure 3(a) with the critical point (blue) at J c .
Interestingly, one finds that the critical hopping J c is modulated as a function of the detuning U D and exhibits a series of lobes, see figure 3(b). The lobe structure is a manifestation of a quantum commensuration effect which favors the hysteretic transition over a smooth crossover whenever the drive frequency corresponds to a mphoton resonance at U m 1 2 + D = [29].
Unfortunately, the lobe structure is particularly hard to calculate with exact numerical methods, because it requires a high single-cavity photon number cutoff n p to capture the physics of multi-photon resonances. This is why quantum trajectory simulations in [29] were initially limited to 6 sites. However, despite the small system size, these simulations strongly substantiate the mean-field prediction: below the critical point (J J c < ), trajectories of each cavity switch independently and at random times between gas and liquid states; this behavior changes drastically beyond the mean-field critical point (J J c > ), where all cavities of the array switch synchronously between gas and liquid phases. In the following, we take a closer look at the gas-liquid transition and analyze compressibility and spatial correlations of the steady-state beyond mean-field using the z 1 expansion described in the previous section. First, we study density fluctuations via the compressibility  figure 3(a). At weak drive f U 1  (gas phase), we find K 1 » as predicted by the mean-field approximation (solid line) and the z 1 expansion. Consequently, the gas phase is well described by a spatially uncorrelated, coherent state with g 1 for all sites j. At large drive f U 1  (liquid phase), the prediction of the mean-field approximation, i.e., K 1 2 » , also agrees well with the results obtained from the z 1 expansion. In fact, the value K 1 2 » can be derived analytically from the single-cavity limit (J = 0), where g m 1 1 00 2 » - ( ) and n m 2 » at the m-photon resonance [50]. We conclude that the effect of the lattice dimension is marginal deep in the gas and liquid phase, where the physics is well described by mean-field theory. However, the crossover region is characterized by strongly enhanced density fluctuations beyond mean-field. In particular, quantum fluctuations due to z 1 corrections in 1D as well as 2D strongly increase the compressibility with respect to the mean-field result. We attribute these enhanced fluctuations to the impending bistable behavior, see also [34]. Our z 1 results are also consistent with the quantum trajectory calculations in [29], which show that synchronization effects already appear below the critical mean-field value J c .
Making use of the z 1 expansion we also calculate the spatial correlation functions of the NESS for J J c < . ) extend further out in the lattice with a larger correlation length, signaling the crossover between gas and liquid phases. This clustering of excitations is consistent with the coherent super-cavity formation as revealed by the quantum trajectory simulations in [29]. Away from the compressibility peak (A and C) photons at different sites are mostly uncorrelated. The symbols at j=0 indicate the mean-field values of the on-site correlator g 00 2 ( ) , which significantly differ from the 1D results only at f U 0.35 = (B, square). Similar outcomes are obtained for the 2D lattice. We note, that the local Hilbert space cutoff n p (maximum photon number per cavity) required in figure 4 is n p =6, which would imply a huge Liouvillian operator of size M 10 21 » in ED for the 2D case with 5×5 sites.
In summary, in this section we have shown with the z 1 method that bunched site-site correlations extend over many lattice sites and largely enhance density fluctuations in the gas-liquid crossover regime of the drivendissipative BHM (1). The low computational cost of the method allowed us to obtain insight for large lattices in 1D and 2D also in a regime of large photon numbers. Unfortunately, it is difficult to analyze the hysteretic transition within the z 1 expansion since the self-consistent approach does not always converge in this region of the phase diagram. In the next section, we will rather focus on the strongly-correlated regime U  ¥ where the BHM (1) is mapped to the spin-1/2 XYM (4).

Spin-1/2 XY model: antibunching-bunching transition
In this section, we investigate the driven-dissipative spin-1/2 XYM in (4). In particular, we study the antibunching-bunching transition of the nearest neighbor correlator as a function of the detuning Δ, which was recently predicted in [52] using large scale MPS simulations. In the following, we (i) provide a simple and analytic explanation of the transition based on a minimal model of two coupled spins (dimer), (ii) reproduce exact numerical results with the self-consistent z 1 method to high accuracy and (iii) go beyond the MPS method by also studying the 2D case.
Before considering large lattices in 1D and 2D, it is instructive to focus on a simpler model consisting of a dimer of two coupled, driven-dissipative spins, i.e., a system described by the XYM in (4) with N=2 sites and the associated four basis states gg ge eg ee , , , ñ ñ ñ ñ {| | | | }. Figure 5 shows the photon amplitude (homodyne signal)  deep in the gas phase, see [29].