Research ArticlePHYSICS

Coupling ultracold matter to dynamical gauge fields in optical lattices: From flux attachment to ℤ2 lattice gauge theories

See allHide authors and affiliations

Science Advances  11 Oct 2019:
Vol. 5, no. 10, eaav7444
DOI: 10.1126/sciadv.aav7444


From the standard model of particle physics to strongly correlated electrons, various physical settings are formulated in terms of matter coupled to gauge fields. Quantum simulations based on ultracold atoms in optical lattices provide a promising avenue to study these complex systems and unravel the underlying many-body physics. Here, we demonstrate how quantized dynamical gauge fields can be created in mixtures of ultracold atoms in optical lattices, using a combination of coherent lattice modulation with strong interactions. Specifically, we propose implementation of ℤ2 lattice gauge theories coupled to matter, reminiscent of theories previously introduced in high-temperature superconductivity. We discuss a range of settings from zero-dimensional toy models to ladders featuring transitions in the gauge sector to extended two-dimensional systems. Mastering lattice gauge theories in optical lattices constitutes a new route toward the realization of strongly correlated systems, with properties dictated by an interplay of dynamical matter and gauge fields.


Gauge fields play a central role in a wide range of physical settings: The interactions in the standard model are mediated by gauge bosons, and everyday phenomena related to electromagnetism are governed by Maxwell’s equations featuring a gauge symmetry. The presence of strong magnetic fields can lead to strong alterations of the behavior of interacting many-body systems; for example, in the fractional quantum Hall (FQH) effect, the statistics of elementary excitations can be transmuted from fermionic to bosonic or to neither of both (anyonic) (1). Last, gauge theories even play a role in strongly correlated quantum systems, where local constraints lead to emergent gauge symmetries at low energies; for example, frustrated quantum spin liquids can be classified by their corresponding gauge theories.

The realization of artificial gauge fields in ultracold gases is an important milestone, enabling studies of the interplay between gauge fields and strong interactions in quantum many-body systems. This feat has further promoted these quantum-engineered systems as versatile quantum simulators (2, 3). While a synthetic magnetic field can be simply introduced by rotating atomic clouds (4, 5), more sophisticated schemes were developed to generate a wide family of gauge field structures, including spin-orbit couplings (6) or patterns featuring staggered magnetic fluxes with alternating signs on a length scale given by the lattice constant (79). The design of magnetic fluxes for ultracold atoms in optical lattices, through laser-induced tunneling or shaking methods, has been recently exploited in view of realizing topological states of matter (10, 11) and frustrated magnetism (8). In the settings described above, artificial gauge fields are treated as classical and nondynamical; however, in the sense that they remain insensitive to the spatial configuration and motion of the atomic cloud, these engineered systems do not aim to reproduce a complete gauge theory, where particles and gauge fields influence each other.

To be able to use ultracold atoms and simulate a wider range of physical problems, including those from high-energy physics, two major steps need to be taken. First, the synthetic gauge fields need to be made intrinsically dynamical, allowing back-actions of the particles on the gauge field. For example, the strength of the synthetic magnetic field may depend explicitly on the local particle density. In a second step, the dynamics of the synthetic gauge fields needs to be constrained to fulfill certain local symmetries. Therefore, the synthetic gauge field interacts with the matter particles, but each lattice site is also associated with a separately conserved charge. Theories of this type are called lattice gauge theories (LGTs), and the detailed conservation laws they satisfy depend on the respective gauge group (12). The simplest instant of an LGT has a ℤ2, or Ising, gauge group, but in the presence of fermionic matter, even this model poses a substantial theoretical challenge (13).

Various theoretical works have already suggested several methods by which synthetic gauge fields can be made intrinsically dynamical. A first approach builds on the rich interplay between laser-induced tunneling and strong on-site interactions, which can be both present and finely controlled in an optical lattice (3): Under specific conditions, the tunneling matrix elements, which not only describe the hopping on the lattice but also capture the presence of a gauge field, can become density dependent (1418); see (19) for an experimental implementation of these density-dependent gauge fields. While these settings include rich physics, they lack local conservation laws and thus still differ significantly from problems relevant to, e.g., high-energy physics.

A second approach aims at directly implementing genuine LGTs with local conservation laws, such as the Kogut-Susskind or quantum link models. This can be achieved, in principle, by engineering specific model Hamiltonians through elaborate laser-coupling schemes involving different atomic species and well-designed constraints; see (2022) for reviews and (23) for an ion trap realization of the Kogut-Susskind Hamiltonian. A digital implementation of ℤ2 LGTs, including couplings to fermionic matter, was proposed in (24). These quantum simulations of LGTs aim to deepen our understanding of fundamental concepts of gauge theories (2023, 25), such as confinement and its interplay with dynamical charges, which are central in high-energy (26) and condensed-matter physics (13, 27, 28) and go beyond a mere density dependence of synthetic gauge fields. While important first steps have been taken, the direct quantum simulation of LGTs is still in its infancy. The implementations proposed so far still contain significant technical challenges that need to be overcome, making alternative implementation schemes desirable.

In this work, we demonstrate how ℤ2 LGTs can be realized in ultracold gases through the use of specifically designed density-dependent gauge fields. Our approach combines the experimental advantages afforded by settings with density-dependent synthetic gauge fields and the additional physical structure added by the presence of local conservation laws. We demonstrate how existing ultracold atom technology can be used to implement toy models relevant to both high-energy and condensed-matter physics, and describe how the procedure can be scaled up to transition from simplistic two-site models to two-dimensional (2D) systems with direct relevance to studies of, e.g., high-temperature superconductivity (13).

As a central ingredient, we devise a scheme to engineer flux attachment for cold atoms moving in an optical lattice. Originally introduced by Wilczek (29, 30), and then widely exploited in the context of the FQH effect (1), flux attachment is a mathematical construction according to which a certain amount of magnetic flux quanta is attached to a particle (e.g., an electron). The resulting composite “flux tube particle” can change its quantum statistics from bosonic to fermionic, or vice versa (30), and naturally appears in field theoretical formulations of FQH states (1). Specifically, we show that an optical lattice loaded with two atomic species (a and f) can be configured in a way that a localized f-particle becomes a source of magnetic flux Φ for the a-particle: The magnetic flux can thus effectively be attached to the f-particles, which are also allowed to move around the lattice (see Fig. 1A). (Cases with only one species, more closely related to the FQH effect, can also be considered.) The flux attachment scheme is our starting point for implementing ℤ2 LGTs using ultracold atoms.

Fig. 1 Flux-attachment and dynamical gauge fields with ultracold atoms.

(A) We propose a setup where one atomic species f becomes a source of magnetic flux Φ (red) for a second species a. Both types of atoms undergo coherent quantum dynamics, described by NN tunneling matrix elements ta and tf, respectively. (B) When realized in a ladder geometry, the flux attachment setup has a ℤ2 lattice gauge structure. By tuning the ratio of the tunneling elements ta/tf, we find that the system undergoes a phase transition. The two regimes can be understood in terms of the elementary ingredients of a ℤ2 LGT, summarized in (C). The matter field a^ has a ℤ2 charge given by the parity of its occupation numbers n^a. It couples to the ℤ2 gauge field τ^i,jz, defined as the number imbalance of the f-particles between different ends of a link. When ∣ta∣≪∣tf∣, the ground state is dominated by tunneling of the f-particles, realizing that eigenstates of the ℤ2 electric field τ^i,jx delocalized over the link. In the opposite limit, ∣ta∣≫∣tf∣, the tunneling dynamics of the a-particles prevails and the system realizes eigenstates of the ℤ2 magnetic field B^p, defined as a product of the gauge field τ^z over all links ℓ ∈ ∂p along the edge of a plaquette p. The ℤ2 magnetic field introduces Aharonov-Bohm phases for the matter field, which are 0 (π) when the f particles occupy the same (different) leg of the ladder, i.e., if Bp = 1 (Bp = −1). The quantized excitations of the dynamical gauge field correspond to ℤ2 vortices of the Ising gauge field, so-called visons.

For specific choices of parameters and carefully designed lattice geometries, we first show that this appealing setting can be readily used to implement interacting quantum systems with ℤ2 link variables and global ℤ2 gauge symmetries. Then, we demonstrate that our method can also be extended to systems with local symmetries, realizing genuine ℤ2 LGTs (26) in various lattice geometries. These latter types of models, where the matter field couples to a ℤ2 lattice gauge field, are especially relevant in the context of high-Tc superconductivity (13, 31) and, more generally, strongly correlated electrons (3234). A central question in this context concerns the possibility of a confinement-deconfinement transition in the LGT (12), which would indicate electron fractionalization (13, 35, 36). The proposed model will allow us to explore the interplay of a global U(1) symmetry with local ℤ2 symmetries, which has attracted particular attention in the context of high-Tc cuprate compounds (34, 37, 38).

We discuss in detail the physics of a toy model characterized by a global U(1) × ℤ2 symmetry, which consists of a two-leg ladder geometry and can be directly accessed with state-of-the-art cold atom experiments. We demonstrate that the toy model features an intricate interplay of matter and gauge fields, as a result of which the system undergoes a phase transition in the ℤ2 sector depending on the ratio of the species-dependent tunnel couplings ta/tf (see Fig. 1B). While this transition can be characterized by the spontaneously broken global ℤ2 symmetry, we argue that an interpretation in terms of the constituents of a ℤ2 LGT (see Fig. 1C) is nevertheless useful to understand its microscopic origin. We also predict a phase transition of the matter field from an insulating Mott state to a gapless superfluid (SF) regime, associated with the spontaneously broken global U(1) symmetry. For appropriate model parameters, an interplay of both types of transitions can be observed.

The paper is organized as follows. We start by introducing the flux attachment scheme, which is at the heart of the proposed experimental implementation of dynamical gauge fields. Particular attention is devoted to the case of a double-well system, which forms the common building block for realizing ℤ2 LGTs coupled to matter. Next, we study the phase diagram of a toy model with a two-leg ladder geometry, consisting of a matter field coupled to a ℤ2 gauge field on the rungs. Realistic implementations of the considered models are proposed afterward, along with a scheme for realizing genuine ℤ2 LGTs with local instead of global symmetries in two dimensions. This paves the way for future investigations of strongly correlated systems, as discussed in the summary and outlook section.

The minimal model of a ℤ2 LGT coupled to matter proposed here has been realized experimentally in a double-well system (39). Besides, density-dependent Peierls phases have been realized with two-component fermions in (40), based on a two-frequency driving scheme, which is proposed below as an ingredient to implement ℤ2 LGTs in extended lattices.


Flux attachment

The recent experimental implementations of classical gauge fields for ultracold atoms (4145) combine two key ingredients (46): First, the bare tunnel couplings t are suppressed by large energy offsets ∣Δ∣≫ t, realized by a magnetic field gradient or a superlattice potential. Second, tunneling is restored with complex phases ϕ by proper time modulation of the optical lattice (47, 48) at the resonance frequency ω = Δ (with ħ = 1 throughout). The phase of the lattice shaking directly determines the value ϕ of the complex hopping element.

Flux attachment operates in a strongly correlated regime, where the energy offsets Δ = ω from an external potential are supplemented by interspecies Hubbard interactions of the same magnitude, U = ω (49). This provides coherent control over the synthetic gauge fields induced by the lattice modulation at frequency ω [see also (1418)].

We consider a situation where atoms of a first species, with annihilation operators a^, represent a matter field. The atoms of the second type, associated with annihilation operators f^, will become the sources of synthetic magnetic flux for the matter field (see Fig. 1A). Namely, the magnetic flux felt by the a-particle, as captured by its assisted hopping over the lattice, is only effective in the presence of an f-particle. To avoid that—vice versa—the f-particles become subject to magnetic flux created by the a-particles, static potential gradients affecting only the f-particles are used. In the following, we assume that both atomic species are hard-core bosons, although generalizations are possible, for instance, when one or both of them are replaced by fermions.

Model. The largest energy scale in our problem is set by strong interspecies Hubbard interactionsH^int=Ujn^jan^jf(1)where n^ja,f denote the density operators of a- and f-particles on lattice site j. To break the symmetry between a- and f-particles, we introduce state-dependent static potentials Vα(j), where α = a, f. We assume that the corresponding energy offsets between nearest-neighbor (NN) lattice sites i and j are integer multiples mi,jα of the large energy scale U, up to small corrections δVi,jαU, which are acceptable, namelyΔi,jαVα(i)Vα(j)mi,jαU(2)

A minimal example is illustrated in Fig. 2A.

Fig. 2 2 LGT in a two-well system.

(A) We consider a double-well setup with one atom of each type, a and f. Coherent tunneling between the two orbitals at j1 and j2 = j1 + ey is suppressed for both species by strong Hubbard interactions U = ω, and for f-particles by the energy offset Δf = ω indicated by the blue triangle. (B) Tunnel couplings can be restored by resonant lattice modulations with frequency ω. The sign of the restored tunneling matrix element is different when the a-particle gains (top) or loses (bottom) energy. (C) This difference in sign gives rise to a ℤ2 gauge structure and allows the implementation of ℤ2 minimal coupling of the matter field a^ to the link variable defined by the f-particles. The action of this term in the effective Hamiltonian acting on a basis state (left) is illustrated. The strength of the ℤ2 electric field is indicated by the thickness of the blue line connecting the two sites of the double-well system. The minimal coupling term is the common building block for realizing larger systems with a ℤ2 gauge structure. (D) These systems are characterized by a symmetry G^j associated with each lattice site j. Here, G^j commutes with the Hamiltonian and consists of the product of the ℤ2 charge, Q^j=(1)n^ja, and all electric field lines—for which τx = −1—emanating from a volume around site j.

Coherent dynamics of both fields are introduced by NN tunneling matrix elements in the μ = x, y directions, tμα, respectively. Thus, the free part of the Hamiltonian isH^0=μ=x,yi,jμ[tμaa^ja^i+tμff^jf^i+h.c.]+j[Va(j)n^ja+Vf(j)n^jf](3)where 〈i, jμ denotes a pair of NN sites along direction μ. Tunnel couplings are initially suppressed by the external potentials Δα = mαU and the strong Hubbard interactionsUtx,yα(4)

To restore tunnel couplings with complex phases, we include a time-dependent lattice modulationH^ω(t)=jVω(j,t)(a^ja^j+f^jf^j)(5)

It acts equally on both species and is periodic in time, Vω(j, t + 2π/ω) = Vω(j, t), with frequency ω = U resonant with the interspecies interactions. In summary, our Hamiltonian isH^(t)=H^0+H^int+H^ω(t)(6)

Effective hopping Hamiltonian. From now on, we consider resonant driving, U=ωtμα, where the lattice modulation H^ω(t) in Eq. 5 restores, or renormalizes, all tunnel couplings of a- and f-particles. As derived in Materials and Methods, we obtain an effective hopping Hamiltonian to lowest order in 1/ωH^eff=μ=x,yi,jμ[tμaa^ia^jλ^i,jμμ eiφ^i,jμμ+h.c.+tμff^if^jΛ^i,jμμ eiθ^i,jμμ+h.c.] (7)

The Hermitian operators λ^i,jμμ and φ^i,jμμ (Λ^i,jμμ and θ^i,jμμ, respectively) in Eq. 7 describe the renormalization of the tunneling amplitudes and phases for a (respectively f) particles; they are mutually commuting and depend only on the number imbalance n^jfn^if (respectively n^jan^ia) associated with the respective complementary species. Our result in Eq. 7 is reminiscent of the models discussed in (17).

Explicit expressions for λ^, φ^, Λ^, and θ^ can be obtained by considering their matrix elements on the relevant many-body states ∣ψr〉 and ∣ψs〉 in the Fock basis that are involved in the various hopping processes. For an a-particle transitioning from state ∣ψs〉 to ∣ψr〉, corresponding to a relative potential and/or interaction energy offset Δrs = nrsω with integer nrs ∈ ℤ, the matrix elements are given byψra^ia^jλ^i,jμμψs=Jnrs(x)(8)

Here 𝒥n denotes the Bessel function of the first kind, x = Ai,j/ω is the dimensionless driving strength, andVω(i,t)Vω(j,t)=Ai,jcos (ωt+ϕi,j)(9)

Without loss of generality, we assume ω, Ai,j > 0 throughout the paper.

The complex phases of the restored tunnelings are also determined by the many-body energy offsets Δrs = nrsω. If nrs ≥ 0, the particle gains energy in the hopping process andψra^ia^jφ^i,jμμψs=nrsϕi,j(10)

In contrast, if nrs < 0, the particle loses energy andψra^ia^jφ^i,jμμψs=nrs(πϕi,j)(11)

In this case, there is an additional nrsπ phase shift due to the reflection properties of the Bessel function, Jn(−x) = (−1)nJn(x) (see Fig. 2B). This nrs π phase shift is at the core of the LGT implementations discussed below. Similar results are obtained for Λ^i,jμμ and θ^i,jμμ by exchanging the roles of a and f (see the “Two-particle two-site problem” section in Materials and Methods). Note, however, that the symmetry between a and f can be broken by a careful design of the potentials Va and Vf, and this will be exploited in the next paragraph.

As illustrated in Fig. 1A, our scheme allows us to implement effective Hamiltonians (Eq. 7) describing a mixture of two species, where one acts as a source of magnetic flux for the other [see also (17)]. A detailed discussion of the resulting Harper-Hofstadter model with dynamical gauge flux is provided in section S1. By analogy with the physics of the FQH effect (50, 51), we expect that this flux attachment gives rise to interesting correlations and possibly to quasiparticle excitations with nontrivial statistics.

2 LGT in a double well. Now, we apply the result in Eq. 7 and discuss a minimal setting, where one a-particle and one f-particle each tunnel between the two sites j1 and j2 = j1 + ey of a double-well potential (see Fig. 2A); ey denotes the unit vector along y. This system forms the central building block for the implementation of ℤ2 LGTs in larger systems, proposed below. We assume that Va(ji) ≡ 0 for i = 1,2 but introduce a potential offset Vf(j2) = Δf + Vf(j1) for the f species, breaking the symmetry between a- and f-particles.

Effective Hamiltonian

For Δf = U = ω and lattice modulations with a trivial phase ϕj1,j2 = 0, the effective Floquet Hamiltonian in Eq. 7 becomesH^eff2well=tyaλy τ^j2,j1z(a^j2a^j1+h.c.)tyfΛ^τ^j2,j1x(12)with notations defined as follows. We describe the f-particle by a pseudo–spin-12τ^j2,j1z=n^j2fn^j1f,  n^j2f+n^j1f=1(13)which becomes a link variable in a ℤ2 LGT (see Fig. 1C). The Pauli matrix τ^j2,j1x=(f^j2f^j1+h.c.) describes tunneling of the f-particle.

As shown in Fig. 2B, the interaction energy of the matter field changes by ±U in every tunneling event. As a result, the amplitude renormalization in Eq. 12 is λy = ∣ 𝒥1(Aj2,j1/ω)∣ (see Eq. 8), and the phase of the restored tunnel couplings is eiφ^=τ^j2,j1z by Eqs. 10 and 11. Because the f-particle is subject to an additional potential offset Δf = U between the two sites, its energy can only change by 0 or 2U in a tunneling event. Hence, the phase of the restored tunneling in Eq. 12 is trivial, θ^=0 as in Eqs. 10 and 11, but the amplitude renormalizationΛ^=J0(Aj2,j1/ω)n^j1a+J2(Aj2,j1/ω)n^j2a(14)depends on the configuration of the a-particle in general.

The effective Hamiltonian (Eq. 12) realizes a minimal version of a ℤ2 LGT: The link variable τ^j2,j1zeiπ𝒜^ provides a representation of the dynamical ℤ2 gauge field 𝒜^, which is quantized to 0 and 1. The corresponding ℤ2 electric field is given by the Pauli matrix τ^j2,j1x, defining electric field lines on the link. The ℤ2 charges Q^ji, defined on the two sites ji with i = 1, 2, are carried by the a-particle, Q^ji=exp(iπn^jia). These ingredients are summarized in Fig. 1C and justify our earlier notion that the a- and f-particles describe matter and gauge fields, respectively. The Hamiltonian in Eq. 12 realizes a minimal coupling (12) of the a-particles to the gauge field (see Fig. 2C).


Each of the two lattice sites ji is associated with a ℤ2 symmetry. The operators generating the ℤ2 gauge group in the double-well systemg^i=Q^jiτ^j2,j1x,  i=1,2(15)both commute with the effective Hamiltonian in Eq. 12, [g^i,H^eff2well]=0 for i = 1,2. This statement is not entirely trivial for the first term in Eq. 12: While τ^j2,j1z and a^j2a^j1 do not commute with τ^j2,j1x and Q^ji individually, their product commutes with g^i. The second term in Eq. 12 trivially commutes with g^i because [Λ^,Q^ji]=0 (see Eq. 14).

Physically, Eq. 15 establishes a relation between the ℤ2 electric field lines, τj2,j1x=1, and the ℤ2 charges from which they emanate (see Fig. 2D). Note that the eigenvalues of g^1 and g^2 are not independent, because g^1g^2=1 for the considered case with a single a-particle tunneling in the double-well system.

The model in Eq. 12 is invariant under the gauge symmetries g^i for all values of the modulation strength Aj2,j1. In general, both terms in the effective Hamiltonian couple the ℤ2 charge to the gauge field. An exception is obtained for lattice modulation strengths Aj2,j1/ω = x02 for whichJ0(x02)=J2(x02)(16)

In this case, neither of the amplitude renormalizationsΛ^Λ02=J0(x02)0.32(17)λy=λ02=J1(x02)0.58(18)is operator valued, and the second term in the Hamiltonian only involves the ℤ2 gauge field. The weakest driving for which Eq. 16 is satisfied has x02 ≈ 1.84.


We take this opportunity to explain in a mechanistic way the meaning of the ℤ2 gauge field in the double-well system and its relation to more general LGTs. As a starting point, consider the situation when tyf=0 and the f-particle is localized. Depending on the position of the f atom (left or right), the restored tunneling amplitude of the a-atom between the two sites has a sign ±1. Formally, this corresponds to the appearance of the factor τ^j2,j1z=n^j2fn^j1f in the first term on the left hand side of Eq. 12. At this point, we have realized a synthetic gauge field 𝒜j2,j1 for the a-atoms, which depends on the density of the f-atoms. That is, the general tunneling matrix element of the a-particles, tyaeiπ𝒜j2,j1, has a phase π𝒜j2,j1〉 depending on the f-density.

The two possible states of the link variable, corresponding to the two positions j1 and j2 of the f atom, define a 2D Hilbert space on the link 〈j2, j1〉. This Hilbert space is equivalent to the Hilbert space of a ℤ2 lattice gauge field, with two orthogonal states on each link of the lattice. Using this language, we can identify the operator τ^j2,j1z with a ℤ2 gauge field. It does not commute with the corresponding ℤ2 electric field operator, τ^j2,j1x, which corresponds physically to the coherence of the f-atom between the two sites j2 and j1. This noncommutativity is related to the noncommutativity of the conjugate electric and magnetic fields e^ and B^ in quantum electrodynamics. Because of the small ℤ2 gauge group, the ℤ2 electric field only takes two possible quantized values: τ^j2,j1x=±1. The eigenstates correspond to even and odd superpositions of the f-atom on the two lattice sites: (f^j2±f^j1)0/2.

If we allow to add arbitrary perturbations to the tunneling Hamiltonian, e.g., terms δyτ^j2,j1y or δzτ^j2,j1z, we can introduce nontrivial dynamics of the synthetic gauge field. While this renders the density-dependent synthetic gauge field dynamical, it does not correspond to a ℤ2 LGT in the strict sense, since local conservation laws are generically absent. These situations, without local conservation laws, lead to interesting physics nonetheless and have been studied for example in the context of the so-called ℤ2 Bose-Hubbard model (52, 53).

A genuine ℤ2 LGT is obtained if only terms are included in the Hamiltonian, which commute with the ℤ2 gauge operators g^i in Eq. 15. In particular, this includes the term tyfΛ^τ^j2,j1x on the right hand side of Eq. 12 or simpler terms such as tyfτ^j2,j1x. These are the ℤ2 analogs of terms ∼E2 in the Hamiltonian of quantum electrodynamics, without the square due to the simpler nature of the ℤ2 gauge group. In the double-well system, the ℤ2 gauge symmetry leads to two decoupled sectors of the Hamiltonian with g1 = −g2 = 1 and g1 = −g2 = −1. As will be shown later, however, in extended systems, each lattice site is associated with its own conserved charge. This has important consequences for the possible many-body phases (12).

Matter gauge field coupling in two-leg ladders

In the following, we study the physics of coupled matter and gauge fields in a two-leg ladder, accessible with numerical density matrix renormalization group (DMRG) simulations (54). Our starting point is a model with minimal couplings to the ℤ2 gauge field on the rungs of the ladder, which is characterized by a global U(1) × ℤ2 symmetry (see Fig. 3A). Here, we study its phase diagram. As explained later, the model can be implemented relatively easily in existing ultracold atom setups by coupling multiple double-well systems, which is our main motivation for studying its phase diagram.

Fig. 3 Coupling matter to a ℤ2 gauge field in a two-leg ladder.

(A) We consider the Hamiltonian (Eq. 24) describing a-particles that are minimally coupled to the ℤ2 gauge field τ^i,jyz on the rungs of a two-leg ladder. (B) The phase diagram, obtained by DMRG simulations at tyf/txa=0.54, contains an SF-to-Mott transition in the charge sector at a commensurate density of the matter field, Na = Lx. In addition, we find a transition in the gauge sector, from an ordered region with a broken global ℤ2 symmetry where the ℤ2 magnetic field dominates and the vison excitations of the gauge field are gapped (red) to a disordered regime where the ℤ2 electric field is dominant and visons are strongly fluctuating in a condensed state (blue). Along the hatched lines at commensurate fillings Na/Lx = 1/2,1, 3/2, insulating CDW states could exist, but conclusive numerical results are difficult to obtain. (C) The conjectured schematic phase diagram of Eq. 24 is shown in the μtyf plane, where μ denotes the chemical potential for a^ particles and 2tyf corresponds to the energy cost per ℤ2 electric field line along a rung. Two scenarios are realized in different parameter regimes: In scenario I, the interplay of gauge and matter fields prevents a fully disordered Mott phase, whereas the latter exists in scenario II. The behavior in scenario I resembles the phase diagram of the more general 2D ℤ2 LGT (13, 2628) sketched in (D). In our DMRG simulations here, as well as in the following figures, we keep up to 1400 DMRG states with five finite-size sweeps; the relative error on the energies is kept smaller than 10−7.

The model. We combine multiple double-well systems (Eq. 12) to a two-leg ladder by introducing tunnelings txa of the matter field along x. Furthermore, we impose that the f-particles can only move along the rungs, txf=0, and each rung is occupied by one f-particle. Thus, we can continue describing the f degrees of freedom by link variables τ^i,jy, as defined in Eq. 13. The number of a-particles Na will be freely tunable.

Effective Hamiltonian

For a properly designed configuration of lattice gradients and modulations, presented in detail later, we obtain an effective HamiltonianH^2leg=i,jx(txaλ^i,jxxa^ja^i+h.c.)i,jy[tyaλy(a^ja^iτ^i,jyz+h.c.)+tyfΛ^i,jyyτ^i,jyx](19)

Expressions for the amplitude renormalizations λy ∈ ℝ and λ^x, Λ^y are provided in section S2B.

For the specific set of driving strengths x = x02 that we encountered already in the double-well problem (see Eq. 16), we find that Λ^y only has a weak dependence on the ℤ2 charges, Q^j=(1)n^ja. Similarly, the amplitude renormalization λ^x depends weakly on the ℤ2 magnetic field B^p only; hereB^p=i,jypτ^i,jyz(20)

is defined as a product over all links 〈i, jy on the rungs belonging to the edge ∂p of plaquette p. Hence, for these specific modulation strengths[Λ^i,jyy,Q^l]=[Λ^i,jyy,τ^k,l]=0(21)[λ^i,jxx,Q^p]=[λ^i,jxx,a^l()]=0(22)


Now, we discuss the symmetries of the effective Hamiltonian (Eq. 19) at the specific value of the driving strengths x02. In the case of decoupled rungs, i.e., for txa=0, every double well commutes with g^i, i = 1, 2 from Eq. 15. These symmetries are no longer conserved for txa0; in this general case, a global ℤ2 gauge symmetry remainsV^i=j=1Lxg^i(jex),  i=1,2(23)

with g^i(jex)=(1)Q^jex+(i1)eyτ^jex+ey,jexyx and for which V^i2=1. Using Eqs. 21 and 22, one readily confirms that [H^2leg,V^i]=0 for i = 1, 2.

In summary, the effective model is characterized by the global U(1) symmetry associated with the conservation of the number of a-particles and the global ℤ2 symmetry V^1. Note that the second ℤ2 symmetry, V^2, follows as a consequence of combining V^1 with the global U(1) symmetry: By performing the global U(1) gauge transformation a^ja^j for all sites j, V^2 is obtained from V^1. Thus, the overall symmetry is U(1) × ℤ2.

Physical constituents

In the following, we will describe the physics of the ladder models using the ingredients of ℤ2 LGTs (see Fig. 1C). The quantized excitations of the ℤ2 lattice gauge field are vortices of the ℤ2 (or Ising) lattice gauge field, so-called visons (13). They are defined on the plaquettes of the ladder: If the plaquette term in Eq. 20 is Bp = 1, there is no vison on p; the presence of an additional ℤ2 flux, Bp = −1, corresponds to a vison excitation on plaquette p. Since the matter field a^ couples to the ℤ2 gauge field, the resulting interactions with the visons determine the phase diagram of the many-body Hamiltonian.

Quantum phase transitions of matter and gauge fields. We start from the microscopic model in Eq. 19 and simplify it by making a mean field approximation for the renormalized tunneling amplitudes, which depend only weakly on Q^j and B^p. Replacing them by ℂ-numbers, txa=txaλ^x, tyf=tyfΛ^y, and tya=tyaλy leads to the conceptually simpler HamiltonianH^2legsimp=i,jxtxa(a^ja^i+h.c.)i,jy[tya(a^ja^iτ^i,jyz+h.c.)+tyfτ^i,jyx](24)

illustrated in Fig. 3A. Later, by introducing a more sophisticated driving scheme, we will show that this model can also be directly implemented using ultracold atoms. The simpler Hamiltonian (Eq. 24) has identical symmetry properties as Eq. 19. Now, we analyze Eq. 24 by means of the DMRG technique. In the phase diagram, we find at least three distinct phases, resulting from transitions in the gauge and matter field sectors (see Fig. 3B). Here, we describe their main features; for more details, the reader is referred to section S4.

Transition in the matter sector

First, we concentrate on the conceptually simpler phase transition taking place in the charge sector. When the tunneling along the legs is weak, txa[(tyf)2+(tya)2]1/2tyf, and the number Na of a-particles is tuned, we observe a pronounced rung-Mott phase (55) at the commensurate filling Na = Lx, where Lx denotes the total number of rungs in the system. Similar to the analysis in (5658), this phase can be characterized by the parity operatorOp(l)=exp [iπj<l(n^jexa+n^jex+eyaNaLx)](25)

In the limit lLx and Lx → ∞, the observable Op(l) remains finite only in the Mott insulating regime. Our results in Fig. 4A confirm that Op takes large values with a weak size dependence for Na = Lx. On the other hand, when Na/Lx ≠ 1 is slightly increased or decreased, the parity Op immediately becomes significantly smaller and a more pronounced Lx dependence is observed, consistent with a vanishing value in thermodynamic limit as expected for an SF phase.

Fig. 4 Characterizing phase transitions of matter coupled to a ℤ2 gauge field in a two-leg ladder.

We consider the Hamiltonian from Eq. 24. (A) In the charge sector, we observe transitions from an SF state, characterized by a vanishing parity correlator Op(Lx/2 → ∞ ) → 0 in the thermodynamic limit, to an insulating rung-Mott state at the commensurate filling Na = Lx, characterized by Op(Lx/2 → ∞ ) > 0 and exponentially decaying correlations. Here, we present exemplary results for tya/txa=3 and tyf/txa=0.54. (B) In the gauge sector, we find a transition from a disordered phase, where the ℤ2 electric field dominates, to a phase where the ℤ2 magnetic field dominates. In the second case, the order parameter τ^i,jyz0 corresponds to a spontaneously broken global ℤ2 symmetry (23). In the two phases, the corresponding vison excitations of the ℤ2 gauge field (C) have different characteristics. The numerical results in (A) [respectively (B)] are obtained by considering periodic boundary conditions (respectively Lx = 96 rungs with open boundaries). (D) Analogs of Wilson loops W^(d) in the two-leg ladder are string operators of visons.

For larger values of txa and Na = Lx, no clear signatures of a rung-Mott phase are found (see section S4B): The parity operator Op takes significantly smaller values, the calculated Mott gap becomes a small fraction of txa, consistent with a finite-size gap, and we checked that the decay of two-point correlations follows a power-law at long distances until edge effects begin to play a role. These signatures are consistent with an SF phase.

Because of the global U(1) symmetry of the model, an SF-to-Mott transition at Na = Lx in the quasi-1D ladder geometry would be of Berezinskii-Kosterlitz-Thouless (BKT) type (59). Hence, the gap is strongly suppressed and the correlation length is exponentially large close to the transition point, making it impossible to determine conclusively from our numerical results whether the ground state is a gapped Mott state in this regime. For hard-core bosons on a two-leg ladder, this scenario is realized: It has been shown by bosonization that an infinitesimal interleg coupling is sufficient to open up an exponentially small Mott gap (55, 60).

Similar considerations apply at the commensurate fillings Na/Lx = 1/2, 3/2, where previous work on single-component bosons in a two-leg ladder (60) pointed out the possible emergence of an insulating charge density wave (CDW) for large enough tya/txa (hatched areas in Fig. 3B). Our numerical results indicate that these CDWs may exist in this regime also in our model (Eq. 24). But since the numerical analysis is plagued by a potentially even larger correlation length and a corresponding strongly suppressed Mott gap, it is difficult to pinpoint the exact location of this BKT transition.

The SF phase observed at incommensurate filling fractions is characterized by a power-law decay of the Green’s function a^dexa^0d1/4K in the charge sector. The exponent is related to the Luttinger parameter K, which approaches K → 1/2 at a transition to the rung-Mott phase for commensurate filling. We confirm this behavior and obtain qualitatively identical results as in the case of a static gauge field (see section S4A for details) (55).

Transition in the gauge sector

In the gauge sector, described by f-particles, we observe a phase transition when the ratios of the tunnel couplings are tuned. In Fig. 4B, we tune tya/txa while keeping tyf/txa fixed. We find a transition from a symmetric regime where τ^i,jyz=0 to a region with a nonvanishing order parameter τ^i,jyz0. Similar behavior is obtained when tuning tyf/tya while keeping txa/tya fixed (see section S4E).

The observed transition is associated with a spontaneous breaking of the global ℤ2 symmetry (Eq. 23) of the model. The f-particles go from a regime where they are equally distributed between the legs, τ^i,jyz=0, to a twofold degenerate state with population imbalance, τ^i,jyz0. This behavior occurs in the insulating and SF regimes of the charge sector, and it is only weakly affected by the filling value Na/Lx (see Fig. 3B).

The two phases are easily understood in the limiting cases. When tya=0, the ground state is an eigenstate of the ℤ2 electric field τ^i,jyx on the rungs, with eigenvalues 1. The ℤ2 magnetic field is strongly fluctuating, and no ℤ2 electric flux loops exist. Thus, also the vison number is strongly fluctuating, and the state can be understood as a vison condensate. In the opposite limit, when tya, the kinetic energy of the matter field dominates. In this case, the ℤ2 magnetic field is effectively static, and its configuration is chosen to minimize the kinetic energy of the a-particles. This is achieved when the effective Aharonov-Bohm phases on the plaquettes vanish, i.e., for B^p=1 (see Fig. 1C). In this case, vison excitations with B^p=1 (see Fig. 4C) correspond to localized defects in the system, which cost a finite energy corresponding to the vison gap.

LGTs are characterized by Wilson loops (12). Their closest analogs in our two-leg ladder model are string operators of visonsW(d)=j=1dB^pj=τ^i,jyzτ^i+dex,j+dexyz(26)

see Fig. 4D. In the disordered phase (electric field dominates), we found numerically that W(d) → 0 when d → ∞. Sufficiently far from the transition, where finite-size effects are small, our data (see fig. S5) show the expected exponential decay. In the ordered phase (magnetic field dominates), W(d) → W quickly converges to the nonzero value W=τ^i,jz2>0 when d → ∞. See section S3 (D and E) for more details.

The qualitatively different behavior of the Wilson loop in the two phases is reminiscent of the phenomenology known from the ubiquitous confinement-deconfinement transitions found in (2 + 1) dimensional LGTs (12, 13, 26): There, visons are gapped in the deconfined phase and the Wilson loop decays only weakly exponentially with a perimeter law; in the confining phase, visons condense and the Wilson loop decays much faster with an exponential area law.

Our numerical results in Fig. 4B indicate that the phase transition in the gauge sector is continuous. Beyond this fact, its nature is difficult to determine. The interactions between the ℤ2 link variables are mediated by the matter field, which has correlations extending over many sites following either a power-law (in the SF regime) or featuring exponential decay with a correlation length ξ ≫ 1 (in the considered Mott-insulating regions). On the one hand, this leads to relatively large finite-size effects in our numerical simulations, which explains the continuous onset of the transition in Fig. 4B. On the other hand, when nonlocal Ising interactions mediated by the gapless matter field compete with the transverse ℤ2 electric field term tyfτ^i,jx in the Hamiltonian, a rich set of critical exponents can be expected (61). For more details, see section S4E.

Interplay of matter and gauge fields

Last, we discuss the interplay of the observed phase transitions in the gauge and matter sectors. To this end, we find it convenient to consider the phase diagram in the μtyf plane, where μ denotes a chemical potential for the a-particles and tyf controls fluctuations of the ℤ2 electric field. We collect our result in the schematic plots in Fig. 3C: Deep in the SF phase, realized for small μ and Na/Lx < 1, tyf drives the transition in the gauge sector. Because of a particle-hole symmetry of the hard-core bosons in the model, similar results apply for large μ and Na/Lx > 1. On the other hand, when tyf is small, permitting a sizable Mott gap at commensurate fillings, μ drives the SF-to-Mott transition.

More interesting physics can happen at the tip of the Mott lobe for commensurate fillings Na = Lx. This corresponds to the hatched regime in Fig. 3B, where we cannot say conclusively if the system is in a gapped Mott phase. To obtain better understanding of the commensurate regime, we first argue that an SF cannot coexist with the ordered phase of the gauge field at commensurate fillings: In this regime, the ℤ2 gauge field acquires a finite expectation value, τi,jyz0. This leads to a term in the Hamiltonian tyaτi,jyza^ia^j, which is expected to open a finite Mott gap, following the arguments in (55, 60). Therefore, only the two scenarios shown in Fig. 3C are possible: In the first case, the Mott insulator coexists only with the ordered phase of the gauge field; in the second scenario, the Mott state coexists with the disordered phase of the gauge field.

To shed more light on this problem, we consider the case when the Mott gap Δ is much larger than the tunneling txa. When txa=0, every rung represents an effective localized spin-1/2 degree of freedom. As shown in section S4C, finite tunnelings txaΔ introduce antiferromagnetic couplings between these localized moments, and in this limit, our system can be mapped to an XXZ chain. It has an Ising anisotropy, and the ground state has a spontaneously broken ℤ2 symmetry everywhere, except when tyf/tya, where an isotropic Heisenberg model is obtained and the ground state has power-law correlations. The transition from the gapped Mott state, corresponding to the ordered phase of the ℤ2 gauge field, to a symmetric state of two decoupled SFs with a disordered gauge field is of BKT type (62).

Our last argument demonstrates that scenario I in Fig. 3C is realized deep in the Mott phase. In this limit of small txa and strong couplings tyf of the gauge field, our analysis proves that an intricate interplay of the phase transitions in the gauge and matter sectors exists. This behavior, characteristic for scenario I in Fig. 3C, is reminiscent of the phase diagram of the 2D ℤ2 LGT (see Fig. 3D) (26, 28). In that case, the phase at weak couplings has topological order as in Kitaev’s toric code (28), and the disordered phases are continuously connected to each other at strong couplings.

On the other hand, a detailed analysis of the Luttinger-K parameter for larger values of txa shows that the ground state at commensurate filling Na = Lx is characterized by K = 1/2, in both the ℤ2 symmetric and ℤ2 broken regimes (see section S4, A and B). This behavior is indicative of scenario II in Fig. 3C, since the BKT transition is characterized by a value K = 1 of the Luttinger parameter (62).

Implementation: Coupled double-well systems

Now, we describe how the models discussed above, and extensions thereof, can be implemented in state-of-art ultracold atom setups. The double-well system introduced around Eq. 12 constitutes the building block for implementing larger systems with a ℤ2 gauge symmetry, or even genuine ℤ2 LGTs, because it realizes a minimal coupling of the matter field to the gauge field (see Fig. 2C) (12). We start by discussing the two-leg ladder Hamiltonian H^2leg (Eq. 19); then, we present a scheme, based on flux attachment, for implementing a genuine ℤ2 LGT coupled to matter in a 2D square lattice.

Two-leg ladder geometry. The ladder system shown in Fig. 3A can be obtained by combining multiple double wells (Eq. 12) and introducing tunnelings txa of the matter field along x, while txf=0. The lattice potential is modulated along y with amplitude Aj2,j1=Vωy, as in the case of a single double well. As described in Fig. 5, we introduce an additional static potential gradient with strength Δxa=U=ω per lattice site along x and modulate it with frequency ω and amplitude Vωx.

Fig. 5 Implementing matter-gauge field coupling in a two-leg ladder.

Multiple double-well systems as described in Fig. 2 are combined to form a two-leg ladder by including hopping elements txa of the a-particles along the x direction. Coherent tunneling is first suppressed by strong interspecies Hubbard interactions U and static potential gradients: Δxa=U for a-particles along x, and Δyf=U for f-particles along y. These gradients are indicated by triangles whose colors refer to the respective atomic species. The tunnel couplings are restored by a resonant lattice shaking with frequency ω = U, realized by a modulated potential gradient Vω(j,t)=(jxVωx+jyVωy)cos (ωt) seen by both species. The modulated gradients are indicated by light-colored triangles. We assume that each rung is occupied by exactly one f-particle, which can thus be described by a link variable, while the number Na of a-particles is freely tunable. As shown in section S2, the special choice for the driving strengths Vωx/ω=Vωy/ω=x02 leads to an effective Hamiltonian with matter coupled to ℤ2 lattice gauge fields on the rungs. The gradient Δxa=U guarantees that the a-particles pick up only trivial phases φ^x=0 while tunneling along the legs of the ladder. Hence, the Aharonov-Bohm phases (red) associated with the matter field become 0, or π corresponding to a vison excitation. They are determined by the plaquette terms B^p defined in Eq. 20, reflecting the configuration of f-particles.

As shown in section S2, this setup leads to the effective Hamiltonian (Eq. 19). For the specific set of driving strengths Vωx/ω=Vωy/ω=x02 (see Eq. 16), the amplitude renormalizations are λy = λ02 andλ^i,jxx=12(1τ^i±ey,izτ^j±ey,jz) J0(x02)+12(1+τ^i±ey,izτ^j±ey,jz)J1(x02)(27)

Simplified model

Now, we discuss a further simplification of the model in Eq. 19, leaving its symmetry group unchanged. We note that, even for the specific choice of the driving strengths Vωx/ω=Vωy/ω=x02, the renormalized tunnel couplings of the a-particles along x still depend explicitly on the ℤ2 gauge fields on the adjacent rungs (see Eq. 27). This complication can be avoided, by simultaneously modulating the gradient along x at two frequencies, ω and 2ω, with amplitudes Vωx and V2ωx; i.e., we consider the following driving term in Eq. 5Vω(j,t)=[jxVωx+jyVωy]cos(ωt)+jxV2ωxcos(2ωt)(28)

Following (48), we obtain expressions for the restored tunnel couplings along x for an energy offset nω introduced by the Hubbard interactions U = ω between a- and f-particles; λn==Jn2(x(1))J(x(2)/2), where x(1)=Vωx/ω and x(2)=V2ωx/ω (see section S2). By imposing the conditions λ0 = λ1 = λ2, we obtain a simplified effective Hamiltonian where λ^i,jxxλx is no longer operator valued, and thus completely independent of the ℤ2 gauge field τ^z. The weakest driving strengths for which this condition is met are given byx(1)=x012(1)1.71,  x(2)=x012(2)1.05(29)

where λx = λ012 ≈ 0.37. A similar approach can be used to make Λ^y independent of the ℤ2 charges, which allows the implementation of H^2legsimp from Eq. 24.

Realizing a2 LGT in a 2D square lattice. Now, we present a coupling scheme of double-wells, which results in an effective 2D LGT Hamiltonian with genuine local symmetries, in addition to the global U(1) symmetry associated with a-number conservation. We will derive a model with ℤ2 gauge-invariant minimal coupling terms τ^i,jza^ja^i along all links of the square lattice.


We consider the setup shown in Fig. 6A in a layered 2D optical lattice, which is a particular type of brick-wall lattice. The a-particles tunnel vertically between the layers in the z direction, with coupling matrix element tza, and along the links indicated in the figure with tunnel couplings txa and tya. Every tube consisting of four lattice sites with coordinates x, y, and nez for n = 1, 2, 3, 4 defines a supersite j = xex + yey in the effective 2D lattice shown in Fig. 6B. The four links connecting every supersite to its nearest neighbors i : 〈i, j〉 are realized by double-well systems, with exactly one f-particle each, in different layers of the optical lattice. The f-particles are only allowed to tunnel between the sites of their respective double-wells in the x-y plane, with amplitudes txf and tyf, while tunneling along the z direction is suppressed, tzf=0.

Fig. 6 Realizing ℤ2 LGT coupled to matter in 2D.

(A) Multiple double-well systems as described in Fig. 2 are combined in the shown brick-wall lattice. Each of its four layers along the z direction is used to realize one of the four links connecting every lattice site of the 2D square lattice (B) to its four nearest neighbors. The double-well systems are indicated by solid lines (colors), and they are only coupled by tunnelings of a-particles along the z direction, with amplitudes tza. The required lattice gradients (their modulations) are indicated by (light) colored triangles. (B) The restored hopping Hamiltonian H^2DLGT in the 2D lattice has local symmetries G^j associated with all lattice sites j, i.e., [H^2DLGT,G^j]=0.

For the realization of the individual double-well systems, we consider a modulated potential gradient along x and y, seen equally by the matter and gauge fields. The modulation amplitudes Vωx/ω=Vωy/ω=x02 are chosen to simplify the amplitude renormalization of f-particle tunneling. As previously, we consider static potential gradients along x and y directions of Δxf=Δyf=U per site, seen only by the f-particles, and work in a regime where U=ωtμν, with μ = x, y, z and ν = a, f.

To realize a-particle tunneling along z, which is independent of the ℤ2 gauge fields τ^z on the links in the x-y plane, we add a static potential gradient of Δza per site along the z direction. It is modulated by two frequency components ω and 2ω, with amplitudes Vωz and V2ωz. These driving strengths are chosen as in Eq. 29, i.e., Vωz/ω=x012(1) and V2ωz/ω=x012(2), such that the restored tunnel couplings with amplitude tzaλ012 become independent of the f-particle configuration.

Effective Hamiltonian

Combining our results from the previous section, we obtain the effective hopping Hamiltonian H^2DLGT for the setup described in Fig. 6H^2DLGT=txyfi,jΛ^i,jτ^i,jxtxyaλ02i,j(τ^i,jza^i,mi,ja^j,mi,j+h.c.)tzaλ012jn=13(a^j,n+1a^j,n+h.c.)(30)

using the same notation as introduced earlier. Here, we treat the z-coordinate nez, with n = 1, …,4, as an internal degree of freedom, while j is a site index in the 2D square lattice; mi,j ∈ {1, 2, 3, 4} denotes the z-coordinate corresponding to double well 〈i, j〉. For simplicity, we assumed that txa=tya=txya and txf=tyf=txyf. The amplitude renormalization for f-particles in the x-y plane depends on the ℤ2 charges Q^j,nΛ^i,j=12[J0(x02)+J1(x02)]+Q^i,mi,jQ^j,mi,j12[J1(x02)J0(x02)](31)

Using the multifrequency driving scheme explained around Eq. 28, a situation where Λ^i,j becomes independent of the ℤ2 charges can be realized.

A simplified effective Hamiltonian, where the internal degrees of freedom are eliminated, can be obtained when U=ωtza and λ012tzatxya; the first inequality is required by the proposed implementation scheme. In this limit, the tunneling of a-particles along z can be treated independently of the in-plane tunnelings txya. The ground state with a single a-particle tunneling along z at supersite j is a^j0, where a^j=n=14ϕna^j,n, with ϕ1=ϕ4=(5+5)1/2 and ϕ2=ϕ3=(1+1/5)1/2/2. It is separated by an energy gap Δε=λ012tzatxya from the first excited state, which justifies our restriction to this lowest internal state.

The ground state energy ε2a with two hard-core a-particles tunneling along z in the same supersite is larger than twice the energy εa of a single a-particle, by an amount Ueff, i.e., ε2a = 2εa + Ueff. By solving the one- and two-particle problems exactly, we find Ueff=λ012tza. In the effective model restricted to the lowest internal state, this offset corresponds to a repulsive Hubbard interaction on the supersites j. Because Uefftxya, double occupancy of supersites is strongly suppressed, and we can treat the new operators a^j() as hard-core bosons.

By projecting the Hamiltonian (Eq. 30) to the lowest internal state on every supersite, we arrive at the following simplified modelH^2DLGTsimp=εaja^ja^jtxyfi,jΛ^i,jτ^i,jxtxyaλ02ϕ12i,jE(τ^i,jza^ia^j+h.c.)txyaλ02ϕ22i,jB(τ^i,jza^ia^j+h.c.)(32)

Here, we distinguish between two sets of links, 〈i,j〉 ∈ E or B, which are realized in layers at the edge n = 1,4 (E) and in the bulk n = 2,3 (B) in the 3D implementation (see Fig. 6A). Because the internal state has different weights ∣ϕ12 ≈ 0.14 and ∣ϕ22 ≈ 0.36, they are associated with different tunneling amplitudes. This complication can be avoided by realizing bare tunnelings of a-particles with different strengths on E- and B-type bonds.


In contrast to the two-leg ladder (Eq. 19), the models in Eqs. 30 and 32 are both characterized by local ℤ2 gauge symmetries. The ℤ2 charge on a supersite is defined as Q^j=exp [iπn=14n^j,na], which becomes Q^j=exp [iπa^ja^j] when projected to the lowest internal state. The ℤ2 gauge group is generated byG^j=Q^ji:j,iτ^j,ix(33)

where the product on the right includes all links 〈j, i〉 connected to site j.

It holds [H^2DLGT,G^j]=0 and [H^2DLGTsimp,G^j]=0 for all j, using the respective ℤ2 charge operators. These results follow trivially for the first line of Eqs. 30 and 32, which contain only the operators τ^i,jx and n^j,na (n^j) (see also Eq. 31). For the last two lines in the effective Hamiltonians, it is confirmed by a straightforward calculation.

In addition to the local ℤ2 gauge invariance, the models (Eqs. 30 and 32) have a global U(1) symmetry associated with the conservation of the a-particle number. Very similar Hamiltonians have been studied in the context of strongly correlated electrons, where fractionalized phases with topological order have been identified (35). When the a-particles condense, effective models without the global U(1) symmetry can also be realized. These are in the same symmetry class as Kitaev’s toric code (28).


We have presented a general scheme for realizing flux attachment in 2D optical lattices, where one species of atoms becomes a source of magnetic flux for a second species. For a specific set of parameters, we demonstrated that the effective Floquet Hamiltonian describing our system has a ℤ2 gauge structure. This allows us to implement experimentally a dynamical ℤ2 gauge field coupled to matter using ultracold atoms, as we have shown specifically for a double-well setup, two-leg ladders, and in a 2D geometry. Because our scheme naturally goes beyond one spatial dimension, the ℤ2 magnetic field—and the corresponding vison excitations—plays an important role in our theoretical analysis of the ground-state phase diagram. Moreover, the link variables in our system are realized by particle number imbalances on neighboring sites, making experimental implementations of our setup feasible using existing platforms [as described, e.g., in (42, 43, 45)].

Our theoretical analysis of hard-core bosons coupled to ℤ2 link variables on the rungs of a two-leg ladder revealed an SF-to-Mott transition in the charge sector as well as a transition in the gauge sector. The latter is characterized by a spontaneously broken global ℤ2 symmetry, but we argued that it can be considered as a precursor of the confinement-deconfinement transitions, which are ubiquitous in LGTs, high-energy physics, and strongly correlated quantum many-body system. Leveraging the powerful toolbox of quantum gas microscopy, our approach paves the way for new studies of LGTs with full resolution of the quantum mechanical wave function. This is particularly useful for analyzing string (57, 63) and topological (64) order parameters, which are at the heart of LGTs but difficult to access in more conventional settings.

As we have demonstrated, extensions of our LGT setting to 2D systems with local rather than global symmetries are possible. Here, we propose a realistic scheme to implement a genuine ℤ2 LGT with minimal coupling of the matter to the gauge field on all links of a square lattice. On the one hand, this realizes one of the main ingredients of Kitaev’s toric code (28, 65, 66)—a specific version of an LGT coupled to matter, which displays local ℤ2 gauge symmetry and hosts excitations with non-Abelian anyonic statistics. On the other hand, the systems that can be implemented with our technique are reminiscent of models studied in the context of nematic magnets (27, 33, 67) and strongly correlated electron systems (13, 35, 36). Other extensions of our work include studies of more general systems with flux attachment, which are expected to reveal physics related to the formation of composite fermions in the FQH effect.

Another application of our work is the realization of the recently suggested ℤ2 Bose-Hubbard Hamiltonian (68) using ultracold atoms in optical lattices. This model contains ℤ2 link variables on a 1D chain, similar to our case, but includes terms in the Hamiltonian, which explicitly break the local ℤ2 gauge symmetry. In contrast to the models studied in this paper, not only the tunneling phases but also the tunneling amplitudes in the ℤ2 Bose-Hubbard Hamiltonian depend on the ℤ2 link variables. The ℤ2 Bose-Hubbard model features bosonic Peierls transitions (68), which can lead to an interesting interplay of symmetry breaking and symmetry-protected topological order (52, 53).

In terms of experimental implementations, we restricted our discussion in this article to ultracold atom setups. However, other quantum simulation platforms, such as arrays of superconducting qubits (69), provide promising alternatives. Generalizations of our scheme to these systems are straightforward, and a detailed analysis of the feasibility of our proposal in such settings will be devoted to future work.


Implementing dynamical gauge fields

Here, we describe in detail how synthetic gauge fields with their own quantum dynamics can be realized and implemented using ultracold atoms. We begin by quickly reviewing results for the case of a single particle in a double-well potential, which we use later on to derive the effective Hamiltonian in a many-body system.

Single-particle two-site problem. We consider the following Hamiltonian describing a single particle hopping between sites ∣1〉 and ∣2〉H^2=t(21+12)+(Δ2,1+Δ2,1ω(t))22(34)

Here, t > 0 denotes the bare tunnel coupling, which is strongly suppressed by the energy offset ∣Δ2,1 ∣ ≫ t. Tunneling is then restored by a modulationΔ2,1ω(t)=A2,1cos (ωt+ϕ2,1)(35)

For resonant shaking, ω = Δ2,1, it has been shown in (47) that the dynamics of Eq. 34 can be described by the following effective HamiltonianH^2,eff=t(21eiϕ2,1+12eiϕ2,1)(36)

The amplitude of the restored tunneling is given byt=tJ1(A2,1ω)(37)

and the complex phase ϕ2,1 is determined directly from the modulating potential Δ2,1ω(t).

More generally, when the offset Δ2,1 = nω is a positive integer multiple n = 0,1,2,3,4, ... of the driving frequency ω, tunneling can also be restored. As shown by a general formalism in (48), the effective Hamiltonian in this case becomesH^2,eff=tn(21einϕ2,1+12einϕ2,1)(38)

For n = 0, the result is independent of the phase ϕ2,1 of the modulation. The tunneling matrix element is renormalized bytn=tJn(A2,1ω)(39)

The first three Bessel functions, n = 0, 1, 2, are plotted in fig. S1 as a function of x = A2,1/ω.

Last, we consider the case when Δ2,1 = −nω, for a positive integer n = 1, 2, 3, …. In this case, we can rewrite the modulation (Eq. 35) asΔ2,1ω(t)=A2,1cos (ωtϕ2,1)(40)

i.e., effectively ω → −ω and ϕ2,1 → −ϕ2,1. By applying the results from Eqs. 38 and 39 for the system with −ω, we obtainH^2,eff=tn(21ein(πϕ2,1)+12ein(πϕ2,1))(41)

The complex phase of the restored hopping in the effective Hamiltonian changes sign, because −ϕ2,1 appears in Eq. 40. In addition, it contains a π phase shift, which takes into account the sign change of the renormalized tunneling matrix element Jn(A2,1ω)=eiπnJn(A2,1ω) if n is odd.

Multiple driving frequencies. Even more control over the restored tunnel couplings can be gained by using lattice modulations with multiple frequency components. Here, we summarize results for the single-particle two-site problem from above, for the case of driving with frequency components ω and 2ω. To do so, we modify our Hamiltonian in Eq. 35 asH^2=t(21+12)+Δ2,122+Δ2,1ω(t)11(42)where the 2π/ω-periodic driving term takes the following formΔ2,1ω(t)=A2,1(1)cos (ωt+ϕ2,1(1))+A2,1(2)cos (2ωt+ϕ2,1(2))(43)

To calculate the effective Hamiltonian, we rewrite the time-dependent Hamiltonian (Eq. 42) in a moving frame by performing a time-dependent unitary transformation realized by the operator (48)R^(t)=exp (iΔ2,1tP^2)exp {i(A2,1(1)ω)sin (ωt+ϕ2,1(1))P^1}×exp {i(A2,1(2)2ω)sin (ωt+ϕ2,1(2))P^1}(44)where we introduced the projectors P^1=11 and P^2=22. In this moving frame, the time-dependent Hamiltonian in Eq. 42 takes the formH2=t12eiΔ2,1t×ei(A2,1(1)ω)sin (ωt+ϕ2,1(1))ei(A2,1(2)2ω)sin (ωt+ϕ2,1(2))+h.c.(45)

Using the Jacobi-Anger identityeiαsin (ωt+ϕ)=k=Jk(α)eik(ωt+ϕ)(46)and time-averaging the time-dependent Hamiltonian in Eq. 45 over a period T = 2π/ω of the drive, we obtain an effective Hamiltonian of the formH^2,eff=12=tn,ei(n2)ϕ2,1(1)+iϕ2,1(2)+h.c.(47)

While this effective Hamiltonian is similar to Eq. 38, the amplitude renormalization now involves a product of two Bessel functionstn,=tJn2(A2,1(1)ω)J(A2,1(2)2ω)(48)

Two-particle two-site problem. Now, we apply the results from the first paragraph (Eqs. 34 to 41)] to the problem of a pair of a- and f-particles in a double-well potential (see Fig. 2). In contrast to the main text, we consider general parameters in our derivation of the effective Hamiltonian. Our starting point is the model in Eqs. 1 to 6 for two sites j1 and j2 = j1 + ey. We assume Va(j1,2) ≡ 0 but introduce a static energy offset Δf = U between the two lattice sites for the f-particles, Vf(j2) = Δf + Vf(j1). Because our analysis is restricted to the subspace with one a-particle and one f-particle, the hard-core constraint assumed in the main text is not required in this case and the statistics of the two species are irrelevant.

The two-site problem has four basis states, f^jma^jn0 with m, n = 1,2. Their corresponding on-site energies are 0, Δf = U, U, Δf + U = 2U (see fig. S2A), which suppress most coherent tunneling processes because Δf=Utya,tyf. When the resonant lattice modulation H^ω(t) with frequency ω = U is included, all tunnel couplings are restored. Now, we will show that the effective Floquet Hamiltonian is given byH^eff2well=tyaλeiφ^a^j2a^j1tyfΛ^eiθ^τ^j2,j1++h.c.(49)where τ^j2,j1+=f^j2f^j1 andφ^=ϕj2,j1+(1τ^j2,j1z)(π2ϕj2,j1)(50)λ=J1(Aj2,j1/ω)(51)θ^=2ϕj2,j1n^j2a(52)Λ^=J0(Aj2,j1/ω)n^j1a+J2(Aj2,j1/ω)n^j2a(53)

To derive Eqs. 49 to 53, we first consider the effect of the coherent driving H^ω(t), characterized by Eq. 9, on the matter field a^. Because the HamiltonianH^a=tya(a^j2a^j1+h.c.)+(n^j2an^j1a)×12(Uτ^j2,j1z+Vω(j2,t)Vω(j1,t))(54)

governing the dynamics of a^, commutes with the link variable, characterizing the gauge field, [H^a,τ^j2,j1z]=0, we can treat τ^j2,j1z as a ℂ-number with two possible values, ±1.

When expressed in terms of the two states 1a=a^j10 and 2a=a^j20, the Hamiltonian H^a is of the same form as H^2 in Eq. 34. It has an energy difference of Δ2,1=Uτ^j2,j1z between the two states, which is caused microscopically by the interspecies Hubbard interaction U (see fig. S2A).

According to Eqs. 38 and 41, the restored tunnel coupling between ∣1〉a and ∣2〉a has a complex phase given by φ = ϕj2, j1 if Δ2,1> 0, i.e., for τj2,j1z=1, and it is φ = π − ϕj2, j1 if Δ2,1 < 0, i.e., for τj2,j1z=1. Because in both cases the magnitude of the energy mismatch between the two sites is ∣Δ2,1 ∣ = ω, the tunneling is renormalized by λ = 𝒥1(Aj2, j1/ω). These results confirm Eqs. 50 and 51.

Next, we consider the dynamics of the f-particles or, equivalently, the link variable τ^j2,j1z. It is governed by the following HamiltonianH^f=tyf(τ^j2,j1++h.c.)+τ^j2,j1z×12(Δf+Uδn^a+Vω(j2,t)Vω(j1,t))(55)

Because H^f commutes with the matter field, [H^f,n^j1a]=[H^f,n^j2a]=0, we can treat the particle number imbalanceδn^a=n^j2an^j1a(56)as a ℂ-number now, which can take two values ±1.

When expressed in terms of the two states 1f=f^j10 and 2f=f^j20, the Hamiltonian H^f is of the same form as H^2 in Eq. 34. It has an energy difference of Δ2,1=Δf+Uδn^a between the two states, which is caused microscopically by the interspecies Hubbard interaction U and the potential gradient Δf, which the f-particles are subject to (see fig. S2A).

In the case of f-particles, the energy offset Δ2,1 can only take positive values 0 and 2ω if Δf = U = ω. From Eq. 38, it follows that the restored tunnel coupling between ∣1〉f and ∣2〉f has a complex phase given by θ = 0 if Δ2,1 = 0, i.e., for δna = − 1, and by θ = 2ϕj2, j1 if Δ2,1 = 2ω, i.e., for δna = 1. Expressed in terms of n^j2a, in a subspace where n^j1a+n^j2a=1, this result confirms Eq. 52.

The magnitudes of the restored tunneling couplings of f-particles in the two-particle Hilbert space depend on the energy offset Δ2,1. In the case when Δ2,1 = 0, i.e., for δna = −1, it becomes Λtyf=tyfJ0(Aj2,j1/ω). When Δ2,1 = 2ω, i.e., for δna = 1, it is given another Bessel function, Λtyf=tyfJ2(Aj2,j1/ω). This result, summarized in fig. S2B, confirms Eq. 53.

Realizations with ultracold atoms. Next, we discuss realizations of the two-particle two-site problem with ultracold atoms. The proposed implementation needs two distinguishable particles with strong interspecies on-site interaction energy Uty. The particles occupy a double well with both species-dependent and species-independent on-site potentials. For the species-dependent contribution, a static potential is sufficient, which introduces a tilt Δf = U between neighboring sites for the f-particles but leads to zero tilt for the a-particles. On the other hand, the species-independent contribution must be time-dependent Vω(t) to restore resonant tunneling for both particles.

For ultracold atoms, a cubic array of lattice sites with period ds can be created by three mutually orthogonal standing waves with wavelengths λ = 2ds. When extending this simple cubic lattice along one axes by an additional lattice with twice the period dl = 2ds, a superlattice of the form Vssin2(πy/ds+π/2)+Vlsin2(πy/dl+ϕSL/2) arises. In the limit VlVs, the superlattice potential resembles a chain of double wells, where tunneling between each double well is suppressed and all dynamics is restricted to two sites. Tuning the relative phase ϕSL allows dynamic control of the on-site potentials. In principle, the time-dependent modulation Vω(t) can be implemented by a fast modulation of ϕSL; however, the modulation frequency may be limited to small values depending on the implementation of the lattices. For a superlattice with a common retroreflector, for instance, the phase ϕSL can only be varied by changing the frequency of the laser. Alternatively, a second lattice Vmodsin2(πy/dl+ϕmod/2) with period dl and phase ϕmod = ±π/2 can be introduced such that it only affects the on-site potential of a single site of the double well. Therefore, amplitude modulation Vmod(t) of this additional lattice induces a relative modulation of the on-site energies. This leads to a nonzero species-independent, time averaged energy offset, which can be compensated by the static phase degree of freedom ϕSL of the superlattice.

The two distinguishable particles can be encoded in different hyperfine sublevels with different magnetic moments, enabling the direct implementation of the static species-dependent potentials by a magnetic field gradient. This is especially appealing for bosonic atoms having a hyperfine sublevel with zero magnetic moment, which directly results in a vanishing, magnetic field–independent tilt for the a-particles in first order. Nevertheless, this is not essential as tilts for the a-particles can be compensated by the present species-independent potentials.


Supplementary material for this article is available at

Section S1. Flux attachment in 2D

Section S2. Implementing matter coupled to a ℤ2 gauge field in the two-leg ladder geometry

Section S3. Gauge structure of two-leg ladders

Section S4. Phase transitions of gauge and matter fields

Fig. S1. Renormalized tunneling amplitudes determined by Bessel functions.

Fig. S2. Two-site two-particle problem.

Fig. S3. Flux attachment in a 2D Hofstadter model.

Fig. S4. Derivation of the effective Hamiltonian.

Fig. S5. Wilson loop scaling.

Fig. S6. The Green’s function in the charge sector.

Fig. S7. The Luttinger-K parameter.

Fig. S8. Rung-Mott state at commensurate filling.

Fig. S9. Phase diagram of the ℤ2 LGT on a two-leg ladder for commensurate filling.

Fig. S10. Transition in the gauge sector.

References (7074)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank I. Bloch and M. Lohse for fruitful discussions. We also acknowledge discussions with P. Hauke, P. Zoller, V. Kasper, A. Bermudez, L. Santos, I. Carusotto, and M. Hafezi. Funding: The work in Brussels was supported by the FRS-FNRS (Belgium) and the ERC Starting Grant TopoCold. The research in Munich was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via Research Unit FOR 2414 under project number 277974659 and under Germany's Excellence Strategy - EXC2111 - 390814868, the European Commission (UQUAM grant no. 5319278), and the Nanosystems Initiative Munich (NIM) grant no. EXC4. The work at Harvard was supported by the Gordon and Betty Moore Foundation through the EPiQS program, Harvard-MIT CUA, NSF grant no. DMR-1308435, AFOSR MURI Quantum Phases of Matter (grant FA9550-14-1-0035), and AFOSR MURI: Photonic Quantum Matter (award FA95501610323). F.G. also acknowledges support from the Technical University of Munich—Institute for Advanced Study, funded by the German Excellence Initiative and the European Union FP7 under grant agreement 291763, from the DFG grant no. KN 1254/1-1, and DFG TRR80 (Project F8). Author contributions: F.G. and N.G. devised the initial concepts. F.G. performed the main analytical calculations, with inputs from N.G. and E.D. All DMRG simulations were performed by L.B. The proposed experimental implementation was devised by C.S., M.A., N.G., and F.G. All authors contributed substantially to the analysis of the theoretical results. The manuscript was prepared by F.G., N.G., L.B., and C.S., with inputs from all other authors. Competing interests: The authors declare that they have no competing interests. Data materials and availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article