Non-local order in Mott insulators, Duality and Wilson Loops
Steffen Patrick Rath, Wolfgang Simeth, Manuel Endres, Wilhelm Zwerger
NNon-local order in Mott insulators, Duality and Wilson Loops
Steffen Patrick Rath a, ∗ , Wolfgang Simeth a , Manuel Endres b , Wilhelm Zwerger a a Technische Universität München, James-Franck-Straße 1, 85748 Garching, Germany b Max-Planck-Institut für Quantenoptik, 85748 Garching, Germany
Abstract
It is shown that the Mott insulating and superfluid phases of bosons in an optical latticemay be distinguished by a non-local ’parity order parameter’ which is directly accessiblevia single site resolution imaging. In one dimension, the lattice Bose model is dual toa classical interface roughening problem. We use known exact results from the latter toprove that the parity order parameter exhibits long range order in the Mott insulatingphase, consistent with recent experiments by Endres et al. [Science , 200 (2011)]. Intwo spatial dimensions, the parity order parameter can be expressed in terms of an equaltime Wilson loop of a non-trivial U (1) gauge theory in dimensions which exhibitsa transition between a Coulomb and a confining phase. The negative logarithm of theparity order parameter obeys a perimeter law in the Mott insulator and is enhanced bya logarithmic factor in the superfluid. Keywords:
Mott insulators, optical lattices, non-local order, Wilson loops, duality
1. Introduction
A central assumption underlying both classical and quantum physics is that its basiclaws are local , i.e., that the fundamental equations can be expressed in terms of relationsbetween physical observables at a given point in space and time [1]. Symmetries anddifferent types of order are thus related to the behavior of correlation functions of localobservables. In particular, a qualitative change in macroscopic properties is typicallyassociated with the appearance of long range order in some local observable ˆ O ( x ) . Itstwo-point correlation (cid:104) ˆ O ( x ) ˆ O ( y ) (cid:105) → O ( ∞ ) (cid:54) = 0 (1)thus approaches a finite constant as | x − y | goes to infinity. In recent years, a lot ofinterest has focussed on systems where this standard characterization of different phasesof matter by finite values of some local order parameter ˆ O ( x ) fails. This is the case, e.g.,in the Quantum Hall Effect, a paradigmatic example of a topological insulator. It is anincompressible state described by a Chern-Simons theory which has gapless excitationsassociated with chiral edge currents [2, 3]. ∗ Corresponding author, phone +49 89 289 12365, fax +49 89 289 12638
Email addresses: [email protected] (Steffen Patrick Rath), [email protected] (WilhelmZwerger)
Preprint submitted to Elsevier June 4, 2018 a r X i v : . [ c ond - m a t . qu a n t - g a s ] M a y ur aim in the present paper is to study non-local orders that may be used to char-acterize Mott insulating phases of lattice bosons. Quite generally, Mott insulators aredefined by their incompressibility κ = ∂n/∂µ ≡ , a response function that is not as-sociated with long range order in any local observable ˆ O ( x ) . Specifically, we focus onMott insulators which do not break any lattice symmetries due to, for instance, theformation of a commensurate charge density wave [4] or due to Neel order in the spindegrees of freedom, as happens in Mott-Heisenberg insulators [5] realized in undopedhigh-temperature superconductors [6]. To investigate possible non-local orders that mayexist in otherwise featureless Mott insulators, the particular case of bosons in an opticallattice is of special interest. By a simple change of the lattice depth, they may be tunedto undergo a superfluid (SF) to Mott insulator (MI) transition [7]. Moreover, thanksto the direct accessibility of the atomic distribution by optical imaging methods whichallow measuring in situ density distributions [8] and even arbitrary density correlationsat the single-atom level [9, 10], they also provide a new perspective on the microscopicdetails of the involved states. The question of whether some non-local order exists in theMI phase of bosons in an optical lattice has been addressed before by Berg et al. [11]in the particular case of one dimension. While their main focus has been the study ofhidden ’string-order’ that appears in the presence of repulsive interactions of finite range,they have shown via Bosonization that even in the much simpler situation of pure on-siterepulsion, there is a non-local observable which takes finite values in the MI and is zeroin the SF. Using single site imaging, the associated ’parity order parameter’ (POP) hasbeen measured by Endres et al. [12]. For lengths of up to eight lattice spacings theirresults were consistent with the theoretical expectation of a non-vanishing parity orderin the MI and an algebraic decay to zero in the superfluid phase.Our aim in the following is a detailed study of parity order in d = 1 and also in d = 2 dimensions using duality transformations. Specifically, in d = 1 , a number ofexact results for the parity order parameter are obtained from a systematic expansionaround the atomic limit and from analytical results for the roughening transition of theequivalent classical interface model in two dimensions [13, 14, 15]. In the d = 2 case, thelattice bosons are dual to a three-dimensional U (1) lattice gauge theory. This mappinghas first been discussed by Peskin [16] starting from the classical three-dimensional XYmodel and has been extended to the two-dimensional quantum XY model by Fisherand Lee [17]. As will be shown here, the parity order parameter is mapped in the dualmodel on a quantity which is ’more local by one dimension’: in d = 1 , it is mapped ona two-point correlation function for a local operator of the type in Eq. (1) which showsan algebraic decay in the SF phase and converges to a constant in the MI. In the d = 2 case, the POP is mapped onto an equal time Wilson loop in the dual gauge theory,which serves to distinguish the superfluid and Mott insulating phases according to thedependence on the system size L . In particular, we find a perimeter law dependence inthe MI phase which leads to an exponential decay of the POP while in the SF phase alogarithmic correction to the perimeter law leads to a super-exponential decay.
2. Bose-Hubbard Model and Parity Order
Ultracold bosons in an optical lattice provide a realization of the Bose-Hubbard model,as suggested theoretically by Jaksch et al. [18] and first realized experimentally by Greiner2t al. [7]. The associated many-body Hamiltonian [19] ˆ H BH = − J (cid:88) (cid:104) i,j (cid:105) ˆ a † i ˆ a j + U (cid:88) i ˆ n i (ˆ n i − , (2)describes the competition between a kinetic energy which involves hopping to nearestneighbor lattice sites with amplitude J > and an on-site repulsive interaction U > which leads to an increase in energy if atoms hop to sites which are already occupied.Specifically, the operator ˆ a i destroys a boson in a single particle Wannier state localizedat lattice site i and the associated occupation number operator ˆ n i = ˆ a † i ˆ a i has eigenvalues , , , . . . . We consider this model in both one and two spatial dimensions, specializingto the case of a simple quadratic lattice in the latter case. In either dimension, the groundstate of the Bose Hubbard model exhibits a continuous transition between a superfluidand a Mott insulating state which may be tuned either by changing the ratio J/U atfixed density or by changing density at fixed
J/U via the dimensionless chemical potential µ/U [19]. The universality class of this quantum phase transition is different in bothcases. For the density driven Mott transition, the associated dynamical scaling exponentis z = 2 [20]. By contrast, changing the ratio J/U at fixed integer density ¯ n = 1 , , . . . ,the transition at the tip of the Mott lobes has z = 1 , i.e., it is described by an O (2) -modelin D = d + 1 dimensions which possesses a formal relativistic invariance [20]. It is thistype of transition which will be considered in the following.The superfluid phase of the Bose Hubbard model is characterized by a conventional,local observable ˆ O ( x ) = ˆ a i associated with long range order in the one-particle densitymatrix, with O ( ∞ ) ≡ n as the condensate density . This order can be observed ina direct manner in time-of-flight images [7], which measure the momentum distributionas the Fourier transform of the one particle density matrix [21]. The excellent quantita-tive agreement between the measured absorption images after time-of-flight and precisequantum Monte Carlo calculations which include both the effects of finite temperatureand the harmonic trap potential [22], shows that the Bose-Hubbard Hamiltonian (2)provides a faithful description of cold atoms in an optical lattice. In the Mott insulatingphase, the correlation function of the local bosonic field operator vanishes exponentiallylike (cid:104) ˆ a † ( x )ˆ a (0) (cid:105) ∼ exp( −| x | /ξ ) with a correlation length ξ ∼ / ∆ that diverges like theinverse of the Mott gap ∆ . The Mott phase is characterized by its incompressibilityand has no associated order parameter, evolving in a continuous manner from a ther-mally disordered state as the temperature is lowered below the Mott gap ∆ . Note thatthe observation of peaks in the noise correlations (cid:104) ˆ n ( x )ˆ n ( x (cid:48) ) (cid:105) in the Mott phase aftertime–of–flight by Fölling et al. [23] does not reflect any long range order: they are a con-sequence of the fact that the Fourier components (cid:80) R n R exp[ i ( k − k (cid:48) ) · R ] of the average density n R = (cid:104) ˆ n R (cid:105) are of order N at wave vector differences k − k (cid:48) = m ( x − x (cid:48) ) / (cid:126) t which are equal to a reciprocal lattice vector G .As will be discussed in the following, the fundamental difference between the super-fluid and Mott insulating phase in terms of their compressibility implies that there is This would not be true if the bosons were charged as, e.g., the Cooper pairs of a conventionalsuperconductor, where the role of ˆ O ( x ) is expected to be taken by the bi-Fermion operator ˆ ψ ↑ ˆ ψ ↓ ( x ) .The associated correlation function (cid:104) ˆ O † ( x ) ˆ O ( y ) (cid:105) , however, is not gauge invariant and therefore does notconstitute a proper order parameter, as noted by Wen [2]. (cid:104)O ( D ) (cid:105) = (cid:68) e iπ (cid:80) i ∈D (ˆ n i − ¯ n ) (cid:69) = (cid:42)(cid:89) i ∈D ( − ¯ n ˆ p i (cid:43) , (3)where D is a spatial domain, i.e., an interval in d = 1 and an area in d = 2 , and ˆ p i = ( − ˆ n i is the parity operator on lattice site i . There are two important proper-ties of this observable which should be noted right away: First of all, the observable iseasily accessible in experiments since, due to light-induced collision losses, quantum gasmicroscopes directly measure the on-site parity rather than the actual occupation num-bers [10, 9]. As a second point, the observable must be calculated and measured in an open domain D which is part of a larger system, otherwise (cid:104)O ( D ) (cid:105) ≡ would triviallybe equal to one due to conservation of particle number. To study the dependence on thesize, we shall characterize the domain D by its linear extension L measured in units ofthe lattice spacing. In analogy to the standard definition (1) of long range order, theparity order parameter (POP) (cid:104)O ( L ) (cid:105) exhibits long range order if (cid:104)O ( ∞ ) (cid:105) is finite.This kind of order parameter is analogous to those studied in Ising models with alocal gauge invariance, which have no conventional phase transitions to states with longrange order, yet may exhibit different phases distinguished by non-local order parameters[24, 25, 26] (see also section 8). In the context of cold atoms, a more complicated ’stringorder’ parameter was introduced, which characterizes a Haldane insulator that can formin one dimensional systems with longer range interactions [27]. Our focus is on thebehavior of the parity order parameter at the conventional SF–MI transition, not only in1d [11, 12] but also in 2d. In particular, we will use a duality transformation to show thatin two dimensions parity order is related to an equal time Wilson loop in a non-trivial U (1) gauge theory which exhibits different behavior as a function of system size L in theMI and SF phases.
3. Number fluctuations and area law
To obtain a qualitative understanding of the dependence of the non-local order definedin (3) on the size L of the domain we start by giving some qualitative arguments for theexpected scaling behavior deep in the MI, where J/U (cid:28) . Starting with the atomiclimit, it is obvious that (cid:104)O ( L ) (cid:105)| J =0 = 1 in any dimension since particle fluctuationsare then completely frozen. For small but finite J/U and for the simple case of a MIwith average density ¯ n = 1 , particle number fluctuations appear as pairs of empty anddoubly occupied sites. Now, for an arbitrary number of pairs which are completelyinside the domain D , the associated two minus signs in the product in Eq. (3) cancel.Only those pairs which are separated by the domain boundary lead to a reduction of (cid:104)O ( L ) (cid:105) [12]. Intuitively one thus expects the parity order (cid:104)O ( L ) (cid:105) ∼ exp( − L d − ) toscale exponentially with the area of the domain’s boundary. This intuitive expectationis supported by a systematic perturbative calculation in an expansion around the atomiclimit J/U = 0 which—as will be shown in section 4 below—yields (cid:104)O ( L ) (cid:105) = 1 − n (¯ n + 1) dL d − (cid:18) JU (cid:19) + · · · (4)4 igure 1: (Color online) Illustration of how the order parameter is reduced from unity for d = 2 . Thegrid lines indicate the lattice, bosons are represented by small circles. The domain D , here taken tobe a square, is shaded in gray. The minus signs from pairs which are completely inside or outside thedomain (yellow ellipses) cancel out so that there is no contribution while pairs which are separated bythe domain boundary (red ellipse) contribute a minus sign which leads to a reduction of (cid:104)O ( L ) (cid:105) . up to second order in J/U (cid:28) . Clearly the expansion is well defined only in d = 1 ,while in higher dimensions the effective expansion parameter ( J/U ) L d − is small onlyup to system sizes of order ( U/J ) / ( d − .Further insight into the origin of the perimeter law for the decay of the parity ordercan be gained by assuming that the expectation value in (3) may be calculated within aGaussian approximation such that (cid:104)O ( L ) (cid:105) ≈ e (cid:104) ( iπ (cid:80) i ∈D δ ˆ n i ) (cid:105) / = e − π (cid:104) δ ˆ N (cid:105) / , (5)where δ ˆ n i = ˆ n i − ¯ n . Within this approximation, − ln (cid:104)O ( L ) (cid:105) is simply a measure of thetotal number fluctuations (cid:104) δ ˆ N (cid:105) in a domain of size L as part of an infinite system. Now,the standard thermodynamic relation (cid:104) δ ˆ N (cid:105) = k B T ∂N ( µ ) /∂µ in this effectively grandcanonical situation seems to indicate that these fluctuations vanish at zero temperaturewhich would imply a trivial result (cid:104)O ( L ) (cid:105) ≡ in the Gaussian approximation. Thisis not true, however, because the relation only applies in the thermodynamic limit andneglects boundary terms. For a careful calculation of (cid:104) δ ˆ N (cid:105) at zero temperature and in afinite system, we generalize the analysis of Giorgini et al. ([28], see also [29, 30]) for Bosegases in a d = 3 continuum to arbitrary dimensions. The particle number fluctuations (cid:104) δ ˆ N (cid:105) = S d (cid:90) L d r r d − τ ( r )¯ nν ( r ) , (6)in a spherical domain of radius L can then be calculated from the pair distributionfunction ν ( r ) = δ ( r ) + n (cid:16) g (2) ( r ) − (cid:17) = (cid:90) d d q (2 π ) e i q · r S ( q ) (7)5nd the volume τ ( r ) of the intersection between two d -dimensional balls of radius L separated by the distance r . Here, S d is the surface of a unit sphere in d dimensionswhile S ( q ) is the standard static structure factor.In the superfluid, the zero temperature static structure factor has the non-analyticbehavior S ( q ) (cid:39) α | q | at small wave numbers, characteristic for any compressible phase.Since collective excitations exhaust the f sum rule for long wavelengths, the prefactor α = (cid:126) / mc s is, moreover, completely fixed by the exact sound velocity c s . As a result ofthe non-analytic behavior of S ( q ) , the associated pair distribution function exhibits analgebraic decay ν ( r ) (cid:39) − α/πr in d = 1 and ν ( r ) (cid:39) − α/ πr in d = 2 at long distances.The behavior of the number fluctuations for large L is then found to be (cid:104) δN (cid:105) ∼ α ¯ nL d − ln( L/ξ h ) ( SF ) , (8)where the effective healing length ξ h = (cid:126) /mc s serves as a short-distance cutoff. Withinthe Gaussian approximation, therefore, the parity order (cid:104)O ( L ) (cid:105) ∼ (cid:40) L − πα ¯ n ( d = 1) L − π α ¯ nL ( d = 2) ( SF ) (9)decays to zero with a power law in the superfluid phase in one dimension while in twodimensions, the decay is super-exponential.In the MI phase, the quadratic number fluctuations are much smaller and the scalingof (cid:104) δN (cid:105) with system size L differs in a qualitative manner from that in the SF. Indeed, atzero temperature, incompressibility of the MI implies that the structure factor vanishesanalytically for q → and thus S ( q ) = γ ( qξ ) + · · · to leading order in an expansionaround q = 0 . Here, ξ is the characteristic length of the exponential decay of theone-particle density matrix mentioned in the previous section, while γ is a numericalconstant which relates the scales appearing in the density correlation and in the one-particle density matrix. From perturbation theory, one finds that in one dimension γ isproportional to ( J/U ) to leading order [31].Since the static structure factor is analytic around q = 0 , the pair correlation function ν ( r ) decays exponentially on the characteristic scale ξ , effectively cutting off the integralin Eq. (6). As a result, one finds that (cid:104) δN (cid:105) → bL d − ( MI ) (10)scales like the area of the boundary of the domain with a coefficient b = b γ ¯ nξ , where b is a numerical constant depending on the geometry and the precise functional form of ν ( r/ξ ) .The results obtained within this Gaussian approximation are—of course—not exact.Yet, as will be shown in sections 6 and 7 below, they give the correct qualitative behaviorof the parity order parameter as a function of the size L of the domain. On a basiclevel, therefore, the different behavior of the non-local parity order in the MI and theSF is a simple consequence of the fundamental difference between number fluctuationsin an incompressible versus a compressible system [32, 33]. In this context, it is alsoinstructive to note that the underlying scaling (cid:104) δN (cid:105) ∼ α ¯ nL d − ln( L/ξ ) of the numberfluctuations in the compressible and gapless superfluid compared to (cid:104) δN (cid:105) ∼ L d − in theincompressible and gapped MI are reminiscent of similar results obtained for the scaling6f the entanglement entropy S ( L ) . The fact that the fluctations of conserved quantitiesare closely related to the latter, exhibiting a simple area law for gapped phases, hasbeen noted by Swingle and Senthil [34]. There are, however, a number of importantdifferences: the entanglement entropy typically scales like S ( L ) = b L d − with a non-universal prefactor b even in gapless phases while the corresponding number fluctuations (cid:104) δN (cid:105) have an additional logarithmic enhancement factor for any phase with a finitecompressibility. For the entanglement entropy, violations of the area law by a logarithmicfactor of the form S ( L ) ∼ L d − ln( L ) appear, e.g., in free fermions with a Fermi surfaceand also for Landau Fermi liquids, not in a gapless superfluid, however. Indeed, as notedby Metlitski and Grover [35], its entanglement entropy S ( L ) = b L d − + ∆ S exhibits an additive —not multiplicative—logarithmic contribution ∆ S = ln (cid:0) ρ s L d − /c s (cid:1) / which isuniversal, i.e., it only depends on the superfluid stiffness ρ s and the sound velocity c s aseffective low energy constants.The simplicity of the Gaussian approximation for studying the parity order also allowsfor a straightforward discussion of how the above results are affected by a finite temper-ature. For a neutral SF, the structure factor at temperatures k B T (cid:28) mc s reads [36] S ( q ) (cid:39) α | q | coth( c s | q | / k B T ) . The characteristic length scale at which the zero temper-ature result S ( q ) (cid:39) α | q | is equal to the q = 0 thermal value S ( q = 0) = 2 k B T α/c s istherefore given by r T = (cid:126) c s / πk B T = λ T /ξ h , where λ T is the thermal wavelength and ξ h = (cid:126) /mc s the effective healing length. Since the parity order effectively probes numberfluctuations at a finite wave vector q (cid:39) π/L the zero temperature results remain validas long as r T (cid:29) L . Since typical values of the healing length are of order ξ h ∼ µ m, thecharacteristic scale of r T becomes larger than a lattice spacing at temperatures below T ≈ nK.An equivalent reasoning can be applied in the Mott insulating phase, where thethermal behavior of the structure factor S ( q = 0) = ¯ nk B T κ T involves the compressiblity κ T . For dimensional reasons, the latter must be of the form κ T = (¯ n ∆) − f ( β ∆) withthe Mott gap ∆ and a function f ( β ∆) ∼ e − β ∆ which has a thermally activated form.The zero temperature results are then valid as long as the domain size L satisfies L (cid:46) λ T ξ (cid:115) m ∆ f ( β ∆) . (11)In the interesting regime T (cid:28) ∆ this is easily satisfied due to the exponential form of f ( β ∆) . The actual constraint is rather the condition β ∆ (cid:29) which becomes increasinglyhard to satisfy as one approaches the critical point.
4. Perturbative analysis in the MI
Deep in the MI, where
J/U (cid:28) , exact results for the parity order can be obtainedby a systematic perturbation theory around the atomic limit J = 0 of the Bose-Hubbardmodel. Here we make use of a method developed by van Dongen [37], which is anextension of a formalism due to Harris and Lange [38]. The basic idea is to constructa canonical transformation of the creation and annihilation operators, ˆ a † = e ˆ S ˆ b † e − ˆ S with an anti-hermitian generator ˆ S such that in the new basis, the occupation numbersand hence the interaction part remain invariant under hopping. Explicitly, one writes7he Hamiltonian in the form ˆ H = U ˆ D + ˆ K and defines the hopping term for the newparticles ˆ T = − J (cid:88) (cid:104) i,j (cid:105) ˆ b † i ˆ b j ⇔ ˆ K = e ˆ S ˆ T e − ˆ S . (12)The requirement that ˆ D is invariant under hopping is obeyed provided it is a constantof motion, i.e., [ ˆ H, ˆ D ] = 0 . This condition fixes the transformation ˆ S , which may becalculated order by order in a systematic expansion ˆ S = (cid:88) n ≥ U n ˆ S n , (13)in powers of /U . By expanding the exponential and substituting the expansion of theoperator ˆ S , the POP takes the form (cid:104)O ( L ) (cid:105) = 1 + e − iπN (cid:26) (cid:88) k ( iπ ) k k ! (cid:104) Φ | (cid:20) U ˆ S + 1 U ˆ S + · · · , (cid:16) ˆ N ( L ) (cid:17) k (cid:21) | Φ (cid:105) + (cid:88) k ( iπ ) k k ! (cid:104) Φ | (cid:20) U ˆ S , (cid:20) U ˆ S , (cid:16) ˆ N ( L ) (cid:17) k (cid:21)(cid:21) | Φ (cid:105) + · · · (cid:27) , (14)where ˆ N ( L ) is the total number operator within the domain D and N is the eigenvalueof ˆ N ( L ) in the atomic limit ground state | Φ (cid:105) . Since | Φ (cid:105) is an eigenstate of ˆ N ( L ) , thefirst term in the curly bracket vanishes at all orders. As a consequence, to calculate thePOP at order n in J/U , one needs to determine the operators ˆ S , . . . , ˆ S n − .Specifically, substituting ˆ S from equation (A.1) (cf. Appendix A) and evaluatingthe commutators, one finds the result (4) stated already in section 3. In one dimension,we have continued the perturbative expansion up to fourth order in J/U . For reasons ofspace, the details of this calculation are deferred to Appendix A and we only give theresult here: (cid:104)O ( L ) (cid:105) = 1 − n (¯ n + 1) (cid:18) JU (cid:19) −
49 ¯ n (¯ n + 1)[¯ n (473¯ n + 217) − (cid:18) JU (cid:19) + · · · . (15)Just as the leading order term, the next-to-leading correction is independent of thedomain size L as long as the latter is larger than the order of perturbation. Hence, (cid:104)O ( L ) (cid:105) is independent of L for sufficiently small J/U , in agreement with a DMRGcalculation in [12] where (cid:104)O ( L ) (cid:105) was found to be essentially independent of L for J/U (cid:46) . for domain sizes ranging from 1 to 60.
5. Duality transformation
In the following, we show that within a reduced description of the SF to MI transitionat fixed density in terms of a quantum rotor model [20], exact results for the parity orderparameter can be obtained from analyzing a ( d + 1) -dimensional classical lattice model.This is particularly interesting in the d = 2 case where the parity order (cid:104)O ( L ) (cid:105) can beexpressed in terms of an equal time Wilson-loop in a nontrivial U (1) gauge field which8s dual to the original lattice boson model. Our aim is twofold: first of all, unlike thequalitative arguments presented in section 3, the duality transformation permits us toderive exact results for the large L behavior of (cid:104)O ( L ) (cid:105) in both the MI and SF phaseswithin one unified framework. From a different point of view, however, the duality of thelattice boson model to a dynamical gauge field in dimensions can also be read theother way, where ultracold atoms in an optical lattice provide a quantum simulator of anontrivial lattice gauge theory.The starting point for the mapping is based on the realization that the SF to MItransition at fixed density is driven by phase fluctuations only [19]. To isolate these fromthe full Bose-Hubbard Hamiltonian, it is convenient to consider the limit of large filling ¯ n (cid:29) , where the boson operators may be rewritten in a density-phase representation, ˆ a j (cid:39) e i ˆ φ j (cid:112) ¯ n + δ ˆ n j [19]. In the limit ¯ n (cid:29) , the number fluctuations δ ˆ n j can be expandedup to second order and the Bose-Hubbard Hamiltonian is thus transformed into theHamiltonian H J = U (cid:88) x ˆ n x + E J (cid:88) x , u [1 − cos( ˆ φ x + u − ˆ φ x )] (16)of a system of quantum rotors at each lattice site x with a discrete eigenvalue spectrum , ± , ± , . . . which are coupled to nearest neighbors x + u by a Josephson energy E J =2¯ nJ . (Note that we have redefined δ ˆ n x → ˆ n x , i.e., ˆ n x is now the deviation from theaverage local occupation number).Following a standard procedure, the partition function associated with the quantumrotor model (16) can now be expressed in terms of a path integral by dividing the interval [0 , β ] into N τ subintervals of width ε = β/N τ . Inserting a complete set of number statesat each step, one thus obtains a discretized path integral representation of the form Z J = (cid:88) n x ,j ∈ Z (cid:104){ n x , }| e − ε ˆ H J |{ n x , }(cid:105) · · · (cid:104){ n x ,N τ }| e − ε ˆ H J |{ n x , }(cid:105) , (17)where the notion { n x ,j } is a reminder of the fact that—at given j = 1 , . . . , N τ —thereare N d variables n x ,j ∈ Z for x ∈ (1 , . . . , N ) d that have to be summed over all integers Z . In the limit ε → , the non-commuting terms in exp( − ε ˆ H J ) can be factorized to give Z J = (cid:88) n x ,j ∈ Z (cid:89) x ,j, u (cid:16) e − εUn x ,j / (cid:104){ n x ,j }| e − εE J [1 − cos( ˆ φ x + u − ˆ φ x )] |{ n x ,j +1 }(cid:105) (cid:17) . (18)The matrix elements of the Josephson coupling terms can now be simplified by using theso called Villain approximation [39] exp {− εE J [1 − cos( ˆ φ x + u − ˆ φ x )] } (cid:39) (cid:88) m x , u exp (cid:32) − m x , u εE J − im x , u ( ˆ φ x + u − ˆ φ x ) (cid:33) . (19)Since the fundamental symmetry φ → φ +2 π due to the discreteness of the boson numberis retained, this leaves the physics qualitatively unchanged. In one spatial dimension, this a constant prefactor in (19) is suppressed because it only gives an irrelevant overall shift of the freeenergy. m lj per lattice site, in two spatial dimensions, itintroduces two integers m x ,j = ( m x ,x,j , m x ,y,j ) at each lattice site.For each given j , one now uses the fact that (cid:89) x exp (cid:16) − im x , u ,j ( ˆ φ x + u − ˆ φ u ) (cid:17) = (cid:89) x exp (cid:16) i ( m x , u ,j − m x − u , u ,j ) ˆ φ x (cid:17) (20)for a periodic chain with ˆ φ x + N u = ˆ φ x and the identity (cid:104) n (cid:48) | e im ˆ φ | n (cid:105) = δ n (cid:48) ,n + m since theoperator e im ˆ φ shifts the particle number by m . The resulting partition function Z J = (cid:88) n x ,j ∈ Z m x , u ,j ∈ Z exp − εU (cid:88) x ,j n x ,j − εE J (cid:88) x , u ,j m x , u ,j (cid:89) x , u ,j δ ∇ x · m x ,j , −∇ τ n x ,j (21)has a Gaussian form, however the variables n x ,j and m x ,j are integer-valued and areconnected by the constraint ∇ x · m + ∇ τ n = 0 , where ∇ x ,τ denotes the discrete derivativeon the dual lattice of links along the physical directions x and the ’time’ direction τ .Thus, the variables n x ,j and m x ,j together form a divergenceless ( d + 1) -dimensionalinteger vector field n ≡ ( n, m ) . In d = 1 , this constraint may be resolved by introducinga single integer field h x,τ such that n = ∇ x h and m = −∇ τ h . The partition functionthen becomes that of the discrete Gaussian model with height variable h which describesthe roughening transition of a 2d interface. Moreover, the POP is mapped on a two-pointcorrelation function of the local variable O ( x ) = exp( iπh ( x )) . This model is discussed insection 6.In d = 2 , the constraint is automatically satisfied if one introduces a three-componentvector potential a such that n = ∇ ∧ a , where the lattice curl is defined as ( ∇ ∧ a x ) i = (cid:80) j,k ε ijk ( a x − ˆ k ,k − a x − ˆ − ˆ k ,k ) . The partition function then becomes that of a (2 + 1) -dimensional U (1) gauge theory, as will be discussed in detail in section 7. Remarkably,under the duality transformation the parity order parameter is mapped onto an equaltime Wilson loop [40] (cid:104)O ( L ) (cid:105) = (cid:42) exp (cid:34) iπ (cid:88) x ∈D d x ( ∇ ∧ a ) τ (cid:35)(cid:43) = (cid:42) exp (cid:34) iπ (cid:88) x ∈ ∂ D (∆ x ) · a (cid:35)(cid:43) , (22)where we have used the discrete version of Stokes’s theorem on the last line, ∂ D is theboundary of the domain D and ∆ x is a unit vector directed along the boundary in thepositive mathematical sense. The Wilson loop is a gauge-invariant quantity which isoften used to characterize phases in gauge theories [41]. As will be shown in section 7,in the present case, the transition between a phase with massless and one with massive’photons’ in the underlying U (1) -gauge theory predicts qualitatively different behavior ofthe parity order in the SF and MI phases of the original lattice Boson model, consistentwith the qualitative considerations discussed in section 3.
6. Discrete Gaussian interface model
As shown in the previous section, in the d = 1 case the partition function can be repre-sented in terms of a single integer h on each lattice site such that ( n, m ) = ( ∇ x h, −∇ τ h ) .10hoosing ε = 1 /U , the resulting partition function Z DG = (cid:88) { h lj } exp − (cid:88) l,j (cid:26) ( ∇ x h lj ) + UE J ( ∇ τ h lj ) (cid:27) (23)defines an anisotropic discrete Gaussian (DG) model for an integer valued height variable h lj above a two-dimensional, perfectly flat interface h lj ≡ (note that an overall shift h lj → h lj + Z of this reference plane is irrelevant). This mapping has been used earlier inthe context of the SF-MI transition in one dimension by one of the present authors [42].The DG model is a classical model which exhibits a phase transition from a smoothinterface in the regime where U dominates to a rough phase in the limit E J (cid:29) U . Thesmooth phase, which corresponds to the MI in the original quantum rotor model, ischaracterized by a finite dimensionless step free energy f s [43] which is related to theMott gap ∆ µ of the dual model (16) by f s = ∆ µ/U [42]. The dimensionless step freeenergy is a decreasing function of E J /U and reaches f s ≡ / at E J = 0 . In this limitan additional boson is described by a step of unit height which is parallel to the τ axis,i.e., the boson world lines exhibit no quantum fluctuations.When E J /U ∼ , the model is essentially isotropic. To render this manifest, it isconvenient to choose ε = 1 / √ E J U so that Z DG = (cid:88) { h lj } exp − (cid:114) UE J (cid:88) l,j (cid:8) ( ∇ x h lj ) + ( ∇ τ h lj ) (cid:9) . (24)In this representation, the ratio ˜ T DG = (cid:112) E J /U of kinetic and interaction energy inthe underlying quantum rotor Hamiltonian (16) plays the role of an effective tempera-ture. The discrete Gaussian model (24) is known to have a roughening transition of theKosterlitz-Thouless type [44] at a critical temperature ˜ T R (cid:39) . . In the smooth phase,the mean square surface displacement ∆ h ( L ) = (cid:104) ( h lj − h l (cid:48) j (cid:48) ) (cid:105) (25)remains finite as the distance L = | ( l, j ) − ( l (cid:48) , j (cid:48) ) | between two points on the surfaceapproaches infinity. By contrast, the rough phase of the discrete Gaussian model at ˜ T DG > ˜ T R is characterized by a logarithmically divergent ∆ h ( L ) ∼ ln( L ) .These results on the discrete Gaussian model, for which the existence of a Kosterlitz-Thouless transition has been proven rigorously by Fröhlich and Spencer [13], can nowbe translated back to understand the nature of non-local order in the original Bose-Hubbard or quantum rotor model. In particular, the qualitative results that were derivedin section 3 within a Gaussian approximation from considering the number fluctuationswithin a domain of linear size L can now be put on a rigorous footing. This relies on thefact that number fluctuations in the original model of bosons hopping on a lattice aretransformed, via the duality, to fluctuations of the normal vector of the 2d interface by Note that the Trotter decomposition in Eq. (17) is usually performed for finite β and thus N τ = β/ε → ∞ requires ε → . Here we keep ε finite but consider the limit of zero temperature β → ∞ . = ∇ x h . As a result, the nonlocal order parameter defined in Eq. (3) translates intothe characteristic function (cid:104)O ( L ) (cid:105) = (cid:104) exp( iπ [ h ( L ) − h (0)]) (cid:105) (26)of the probability distribution p ( h, | ( l, j ) − ( l (cid:48) , j (cid:48) ) | ) = (cid:104) δ ( h lj − h l (cid:48) j (cid:48) − h ) (cid:105) that the heightvariables at two sites at a distance L differ by h .For a detailed understanding of the behavior of the parity order parameter, we cannow use exact results on the classical roughening transition of two-dimensional interfaces.Specifically, within the smooth phase, Forrester [14] has obtained the exact probabilitydistribution for the height of a lattice site near the center in the body-centered solid-onsolid (BCSOS) model for fixed boundary conditions and the thermodynamic limit: itturns out to be a discrete Gaussian distribution p ( h ) = N e − ( h − / / σ , where N is anormalization constant, and with σ = C/ (cid:112) − T /T R where C = (cid:112) / ln 2 . The non-zero expectation of this distribution stems from the fact that in the BCSOS model, onehas to deal with two sublattices where the sites of one take only even values while thesites of the other take only odd values. Since the outermost sites are fixed at values and (according to the sublattice), the height in the center may equivalently be seen as theheight difference to the border, or, in the thermodynamic limit, as the infinite-distancelimit of this difference. Hence, the characteristic function of this probability distributionis just the two-point correlation function discussed above within the DG model. Sincethe DG and the BCSOS model are in the same universality class, we may conclude thatin the limit of infinite distance, the probability distribution for the height differencebetween two points in the DG model equally becomes a discrete Gaussian distribution p (∆ h = n ) = N e − n / σ with σ = C DG / (cid:113) − ˜ T DG / ˜ T R . Using the Poisson formula,one then obtains the characteristic function (cid:88) n ∈ Z p ( n ) e ikn = √ πσ (cid:80) n ∈ Z e − n / σ (cid:88) m ∈ Z e − σ ( k − πm ) / , (27)i.e., the characteristic function is a periodic sum of Gaussians centered on , ± π, ± π, . . . .The prefactor is just the ratio of the normalization constants of the continuous and thediscrete Gaussian distribution and tends to unity for σ (cid:29) . As ˜ T DG approaches ˜ T R from below, σ diverges and the individual peaks of the characteristic function becomevery narrow. There are then only two contributions to the POP (from the peaks centeredon and π ) which goes to zero as (cid:104)O ( ∞ ) (cid:105) = p ( k = π ) ∼ − π C DG (cid:113) − ˜ T DG / ˜ T R , (28)confirming the conjecture (cid:104)O ( L ) (cid:105) ∼ exp {− A [( J/U ) c − ( J/U )] − / } made in [12] on howthe POP reaches zero as one approaches the critical point from the smooth phase (close Note that a priori, C DG (cid:54) = C [14]. The numerical values of the roughening temperatures are model-specific as well.
12o the critical point, writing (cid:104)O ( L ) (cid:105) as a function of J/U rather than ˜ T DG ∝ (cid:112) J/U only affects the non-universal constant A ).To understand the behavior of the parity order within the superfluid, which corre-sponds to the rough phase of the associated 2d classical interface model, it is convenientto use the equivalence between the DG model and the classical two-dimensional Coulombgas (CG) derived by Chui and Weeks [44]. Within the 2d CG picture, the underlying SF-MI transition is translated into a phase transition of the Kosterlitz-Thouless type betweenan insulating phase where charges are bound—the SF phase of the original model—anda metallic phase of effectively free charges which describes the MI. The metallic phaseis characterized by a divergent polarizability (cid:82) d r r p ( r ) , where p ( r ) is the probabilitydistribution function for the distance between a pair of opposite unit charges added tothe system. Indeed, as shown by Chui and Weeks [44], for < ξ < π the two-pointcorrelation function of the DG model maps on the partition function of the Coulomb gasin the presence of two opposite charges ± ξ at a distance r = | r − r | , (cid:68) e iξ [ h ( r ) − h ( r )] (cid:69) = e − β CG F ( r,ξ ) = p ( r ) , (29)where F ( r, ξ ) is the free energy of the neutral Coulomb gas in the presence of the addedcharges. It is well known [45] that when these added charges have unit strength [withinthe mapping by Chui and Weeks, a unit charge corresponds to ξ = 2 π in F ( r, ξ ) ], e − β CG F ( r, π ) = exp (cid:18) − e ε ln r (cid:19) = r − e /ε , (30)with the dielectric constant ε which characterizes the insulating phase of the Coulombgas. At the unbinding transition, the polarizability diverges, i.e., the critical value of thedielectric constant obeys ( e /ε ) c = 4 . In our case of ξ = π , we are dealing with a pairof half-unit charges . As a result the POP (cid:104)O ( L ) (cid:105) ∼ r − e / ε , (31)decays algebraically in the SF phase of the original lattice Bose model, and the mappingto the Coulomb gas gives us the exact value of the universal jump of the exponent at theroughening temperature: for ˜ T DG → ˜ T R from above, the exponent e / ε approaches onebefore jumping to zero, a consequence of the universal jump of the superfluid stiffnessat the SF-MI transition of the 1d quantum rotor model (for a more detailed discussionsee [46]). U (1) gauge theory in dimensions In d = 2 , the duality mapping discussed in section 5 and the introduction of the three-component vector potential a with n = ∇ ∧ a lead to a (2 + 1) -dimensional U (1) gaugetheory [2, 16, 17, 47, 48] (for an alternative approach to this type of duality mapping To avoid possible confusion, we stress that the right-hand side of Eq. (29) is well-defined for arbitrary ξ , but the equality to the left-hand side (i.e., the existence of the mapping) is restricted to values <ξ < π , and thus to fractional charges. For ξ = 2 πn with integer n , one trivially has e iξ [ h ( r ) − h ( r )] = 1 . ε = 1 / √ U E J , which results in a partition function of the form Z = (cid:88) { a } exp (cid:32) − (cid:114) UE J (cid:88) x ,τ ( ∇ ∧ a ) (cid:33) . (32)Since ( ∇ ∧ a ) = e + b is the energy density associated with a two-component ’electricfield’ e i = ∇ τ a i − ∇ i a τ and a scalar ’magnetic field’ b = ∇ x a y − ∇ y a x , this partitionfunction appears to describe the free field theory of pure electrodynamics in (2 + 1) di-mensions [2]. It is gauge invariant since the action only depends on the curl of a and isthus unchanged if the lattice gradient of an arbitrary function of x and τ is added to a .What makes the model in Eq. (32) nontrivial is the fact that all fields take only integervalues on a discrete space-time lattice. The SF to MI transition of the underlying latticeBose model shows up at the level of the dual U (1) gauge theory as a transition betweena phase in which a really is a free field and the photon is massless and one with massivephotons.In order to understand the physical meaning of this transition in terms of the gaugefield degrees of freedom, it is convenient to consider the equal-time one-body densitymatrix in the quantum rotor model (16), which is mapped onto a ratio of two gauge fieldpartition functions (cid:104) ˆ a † ( x )ˆ a (cid:105) ¯ n = (cid:104) e i ˆ φ ( x ) e − i ˆ φ (0) (cid:105) = Z [ x , Z = exp( − ∆ F [ x , , (33)where Z [ x , differs from Z in that the constraint ∇ · n ( y ) = 0 is replaced by ∇ · n ( y ) = δ y , x − δ y , [50]. Physically, this corresponds to a configuration with a pair of oppositelycharged magnetic monopoles situated at x and . The fact that the one-body densitymatrix approaches a constant in the SF and decays exponentially in the MI thus leadsto a fundamentally different behavior of the dimensionless free energy increase ∆ F [ x ,
0] = (cid:40) | x | /ξ ( MI ) const. − c s πρ s | x | ( SF ) (34)associated with the introduction of a monopole–antimonopole pair at distance x in thedual gauge theory (here, ρ s and c s denote the superfluid stiffness and sound velocity,respectively). In particular, the exponential decay of the one-body density matrix in theMI phase translates to a linear confinement while in the SF phase the monopoles interactvia a 3d Coulomb potential.An effective low energy description of the gauge theory which properly accounts forthe two different phases is provided by a Gaussian model of the form S = 12 ˜ T (cid:90) d x (cid:26) [ ∇ ∧ a ( x )] + 1 ξ [ a ( x )] (cid:27) , (35)where a is now treated as a continuous variable and ˜ T (cid:39) (cid:112) E J /U is a renormalizeddimensionless temperature or coupling constant, similar to the one in the previous section.In the SF regime, where ˜ T is above a critical value of order one, the gauge field is in14ts Coulomb phase where ξ = ∞ . Elementary excitations are then massless photons,which are just the phonons of the superfluid with linear dispersion ω = c s | q | (notethat ’photons’ in (2 + 1) -dimensional electrodynamics have no polarization degrees offreedom). By contrast, in the MI for small values of ˜ T , the gauge field is in a confiningphase, with ξ = 1 /m finite. The photons thus acquire a mass m and now represent theelementary particle–hole excitations of the MI with dispersion ω = c s (cid:112) m c s + q (atthe transition, one has c s = 4 . J for the Bose–Hubbard model with ¯ n = 1 [51]). Whenevaluating expectations for the massless phase, it is convenient to keep a finite ξ duringthe calculation and only take the limit ξ → ∞ at the end of the calculation, which avoidsthe necessity of introducing an explicit gauge-fixing term [52]. This can be seen explicitlyin the correlation function of the vector potential, which reads (cid:104) a ( q ) j a ( q (cid:48) ) k (cid:105) = ˜ T (2 π ) δ ( q + q (cid:48) ) (cid:20) [ P t ( q )] jk q + ξ − + [ P l ( q )] jk ξ − (cid:21) , (36)where [ P t ( q )] jk = δ jk − q j q k / q and [ P l ( q )] jk = q j q k / q are the components of thetransverse and longitudinal projector with respect to q , respectively. Since the effectivemodel (35) is Gaussian, the calculation of the parity order parameter from the Wilsonloop in Eq. (22) is now easy and gives (cid:104)O ( L ) (cid:105) = exp (cid:32) − π (cid:42)(cid:18)(cid:90) D d x [ ∇ ∧ a ( x )] τ (cid:19) (cid:43)(cid:33) . (37)Using (36), the expectation appearing in the exponent in Eq. (37) reduces to (cid:42)(cid:18)(cid:90) D d x [ ∇ ∧ a ( x )] τ (cid:19) (cid:43) = ˜ T (cid:90) d q (2 π ) q x + q y q + ξ − (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) D d x e i q ⊥ · x (cid:12)(cid:12)(cid:12)(cid:12) , (38)where q ⊥ = ( q x , q y ) . Note that q has three finite components whereas x is restricted tothe τ = 0 plane. For a circular disk of radius R , the spatial integral appearing in (38)gives (cid:90) D d x e i q ⊥ · x = 2 πRq ⊥ J ( q ⊥ R ) . (39)In the deconfined phase, which corresponds to the superfluid, the mass parameter /ξ isequal to zero and the resulting parity order parameter is given by − ln (cid:104)O ( R ) (cid:105) ∼ π ˜ T R ln( πR ) . (40)in agreement with the qualitative result obtained in Eq. (9) within the Gaussian approx-imation. Conversely, in the confined phase corresponding to the MI, The q is negligiblecompared with /ξ , leading to a perimeter law − ln (cid:104)O ( R ) (cid:105) ∼ ˜ T ξ πR ≡ R(cid:96) O . (41) The resulting scaling of the parity order will thus turn out to be the same as found in section 3 upto numerical factors, since both are obtained from a Gaussian model. Note, however, that the mappingof (cid:104)O ( L ) (cid:105) to a Wilson loop in the dual gauge theory provides a proper justification for these scalings. The momentum integrals in (38) are cut off at Λ ∼ π since the original problem is defined on alattice with lattice constant set equal to unity. (cid:96) O ∼ ∆ (cid:112) U/J close to thecritical point. Remarkably, a measurement of the parity order, which only involves thestatistics of number fluctuations, therefore allows to extract the Mott gap.
8. Conclusion & Outlook
We have presented a detailed study of parity order (cid:104)O ( L ) (cid:105) for lattice bosons inboth one and two spatial dimensions, using duality transformations. Consistent withprevious theoretical work [11] and recent experiments [12], it has been shown that theMott insulating phase in one dimension exhibits long range parity order. An intuitiveunderstanding of this result relies on the observation that for any incompressible phase,the number fluctuations (cid:104) δ ˆ N (cid:105) ∼ L d − in a domain of size L scale with the area of theboundary ∼ L d − . Using the duality to a discrete, classical interface roughening problemin two dimensions, these results have been put on a rigorous footing. In two dimensions,the parity order again shows qualitatively distinct behavior in the MI and SF phase.In this case, the variable (cid:104)O ( L ) (cid:105) can be expressed in terms of an equal time Wilsonloop of a nontrivial U (1) gauge theory in dimensions. This is related to the factthat the density fluctuations in the original lattice Bose model are mapped to the scalarmagnetic field in the dual gauge theory. A quite interesting result obtained from thismapping is the fact that the decay of parity order in the MI allows to measure the Mottgap from the statistics of number fluctuations. Since experimental measurements of theparity order in two dimensions are straightforward in principle, this seems a promisingroute to infer microscopic information from single site resolution imaging which is verydifficult to obtain otherwise. The detailed numerical factors in the scaling of (cid:104)O ( L ) (cid:105) can unfortunately not be predicted from our effective long wavelength description of thequantum rotor model. Since we are dealing with a bosonic system, however, effectivenumerical methods are available to obtain precise results in a realistic experimentalsetup, as has been shown for thermodynamic properties and excitation energies [51].In particular, numerical simulations directly deal with the Bose Hubbard model whichapplies in the relevant case of low filling ¯ n = 1 instead of the qualitatively similar situation ¯ n (cid:29) that has been studied here within the quantum rotor model.An important open question is, of course, to which extent our results for Bose Mottinsulators can be generalized to the fermionic case, which has also been realized experi-mentally with ultracold atoms [53, 54]. Based on the qualitative description in terms ofthe scaling of number fluctuations in section 3, we expect that the results obtained herecarry over to the fermionic case despite the fact that no duality transformations exist inthis case which allow to connect the parity order with a Wilson loop in a dual gauge the-ory. The presence of long range parity order for 1d Mott insulators also in the fermioniccase is consistent with the results of a recent DMRG study of the fermionic Hubbardmodel by Montorsi and Roncaglia [55]. A different kind of non-local order characterizedby sub-lattice parity was found in this Model by Kruis et al. several years earlier [56, 57].Finally, a quite interesting direction of further research on non-local orders for coldatoms in optical lattices is connected with the recent realization of a 1d transverse Isingmodel using a tilted optical lattice [58]. Extending this setup to the case of two dimen-sions, a number of complex phases may appear depending on the type of lattice and thedirection of the tilt [59]. A quite intriguing perspective would appear in a setup that16llows to realize a ferromagnetic version of this standard model for quantum phase transi-tions [20]. In the ferromagnetic case, the transverse Ising model is self-dual in one dimen-sion, while in two dimensions the dual theory is given by an Ising gauge theory [41]. Forthe latter, one can define non-local correlation functions as C ( D ) = (cid:104) (cid:81) m ∈D ˆ S xm (cid:105) , where ˆ S xm is the x component of the spin operator at site m and D is an area in two dimensions,in complete analogy with our definition of the parity order in equation (3). After theduality transformation, this observable is transformed into C ( D ) = (cid:104) (cid:81) k ∈ ∂ D ˆ σ zk (cid:105) , where ˆ σ zk is the z component of the spin operator at site k in the dual theory and ∂ D arethe two sites at both ends of the string in one dimension or the border of the area intwo dimensions. As a result, one has in one dimension that lim L →∞ C ( L ) > for theparamagnetic phase. In two dimensions, C ( D ) is a Wilson-Wegner loop around a closedpath [24, 41], which shows an exponential scaling with the perimeter of the loop in theparamagnetic phase and with the enclosed area in the ferromagnetic phase of the orig-inal model [41, 60]. However, the detection of the non-local order parameter requiresthe measurement of the x component of the spin operator, which in turn requires thesingle-site resolved detection of the phase coherence between superposition states withdifferent on-site occupation numbers, a technique that is so far not available. Acknowledgements
We acknowledge useful comments by Nigel Cooper and Senthil Todadri. This workhas been supported by the DFG Forschergruppe 801.
Appendix A. Details about the perturbative calculation
According to the program outlined in Section 4, calculting the POP to n th orderin J/U amounts to calculating the operators ˆ S , . . . , ˆ S n − . Hence, for the second orderresult, we only need to calculate one operator. One finds ˆ S = (cid:88) D ,D (cid:48) D − D ˆ P D ˆ T ˆ P D , (A.1)where the sums go over the eigenvalues D of the operator ˆ D and ˆ P D is the projector onthe eigenspace corresponding to this eigenvalue. Moreover, here and in the following, aprime on a sum indicates that values which make the denominator vanish are excludedfrom the sum.One may imagine the calculation of the expectation (cid:104)O ( L ) (cid:105) as summing over all pos-sible hopping processes where the order in J/U indicates the number of occurring hops.For example, at second order, starting from a situation with uniform filling correspondingto | Φ (cid:105) , all possible processes consist of a single particle hopping to a neighboring siteand back again. Different contributions stem from the position of the starting lattice siteand its neighbor relative to the domain boundary and from the position of the operator ˆ N ( D ) in the product of operators. 17he next non-vanishing contribution to the POP is of fourth order, so we need tocalculate ˆ S , . The former turns out to be given by ˆ S = (cid:88) D ,D ,D (cid:48)
12 ( D − D ) (cid:18) D − D − D − D (cid:19) ˆ P D ˆ T ˆ P D ˆ T ˆ P D + (cid:88) D ,D (cid:48) D − D ) (cid:16) ˆ P D ˆ T ˆ P D ˆ T ˆ P D − ˆ P D ˆ T ˆ P D ˆ T ˆ P D (cid:17) . (A.2)The expression for ˆ S is too unwieldy to reproduce here, so we limit ourselves to adescription of the involved hopping processes. They can be grouped into four types: • Back- and forth hoppings of two particles separated by more than two lattice sites,i.e., independent second order processes. • Processes where the same particle or hole hops twice before returning to its originalsite. • Processes where a particle hops twice and is then followed by the created hole orvice versa. • Processes where two particles start or end on the same lattice site before hoppingback (cf. Fig. A.2).
Figure A.2: (Color online) Possible processes at fourth order where two particles start or end on thesame lattice site. The images show the configuration after two hopping events (the remaining two restoreuniform filling) for ¯ n = 4 . Summing over all contributions finally yields Eq. (15).
References [1] R. Haag,
Local Quantum Physics (Springer, 1996).[2] X.-G. Wen,
Quantum Field Theory of Many-body Systems (Oxford University Press, 2004).[3] J. Fröhlich and U. M. Studer, Rev. Mod. Phys. , 733 (1993).[4] D. H. Lee and R. Shankar, Phys. Rev. Lett. , 1490 (1990).[5] F. Gebhard, The Mott Metal-Insulator Transition - Models and Methods , no. 137 in Springer Tractsin Modern Physics (Springer, 1997).[6] P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. , 17 (2006).[7] M. Greiner, O. Mandel, T. Esslinger, T. Hänsch, and I. Bloch, Nature , 39 (2002).[8] N. Gemelke, X. Zhang, C.-L. Hung, and C. Chin, Nature , 995 (2009).[9] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Nature , 68(2010).
10] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet, and M. Greiner,Science , 547 (2010).[11] E. Berg, E. G. Dalla Torre, T. Giamarchi, and E. Altman, Phys. Rev. B , 245119 (2008).[12] M. Endres, M. Cheneau, T. Fukuhara, C. Weitenberg, P. Schauß, C. Gross, L. Mazza, M. C. Bañuls,L. Pollet, I. Bloch, et al., Science , 200 (2011).[13] J. Fröhlich and T. Spencer, Commun. Math. Phys. , 527 (1981).[14] P. J. Forrester, J. Phys. A , L143 (1986).[15] D. B. Abraham, in Phase transitions and critical phenomena , edited by C. Domb and J. L. Lebowitz(Academic Press, London, 1986), vol. 10.[16] M. E. Peskin, Ann. Phys. , 122 (1978).[17] M. P. A. Fisher and D. H. Lee, Phys. Rev. B , 2756 (1989).[18] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. , 3108 (1998).[19] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B , 546 (1989).[20] S. Sachdev, Quantum phase transitions (Cambridge University Press, Cambridge, 1999).[21] I. Bloch, J. Dalibard, and W. Zwerger, Rev. of Mod. Phys. , 885 (2008).[22] S. Trotzky, L. Pollet, F. Gerbier, U. Schnorrberger, I. Bloch, N. Prokof’ev, B. Svistunov, andM. Troyer, Nature Physics , 998 (2010).[23] S. Fölling, F. Gerbier, A. Widera, O. Mandel, T. Gericke, and I. Bloch, Nature , 481 (2005).[24] F. Wegner, J. Math. Phys. , 2259 (1971).[25] E. Fradkin and L. Susskind, Phys. Rev. D , 2637 (1978).[26] H. Tasaki, Phys. Rev. Lett. , 798 (1991).[27] E. G. Dalla Torre, E. Berg, and E. Altman, Phys. Rev. Lett. , 260401 (2006).[28] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. , 5040 (1998).[29] G. E. Astrakharchik, R. Combescot, and L. P. Pitaevskii, Phys. Rev. A , 063616 (2007).[30] M. Klawunn, A. Recati, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A , 033612 (2011).[31] S. Ejima, H. Fehske, and F. Gebhard, Europhys. Lett. , 30002 (2011).[32] S. Rachel, N. Laflorencie, H. F. Song, and K. Le Hur, Phys. Rev. Lett. , 116401 (2012).[33] H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Phys. Rev. B , 035409(2012).[34] B. Swingle and T. Senthil, ArXiv e-prints (2011), .[35] M. A. Metlitski and T. Grover, ArXiv e-prints (2011), .[36] L. D. Landau and E. M. Lifshitz, in Statistical Physics, Part 2 (Butterworth-Heinemann, 1980),vol. IX.[37] P. G. J. van Dongen, Phys. Rev. B , 7904 (1994).[38] A. B. Harris and R. V. Lange, Phys. Rev. , 295 (1967).[39] J. Villain, J. Phys (Paris) , 581 (1976).[40] K. G. Wilson, Phys. Rev. D , 2445 (1974).[41] J. B. Kogut, Rev. Mod. Phys. , 659 (1979).[42] W. Zwerger, Europhys. Lett. , 421 (1989).[43] D. S. Fisher and J. D. Weeks, Phys. Rev. Lett. , 1077 (1983).[44] S. T. Chui and J. D. Weeks, Phys. Rev. B , 4978 (1976).[45] P. Minnhagen and G. G. Warren, Phys. Rev. B , 2526 (1981).[46] W. Zwerger, Z. Phys. B , 111 (1990).[47] C. Dasgupta, B. I. Halperin, Phys. Rev. Lett. , 1556 (1981).[48] I. F. Herbut, Phys. Rev. B , 13729 (1998).[49] L. Balents, L. Bartosch, A. Burkov, S. Sachdev, and K. Sengupta, Progress of Theoretical PhysicsSupplement , 314 (2005).[50] I. Herbut, A Modern Approach to Critical Phenomena (Cambridge University Press, 2007).[51] B. Capogrosso-Sansone, S. G. Söyler, N. Prokof’ev, and B. Svistunov, Phys. Rev. A , 015602(2008).[52] A. Zee, Quantum Field Theory in a Nutshell (Princeton University Press, 2010), 2nd ed.[53] R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, Nature , 204 (2008).[54] U. Schneider, L. Hackermüller, S. Will, T. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch,and A. Rosch, Science , 1520 (2008).[55] A. Montorsi and M. Roncaglia, Phys. Rev. Lett. , 236404 (2012).[56] H. V. Kruis, I. P. McCulloch, Z. Nussinov, and J. Zaanen, Europhys. Lett. , 512 (2004).[57] H. V. Kruis, I. P. McCulloch, Z. Nussinov, and J. Zaanen, Phys. Rev. B , 075109 (2004).[58] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner, Nature , 307 (2011).[59] S. Pielawa, T. Kitagawa, E. Berg, and S. Sachdev, Phys. Rev. B , 205135 (2011).
60] R. Savit, Rev. Mod. Phys. , 453 (1980)., 453 (1980).