Effective spin physics in two-dimensional cavity QED arrays

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 [50] and an active experimental effort to combine the waveguides with atomic microtraps on a single device [51, 52] 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 ith row (νth column) of the array, as shown in figure 1(a). 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\rangle ,| 1\rangle$ denote the ground states and $| e\rangle$ the excited state, see figure 1(b). The cavity modes are coupled with strength g₀ to the $| 0\rangle -| e\rangle$ transition, while the $| 1\rangle -| e\rangle$ transition couples to a classical field Ω with frequency ${\omega }_{T}$ which propagates perpendicularly to the plane of the array and which is detuned by ${{\rm{\Delta }}}_{e}={\omega }_{1}-{\omega }_{T}$ with respect to the $| 1\rangle -| e\rangle$ transition. In the limit ${{\rm{\Delta }}}_{e}\gg {g}_{0},{\rm{\Omega }}$ one can eliminate the excited state, which we described in detail in our previous publication [53]. Furthermore, under the condition of strong driving

$$\begin{eqnarray}&&{\rm{\Omega }}\gg {g}_{0}\end{eqnarray} \tag{ 1 }$$

the resulting Hamiltonian is of the TC type and reads (see appendix A for the details of the derivation and of the full Hamiltonian)

$$\begin{eqnarray}&&{H}_{\mathrm{TC}}=\displaystyle \sum _{i}{{\rm{\Delta }}}_{i}{a}_{i}^{\dagger }{a}_{i}+\displaystyle \sum _{\nu }{{\rm{\Delta }}}_{\nu }{b}_{\nu }^{\dagger }{b}_{\nu }+\displaystyle \frac{{\omega }_{\mathrm{at}}}{2}\displaystyle \sum _{i,\nu }{\sigma }_{i\nu }^{z}+g\displaystyle \sum _{i,\nu }({\sigma }_{i\nu }^{+}{a}_{i}+{a}_{i}^{\dagger }{\sigma }_{i\nu }^{-})+g\displaystyle \sum _{i,\nu }({\sigma }_{i\nu }^{+}{b}_{\nu }+{b}_{\nu }^{\dagger }{\sigma }_{i\nu }^{-}).\end{eqnarray} \tag{ 2 }$$

with the effective atomic transition frequency and the effective coupling strength

$$\begin{eqnarray}&&{\omega }_{\mathrm{at}}=-\displaystyle \frac{{{\rm{\Omega }}}^{2}}{{{\rm{\Delta }}}_{e}},\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&g=-\displaystyle \frac{{g}_{0}{\rm{\Omega }}}{{{\rm{\Delta }}}_{e}}.\end{eqnarray} \tag{ 3b }$$

Here, ${\sigma }_{i\nu }^{z,\pm }$ are the usual Pauli matrices written in the $\{| 1\rangle ,| 0\rangle \}$ basis, ${{\rm{\Delta }}}_{i(\nu )}$ the effective cavity frequencies, ${\omega }_{\mathrm{at}}$ the effective frequency of the $| 0\rangle -| 1\rangle$ transition and $i=1..{L}_{y},\nu =1..{L}_{x}$, where ${L}_{y}({L}_{x})$ is the number of rows (columns) of the array.

Zoom In Zoom Out Reset image size

In the large detuning (dispersive) limit

$$\begin{eqnarray}&&g\ll | {\omega }_{\mathrm{at}}-{{\rm{\Delta }}}_{i(\nu )}| ,\end{eqnarray} \tag{ 4 }$$

one can further perform a unitary transformation to perturbatively eliminate the cavity fields in order to obtain an effective spin Hamiltonian, where a given spin is coupled to all other spins belonging to the same cavity mode (see figure 1(c) and appendix A),

$$\begin{eqnarray}&&{H}_{\mathrm{spin}}={H}_{\mathrm{spin},0}+\delta H+{H}_{\mathrm{spin},\mathrm{int}},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{H}_{\mathrm{spin},0}=\displaystyle \sum _{i,\nu }\left(\displaystyle \frac{{\omega }_{\mathrm{at}}}{2}+{\lambda }_{i}+{\lambda }_{\nu }\right){\sigma }_{i\nu }^{z},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&\delta H=\displaystyle \sum _{i,\nu }\delta {\omega }_{\mathrm{at},i\nu }{\sigma }_{i\nu }^{z},\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&{H}_{\mathrm{spin},\mathrm{int}}=\displaystyle \sum _{i,\nu \ne \mu }{\lambda }_{i}({\sigma }_{i\nu }^{-}{\sigma }_{i\mu }^{+}+{\sigma }_{i\nu }^{+}{\sigma }_{i\mu }^{-})+\displaystyle \sum _{i\ne j,\nu }{\lambda }_{\nu }({\sigma }_{i\nu }^{-}{\sigma }_{j\nu }^{+}+{\sigma }_{i\nu }^{+}{\sigma }_{j\nu }^{-}),\end{eqnarray} \tag{ 6c }$$

where

$$\begin{eqnarray}&&{\lambda }_{i}=-\displaystyle \frac{{g}^{2}}{2({{\rm{\Delta }}}_{i}-{\omega }_{\mathrm{at}})},\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{\lambda }_{\nu }=-\displaystyle \frac{{g}^{2}}{2({{\rm{\Delta }}}_{\nu }-{\omega }_{\mathrm{at}})}\end{eqnarray} \tag{ 7b }$$

are the effective spin–spin couplings along the rows (columns) and

$$\begin{eqnarray}&&\delta {\omega }_{\mathrm{at},i\nu }={\lambda }_{i}(2{a}_{i}^{\dagger }{a}_{i}+{a}_{i}^{\dagger }{b}_{\nu }+{b}_{\nu }^{\dagger }{a}_{i})+{\lambda }_{\nu }(2{b}_{\nu }^{\dagger }{b}_{\nu }+{a}_{i}^{\dagger }{b}_{\nu }+{b}_{\nu }^{\dagger }{a}_{i}).\end{eqnarray} \tag{ 8 }$$

Dispersively eliminating photonic or phononic fields leading 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 [39, 40, 54].

In order to simplify the parameter space, in what follows we choose all the couplings to be the same along rows (columns): ${\lambda }_{a}\equiv {\lambda }_{i},\,\forall i$ (${\lambda }_{b}\equiv {\lambda }_{\nu },\,\forall \nu$). Schematically, the parameter regime of the spin Hamiltonian (5) is summarized in figure 1(d).

First we note, that the parameters ${\omega }_{\mathrm{at}}/2$, λ and $\delta {\omega }_{\mathrm{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, ${\lambda }_{a(b)}\lt 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 TC Hamiltonian (2) is the condition (4), $g\ll | {\omega }_{\mathrm{at}}-{{\rm{\Delta }}}_{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)

$$\begin{eqnarray}&&| {\lambda }_{a(b)}| \ll | {\omega }_{\mathrm{at}}| .\end{eqnarray} \tag{ 9 }$$

Here, we have neglected the $\delta {\omega }_{\mathrm{at}}$ term contributing to the atomic transition frequency since $\delta {\omega }_{\mathrm{at}}\propto \lambda$ ⁴ . One should verify that reaching the weakly interacting regime is compatible with the conditions (1) and (4) used in the derivation of the Hamiltonians (2) and (5). It is easy to show that it is indeed the case: substituting (7) for λ in (9), we get $| {g}^{2}/{\omega }_{\mathrm{at}}| \ll | {\rm{\Delta }}-{\omega }_{\mathrm{at}}|$. This implies, together with (4), that $| g| \lesssim | {\omega }_{\mathrm{at}}|$. Finally, substituting for g from (3b), we get $| {g}_{0}| \lesssim | {\rm{\Omega }}|$. This is enforced by the stronger condition (1), which completes our consistency check.

(ii) Strongly interacting regime. We refer to the strongly interacting regime as the regime where $| {\lambda }_{a(b)}| \gtrsim | {\omega }_{\mathrm{at}}/2+\langle \delta {\omega }_{\mathrm{at}}\rangle |$. Here, the cavity fields dependence of the $\delta {\omega }_{\mathrm{at}}$ term plays an essential role. We leave this interesting possibility for section 5 and focus first on the scenario where the dynamics of the cavity fields decouples from the spins leading to a pure spin Hamiltonian.

While cavity QED has become a well established experimental research direction, a brief discussion of whether the ground state physics studied in this paper can be accessed in a realistic experiment is in order. We note that since the excited state $| e\rangle$ has been eliminated, the atoms are assumed to be trapped at the positions of the lattice sites and the levels $| 0\rangle ,| 1\rangle$ correspond to the long-lived hyperfine states, the main decoherence mechanism is due to the decay κ of the cavity modes (we consider the same decay rate for all modes). In order to estimate whether one can reach the ground state of the spin Hamiltonian (5), we consider an adiabatic preparation scheme, i.e. a situation where the parameters of the Hamiltonian are tuned sufficiently slowly so that the state of the system at a given time is also its ground state [55]. Specifically, the rate of change of the Hamiltonian parameters r should be much smaller than the energy gap ${\rm{\Delta }}E$ between the ground and the first excited state in order to avoid level mixing, $r\ll {\rm{\Delta }}E$. Provided ${\omega }_{\mathrm{at}}$ is the largest energy scale in (5), the gap corresponds to the energy cost of exciting a single spin, i.e. ${\rm{\Delta }}E\approx {\omega }_{\mathrm{at}}$. At the same time, r should be much faster than the cavity decay in order to be able to reach the desired ground state without being affected by the decay, which yields the following condition for the timescales

$$\begin{eqnarray}&&| {\omega }_{\mathrm{at}}| \gg r\gg \kappa \,\Rightarrow \,| {\omega }_{\mathrm{at}}| \gg \kappa .\end{eqnarray} \tag{ 10 }$$

Using (3a) and the fact that ${{\rm{\Delta }}}_{e}\gg {\rm{\Omega }}\gg {g}_{0}$ (see equation (1) and the discussion above it) we arrive at $| {\omega }_{\mathrm{at}}| \gg {g}_{0}$ provided ${\rm{\Omega }}/{g}_{0}\gg {{\rm{\Delta }}}_{e}/{\rm{\Omega }}$. Therefore one can see, that the condition (10) is automatically fulfilled in the strong coupling regime, where ${g}_{0}\gtrsim \kappa$. This regime has been realized in a number of platforms, such as optical microcavities [56], fibre-based cavity on atomic chip [57] or open-access on-chip microcavities [52] all combined with ⁸⁷Rb atoms with $({g}_{0},\kappa )\approx 2\pi \times (25,2.5)$ MHz, $2\pi \times (250,50)$ MHz and $2\pi \times (1,6.5)$ GHz for [56, 57] and [52] respectively.

Other elements required to realize the proposed setup have been also achieved experimentally, such as a creation of optical lattices with unit filling using optical tweezer arrays [58, 59] or in principal arbitrary geometries of silica-engraved waveguide arrays [30–32]. While those elements, together with the condition (10), are within reach of current technology, their combination on a single device remains a challenge. Fortunately, there are active experimental efforts to achieve this goal, for example to combine the waveguide arrays with atomic microtraps [51] or atoms with photonic-crystal nanocavities [60].

Symmetry considerations. We start our analysis in this section by noting that the total number of spin excitations, ${\hat{N}}_{\mathrm{exc}}={\sum }_{i,\nu }\tfrac{1}{2}({\sigma }_{i\nu }^{z}+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)

$$\begin{eqnarray}&&{H}_{\mathrm{spin},\mathrm{int}}| {\psi }_{\mathrm{GS}}\rangle ={E}_{\mathrm{GS}}| {\psi }_{\mathrm{GS}}\rangle \end{eqnarray} \tag{ 23 }$$

in the given excitation number sector ${N}_{\mathrm{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 [68]. Considering the most symmetric situation, i.e. a square array with equal couplings (L_x = L_y, ${\lambda }_{a}={\lambda }_{b}$), the discrete symmetry group of the Hamiltonian (5) is ${G}_{{L}_{x}\times {L}_{x}}=\{{\mathfrak{R}},{\mathfrak{C}}\}\cup {D}_{4}$, where ${\mathfrak{R}},{\mathfrak{C}}$ and D₄ stand for permutation of rows, permutation of columns and the dihedral group of the square array (i.e. successive rotations ${r}_{\pi /2},{r}_{\pi },{r}_{3\pi /2}$ by $\pi /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={L}_{x}{L}_{y}}$ [69], the subduced representations of the subgroup $G\subset {S}_{N}$ are in general reducible [70, chapter 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 $\lambda \lt 0$ in the simplest non-trivial example, a plaquette (i.e. 2 × 2 array) with ${N}_{\mathrm{exc}}=2$ excitations. The vertices of the plaquette are labeled 1–4, see figure 4. The ground state can be written as

$$\begin{eqnarray}&&| {\psi }_{\mathrm{GS}}\rangle =\displaystyle \frac{1}{\sqrt{2}}(| {\theta }_{1}\rangle +| {\theta }_{2}\rangle ),\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&| {\theta }_{1}\rangle =\displaystyle \frac{1}{\sqrt{4}}(| 1100\rangle +| 1010\rangle +| 0101\rangle +| 0011\rangle ),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&| {\theta }_{2}\rangle =\displaystyle \frac{1}{\sqrt{2}}(| 1001\rangle +| 0110\rangle ),\end{eqnarray} \tag{ 26 }$$

which we symbolically write as

$$\begin{eqnarray}&&| {\theta }_{i}\rangle =\displaystyle \frac{1}{\sqrt{| {\theta }_{i}| }}\displaystyle \sum _{j\in {\theta }_{i}}| {s}_{j}\rangle ,\end{eqnarray} \tag{ 27 }$$

where $| {s}_{j}\rangle$ is a specific spin configuration and $| {\theta }_{i}|$ the number of such configurations belonging to a given set ${\theta }_{i}$. This seemingly artificial decomposition of the ground state into $| {\theta }_{1}\rangle$ and $| {\theta }_{2}\rangle$ is in fact directly related to the colouring of a graph as we now discuss.

Zoom In Zoom Out Reset image size

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

• 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\in G$ can be written as a product of j-cycles x_j, $g\to {({x}_{1}^{{b}_{1}}{x}_{2}^{{b}_{2}}...{x}_{| S| }^{{b}_{| S| }})}_{g}$, where b_j counts how many j-cycles appear in the decomposition of g. The product ${({x}_{1}^{{b}_{1}}{x}_{2}^{{b}_{2}}...{x}_{| S| }^{{b}_{| S| }})}_{g}$ is thus a monomial representing the cycle structure of the element g
• Lets consider m colours ${c}_{1},\ldots ,{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, ${\mathrm{orb}}_{G}(C)=\{g(C),g\in G\}$

An orbit is thus an equivalence class of all colourings belonging to the orbit.

Definition 3. A stabilizer of a colouring C is a set of all group elements g which leave C invariant, ${\mathrm{stab}}_{G}(C)=\{g\in G,g(C)=C\}$

Definition 4. A generating function (or pattern inventory or cycle index) is a polynomial given by the sum of all monomials of elements of G acting on S

$$\begin{eqnarray}&&{P}_{G}({x}_{1},\ldots {x}_{| S| })=\displaystyle \frac{1}{| G| }\displaystyle \sum _{g\in G}{({x}_{1}^{{b}_{1}}{x}_{2}^{{b}_{2}}...{x}_{| S| }^{{b}_{| S| }})}_{g}.\end{eqnarray} \tag{ 28 }$$

With the definitions above we now introduce two theorems:

Theorem 1. Pólya's enumeration theorem [72]. Let ${ \mathcal C }=\{C\}$ be a set of all colourings of $S$ using colours ${c}_{1},\ldots ,{c}_{m}$. Let $G$ induce an equivalence relation on ${ \mathcal C }$. Then

$$\begin{eqnarray}&&{P}_{G}\left(\displaystyle \sum _{i=1}^{m}{c}_{i},\displaystyle \sum _{i=1}^{m}{c}_{i}^{2},\ldots ,\displaystyle \sum _{i=1}^{m}{c}_{i}^{| S| }\right)\end{eqnarray} \tag{ 29 }$$

is the generating function for the number of non-equivalent colourings of $S$ in ${ \mathcal C }$.

Theorem 2. Orbit-stabilizer theorem [73, 74, chapter 7].

$$\begin{eqnarray}&&| {\mathrm{orb}}_{G}(C)| =\displaystyle \frac{| G| }{| {\mathrm{stab}}_{G}(C)| }.\end{eqnarray} \tag{ 30 }$$

Equipped with the necessary notions, we return back to the example of the ground state (24)⁸ . In order to find the structure of the ground state corresponding to a given excitation number sector ${N}_{\mathrm{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 $| \theta |$ is thus given by the orbit-stabilizer theorem. We now demonstrate the use of the above theorems on our example of (24). The generating functional (28) of the ${G}_{2\times 2}=\{{\mathfrak{R}},{\mathfrak{C}}\}\cup {D}_{4}={D}_{4}$ group of the plaquette reads (see footnote 7)

$$\begin{eqnarray}{P}_{{G}_{2\times 2}} & = & \displaystyle \frac{1}{8}({x}_{1}^{4}+2{x}_{1}^{2}{x}_{2}+3{x}_{2}^{2}+2{x}_{4})\\ & = & {b}^{4}+{b}^{3}r+2{b}^{2}{r}^{2}+{{br}}^{3}+{r}^{4}.\end{eqnarray} \tag{ 31 }$$

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² in (31), the numerical prefactor 2 means there are two equivalence classes ${\theta }_{1},{\theta }_{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 (30) so that $| {\theta }_{1}| =4$ and $| {\theta }_{2}| =2$ with the states written explicitly in (26).

The above results can be generalised straightforwardly to larger arrays. In that case the ground state can be written as

$$\begin{eqnarray}&&| {\psi }_{\mathrm{GS}}\rangle =\displaystyle \sum _{i}{\psi }_{i}| {\theta }_{i}\rangle ,\end{eqnarray} \tag{ 33 }$$

where the orbit states $| {\theta }_{i}\rangle$ are orthonormal, $\langle {\theta }_{i}| {\theta }_{j}\rangle ={\delta }_{{ij}}$. To give an explicit example going beyond the plaquette, we consider a 3 × 3 array and we choose ${N}_{\mathrm{exc}}=4$ sector. We find that the total of $\left(\genfrac{}{}{0em}{}{9}{4}\right)=126$ spin basis states form five equivalence classes depicted in figure 4.

The Hamiltonian in the orbit states basis $\{| {\theta }_{1}\rangle ,| {\theta }_{2}\rangle \,| {\theta }_{3}\rangle ,| {\theta }_{4}\rangle ,| {\theta }_{5}\rangle \}$ reads

$$\begin{eqnarray}{H}_{\mathrm{spin},\mathrm{int}}=\left(\begin{array}{ccccc}3 & 0 & 2 & 2 & 2\\ 0 & 0 & 4 & 0 & 0\\ 2 & 4 & 2 & 4 & 0\\ 2 & 0 & 4 & 4 & 4\\ 2 & 0 & 0 & 4 & 0\end{array}\right).\end{eqnarray} \tag{ 34 }$$

We note that $\langle {\theta }_{2}| {H}_{\mathrm{spin},\mathrm{int}}| {\theta }_{2}\rangle =\langle {\theta }_{5}| {H}_{\mathrm{spin},\mathrm{int}}| {\theta }_{5}\rangle =0$. This can be easily understood when inspecting the structure of ${\theta }_{2},{\theta }_{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 figure 4 that such operations necessarily take a state belonging to ${\theta }_{2}$ or ${\theta }_{5}$ out of its equivalence class. For instance for the example of ${\theta }_{2}$ in figure 4, anihilating the excitation at position 5 and creating an excitation at position 6 results in the state ${\theta }_{3}$ shown and the corresponding non-zero matrix element $\langle {\theta }_{2}| {H}_{\mathrm{spin},\mathrm{int}}| {\theta }_{3}\rangle$ in (34) (in fact an operation displacing two excitations at once would be required for the state to remain in ${\theta }_{2}$).

The main result of the present discussion is that the set of states (27), which has a clear group-theoretical interpretation in terms of equivalence classes of a coloured graph, constitutes a natural basis for the ground state of the system. While it provides a useful insight into the structure of the ground state, it does not represent a clear computational advantage as we did not provide a prescription for obtaining the Hamiltonian (6c) in the $\{| {\theta }_{i}\rangle \}$ 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.

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 $\langle {\sigma }_{i\nu }^{+}{\sigma }_{j\mu }^{-}\rangle$. Due to long (infinite) range connectivity along the rows and the columns, the system features only two length scales corresponding to spins belonging to the same cavity mode (intra-cavity spins, IC) and to different cavity modes (extra-cavity spins, EC) respectively as there are at most two different cavity modes connecting any two spins of the array. We thus define two types of correlation functions, ${{\rm{\Sigma }}}_{\mathrm{IC}}\equiv \{\langle {\sigma }_{i\nu }^{+}{\sigma }_{i\mu }^{-}\rangle ,\nu \ne \mu \}\cup \{\langle {\sigma }_{i\nu }^{+}{\sigma }_{j\nu }^{-}\rangle ,i\ne j\}$ and ${{\rm{\Sigma }}}_{\mathrm{EC}}\equiv \{\langle {\sigma }_{i\nu }^{+}{\sigma }_{j\mu }^{-}\rangle ,i\ne j,\nu \ne \mu \}$, where we have excluded self-correlations of the type $\langle {\sigma }^{+}{\sigma }^{-}\rangle =\langle | 1\rangle \langle 1| \rangle$. This situation is schematically depicted in the inset of figure 5(a). Here, ${{\rm{\Sigma }}}_{\mathrm{IC}}$ corresponds to correlations between the spin in the green box and the spins belonging to the same row and column (red-shaded region). Similarly, ${{\rm{\Sigma }}}_{\mathrm{EC}}$ corresponds to correlations between the spin in the green box and the spins belonging to the blue-shaded region. In figure 5 we plot the ratio ${{\rm{\Sigma }}}_{\mathrm{EC}}/{{\rm{\Sigma }}}_{\mathrm{IC}}$ as a function of the number of spin excitations ${N}_{\mathrm{exc}}$ in an $N=3\times 3$ (figure 5(a)) and $N=5\times 5$ array (figure 5(b)). 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⁹ . Moreover we find a good agreement also with the QMC simulation of the TC model (2), which is improving with increasing value of ${\rm{\Delta }}/{\omega }_{\mathrm{at}}$ as it should. Similar agreement between the QMC simulations of the spin and the TC model is observed for the 5 × 5 array.

Zoom In Zoom Out Reset image size

In analogy to section 4 we seek to diagonalize the interaction part of the spin Hamiltonian (6c)

$$\begin{eqnarray*}&&H={\lambda }_{a}\displaystyle \sum _{j=1}^{{L}_{x}}\displaystyle \sum _{\mu \ne \nu =1}^{{L}_{y}}{\sigma }_{j\mu }^{+}{\sigma }_{j\nu }^{-}+{\lambda }_{b}\displaystyle \sum _{\alpha =1}^{{L}_{y}}\displaystyle \sum _{k\ne l=1}^{{L}_{x}}{\sigma }_{k\alpha }^{+}{\sigma }_{l\alpha }^{-}+{\rm{h.c.}},\end{eqnarray*}$$

where we have now explicitly written the summation limits.

0-excitation sector. Here the situation is trivial, the unique ground state being simply $| {\psi }^{0}\rangle =| 00..0\rangle$, i.e. all spins down and correspondingly ${s}_{l\mu }^{z}\equiv \langle {\psi }^{0}| {\sigma }_{l\mu }^{z}| {\psi }^{0}\rangle =-1,\,\forall l,\mu$

1-excitation sector. In the basis $\{| 100..0\rangle ,| 010..0\rangle ,\ldots ,| 000..1\rangle \}$ of single particle states, the interaction Hamiltonian (6c) takes a simple form

$$\begin{eqnarray}&&{H}_{1\mathrm{exc}}=2{\lambda }_{a}{M}_{a}\otimes {{\mathbb{1}}}_{b}+2{\lambda }_{b}{{\mathbb{1}}}_{a}\otimes {M}_{b},\end{eqnarray} \tag{ 39 }$$

where the a (b) matrices have dimensions ${L}_{x}\times {L}_{x}({L}_{y}\times {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-{\delta }_{{ij}}$. The corresponding eigenvalues and multiplicities are

$$\begin{eqnarray}\boxed{\begin{array}{c}\begin{array}{ccc}{\rm{M}}\mathrm{atrix} & {\rm{E}}\mathrm{igenvalue} & {\rm{M}}\mathrm{ultiplicity}\\ {M}_{a} & -1 & {L}_{x}-1\\ & {L}_{x}-1 & 1\\ {M}_{b} & -1 & {L}_{y}-1\\ & {L}_{y}-1 & 1\end{array}\end{array}}\end{eqnarray} \tag{ 40 }$$

Since ${\lambda }_{a}\gt 0$ and ${\lambda }_{b}\lt 0$, the minimum energy is

$$\begin{eqnarray}&&{E}_{\min }^{1\mathrm{exc}}=-2{\lambda }_{a}+({L}_{y}-1)2{\lambda }_{b}\end{eqnarray} \tag{ 41 }$$

with multiplicity ${L}_{x}-1$. The corresponding eigenvectors are

$$\begin{eqnarray}&&{E}_{a}=-1\to | {v}_{a}^{j}\rangle =\displaystyle \frac{1}{\sqrt{2}}(| {1}_{1}\rangle -| {1}_{j}\rangle ),\,\,\,j=2..{L}_{x},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{E}_{b}={L}_{y}-1\to | {v}_{b}\rangle =\displaystyle \frac{1}{\sqrt{{L}_{y}}}\displaystyle \sum _{j=1}^{{L}_{y}}| {1}_{j}\rangle ,\end{eqnarray} \tag{ 43 }$$

where $| {1}_{j}\rangle \equiv | 0..{01}_{j}0..0\rangle$. The ground state eigenvectors of ${H}_{1\mathrm{exc}}$ are then $| {\psi }_{j}^{0}\rangle =| {v}_{a}^{j}\rangle \otimes | {v}_{b}\rangle$ and the spin expectation values become

$$\begin{eqnarray}&&\langle {\psi }_{j}^{0}| {\sigma }_{l\mu }^{z}| {\psi }_{j}^{0}\rangle =\displaystyle \frac{1}{2{L}_{y}}-\left(1-\displaystyle \frac{1}{2{L}_{y}}\right)=-1+\displaystyle \frac{1}{{L}_{y}}\end{eqnarray} \tag{ 44 }$$

if $| {\psi }_{j}^{0}\rangle$ contains an excitation at site $l\mu$ or −1 otherwise. We note that ${s}_{l\mu }^{z}\to -1$ in the thermodynamic limit as one would expect.

In summary, we have for the ground state energies

$$\begin{eqnarray}&&{E}_{\mathrm{GS}}^{0\mathrm{exc}}=\displaystyle \frac{{\omega }_{\mathrm{at}}^{\prime }}{2}(-N),\end{eqnarray} \tag{ 45a }$$

$$\begin{eqnarray}&&{E}_{\mathrm{GS}}^{1\mathrm{exc}}=\displaystyle \frac{{\omega }_{\mathrm{at}}^{\prime }}{2}(-N+2)-2{\lambda }_{a}+({L}_{y}-1)2{\lambda }_{b},\end{eqnarray} \tag{ 45b }$$

where

$$\begin{eqnarray}&&{\omega }_{\mathrm{at}}^{\prime }={\omega }_{\mathrm{at}}+2({\lambda }_{a}+{\lambda }_{b}).\end{eqnarray} \tag{ 46 }$$

We are now in position to diagonalize the photonic quadratic form (38). 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}\equiv {s}_{l\mu }^{z}=-1,\,\forall l,\mu$, so that the photonic Hamiltonian becomes

$$\begin{eqnarray}&&{H}_{\mathrm{ph}}={p}^{\dagger }{M}_{\mathrm{ph}}p,\end{eqnarray} \tag{ 47 }$$

where we have introduced $p={({a}_{1},..,{a}_{{L}_{y}},{b}_{1},..,{b}_{{L}_{x}})}^{T}$. The matrix ${M}_{\mathrm{ph}}$ can be written as

$$\begin{eqnarray}{M}_{\mathrm{ph}}=\left(\begin{array}{cc}{W}_{a} & G\\ {G}^{T} & {W}_{b}\end{array}\right)\end{eqnarray} \tag{ 48 }$$

with

$$\begin{eqnarray}&&{({W}_{a})}_{{ij}}={\delta }_{{ij}}\left[{{\rm{\Delta }}}_{a}+2{\lambda }_{a}\displaystyle \sum _{\nu =1}^{{L}_{x}}{s}_{i\nu }^{z}\right],\,\,\,i,j=1..{L}_{y},\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{({W}_{b})}_{{ij}}={\delta }_{{ij}}\left[{{\rm{\Delta }}}_{b}+2{\lambda }_{b}\displaystyle \sum _{l=1}^{{L}_{y}}{s}_{{lj}}^{z}\right],\,\,\,i,j=1..{L}_{x},\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{G}_{{ij}}=({\lambda }_{a}+{\lambda }_{b}){s}_{{ij}}^{z},\,\,\,i=1..{L}_{y},j=1..{L}_{x}.\end{eqnarray} \tag{ 51 }$$

${M}_{\mathrm{ph}}$ has the following eigenvalues and multiplicities

$$\begin{eqnarray}\boxed{\begin{array}{c}\begin{array}{cc}{\rm{E}}\mathrm{igenvalue} & {\rm{M}}\mathrm{ultiplicity}\\ {E}_{a}^{\mathrm{ph}}\equiv {{\rm{\Delta }}}_{a}+2{\lambda }_{a}{L}_{x}{s}^{z} & {L}_{y}-1\\ {E}_{b}^{\mathrm{ph}}\equiv {{\rm{\Delta }}}_{b}+2{\lambda }_{b}{L}_{y}{s}^{z} & {L}_{x}-1\\ {E}_{+}^{\mathrm{ph}}\equiv \tfrac{1}{2}(\epsilon +\xi ) & 1\\ {E}_{-}^{\mathrm{ph}}\equiv \tfrac{1}{2}(\epsilon -\xi ) & 1\end{array}\end{array}}\end{eqnarray} \tag{ 52 }$$

where

$$\begin{eqnarray}&&\epsilon ={{\rm{\Delta }}}_{a}+{{\rm{\Delta }}}_{b}+2{s}^{z}({\lambda }_{a}{L}_{x}+{\lambda }_{b}{L}_{y}),\end{eqnarray} \tag{ 53a }$$

$$\begin{eqnarray}&&\xi =\sqrt{{({{\rm{\Delta }}}_{b}-{{\rm{\Delta }}}_{a}+2{s}^{z}({\lambda }_{a}{L}_{x}-{\lambda }_{b}{L}_{y}))}^{2}+4{L}_{x}{L}_{y}{({s}^{z})}^{2}{({\lambda }_{a}+{\lambda }_{b})}^{2}}.\end{eqnarray} \tag{ 53b }$$

First we note, that due to (36), ${{\rm{\Delta }}}_{b}$ is not independent and can be expressed as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{b}=\displaystyle \frac{{{\rm{\Delta }}}_{a}-{\omega }_{\mathrm{at}}}{\eta }+{\omega }_{\mathrm{at}}\end{eqnarray} \tag{ 54 }$$

(here ${\omega }_{\mathrm{at}}$ is the bare atomic frequency, not ${\omega }_{\mathrm{at}}^{\prime }$). 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

$$\begin{eqnarray}&&Q=\min \left(\displaystyle \frac{| {{\rm{\Delta }}}_{a}-{\omega }_{\mathrm{at}}| }{{g}_{c}^{\mathrm{spin}}},\displaystyle \frac{| {{\rm{\Delta }}}_{b}-{\omega }_{\mathrm{at}}| }{{g}_{c}^{\mathrm{spin}}}\right),\end{eqnarray} \tag{ 55 }$$

where ${g}_{c}^{\mathrm{spin}}$ is given by the critical value of ${\lambda }_{a}$, ${g}_{c}^{\mathrm{spin}}=\sqrt{-2{\lambda }_{a,c}({{\rm{\Delta }}}_{a}-{\omega }_{\mathrm{at}})}$, see below.

We start by determining the critical value of ${\lambda }_{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}_{\mathrm{GS}}^{0\mathrm{exc}}={E}_{\mathrm{GS}}^{1\mathrm{exc}}$ and with the help of (45) we get

$$\begin{eqnarray}&&{\lambda }_{a,c}^{\mathrm{spin}}=-\displaystyle \frac{{\omega }_{\mathrm{at}}}{2\eta {L}_{y}}.\end{eqnarray} \tag{ 56 }$$

Next, the breakdown of the effective model is indicated when any of the photonic eigenvalues (52) becomes negative, corresponding to infinitely many photons in the ground state. Substituting ${s}^{z}=-1$ in (52), the only two candidates for the minimum eigenvalue are ${E}_{a}^{\mathrm{ph}}$ and ${E}_{-}^{\mathrm{ph}}$ (we recall that ${\lambda }_{a}\gt 0,{\lambda }_{b}\lt 0$). The value of ${\lambda }_{a}$ where ${E}^{\mathrm{ph}}$ becomes negative is determined as

$$\begin{eqnarray}&&{\lambda }_{a,c}^{\mathrm{ph}}\equiv \min ({\lambda }_{a}\,:\,{E}_{a}^{\mathrm{ph}}({\lambda }_{a})=0,{E}_{-}^{\mathrm{ph}}({\lambda }_{a})=0),\end{eqnarray} \tag{ 57 }$$

where only the positive branch of ${\lambda }_{a}$ in the solutions of ${E}_{-}^{\mathrm{ph}}({\lambda }_{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

$$\begin{eqnarray}&&R\equiv \displaystyle \frac{{\lambda }_{a,c}^{\mathrm{ph}}}{{\lambda }_{a,c}^{\mathrm{spin}}}\gt 1\,\,\wedge \,\,Q\gg 1.\end{eqnarray} \tag{ 58 }$$

In figure 6 we plot the contours of constant R (a) and Q (b) respectively in the $\eta -{L}_{x}/{L}_{y}$ plane. It is apparent from the figure that increasing both ${L}_{y}/{L}_{x}$ and $| \eta |$ 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

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{L}_{y}\to \infty }R=\infty ,\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\eta \to -\infty }R={R}_{\mathrm{asymptotic}}\end{eqnarray} \tag{ 60 }$$

when keeping all the other parameters fixed, as expected from the contours in figures 6(a) and (b) (here ${R}_{\mathrm{asymptotic}}$ is some finite asymptotic value).

Zoom In Zoom Out Reset image size

So far we have concentrated 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 it is possible to extend the above analysis to the ${N}_{\mathrm{exc}}\to {N}_{\mathrm{exc}}+1$ spin transitions for ${N}_{\mathrm{exc}}\geqslant 1$ and for each of the transitions evaluate the critical ratio R. We note that the fully analytical approach is unfortunately obscured by the fact that for excitation sectors ${N}_{\mathrm{exc}}\gt 1$, the spin interaction matrix (6c) does not take the simple structure of (39) and the evaluation of eigenvalues in principle amounts to solving higher order polynomials. Specifically, for a ${N}_{\mathrm{exc}}\to {N}_{\mathrm{exc}}+1$ spin transition, the critical value of the coupling ${\lambda }_{a}$ can be obtained from the condition ${E}_{\mathrm{GS}}^{{N}_{\mathrm{exc}}}={E}_{\mathrm{GS}}^{{N}_{\mathrm{exc}}+1}$, where ${E}_{\mathrm{GS}}^{{N}_{\mathrm{exc}}}=\tfrac{{\omega }_{\mathrm{at}}^{\prime }}{2}(-N+2{N}_{\mathrm{exc}})+\delta {E}^{{N}_{\mathrm{exc}}}$. Here, $\delta {E}^{{N}_{\mathrm{exc}}}$ is obtained from exact diagonalization of (6c). Furthermore, in order to evaluate the critical value of ${\lambda }_{a}$ for the photonic transition (57) (needed for the evaluation of R, (58)), we assume ${s}_{l\mu }^{z}={s}^{z},\,\forall l,\mu$ in (48), where we take ${s}^{z}({N}_{\mathrm{exc}})=-1+\tfrac{2{N}_{\mathrm{exc}}}{N}$.

In practice, the restriction of the Hamiltonian to a given ${N}_{\mathrm{exc}}$ sector clearly simplifies the analysis, however the number of basis states still grows rapidly as $\left(\genfrac{}{}{0em}{}{N}{{N}_{\mathrm{exc}}}\right)$ and an extensive numerical investigation of higher excitation number sectors goes beyond the scope of the present article. For that reason, and in order to unambiguously show that the frustrated regime can be reached in the relevant ${N}_{\mathrm{exc}}\gt 1$ sector, we focus on the simpliest such situation, namely the $1\to 2$ spin transition. In figure 6 we plot R (c) and Q (d) for two distinct values of η ($\eta =-1$ and $\eta =-0.5$) for $0\leftrightarrow 1$ and $1\leftrightarrow 2$ spin transition. The solid lines correspond to the analytical result (56) and coincide with the exact diagonalization results for the $0\leftrightarrow 1$ (black crosses and red stars for $\eta =-1$ and $\eta =-0.5$ respectively)¹⁰ . It is apparent from figures 6(c) and (d) that the desired parameter range $R\gt 1$, $Q\gg 1$ can be achieved also for the $1\leftrightarrow 2$ spin transition and therefore the frustrated regime can be indeed reached (this is also a confirmation of the intuitive expectation based on the analysis of the $0\to 1$ transition alone, namely that the arbitrarily large values of R are a strong indication that higher excitation number sectors can be reached without breaking the validity of the spin Hamiltonian). The smallness of the change of R and Q between the $0\leftrightarrow 1$ and $1\leftrightarrow 2$ transitions is a consequence of the fact that the difference between the critical values of ${\lambda }_{a}$ for two adjacent excitation numbers spin transitions (${N}_{\mathrm{exc}}\to {N}_{\mathrm{exc}}+1$ and ${N}_{\mathrm{exc}}+1\to {N}_{\mathrm{exc}}+2$) decreases rapidly with the system size (another manifestation of this is the finite size scaling of the superradiant phase transition, see the insets of figure 2, where the critical coupling tends to 0 in the thermodynamic limit.)

To recap, we have shown that the spin model (5) with both frustrated and non-frustrated interactions can emerge as an effective description of the parent TC Hamiltonian (2). This can be achieved when considering an elongated geometry of the square array, ${L}_{y}\gg {L}_{x}$. On one hand this circumvents the limitations related to realizing effective spin Hamiltonians with frustrated interactions using optical setups governed by TC Hamiltonians [63]. Finally, we note that one can simulate the parent TC 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.