Skip to main content



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.


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 [911].

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. [1215]).

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 [1624]. 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 [2527]. As of today, the most concrete possibility to realize controllable and scalable quantum simulators with cavity arrays involves circuit-QED cavities [2830].

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 [3337]. 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 phase-separation 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 [4547]. 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 nonlinearity. 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 (Section 2) and the quantities we are going to address, namely the energy gaps, and the staggered number-number correlations (Section 3). In Section 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 Section 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 Figure 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 \(\{ \vert i \rangle \}_{i = 1 ,\ldots,4}\), and are depicted in Figure 1(b). An external laser with frequency Ω resonantly drives the transition \(\vert 3 \rangle \leftrightarrow \vert 2 \rangle\). The transition \(\vert 1 \rangle \leftrightarrow \vert 3 \rangle\) is resonantly coupled to the cavity mode of the same site with strength \(g_{1}\), while the transition \(\vert 2 \rangle \leftrightarrow \vert 4 \rangle\) couples to the cavity mode of its right nearest-neighbor site with strength \(g_{2}\), and a detuning Δ.

Figure 1

Scheme for nonlinearly coupled QED cavities. (a) An array of QED cavities nonlinearly coupled by N-type atoms. The photon hopping between nearest-neighbor cavities has a strength t. Each effective site is composed of a cavity and an atom (dashed box). (b) Level structure of the N-type atoms. The transition \(\vert 1 \rangle \leftrightarrow \vert 3 \rangle \) is resonantly coupled to the cavity mode of its own site with strength \(g_{1}\), while \(\vert 2 \rangle \leftrightarrow \vert 4 \rangle\) is coupled to the cavity mode of its right nearest-neighbor site with strength \(g_{2}\), and has a detuning Δ. The transition \(\vert 2 \rangle \leftrightarrow \vert 3 \rangle \) is in resonance with an external laser field of strength Ω.

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\rangle\leftrightarrow|3\rangle\) with left cavity mode, and \(|2\rangle\leftrightarrow|4\rangle\) 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 four-level 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

$$ \mathcal{H}= \sum_{i} \bigl[ \Delta\sigma^{44}_{i} + \bigl( \Omega\sigma^{23}_{i} + g_{1} \sigma_{i}^{13} a_{i}^{\dagger}+ \mathrm{H.c.} \bigr) + \bigl(-t a_{i} a_{i+1}^{\dagger}+ g_{2} \sigma_{i}^{24} a_{i+1}^{\dagger}+ \mathrm{H.c.} \bigr) \bigr] , $$

where \(\sigma^{mn} = \vert m \rangle \langle n\vert \) (\(m,n = 1,2,3,4\)), and a (\(a^{\dagger}\)) 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 \(\hbar= 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.Footnote 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^{\mathrm{pol}}_{i}=2\sigma_{i}^{44}+\sigma^{33}_{i}+\sigma ^{22}_{i}+a^{\dagger}_{i}a_{i}\). For the closed system described by the Hamiltonian (1), the total number \(N^{\mathrm{pol}} = \sum_{i} n^{\mathrm{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.

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

$$ \Delta E_{c} = \Delta E^{+} - \Delta E^{-} , $$

where, in the canonical ensemble, \(\Delta E^{+}\) (\(\Delta 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, \(\Delta E^{+} = E_{L+1}-E_{L}\) and \(\Delta E^{-} = E_{L}-E_{L-1}\), 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 \(\Delta E_{c}\) by running three different DMRG simulations with fixed number of polaritons \(N^{\mathrm{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,

$$ \Delta E_{n} = E_{L}^{1} - E_{L} , $$

where, again working in the canonical ensemble, \(E_{L}^{1}\) 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

$$ C_{\mathrm{DW}}(r)=(-1)^{r} \bigl\langle \delta n_{i}^{\mathrm{pol}} \delta n_{i+r}^{\mathrm{pol}} \bigr\rangle , $$

where \(\delta n_{i}^{\mathrm{pol}} = n^{\mathrm{pol}}_{i}-\bar{n}\) denotes the polariton fluctuation from the average filling \(\bar{n}\). The order parameter identifying the DW phase is thus given by: \(\mathcal{O}_{\mathrm{DW}}\equiv\lim_{r\rightarrow\infty }{C_{\mathrm{DW}}(r)}\). A finite value of \(\mathcal{O}_{\mathrm{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_{\mathrm{DW}}(r)\) vanishes exponentially with increasing distance r.

Phase diagram

The zero-temperature phase diagram of model (1), at unit polariton filling \(\bar{n} = 1\) and in the \(g_{2}-t\) plane, is summarized in Figure 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^{\mathrm{phot}}_{i} \leq3\). We also truncated the effective Hilbert space dimension to a value \(m=80\) in all the simulations, except for those shown in Figure 6 for the neutral energy gap (see the discussion in Section 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.

Figure 2

Phase diagram. Ground-state phase diagram for a 1D system of coupled cavities described by Hamiltonian (1) in the \(g_{2}-t\) parameter space. Here and in the subsequent figures we choose \(g_{1} / \Omega= 0.8\), \(\Delta/ \Omega= -2\). The various symbols correspond to points belonging to different phases (triangles = MI, circles = DW, squares = SF). Filled points lie close to a phase transition, and are used to draw the interpolating curves (dashed lines). The DW-to-MI/SF boundary (filled circles) has been obtained by analyzing the density-wave order parameter, while the MI-SF boundary (filled squares) through the charge gap. The two vertical dotted lines denoted two cuts in the phase diagram which will be analyzed in details below. The error bars in the points characterizing the DW-to-MI/SF boundary take into account the discretization of the \(g_{2}\) values that we adopted in our numerical simulations (see the text).

For small photon hopping (\(t / \Omega\lessapprox0.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 / \Omega\gtrapprox0.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 (\(\bar{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 Figure 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 \to\infty\).

In Figure 3, left panel, we highlight the size-dependence of \(\Delta 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 Figure 2). We expect to see a quadratic dependence \(\Delta E_{c} \sim 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 \(\Delta 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 Figure 2. In the specific case of Figure 3, we fixed \(t / \Omega= 0.25\) and chose two different values of \(g_{2} / \Omega\) corresponding to configurations in the gapped MI (\(g_{2} / \Omega= 1.35\), triangles) and in the gapless SF phase (\(g_{2} / \Omega= 1.5\), squares). The MI is signaled by an extrapolated finite value of \(\lim_{L \to \infty} \Delta E_{c} > 0\), while in the SF this is zero.

Figure 3

Analysis of the MI-SF boundary. Left panel: system-size dependence of the charge gap \(\Delta E_{c}\) per site, for fixed \(t/\Omega= 0.25\) and two values of \(g_{2}\) in the MI (\(g_{2}/\Omega= 1.35\)) and in the SF (\(g_{2}/\Omega= 1.5\)) phase. Symbols denote the DMRG results. Solid and dashed lines are linear and quadratic fitting curves, respectively. The difference between the extrapolated values in the two fits \(\Delta E_{c}^{\infty} = \lim_{L \to\infty} \Delta E_{c}^{(L)}\) is negligible. Right panel: determination of the critical \(g_{2}^{*}\) value for the quantum phase transition. The triangles denote the charge gap per site \(\Delta E_{c}^{\infty}\) at the thermodynamic limit, as extrapolated in the left panel. The dashed line is a best linear fit of the data vs. \(g_{2}\). The critical point is obtained when \(\Delta E_{c}\) vanishes. For \(t / \Omega=0.25\), we get \(g_{2}^{*} / \Omega\approx1.379\).

In order to locate the critical \(g_{2}\) for a given value of t (filled squares in Figure 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 Figure 3, where we plot \(\Delta E_{c}\) as a function of \(g_{2}\), when this is close to the phase transition (in the specific, here we set \(t / \Omega= 0.25\)). After a linear extrapolation, we get a critical \(g_{2}\) value corresponding to \(g_{2}^{*} / \Omega\approx1.379\). An analogous procedure is repeated for all the filled squares shown in Figure 2, thus identifying the MI/SF boundary.

Boundary of the DW phase

The DW phase is characterized by a finite order parameter \({\mathcal{O}}_{\mathrm{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.Footnote 2 The left panel of Figure 4 shows how differently such polariton correlations behave when the system goes from the MI to DW phase, for a fixed system size.

Figure 4

Determination of the DW boundaries. The two-point correlation function \(C_{\mathrm{DW}}(r)\) for the polariton number and its asymptotic value near the MI-DW quantum phase transition. Here we fix \(t/\Omega= 0.05\) and vary \(g_{2}/\Omega\). Left panel: behavior at fixed system size \(L=300\), as a function of the distance r and for different values of \(g_{2}/\Omega= 1.3\), 1.35, 1.4, and 1.6. To minimize boundary effects, we chose the two points \((i, i+r)\) symmetrically with respect to the center of the array. Right panel: finite-size scaling close to the transition. Empty circles, filled circles, and triangles respectively are for \(g_{2} / \Omega= 1.4\) (DW), 1.35 (near the critical point), and 1.3 (MI phase). In the MI phase, \(C_{\mathrm{DW}}\) vanishes exponentially with L, according to the fits: \(C_{\mathrm{DW}}^{g_{2}=1.35} \approx0.241 \times e^{-0.032 L}\) and \(C_{\mathrm{DW}}^{g_{2}=1.3} \approx0.253 \times e^{-0.065 L}\). In the DW phase, \(C_{\mathrm{DW}}\) converges to a finite value. The inset shows such obtained asymptotic value \(\mathcal{O}_{\mathrm{DW}}\), as a function of \(g_{2}\).

To be more accurate, in the right panel we performed a finite-size scaling and showed that the staggered correlation \(C_{\mathrm{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/\Omega=0.05\) and show that for \(g_{2}/\Omega= 1.3, 1.35\) the DW order is exponentially suppressed with L, while for \(g_{2}/\Omega= 1.4\) it remains finite. The \({\mathcal{O}}_{\mathrm{DW}}\) order parameter reached for \(L \to \infty\) is displayed in the inset as a function of \(g_{2}\).

In order to determine the DW boundary in the phase diagram of Figure 2, we adopted the following protocol. For a fixed value of \(t/\Omega\), we start increasing \(g_{2}\) from zero up to a finite value, with a fixed increment \(\delta g_{2} = 0.05 \Omega\), 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 Figure 2), for any fixed t, is located by the \(g_{2}^{*}(t)\) that gives the first non-vanishing order parameter \(\mathcal{O}_{\mathrm{DW}}\).

Here we stress that, because of the arrangement of our 1D array [see Figure 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 \(\vert 4 \rangle\) of the L-th atom will be never occupied, since the transition \(\vert 2 \rangle \leftrightarrow \vert 4 \rangle\) 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.

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 Figure 5). Quite remarkably, we note that a positive Δ will never stabilize an antiferromagnetic DW ordering.

Figure 5

DW phase order parameter as a function of detuning. For a positive detuning, the system will never present as a DW phase, meaning that \(C_{\mathrm{DW}}(L/2)=0\). Taking a negative value of the detuning, we observe that \(C_{\mathrm{DW}}(L/2)\) increases with \(\vert \Delta \vert \). The parameters in this figure are \(t/\Omega=0.05\), \(g_{1}/\Omega=0.8\), \(g_{2}/\Omega=1.6\) and the system size is \(L=100\).

Neutral gap

The analysis leading to the phase diagram in Figure 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 Figure 6 show the behavior of \(\Delta E_{n}\) as a function of \(g_{2}\), for a fixed value of \(t / \Omega\). In particular we analyzed a vertical cut in the phase diagram of Figure 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 \(\Delta E_{n}\), for \(t/\Omega\lesssim1.45\)) to the SF phase (zero \(\Delta E_{n}\), for \(1.45 \lesssim t/\Omega \lesssim1.8\)), and then to the DW phase (nonzero \(\Delta E_{n}\), for \(t/\Omega\gtrsim1.8\)).

Figure 6

Analysis of the neutral gap. Neutral energy gap as a function of \(g_{2}\) for \(t/\Omega=0.25\), i.e., along the vertical cut depicted in Figure 2. In the left panel, the different curves are for various system sizes according to the legend, and for a fixed number of kept states \(m=80\) in the DMRG algorithm. The right panel evidences the convergence of the data, at fixed \(L=80\), by increasing m (see also the inset, where we show the behavior of \(\Delta E_{n}\) as a function of m, for three different values of \(g_{2}\)).

While we cannot see a clear signature of a finite gap for \(g_{2} = 1.8 \Omega\), the scaling with the size displayed in the left panel of Figure 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 Figure 6. We observe that the non monotonic features that are visible in the region \(1.45 \lesssim t/\Omega\lesssim1.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).


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.


  1. 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.

  2. 2.

    The two points of \(\langle\delta n_{i}^{\mathrm{pol}} \delta n_{j}^{\mathrm{pol}} \rangle\), with \(\vert i - j\vert = 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).


  1. 1.

    Bloch I, Dalibard J, Zwerger W. Many-body physics with ultracold gases. Rev Mod Phys. 2008;80:885.

    ADS  Article  Google Scholar 

  2. 2.

    Fisher MPA, Weichman PB, Grinstein G, Fisher DS. Boson localization and the superfluid-insulator transition. Phys Rev B. 1989;40:546.

    ADS  Article  Google Scholar 

  3. 3.

    Jaksch D, Bruder C, Cirac JI, Gardiner CW, Zoller P. Cold bosonic atoms in optical lattices. Phys Rev Lett. 1998;81:3108.

    ADS  Article  Google Scholar 

  4. 4.

    Greiner M, Mandel O, Esslinger T, Hänsch TW, Bloch I. Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms. Nature. 2002;415:39.

    ADS  Article  Google Scholar 

  5. 5.

    Lewenstein M, Sanpera A, Ahufinger V, Damski B, Sen(De) A, Sen U. Ultracold atomic gases in optical lattices: mimicking condensed matter physics and beyond. Adv Phys. 2007;56:243.

    ADS  Article  Google Scholar 

  6. 6.

    Lahaye T, Menotti C, Santos L, Lewenstein M, Pfau T. The physics of dipolar bosonic quantum gases. Rep Prog Phys. 2009;72:126401.

    ADS  Article  Google Scholar 

  7. 7.

    Schmidt H, Imamoǧlu A. Giant Kerr nonlinearities obtained by electromagnetically induced transparency. Opt Lett. 1996;21:1936.

    ADS  Article  Google Scholar 

  8. 8.

    Imamoǧlu A, Schmidt H, Woods G, Deutsch M. Strongly interacting photons in a nonlinear cavity. Phys Rev Lett. 1997;79:1467.

    ADS  Article  Google Scholar 

  9. 9.

    Hartmann MJ, Brandão FGSL, Plenio MB. Strongly interacting polaritons in coupled arrays of cavities. Nat Phys. 2006;2:849.

    Article  Google Scholar 

  10. 10.

    Greentree AD, Tahan C, Cole JH, Hollenberg LCL. Quantum phase transitions of light. Nat Phys. 2006;2:856.

    Article  Google Scholar 

  11. 11.

    Angelakis DG, Santos MF, Bose S. Photon-blockade-induced Mott transitions and XY spin models in coupled cavity arrays. Phys Rev A. 2007;76:031805(R).

    ADS  Article  Google Scholar 

  12. 12.

    Hartmann MJ, Brandão FGSL, Plenio MB. Quantum many-body phenomena in coupled cavity arrays. Laser Photonics Rev. 2008;2:527.

    Article  Google Scholar 

  13. 13.

    Tomadin A, Fazio R. Many-body phenomena in QED-cavity arrays. J Opt Soc Am. 2010;27:A130.

    Article  Google Scholar 

  14. 14.

    Houck AA, Türeci HE, Koch J. On-chip quantum simulation with superconducting circuits. Nat Phys. 2012;8:292.

    Article  Google Scholar 

  15. 15.

    Schmidt S, Koch J. Circuit QED lattices: towards quantum simulation with superconducting circuits. Ann Phys. 2013;525:395.

    Article  Google Scholar 

  16. 16.

    Carusotto I, Gerace D, De Liberato S, Ciuti C, Imamoǧlu A. Fermionized photons in an array of driven dissipative nonlinear cavities. Phys Rev Lett. 2009;103:033601.

    ADS  Article  Google Scholar 

  17. 17.

    Tomadin A, Giovannetti V, Fazio R, Gerace D, Carusotto I, Türeci HE, Imamoǧlu A. Signatures of the superfluid-insulator phase transition in laser-driven dissipative nonlinear cavity arrays. Phys Rev A. 2010;81:061801(R).

    ADS  Article  Google Scholar 

  18. 18.

    Hartmann MJ. Polariton crystallization in driven arrays of lossy nonlinear resonators. Phys Rev Lett. 2010;104:113601.

    ADS  Article  Google Scholar 

  19. 19.

    Nunnenkamp A, Koch J, Girvin SM. Synthetic gauge fields and homodyne transmission in Jaynes-Cummings lattices. New J Phys. 2011;13:095008.

    Article  Google Scholar 

  20. 20.

    Nissen F, Schmidt S, Biondi M, Blatter G, Türeci HE, Keeling J. Nonequilibrium dynamics of coupled qubit-cavity arrays. Phys Rev Lett. 2012;108:233603.

    ADS  Article  Google Scholar 

  21. 21.

    Grujic T, Clark SR, Angelakis DG, Jaksch D. Non-equilibrium many-body effects in driven nonlinear resonator arrays. New J Phys. 2012;14:103025.

    Article  Google Scholar 

  22. 22.

    Grujic T, Clark SR, Jaksch D, Angelakis DG. Repulsively induced photon superbunching in driven resonator arrays. Phys Rev A. 2013;87:053846.

    ADS  Article  Google Scholar 

  23. 23.

    Jin J, Rossini D, Fazio R, Leib M, Hartmann MJ. Photon solid phases in driven arrays of nonlinearly coupled cavities. Phys Rev Lett. 2013;110:163605.

    ADS  Article  Google Scholar 

  24. 24.

    Jin J, Rossini D, Leib M, Hartmann MJ, Fazio R. Steady-state phase diagram of a driven QED-cavity array with cross-Kerr nonlinearities. Phys Rev A. 2014;90:023827.

    ADS  Article  Google Scholar 

  25. 25.

    Underwood DL, Shanks WE, Koch J, Houck AA. Low-disorder microwave cavity lattices for quantum simulation with photons. Phys Rev A. 2012;86:023837.

    ADS  Article  Google Scholar 

  26. 26.

    Abbarchi M, Amo A, Sala VG, Solnyshkov DD, Flayac H, Ferrier L, Sagnes I, Galopin E, Lemaître A, Malpuech G, Bloch J. Macroscopic quantum self-trapping and Josephson oscillations of exciton polaritons. Nat Phys. 2013;9:275.

    Article  Google Scholar 

  27. 27.

    Toyoda K, Matsuno Y, Noguchi A, Haze S, Urabe S. Experimental realization of a quantum phase transition of polaritonic excitations. Phys Rev Lett. 2013;111:160501.

    ADS  Article  Google Scholar 

  28. 28.

    Lucero E, Barends R, Chen Y, Kelly J, Mariantoni M, Megrant A, O’Malley P, Sank D, Vainsencher A, Wenner J, White T, Yin Y, Cleland AN, Martinis JM. Computing prime factors with a Josephson phase qubit quantum processor. Nat Phys. 2012;8:719.

    Article  Google Scholar 

  29. 29.

    Steffen L, Salathe Y, Oppliger M, Kurpiers P, Baur M, Lang C, Eichler C, Puebla-Hellmann G, Fedorov A, Wallraff A. Deterministic quantum teleportation with feed-forward in a solid state system. Nature. 2013;500:319.

    ADS  Article  Google Scholar 

  30. 30.

    Chen Y, Roushan P, Sank D, Neill C, Lucero E, Mariantoni M, Barends R, Chiaro B, Kelly J, Megrant A, Mutus JY, O’Malley PJJ, Vainsencher A, Wenner J, White TC, Yin Y, Cleland AN, Martinis JM. Emulating weak localization using a solid-state quantum circuit. Nat Commun. 2014;5:5184.

    ADS  Article  Google Scholar 

  31. 31.

    Zueco D, Mazo JJ, Solano E, García Ripoll JJ. Microwave photonics with Josephson junction arrays: negative refraction index and entanglement through disorder. Phys Rev B. 2012;86:024503.

    ADS  Article  Google Scholar 

  32. 32.

    Peropadre B, Zueco D, Wulschner F, Deppe F, Marx A, Gross R, García Ripoll JJ. Tunable coupling engineering between superconducting resonators: from sidebands to effective gauge fields. Phys Rev B. 2013;87:134504.

    ADS  Article  Google Scholar 

  33. 33.

    Sowiński T, Dutta O, Hauke P, Tagliacozzo L, Lewenstein M. Dipolar molecules in optical lattices. Phys Rev Lett. 2012;108:115301.

    ADS  Article  Google Scholar 

  34. 34.

    Dalla Torre EG, Berg E, Altman E. Hidden order in 1D Bose insulators. Phys Rev Lett. 2006;97:260401.

    ADS  Article  Google Scholar 

  35. 35.

    Berg E, Dalla Torre EG, Giamarchi T, Altman E. Rise and fall of hidden string order of lattice bosons. Phys Rev B. 2008;77:245119.

    ADS  Article  Google Scholar 

  36. 36.

    Rossini D, Fazio R. Phase diagram of the extended Bose-Hubbard model. New J Phys. 2012;14:065012.

    Article  Google Scholar 

  37. 37.

    Deng X, Citro R, Orignac E, Minguzzi A, Santos L. Polar bosons in one-dimensional disordered optical lattices. Phys Rev B. 2013;87:195101.

    ADS  Article  Google Scholar 

  38. 38.

    Wikberg E, Larson J, Bergholtz EJ, Karlhede A. Fractional domain walls from on-site softening in dipolar bosons. Phys Rev A. 2012;85:033607.

    ADS  Article  Google Scholar 

  39. 39.

    Batrouni GG, Scalettar RT, Rousseau VG, Grémaud B. Competing supersolid and Haldane insulator phases in the extended one-dimensional bosonic Hubbard model. Phys Rev Lett. 2013;110:265303.

    ADS  Article  Google Scholar 

  40. 40.

    Batrouni GG, Rousseau VG, Scalettar RT, Grémaud B. Competing phases, phase separation, and coexistence in the extended one-dimensional bosonic Hubbard model. Phys Rev B. 2014;90:205123.

    ADS  Article  Google Scholar 

  41. 41.

    Schiró M, Bordyuh M, Öztop B, Türeci HE. Phase transition of light in cavity QED lattices. Phys Rev Lett. 2012;109:053601.

    ADS  Article  Google Scholar 

  42. 42.

    Kumar B, Jalal S. Quantum Ising dynamics and Majorana-like edge modes in the Rabi lattice model. Phys Rev A. 2013;88:011802(R).

    ADS  Article  Google Scholar 

  43. 43.

    Schollwöck U. The density-matrix renormalization group. Rev Mod Phys. 2005;77:259.

    ADS  MathSciNet  Article  MATH  Google Scholar 

  44. 44.

    De Chiara G, Rizzi M, Rossini D, Montangero S. Density matrix renormalization group for dummies. J Comput Theor Nanosci. 2008;5:1277.

    Article  Google Scholar 

  45. 45.

    Rossini D, Fazio R. Mott-insulating and glassy phases of polaritons in 1D arrays of coupled cavities. Phys Rev Lett. 2007;99:186401.

    ADS  Article  Google Scholar 

  46. 46.

    Rossini D, Santoro GE, Fazio R. Photon and polariton fluctuations in arrays of QED-cavities. Europhys Lett. 2008;83:47011.

    ADS  Article  Google Scholar 

  47. 47.

    D’Souza AG, Sanders BC, Feder DL. Fermionized photons in the ground state of one-dimensional coupled cavities. Phys Rev A. 2013;88:063801.

    ADS  Article  Google Scholar 

  48. 48.

    Hu Y, Ge GQ, Chen S, Yang XF, Chen YL. Cross-Kerr-effect induced by coupled Josephson qubits in circuit quantum electrodynamics. Phys Rev A. 2011;84:012329.

    ADS  Article  Google Scholar 

  49. 49.

    Rebić S, Twamley J, Milburn GJ. Giant Kerr nonlinearities in circuit quantum electrodynamics. Phys Rev Lett. 2009;103:150503.

    ADS  Article  Google Scholar 

  50. 50.

    Verstraete F, García-Ripoll JJ, Cirac JI. Matrix product density operators: simulation of finite-temperature and dissipative systems. Phys Rev Lett. 2004;93:207204.

    ADS  Article  Google Scholar 

  51. 51.

    Zwolak M, Vidal G. Mixed-state dynamics in one-dimensional quantum lattice systems: a time-dependent superoperator renormalization algorithm. Phys Rev Lett. 2004;93:207205.

    ADS  Article  Google Scholar 

  52. 52.

    Kühner TD, Monien H. Phases of the one-dimensional Bose-Hubbard model. Phys Rev B. 1998;58:R14741.

    ADS  Article  Google Scholar 

  53. 53.

    Kühner TD, White SR, Monien H. One-dimensional Bose-Hubbard model with nearest-neighbor interaction. Phys Rev B. 2000;61:12474.

    ADS  Article  Google Scholar 

Download references


We would like to acknowledge our previous collaboration with M. Hartmann and M. Leib which was inspiring for the present work. This work was supported by Italian MIUR via FIRB Project RBFR12NLNA and PRIN Project 2010LLKJBX, by EU through IP-SIQS, and by National Natural Science Foundation of China under Grant No. 11175033 and No. 11305021.

Author information



Corresponding author

Correspondence to Jiasen Jin.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

All the authors participated in the design of the research, analysis of the results, and writing of the paper. The DMRG code used to run all the simulations of this research has been developed and written by DR and coworkers [44] (see also The DMRG simulations were performed by JJ.

Rights and permissions


Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Jin, J., Fazio, R. & Rossini, D. Phase diagram of a QED-cavity array coupled via a N-type level scheme. EPJ Quantum Technol. 2, 5 (2015).

Download citation

PACS Codes

  • 42.50.Pq
  • 85.25.Cp
  • 64.70.Tg


  • cavity quantum electrodynamics
  • strongly correlated polariton systems
  • quantum phase transitions