Table Of Content

Merging and alignment of Dirac points in a shaken honeycomb optical lattice Selma Koghee1,2, Lih-King Lim2, M.O. Goerbig2, and C. Morais Smith1 1Institute for Theoretical Physics, Utrecht University, Leuvenlaan 4, 3584 CE Utrecht, The Netherlands and 2Laboratoire de Physique des Solides, CNRS UMR 8502, Université Paris-Sud, 91405 Orsay, France (Dated: February 1, 2012) Inspired by the recent creation of the honeycomb optical lattice and the realization of the Mott insulatingstateinasquarelatticebyshaking,westudyheretheshakenhoneycombopticallattice. Foraperiodicshakingofthelattice,aFloquettheorymaybeappliedtoderiveatime-independent Hamiltonian. In this effective description, the hopping parameters are renormalized by a Bessel 2 function,whichdependsontheshakingdirection,amplitudeandfrequency. Consequently,thehop- 1 ping parameters can vanish and even change sign, in an anisotropic manner, thus yielding different 0 band structures. Here, we study the merging and the alignment of Dirac points and dimensional 2 crossoversfromthetwodimensionalsystemtoonedimensionalchainsandzerodimensionaldimers. n Wealsoconsidernext-nearest-neighborhopping,whichbreakstheparticle-holesymmetryandleads a to a metallic phase when it becomes dominant over the nearest-neighbor hopping. Furthermore, J we include weak repulsive on-site interactions and find the density profiles for different values of 1 the hopping parameters and interactions, both in a homogeneous system and in the presence of 3 a trapping potential. Our results may be experimentally observed by using momentum-resolved Raman spectroscopy. ] s a g I. INTRODUCTION theatomicspecies,aswellastheinteractionscanbeengi- - neeredwithahighdegreeofprecision. Themoreinvolved t n triangularandhoneycombgeometrieswererecentlyreal- a ized experimentally and exotic correlated states of mat- u The study of Dirac points, i.e. the contact points be- ter have been observed experimentally [8] or predicted q tween different energy bands with an approximate linear theoretically [9–12]. . dispersion relation, has become a major issue since the t a experimental breakthrough in graphene-based electron- Theapplicationofatime-periodicperturbationonthe m optical lattice introduces yet another parameter scale ics [1, 2]. Indeed, the low-energy electronic properties of - graphene are governed by a pseudo-relativistic 2D Dirac intothesystem. Aperiodicshakingoftheopticallattice, d equation for massless fermions situated at the K and K(cid:48) uptothekHzfrequencyrange,hasbeenimplementedby n placing one of the mirrors used to create the optical lat- corners of the Brillouin zone [3]. The Dirac points are o tice on a piezoelectric material, such that the mirror can c topologically protected and a gap is opened only when [ theinversionsymmetryofthelatticeorthetime-reversal be moved back and forth in the direction of the beam symmetry are broken. [13, 14]. The Floquet formalism shows that the hopping 2 energy of the atoms in the shaken lattice is renormalized v The possibility to generate topological phase transi- by a Bessel function, as a function of the shaking fre- 0 tions in graphene-like systems has recently attracted a quencyandamplitude,thusallowingboththemagnitude 1 great deal of attention. Within a tight-binding descrip- 2 and the sign of the hopping parameter to change. This tion, an anisotropy in the nearest-neighbor hopping pa- 2 rathercounter-intuitivephenomenon,ascomparedtothe rameters makes the Dirac points move away from the . standard tight-binding physics, has been experimentally 1 high-symmetry K and K(cid:48) points and, under appropri- 1 ateconditions,mergeattime-reversalinvariantpointsin observed in a one-dimensional cold-atomic system [14]. 1 the first Brillouin zone [4–6]. Most saliently, this merg- In this paper, we consider ultracold fermions trapped 1 ingofDiracpointsisassociatedwithatopologicalphase in a shaken honeycomb optical lattice. Within the Flo- : v transition between a semimetallic phase and a gapped quet formalism, we derive an effective Hamiltonian that Xi band-insulating phase. An experimental investigation of generalizesthatofagraphene-likematerialunderstrain. the mergingtransition in grapheneturns outto beprob- In particular, we find that the alignment and merging r a lematic, since in order to appropriately modify the hop- of Dirac points in momentum space are now accessible pingparameters,anunphysicallylargestrainneedstobe with ultracold fermions in the shaken optical lattice and applied to the graphene sheet [7]. thephasediagramconsistsofvariousphasesofthecorre- spondingsolid-statesystemthatareotherwisedifficultto An alternative system for the study of such topolog- realize. Furthermore, bytakingintoaccountaHubbard- ical transitions is that of ultracold atoms trapped in a likeinteractionforspinfulfermions,westudythedensity honeycomb optical lattice. Since the seminal realization profiles for the homogenous and the trapped gas within of the superfluid-Mott-insulator transition in the Bose- a Hartree-Fock theory. Hubbard model, ultracold atoms in optical lattices have become promising systems to emulate condensed-matter The outline of this paper is the following: in Sec.IIA physics. Indeed,thelatticegeometry,thedimensionality, we introduce the time-dependent Hamiltonian and in 2 Sec.IIB we derive the time-independent one, by apply- energy gained in hopping to the nn and next-nearest- ing the Floquet formalism. In Sec.III we investigate the neighbor (nnn) sites, respectively, and µ is the on-site merging and alignment of Dirac points, when the optical energy. We remark that the nnn hopping is taken into lattice is shaken along specific directions. The descrip- accountbecausethennhoppingmayberenderedvanish- tion is extended to include interactions in Sec.IV, where inglysmallintheeffectivetime-independentdescription. we derive the dependence of the density on the chemical Inthisregime,thennnmaybecomethedominantkinetic potential. Implicationsofourresultsforexperimentsare term. Inasquarelattice,wherethepotentialisseparable discussedinSec.V.Finally,ourconclusionsarepresented in independent eˆ and eˆ components, the nnn hopping x y in Sec.VI. isidenticallyzero[15]. However,thennnhoppingcanbe nonzero in the honeycomb lattice, since its potential is not separable in eˆ and eˆ components. Nevertheless, it x y II. THE SHAKEN HONEYCOMB LATTICE may be expressed as the sum of two triangular lattices. In this section, we derive a time-independent effective description for ultracold atoms trapped in a periodically shaken honeycomb optical lattice by utilizing a d2 Foflosqinugelte-tchoemorpyo.neFnotrfesrimmipolnicicitya,towmesfoacnudscoonnsaidesrysotnemly d3 d1 single-particle terms in this section. The results from theFloquettheoryarevalidforfermionicatomswithin- ternaldegreesoffreedomaswellasforbosonicatoms. In particular, the hyperfine state of fermionic atoms, play- ingtheroleofaneffectivespin-1/2degreeoffreedomfor electrons, will be considered when interaction effects are taken into account in Sec. IV. Shaking direction 1 Shaking direction 2 A. Time-dependent Hamiltonian FIG.1: (Color online) Laserconfiguration tocreatethe honeycomb lattice, which consists of two triangular sublattices In the tight-binding limit, the system of ultracold (A, red dots, and B, blue dots). The vectors d , d , and d fermionic atoms trapped in a 2D shaken honeycomb op- 1 2 3 connectasiteontheAsublatticetoitsnearestneighborson tical lattice can be described by the Hamiltonian the B sublattice. H(t)=H +W(t), (1) 0 The time-dependent part of the Hamiltonian (1), which consists of two distinct parts. The static part (cid:32) (cid:33) (cid:88) (cid:88) W(t)=mΩ2cos(Ωt) r·ρa†a + r·ρb†b , 3 r r r r (cid:88)(cid:88)(cid:16) (cid:17) H =−γ a†b +b† a r∈A r∈B 0 r r+dj r+dj r (4) j=1r∈A describestheharmonicshakingofthelatticeinthedirec- (cid:88)3 (cid:88)3 (cid:18)(cid:88) (cid:88) (cid:19) tionρwithadrivingfrequencyΩintheco-movingframe −γ(cid:48) a†a + b†b r r+di−dj r r+di−dj of Ref. [16]. As a consequence of the transformation to i=1j=1,j(cid:54)=i r∈A r∈B the co-moving frame, W(t) describes atoms of mass m (cid:18) (cid:19) (cid:88) (cid:88) experiencing a position-dependent sinusoidal force. −µ a†a + b†b (2) r r r r r∈A r∈B B. Effective Hamiltonian is simply the tight-binding Hamiltonian in the honeycomb lattice, where a† (b†) and a (b ) are, respectively, r r r r fermionic creation and annihilation operators on the lat- The unavoidable complexity that arises when dealing tice site r in the A (B) sublattice. The three vectors with a quantum many-body system out of equilibrium has recently motivated the development of new theoret- d √ d √ ical tools, for example time-dependent density matrix d =deˆ , d = (−eˆ + 3eˆ ), d = (−eˆ − 3eˆ ), 1 x 2 2 x y 3 2 x y renormalization group [17], time-dependent dynamical (3) mean field theory [18], and exact diagonalization [19]. connect an A-lattice site with its three nearest-neighbor However, for a periodically driven quantum system, (nn)B-latticesitesandaregivenintermsofthedistance the Floquet theory offers a simplified description of √ d=8π/3 3k betweennnsites,wherek isthelaserwave the system, in the form of a time-independent effective number (see Fig.1). Here, γ,γ(cid:48) > 0 characterize the Hamiltonian,iftheperiodT =2π/Ωistheshortesttime 3 scaleintheproblem[20]. Inthislimit, theatomscannot states in the first Brillouin zone, i.e. states with quasi- follow the shaking motion adiabatically and remain thus energies −(cid:126)Ω/2<(cid:15) ≤(cid:126)Ω/2. φ attheiraveragelatticeposition,albeitwithrenormalized The space in which the states |φ(q,t)(cid:105) are defined is hoppingparameters. Thesystemisthusconsideredtobe the composite of the Hilbert space spanned by square in a stationary state and the knowledge of equilibrium integrable functions on configuration space, |α(q)(cid:105), and physics can be employed. thespaceofT-periodicfunctions. Thestate|φ(q,t)(cid:105)may Let us consider the Floquet Hamiltonian defined by bewrittendowninanorthonormalbasisinthecomposite HF =H(t)−i(cid:126)∂t, where H(t+T)=H(t) is periodic in space according to time [20]. The eigenvalue equation is then given by ∞ HF|φ(q,t)(cid:105)=(cid:15)φ|φ(q,t)(cid:105), (5) |φ(q,t)(cid:105)= (cid:88)(cid:88)c exp[−iFˆ(t)+inΩt]|α(q)(cid:105), (6) n,α n=0 α where(cid:15) isthequasienergydefineduniquelyuptoamul- φ tipleof(cid:126)Ω. Anysolution|φ(q,t)(cid:105)ispartofasetofsolu- where c are coefficients to normalize |φ(q,t)(cid:105) and the tions exp(inΩt)|φ(q,t)(cid:105) with integer n, which all corre- n,α operatorFˆ(t)canbeanyT-periodicHermitianoperator. spondtothesamephysicalsolution. Hence,thespectrum of the Floquet Hamiltonian possesses a Brillouin-zone- Therefore, wecanconvenientlychoose Fˆ(t)tobeFˆ(t)= like structure [20]. The interest therefore lies with the (cid:126)−1(cid:82)tdt(cid:48)W(t(cid:48)), such that H(t)−(cid:126)∂ Fˆ(t)=H . 0 t 0 If the condition (cid:104)α(cid:48)(q)|(cid:104)exp[iFˆ(t)] exp[i(n−n(cid:48))Ωt]H exp[−iFˆ(t)](cid:105) |α(q)(cid:105)(cid:28)(cid:126)Ω (7) 0 T is satisfied for any two states |α(q)(cid:105) and |α(cid:48)(q)(cid:105), then the eigenvalues (cid:15) are approximately φ (cid:15) =(cid:104)φ(q,t)|(cid:104)H (cid:105) |φ(q,t)(cid:105)≈ (cid:88)c c (cid:104)α(cid:48)(q)|(cid:104)exp[iFˆ(t)]H exp[−iFˆ(t)](cid:105) |α(q)(cid:105). (8) φ F T 0,α(cid:48) 0,α 0 T α,α(cid:48) Here,(cid:104)O(t)(cid:105) =T−1(cid:82)T dtO(t)denotesthetimeaverage condition (7) is satisfied if γ (cid:28) (cid:126)Ω, and the effective T 0 of the operator O(t) over the period T. The condition Hamiltonian becomes (7) will hold for n (cid:54)= n(cid:48), if H is nearly constant during 0 3 the period T, which is small if Ω is large. In this case, (cid:88)(cid:88) (cid:16) (cid:17) H =− γ a†b +b† a states with different n do not mix. If Ω is large enough, eff j r r+dj r+dj r suchthatthecondition(7)alsoholdsforn=n(cid:48),thenthe j=1r∈A energy spectrum will split up into energy bands labelled (cid:88)3 (cid:88)3 (cid:18)(cid:88) (cid:88) (cid:19) − γ(cid:48) a†a + b†b by an index n, where the details within the energy band i,j r r+di−dj r r+di−dj aredeterminedbyH . Becausethestateswithadifferent i=1j=1,j(cid:54)=i r∈A r∈B 0 (cid:18) (cid:19) indexnareseparatedbyanenergywhichisamultipleof (cid:88) (cid:88) −µ a†a + b†b , (10) (cid:126)ΩandbecausethespectrumpossessesaBrillouin-zone- r r r r likestructure,onlythetermswithn=0needtobetaken r∈A r∈B intoaccount. TheeffectiveHamiltonianH ,whichgives eff where the renormalized nn hopping parameters γ are j risetothesamespectrumastheFloquetHamiltonian,is given by then defined by [21] (cid:18)(cid:12) (cid:12)(cid:19) H =(cid:28)exp[iFˆ(t)]H exp[−iFˆ(t)](cid:29) γj =γJ0 (cid:12)(cid:12)(cid:12)dj ·ρm(cid:126)Ω(cid:12)(cid:12)(cid:12) , (11) eff 0 T =(cid:28)(cid:88)∞ inn![Fˆ(t),H0]n(cid:29) . (9) abnyd the renormalized nnn hopping parameters are given n=0 T (cid:18)(cid:12) (cid:12)(cid:19) Here, [Fˆ,Gˆ]n denotes the multiple commutator, which is γi(cid:48),j =γ(cid:48)J0 (cid:12)(cid:12)(cid:12)(di−dj)·ρm(cid:126)Ω(cid:12)(cid:12)(cid:12) (12) defined by [Fˆ,Gˆ] =[Fˆ,[Fˆ,Gˆ] ] and [Fˆ,Gˆ] =Gˆ. n+1 n 0 Effective Hamiltonians corresponding to Eq.(9) have (see Appendix A for detailed calculations). In these ex- been derived for linear shaking of a one-dimensional lat- pressions, J (x) denotes the zeroth order Bessel func- 0 tice [22] and for elliptical shaking of a triangular lattice tion of the first kind, which shows a damped oscillation [23]. For the shaken honeycomb lattice studied here, the around zero. 4 In terms of the renormalized nn and nnn hopping pa- zone, start to move in the q -direction along the verti- y rameters,thediaginalizationoftheeffectiveHamiltonian cal edges of the latter. This motion is depicted by the (10) yields the dispersion relation arrows in Fig.3(b). Even if the two Dirac points are no longer located at the high-symmetry points K and (cid:15)λ(q)=h(q)+λ|f(q)|, (13) K(cid:48),theyremainrelatedbytime-reversalsymmetry,such that their Berry phases π and −π are opposite. This where λ=± is the band index, and we have defined the non-zero Berry phase topologically protects each of the functions Dirac points and thus the semimetallic phase remains (cid:88) robust until γ2,3 = γ1/2, where the two points merge at f(q)= γ exp(−iq·d ) (14) j j a time-reversal invariant momentum, i.e. half of a re- j ciprocal lattice vector [6]. In the present example, this point is situated at the center of the vertical edges of and thefirstBrillouinzone,andthebanddispersionbecomes h(q)=2(cid:88)γ(cid:48) cos[q·(d −d )]. (15) parabolicinthey-directionwhileremaininglinearinthe i,j i j x-direction [see Fig.2(b)]. The merged Dirac points are i<j no longer topologically protected due to the annihilation oftheoppositeBerryphases. Consequently,afurtherin- creaseoftheshakingamplitude,whichresultsinafurther III. MERGING AND ALIGNMENT OF DIRAC POINTS decrease of γ2,3, leads to the opening of a gap between thetwobands. Thus,thesystemundergoesatopological phase transition from a semimetal to a band insulator. Inthissection,thehoneycomblatticewithanisotropic Thismergingtransitionwasalsostudiedinastaticsetup hopping is studied. In the first two subsections, only nn in Ref. [24], where the hopping amplitudes γ were pro- hopping is considered for illustration reasons. Indeed, posed to be modified by a change in the intensity of one this allows for a simple understanding of the main con- ofthelasersusedtocreatetheopticallattice. Incontrast sequences of shaking on Dirac-point motion and dimen- tothisstaticsetup,shakingthehoneycomblatticeallows sional crossover. In Subsec. IIIC, we discuss how the one to completely annihilate some of the nn hopping pa- picture evolves when nnn hopping is included, and the rametersandtoevenchangetheirsign. Thissignchange signofthehoppingparametersisinvestigatedinSubsec. occurs at the zeros of the Bessel function [see Eq.(11)]. IIID.Sincethesystemwithtwonnequalhoppingparam- For an example system of 40K atoms in a lattice created eters and a single independent one captures the essential bylaserswithawavelengthof830nm,whichisshakenin features of the systems with three independent nn hop- the direction perpendicular to d , the situation γ =0 pingparameters,wewillfocusonthissystem. Thenum- 1 2,3 is encountered for bering of the γ s is chosen such that |γ | = |γ | = γ , j 2 3 2,3 which can be achieved by shaking in a direction parallel ρ=180nm; Ω/2π =6kHz, (16) or perpendicular to d . 1 Although the atoms in the optical lattice are charge- which corresponds to the first zero of the Bessel func- neutral objects, we shall adopt the language from tion. At this particular point, and if γ(cid:48) = 0 in addition, condensed-matterphysicsandcallazero-gapphasewith the system consists of a set of effectively decoupled hor- a pair of Dirac cones and a vanishing density of states at izontal bonds along which the atoms are solely allowed the band-contact points a semimetal, whereas a gapped to hop. This yields two flat bands at ±γ [see Fig.2(c)] 1 phase is called band insulator. Furthermore, nnn hop- that may be viewed as the extreme limit of the band- ping induces a metallic phase for small values of γ be- insulating phase. Alternatively, one may view this situ- j cause of an overlap between the two bands that yields a ation upon decreasing the value of γ as a dimensional 2,3 non-vanishing density of states at the energy level of the crossoverfroma2Dbandinsulatortoazero-dimensional band-contact points. (0D) system. A small non-zero value of γ simply pro- 2,3 vides a weak dispersion of these decoupled bands (not shown). A. Merging of Dirac points Ifthelaticeisshakeninthedirectionperpendicularto B. Alignment of Dirac points d (direction 1 in Fig.1), γ remains equal to γ, whereas 1 1 γ and γ are renormalised to a smaller value. An in- Anotherdimensionalcrossover,from2Dto1D,maybe 2 3 crease in the shaking amplitude results in a decrease in obtained if the lattice is shaken in the direction parallel γ = γ = γ , which is depicted by the arrow M in to one of the nn vectors (direction 2 in Fig.1). Here, 2,3 2 3 Fig.3(a). When the hopping parameters change accord- we choose d to maintain the symmetry γ = γ = 1 2,3 2 ing to this arrow M, the energy spectrum evolves from γ . In this case, both γ and γ are renormalized by 3 1 2,3 Fig.2(a) to Fig.2(b). The Dirac points, originally sit- Bessel functions, albeit with different arguments. Since uated at the corners K and K(cid:48) of the first Brillouin all hopping parameters are renormalized, the trajectory 5 (a) (b) forthesamesystemof40Katomsmentionedabove. Here, for illustrative purposes, a simplified trajectory of the system is depicted in Fig.2 by arrow A, which corresponds to the motion of the Dirac points in reciprocal space as shown in Fig.3(c). As γ approaches zero, 1 the Dirac points align in lines parallel to the x-axis at √ q =±π/ 3dandtheenergybarriersbetweenthealign- y ing points are lowered. Consequently, when γ = 0, 1 (c) (d) theenergyspectrumcontainslineswherethetwoenergy band meet and the dispersion is linear, as is shown in Fig.2(d). The dispersion relation (13) then reads simply (cid:12) (cid:32)√ (cid:33)(cid:12) (cid:12) 3 (cid:12) (cid:15) (q)=2λγ (cid:12)cos q d (cid:12), (18) λ 2,3(cid:12) 2 y (cid:12) (cid:12) (cid:12) and one clearly sees the 1D character. Indeed, there is FIG.2: (Coloronline)Energydispersionfortheshakenhon- no dispersion in the q -direction, as is also evident from eycomb optical lattice, with k = 1, γ = 1, and γ(cid:48) = 0. The x Fig.2(d), and the system may be viewed as completely labels q and q represent the x and y components of the x y decoupled 1D chains in which the zig-zag arrangement momentum,respectively. Thexandy axeshavebeenchosen suchthatthennvectorsaregivenbyEq.(3). (a)Theisotropic is of no importance. In this particular limit, the sites A case, where γ = γ . (b) The merged Dirac points, where and B are therefore no longer inequivalent such that the 1 2,3 γ = γ /2. (c) The zero dimensional case, where γ = 0. unit cell is effectively divided by two, and the size of the 2,3 1 2,3 (d) The aligned Dirac points, where γ =0. first Brillouin zone is consequently doubled. The aligned 1 Diracpointsmaythus,alternatively,beviewedasdueto an artificial folding of the second (outer) half of the first (b) M (a) Brillouin zone into its inner half. However, this aspect M is very particular in that the Brillouin zone immediately Band retrieves its original size when γ is small, but non-zero, insulator 1 orifnnnhoppingsaretakenintoaccount. Inbothcases, oneneedstodistinguishthetwodifferentsublatticesand A one obtains a dispersion in the q -direction. x The actual behavior of the system for an increasing Zero-gap (c) A semimetal shaking amplitude is discussed in Sec.V. This behavior is more complicated because all three nn hopping pa- rametersarerenormalized,which,beyondthealignment, 0D leads to the merging of Dirac points and the opening of Merging of Dirac points agapalsointhecaseofshakingparalleltooneofthenn Isotropic system vectors. In the absenceof nnn hopping, the 0D limitcan Alignment of Dirac points be reached in addition. FIG.3: (Coloronline)(a)Phasediagram,showingthephase transition between the zero-gap semimetallic phase and the C. Next-nearest-neighbor hopping insulatingphase,whichhappensatγ =2γ . Here,wehave 1 2,3 chosen γ(cid:48) = 0. (b) Dirac-point motion in the first Brillouin zoneforashakingdirectionperpendiculartod [directionM The major consequence of nnn hopping is to break 1 inthephasediagram(a)]. (c)Dirac-pointmotioninthefirst particle-hole symmetry, as may be seen from Eq. (13), Brillouinzoneforashakingdirectionparalleltod [direction where non-zero values of γ(cid:48) yield (cid:15) (q)(cid:54)=−(cid:15) (q). Its 1 i,j λ −λ Ainthephasediagram(a)]. Thecontourplotsdepictthedis- relevance depends sensitively on the shaking direction, persion of the isotropic system with an arbitrary color scale. because of the different renormalization of the nn hop- The area with higher contrast is the first Brillouin zone. ping parameters. The band structure with nnn hopping included is depicted in Fig.4 for different shaking directions. of the system in the phase space upon increasing the shaking amplitude is not a straight line, as was the case for shaking perpendicular to a hopping direction, and 1. Shaking in the direction perpendicular to d 1 has a new feature: the alignment of Dirac points, which occurs for γ1 =0. The first zero of γ1 is found at In the case of a shaking perpendicular to d1, only γ2,3 are decreased, whereas γ = γ remains the leading en- 1 ρ=92nm; Ω/2π =6kHz, (17) ergy scale in the band structure [32]. The band struc- 6 ture for the unshaken lattice is depicted in Fig.4(a) for (a) (b) γ(cid:48)/γ = 0.1, and one notices that the main features of the band structure, namely the Dirac points, are unal- tered with respect to the case γ(cid:48) = 0 in Fig.2(a), apart from the flattening of the upper band as compared to the lower one. When approaching the merging transition γ =γ /2, thevalueof whichisdeterminedbythe 2,3 1 zeros of f(q) in Eq. (14) and that therefore does not (c) (d) depend on the nnn hopping parameters, the band width remainsdominatedbythelargesthoppingparameterγ , 1 such that the band structure [Fig.4(b)] at the transition is essentially the same as in Fig.2(b) for γ(cid:48) = 0. In the 0D limit, with γ = 0 the originally flat bands 2,3 [Fig.2(c)]acquiretheweakdispersionofatriangularlat- ticeasaconsequenceofthenon-zeronnnhoppingparam- eters. However, as expected from the above arguments, (e) (f) the dispersion is on the order of γ(cid:48) and thus small as compared to the energy separation ∼ 2γ = 2γ between 1 the two bands. 2. Shaking in the direction parallel to d1 FIG.4: (Coloronline)Energydispersionfortheshakenhon- eycomb optical lattice, with k = 1, γ = 1, and γ(cid:48) = 0.1. The labels q and q represent the x and y components of x y Incontrasttoashakingdirectionperpendiculartod1, the momentum, respectively. The x and y axes have been nnn hopping has more drastic consequences if the lat- chosen such that the nn vectors are given by Eq.(3). (a) tice is shaken in the direction parallel to d . In this The homogeneous case, where γ = γ . (b) The merged 1 1 2,3 case, all nn hopping parameters are decreased, and the Dirac points, where γ2,3 = γ1/2. (c) The zero dimen- relative importance of nnn hopping is enhanced. No- sional case, where γ2,3 = 0. (d) The aligned Dirac points, tice further that the nnn lattice vectors ±(d −d ) are where γ1 = 0. (e) An example of the metallic phase with 2 3 ρ=5.2((cid:126)/mΩd)eˆ . (f)Anotherexampleofthemetallicphase now perpendicular to the shaking direction such that x with ρ=4.8((cid:126)/mΩd)eˆ . γ(cid:48) = γ(cid:48) = 0.1γ remains unrenormalized. Also in this x 2,3 case,thesystemisapproachingthe1Dlimit,withγ =0 1 [see Fig.4(d)]. However, in contrast to Fig.2(d), the chains remain coupled by nnn hopping that yields a dispersion in the q -direction. Furthermore, as mentioned x above, the A and B sites are now not equivalent from a crystallographic point of view, such that the outer parts D. The signs of the hopping parameters of the first Brillouin zone cannot be folded back into the inner one, as may be seen from Fig.4(d). Finally, for particular values of the shaking amplitude in the direction parallel to d , the nn hopping As already alluded to in the previous sections, the 1 parameters can be decreased in such a manner as to shakingofahoneycomblatticecanleadtoasignchange render γ more relevant. In this case, the two bands of the hopping parameters. Quite generally, Fig.5 shows 2,3 can overlap in energy, as depicted in Figs.4(e) and 4(f) that changing the relative signs of the nn hopping pa- for ρ = 5.2((cid:126)/mΩd)eˆ (in which case γ ≈ γ ) and rameters results in a translation of the energy spectrum x 1 2,3 ρ = 4.8((cid:126)/mΩd)eˆ (with γ ≈ 0), respectively. In the in momentum space. This effect was also mentioned in x 2,3 latter example there are no band contact points, in spite Ref. [4]. Indeed, the relative signs determine at which of of the overlap between the two bands, and the system thefourtime-reversalinvariantmomentainthefirstBril- would be in an insulating phase if nnn hopping terms louinzonethemergingofDiracpointsandthesemimetal- were not taken into account. This overlap in energy be- insulator transition take place when |γ |=2|γ |. How- 1 2,3 tween the two bands yields a non-zero density of states ever, the sign change of the nn hopping parameters can at any energy, such that the semi-metallic (or insulat- betransformedawaybyagaugetransformation[4]. Nev- ing) phase vanishes and yields, at half-filling, a metallic ertheless, the sign of the nnn hopping parameter is im- phase with particle and anti-particle pockets in the first portant, since it determines whether the upper or the Brillouin zone. lower band is flattened. 7 a b small-U large-frequencylimit,criticalresonancesmayoc- cur for m (cid:29) n [26]. Nevertheless, it has been shown in a b Ref. [26] that these resonances, which occur in higher- orderperturbationtheory,arestronglysuppressedinthe e f c d large-m limit. In the weakly-interacting regime considered here, the ground state is adiabatically connected to that of the e f non-interactingsystem,withnobrokensymmetry. First, we use the Fourier transform of the creation and annihi- g h lationoperators,a =N−1/2(cid:80) exp(iq·r)a ,tofind r,σ q q,σ theHamiltonianinmomentumspace. WithinaHartree- g h c d Fock theory, we introduce a mean-field decoupling of the Band insulator interaction terms, Zero-gapsemimetal a† a† a a ≈ (20) q1,σ q2,σ(cid:48) q3,σ(cid:48) q4,σ FIG. 5: (Color online) Phase diagram and contour plots of (cid:10)a† a (cid:11)a† a −(cid:10)a† a (cid:11)a† a the energy bands, showing the effects of the renormalized nn q2,σ(cid:48) q3,σ(cid:48) q1,σ q4,σ q2,σ(cid:48) q4,σ q1,σ q3,σ(cid:48) herogpypibnagndpsa,rawmheetreertshγejc.ol(oar)stcoal(ihn)gCisoanrtboiutrraprylotasndoftthheefiernst- +a†q2,σ(cid:48)aq3,σ(cid:48)(cid:10)a†q1,σaq4,σ(cid:11)−a†q2,σ(cid:48)aq4,σ(cid:10)a†q1,σaq3,σ(cid:48)(cid:11) Brillouin zone is the area with higher contrast. The value of −(cid:10)a†q2,σ(cid:48)aq3,σ(cid:48)(cid:11)(cid:10)a†q1,σaq4,σ(cid:11)+(cid:10)a†q2,σ(cid:48)aq4,σ(cid:11)(cid:10)a†q1,σaq3,σ(cid:48)(cid:11), thennhoppingparametersforeachcontourplotaregivenby thepositionofthecorrespondingletterinthephasediagram such that the expectation values of both sides are equal. and γ(cid:48) =0. The dark regions indicate energies close to zero, For the mean value we take whereas brighter regions are further away in energy from the Fermi level at half filling. (cid:104)a† a (cid:105)=(cid:104)b† b (cid:105)=Nn δ δ , (21) q,σ q(cid:48),σ(cid:48) q,σ q(cid:48),σ(cid:48) q,σ q,q(cid:48) σ,σ(cid:48) where N is the number of sites per sublattice, n is q,σ IV. INTERACTIONS the density of atoms with momentum q and spin index σ, and δ is the Kronecker delta. We then obtain the α,α(cid:48) Until now, we have considered single-component mean-field Hamiltonian fermionic atoms and, due to the Pauli principle, the absence of s-wave interaction naturally results in an ideal H =H − UNn2 + Un(cid:88)(cid:0)a† a +b† b (cid:1), Fermi lattice gas, albeit with an unusual band structure. MF eff 8 4 q,σ q,σ q,σ q,σ σ,q By trapping two hyperfine states of the fermionic atoms, (22) Hubbard-like interaction terms arise, where the total density is defined by n=(cid:80) n . q,σ q,σ The Hamiltonian (22) may be rewritten in a matrix (cid:88)(cid:88)U H = a† a† a a form: int 2 r,σ r,σ(cid:48) r,σ(cid:48) r,σ r∈Aσ,σ(cid:48) UNn2 (cid:88)(cid:88)U H =− (23) + 2 b†r,σb†r,σ(cid:48)br,σ(cid:48)br,σ, (19) MF 8 (cid:18) (cid:19)(cid:18) (cid:19) r∈Bσ,σ(cid:48) +(cid:88)(cid:0)a† b† (cid:1) h(µ,q) f(q) aq,σ , q,σ q,σ f∗(q) h(µ,q) b where the fermionic operators now acquire an additional σ,q q,σ spin index σ = {↑,↓}, which is summed over, and U is where we have introduced the functions the interaction energy. Naturally, in order to be able to apply the Floquet theory in the presence of interactions, 3 3 the Hamiltonian H0 in Eq. (7) must now be replaced by h(µ,q)= Un−µ−γ(cid:48)(cid:88) (cid:88) exp[−iq·(d −d )], (24) H =H +H . This is the case in our study because we 4 i j 0 int i=1j=1,j(cid:54)=i investigate the weak-coupling limit with U (cid:28) γ. Since theinteractiontermcommuteswiththeshaking,itisnot and f(q) is defined in Eq. (14) The Hamiltonian (23) renormalized,similarlytotheon-siteenergytermpropor- can then be diagonalized by the unitary operator tional to µ in Eq.(2). However, it has been shown that (cid:18) (cid:19) complications may arise when a multiple of the energy 1 1 if(q)/|f(q)| U is in resonance with a harmonic of (cid:126)Ω, m(cid:126)Ω = nU, Uˆ = √2 f∗(q)/|f(q)| −i , (25) for integer m and n. Whereas the limit [25] m (cid:28) n is notconsideredherebecauseitisincontradictionwiththe which yields 8 H =−UNn2 +(cid:88)(cid:0)c† d† (cid:1)(cid:18)h(µ,q)−|f(q)| 0 (cid:19)(cid:18)cq,σ(cid:19). (26) MF 8 q,σ q,σ 0 h(µ,q)+|f(q)| dq,σ σ,q Because the c and d quasiparticles are free, the partition function corresponding to the Hamiltonian (26) reads (cid:20) (cid:18) (cid:26) (cid:27) (cid:26) (cid:27)(cid:19)(cid:21) (cid:88) (cid:2) (cid:0) (cid:1)(cid:3) (cid:2) (cid:0) (cid:1)(cid:3) Z =exp log 1+exp −β h(µ,q)−|f(q)| +log 1+exp −β h(µ,q)+|f(q)| , (27) σ,q where β =(k T)−1 with k denoting the Boltzmann’s constant and T the temperature. B B The total number of particles N is given by (1/β)∂logZ/∂µ, and one obtains (cid:18) (cid:19) (cid:88) 1 1 N = + . (28) (cid:2) (cid:0) (cid:1)(cid:3) (cid:2) (cid:0) (cid:1)(cid:3) 1+exp β h(µ,q)−|f(q)| 1+exp β h(µ,q)+|f(q)| σ,q Since the expression inside the sum does not depend on spin, summing over σ yields a factor 2. One recognizes in Eq.(28) the Fermi-Dirac distribution function N (x)=[1+exp(x)]−1. The number of particles N is related to the FD density n, which is defined here as the number of particles per lattice site, i.e. n=N/2N. Converting the sum over q into an integral, the following self-consistent equation for the density is derived n(µ)= 1 (cid:90) d2q(cid:26)N (cid:2)β(cid:0)h(µ,q)−|f(q)|(cid:1)(cid:3)+N (cid:2)β(cid:0)h(µ,q)+|f(q)|(cid:1)(cid:3)(cid:27), (29) V FD FD 1BZ 1BZ where the integral is restricted to the first Brillouin zone, the surface of which is V . 1BZ In Fig.6(a), the density n(µ) is plotted for several val- honeycombopticallatticehasbeenevaluatedinRef. [24], ues of γ /γ . For the isotropic case, γ /γ = 1, the result of2Z,3hu1et al. is reproduced [27]. Fo2r,3the10D limit, γ ≈1.861E (cid:18)V0 (cid:19)3/4exp(cid:32)−1.582(cid:114)V0 (cid:33), (30) theflatlineduetothegapinthespectrumisclearlyvis- R E E R R ible at the chosen temperature. Fig.6(b) confirms that repulsive interactions lead to a lower density than in a intermsoftherecoilenergyE =(cid:126)2k2/2mandthemag- R system without interactions for the same chemical po- nitude of the potential barrier between nearest-neighbor tential. Fig.6(c) agrees with the observation that the lattice sites V . The magnitude of the nnn hopping pa- 0 nnn hopping breaks the particle-hole symmetry. This ef- rameter γ(cid:48) is not yet known, but could be determined fect is also visible in Fig.6(d), where the dependence of from numerical band structure calculations. In a typical the density on the chemical potential is calculated for a experimental situation, we expect the ratio γ(cid:48)/γ to be shaking vector where the system is in the zero-gapped in the 5−10% range, in agreement with the parameter semimetallic phase for γ(cid:48) = 0 and in the metallic phase chosen in the discussion of Sec.IIIC. for γ(cid:48) =0.1. In a typical experiment, the shaking amplitude would be increased from zero to a finite value. Fig.7 shows in whichorderthesystemgoesthroughthedifferentphases and dimensions upon increasing the shaking amplitude. V. POSSIBILITIES FOR EXPERIMENTAL Here, also the values of the shaking amplitude required OBSERVATION forthedimensionalcrossoversaregivenforthesamesys- tem as discussed in Sec.IIIC and γ(cid:48) =0.1γ. If the shak- Honeycombopticallatticeshaverecentlybeenrealized ing direction is perpendicular to one of the nn vectors, experimentally, although the existing setups have only thesystemwillbeinthegappedinsulatingphasebeyond beenusedtoinvestigatebosonicatoms[8,28]. Shakingof acertainvalueoftheshakingamplitude,sincetheBessel a lattice has been experimentally implemented in a one- functioncrossesthevalue0.5onlyonceandneverobtains dimensionalonebyaperiodicmodulationontheposition the value -0.5. If the shaking is parallel to one of the nn of the reflecting mirrors [14]. For a honeycomb lattice, vectors, the system will be in the metallic phase beyond the shaking could be realized by means of an acousto- a certain value of the shaking amplitude. Nevertheless, opticaldevice,asproposedforatriangularlatticeinRef. it is still possible to induce a merging of Dirac points [23]. and to open up a gap in the spectrum, since the nn hop- The magnitude of the nn hopping parameter γ in a ping parameters are renormalized such that the value of 9 (a) (b) is imposed to confine the atoms. It is described by Isotropic U = 0 Merging U = 0.1 γ 0D U = 0.2 γ 1 Alignment U = 0.3 γ V (r)= mω2 r2, (31) trap 2 trap where ω is the trapping frequency, and r is the posi- trap tion measured from the center of the trap. By apply- ing the local density approximation (LDA), one finds (c) (d) that the chemical potential evolves radially according to Isotropic Isotropic µ→µ−V (r2). Merging Metallic iso. trap 0D 0D Fig.8(a) shows the density profile for several ratios of Alignment Metallic 0D γ /γ , without nnn hopping or interactions. The case 2,3 1 withγ /γ =0,whenthesystemisintheextremelimit 2,3 1 of the band insulating phase, can be well distinguished from the other cases. Fig.8(b) shows that stronger interactions lead to a higher density away from the center of the trap. This effect becomes visible when the FIG.6: (Coloronline)Densitynasafunctionofthechemical density starts to deviate from one particle per lattice potential µ. Unless specified otherwise in the figure, the nn hopping parameters γ = γ = γ = 1, the nnn hopping site. Next-nearest-neighbor hopping leads to a higher 2,3 1 parameter γ(cid:48) = 0, the interaction strength U = 0, and the density at the edge of the cloud compared to the case inversetemperatureβ =20. (a)Effectoftherenormalization withoutnnnhopping,whichcanbeseenfromcomparing of the nn hopping parameters. (b) Effect of the interaction Figs.8(a) and (c) and from Fig.8(d). The latter shows strengthU fortheisotropiccase. (c)Effectofthennnhopping the effect of nnn hopping on the density profile for the parameter γ(cid:48) = 0.1 in the shaken lattice. (d) The metallic case where the nnn hopping gives rise to the metallic phase. For the isotropic cases, ρ = 5.2((cid:126)/mΩd)eˆ , whereas x phase for two different shaking vectors. In the first case, forthe0Dcasesρ=4.8((cid:126)/mΩd)eˆx. Thesesystemsareinthe ρ = 5.2((cid:126)/mΩd)eˆ , which gives γ ≈ γ , such that metallicphaseforγ(cid:48) =0.1,whereasforγ(cid:48) =0theyareinthe x 1 2,3 without nnn hopping, the system is in the zero-gapped semi-metallic and the insulating phase, respectively. semi-metallicphaseandtheDiracpointsarelocatedvery closetothecornersofthefirstBrillouinzone. Inthesec- ond case, ρ = 4.8((cid:126)/mΩd)eˆ , which results in γ ≈ 0, one of them will in general differ from that of the other x 2,3 such that without nnn hopping the system is in the in- two. However,whetherthesystemisactuallydriveninto sulating phase and the two energy bands are almost flat. an insulating phase or remains metallic depends on the precise value of the ratio γ(cid:48)/γ. We emphasize that, in the present paper, we only discuss weak correlations that adiabatically affect the density. However, when increasing further the onsite inter- (a) Perpendicular action, one may expect correlated phases with inhomo- 135 210 480  (nm) geneous density, even at half-filling. A detailed study of thesecorrelatedphasesisavastresearchissuethatisyet   ongoing and that is beyond the scope of the present pa- md per. Here,weonlyprovideaglimpseonhowthedensity, whichwediscussedaboveintheweak-couplinglimit,may Gapped Metallic 0D 1D evolve in view of some phases studied in the literature. Mean-fieldcalculationsindicateatransitiontoanantifer- romagneticstateaboveavalueofU/γ (cid:39)2.2[9], whereas (b) Parallel more sophisticated quantum Monte-Carlo calculations 180 278 360 412  (nm) indicate an intermediate spin-liquid phase between the semimetal and the anti-ferromagnetic phase [10]. The   spin-liquidphasemaybeviewedasaMottinsulatorwith md achargelocalizationonthelatticesites,andrecentslave- rotor calculations indicate that such spin-liquid phases FIG. 7: (Color online) Overview of the different phases as a dominate the phase diagram for γ > γ [11], which is function of the shaking amplitude. The bottom scale gives 1 2,3 theparameterrangewheretheDiracpointswouldmerge the size of the argument of the Bessel function, whereas the in the absence of interactions. The precise transition values for ρ in nm correspond to the optical lattice discussed in Sec.IIIC with γ(cid:48) = 0.1γ, Ω/2π = 6kHz, k = 830nm, and between the weakly-interacting liquid phases and these containing 40K atoms. (a) Shaking perpendicular to one of strongly-correlated Mott insulators could in principle be the nn hopping directions. (b) Shaking parallel to one of the determinedwiththehelpoftheabove-mentioneddensity nn hopping directions. measurements. A more promising technique for detecting Dirac-point Inexperiments,anoverallharmonictrappingpotential motionismomentum-resolvedRamanspectroscopy. This 10 thermore, Raman spectroscopy yields better results for (a) (b) a system with strong interactions compared to standard time-of-flight measurements [29]. Isotropic U = 0 Notice that momentum resolved Raman spectroscopy Merging U = 0.1 γ 0D U = 0.2 γ was originally proposed to be applied to a gas of ultra- Alignment U = 0.3 γ coldfermionicatomsatequilibriumandnotforashaken lattice. Wethereforediscuss,inthisfinalparagraph,why we think that this technique may also be applied to the (c) (d) present case. Naturally, as long as the frequencies of the additional lasers in the Raman-spectroscopy setup are small with respect to the shaking frequency, ω ,ω (cid:28)Ω, Isotropic Isotropic 1 2 Merging Metal iso. eventhefullsystemsatisfiesthecondition(7)fortheva- 0D 0D Alignment Metal 0D lidity of Floquet theory. As in the case of interactions, one needs, however, to avoid resonances between the different laser frequencies that could become critical [26]. The opposite limit, in which the laser frequencies and FIG. 8: (Color online) Density n as a function of the dis- that of the hyperfine transition are larger than the shak- tance from the trap’s centre r = |r|, which is expressed in ing frequency, is more delicate. However, even then, the units of the nearest-neighbor distance d. The trapping fre- shaken system remains at quasi-equilibrium as long as quency has been chosen such that the trapping potential is the intensities of the lasers used in Raman spectroscopy given by V (r) = 0.001γr2/d2. The chemical potential µ trap are weak, such that they only constitute a small per- for each case has been chosen such that the density at the turbation. The atomic dynamics probed even at high trap’s centre n is one particle per site. This corresponds to half-filling, since we consider 2 species of fermions. Unless frequencies is therefore still that of the atoms at quasi- specified otherwise in the figure, the nn hopping parameters equilibrium, with the band strucure obtained from Flo- γ = γ = γ = 1, the nnn hopping parameter γ(cid:48) = 0, quet theory. Furthermore, in the experimental studies 2,3 1 the interaction strength U = 0, and the inverse temperature by Zenesini et al. [14], time-of-flight measurements were β = 20. (a) Effect of the renormalization of the nn hop- usedtodeterminethemomentumdistributionofbosonic ping parameters. (b) Effect of the interaction strength U for atomsinashakenlattice. Apartfromthetime-scalecon- the isotropic case. (c) Effect of the nnn hopping parameter siderations, there are also some length scales that need γ(cid:48) =0.1intheshakenlattice. (d)Themetallicphase. Forthe to be taken into account. There are indeed two require- isotropic cases, ρ=5.2((cid:126)/mΩd)eˆ , whereas for the 0D cases x mentsforthecorrectsizeofthefocusofthelaserbeams. ρ = 4.8((cid:126)/mΩd)eˆ . These systems are in the metallic phase for γ(cid:48) =0.1, wherxeas for γ(cid:48) =0 they are in the semi-metallic On the one hand, it needs to be larger than the lattice spacing, such that sufficiently many atoms can be ex- and the insulating phase, respectively. cited, while on the other hand the focus of the beams should be small enough, in order to have an approximately flat trapping potential inside the focus area. In technique has been proposed as an equivalent of angle- addition,choosingthelengthofthepulsescouldpossibly resolved photoemission spectroscopy for cold atom sys- be a problem, since for shorter pulses, the excited atoms tems [29]. It has not yet been realized experimentally, will be less affected by the lattice potential, whereas for while momentum-resolved radio-frequency spectroscopy, longer pulses more atoms can be excited, leading to a which is a very similar technique, has already been im- stronger signal. plemented [30]. Notice further that another very similar technique, momentum-resolved Bragg spectroscopy, has been applied to ultracold bosonic atoms in a static optical lattice by Ernst et al. [31]. Momentum-resolved VI. CONCLUSIONS spectroscopycanallowustoindirectlyvisualizetheband structure. In momentum-resolved Raman spectroscopy, Inconclusion, wehaveinvestigatedthebandengineer- two laser pulses with frequencies ω and ω are irradi- ing of fermionic atoms in an optical honeycomb lattice 1 2 ated upon the system. If the frequency difference is in with the help of a periodic shaking of the lattice. If the resonance with a transition ω between atomic hyper- shakingfrequencyΩislargeenough,i.e. if(cid:126)Ωconstitutes hf fine states, ω −ω = ω , some atoms are excited in a thelargestenergyscaleinthesystem,theFloquettheory 1 2 hf second-orderprocesstothehigherhyperfinestate. Then, maybeappliedandthesystemisatquasi-equilibriumin with state-selective time-of-flight measurements, the dis- the sense that the atoms cannot follow the rapid motion persionoftheatomsinthenewstatearemeasured,from associated with the shaking. Depending on the direction whichthedispersionoftheoriginalatomscanbederived. of the shaking, one may render the hopping amplitudes Whentheatomsareconfinedinatrappingpotentialand in the quasi-static lattice anisotropic, due to a renormal- thelaserpulsesarefocussedonthecenterofthetrap,the ization of the nn and nnn hopping parameters by Bessel qualityoftheresultsobtainedbyRamanspectroscopyis functions that go through zero and change sign. As a comparable to those of a homogeneous system [29]. Fur- consequence, dimensional crossovers can be induced in