Phase diagram of incoherently driven strongly correlated photonic lattices

TL;DR

This paper theoretically investigates the nonequilibrium phases of a driven-dissipative photonic lattice, revealing a phase transition from a Mott-like state to a coherent phase driven by hopping rates and symmetry breaking.

Contribution

It introduces a phase diagram for incoherently driven strongly correlated photonic lattices using a Gutzwiller approach and finite-size simulations, highlighting a novel nonequilibrium phase transition.

Findings

Identifies a second-order phase transition with symmetry breaking.

Shows the transition is driven by commensurability, not nonlinearity.

Validates mean-field results with matrix product and corner-space methods.

Abstract

We explore theoretically the nonequilibrium photonic phases of an array of coupled cavities in presence of incoherent driving and dissipation. In particular, we consider a Hubbard model system where each site is a Kerr nonlinear resonator coupled to a two-level emitter, which is pumped incoherently. Within a Gutzwiller mean-field approach, we determine the steady-state phase diagram of such a system. We find that, at a critical value of the inter-cavity photon hopping rate, a second-order nonequilibrium phase transition associated with the spontaneous breaking of the $U(1)$ symmetry occurs. The transition from an incompressible Mott-like photon fluid to a coherent delocalized phase is driven by commensurability effects and not by the competition between photon hopping and optical nonlinearity. The essence of the mean-field predictions is corroborated by finite-size simulations obtained with matrix product operators and corner-space renormalization methods.

Equations

$$\hat{H}_{ph}=\sum_{i}\Bigl{(}\omega_{c}\ \hat{a}^{\dagger}_{i}\hat{a}_{i}+U\ \hat{a}^{\dagger}_{i}\hat{a}^{\dagger}_{i}\hat{a}_{i}\hat{a}_{i}\Bigr{)}-J\sum_{\braket{i,j}}\hat{a}_{i}\hat{a}^{\dagger}_{j},$$

$${\hat{H}}_{at}$$

$${\hat{H}}_I$$

$$\dot{\rho}=-{\rm i}\bigl{[}\hat{H},\rho\bigr{]}+\mathcal{L}[\rho],$$

$$\mathcal{L}[\rho ]=\sum _i\left(\frac{{\Gamma }_l}{2}\mathcal{D}[{\hat{a}}_i;\rho ]+\frac{\gamma }{2}\mathcal{D}[{\hat{\sigma }}_i^{-};\rho ]+\frac{{\Gamma }_p}{2}\mathcal{D}[{\hat{\sigma }}_i^{+};\rho ]\right),$$

$${\omega }_N=N{\omega }_c+N(N-1)U.$$

$$\frac{{\Gamma }_{em}^0}{{\Gamma }_l}\gg 1,\frac{{\Gamma }_{em}^0{\Gamma }_p^2}{{\Gamma }_l{U}^2}\ll 1,$$

$$n=\frac{4{\Gamma }_p{\Omega }_R^2}{({\Gamma }_p+\gamma +{\Gamma }_l)\left({\Gamma }_l({\Gamma }_p+\gamma )+4{\Omega }_R^2\right)}.$$

$$n=\frac{{\Gamma }_p}{{\Gamma }_{tot}}-\frac{{\Gamma }_p+\gamma }{{\Gamma }_{tot}}\eta +\mathcal{O}({\eta }^2),$$

$$\Delta n^2=\frac{{\Gamma }_p(\gamma +{\Gamma }_l)}{{\Gamma }_{tot}^2}-\frac{({\Gamma }_p+\gamma )({\Gamma }_p-\gamma -{\Gamma }_l)}{{\Gamma }_{tot}^2}\eta +\mathcal{O}({\eta }^2).$$

$${\rho }_{\text{MF}}=\underset{i}{⨂}{\rho }_i,$$

$$\dot{\rho}_{i}=-{\rm i}\bigl{[}\hat{H}_{\text{MF}},\rho_{i}\bigr{]}+\mathcal{L}_{i}[\rho_{i}],$$

$${\rho }_i^{SS}\simeq |1,\uparrow ⟩⟨1,\uparrow |.$$

$$⟨\hat{a}⟩=|⟨\hat{a}⟩|{\mathrm{e}}^{-\mathrm{i}{\omega }_{L}t},$$

$$W(\alpha )=\frac{2}{\pi }\mathrm{Tr}\left[{\rho }_{ph}^{SS}\hat{D}(\alpha ){\mathrm{e}}^{\mathrm{i}\pi {\hat{a}}^{†}\hat{a}}{\hat{D}}^{†}(\alpha )\right],$$

$${\Gamma }_{em}^0\frac{{\left({\Gamma }_p/2\right)}^2}{{\left({\omega }_{at}-{\omega }_c+z{J}_c^{\mathrm{HC}}\right)}^2+{\left({\Gamma }_p/2\right)}^2}\simeq {\Gamma }_l.$$

$$\frac{{\Gamma }_{em}^0}{{\Gamma }_l}=\frac{4{\Omega }_R^2}{{\Gamma }_p{\Gamma }_l}\gg 1,$$

$$z{J}_c^{\mathrm{HC}}\simeq \sqrt{\frac{{\Gamma }_p}{{\Gamma }_l}}{\Omega }_R.$$

$$c(i,j)=⟨{\hat{a}}_i^{†}{\hat{a}}_j⟩.$$

$$|c(i,i+r)|\propto {\mathrm{e}}^{-r/\lambda },$$

$$\lambda \propto \frac{1}{|\ln (J/{\Gamma }_p)|},$$

$$\hat{a}$$

$${\hat{a}}^{†}$$

$$\hat{H}=\frac{{\omega }_c}{2}({\hat{\Sigma }}^z+1)+\frac{{\omega }_a}{2}({\hat{\sigma }}^z+1)+\frac{{\Omega }_R}{2}({\hat{\Sigma }}^x{\hat{\sigma }}^x+{\hat{\Sigma }}^y{\hat{\sigma }}^y).$$

$$\dot{\theta }=\mathrm{i}[\hat{H},\hat{\theta }]+\tilde{\mathcal{L}}[\rho ],$$

$$\tilde{\mathcal{L}}[\rho ]=\sum _i\left(\frac{{\Gamma }_l}{2}\tilde{\mathcal{D}}[{\hat{a}}_i;\rho ]+\frac{\gamma }{2}\tilde{\mathcal{D}}[{\hat{\sigma }}_i^{-};\rho ]+\frac{{\Gamma }_p}{2}\tilde{\mathcal{D}}[{\hat{\sigma }}_i^{+};\rho ]\right),$$

$$⟨{\hat{\Sigma }}^x⟩$$

$$⟨{\hat{\Sigma }}^z⟩$$

$$⟨{\hat{\sigma }}^z⟩$$

$$⟨{\hat{\Sigma }}^x{\hat{\sigma }}^z⟩$$

$$⟨{\hat{\Sigma }}^x{\hat{\sigma }}^x⟩$$

$$⟨{\hat{\Sigma }}^x{\hat{\sigma }}^y⟩$$

$$⟨{\hat{\Sigma }}^z{\hat{\sigma }}^z⟩$$

$$\hat{H}\to \hat{H}-\frac{zJ}{2}({\hat{\Sigma }}^x⟨{\hat{\Sigma }}^x⟩+{\hat{\Sigma }}^y⟨{\hat{\Sigma }}^y⟩).$$

Full text

I Introduction

Interacting many-particle systems driven away from the thermodynamic equilibrium exhibit a number of interesting features in both the classical cross1993 ; vicsek ; zia and quantum petruccione ; eisert2015 regimes. Our understanding of such systems is limited by the fact that nonequilibrium can emerge in a broad variety of forms, making it difficult to develop a unitary picture. In contrast, in equilibrium situation, very general paradigms allow us to predict the behavior of large ensembles of particles fisher1974 ; fisher1998 ; pathria ; sachdev .

Nonequilibrium dynamics can be realized both in isolated quantum systems subject to quench dynamics polko2011 as well as in open systems under the effect of driving and dissipation. The investigation of nonequilibrium phenomena in extended open quantum systems is a complex problem: theoretical efforts to develop analytical techniques weimer2015 ; sieberer2016 and numerical methods orus2008 ; mascarenhas2015 ; cui2015 ; finazzi2015 ; montangero2014 ; ksh2016 have been carried on over the years. Another strategy involves the experimental realization of a tunable and well controllable quantum simulator feynman ; buluta2009 which mimics the behavior of the real system under consideration.

The quantum simulation of many-body systems in equilibrium conditions have been implemented in very different contexts ranging from ultracold atoms in optical lattices bloch2008 ; blatt2012 to trapped ions monroe2013 . The impressive experimental advances of the last decade allowed to extend this idea to the nonequilibrium realm using Rydberg atoms and trapped ions muller2012 , exciton-polariton condensates carusotto2013 , cold atoms in cavities ritsch2013 , and arrays of coupled QED-cavities  houck2012 ; tomadin2010rev . Among them, coupled-cavities arrays are particularly appealing from the condensed matter physics perspective since they allow the simulation of archetypal (interacting) lattice models under nonequilibrium conditions reviews .

The emergence of collective behaviors and critical phenomena in these platforms have attracted an increasing interest over the years. Their phase diagram can be remarkably rich (see for example Refs. hartmann2010 ; lee2011 ; umucalilar2012 ; jin2013 ; lee2013 ; hoening2014 ; chan2015 ; wilson2016 ; ff2017 ): Exotic ordering schiro2016 and phases without an equilibrium counterpart may appear jin2016 . Only very recently, the build up of quantum correlations and the peculiar behavior at criticality started to be investigated torre2012 ; sieberer2013 ; marino2016 ; rota2016 . As it happened for their equilibrium counterpart, the possibility to engineer and manipulate complex many-body states would allow to study fundamental questions and obtain deep insight about the nature of phase transitions and spontaneous out-of-equilibrium ordering.

An important ingredient determining the dynamics in these systems is the competition between the on-site photon-photon interaction (mediated by the atoms) and the photon-hopping between neighbouring cavities: A large local interaction favours the formation of states with a fixed number of particles per site, while a large photon-hopping allows delocalization and enhances the density fluctuations. The early works in the literature hartmann2006 ; greentree2006 ; angelakis2007 explored this physics considering essentially the case of negligible photon losses. The large number of studies in this regime clearly showed the striking resemblance of the thermodynamic phases greentree2006 ; zhao2008 ; aichhorn2008 ; schmidt2009 as well as of the critical properties rossini2007 ; koch2009 ; hohenadler2011 with the celebrated Bose-Hubbard (BH) model (see also the reviews reviews and the references therein).

However, the unavoidable presence of photon dissipation and atomic relaxation affects the thermodynamic properties of these quantum simulators. The dissipative processes are counteracted by an external (coherent or incoherent) driving source which makes the long-time dynamics nontrivial: As a matter of fact, this is determined by the simultaneous interplay between Hamiltonian dynamics, dissipative processes and external driving. As a result, the scenario considerably enriches and the intrinsic nonequilibrium nature of these platforms emerges in different aspects, ranging from their dynamical response tomadin2010 and transport properties biella2015 ; angelakis2015 to the steady-state behavior nissen2012 ; leboite2013 ; biondi2016 . Furthermore, nonequilibrium effects have been highlighted in single planar cavities, where new complex shapes for the condensate wavefunction have been found and explained wouters2008 and novel dispersions for the Goldstone mode predicted wouters2007 . Even if the competition between the photon-hopping and on-site nonlinearity always plays an important role, the analogies with the underlying BH model are in general difficult to be found.

In this work we study the steady-state phases of an incoherently driven photonic lattice (see sketch in Fig. 1). Each site is represented by a lossy nonlinear cavity, where the photon mode is coupled to a two-level system, which is pumped incoherently lebreuilly2015 ; rivas2014 . A recent study jose2017 has shown that, via a non-Markovian pump scheme based on a reservoir of two-level systems with engineered spectral features, it is possible to stabilize Mott-like states of photons in spite of the losses. However, a nonequilibrium phase diagram is unknown for this class of systems where a photonic lattice is incoherently-driven. Here we determine the phase diagram within a Gutzwiller mean-field approach including both photon modes and two-level systems, showing the emergence of a second-order phase transition with a $U(1)$ symmetry breaking from a Mott-like incompressible fluid of light to a coherent delocalized phase. Our Gutzwiller theoretical predictions are consistent with finite-size numerical simulations obtained with two different methods (matrix product operators and corner-space renormalization).

The paper is organized as follows. In Sec. II we introduce our model, highlighting the main features of the incoherent driving scheme. In Sec. III we discuss the steady-state phase diagram of the system under the Gutzwiller mean-field approximation. In Sec. III.3 we go beyond the Gutzwiller approximation and we show the signatures of our findings in cavity arrays, comparing the mean-field results with finite-size numerical simulations. Finally, in Sec. IV we draw our conclusions and discuss some future perspectives.

II The Model

We consider a driven-dissipative BH model for photons in a $d$-dimensional array of QED-cavities (setting $\hbar=1$)

[TABLE]

where $\hat{a}_{i}$ ($\hat{a}^{\dagger}_{i}$) are bosonic photon annihilation (creation) operators associated with the $i$-th cavity of the chain with natural frequency $\omega_{c}$, which obey the canonical commutation relations (\bigl{[}\hat{a}_{i},\hat{a}^{\dagger}_{j}\bigr{]}=\delta_{ij},\bigl{[}\hat{a}_{i},\hat{a}_{j}\bigr{]}=0), $J$ is the hopping rate and $U$ sets the scale of the Kerr nonlinearity.

Each cavity is coupled to a two-level emitter which is pumped incoherently and provides a driving source for the array. The atomic evolution and the coupling to the cavities are ruled by

[TABLE]

where $\hat{\sigma}_{i}^{\pm}=(\hat{\sigma}^{x}_{i}\pm{\rm{i}}\hat{\sigma}^{y}_{i})/2$ and $\{\hat{\sigma}_{i}^{\alpha}\ |\ \alpha=x,y,z\}$ are the Pauli matrices acting on the $i$-th site.

The photon leakage from the cavities, atomic relaxation and pumping processes are taken into account by means of a master equation for the density matrix in the Lindblad form

[TABLE]

where $\hat{H}=\hat{H}_{ph}+\hat{H}_{at}+\hat{H}_{I}$ and

[TABLE]

with $\mathcal{D}[\hat{O};\rho]=[2\hat{O}\rho\hat{O}^{\dagger}-\{\hat{O}^{\dagger}\hat{O},\rho\}]$. A sketch of the system is provided in Fig. 1. In what follows we will be interested in the nonequilibrium steady-state (NESS) of this model, as determined by computing the long-time limit of Eq. (5) $\rho(t\to\infty)=\rho^{SS}$.

III Results

In this section, we present our main results for the steady-state properties of the considered system obtained by solving the the master equation (5) via several approaches. In Sec. III.1, we first consider the single-cavity case and discuss the photon number selection mechanism that can be achieved with the considered incoherent pump scheme for the two-level emitters. In particular, in the limit of large nonlinearity (hard-core photons), we provide the exact analytical solution for the steady-state density matrix. In Sec. III.2, we study the many-cavity case and explore its steady-state phase diagram within a Gutzwiller mean-field approach. Finally, in Sec. III.3 we compare our findings with exact finite-size simulations using two different techniques.

III.1 Photon number selection in a single cavity

The eigenvalues of the Hamiltonian in Eq. (1) for $J=0$ (labelled by their photon number) are

[TABLE]

This implies that the $N\to N+1$ transition has a frequency $\omega_{N+1,N}=\omega_{c}+2NU$. Choosing the emitter transition frequency to be resonant with the $N\to N+1$ transition (i.e., $\omega_{at}=\omega_{N+1,N}$), it is possible to obtain a NESS that is a mixed state (with no coherences) dominated by the $N+1$ photon state (and the atom in the excited state). While the specific conditions to obtain this for generic integer $N$ require a fine tuning of parameters, as discussed in full detail in lebreuilly2015 , the ones for the $N=1$ photon state have the simple form

[TABLE]

where $\Gamma_{em}^{0}=4\Omega_{R}^{2}/\Gamma_{p}$ . We checked this (numerically) focusing on the $0\to 1$ transition by solving the single-cavity master equation via diagonalization of the corresponding Liouvillian. In the following we will work in units of $\Gamma_{p}$. In Fig. 2 we show the steady-state value of the photon density $n=\braket{\hat{a}^{\dagger}\hat{a}}$ (where $\braket{\hat{O}}={\rm Tr}[\rho^{SS}\hat{O}]$ and ${\rm Tr}[\rho^{SS}]=1$) and its variance $\Delta n$ as a function of $U/\Gamma_{p}$ for different values of the cavity dissipation rate $\Gamma_{l}/\Gamma_{p}$. Moreover, we also show the purity of the density matrix $\mathcal{P}={\rm Tr}[(\rho^{SS})^{2}]$ and the population $\Pi_{1}=\bra{1,\uparrow}\rho^{SS}\ket{1,\uparrow}$, where $\ket{1,\uparrow}$ denotes the state with one photon in the cavity mode and the two-level system into its excited state. As highlighted in the right bottom panel of Fig. 2, it is possible to prepare the desired Fock state with arbitrary precision for large enough nonlinearity and small photon leakage.

When the nonlinearity is large enough to forbid the double occupation of a lattice site, Eq. (5) can be analytically solved for the steady-state (see Appendix A). At resonance ($\omega_{at}=\omega_{c}$) the photon density reads as

[TABLE]

Expanding Eq. (9) for small effective loss/gain ratio $\eta=\Gamma_{l}/\Gamma_{em}^{0}$ we obtain

[TABLE]

where $\Gamma_{tot}=\Gamma_{p}+\gamma+\Gamma_{l}$. The expression above for the photon population tells us that, in the regime of large nonlinearity, it is possible to stabilize in the steady-state a single-photon Fock state for small enough cavity and atomic dissipation rates. The number of photons fluctuates as

[TABLE]

What we want to do next is to study what the effect of a finite coupling between neighbouring cavities ($J\neq 0$) is.

III.2 Gutzwiller phase-diagram for cavity lattices

In this section, we consider a lattice of cavities, using the same kind of incoherent driving that we have analyzed for the single-cavity case. Due to the complexity of the problem, we perform a Gutzwiller mean-field approximation leboite2013 ; biondi2016 assuming a factorized ansatz for the global density matrix

[TABLE]

where $\rho_{i}$ is the density matrix of the $i$-th site. Inserting such ansatz into Eq. (5) and assuming the translational invariance ($\rho_{i}=\rho_{j}$ , $\forall i,j$) we get an effective master equation of the form

[TABLE]

where $\hat{H}_{\text{MF}}=\hat{H}_{i}+\hat{H}_{\mathcal{B}}$. Here $\hat{H}_{i}$ and $\mathcal{L}_{i}$ contain all the local terms of the full Hamiltonian $\hat{H}$ and of the superoperator $\mathcal{L}$ (see Eq. (6)) respectively, acting on the $i$-th site. The term $\hat{H}_{\mathcal{B}}=-zJ(\hat{a}_{i}^{\dagger}\braket{\hat{a}}+\rm{H.c.})$ (where $z$ is the coordination number of the lattice), takes into account the mean-field interactions of the site $i$ with its neighbours.

The steady-state is reached dynamically after integrating Eq. (13) by means of a fourth-order Runge-Kutta method. The mean-field is computed dynamically $\braket{\hat{a}}(t)={\rm Tr}[\hat{a}\ \rho_{i}(t)]$ and used to self-consistently evolve $\rho_{i}(t)$ until the steady-state is reached. This approach has been proven to be very effective to determine the phase diagram of the coherently driven BH model leboite2013 ; biondi2016 ; wilson2016 and similar systems jin2013 since the local pumping and decay drastically restrict the range of the correlations. We therefore expect that the Gutzwiller ansatz (12) is well suitable to capture the physics of our system.

In Fig. 3 we show the value of the mean-field order parameter $|\braket{\hat{a}}|$, the photon density $n$, its fluctuations $\Delta n$ and the compressibility $\mathcal{K}=\Delta n^{2}/n$ note00 in the steady-state as a function of the inter-cavity hopping rate in the hard-core limit ($U/J=+\infty$). The main findings of the Gutzwiller mean-field theory can be summarized as follows.

For $0\leq J<J_{c}^{\rm HC}$ (where $J_{c}^{\rm HC}$ denotes the critical hopping rate in the in the hard-core limit), we get a vanishing value of the mean-field order parameter $|\braket{\hat{a}}|=0$. As a consequence, each local steady-state density matrix can be approximately written as

[TABLE]

By construction, a vanishing value of the MF order parameter forces $\rho_{i}^{SS}$ to be the steady-state solution of the ME (5) for a single cavity. However, the short-range coupling induced by the photon hopping may play an important role. Consequently, in order to characterize the phase with the unbroken symmetry it is necessary to go beyond the mean-field theory. This has been done in Sec. III.3 where we show that in a range of $zJ/\Gamma_{p}$ compatible with the Gutzwiller prediction the number of photons remains very close to the unity with very small fluctuations. This Mott-like phase, is also characterized by an (almost) vanishing compressibility $\mathcal{K}$, analogously to what happens in equilibrium situations note01 .

At $J=J_{c}^{\rm HC}$, a second-order phase transition takes place. For $J>J_{c}^{\rm HC}$, the system enters a coherent delocalized phase characterized by the emergence of limit cycles, namely

[TABLE]

where $\omega_{L}$ depends on the system parameters. This transition is associated with the spontaneous breaking of the $U(1)$ symmetry possessed by Eq. (5) note02 . The steady-state density matrix becomes mixed (not shown) and the photonic population feels the transition: The number of photons starts to be significantly different from $1$ and its fluctuations become relevant. The transition from a Mott-like to a coherent phase can be measured by monitoring the behavior of the compressibility which becomes finite in the symmetry broken phase. For large values of $zJ/\Gamma_{p}$ the photon density remains finite approaching asymptotically $n=1/2$ and the steady-state density matrix becomes maximally mixed. We also note that the use of the translationally invariant ansatz (12) imposes that the condensate appears in the $\mathbf{k}=0$ mode note03 . This restriction is however fully justified as the pumping conditions chosen for the numerics of Sec. III.3 (see Fig. 7 and Fig. 8) explicitly favor condensation into the $\mathbf{k}=0$ mode.

Further information about the nature of the steady-state can be extracted by looking at the Wigner quasi-probability distribution

[TABLE]

where $\rho^{SS}_{ph}={\rm Tr}_{at}\left[\rho^{SS}\right]$ is the photonic reduced density matrix, $\hat{D}(\alpha)={\rm{e}}^{\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}}$ is the displacement operator and the prefactor ensures that $\int_{\mathbb{C}}d^{2}\alpha\ W(\alpha)=1$. As it is shown in Fig. 4, for $0\leq J<J_{c}^{\rm HC}$ the steady-state is approximately a one-photon Fock state and then the Wigner distribution is negative around $\alpha=0$ (see top panel) indicating strong nonclassicality haroche_book . For $J>J_{c}^{\rm HC}$ the $U(1)$ symmetry is spontaneously broken in the steady-state. The Wigner function thus becomes asymmetric and rotates at a frequency $\omega_{L}$ around $\alpha=0$ (see bottom panel). Let us note that also in this case $W(\alpha)$ is negative around $\alpha=0$ again indicating the nonclassical nature of the steady-state at low photon number per site.

The phase transition we observe is not related to the competition between photon hopping and nonlinearity since we explored a range of parameters such that the $zJ/U$ ratio remains very small. Increasing $J$ we increase the bandwidth ($2zJ$), therefore the photons start to be off-resonant with respect to the incoherent driving source provided by the atoms. As a consequence, the photons cannot be efficiently pumped into the array and their number is no longer commensurate to the system size. The presence of a significant number of empty sites allows photons to move along the lattice and to establish a long-distance coherence. As we will discuss in Sec. III.3, the transition can be characterized also in terms of spatial correlation functions: it is associated to the appearance of a long-range order in the first-order coherence function. Very recently, the transition between a Mott-like regime to a coherent one driven by the $J-U$ competition has been theoretically investigated in jose2017 by making use of specifically designed emission spectra based on a non-Markovian reservoir.

Also for the configuration under exam in the present work, in spite of its nonequilibrium nature it is possible to trace an analogy with the phase diagram of hard-core BH model at zero temperature schmid2002 . The ground-state of this model is a superfluid (in two dimensions) when the band is not completely filled (i.e., $|\mu/J|<z$, where $\mu$ is the chemical potential). Here, this commensurability driven transition is extended to the nonequilibrium realm and the emergent coherent phase ($|\braket{a}|\neq 0$) has a mixed nature.

In the hard-core regime, it is also possible to exploit the single-cavity exact solution to infer the structure of the phase diagram. In the unbroken symmetry phase ($\braket{\hat{a}}=0$), the MF master equation (13) always admits a stationary solution which corresponds to the single-cavity steady-state. For certain values of the parameters such solution becomes unstable and the system approaches a stationary state with $\braket{\hat{a}}\neq 0$. The stability analysis reveals that, for the typical values of the parameters we are interested in, the single-cavity fixed point gets unstable as the hopping rate is increased (see Fig. 5). The scaling of the instability point with the system parameters can be estimated by comparing the effective emission rate at the band boundary with the photon leakage rate

[TABLE]

In the regime where the effective pumping rate dominates over the photon losses, that is

[TABLE]

at resonance ($\omega_{at}=\omega_{c}$) we get

[TABLE]

The prediction of Eq. (19) (solid black lines in Fig. 5) well reproduces the phase boundaries in the regime (18). The dashed horizontal lines in Fig. 5 denote the $\Gamma_{em}^{0}/\Gamma_{l}=1$ threshold. The lasing condition $\Gamma_{em}^{0}/\Gamma_{l}\geq 1$ is necessary in order to have a significant population in the symmetry broken phase. For the details about the Gutzwiller mean-field stability analysis see App. B.

It is interesting to extend our study considering finite values of the nonlinearity. In Fig. 6 we show the mean-field steady-state phase diagram in the $U/\Gamma_{p}-zJ/\Gamma_{p}$ plane. For each value of $U$ we explored, we found a critical value of the tunneling rate $J_{c}(U)$ such that for $J>J_{c}(U)$ the $U(1)$ symmetry is spontaneously broken ($|\braket{\hat{a}}|\neq 0$) and the system exhibits limit cycles. The value of the critical hopping rate increases as $U$ is increased, approaching the value predicted in the hard-core limit $zJ_{c}^{\rm HC}/\Gamma_{p}=1.51$ (predicted by the Gutzwiller stability analysis).

III.3 Beyond the Gutzwiller approximation

In this section, we go beyond the Gutzwiller mean-field approximation employed in Sec. III.2 and present finite-size simulations using different numerical methods. The goal is to show what is the fate of the phases predicted by the Gutzwiller mean-field theory. Since the Hilbert space dimension increases exponentially with the system size, the exact diagonalization of the full Liouvillian becomes soon impracticable. Numerically, the present problem is also particularly challenging since there are several time scales which differ by orders of magnitude (in particular, the two-level and cavity dissipation rates compared to the incoherent driving rate). In order to compute the steady-state $\rho^{SS}$ of Eq. (5) we use two methods: an algorithm based on the matrix product operators (MPO) formalism vestraete2004 ; zwolak2004 and the corner-space renormalization method finazzi2015 . The MPO algorithm enables us to simulate relatively large chains of cavities (up to $20$ sites) but is unable to explore the region where long-range spatial correlations develop (convergence with respect to the bond dimension was not achieved in that region). The MPO approach has been therefore used extensively to investigate the nature of the strongly localized Mott-like phase in large arrays. The corner-space renormalization approach allows us to compute the global density matrix with arbitrary accuracy across the critical region (and therefore can be used to access information as the entropy of the steady state), but it is more limited by the system size (we report here results across the critical region for chains with size up to $8$ sites).

In Fig. 7 we show the average photon density $n=\sum_{i=1}^{M}\braket{\hat{a}^{\dagger}_{i}\hat{a}_{i}}/M$, its fluctuations $\Delta n^{2}=\sum_{i=1}^{M}\Delta n_{i}^{2}/M$, and the compressibility $\mathcal{K}=\sum_{i=1}^{M}\Delta n_{i}^{2}/(n_{i}M)$ for different system sizes ($M$) as a function of the hopping rate $zJ/\Gamma_{p}$ (here $z=2$) in the hard-core limit. In these numerical calculations, we have chosen a slightly shifted atomic frequency $\omega_{at}=\omega_{c}-zJ$. Under this condition, gain is strongest at the bottom of the photon band, so that condensation into the $\mathbf{k}=0$ mode is explicitly favored from the outset. This ensures that the spatially homogeneous condensation process that we have assumed in the mean-field calculation is not disturbed by mode competition phenomena between pairs of modes with opposite wavevector which display the same gain and could give rise to condensate fragmentation effects frag . This fragmentation mechanism was likely the reason of the reduced coherence numerically found in rivas2014 .

We observe that the data collapse on the single-cavity predictions ($J=0$) for $J\lesssim J_{c}^{\rm HC}$. This result is in agreement with the findings of the mean-field analysis which predicts a region where local interactions dominate over cooperative effects. Such Mott-like phase (highlighted by the logarithmic scale) is thus characterized by an almost integer local density and almost vanishing fluctuations and compressibility. This means that in this parameter range, the correlations among different cavities are very small and cooperative effects are suppressed. As the photon hopping rate is increased ($J\gtrsim J_{c}^{\rm HC}$), the system enters a novel regime with finite (significant) local density fluctuations and consequently a large value of the compressibility. Our numerical results in one-dimensional systems are consistent with the expectation that in 1D a true phase transition to a lasing regime, as predicted by the Gutzwiller theory, is replaced by a crossover from a Mott-like phase to a mixture of extended Tonks-Girardeau states with different number of photons girardeau1960 ; carusotto2009 . We remark that the hard-core regime forbids the multiple occupation of a lattice site but allows for the multiple (eventually macroscopic) occupation of a given $\mathbf{k}$-mode of the lattice.

The transition between these two phases is also signaled by a sharp variation of the entropy of the steady-state density matrix ($S={\rm Tr}[\rho^{SS}\ln(\rho^{SS})]$, lower right panel of Fig. 7). As predicted at mean-field level, in the Mott-like phases the system is an almost pure state while the symmetry broken phase has a strongly mixed character. This appears as a common feature of phase transitions in driven-dissipative lattices jin2016 ; rota2016 . In the inset, we show $S/M$, the entropy rescaled by the number of cavities. In the Mott-like phase, where the correlations among the cavities are very weak, the data referring to different system sizes collapse, indicating the extensive nature of the entropy. Indeed, when the steady-state density matrix is factorizable one gets $S=M\ S_{\rm sc}$, where $S_{\rm sc}$ is the single cavity entropy.

In Fig. 8, we show $\partial S/\partial(zJ/\Gamma_{p})$, the derivative of the entropy with respect to the hopping rate: it displays a peak that, for the considered range of sizes, becomes more pronounced and closer to the Gutzwiller mean-field critical coupling ($zJ_{c}^{\rm HC}/\Gamma_{p}=3.3$) as the number of sites is increased. However, it is expected that for one-dimensional systems in the thermodynamic limit no divergence of the entropy derivative occurs since the phase transition predicted by the mean-field theory should be replaced by a crossover between the two phases wouters2006 ; szy2006 .

To better understand the build up of the quantum correlations among the different cavities and to further characterize the Mott-like phase we also computed the one-body correlation function in the steady-state

[TABLE]

To do so we exploited an MPO algorithm which allows to dynamically reach the steady-state for large systems in the Mott-like regime. In Fig. 9 we show the behavior of $c(i,i+r)$ for $M=20$ as $zJ/\Gamma_{p}$ is varied. Since $c(i,j)$ has an oscillatory behaviour due to the very large nonlinearity, we considered its absolute value. The one-body correlation function decays exponentially with the distance

[TABLE]

where the correlation length $\lambda$ increases as the photon hopping increases. In particular, it is possible to show that $\lambda$ scales as

[TABLE]

in the $J/\Gamma_{p}\ll 1$ limit. The details about the derivation of Eq. (22) can be found in App. C.

IV Conclusions

In this work, we have investigated the steady-state phases of a photonic lattice in presence of incoherent driving, dissipation and strong photon-photon interactions. We have explored a general model where each lattice site (a nonlinear Kerr resonator) is coupled to a two-level emitter, which is pumped incoherently. Via a Gutzwiller decoupling theory, we have determined the non-equilibrium phase diagram. We have found that the interplay between on-site interactions, photon hopping and driven-dissipative processes, lead to a second-order nonequilibrium phase transition related to the spontaneous breaking of the $U(1)$ symmetry possessed by the model. Furthermore, we have shown that such incoherent driving scheme allows to stabilize Mott-like phases of light characterized by an almost integer local density and almost vanishing compressibility. The picture predicted by the Gutzwiller mean-field has been validated by numerical finite-size simulations using matrix product operators and the corner-space renormalization method. By driving the system across the critical point, we have characterized the transition both in terms of one-body correlation functions (which display long-range order) and of the entropy of the system (which increases in the symmetry broken phase). The phase transition from an (almost) incompressible Mott-like photon fluid to a coherent delocalized phase is driven by commensurability effects. The control parameter of this transition can be thus deduced by comparing the effective emission rate at the band boundary with the photon leakage rate. The mixed character of the steady-state density matrix reflects the intrinsic nonequilibrium nature of the coherent phase. Remarkably, signatures of the phase-transition are present already in small arrays.

The strongly correlated photon phases proposed here could be explored by using photonic quantum simulators based on circuit QED lattices houck2012 . In particular, in these systems thanks to superconducting quantum resonators and Josephson junctions it is possible to engineer large Kerr photon-photon interactions and to tailor the interaction with two-level emitters wallraff2004 ; fink2008 , paving the way to the study of strongly correlated quantum fluids of light and non-equilibrium phase transitions. The incoherent pump scheme we exploit can be implemented by coherently driving the emitter into a third metastable level from which it fast decays into the excited state of the active transition, thus resulting into an effective incoherent pump ma2017 . The role of dimensionality, disorder and criticality in these nonequilibrium quantum phases are intriguing topics that need to be explored experimentally and theoretically in the future.

V Acknoledgments

Discussions with Leonardo Mazza are warmly acknowledged. AB, FS and CC acknowledge support from ERC (via Consolidator Grant CORPHO No. 616233). AB, CC, DR, FS, IC and RF acknowledge the Kavli Institute for Theoretical Physics, University of California, Santa Barbara (USA) for the hospitality and support during the early stage of this work. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. JL and IC are supported by the EU-FET Proactive grant AQuS, Project No. 640800, and by the Autonomous Province of Trento, partially through the project “On silicon chip quantum optics for quantum computing and secure communications” (“SiQuro”). RF acknowledges Oxford Martin School Technologies and the Singapore Ministry of Education and National Research Foundation (QSYNC) for financial support.

Appendix A Single-cavity solution (hard-core limit)

When the nonlinearity is large enough to forbid the double occupation of a lattice site, one can map the bosonic degree of freedom into an effective spin

[TABLE]

Doing so, the single-cavity Hamiltonian becomes

[TABLE]

In the Heisenberg representation the evolution of a given operator $\hat{\theta}$ is ruled by

[TABLE]

where

[TABLE]

with $\mathcal{\tilde{D}}[\hat{O};\rho]=[2\hat{O}^{\dagger}\rho\hat{O}-\{\hat{O}^{\dagger}\hat{O},\rho\}]$. Writing the equation of motion for $\braket{\hat{\Sigma}^{\alpha}}$, $\braket{\hat{\sigma}^{\alpha}}$, $\braket{\hat{\Sigma}^{\alpha}\hat{\sigma}^{\beta}}$ ($15$ equations) and solving for the stationary state one always find a single (stable) solution given by

[TABLE]

where $\Delta\omega=\omega_{at}-\omega_{c}$, $\Gamma_{tot}=\Gamma_{p}+\gamma+\Gamma_{l}$ and $\Lambda=\Gamma_{tot}^{2}+4\Delta\omega^{2}$.

Appendix B Mean-field stability analysis (hard-core limit)

At mean-field level the dynamics in the hard-core regime is ruled by the master equation (27) where the single-cavity Hamiltonian (26) is replaced by

[TABLE]

The single-cavity steady-state solution (29) is always a fixed point of the single-site MF equations with $\braket{\hat{a}}=0$ ($\braket{\hat{\Sigma}^{x}}=\braket{\hat{\Sigma}^{y}}=0$) which becomes unstable for certain values of the parameters and allows the emergence of a phase characterized by $\braket{\hat{a}}\neq 0$. In order to study the stability of this solution we computed the Jacobian matrix of the MF equations leboite2013 ; lee2013 ; wilson2016 evaluated with respect to the single-cavity solution (29). If the real part of one of its eigenvalues becomes positive for certain values of the parameters then the single-cavity fixed point is unstable. In Fig. 5 we show the Heaviside step function of the real part of the most unstable eigenvalue of the Jacobian $\Theta\left[{\rm Re}\left(\lambda_{u}\right)\right]$ as a function of the typical values of the system parameters. The result is a plot which highlights the stable ($\braket{\hat{a}}=0$) and the unstable ($\braket{\hat{a}}\neq 0$) regions in yellow and light blue respectively.

Appendix C Spatial correlation functions (1D)

In this appendix we provide a scaling law for the spatial decay of the one-body correlations in the $J/\Gamma_{p}\ll 1$ regime by using an ansatz for the photonic steady-state density matrix in the hard-core limit. Even if the derivation might at first appear somewhat heuristic, it will provide a good approximation of how the correlation length depends on the system parameters (see Sec. III.3). In analogy with the equilibrium physics of hard-core bosons in one dimension, we suppose that the steady-state is fermionized, i.e., that the photonic density matrix $\rho$ is diagonal in the fermionic momentum basis, up to a unitary Jordan-Wigner transformation $\hat{U}$ which anti-symmetrizes the bosonic density matrix:

[TABLE]

The bosonic and fermionic annihilation operators are related through the unitary relation

[TABLE]

where the local particle number operator $\hat{n}_{l}=\hat{n}^{B/F}_{l}$ is left unchanged by the anti-symmetrization process.

In the simple case of free bosons, the steady state momentum distribution can be exactly calculated analytically:

[TABLE]

where $\tilde{\delta}=-2\Delta\omega/\Gamma_{p}$ and $\epsilon_{k}=-2J\cos(k)$. So that, a natural ansatz for the Fermi non-equilibrium distribution would be

[TABLE]

The analytical function

[TABLE]

of the complex variable $z$ can be expanded as

[TABLE]

where $r_{c}$ is the convergence radius of this power serie, and $\alpha_{n}$ is a sub-exponential function of the variable $n$ which is dominated by any geometric function. The convergence radius is given by the modulus of the complex pole of $g(z)$ which is closest to the origin $z=0$, and thus has the following expression $r_{c}=\sqrt{1+\tilde{\delta}^{2}+\Gamma_{em}^{0}/\Gamma_{l}}$.

From the ansatz Eq. (53), it is possible to calculate the long-range properties of the fermionic one-body correlation function

[TABLE]

Exploiting the expansion Eq. (55) we obtain

[TABLE]

as the integral in the right side is non-zero only for $n\geq j$. In the limit, $2J/(r_{c}\Gamma_{p})\ll 1$ we can keep only the lowest power of $J$ in the serie expansion, which corresponds to $n=j$

[TABLE]

where $\alpha_{j}$ is sub-exponential. The fermionic autocorrelation thus spatially decays exponentially with a correlation length $\lambda_{F}$ given by

[TABLE]

and scales as $1/\text{ln}(\Gamma_{p}/J)$ for small $J/\Gamma_{p}$.

The true photonic correlation function is related to the fermionic one as follows

[TABLE]

For a vanishing $J$, the fermionic distribution is nearly momentum independent, and the particles are thus fully localized with the uniform spatial density $n=1/(\Gamma_{l}/\Gamma_{em}^{0}+1)$. The correlation function can be thus factorized as $c(j)\simeq\braket{\hat{a}^{F}_{j}\hat{a}_{0}^{F\dagger}}\prod_{0<l<k}\braket{e^{i\pi\hat{n}_{l}}}$. The local expectation value $\braket{e^{i\pi\hat{n}_{l}}}\simeq(1-\Gamma_{em}^{0}/\Gamma_{l})/(1+\Gamma_{em}^{0}/\Gamma_{l})$, is positive for hole-dominated statistics ($\Gamma_{em}^{0}/\Gamma_{l}<1$) and negative for particle dominated statistics ($\Gamma_{em}^{0}/\Gamma_{l}>1$).

The expression for the photonic one-body correlation function in the $J/\Gamma_{p}\ll 1$ limit reads as

[TABLE]

for particle-dominated statistics. The corresponding correlation length is given by

[TABLE]

It maintains the logarithmic scaling, and is slightly shorter than the fermionic one, due to a scrambling induced by sign changes when crossing intermediary particles.