Quantum phase diagrams of the Jaynes–Cummings Hubbard models in non-rectangular lattices

In this paper, we investigate systematically the quantum phase transition between the Mott-insulator and superfluid states of the Jaynes–Cummings Hubbard model in triangular, square, honeycomb and kagomé lattices. With the help of Green’s function method, by treating the hopping term in the Jaynes–Cummings Hubbard model as perturbation, we calculate the phase boundaries of Jaynes–Cummings Hubbard models on different geometrical lattices analytically up to second order for both detuning Δ=0 and Δ≠0. It is shown that the Mott phase becomes more unstable when increasing the polariton number n. Meanwhile, the Mott lobe of n  =  0 is enlarged while other Mott lobes shrink as the detuning parameter Δ increases.


Introduction
As a new platform for testing and simulating quantum many body physics, systems of cavity quantum electrodynamics (cavity QED) provide new insights about the quantum aspect of the interaction between photons and atoms [1][2][3][4]. The cavity QED systems consisting of an array of coupled optical cavities with each cavity containing a two-level atom are often described by the Jaynes-Cummings Hubbard model (JCHM) [2,5,6]. The coupling between the atom and photons results in an effective photon-photon repulsion interaction [7][8][9][10][11][12]. The delicate balance between the hopping of photons and photon-atom on-site interaction leads to a quantum phase transition between superfluid (SF) and Mott-insulator (MI) of polaritons [5].
Although the experimental realization of JCHM with trapped ions was performed recently [13], the JCHM has been widely studied theoretically in the past few years [5,[14][15][16][17][18][19][20][21][22][23][24][25][26]. The quantum phase diagrams of the JCHM has been extensively studied via a wide range of methods, such as the mean field theory [5,20], Ginzburg-Landau approach [22], strong coupling theory [19,21] and numerical methods [23][24][25][26]. Despite so many efforts, these outstanding works have mainly focused on mean-field level or systems with simple geometrical lattices. Due to the potential applications in quantum simulation, it is thus desirable to investigate the quantum phase transitions of JCHM in more complex lattice structures beyond mean-field level in a systematic way.
As is well known, the one-particle Green's function diverges at the phase transition point of second order phase transitions [27]. Based on strong coupling perturbation theory [28,29], the perturbative calculation of one-particle Green's function has been used to determine quantum phase diagrams of JCHM systems in square and cubic lattices [19]. In this In this paper, we investigate systematically the quantum phase transition between the Mottinsulator and superfluid states of the Jaynes-Cummings Hubbard model in triangular, square, honeycomb and kagomé lattices. With the help of Green's function method, by treating the hopping term in the Jaynes-Cummings Hubbard model as perturbation, we calculate the phase boundaries of Jaynes-Cummings Hubbard models on different geometrical lattices analytically up to second order for both detuning 0 ∆ = and 0 ∆ ≠ . It is shown that the Mott phase becomes more unstable when increasing the polariton number n. Meanwhile, the Mott lobe of n = 0 is enlarged while other Mott lobes shrink as the detuning parameter ∆ increases.
Keywords: Jaynes-Cummings Hubbard model, quantum phase transition, Bose-Einstein condensates, polariton (Some figures may appear in colour only in the online journal)

Astro Ltd
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. paper, by treating the hopping part of the JCHM as perturbation, with the help of Green's function technique and linked cluster expansion [30], the SF-MI quantum phase boundaries of JCHM in square, triangular, honeycomb and kagomé lattices are investigated systematically for various polariton numbers and detuning parameters. By taking the leading quantum correction into account, the analytical results beyond mean-field and the corresponding phase diagrams are presented.

The model
The cavity QED systems can often be described by the following Jaynes-Cummings-Hubbard Hamiltonian (setting which consists of two parts [2,5,6]: the local one H 0 which is a sum over localized Jaynes-Cummings Hamiltonians, and the hopping part H 1 describing the tunneling of photons between different neighbor cavities, i.e. where a ĵ (a ĵ † ) is photon annihilation (creation) operator on lattice site j and ĵ σ − ( ĵ † σ ) is lowering (raising) operator between the ground state g⟩ | and excitation state e⟩ | of the two-level atom in the cavity. The atom-photon coupling g gives rise to the formation of polaritons with n a a ĵˆˆˆ † † σ σ = + − denoting the number operator of the polaritons on site j [19,20]. ω is the frequency of the one-mode electromagnetic field with its polarization being fixed. ω ∆ = − ε represents the detuning parameter with ε being transition frequency between the two energy levels g⟩ | and e⟩ | of an atom in optical cavities. μ is chemical potential and ij κ describes the hopping strength of photons between sites i and j. Here we only consider the nearest neighbor tunneling (which is denoted by 〈 〉 i j , ), so we define the tunneling parameter as follows Due to the non-commutativity of H 0 and H 1 , the exact eigenstates of Ĥ can hardly be determined. In the atomic limit ( 0 → κ ), the Hamiltonian decouples at the site index and all the excitations remain at their respective lattice sites, thus the ground state of the lattice system is a MI phase which is labeled by polaritons number. While in the hopping limit (g 0 → ), the tunneling of photons between the nearest neighbor cavities overwhelms the coupling between cavity mode and two-level atom, the ground state becomes photons SF, where the photons are delocalized over the lattices in real space. With the increase of g / κ , a phase transition from MI to SF is thus expected. In the following, by treating the hopping term H 1 as a perturbation, we will investigate the MI-SF quant um phase transitions for JCHM systems with various lattice structures.
Before performing our perturbation calculations, let us clarify the ground states first, which can be determined by diagonalizing the local Jaynes-Cummings Hamiltonian [22,31]. Note that the site index j will be ommited to simplify notation for the reminder of the paper. From the expression of H 0 , it can be easily found that the occupation number operator n commutes with H 0 , i.e. n H , 0 0 [ˆˆ] = . Hence we know that i) the occupation number n will be a conserved quantity; ii) H 0 and n share a common set of eigenstates. The number-conserved particles are named polariton in unperturbed systems H 0 . In order to find the eigenstates, it is necessary to examine the occupation number operator n n n p âˆ= + further with n a a pˆ † = and n âˆ † σ σ = − . Since operators n p and n â commute with each other, they also share a set of common eigenstates. These eigenstates are often called bare basis states which can be constructed with the photonic eigenstates n p ⟩ | and two-level states of the atom s⟩ | as The eigenstates n s , p ⟩ | satisfy the orthogonality and completeness conditions, i.e.
According to the eigenvalue equation n n s n n s n n n s n n s n , , , 0 , 1, 2 , , , , 0 , 1.
we have n n s n n s , , , with n n n p a ( ) = + . From equation (8), we know that for a fixed number n of polaritons there exist two possible micro states, i.e. n, 0⟩ | and n 1, 1⟩ | − . With these bare basis states, the Jaynes-Cummings Hamiltonian can be rewritten as

H ns n s H n s n s
n n e n e n e n e n n g n g g n n e n g g n n g n e , , The above expression can be transformed to a simpler matrix form in the bare basis which consists of 2 2 × block matrices Now the task of calculating the characteristic equation of H 0 is converted to solve the eigenvalues of each 2 2 × matrix labeled by the number of polaritons After straightforward calculations, we get the following eigenstates and eighenvalues of H 0 , i.e. for n = 0 with the generalized Rabi frequency R ng 4 for fixed polaritons number. By adjusting the chemical potential, one eventually reaches a point where adding an excitation to the system becomes energetically favorable. It imposes constraints on the chemical potential μ in each Mott region, i.e.
The above equation leads to for n = 0 and

Cumulants expansion method
With the configurations of the ground states of the unperturbed Hamiltonian H 0 , we can now tune on the value of κ gradually. Under the influence of quantum fluctuations, the localized states will become more and more softened. When g / κ goes beyond some critical value, a quantum phase transition is expected to happen, and the system will enter a phase of polariton superfluid. To study the quantum phase diagrams of JCHM, we are going to treat H 1 as a perturbation and calculate the one-particle Green's function of the systems in Dirac picture. In such a case, the one-particle Green's function is represented as [30,32] Tr e ,0 , By expanding the evolution operator in equation (18) perturbatively, the expansion expression of the Green's function being the (n + 1)-particle Green's function with respect to the unperturbed Hamiltonian H 0 . Now the problem of calculating the one-particle Green's function with respect to the full Hamiltonian (1) is reduced to calculating the (n + 1)-particle Green's function G n 1 0 ( ) + . Due to the fact that the Hamiltonian H 0 is diagonal in the site representation, each annihilation (creation) operator in the expression of G n 1 0 + has to be paired with a creation (annihilation) operator on the same site. Therefore, G n 0 ( ) can be decomposed in terms of the cumulants C ; ; ; ; in which all the operators are at the same lattice site [30,33,34]. For example, the 1-particle Green's function G 1 0 ( ) of the unperturbed system can be expressed as with the 1st-order cumulant C T a a . In fact, the cumulants can in principle be calculated explicitly. In order to illustrate this point clearly, as an example, the first and second order cumulants will be calculated in detail. According to the definition (22), we have , We need know how the ladder operator â and â † act on the the dressed state. The relationship between the dressed state and the bare state is where the matrix n R is , cos , sin , .
where T n is the transition matrix. On the other hand, we have the relation a n n n g n n e a n n n g n n e Therefore, we get the relation n n n a a n b a n a b n b b n a a n b a n a b n b b In fact, the calculation of the one-particle Green's function can be conducted conveniently in Matsubara space [30,34]. This is mainly because the convolution of cumulants over the imaginary time variables amounts to simple multiplication of cumulants. With the help of the following integral formula After a complicated yet straightforward calculation, we get All other higher order cumulants can be calculated in the same manner, yet extremely complicated and tedious.
With the cumulants in hand, the Green's function can then be calculated out. To avoid writing the tedious complex expression, a diagrammatic representation may be adopted. We denote an nth-order cumulant at a lattice site by a vertex with n incoming lines and n outgoing lines with imaginary time variables. For example, the 1st-and 2nd-order cumulants are represented diagrammatically as Meanwhile, the hopping parameter ij κ between the nearest neighbor sites (i and j) is represented as Now the (n + 1)-particle Green's function (20) can be diagrammatically expressed. Taking into account of the cancelation effect of the numerator and the denominator in one-particle Green's function (17), one can find that the one-particle Green's function only consists of connected diagrams [30,33]. The one-particle Green's function can then be expanded as a geometric series of building block which is the sum of all one-particle irreducible diagrams with their respective weights [27,35], i.e.
with the building block C k , being represented diagrammatically as

The phase diagram of the JCHM
Since the one-particle Green's function of the systems diverges at the phase transition point of a second order phase transition [27,35], it can be used to locate the MI-SF phase boundary. As shown in equation (46), the one-particle Green's function is a series of hopping parameter κ. Thus the radius of convergence of the series will reveal the location of the phase boundary. Combining the fact that the phase transitions are governed by long-wavelength fluctuations [27,35], i.e. k 0 = and 0 m ω = , according to equation (46) we have the phase boundary to the first order  where z is the coordination number and G 1 loop ( ) is the lowest quantum correction which can be represented diagrammatically as . The expression of the loop is where z is the coordination number. The zero temperature phase diagrams of JCHM on square, triangular, honeycomb and Kogomé lattices for different polaritons numbers at zero detuning ( 0 ∆ = ) are shown in figures 1 (log-linear plot) and 2 (linear-linear plot). From the pictures we see that the Mott lobes becomes more and more smaller with the increase of the polariton number n.  This is because effective polaritons repulsion decreases with n increasing. Moreover, we find that the quantum fluctuation correction is the strongest in the case with n = 1. With the increase of coordination number, the correction of the lowest order quantum fluctuation decreases gradually which is consistent with the numerical results [21].
In figure 3, for the detuning parameter 0 ∆ ≠ , we show the second order perturbation results of the quantum phase diagrams of JCHM on square, triangular, kagomé and honeycomb lattices. We can see clearly that the size of the Mott lobe with n = 0 enlarges while the Mott lobe with n = 1 shrinks as the detuning parameter ∆ increases.

Summary
In this work, by treating the hopping parameter of photons in the Jaynes-Cummings-Hubbard model (JCHM) as a perturbation, we have calculated the one-particle Green's function of JCHM analytically with the help of cumulant expansion. By use of the Green's function d'Alembert ratio test method, we systematically determined the quantum phase diagrams of JCHM on different geometrical lattices up to the second order for both detuning 0 ∆ = and 0 ∆ ≠ with the explicit expression of the lowest quantum correction. We found that the lowest quantum correction is the strongest in the case with n = 1. Meanwhile, the Mott lobes become more unstable with the increase of polariton number n. On the other hand, when the detuning 0 ∆ ≠ , the Mott lobe for n = 1 shrinks as the detuning parameter ∆ increases.