Effective Three-Body Interactions in Jaynes-Cummings-Hubbard Systems

A generalisation of the Jaynes-Cummings-Hubbard model for coupled-cavity arrays is introduced, where the embedded two-level system in each cavity is replaced by a Ξ-type three-level system. We demonstrate that the resulting effective polariton-polariton interactions at each site are both two-body and three-body. By tuning the ratio of the two transition dipole matrix elements, we show that the strength and sign of the two-body interaction can be controlled whilst maintaining a three-body repulsion. We then proceed to demonstrate how different two-body and three-body interactions alter the mean field superfluid-Mott insulator phase diagram, with the possible emergence of a pair superfluid phase in the two-body attractive regime.

Pfaffian states via cavity QED 14 . As we assume that the three-level system is in the Ξ-configuration, we allow for photon-mediated atomic transitions only between the |0 a 〉 → |1 a 〉 and |1 a 〉 → |2 a 〉 atomic levels; we denote the corresponding transition dipole moments as β 01 and β 12 respectively, with β 02 = 0. We also introduce the photonic creation and annihilation operators ˆ † a and â, along with the atomic operators λ λ = ↓ + ↓ − † ( ) and λ λ = ↑ + ↑ − † ( ), where λ ↓ + and λ ↑ + excite the atom from |0 a 〉 → |1 a 〉 and |1 a 〉 → |2 a 〉 respectively. Subsequently, with β β β =  / 12 01 , the application of the rotating wave approximation results in the Hamiltonian for a single site taking the form: )/ 01 represent the rescaled chemical potential and cavity detuning respectively. Since n commutes with Ĥ JC3 , the Hamiltonian for this system may be diagonalised by decomposing it into a direct sum of operators Ĥ n JC3 , each acting on the Fock subspace spanned by basis states | 〉 ≡ | 〉 ⊗ | 〉 n n n n , a ph a ph where n a + n ph = n. For each n-particle subspace, we work in the ordered basis | 〉 | − 〉  n i n i { 0, , , , }, where i = 0 for n = 0, i = 1 for n = 1 and i = 2 for n ≥ 2. For the vacuum, where n = 0, Ĥ JC 0 3 is trivially zero. For the one-particle subspace (n = 1), λ ↑ ± is undefined and the usual Pauli matrix-based representation applies to the λ ↓ ± operators 26 . However, for n ≥ 2, these new atomic operators are representable via the Gell-Mann matrices, generators of su (3) 29 , given by This yields the matrix blocks to be diagonalised for each n, given by We initially detail the nature of the spectrum of Eq. (1). The spectra for n = 0, 1 are necessarily identical to those in the Jaynes-Cummings model 27 with the existence of two distinct eigenstates for n = 1. However, it is readily seen from Eq. (2) that there must exist three independent eigenstates for n > 1. The energies of each of these branches is 30,31 given by n n n 3 The Fock state representations for the corresponding polariton branches are provided by 30 where for n ≥ 2 the normalisation factor ζ k n ( ) is defined as It is necessary to be able to identify which branch (k) has the lowest energy for a given polariton number n. If the arccosine in Eq. (7) is restricted to its principle branch, i.e. [0, π], then θ ∈       ∀ ≥ π n 0, 2 n 3 and thus Therefore, the ground state for any polariton number n ≥ 2 may be taken to be the k = 1 branch, or equivalently the '−' branch for n = 1, and for future reference we denote the dimensionless energies corresponding to these branches as E g (n).
An effective two-body repulsion for photons is provided in a Jaynes-Cummings cavity by the interactions between the two-level atom and the photons; this is the cause of the photon blockade effect 32,33 . In the three-level system coupled to a single cavity mode, the existence of the third atomic level, and its interaction with the cavity mode, induces an additional three-body nonlinearity 14 . Furthermore, in this case, the sign of neither the two-body or three-body nonlinearities are immediately evident, but may nonetheless be discerned from inspecting the eigenenergies of the lowest energy polariton branches for each occupation n. For each n, we initially define the quantity E g (n) to be the eigenenergy of the respective lowest polariton branch, i.e.
∀ ≥ E n n ( ) 2 1 , E − (1) for n = 1 and E(0) for n = 0. The quantity = ′ E E n n ( )/ n g is thus the energy of each polariton quasiparticle in the ground state corresponding to an occupation of n polaritons in a cavity. Let us initially consider the case where there is a single polariton in the cavity. If the two-body nonlinearity is repulsive, there must be an energy cost associated with adding a second polariton to the cavity, i.e. > ′ ′ E E 2 1 . Conversely, a two-body attraction is associated with an energy cost associated with removing a polariton from a cavity with two polaritons, a scenario marked by the property < ′ ′ E E 2 1 . We thus plot ′ E n as a function of β  , at ∆ = ∼ 0, in Fig. 2(a), in order to determine the sign of the two-body nonlinearity. From this it is evident that as β  is increased, the proportional energies for one (the dotted curve) and two (the blue curve) polaritons cross, suggesting that a crossover of the two-body interaction, from repulsion to attraction, occurs at a critical value of β  . As for the three-body nonlinearity, it is readily seen that E E , for all β  . Since there is always an energy cost associated with adding an additional polariton to a cavity with two polaritons already present, we conclude that the three-body nonlinearity is not only repulsive, but sufficiently so to ensure the presence of the energy cost even when the two-body nonlinearity is an attraction.
In the Methods section, we show that taking the difference of the proportional energy ′ E n for n = 1, 2 yields the . Note that this result may also be obtained from a mapping onto a two-mode bosonic system, with intermode tunnelling and a three-body constraint 14 . That the two-body repulsion-attraction crossover is independent of ∆ ∼ is numerically demonstrated in Fig. 2 In this regime, the energy per polariton is degenerate for n = 1 and n = 2 regardless of ∆ ∼ , which is evident in the coincidence of the blue and dotted lines in Fig. 2(b). This figure also demonstrates that the three-body nonlinearity is a repulsion, irrespective of β  and ∆ ∼ .

Jaynes-Cummings-Hubbard Mean-Field Phase Diagram.
With the exact single-site solutions in hand we now consider a lattice of such systems with photon hopping permissible between nearest neighbours. For such a system, we wish to compute the resulting zero-temperature mean-field phase diagram. A reasonable first guess would be that the superfluid and Mott insulator phases also exist in this modified system, and thus we proceed to apply the mean-field approximation to the hopping terms to yield a purely local Hamiltonian. Explicitly, this involves the substitution of the expression 34 i where z is the number of nearest-neighbour sites, a quantity that is dependent on the lattice geometry (the superfluid order parameter ψ is assumed to be real). This yields an energy approximation consisting of a sum over all sites of the Hamiltonian To numerically diagonalise the Hamiltonian in Eq. (16), we employ a self-consistent field approach, which we shall briefly detail here. Initially, we truncate the single-site basis to Fock states |n a , n ph 〉 such that + ≤ n n N a p h max , where, in this work, we specify the maximum polariton occupation to be N max = 16. With an initial guess for ψ, we diagonalise Ĥ MF JCH3 with respect to this truncated basis, and check that the resulting ground state respects ψ 〈 〉 = a . If it does not, the obtained 〈 〉 a is used as a new guess for ψ, and the process is repeatedly iterated until adequate convergence is achieved.
Subsequently, the phase diagram may be obtained from the solution space by plotting the self-consistently determined values of ψ in a given region of parameter space. However, analytical expressions for the phase boundaries may also be found via Landau theory, provided that the transitions are second-order continuous. Following Landau's prescription, the free energy (which reduces to the many-body ground state energy in the zero temperature limit) may be expressed as a power series in the order parameter ψ, i.e 5,34 . At a continuous phase transition, minimisation of the free energy with respect to ψ yields the criterion that the phase boundary is demarcated by solutions to the curve 35 a 2 = 0. Solutions to this are especially simple to find if the order parameter ψ is manifest in a perturbative expansion of the ground state energy as a power series, which one finds to be the case in Hubbard-type systems 34 . To construct this expansion, we note that Eq. (16) is purely diagonal except for the term κψ + ˆ † a a ( ) , and so it is partitioned as a 0 2 ) . Then, we can employ Rayleigh-Schrodinger perturbation theory to expand each perturbed polariton eigenstate's energy in a power series of κψ  of the form 36 κψ For the unperturbed reference eigenstates Eqs (9, 10, 11), the quantity 〈n〉 is clearly a good quantum number, implying that, by Eq. (22), the first-order correction E g (1) vanishes identically. Subsequently, a comparison of Eqs (17) and (20) suggests that we make the identification that 5 The calculation of the second-order perturbative correction is somewhat tedious; we emphasise that in the summation over accessible virtual states in Eq. (23), the contributions from every polariton branch must be considered and not merely the lowest-energy branches for a given polariton occupation. For the sake of clarity, the full expressions for E n (2) can be found in the Methods section. Below, we merely compare the predictions for the phase boundary obtained via the numerical self-consistent field and analytical Landau-perturbative methods. For simplicity, the analysis is limited to the system with zero detuning, that is, ∆ = ∼ 0. Nontheless, it is possible to extend these results to a nonzero detuning in the analytical method, as well as in the numerical formalism.
Two-Body Repulsion. We first investigate the regime in which the two-body interaction is predicted to be repulsive, i.e. β <  2 . In Fig. 3, the superfluid-to-Mott insulator phase diagram is presented at ∆ = ∼ 0 for the regimes β = .  {1, 1 35, 2}. It is clear that for β =  1 (the first column in Fig. 3), the phase diagram closely resembles that of the two-level JCH system, i.e. successive Mott lobes corresponding to successive integer 〈 〉 n . However, as β  is increased to 1.35 (the middle column in Fig. 3), the first Mott lobe (〈 〉 = n 1) shrinks. At β =  2 the 〈 〉 = n 1 Mott lobe vanishes, with the 〈 〉 = n 2 Mott lobe being adjacent to the vacuum (〈 〉 = n 0).

Two-Body Attraction.
We now focus on the regime in which we predict the existence of an attractive two-body interaction, that is, when β ≥  2. As β →  2 from below, it was found, see Fig. 3, that the first Mott lobe shrank and subsequently vanished at the threshold of β =  2 , raising the question of what features might be seen in the attractive regime. Retaining ∆ = ∼ 0, we cross the threshold and examine the numerical and analytical predictions for the phase diagram at β = .  1 5 in Fig. 4(a). The analytical prediction is a first Mott lobe in the negative κ  half-plane, which we can reject as unphysical, and is thus not included in the figure. Of more concern is the prediction that the vacuum lobe and second Mott lobe are overlapping, evident in the overlap of the white lines in Fig. 4(a) in the region µ κ ′ ∈ − . ∈ .   { 1 06, 1}, {0, 0 35}, which is also clearly unphysical but not representative of a discardable solution. A plot of the self-consistently determined value of 〈 〉 n , in Fig. 4(b), seems to resolve this question with a prediction that the zeroth and second Mott lobes do not overlap, but abruptly meet each other with a boundary at κ ≥  0.

Discussion
The discrepancy between the numeric and analytical predictions, for β >  2 , is quite unlike what one observes in the Bose-Hubbard system 20,34 , or the JCH system 5 , or even the system in consideration for β ≤  2 . In these systems it was always assumed that the phase transition is second-order, that is, the order parameter ψ is continuous across the transition. A discontinuous, first-order phase transition features a jump in ψ from zero in the disordered phase to a distinct non-zero value in the ordered phase, and the condition a 2 = 0 does not describe the value of the phase boundary for such transitions 37 .
To understand the nature of the transition it is instructive to recall from Fig. 4(b) that the Mott lobe adjacent to the vacuum corresponds to 〈 〉 = n 2. The complete absence of the n = 1 Mott lobe for β ≥  2 must be attributable to the fact that in the two-body attractive regime, there is no value of µ ′  in the limit of zero hopping at which it is energetically favourable for there to be a single polariton at a given site, rather than zero or two of them.  However, the existence of the finite-length boundary, as a function of κ  (see Fig. 4(b)), between these lobes raises the question of what phase exists between them. In general, two Mott lobes meet at the value of µ ′  such that the two corresponding lowest polariton branches are degenerate. A reasonable expectation would therefore be that the quantity ϕ = 〈 〉 ∀â a i (26) i i is nonzero along this boundary. This corresponds to the mean-field approximation of the order parameter for the pair superfluid phase, in which superfluidity occurs due to the condensation of (quasi-)particle pairs 25 . The true pair superfluid phase occurs when ψ = 0 and ϕ ≠ 0, and is associated with the partial breaking of the system's (1) symmetry to U Z Z (1)/( /2 ) 24 . This anomalous behaviour in the phase transitions of the three-level JCH system when β >  2 is not entirely unexpected. While the precise details are different, similar behaviour has been noted in the Bose-Hubbard model with attractive onsite two-body interactions and an additional onsite three-body repulsion 38 . The relevant Hamiltonian, scaled appropriately such that the interaction strengths are dimensionless, is provided by where κ µ   , and ∼ W are all non-negative. This model has been studied in the context of experimental proposals which suggest that accounting for the three-body loss rate in a Bose-Hubbard model with attractive onsite two-body interactions may potentially result in either an effective finite three-body repulsion mechanism 19 or a hardcore three-body constraint 17 . We note that in the Hamiltonian described by Eq. (27), a hardcore three-body constraint may be modelled by taking the limit → ∞ ∼ W . In the case where ∼ W is finite, a zero-temperature phase diagram has previously been obtained for the one-dimensional chain using exact diagonalisation and the Density Matrix Renormalisation Group (DMRG) 25 . Furthermore, in the hardcore limit (W → ∞), a sophisticated procedure of spin-model mapping followed by the derivation of an effective field theory has provided a rich understanding of both the phase diagram and the dynamics of the system in each phase 21,22 . These studies demonstrate that in the presence of an additional three-body repulsion, the two-body-attractive Bose-Hubbard system exhibits both superfluid (ψ ≠ 0, ϕ ≠ 0) and pair superfluid (ψ = 0, ϕ ≠ 0) phases, with the phase transition between the pair superfluid and the Mott phase being first-order 17,[21][22][23][24][25] . Noting the similarities between this system and the three-level JCH system studied here, we predict that when β >  2 , the three-level JCH system exhibits the pair superfluid phase and an associated pair superfluid-Mott insulator first-order phase transition.
However, it is also noted that in these previous studies of similar systems, the phase boundary between the zeroth and second Mott phases is resolved in such a manner that the pair superfluid phase is not merely confined to a straight line at µ = − .  0 5 but in a narrow but nonetheless finite band of parameter space. It is not possible to explore this feature rigorously via the mean-field formalism employed in this work due to the somewhat crude approximation, made vis-a-vis spatially off-diagonal correlations, inherent in the mean-field approximation. It would therefore seem that utilising one or more of the methods employed for the three-body-repulsive, two-body-attractive Bose-Hubbard model-which include an effective potential mechanism with both single and paired source fields 23 , and quantum Monte Carlo methods 24 -would serve to elucidate the possible pair superfluid phase in the three-level JCH system.
The addition of a third energy level to each 'atom' in a Jaynes-Cummings cavity array is seen here to result in an additional effective three-body repulsion between the resonant cavity photons, analogous to existing modified Bose-Hubbard models with an additional onsite three-body repulsion. By uniformly adjusting the ratio of the two transition dipole moments of each atom, we predict that a crossover manifests itself, between the regimes of two-body repulsion and attraction. Strikingly, this is seen to result in a possible first-order phase transition between the superfluid and Mott vacuum phases, and the coexistence of first-and second-order superfluid-Mott insulator phase transitions across the second Mott lobe. Furthermore, the first Mott lobe vanishes in the regime of two-body attraction.
However, the presence of a pair superfluid phase, which our results seem to suggest, is not provable via the methods employed here. Its existence in the three-body constrained Bose-Hubbard model, demonstrated via exact diagonalisation, effective potential and quantum Monte Carlo simulations, motivated our hypothesis that it is also present in this modified Jaynes-Cummings-Hubbard system. Thus it would seem prudent to employ one or more of these to verify the existence of a pair superfluid of polaritons in such a coupled-cavity array.  ∼ ∼ ∼   a n a n a n n ( 0) 0, ( 1   where E k (n) represents the unperturbed energy eigenvalues. We note here that a solution to a 2 (n) = 0, for a given value of µ ′  , may only be utilised to describe the superfluid-Mott insulator phase boundary if, at κ =  0, the ground state is given by the lowest polariton branch corresponding to 〈 〉 = n n. Unphysical solutions such as Mott lobes seemingly extending to negative values of κ  can otherwise arise; these are necessarily discarded.