Effective spin physics in two-dimensional cavity QED arrays

We investigate a strongly correlated system of light and matter in two-dimensional cavity arrays. We formulate a Jaynes-Cummings Hamiltonian for two-level atoms coupled to cavity modes and driven by an external laser field which reduces to an effective spin Hamiltonian in the dispersive regime. In one dimension we provide exact analytical solution. In two dimensions, we perform mean-field study and large scale quantum Monte Carlo simulations of both the Jaynes-Cummings and the effective spin models. We discuss the phase diagram and the parameter regime which gives rise to frustrated interactions between the spins. We provide quantitative description of the phase transitions and correlation properties featured by the system and we discuss graph-theoretical properties of the ground states in terms of graph colorings using P\'{o}lya's enumeration theorem.


I. INTRODUCTION
Strongly coupled light-matter systems are at the heart of much of the effort in modern atomic and optical physics with applications ranging from quantum information processing to quantum simulations.
In this context, the use of cavities plays a prominent role as the strong confinement of the electromagnetic field implies strong interaction with matter coupled to the cavity modes. In particular, it offers possibilities to realize and study a plethora of quantum light-matter many-body Hamiltonians such as the so-called Jaynes-Cummings-Hubbard or Rabi-Hubbard models [1][2][3][4][5][6][7][8][9][10][11], or quantum fluids of light, where the effective interaction between light fields is mediated by a non-linear medium [12][13][14]. This offers ways to study various physical phenomena such as excitation propagation in chiral networks [15][16][17], the physics of spin glasses [18][19][20], self-organization of the atomic motion in optical cavities [21][22][23][24][25] or quantum phase transitions in arrays of nanocavity quantum dots [26] and in Coulomb crystals [27]. Furthermore, modern implementations of optical and microwave cavities allow for the creation of lattices with tunable geometries and dimensionality [28][29][30].
The paradigmatic description of cavity and circuit QED systems is typically in terms of the famous Jaynes-Cummings (JC) [31] or Dicke [32] models, which describe the interaction between the modes of the light field and the matter constituents, typically spin or phononic degrees of freedom of atoms or ions. Importantly, effective spin models emerge in the dispersive limit of the JC or Dicke models [33,34]. Under some circumstances this leads to spin Hamiltonians with frustrated or long-range interactions [35,36], thus offering ways to study rich physics of quantum magnetism. This is a particularly interesting direction allowing e.g. for studies of spin liquids [37][38][39] with optical quantum simulators.
While advanced numerical techniques, such as tensor networks, have been developed for spin Hamiltonians [40][41][42][43], the use of similarly efficient techniques for quantum optical systems, where a system of spins is coupled to the bosonic modes of an electromagnetic field remains limited. In this work we use mean-field (MF) description, exact diagonalization and large-scale quantum Monte Carlo (QMC) algorithm to study arrays of waveguide cavities (we note that in the context of cavity QED, QMC was implemented to study both the Jaynes-Cummings-Hubbard [44] and the Rabi-Hubbard [45] models). This gives us rigorous tools to investigate the emergence of the effective spin physics as a limiting case of the parent JC Hamiltonian for arbitrary lattice geometries and dimensions. Specifically, in this work we focus on square lattice geometry and we study the ground state properties of the JC and the effective spin models for various parameter regimes. We then show, that depending on the parameter regime and the array geometry, spin models with both non-frustrated and frustrated interactions can be engineered.
The paper is organized as follows. In Sec. II we introduce the system and derive the effective spin model from the parent JC Hamiltonian. In Sec. III we present exact analytical solution of the spin model in one dimension. In Sec. IV we present the results of the QMC simulation and MF analysis and discuss different regimes provided by the investigated model. We explain how the present work opens possibilities for simulating frustrated systems in Sec. V and conclude in Sec. VI. The atomic level scheme. The |1 − |e transition is driven by a classical field with Rabi frequency Ω and the excited state |e is adiabatically eliminated, see text for details. (c) Schematic of the emergent spin system after the elimination of the cavity modes. (d) Graphical representation of the parameter regime of the emergent spin Hamiltonian (5). The sign of the effective spin-spin interaction λ determines the nature of the spin-spin interaction: non-frustrated if λ < 0, frustrated if λ > 0. As |λ| is increased, a transition to a superradiant (SR) phase occurs, corresponding to a non-zero spin excitations of the system. While arbitrary rectangular arrays can be considered in the non-frustrated regime, only elongated geometries give rise to a non-trivial spin physics in the frustrated regime, cf. Sec. V for details.

A. Jaynes-Cummings Hamiltonian
Recent advances in integrated optical circuits, where in principle arbitrary waveguide geometries can be created with high precision by laser engraving in the silica substrate [46] and an active experimental effort to combine the waveguides with atomic microtraps on a single device [47,48] motivate us to investigate a system of three-level atoms embeded in waveguide cavities.
We consider a square cavity array, where we denote by a i (b ν ) the modes in the i-th row (ν-th column) of the array, as shown in Fig. 1a. We use the latin (greek) indices to denote the rows (columns) throughout the article. All sites of the array are occupied by identical three-level systems in a Λ configuration, where |0 , |1 denote the ground states and |e the excited state, see Fig. 1b. The cavity modes are coupled with strength g 0 to the |0 − |e transition, while the |1 − |e transition couples to a classical field Ω with frequency ω T which propagates perpendicularly to the plane of the array and which is detuned by ∆ e = ω 1 − ω T with respect to the |1 − |e transition. In the limit ∆ e g 0 , Ω one can eliminate the excited state, which we described in detail in our previous publication [49]. Furthermore, under the condition of strong driving the resulting Hamiltonian is of the JC type and reads (see Appendix A for the details of the derivation and of the full Hamiltonian) with the effective atomic transition frequency and the effective coupling strength Here, σ z,± iν are the usual Pauli matrices written in the {|1 , |0 } basis, ∆ i(ν) the effective cavity frequencies, ω at the effective frequency of the |0 − |1 transition and i = 1..L y , ν = 1..L x , where L y (L x ) is the number of rows (columns) of the array.

B. Spin Hamiltonian
In the large detuning limit one can further eliminate the cavity fields to obtain an effective spin Hamiltonian, where a given spin is coupled to all other spins belonging to the same cavity mode (see Fig.  1c and Appendix A), where where are the effective spin-spin couplings along the rows (columns) and C. Validity of the approximations and frustrated vs. non-frustrated regime The elimination of photonic or phononic fields giving rise to an effective spin physics is a known technique often used in the design of various quantum optical simulators. It can lead to interesting frustrated spin Hamiltonians, e.g. in the context of trapped ions [35,36,50].
In order to simplify the parameter space, in what follows we choose all the couplings to be the same along rows (columns): λ a ≡ λ i , ∀i (λ b ≡ λ ν , ∀ν). Schematically, the parameter regime of the spin Hamiltonian (5) is summarized in Fig. 1d.
First we note, that the parameters ω at /2, λ and δω at of the Hamiltonian (5) given by (3a),(7) and (8) can take both positive or negative sign. In particular the sign of λ determines the kind of physical situation provided by the interaction Hamiltonian (6c): non-frustrated if both couplings are negative, λ a(b) < 0 and frustrated otherwise. This is apparent from the form of the interaction which tends to align each pair of spins antiparallel whenever the corresponding coupling is positive. This then leads to frustration as the antiparallel alignment cannot be satisfied simultaneously for all the spins. Note that while we consider square lattice for concreteness, the presence of frustration in cavity arrays stems from all-to-all interaction between spins belonging to the same cavity mode and, hence, is independent of the lattice geometry.
Next, we discuss the parameter regimes of the Hamiltonian (5). We recall that the only requirement used in the derivation of (5) from the parent JC Hamiltonian (2) is the condition (4), g |ω at − ∆ a(b) |. (i) Weakly interacting regime. We refer to the weakly interacting regime as the regime where (we drop the a, b indices for simplicity) |λ| |ω at |.
Here, we have neglected the δω at term contributing to the atomic transition frequency since δω at ∝ λ [66]. One should verify that reaching the weakly interacting regime is compatible with the conditions (1), (4) used in the derivation of the Hamiltonians (2), (5). It is easy to show that it is indeed the case: substituting (7) for λ in (9), we get |g 2 /ω at | |∆ − ω at |. This implies, together with (4), that |g| |ω at |. Finally, substituting for g from (3b), we get |g 0 | |Ω|. This is enforced by the stronger condition (1), which completes our consistency check.
(ii) Strongly interacting regime. Here we refer to the strongly interacting regime as the regime where |λ| |ω at /2 + δω at |. Here, the cavity fields dependence of the δω at term plays an essential role. We leave this interesting possibility for Sec. V and focus first on the scenario where the dynamics of the cavity fields decouples from the spins leading to a pure spin Hamiltonian.

III. 1D: EXACT SOLUTION OF THE SPIN MODEL
It is illustrative to clarify on a simple example some of the basic properties of the Hamiltonian (5). Specifically, we are interested in the nature of phase transitions featured by (5) and the scalings of critical couplings. To this end we consider a one dimensional limit of (5) by taking a single cavity mode a. The Hamiltonian simplifies to Here, J are the total angular momentum operators We note that in the absence of the cavity fields, (10) is the well-known Lipkin-Meshkov-Glick model [51], which has been recently studied also in the context of cavity QED [52]. The advantage of the model (10) is that it is exactly solvable providing us with useful analytical insights. Using the usual angular momentum algebra where J is the half-integer total angular momentum (J = N/2 for N spins) and m = −J, −J + 1, ..., J, it follows that |J, m, n , where n is the photon number, are the eigenstates of the Hamiltonian (10). The eigenenergies are where in the second line we have regrouped the terms in order to emphasize the dependence on the photon number n.
The implications of the first bracket in the second line of (13) are the following. For the ground state photon number is 0. On the other hand, for E n < 0 the ground state photon number is n = ∞, which invalidates the approximate description in terms of the effective Hamiltonian (10). At this point it is important to note that since both ∆ − ω at in the denominator of λ and m can take positive or negative values, there is always a combination of m and ∆ − ω at where the transition n = 0 ↔ n = ∞ occurs as λ is varied. The situation is summarized in Table I. The main message contained in the Table I is that it is not possible to simulate the frustrated spin system using (2) in one dimension (see also [53]). In Sec. V we will show, how this limitation can be circumvented in two dimensions by exploiting the properties of the δω at term (8).
In what follows we shall investigate this transition and its relation to the parent JC Hamiltonian (2) further. The n = 0 ↔ n = ∞ transition occurs when E n changes sign. From (7) and (14) we get the expression for a critical coupling g c Lets first consider ω at > 0. In the non-frustrated case (λ < 0, ∆ > ω at ), the minimal possible value of g c corresponds to m = N/2 (i.e. all spins excited). On the other hand, for λ > 0 and positive ω at assumed here, we can have either ∆ > 0 or ∆ < 0. For ∆ < 0, we can see immediately from (14) that E n can be made always negative by a suitable choice of m. Specifically, considering the spin ground state m = −N/2, the global ground state would correspond to n = ∞ invalidating the description in terms of (10). On the other hand, for ∆ > 0 the system undergoes the n = 0 ↔ n = ∞ transition as λ is increased. However, it occurs for m = −N/2, i.e. before any spin transition could possibly take place. One could now proceed analogously for ω at < 0 [67]. In summary, the critical coupling at which the n = 0 ↔ n = ∞ transition occurs is given by where we have emphasized by the label "ph", that the transition is in the photon number.
In one dimension, the only non-trivial situation is thus the non-frustrated case, λ < 0, where n = 0. Here, a series of transitions between phases with N exc (m) and N exc (m + 1) = N exc (m) + 1 excited spins takes place as |λ| is increased (here N exc = (2m + N )/2). The corresponding coupling strengths at which these transitions occur are obtained from (13) by solving for E J,m,0 = E J,m+1,0 .
For instance, considering ω at > 0 and ∆ > ω at , third line in Table I, the first spin transition from m = −N/2 to m = −N/2 + 1 occurs at One can also read off from the expression (17) the scaling properties of the critical point of the spin transition with the system size, g c ∝ 1/ √ N and correspondingly for the critical coupling λ, |λ c | ∝ 1/N .

IV. 2D: ANALYTICAL STUDY AND QMC SIMULATIONS
After having analyzed the situation in 1D, we now turn our attention to 2D. It is well known that the JC model features second order superradiant phase transition as  (10). + (−) stand for positive (negative) values respectively. The "(n = ∞)" description in the second and fifth line indicates that in these cases, the ground state corresponds to the infinite photon number independently of the value of λ, see (14) and text for details.
the coupling strength is varied [32,54,55]. We will analyse the scaling properties at this superradiant phase transition and evaluate the two-point spin correlation functions of the cavity array. In order to do so, we employ large scale QMC simulations using the worm algorithm [56,57]. In the following we compare the QMC results with the MF solutions. We emphasize that in the considered square lattice geometry the spins are not all-to-all connected (they are connected only if they belong to the same row/column), i.e. it is not apriori obvious whether the MF solutions provide an accurate quantitative description.
In order to simplify the discussion, in this section we consider equal couplings along all rows and columns, λ ≡ λ a = λ b (i.e. ∆ ≡ ∆ i = ∆ ν , ∀i, ν). Motivated by the findings in the one-dimensional case, we focus only on the non-frustrated case with ω at > 0 and ∆ > ω at . We will address the frustrated case in Sec. V.
Mean-field solutions. In the thermodynamic limit, one can find MF solutions of the JC model (2) which we describe in detail in Appendix B and which we use for the sake of comparison with the QMC data. In particular, for an array of size N = L x ×L y one can find expressions for the critical strength g MF c of the coupling at which the superradiant transition occurs and the number of spin excitations N exc in the superradiant phase, which read and respectively. In the specific case of a square array L ≡ L x = L y and in the limit ∆ → ∞, where the descriptions in terms of (2) and (5) should coincide, we get with the help of (7) and The spin model (5) is an effective description of the parent JC model (2) in the limit of large detuning (4). Therefore, the excitations of the JC model in the superradiant phase result in spin excitations in the effective spin model. Here, QMC provides an efficient numerical tool to study this limit behaviour of the JC model and how well it is described by the effective spin model. The results of the simulations are presented in Fig. 2. Here we show the critical couplings of the superradiant phase transition g c determined using QMC (using the total number of photonic and spin excitations as order parameter, square data points) and the MF prediction (18), solid lines. The red (blue) data correspond to two different system sizes N = 18 × 18 (N = 28 × 28) respectively and the horizontal black lines are the values of the critical couplings λ c obtained from the QMC simulation of the effective spin model (5). As expected, we find that the values of the critical couplings approach asymptotically in the limit ∆ ω at where the two models (2) and (5) coincide. In the insets we show the finite size scaling of the critical couplings for the JC (left inset) and effective spin (right inset) models. As in the main plot, the squares represent the QMC data and the solid lines are the MF predictions (18) and (20). The slight departure of the scaling for the spin model for small system sizes is indeed a finite size effect on which we will comment momentarily. We also note, that the couplings for the present 2D model scale in the same way as the 1D predictions (17), i.e. in the linear extent of the system, g c ∝ 1/ √ L. This is due to the fact that the scaling is determined by the number of cavity modes to which the atoms couple rather than by the system size N (see also Appendix B). After we have verified, that the critical couplings of the JC model coincide with those of the spin model, we now study the number of excitations of the spin model as the coupling strength is varied. This is shown in Fig. 3. The solid line corresponds to the MF prediction (21). The squares correspond to the QMC simulation of the spin model (5) on a N = 20 × 20 array, where we have neglected λ a , λ b in (6a) as |λ| ω at in the studied regime. The discrepancy between the MF prediction and the QMC simulation is precisely the consequence of neglecting λ a , λ b in (6a) and results in an offset −1/(4L 2 ) in the values of λ -the dashed line corresponds to the MF solution corrected for this offset [68].
So far, we concentrated only on one-point observables in our QMC simulations and found a good agreement with the MF predictions. In order to go beyond the MF picture we next consider the correlation functions. Before presenting the results and in order to get a deeper insight in the structure of the spin Hamiltonians emergent in cavity arrays, in the following section we study the properties of the ground states from the group and graph theory perspective.
A. Ground state structure Symmetry considerations. We start our analysis in this section by noting that the total number of spin excitations,N exc = i,ν 1 2 (σ z iν + 1), is the constant of motion of the Hamiltonian (5). This significantly simplifies the problem in that in order to find the ground states of (5), one only needs to solve for the eigenstates of the interaction Hamiltonian (6c) in the given excitation number sector N exc . The problem can be further simplified by exploiting the real space symmetries of the Hamiltonian (5), similarly to the analysis carried out e.g. for the antiferromagnetic Heisenberg chain [58]. Considering the most symmetric situation, i.e. a square array with equal couplings (L x = L y , λ a = λ b ), the discrete symmetry group of the Hamiltonian (5) is G Lx×Lx = {R, C} ∪ D 4 , where R, C and D 4 stand for permutation of rows, permutation of columns and the dihedral group of the square array (i.e. successive rotations r π/2 , r π , r 3π/2 by π/2 and reflections about the horizontal (h), vertical (v) and the two diagonal(p, n) axes of the array) respectively. In order to get use of the symmetries, one has to find the irreducible representations (irreps) of G. While a systematic approach exists for finding irreps of the full symmetric group S N =LxLy [59], the subduced representations of the subgroup G ⊂ S N are in general reducible [60, ch.3]. Motivated by exact diagonalization results, instead of finding the irreps of G, we focus on the graph-theoretical properties of the ground states in what follows.
Let us start with the following observation based on the exact diagonalization results of (6c) in the non-frustrated case λ < 0 in the simplest non-trivial example, a plaquette (i.e. 2 × 2 array) with N exc = 2 excitations. The vertices of the plaquette are labeled 1-4, see Fig. 4. The ground state can be written as where which we symbolically write as where |s j is a specific spin configuration and |θ i | the number of such configurations belonging to a given set θ i . This seemingly artificial decomposition of the ground state into |θ 1 and |θ 2 is in fact directly related to the coloring of a graph as we now discuss. Let us start by introducing the notions necessary for our considerations which we then demonstrate on specific examples of the ground state construction. To this end we follow closely the treatment presented in [61].
• Lets consider a set S and a group G acting on S with ranks |S| and |G| respectively.
• For G a discrete group, each element g ∈ G can be written as a product of , where b j counts how many j−cycles appear in the decomposition of g. The product is thus a monomial representing the cycle structure of the element g • Lets consider m colours c 1 , ..., c m such that a specific colour c j is assigned to each element of S Definition 1: A colouring C is a specific configuration of colours on S.
For example, considering two colours (black and red), is a possible colouring of a plaquette.
Definition 2: An orbit of a colouring C is a set of all colourings produced by the action of the group G on C, orb G (C) = {g(C), g ∈ G} An orbit is thus an equivalence class of all colourings belonging to the orbit.
With the definitions above we now introduce two theorems: Theorem 1: Pólya's enumeration theorem [62]. Let C = {C} be a set of all colourings of S using colours c 1 , ..., c m . Let G induce an equivalence relation on C. Then is the generating function for the number of non-equivalent colorings of S in C.
Equipped with the necessary notions, we return back to the example of the ground state (23) [69]. In order to find the structure of the ground state corresponding to a given excitation number sector N exc , we need to enumerate the number of the sets θ and how many elements belong to each of the set. Here, we are concerned only with two colours, say black and red, which correspond to spins in ground and excited state respectively. In other words, θ is precisely an orbit and |θ| is thus given by the orbit-stabilizer theorem. We now demonstrate the use of the above theorems on our example of (23). The generating functional (27) In the second line, we have used the Pólya's theorem, where we have substituted the black (b) and red (r) colours, x j = b j + r j for j = 1, 2, 4. In our example of two excitations, i.e. the term with r 2 in (30), the numerical prefactor 2 means there are two equivalence classes θ 1 , θ 2 of the colourings. These can be found explicitly and read where e stands for the identity element of the group G. Finally, one can verify that the above relations obey the orbit-stabilizer theorem (29) so that |θ 1 | = 4 and |θ 2 | = 2 with the states written explicitly in (25).
The above results can be generalized straightforwardly to larger arrays. In that case the ground state can be written as where the orbit states |θ i are orthonormal, θ i |θ j = δ ij . To give an explicit example going beyond the plaquette, we consider a 3 × 3 array and we choose N exc = 4 sector. We find that the total of 9 4 = 126 spin basis states form five equivalence classes depicted in Fig. 4 The Hamiltonian in the orbit states basis {|θ 1 , |θ 2 |θ 3 , |θ 4 , |θ 5 } reads We note that θ 2 | H spin,int |θ 2 = θ 5 | H spin,int |θ 5 = 0. This can be easily understood when inspecting the structure of θ 2 , θ 5 and realizing that the action of the Hamiltonian (6c) is to anihilate an excitation at a given site and create an excitation at another site. It can be seen from Fig. 4 that such operations necessarily take a state belonging to θ 2 or θ 5 out of its equivalence class. For instance for the example of θ 2 in Fig. 4, anihilating the excitation at position 5 and creating an excitation at position 6 results in the state θ 3 shown and the corresponding non-zero matrix element θ 2 | H spin,int |θ 3 in (33) (in fact an operation displacing two excitations at once would be required for the state to remain in θ 2 ). While the above considerations offer a useful insight into the structure of the ground state, so far they do not constitute a clear computational advantage as we did not provide a prescription for obtaining the Hamiltonian (6c) in the {|θ i } basis. Such prescription is likely to be equivalent to finding the irreps of G as discussed at the beginning of this section. We leave this investigation for future work and restrict ourselves only to exact diagonalization in the comparative QMC study of the correlations presented in the following section.

B. Correlations
We are now in position to study the correlation functions of the spin model. To this end we consider (connected) two-point correlations of the type σ + iν σ − jµ . Due to long (infinite) range connectivity along the rows and the columns, the system features only two length scales, the nearest-neighbours (NN, spins belonging to the same cavity mode) and next-to-nearest-neighbors (NNN) as there are at most two different cavity modes connecting any two spins of the array. We thus define two types of correlation functions, where we have excluded self-correlations of the type σ + σ − = |1 1| . This situation is schematically depicted in the inset of Fig. 5a. Here, Σ NN corresponds to correlations between the spin in the green box and the spins belonging to the same row and column (red-shaded region). Similarly, Σ NNN corresponds to correlations between the spin in the green box and the spins belonging to the blue-shaded region. In Fig. 5 we plot the ratio Σ NNN /Σ NN as a function of the number of spin excitations N exc in an N = 3 × 3 (Fig. 5a) and N = 5 × 5 array (Fig. 5b). For the 3 × 3 array we find perfect agreement between the exact results obtained by exact diagonalization of the spin Hamiltonian (6c) in each excitation sector and the QMC simulation of that Hamiltonian [70]. Moreover we find a good agreement also with the QMC simulation of the JC model (2), which is improving with increasing value of ∆/ω at as it should. Similar agreement between the QMC simulations of the spin and the JC model is observed for the 5 × 5 array In (a) we find good agreement between the exact diagonalization results of the spin model (blue circles), the QMC simulation of the spin model (red squares) and the QMC simulation of the JC model (green and magenta triangles for ∆/ω at = 20 and 40 respectively). (b) Similar agreement is obtained for the 5 × 5 array. The inset in (a) is a schematic representation of the connectivity in the cavity array: NN is represented by the red-shaded (NNN by the blue-shaded) regions, see text for details.

V. TOWARDS SIMULATION OF FRUSTRATED SPIN SYSTEMS IN CAVITY ARRAYS
In Sec. III we have shown, that in one dimension it is not possible to obtain the effective spin Hamiltonian (5) with frustrated interactions λ > 0. The aim of this section is to show that this limitation can be circumvented in two dimensions by exploiting the properties of the δω term (8).
In order to see this, we first perform a back-of-the-envelope estimation. The reason of the breakdown of the effective model in the frustrated case in one dimension is that as λ is decreased, the term λa † a in (10) is decreasing in a way that the n = 0 ↔ n = ∞ transition occurs before any spin excitation can appear.
In order to simplify the analytical treatment, from now on we assume all the couplings along rows to be the same, λ a ≡ λ i , ∀i and similarly for the columns, λ b ≡ λ ν , ∀ν. Assuming further that the row (column) cavity photon occupation numbers are n a (n b ), the expectation value of the δω term (8) becomes It is apparent from the above expression that the magnitude of δω can be made small when the couplings along rows and columns have opposite signs, λ a = −λ b , so that the amplitude of each bracket in (34) becomes significantly smaller than if we take both λ with the same sign. We will thus consider a scenario with frustrated interactions along one direction (we choose the a cavity modes) and non-frustrated one along the other (b modes) and we paramterize the couplings as In order to show that one can get a non-trivial frustrated-non-frustrated situation without breaking the effective spin description, we will use a self-consistency argument as follows. First we restore the free cavity fields terms we omitted in (5) so that the effective spin Hamiltonian can be written as (see Appendix A) where Ly µ=1 δω at,jµ σ z jµ (37) is now the photonic part (note that we have absorbed the δH term (6b)) in H ph and H spin,int is given by (6c). Exploiting the fact that the σ z term is diagonal in eigenstate basis in all excitation sectors N exc of the spin Hamiltonian (5), we substitute the spin expectation values σ z in (37) and subsequently diagonalize the photonic Hamiltonian, which is a straightforward exercise as it is quadratic in the photonic degrees of freedom. We then compare the values of the critical couplings at which a transition N exc → N exc + 1 occurs with that of 0 → ∞ photon number in analogy to the analysis in Sec. III. We anticipate that a non-trivial frustrated regime can be always obtained by appropriate tuning of the system parameters and in particular its geometry. This is also the regime which fulfills the self-consistency criterion, namely taking δω = 0 in the spin model, using the corresponding solutions in the photonic Hamiltonian (37) and finding that its solutions again yield δω = 0.

A. Diagonalization of the spin Hamiltonian
In analogy to Sec. IV we seek to diagonalize the interaction part of the spin Hamiltonian (6c) where we have now explicitly written the summation limits.

1-excitation sector
In the basis {|100..0 , |010..0 , ..., |000..1 } of single particle states, the interaction Hamiltonian (6c) takes a simple form where the a (b) matrices have dimensions L x × L x (L y × L y ) respectively and the prefactor 2 comes from accounting twice for each spin configuration in (6c). M are matrices with 1 everywhere except the diagonal, where it is 0, M ij = 1 − δ ij . The corresponding eigenvalues and multiplicities are matrix eigenvalue multiplicity Since λ a > 0 and λ b < 0, the minimum energy is with multiplicity L x − 1. The corresponding eigenvectors are where |1 j ≡ |0..01 j 0..0 . The ground state eigenvectors of H 1exc are then |ψ 0 j = |v j a ⊗ |v b and the spin expectation values become if |ψ 0 j contains the excitation at site lµ or -1 otherwise. We note that s z lµ → −1 in the thermodynamic limit as one would expect.
In summary, we have for the ground state energies where

B. Diagonalization of the photonic Hamiltonian
We are now in position to diagonalize the photonic quadratic form (37). In analogy to the previous paragraph, we start our examination in the 0-excitation sector. Here the expectation values of the spin operator is simply s z ≡ s z lµ = −1, ∀l, µ, so that the photonic Hamiltonian becomes where we have introduced p = (a 1 , .., a Ly , b 1 , .., b Lx ) T . The matrix M ph can be written as M ph has the following eigenvalues and multiplicities C. Validity and breakdown of the effective spin model First we note, that due to (35), ∆ b is not independent and can be expressed as (here ω at is the bare atomic frequency, not ω at ). Motivated by the condition (4) needed for the spin model (5) to be valid, we define what we call the quality factor of the approximation as where g spin c is given by the critical value of λ a , g spin c = −2λ a,crit (∆ a − ω at ), cf. below. We start by determining the critical value of λ a at which a transition from 0-to 1-excitation sector occurs in the spin model. This can be simply obtained from the condition E 0exc GS = E 1exc GS and with the help of (44) we get Next, the breakdown of the effective model is indicated when any of the photonic eigenvalues (51) becomes negative, corresponding to infinitely many photons in the ground state. Substituting s z = −1 in (51), the only two candidates for the minimum eigenvalue are E ph a and E ph − (we recall that λ a > 0, λ b < 0). The value of λ a where E ph becomes negative is determined as where only the positive branch of λ a in the solutions of E ph − (λ a ) = 0 is considered. The criterion of having a valid and non-trivial regime in the effective spin Hamiltonian (i.e. finite number of photons and non-zero spin excitations) thus translates into In Fig. 6 we plot the contours of constant R (left pane) and Q (right pane) respectively in the η-L x /L y plane. It is apparent from the figure that increasing both L y /L x and |η| leads to a larger critical ratio R, i.e. we can ensure the presence of the non-trivial region by tuning these parameters. Additionally, one can show analytically that when keeping all the other parameters fixed, as expected from the contours in Fig. 6 (here R asymptotic is some finite asymptotic value).
So far we have concetrated only on the simple 0 and 1-excitation sectors of the spin Hamiltonian. Clearly, for the interactions to become relevant one is interested in sectors with larger number of excitations. To this end one could extend the above analysis to the N exc → N exc + 1 particle transitions for N exc ≥ 1 and for each of the transitions evaluate the critical ratio R. The fully analytical approach is unfortunatelly obscured by the fact that for particle sectors N exc > 1, the spin interaction matrix (6c) does not take the simple structure of (38) and the evaluation of eigenvalues thus amounts to solving higher order polynomials which in general requires a numerical approach. However, the arbitrarily large values of R for the 0 → 1 spin excitation transition are suggestive and it is likely that also higher excitation number sectors can be reached without breaking the validity of the spin Hamiltonian. To recap, we have shown that the spin model (5) with both frustrated and nonfrustrated interactions can emerge as an effective description of the parent JC Hamiltonian (2). This can be achieved when considering an elongated geometry of the square array, L y L x . On one hand this circumvents the limitations related to realizing effective spin Hamiltonians with frustrated interactions using optical setups governed by JC Hamiltonians [53]. Finally, we note that one can simulate the parent JC model in the regime where it yields the effective spin model using QMC avoiding thus a sign problem of the spin model which opens up an interesting perspective on the QMC simulations of Hamiltonians with sign problem.

VI. CONCLUSIONS AND OUTLOOK
In this work we have analysed the ground states of a cavity array where each intersection of cavity modes is occupied by a single atom. We have shown that the system's description in terms of the JC model leads to an effective description -in a suitable parameter regime, where the cavity modes can be adiabatically eliminated -in terms of a two-component spin model. In one dimension, we have provided exact solution of the spin model demonstrating explicitly the need of higher dimensions in order to obtain frustrated spin-spin interactions. Using large-scale QMC simulation of the JC model we have performed a quantitative comparison between the parent JC and the emerging spin model. Specifically, in two dimensions, we have studied the superradiant phase transition and the properties of two-point correlation functions in the cavity array and we have described the graph-theoretical structure of the ground states of the spin model. In all cases we found a firm agreement between the two models in the regime of validity of the approximations used. Finally, we have outlined the possibility, by exploiting the non-linearities of the effective spin model, of studying frustrated spin Hamiltonians using two-dimensional cavity arrays.
In conclusion, the theoretical framework and numerical tools established in this work open ways to address the cavity QED physics in a quantitative way beyond the traditional mean-field or perturbative treatments. The present developments can be exploited in various scenarios, such as the study of the ground state properties in different lattice geometries and dimensions.