Fragmentation and destruction of the superfluid due to frustration of cold atoms in optical lattices

A one-dimensional optical lattice is considered where a second dimension is encoded in the internal states of the atoms giving effective ladder systems. Frustration is introduced by an additional optical lattice that induces tunnelling of superposed atomic states. The effects of frustration range from the stabilization of the Mott insulator phase with ferromagnetic order, to the breakdown of superfluidity and the formation of a macroscopically fragmented phase.


Introduction
One of the most interesting open fields in quantum magnetism is the study of frustration [1]. Frustrated models are described by Hamiltonians with competing local interactions such that the ground state cannot minimize their energy simultaneously. Frustration can appear due to the geometry of a problem, the Ising model on a triangular lattice being a paradigmatic case, or due to the coexistence of ferro-and antiferromagnetic interactions, as in spin glasses. Frustrated models typically have highly degenerate ground states, which can become ordered by increasing the temperature or by quantum fluctuations-i.e. 'order by disorder'. Theoretical and numerical problems, such as the large dimensionality or the sign problem in Monte Carlo simulations make it very difficult to study frustrated Hamiltonians.
Quantum simulation has been put forward as a tool to probe the physics of a variety of manybody systems. In particular, schemes with cold atoms in optical lattices have been suggested to simulate arbitrary spin models [2]- [7] and some interesting frustrated Hamiltonians [8,9]. Compared to the alternative of looking for or inducing low-dimensional behaviour on existing magnetic materials [10,11] and in organic conductors [12,13], cold atoms and molecules offer greater flexibility in terms of variable geometry and interaction strength. However, the majority of current proposals are perturbative and their effective interactions are rather weak. This makes their experimental realization challenging, as it requires very low temperatures. A solution suggested to solve the problem of weak interactions is to replace the atoms with polar molecules [14].
In this paper, we follow a different approach. First of all, spin up and spin down states are identified with particle and hole configurations. This gives naturally an XX coupling which is proportional to the tunnelling amplitude [15] and, thus, it is strong and less experimentally demanding. By using the internal state of the atoms to encode a second virtual dimension and applying a spatially-dependent Raman coupling, we are able to induce in-plane frustration. The result is a new frustrated model, which contrasts with related literature on spin ladders, where the Hamiltonian either uses antiferromagnetic isotropic [16]- [21] or XXZ [22] interactions, or uses purely XY interactions but requires a different lattice geometry [23], or both. More important, the problem we consider here is a full Bose-Hubbard Hamiltonian, and the main result is that the transport properties of the lattice change dramatically, causing the breakdown of superfluidity even for large densities in which the analogy to spin models no longer applies. In particular, in the most interesting case the Mott insulator is replaced by a new dimerized phase, where neighbouring pairs of sites have maximal coherence and there is a fast decay of coherence between pairs of sites as a function of distance. This new phase is reminiscent of a Bose glass [24], with the important difference that it has been induced by frustration and not by disorder. Our statements are supported by a variety of analytical, variational and numerical solutions.
The paper is organized as follows. In section 2, we introduce the two Hubbard models that we are going to use, explain how they are implemented with cold atoms in optical lattices and briefly justify the use of the word 'frustration'. In section 3, we focus on a model with diagonal interactions across each square plaquette. We first discuss the different phases as obtained from density matrix renormalization group (DMRG) calculations, including qualitative arguments about why these phases are expected. Afterwards we explain in more detail how fragmentation happens and what happens in the limit of strong interactions, and what are the correlation properties of all available phases. In section 4, we briefly present a model with frustrated rectangular lattice. In section 5, we study in great detail all issues regarding the implementation of our ideas in current experiments with optical lattices. First we do a microscopic derivation of the parameters in the Hubbard Hamiltonians, including all possible sources of error. The main conclusion is that the requirements for studying these frustrated models appear to be within reach of current experiments and that there are no side effects from introducing a Raman coupling in the experiment. We explain how the different phases studied in this paper could be detected, either from time of flight measurements or more sophisticated correlation measurements. Finally, we discuss possible sources of imperfection, such as temperature or a residual harmonic confinement. The last section (section 6) contains a brief summary of our main results and possible implications and connections to other work.

Bose-Hubbard model
The system that we are studying is that of cold atoms confined by an off-resonance optical lattice that forms a one-dimensional (1D) bosonic lattice gas [15]. As in recent experiments [25], the lattice traps atoms in two internal states. If the confinement is strong, the effective model will be a Hubbard Hamiltonian [26] where J is the hopping amplitude and we assume that the on-site atomic interactions are both repulsive and symmetric, U σσ = U > 0. Note that in this work, the internal state of the atoms are denoted by Greek letters, σ =↑, ↓; the site indices are denoted by italic characters, i, j, and go  (5) replaces the frustrating terms with weak and strong tunnelling. (d)-(f) Similar drawings for the Villain model in equation (2). from 1 to L which is the length of the lattice; finally N := N ↑ + N ↓ , represents the total number of particles.
The new ingredient in our Hamiltonian are two counter-propagating laser beams that induce Raman transitions between the two atomic states. By adjusting the phase, the polarization and the alignment of these beams, we can ensure that the effective Rabi frequency, (x), forms a lattice with twice the period of the confinement (figure 1). We can interpret our 1D system as having two 1D lattices, one for each internal state of the atom. Both chains are coupled forming a ladder thanks to both the interaction, U, and to H R , a Raman term with the spatially-dependent Rabi frequency, (x).
We shall consider two specific configurations of the Raman coupling. When its maxima and minima coincide with those of the confining lattice (figure 1(d)), the result will be on-site Rabi oscillations Alternatively, when they are displaced by half a period we will have diagonal interactions ( figure 1(a))

Origin of frustration
The frustration induced by H R becomes evident in the limit of hard-core bosons, a 2 iσ = 0 [15]. Identifying the bosonic operators, a i↑ and a i↓ , with the spin operators σ − i1 and σ − i2 , the Hamiltonian (1) becomes an XX model on a spin ladder (figures 1(b) and (e)). Thus, while the hopping J translates into a ferromagnetic XX interaction along the legs, the Raman coupling H R is a transverse XX interaction which alternates from ferro-to antiferromagnetic. Take for instance the diagonal interactions (3). The associated spin model is The frustration arises from the competition of positive and negative contribution terms. For instance, looking at the sites (1, j), (1, j + 1), (2, j + 1) and (1, j + 2): the bonds on this triangle have three ferro-and one antiferromagnetic terms, and there is no way in which the spins can be aligned so as to minimize the energy of all terms (see also figure 1(b)). While the previous analogy is somewhat pleasing, one may wonder whether having more than one particle per site or working with a system which is not properly 2D will give rise to trivial physics. Regarding the first point, when we have more than one atom per site our model becomes equivalent to an array of Josephson junctions [9] where the frustration is still present and lays on the choice of phase for each individual site. Regarding the second point, some of the best studied models in frustrated quantum magnetism, such as the Majumdar-Ghosh [27,28] and the zig-zag chain [22] are also quasi-1D models. As we will show in the following sections, the competition between the different terms in equation (1) does indeed produce new and interesting effects with the advantage of being analytically tractable.

Diagonal interactions
In contrast to fully 2D systems, our quasi-1D models can be greatly simplified by rewriting the frustrated coupling using dressed states We will begin by studying the Raman interaction (3). With the previous change of variables, our Hamiltonian becomes where n i := c † i+ c i+ + c † i− c i− is the total population of a single site. Equation (6) models a ladder whose legs are made of junctions with strong, J + = J + J , and weak hoppings, J − = |J − J |, as shown in figure 1c. Notice that the model is symmetric under the exchange of J and J . In the following sections, we will often exchange descriptions between J + , J − and J, J , and we will assume, without loss of generality, that all these parameters are not negative.

Phase diagram
In this section we present zero temperature results for Hamiltonian (6). We have computed the ground states and first excitations of this model using the DMRG method in the matrix product states (MPS) formalism [29,30]. We have conducted accurate simulations for L = 8, 16, 32, 50 and 70 sites, using a cutoff of four particles per spin. These values give us a local dimension of the Hilbert space of 25 states, which is very big and can only be handled with state-of-the art optimizations, such as working in sectors with well-defined number of particles or angular momenta. Due to the large size of the local Hilbert space, we have typically worked with MPS of size 100 and checked convergence for different points with up to D = 200. 6 In our simulations we have computed the ground state and lowest energy states for N/L = 1/2, 1, 5/4, 3/2, 7/4, 2 particles, as well as the energy to add a particle, µ p := E N+1 − E N , or to create a hole, µ h := E N − E N−1 , and verified that these values did not change with a larger cutoff.
These zero temperature simulations reveal a very simple picture that we summarize here. As illustrated in figure 1(d), when we increase U for fixed J + and J − , the system first experiences a phase transition from a superfluid to a fragmented or dimerized phase in which coherence is maximal between pairs of neighbouring sites. As the interaction is further increased these fragments evolve smoothly (i.e. a crossover) to a Mott insulator. While the first phase is gapless and its excitations are phonons, the crossover contains a gapped phase whose lowest excitations are localized spin flips.
In the following sections, we will study in detail the properties of the different coupling regimes, numerically as well as analytically where possible. Nevertheless it is possible to obtain a qualitative picture of the three regimes and the associated ground states by variational wavefunctions. First of all, for weak interactions, U J ± , the ground state is a uniform superfluid that spread over the lattice with the usual rotational symmetry on the (α, β) space. For stronger interactions, J − < U/4 J + , the energy of |ψ SF is larger than a state made of fragmented condensates which reside on the junctions with large hopping, J + (figure 1(c)), Here, for integer filling, N/L = 1, 2, 3, . . ., the ground state is one with spontaneously broken rotational symmetry, a ferromagnetic state with either n + = 0 or n − = 0. This effect arises from the contributions to the energy of both the interaction U and the hopping J − , as it will be explained in section 3.2.
Another symmetry break happens when the interactions become dominant, U J ± . Indeed, in the limit of strongly interacting bosons and N = L, all ground states can be written in the form The quantum fluctuations generated by the hopping select, among all disordered ground states, one particular ferromagnetic state with either n − = 0 or n + = 0. This case is further discussed in section 3.3. In the space of parameters, hopping versus chemical potential, J/U vs µ, we can draw a phase diagram similar to that of the Mott-insulator transition [24]. As figure 2 shows, for J = 0.9J, the curves of particle and hole excitations, µ p and µ h , create lobes (grey areas) containing fragmented phases. Crossing these lines amounts to a quantum phase transition to a superfluid, on the exterior of the lobes. Qualitatively, these lobes are similar to those of the Bose-Hubbard model [24], but about J/|J − J | times larger, so that for J = J the lobes become infinitely large and all ground states are fragmented. In this limit of J = J and for integer and half-integer filling, that is N/L = 1, 3/2, 2, . . ., we have on phase space infinitely many incompressible 'stripes' whose borders, µ p and µ h , can be computed exactly (see section 3.2 and equation (12)). As explained above, within each stripe the state of the atoms changes smoothly from a Mott insulator (small J + ) to fragmented or dimerized phase where atoms delocalized between pairs of neighbouring sites (large J + ). Interestingly, in the macrocanonical ensemble, all phases with other filling factors collapse to the lines between these regions. The large degeneracies of these other fillings has made it impossible for us to study systematically the nature of those states and of their excitations.
To gain further insight on the properties of the ground state, let us focus on the unit filling case. The transition from the Mott phase to the superfluid is accompanied by the change of several excitation gaps. One is related to the energy required to add or remove a particle, the other one is the energy gap required to create a density perturbation for fixed number of particles. Both gaps close at the same point when we jump from the fragmented phase to the uniform superfluid, signalling a phase transition from an incompressible to a compressible phase and also the establishment of algebraically decaying correlations (figure 3). This transition occurs around the critical value of the interaction for a single-component 1D Bose-Hubbard model [31]- [34] with hopping J − , that is U ∼ U 1D := 3.84J − (see section 3.3). A third energy gap is present for excitations that take a particle out of one of the coherent fragments into a different spin state, while preserving the total number of particles. This gap freezes the dynamics of the spin of the atoms. As shown in figures 3(a) and (b), it is maximal in the fragmented phase, with a maximum value of order U/4. For strong interactions it grows as J 2 /U while in the superfluid phase one would expect it to close as U → 0, contrary to what the numerical simulations reveal.
Intuitively, this small gap could be due to the fact that an atom changing the spin means it jumps to the empty superlattice and increases the total repulsive interaction. However, the fact that for U J ± the ground state becomes gapless implies that our DMRG estimates for the spin gap are unreliable and should be confirmed using independent analytical and numerical tools.

Fragmented state
In this subsection, we will comment on the configurations that appear for U ∈ [J − , J + ). For simplicity we will begin with the case of balanced hopping and Raman coupling and no interaction, U = J − = 0. Then it is trivial to see that all ground states will be of the form (8) with atoms delocalized on pairs of sites. In particular, for unit filling N = L we have three degenerate states: a uniform state with all junctions populated n ± = 1 (figure 4(b)) and two states with broken translational symmetry, n − = 2 and n + = 2 ( figure 4(a)). This degeneracy persists for small U, but it is broken as soon as we set a weak imbalance in the hoppings. For any 0 < J − U, terms proportional to J − connects these states off-resonantly to excited states like It follows that to second order in perturbation theory the degeneracy of these states is broken and the energy of the uniform state, n + = n − , is shifted by an amount E = −8(J 2 − /U) × L smaller than the ferromagnetic states with broken symmetry, E = −24(J 2 − /U) × L, which become the true ground state. Once we have established that the system spontaneously selects either the c + or c − order, we can compute the ground state for any value of U while keeping J − = 0. Such state will be a product of ground states of the two-well problem, Heren = N/L = 1 is the density, φ (n) m is the wavefunction of the double-well problem withn atoms and σ = ± is the selected polarization of the atom.
To delimit the incompressible regions in phase space we only need to compute the energies of a state with N = L, L + 1 and L − 1 particles, which are variations of the one in equation (10). Since the ground state energies of a double well with 1, 2 and 3 particles are, respectively, the boundaries of the Mott region with n = 1 particle per site are given by These boundaries are the dash-dot lines we plotted in figure 2. A simple inspection of the previous formulae shows that these curves never cross and become parallel in the limit of J + /U 1. The previous exact results, which are valid for J − = 0, can be extended to J − U, J + by means of perturbation theory. In that case there will be an additional coupling between fragments up to a distance of the order of the correlation length. The wavefunction can be approximated by linear combinations of the terms (10) with different occupations of the double wells. The corrections are of order O(J − /U) to the wavefunction, and of order O(J 2 − /U) to the energy, but the dominant term is always the one with equally populated double wells. This will be further explained in section 3.3, where we analyse the location of the Mott insulator to superfluid transition. That this picture is applicable to a great degree is appreciated in figure 5, where we plot the correlations between different sites. For J = 0.9, sites connected by J + are strongly correlated, while the correlations decrease significantly between different pairs of sites.

Hard-core bosons
To understand the behaviour of the bosons in the strongly interacting limit, U J ± it is useful first to see what happens if there are no frustrating terms, J = 0. We know that on the one hand, if we are below half-filling, N < L, the gas is a Tomonaga-Luttinger liquid with spin-charge separation [35]. In the limit of U → ∞ the spin degrees of freedom become degenerate and the whole system behaves like free fermions or a Tonks gas. On the other hand, for half filling, N = L, the sample becomes a Mott for any strong interaction U J. The atoms cannot move but their spin degrees of freedom interact with a ferromagnetic Heisenberg interaction of strength t = −J 2 /U and a gapless spectrum made of spin waves. The diagonal coupling J changes this landscape slightly. For half filling, N = L, we still have a Mott phase, with one atom per site. Connecting with the previous section, as we increase U the wavefunction of the two-well fragment (10) becomes closer and closer to φ (n) m = δ m,n/2 . This is a smooth evolution which is why we speak of a crossover. Furthermore, compared to the unfrustrated case J = 0, the spin rotational symmetry is now broken, favouring the X-direction or c ± states in the second quantization language. The effective spin model is, up to constants, a ferromagnetic one and spin excitations are prevented by an energy gap = 2J 2 /U [36]. This gap is evident also in figure 3(a), in the region of small hopping, J U. For smaller fillings, i.e. N < L, the previous spin model is no longer valid because the chain has holes. Nevertheless one can easily compute things if J − = 0. The ground state manifold is composed of variations of the fragmented wavefunction (10). In particular, the set of available configurations is shown in figure 6, where we show three kinds of bonds. Those with zero or two particle, which are frozen, and bonds with one particle on the left or right sites, which can be identified with spins. These effective spins form ferromagnetic blocks which are separated by the frozen sites. If the block is surrounded by at most one doubly occupied bond, then it can only adopt one state and it has zero energy ( figure 6(d)). If the block, on the other hand is surrounded by empty blocks it may host at most one domain wall (figure 6(c)) with some small negative energy due to the motion of the domain wall (for instance, the states in figures 6(b) and (c) are connected by hopping). The effective Hamiltonian for a block of size M is where

Scaling of the Mott transition
In this subsection, we attempt to analyse how the superfluid to insulator transition changes due to frustration. We will do it in two different perturbative methods. The first one is based on a strong coupling expansion around J, J = 0, while the second one is based on a strong coupling expansion around J = J. Remarkably, this second method, which is more accurate than the trivial strong coupling expansion [33], may be contemplated as a real space renormalization technique in which the alternating Hubbard model is replaced by an effective one with a renormalized hopping and half the sites. Our first attempt starts with the estimates of the critical value of J/U given in [33]. Applying perturbation theory up to O(J 3 /U 3 ), it is possible to estimate the change in the energy induced by a small amount of hopping. The resulting energies can then be used to compute the particle and hole excitation gaps, and delimit the borders of the incompressible regions in the hopping versus chemical potential parameter space, as in figure 2. To apply the theory in [33] we need the hopping matrix 7 the smallest eigenvalue of (−t ij ) the wavefunction of the single-particle ground state, f i = 1/ √ L, and finally the sums As shown in figure 7, the transition does indeed shift to larger values of J/U as we switch on the diagonal coupling J . However, the prediction is not very consistent between second and third order expansions, and it even fails to provide the right value for J = 0, probably due to the asymptotic nature of these expansions. 7 We assume periodic boundary conditions. 84J computed with DMRG for unfrustrated models [32]. Note the agreement between the two latter methods for J = 0.
Our second attempt focuses on the regime of strong frustration or values of J close to J. We have seen that for J = J the incompressible regions fill the entire phase space. With very simple arguments it is possible also to compute how they shrink for nonzero |J − J| = U. The reasoning begins by noticing that our system is a superlattice, where the weak hopping terms can be treated perturbatively. However, instead of using the Mott states as unperturbed states one must use the exact ones for J + > 0. Basically, we have an unperturbed state with energy E = L 2 /2 and several unperturbed manifolds above it, consisting of holes and excitations in the superlattice, with energy gaps of order 3 − 2 , 1 − 2 . We can compute the corrected µ p and µ h as a function of J − . The equation for the critical point then looks as follows where the energies and the matrix elements of the hopping Hamiltonian, c(U), are functions of U. We have solved numerically the highly nonlinear equation (18) for U as a function of J − and J + . The results are plotted in figure 7. As seen there, this second method interpolates properly between the critical value of U = 3.84 obtained by DMRG for the unfrustrated model [32], and the critical value U = 0 computed for J = J . Furthermore, the whole line deviates very little from U c = U c (J = 0)|J − J |.

Correlators
As it will be useful to identify the different phases experimentally, we have computed the singlebody and two-body correlation functions in the three regimes. The most interesting values are the single-body correlators, a † i a j . In the superfluid limit, due to the 1D character of our problem, we expect that correlations decay algebraically. In the fragmented regime, we have two possible ground states, depending on whether the fragments sit on even or odd junctions, which is equivalent to specify whether the fragments are made of c + or c − particles. In the first case we will have with all other correlators being zero, and in the c − case we will have the same correlation pattern but displaced by one lattice site. Finally, in the Mott insulator regime we have no correlations and c † i± c j± frag ∼ δ i−j . In experiments what is measured is the time of flight images. These pictures are related to the momentum distribution of the sample or the Fourier transform of the single-particle correlations This function is peaked at q = 0 and decays towards the sides of the Brillouin zone. The halfwidth of the peak is related to the correlation length of the sample. In the superfluid region it will be proportional to the size of the lattice and limited by finite temperature effects [37]. As soon as we cross the phase transition towards the fragmented phase, though, the correlation length decays abruptly to ξ = 2 sites and the previous function has a constant width n q ∼ 1 + cos(2πq/L).
From there on we expect a smooth crossover towards n q = 1 which is the Mott-Insulator regime. All this phenomenology is evident in figure 9.

Villain model
The other limiting case in our set-up is that of vertical frustrating interactions, a configuration resembling the odd model or Villain lattice [1]. Once more it is convenient to change basis as in equation (5), so that the frustrating terms become alternating energy shifts The effect of the Raman coupling is now equivalent to a superlattice. This additional potential splits the Wannier band into two effective bands separated by a gap J , each band having a reduced width J := ( √ 4J 2 + J 2 − J )/2. The superlattice localizes the atoms on alternating rungs, forming an antiferromagnetic configuration. Instead of exhibiting fragmentation, the system goes straight from a superfluid to a Mott-insulator, but now the transition happens for weaker interactions, U ∼ J.

Experimental realization
There are several issues that one has to consider when implementing our Hubbard models using neutral atoms. The first one is how to relate U, J and J with experimental parameters such as the intensity of the Raman laser and the strength of the optical lattice. Additionally one has to consider how to prepare the ground state and how cold the sample should be. Finally one has to think about procedures to detect the different phases. We will address all of these questions in the following sections.

Microscopic theory of Hubbard model
In this subsection, we show how to derive the constants U, J and J using a band structure calculation that takes into account the Raman coupling, (x) = 0 cos(kx + φ), and the confining lattice, Our work generalizes the ideas from [26], but now our microscopic Hamiltonian contains a potential term that includes not only the trapping potential, but the coupling and the interaction between atomic states. Following [26], we perform a tight-binding approximation and replace the bosonic operator by an expansion of the form where b = π/k is the period of the lattice and w x,y,z are the Wannier function corresponding to the lowest Bloch band of the lattice, along the X, Y and Z-directions, respectively. By substituting this expansion into the full Hamiltonian (24) and retaining the most important terms, one arrives to a Bose-Hubbard model (1).
To compute the parameters of our model, we will assume that the confinement along the transverse directions, V 0⊥ , is big enough, so that hopping is negligible and our system is made of multiple, disconnected 1D lattices. However, even with this assumption there are still too many 'knobs' to tune. To simplify the problem it is customary to adimensionalize all magnitudes using the energy of the lattice, E R =h 2 k 2 /2m, and its period, b, as units. We replace the trapping potential with V ∼ V 0 cos(2kx), make the change of variables 2kx → 2πx and rescale the wavefunction accordingly, w → w √ k/π. As a result we obtain formulae for the hopping the Raman coupling and the interaction energy Notice that w x (x) is now a Wannier function computed for a single, 1D problem with period 2π and trapping V 0 /E R . The only other free parameters are 0 /E R , and the effective 1D scattering length 8 At this point one has to solve a family of 1D problems with the periodic potential cos(2πx) in order to obtain the relative values of J, J and U.

Experimental parameters
From the computations sketched in section 5.1, it follows that if we choose the Rabi coupling of similar strength to the confinement, 0 = V 0 /2, then the ordinary hopping, J, and the ladder coupling, J , are also of the same order of magnitude ( figure 8). That means that by slightly tuning 0 it should be experimentally feasible to cover all possible regimes, from J − = J + to J − = 0. Now it only remains the question of whether we can make U/J small or big enough to explore the different phases. For a typical experiment with 87 Rb, we have a 3D ∼ 5 nm and the wavelength of the laser is λ ∼850 nm, so that a (1D) σσ ∼ 1.8 × 10 −3 I 2 . Now, for a strong confinement we have that I ranges from 2 to 8, and comparing with the evolution of J and J in figure 8 we see that the value of the interaction can indeed change from a superfluid regime, U < |J − J |, to the Mott insulator regime, |J + J | U. In other words, if we keep J − J not too small, the laser intensities will not differ much from those used in current experiments [15,39] and we will be able to generate all phases.
While the validity of the tight-binding approximation is well established for the ordinary Bose-Hubbard model, we have introduced a coupling between neighbouring sites that might excite the atoms to higher Bloch bands. We have studied this effect and computed the following effective strength of the coupling which is a function of the energy difference between bands, E, and the coupling between Wannier functions of the lowest and first excited bands When small, the value ε measures how much mixing of the excited wavefunctions there is in the ground state and it is a function of 0 . As figure 8 shows, for the diagonal interactions (φ = π/2) the coupling is indeed extremely small, ε < 10 −2 , and we can be sure that our approximations are valid. Similar results are obtained for φ = 0.

Preparation and observation
The remaining issues are related to the possibility of preparing and unambiguously identifying the ground states of our model (1). The state preparation can be done in a two step process. One begins with a 1D superfluid and all atoms in the same internal state, and slowly increasing the lattice depth. Judging from current experiments, this should be a robust way to create a Mott or Tonks gas phase [15]. For U J, J the density excitations will be frozen and the dynamics will be ruled by the effective ferromagnetic model (13). One can then slowly increase J to create the ground state of (13) which will be the ground state of the frustrated model.
Regarding detection, the three phases of model (3) can be identified from time of flight images. First of all, for U J ± we expect a quasi-condensate which will produce strong interference peaks in either of the atomic species. The width of these peaks may be slightly modified by temperature [37] but it should remain approximately constant for U < J ± . For moderately strong interactions, J − < U < J + , the correlation length will decrease substantially as we expect coherences only inside each fragment. Thus according to section 3.5 the interference pattern will be approximately |ψ(tp/m)| 2 ∼ |w(p)| 2 [1 + cos(p)] =: |w(p)| 2 n(p). Herew(p) is the Fourier transform of a Wannier wavefunction, n is the ideal interference pattern and the position x = tp/m is related to the time, t, the mass of the atoms, m, and its momentum in the lattice, p. In practice, for J = J, the remaining coherence between fragments can modify this pattern, giving thinner peaks, as shown in figure 9. Another signature of the transition from quasi-condensate to fragmented phase is a spontaneous symmetry break that makes the ground state be formed of either c + or c − particles. The energy gap in this phase is related to the stiffness of the atoms with respect to a change in the internal state and it could be probed spectroscopically.
A third signature which applies to both the fragmented phases and the ground state of the Villain model is the existence of magnetic order. More precisely, for strong interactions J + < U, the ground states of the diagonal and Villain frustrated models have ferro-and antiferromagnetic order, respectively. This pattern reveals itself in high order correlations between the number of atoms of c + and c − at different sites: c † jσ c jσ c † iσ c iσ ∝ 1 + σσ (−1) i−j , where = 1 or 0 for each model, respectively. To observe these correlations one must first rotate the atoms with a π/2 pulse, as in equation (5), and afterwards analyse the quantum fluctuations in the time of flight images [40]. Alternatively, if the experiment allows for state-dependent lattices and photoassociation, it will be advantageous to use the techniques suggested in [7] to measure the number of correlations. This second method should indeed provide a much stronger signal.
In current experiments, we have to consider two additional sources of imperfection that can influence our observations. One of them is the residual harmonic confinement which is always present due to the difference on intensity along the lasers that trap the atoms. From the diagram of figure 2, the harmonic confinement will cause the atomic cloud to be made of nested insulating shells with different densities [26], while the size of the superfluid regions will substantially decrease and strictly vanish for J = J . This structure can cause noise in the interference patterns, but, as it has been shown in recent experiments [41], it can also be probed using RF knifes and spin-changing collisions.
The second source of imperfections is temperature, which will influence the measurements in two ways. First of all, for T = 0, the weakly interacting region will not be a condensate, but a quasi-condensate, where phase fluctuations are energetically cheap but the system remains superfluid. One might argue that these phase fluctuations will destroy the interference peaks, but as we have seen in 1D experiments with optical lattices [39], this does not seem to happen in practice. Moreover, the correlation length of the quasi-condensate will be much larger than that of the fragmented phase, a fact that will be evident in the interference pattern. It remains the question of whether the fragmented phase itself is robust against temperature and we argue that indeed it is. This phase will be observed if the temperature of the system is smaller or comparable to the energy gap, . But as shown in figure 3 the energy gap can be considerably large, ∼ 0.6J, sufficiently larger than the temperature limitations of the latest experiments [15,41]. We expect that under these conditions the fragmented phase can be unambiguously identified.

Conclusions
In this work, we have introduced several new ideas. The first one is to use internal states of atoms in an optical lattice to simulate additional spatial dimension and to implement ladders. While it may be argued that there are previous works which have actively employed the internal degrees of freedom of the atom for quantum simulation [2,4,6,7,42,43], it is in this work that they are used to effectively increase the dimensionality of the lattice.
The second idea is to couple the atomic degrees of freedom in a spatially dependent way [42], so as to induce a frustrating hopping in the new virtual dimensions. The frustrating nature if this assisted tunnelling arises from the sign of the Raman coupling and it is best appreciated in the hard-core bosons limit, in which the model is tantamount to a spin Hamiltonian (section 2.2). Additionally, this frustrating hopping can also be seen as originating multiple small Mach-Zender interferometers in the lattice, such that destructive interference disturb the motion of atoms along the lattice (figure 4(e)).
The most dramatic consequence of these two kinds of tunnelling is the breakdown of superfluidity. The ground state fragments into a macroscopic number of double-well subsystems which lose coherence among them as we increase the strength of interactions. This situation is reminiscent of a Bose glass where a disordered energy landscape creates islands of superfluid regions with no coherence among them [24]. However, in our case the fragmentation is induced by frustration and not by disorder.
It is possible to qualitatively connect the physics observed here with that of other spin models. On the one hand, the modulation of the hopping leads to an effective superlattice in the internal plus motional degrees of freedom. This superlattice gives rise to a larger coherence between neighbouring lattice sites, similar to what happens in the spin-Peierls effect, where a modulation of a Heisenberg interaction leads to dimerization [44]. Thus, in contrast to other models where dimerization happens spontaneously [22], here it is imposed by the quantum interference of hopping terms.
However, while this analogy to spin models explains the fragmentation, our system goes beyond this simple picture, since we actually have two superlattices which are connected by repulsive interactions. The internal degrees of freedom of the atoms give rise to rich physics both in the strongly repulsive regime and in the insulator to superfluid transition. In particular, the dynamics of the spin deep in the superfluid region deserve tocc be studied in future work.
The methods shown here can be applied to simulating spin ladders as the ones studied in [7]. Additionally, the internal degrees of freedom of the atoms can be used to change the topology of the lattice from square to triangular. Finally, a natural extension would be to combine the Raman coupling used here with a 2D or even a 3D optical lattice, so as to simulate new topologies and other frustrated bilayer magnetic models that may have not been considered so far.