Phase diagram of a QED-cavity array coupled via a N-type level scheme

We study the zero-temperature phase diagram of a one-dimensional array of QED cavities where, besides the single-photon hopping, an additional coupling between neighboring cavities is mediated by an N-type four-level system. By varying the relative strength of the various couplings, the array is shown to exhibit a variety of quantum phases including a polaritonic Mott insulator, a density-wave and a superfluid phase. Our results have been obtained by means of numerical density-matrix renormalization group calculations. The phase diagram was obtained by analyzing the energy gaps for the polaritons, as well as through a study of two-point correlation functions.


Introduction
The recent impressive advances in the field of quantum simulators allowed to probe the many-body physics of strongly correlated systems at the level of the single quantum object. At present cold atoms trapped in optical lattices can be considered among the most promising examples of quantum simulators. By means of ultracold atomic and molecular gases, it is nowadays possible to reach a degree of control and accuracy in engineering the dynamics of many-body systems that were unimaginable in previously. As a consequence, the coherent quantum dynamics emerging from carefully tailored microscopic Hamiltonians can now be tested experimentally [1]. It has been possible, just to recall one example, to implement the Bose-Hubbard (BH) model [2,3] and to detect its zero-temperature superfluid (SF) to Mott insulator (MI) quantum phase transition [4]. Other models involving spinor gases, Fermi systems, Bose-Fermi mixtures, or dipolar gases have been also devised and realized, providing an even richer phase diagram (see for example the review [5]). We mention the stabilization of density-wave (DW) phases for bosons, as well as more peculiar topological or supersolid orderings, which can arise in the presence of finite-range interactions [6].
More recently a novel kind of many-body quantum simulator has been introduced, based on the idea to use single photons as quantum objects. Since photons hardly interact in open space, the most natural way to significantly increase their interactions is to trap them into an optical QED cavity, and couple the field with atoms/molecules inside it in order to create an optical nonlinearity. If the nonlinearity is sufficiently large, the so called photon blockade sets in [7,8], namely, the presence of a single photon inside a cavity prevents a second one to enter it. In the rotating-wave approximation, the simplest light-matter interaction scheme of this type can be accurately described by the Jaynes-Cummings (JC) model. By arranging an array of cavities coupled through the photon hopping, such to generate a competition between the hopping and the on-site nonlinearities, one can devise a setup that is well described by the so called Jaynes-Cummings-Hubbard (JCH) model [9][10][11].
In many respects, if one ignores dissipation, the physics emerging from the JCH Hamiltonian resembles, at low-energies, that of an effective BH model. Probably the main difference between the two systems is that, instead of having neutral bosons as building blocks of the model, in the JCH Hamiltonian one has to think in terms of polaritons, i.e., combined photonic/atomic excitations. Many different works already addressed the JCH equilibrium phase diagram with analytical, as well as numerical methods, leading to a fairly complete theoretical understanding of the nature and the location of the emerging phases and phase transitions in terms of the parameters governing the system (the field has been recently reviewed in, e.g., Refs. [12][13][14][15]).
Additional interest in cavity arrays comes from the fact that these systems can be naturally considered as open-system quantum simulators. Some related features have been recently explored [16][17][18][19][20][21][22][23][24]. In the following we will not touch on this and consider only the "equilibrium" phase diagram.
This intense activity has been very recently boosted by the first experiments on QED cavity arrays [25][26][27]. As of today, the most concrete possibility to realize controllable and scalable quantum simulators with cavity arrays involves circuit-QED cavities [28][29][30].
So far the coupling between cavities has been mostly considered through photon hopping. Only few works started addressing more general schemes, where the cavity coupling can be induced also by means of non-linear elements [23,24,31,32]. Such configurations include cross-Kerr interactions and/or correlated hopping terms, which lead to generalizations of the JCH model in a way similar to the extended BH (EBH) Hamiltonian for atoms with large dipole momentum loaded in optical lattices [33]. The underlying physical model is believed to possess a much richer structure, with the emergence of exotic phases of correlated polaritons. It is particularly interesting to address these schemes in one-dimensional (1D) systems, where interactions become crucial to stabilize exotic phases of matter [33][34][35][36][37]. These notably include a series of nontrivial density-wave (DW) states, which can arise in the strong coupling regime [38], as well as supersolidity and phaseseparation effects [39,40]. Extension to consider also counter-rotating terms in the ultrastrong coupling regime, thus leading to the so called Rabi-Hubbard model [41], have been investigated [42]. However we are not aware of numerical investigations of coupled cavity models beyond the JCH and Rabi-Hubbard model.
In all such situations, non-perturbative, either numerical or analytical calculations are necessary. Here the density-matrix renormalization group (DMRG) algorithm [43,44] has been employed to work out the quantitative zero-temperature phase diagram of the JCH model [45][46][47]. This is a particularly efficient method for the statics of 1D many-body problems. Its key strategy consists in constructing a portion of the system (called block) and then recursively enlarge it. At each step, the basis of the corresponding Hamiltonian is truncated to a given value m, so that one can manage the Hamiltonian in an effective Hilbert space of fixed dimensions, as the physical system grows. This truncation is performed by retaining the eigenstates corresponding to the m highest eigenvalues of the reduced density matrix of the block.
The aim of this paper is to quantitatively study a generalization of the JCH Hamiltonian, aimed at taking into account an effective nearest-neighbor nonlinearity between cavities mediated by an N-type four-level system as discussed for two cavities in Ref. [48]. The presence of this coupling leads to an effective cross-Kerr non-linearity. An analysis at the mean-field level of a dissipative open EBH as an effective model for nonlinearly coupled cavities has been performed, unveiling the emergence of novel photon crystal and supersolid phases [23,24]. Here we do not resort to the effective EBH model and analyze the full model as introduced in [48]. Using the DMRG algorithm, we work out the 1D ground-state phase diagram. We show that a physics similar to the EBH model appears, with a rich phase diagram including gapless SF, as well as MI and DW phases of polaritons. We postpone the analysis of the interplay of driving and dissipation to a future work.
The paper is organized as follows. In the next two sections we introduce the model of coupled cavities of our interest (Sec. 2) and the quantities we are going to address, namely the energy gaps, and the staggered number-number correlations (Sec. 3). In Sec. 4 we discuss the zero-temperature equilibrium phase diagram, focusing on the MI/SF boundary and on the boundary separating the DW from the other phases. Finally, in Sec. 5 we draw our conclusions.

The model
Let us consider a 1D array of QED cavities, where photons can hop between neighboring cavities. Moreover two adjacent resonators are also nonlinearly coupled to each other via a N-type four-level system, as shown in Fig. 1(a). For the sake of clarity in our description, we shall divide the 1D array into coupled effective sites composed of a cavity and an atom. The four levels are denoted by {|i } i=1...4 , and are depicted in Fig. 1 Each effective site is composed of a cavity and an atom (dashed box). (b) Level structure of the N-type atoms. The transition |1 ↔ |3 is resonantly coupled to the cavity mode of its own site with strength g 1 , while |2 ↔ |4 is coupled to the cavity mode of its right nearest-neighbor site with strength g 2 , and has a detuning ∆. The transition |2 ↔ |3 is in resonance with an external laser field of strength Ω.

(b). An external laser
with frequency Ω resonantly drives the transition |3 ↔ |2 . The transition |1 ↔ |3 is resonantly coupled to the cavity mode of the same site with strength g 1 , while the transition |2 ↔ |4 couples to the cavity mode of its right nearest-neighbor site with strength g 2 , and a detuning ∆.
The use of such N-type atom for generating large Kerr nonlinearity has been extensively studied in the literature [7,8,49], however the vast majority of the scenarios only focused on a single-mode cavity. Our work is inspired by the idea of Ref. [48], where the cross-Kerr nonlinearity is generated between two different and neighboring cavities, in circuit-QED systems. In practice, we use the unbalanced couplings of atomic transition |1 ↔ |3 with left cavity mode, and |2 ↔ |4 with right cavity mode respectively, in order to generate the local (g 1 ) and nonlocal (g 2 ) nonlinearities of our many-body system. This kind of fourlevel artificial molecule can be realized using two Josephson transmon qubits coupled by a superconducting quantum interference device.
Using the interaction picture and in the rotating-wave approximation, the system Hamiltonian reads where σ mn = |m n| , (m, n = 1, 2, 3, 4), and a (a † ) is the annihilation (creation) operator of the cavity mode. The subscripts denote the site position along the 1D chain. The first three terms in the r.h.s. of Eq. (1) describe the local terms and the nonlinearities on each site. Inside the latter brackets, the first term is the photon hopping, while the second term describes the coupling of the atom to its right neighboring cavity, which generates an effective nonlocal cross-Kerr nonlinearity between the two cavities.
Hereafter we concentrate on the 1D model in Eq. (1) at zero temperature, specifically addressing the case without dissipation with DMRG. Let us also fix the Hamiltonian quantities in units of Ω, set = 1, and work with open boundary conditions. We recall that, in the presence of dissipation, the problem becomes much more difficult to be handled numerically [1] .
For the system we are considering here, in the strong coupling regime atoms and photons cannot be considered as two separate entities. It is thus natural to investigate the phase diagram in terms of combined atomic/photonic modes, named polaritons. The polaritonic number operator on each site i, representing the number of local excitations, is defined as n pol i = 2σ 44 i + σ 33 i + σ 22 i + a † i a i . For the closed system described by the Hamiltonian (1), the total number N pol = i n pol i of such polaritons is a conserved quantity. In the following we work in the canonical ensemble for polaritons, and focus on the integer filling situation. [1] It is however possible to address the effect of dissipation with a DMRG approach in a 1D chain, when this is described by a master equation within the Lindblad formalism. In the language of tensor networks, one has to generalize the matrix-product-state ansatz to a matrix-product-density-operator ansatz for mixed states, as originally proposed in Refs. [50,51]. The computational complexity is greater than for static computations, and is eventually related to the amount of entanglement in the steady state.

Energy gaps and correlation functions
The different nature of the various phases is sensitive to a number of properties which we are going to focus on. Here we are going to study quantities that resemble those characterizing the various phases of the EBH model [36].
First of all, the ground-state energy gap is an important indicator which characterizes the presence or absence of criticality in the model. In particular, in the critical SF phase, the charge gap vanishes in the thermodynamic limit. On the other side, in the insulating MI and DW phases, such gap remains finite. In order to make connection with a similar notation in the EBH model, below we introduce the so called charge and neutral gaps referring respectively to the gaps corresponding to adding one extra particle ("charge" sector) or remaining with the same number of particles ("neutral" sector). We stress however that in the present model the excitation carry no real charge. This has to be understood only as a convention.
The charge gap is defined as where, in the canonical ensemble, ∆E + (∆E − ) denotes the extra energy needed to add (remove) one particle, i.e. one polariton, in the system. In the specific, focusing on the unit filling, where E L is the ground-state energy per site of an L-sites cavity-array with exactly L excitations, and E L+1 (E L−1 ) is the corresponding energy per site with one excitation more (less). It is therefore possible to extrapolate ∆E c by running three different DMRG simulations with fixed number of polaritons N pol = L − 1, L, L + 1 [52,53]. While the charge gap is able to detect particle-hole excitations, in some circumstances it is possible that the dominant low-energy excitations are of a different type. Their presence can be detected only by the so called neutral gap at a fixed number of particles, where, again working in the canonical ensemble, E 1 L denotes the first excited energy per site of an L-site system with L excitations.
In the following, we also focus on the analysis of the staggered diagonal order for the polaritons, in order to distinguish the DW from the other phases. We do this by investigating the two-point correlation function where δn pol i = n pol i −n denotes the polariton fluctuation from the average fillingn. The order parameter identifying the DW phase is thus given by: O DW ≡ lim r→∞ C DW (r). A finite value of O DW indicates a tendency to establish, in the thermodynamic limit, a staggered occupation of the polaritons. On the other side, in the MI as well as the SF phases, C DW (r) vanishes exponentially with increasing distance r.

Phase diagram
The zero-temperature phase diagram of model (1), at unit polariton fillingn = 1 and in the g 2 − t plane, is summarized in Fig. 2. We observe that three different phases can be stabilized. Their boundaries have been obtained by means of a finite-size scaling of the numerical data, for systems up to L = 300 sites. In our simulations we imposed a cutoff photon number in each cavity, such that n phot i ≤ 3. We also truncated the effective Hilbert space dimension to a value m = 80 in all the simulations, except for those shown in Fig. 6 for the neutral energy gap (see the discussion in Sec. 4.3). We checked that, by increasing m and the local fock-space truncation over the photon number, the results concerning the charge gap and the DW order parameter do not change on the scales shown here.
For small photon hopping (t/Ω 0.2), by increasing the nonlocal nonlinearity g 2 the system exhibits a direct transition from the MI to the DW phase. On the other hand, for t/Ω 0.2, the MI-to-DW transition is mediated by an extended region appearing at intermediate g 2 values, where the system stabilizes into a gapless SF. In the following we are going to elucidate our finite-size scaling procedure and how we were able to distinguish between the different phases.

Boundary between MI and SF phases
In the limit of small g 2 and t values, the dominant presence of the on-site interactions stabilize the system into a MI phase with exactly one polariton per cavity (n = 1), and where the charge energy gap has a finite value. As long as the hopping strength t is progressively increased (and for fixed g 1 , g 2 ), the system eventually enters a SF phase, with a vanishing gap. The filled squares of Fig. 2 denoting the MI/SF boundaries have been obtained by means of a finite-size scaling of the charge gap. We performed simulations up to L = 100 sites and analyzed whether the gap closes or remains finite in the thermodynamic limit L → ∞.
In Fig. 3, left panel, we highlight the size-dependence of ∆E c as a function of 1/L for two points in the phase space close to the MI/SF transition (see points along the dotted line in Fig. 2). We expect to see a quadratic dependence ∆E c ∼ L −2 (dashed line) at large L [52,53], however a linear extrapolation (solid line) is already a good approximation to the scaling, and we can use it to determine ∆E c in the thermodynamics limit. Indeed, we observe that the difference between quadratic and linear extrapolation is tiny ( 10 −3 ) and does not produce any distinguishable modification on the scale of Fig. 2. In the specific case of Fig. 3, we fixed t/Ω = 0.25 and chose two different values of g 2 /Ω corresponding to configurations in the gapped MI (g 2 /Ω = 1.35, triangles) and in the gapless SF phase (g 2 /Ω = 1.5, squares). The MI is signaled by an extrapolated finite value of lim L→∞ ∆E c > 0, while in the SF this is zero.
In order to locate the critical g 2 for a given value of t (filled squares in Fig. 2) we perform a linear extrapolation of the charge gaps in the vicinity of the critical value of g 2 . An example of such procedure is shown in the right panel of Fig. 3, where we plot ∆E c as a function of g 2 , when this is close to the phase transition (in the specific, here we set t/Ω = 0.25). After a linear extrapolation, we get a critical g 2 value corresponding to g * 2 /Ω ≈ 1.379. An analogous procedure is repeated for all the filled squares shown in Fig. 2, thus identifying the MI/SF boundary.

Boundary of the DW phase
The DW phase is characterized by a finite order parameter O DW . Let us therefore look at the two-point staggered correlator in Eq. (4). Since in DMRG simulations we are employing open boundary conditions, to minimize the border effects we analyze the correlations in such a way that the two points are taken symmetrically with respect to the center of the system [2] . The left panel of Fig. 4 shows how differently such polariton correlations behave when the system goes from the MI to DW phase, for a fixed system size.
To be more accurate, in the right panel we performed a finite-size scaling and showed that the staggered correlation C DW (r) approaches the zero value exponentially with L, in the MI phase (a similar behavior occurs in the SF region). On the other hand, in the DW such correlator asymptotically converges to a finite value. In the specific, here we fix t/Ω = 0.05 and show that for g 2 /Ω = 1.3, 1.35 the DW order is exponentially suppressed with L, while for g 2 /Ω = 1.4 it remains finite. The O DW order parameter reached for L → ∞ is displayed in the inset as a function of g 2 .
In order to determine the DW boundary in the phase diagram of Fig. 2, we adopted the following protocol. For a fixed value of t/Ω, we start increasing g 2 from zero up to a finite value, with a fixed increment δg 2 = 0.05Ω, and to compute the DW order parameter for all such values of g 2 . The boundary of DW phase in the g 2 −t plane (filled circles in Fig. 2), for any fixed t, is located by the g * 2 (t) that gives the first non-vanishing order parameter O DW . Here we stress that, because of the arrangement of our 1D array [see Fig. 1(a)] and of the asymmetric coupling between the atom and its right/left cavity, the antiferromagnetic symmetry of the system is spontaneously broken. In particular, the state |4 of the L-th atom will be never occupied, since the transition |2 ↔ |4 does not couple to any cavity mode [see Eq. (1)]. As a consequence, in our simulations we do not need any symmetry-breaking potential. We can observe that the expectation value for the onsite number of polaritons explicitly exhibits a staggered behavior, in that the occupation of the (2n − 1)-th site is always higher than that of the (2n)-th site (for any integer value of n). Finally we notice that such staggering persists at finite size, also for the set of parameters corresponding to the MI phase, although it is extremely tiny and decreases with L. This effect eventually disappears in the thermodynamic limit. [2] The two points of δn pol i δn pol j , with |i− j| = r, have been chosen such that i = (L−r+1)/2, j = (L+r+1)/2 for odd r, and i = (L − r)/2, j = (L + r)/2 for even r (e.g. for L = 100 sites, r = 1 corresponds to i = 50, j = 51; r = 2 corresponds to i = 49, j = 51; r = 3 to i = 49, j = 52, and so on) For a positive detuning, the system will never present as a DW phase, meaning that C DW (L/2) = 0. Taking a negative value of the detuning, we observe that C DW (L/2) increases with |∆|. The parameters in this figure are t/Ω = 0.05, g 1 /Ω = 0.8, g 2 /Ω = 1.6 and the system size is L = 100.
The extension of the DW phase depends on the cavity detuning ∆. In particular, the robustness of the order parameter increases with increasing the modulus of the detuning (see Fig. 5). Quite remarkably, we note that a positive ∆ will never stabilize an antiferromagnetic DW ordering.

Neutral gap
The analysis leading to the phase diagram in Fig. 2 has been corroborated by a study of the neutral gap, which vanishes both in proximity of the phase transitions and in the entire superfluid region. Differently for the charge gap, it is able to detect the presence of excitations other than particle-hole, and thus locates the boundaries of insulating regions (as the DW) beyond the MI.
The data displayed in Fig. 6 show the behavior of ∆E n as a function of g 2 , for a fixed value of t/Ω. In particular we analyzed a vertical cut in the phase diagram of Fig. 2 (see the rightmost vertical dotted line in that figure), where the system can be in three different phases according to the value of g 2 . With increasing g 2 , it goes from the MI phase (nonzero ∆E n , for t/Ω 1.45) to the SF phase (zero ∆E n , for 1.45 t/Ω 1.8), and then to the DW phase (nonzero ∆E n , for t/Ω 1.8).  While we cannot see a clear signature of a finite gap for g 2 = 1.8Ω, the scaling with the size displayed in the left panel of Fig. 6 seems to suggest the scenario depicted above. It is however important to stress that the DMRG simulations needed to compute the neutral gap have to target the two lowest-lying eigenstates in a single run. Thus they generally require a larger dimension m of the effective Hilbert space, as compared to all the other ground-state calculations discussed before. The analysis of the neutral gap requires a careful convergence test of the results with m, which we provide in the right panel of Fig. 6. We observe that the non monotonic features that are visible in the region 1.45 t/Ω 1.8 have to be probably ascribed to the inaccuracy of the method at small m values. This signals the presence of the gapless SF phase there, in agreement with the results provided by the charge gap (MI/SF boundary) and for the DW order parameter (SF/DW boundary).

Summary
Using the density-matrix renormalization group with open boundary conditions, we studied the equilibrium phase diagram of a system of coupled QED cavities in one dimension. We provided results beyond the standard model of couplings through photon hopping, and also considered nearest-neighbor cross-Kerr nonlinearities. Our analysis is based on a finite-size scaling of the ground-state charge and neutral gaps, as well as of the density-wave order parameter, for systems up to 300 sites. We showed that, beyond the conventional Mott insulator and superfluid phases, the presence of a nearest-neighbor nonlinear coupling can also stabilize a density-wave ordering of polaritons.