Simulation of nodal-line semimetal in amplitude-shaken optical lattices
Tanji Zhou, Zhongcheng Yu, Zhihan Li, Xuzong Chen, Xiaoji Zhou
SSimulation of nodal-line semimetal in amplitude-shaken optical lattices
Tanji Zhou,
1, 2
Zhongcheng Yu,
1, 2
Zhihan Li, Xuzong Chen, and Xiaoji Zhou
1, 3, ∗ State Key Laboratory of Advanced Optical Communication Systems and Networks,Department of Electronics, Peking University, Beijing 100871, China School of Physics, Peking University, Beijing 100871, China Collaborative Innovation Center of Extreme Optics,Shanxi University, Taiyuan, Shanxi 030006, China (Dated: August 28, 2020)As the research of topological semimetals develops, semimetals with nodal-line ring come intopeople’s vision as a new platform for studying electronic topology. We propose a method usingultracold atoms in a two-dimensional amplitude-shaken bipartite hexagonal optical lattice to simu-late the nodal-line semimetal, which can be achieved in the experiment by attaching one triangularoptical lattice to a hexagonal optical lattice and periodically modulating the intensity and posi-tion of the triangular lattice. By adjusting the shaking frequency and well depth difference of thehexagonal optical lattice, a transformation from Dirac semimetal to the nodal-line semimetal isobserved in our system. This transformation can be demonstrated by the Berry phase and Berrycurvature, which guides the measurement. Furthermore, the amplitude shaking and C symmetrygenerate a time-reversal-symmetry-unstable mode, and the proportion of the mode and the trivialmode of the hexagonal lattice controls the transformation. This proposal paves a way to studynodal-line semimetals in two-dimensional systems, and it contribute to the application of nodal-linesemimetals. I. INTRODUCTION
Since Hermann Weyl found a massless solution tothe Dirac equation [1] in 1929, the research of Weylfermion has attracted intensive attention. Last century,the neutrino was considered as a strong candidate forWeyl fermion, until the observation of neutrino oscilla-tion shows that neutrino is not massless [2]. Recently abreakthrough about Weyl fermion has been achieved –a gapless semimetal, named Weyl semimetal, was real-ized in photonics, condensed matter, and cold atoms [3–7]. In Weyl semimetal, low-energy excited electrons, atWeyl points, behave as Weyl fermion [8–10]. Two Weylpoints with opposite chirality can be considered as split-ting from a Dirac point by breaking time-reversal symme-try or inversion symmetry [10, 11]. During the research, anew topological semimetal was discovered, which is callednodal-line semimetal [12].As one kind of nontrivial topological semimetals, thenodal-line semimetal is different from Weyl semimetaland Dirac semimetal [13–16]. In Weyl semimetal andDirac semimetal, the bands cross at points. However,the valence and conduction bands cross and form a ring-shaped nodal line in nodal-line semimetal [17–19].In bulkstates, the unusual electromagnetic and transport re-sponse of nodal-line semimetals have attracted intenseattention [20–23]. In terms of surface states, Weylsemimetal has the Fermi arc, which is the intersectionof the 2D surface of the Brillouin zone and stretches be-tween two Weyl points. While in nodal-line semimetal,the so-called drumhead-like states are the flat bands nes-tled inside of the circle on the surface projected from ∗ [email protected] bulk nodal-line[24, 25], which is different from Fermi arcand may demonstrate the presence of interactions[23].Further, nodal-line semimetal is expected to have highfermion density at nodal-line ring due to the crossing oftwo bands, which may contribute to filling the gap be-tween fundamental physics of topological materials andpractical applications in quantum devices[26].On the other hand, the development of the artificialgauge field has paved the way for simulating topologicalmaterials by using ultracold atoms [27–35]. Furthermore,due to their high controllability, ultracold atoms in opti-cal lattices are widely used to simulate unique topologi-cal phenomena, including the measurement of the secondChern number [36], the observation of topologically pro-tected edge states [37, 38], and the chiral interaction ofWeyl semimetals [39–41]. With the development of cor-responding technologies, recently, 3D spin-orbit coupledultracold atoms in optical lattice have been used in thesimulation of nodal-line semimetal [42]. In condensedmatter systems, due to its complexity, detection of thenodal-line ring is always influenced by other bulk bands[23–25]. In comparison, optical lattices are more control-lable and can greatly remove influence from other bulkbands to directly observe nodal-line ring.In this paper, we propose a method to simulate nodal-line semimetal based on ultracold atoms in an amplitude-shaken bipartite hexagonal optical lattice. By periodi-cally changing the lattice depth of the bipartite hexag-onal optical lattice, a time-reversal-symmetry-unstablemode is imported into the system, which correspondsto nodal-line phase. As the frequency and shakingamplitude change, the energy spectrum of the opticallattice transforms from Dirac semimetal to nodal-linesemimetal. Then we discuss symmetry of the systemwhich causes the transformation and analyze the topo- a r X i v : . [ c ond - m a t . qu a n t - g a s ] A ug a b y 𝑥𝑥 FIG. 1.
Schematic diagram of 2D bipartite hexagonal lattice. (a) The lattice depth at point A and B of the bipartitehexagonal lattice are different and can be modulated periodically. (cid:126)v j are lattice vectors, (cid:126)u j are vectors connecting nearest-neighbour points ( j =1,2,3). a is the distance between nearest neighbor points. (b) We change the amplitude shaking of thelattice depth at A and B as the form of cosine function with a phase difference π , where the orange solid line and green dottedline denote points A and B, respectively. d is the difference of the average lattice depth between points A and B. ω and T denote the frequency and period of shaking, respectively. V A and V B are the lattice depth of trivial hexagonal optical lattice,and δ = V B − V A is the well difference between points A and B. ∆ is the shaken amplitude of well depth. logical characteristics of the shaking optical lattice. TheBerry curvature and Berry phase are demonstrated in thetransformation process, which provides a method to de-tect the transformation in further experiments. Finally,the influence of the well depth difference of the hexago-nal optical lattice on the nodal-line phase is specificallydiscussed, which indicates the effect of crystalline sym-metries on the nodal-line phase. Combining with the ad-vantages of ultracold atoms, it is likely to pioneer a newapproach to nodal-line semimetal in two-dimensional sys-tems. Moreover, it may serve as an important basis forfuture studies of the symmetry of nodal-line semimetals.Different from existing works [42], our system has threedistinguishing features: (1) We used amplitude-shakenoptical lattice to simulate nodal-line semimetal, which iseffective and easy to be implemented, while a great ma-jority of previous methods use spin-orbit coupling(SOC)in an optical lattice. (2) Our method can simulate atwo-dimensional nodal-line semimetal, while the existingworks mainly focus on three-dimensional systems. (3)Band structure in our system takes the form of semicon-ductor, where the conduction and valence bands touch atnodal-line ring without crossing each other, while existingworks about nodal-line semimetal are mainly semimetalform, which means the existence of an overlap betweenthe bottom of the conduction band and the top of thevalence band. This difference directly influences the dis- tribution of the density of states, which may cause un-usual electromagnetic and transport response. Here weuse the common name “nodal-line semimetal” to describeour system, while the name “nodal-line semiconductor”may be more suitable in practical terms.The remainder of this paper is organized as follows. InSec.II, we introduce a particular model which is calledthe amplitude-shaken bipartite hexagonal optical latticeand propose its feasible experimental scheme. In Sec.III,we describe the calculation to get effective Hamiltonianand derive the energy dispersion of our system. Theband structure during the transformation process fromDirac semimetal to nodal-line semimetal is shown andexplained by symmetry in Sec.IV. The discussion abouttopological properties by calculating Berry curvature andBerry phase in a certain area and the study of nodal-line phase with well depth difference are shown in Sec.V.Then we give a conclusion in Sec.VI. II. MODEL DESCRIPTION AND FEASIBLEEXPERIMENTAL SCHEMEA. Model description
In this model, the amplitude shaking is applied toa trivial honeycomb optical lattice, and the relation-ship between well depth V i and time t is shown inFig.1, where i = A, B . In these two figures, Aand B are different points which are inequivalent inhoneycomb lattice structure with a difference of welldepth d . t NN and t NNN denote the nearest-neighborhopping coefficient and next-nearest-neighbor hoppingcoefficient, respectively. And (cid:126)v j are lattice vectors,where (cid:126)v = ( − a/ , −√ a/ , (cid:126)v = (0 , √ a ) and (cid:126)v =(3 a/ , −√ a/ (cid:126)u j are vectors connecting nearest-neighbour points, where (cid:126)u = ( −√ a/ , a/ , (cid:126)u = ( a, (cid:126)u = ( −√ a/ , − a/ a = 2 √ λ/ λ is the wavelength of thelattice laser beam.The amplitude-shaken model is shown in Fig.1(b) .Thelattice depths of points A and B shake periodically inthe form of the cosine function at a fixed phase differ-ence π . In the physical image, the hexagonal lattice canbe divided into two sets of amplitude shaking triangularlattices, the amplitude of each shakes with a phase dif-ference π between them. In the figure, ω is the angularfrequency of periodically shaking, and ∆ is the shakingamplitude, which is small enough to be considered as per-turbation comparing with V A,B , where V A,B is the latticedepth of point A and B.Considering ultracold atoms, we neglect the weakatomic interactions as demonstrated in usual experiments[43], and start with a single atom model in the honey-comb lattice. The Hamiltonian ˆ H of the system canbe written as the zeroth-order Hamiltonian of the hon-eycomb lattice ˆ H adding a first-order perturbation ˆ H caused by amplitude shakingˆ H = ˆ H + ˆ H (1)According to the single-particle two-band tight-binding model, ˆ H can be written as:ˆ H = (cid:88) i V i c † i c i + (cid:88) i (cid:54) = j t (cid:48) ij c † i c j (2)where i represents each node of honeycomb lattice, V i equals to V A and V B for A and B site respectively, witha difference V B − V A = d . The operators c † i and c i denote the creation and annihilation operator at node i . t (cid:48) ij denotes the hopping coefficient between point i and j .And (cid:80) i (cid:54) = j denotes summation over all pairs of differentpoints.Next, we give the expression ˆ H when the shaking ofthe lattice depth takes the following form: V i ( t ) = V i + χ ( i ) ∆2 cos( ωt ) (3)where ω denotes the frequency of shaking. χ ( i ) at pointA equals to 1, while at point B equals to −
1. So ˆ H readsˆ H = (cid:88) i χ ( i ) ∆2 cos( ωt ) c † i c i + (cid:88) i (cid:54) = j δt ij ( t ) c † i c j (4) ba c1 c2 d FIG. 2.
The proposed experimental scheme diagram. (a) Using 3 intersecting elliptical polarized laser beams withthe enclosing angle of 120 degrees, the bipartite hexagonal isformed. Each laser beam is formed by combining two lin-early polarized laser beams with polarization directions inthe lattice plane ( V in ) and perpendicular to the lattice plane( V out ). (b) Hexagonal lattice corresponds to V out and trian-gular lattice corresponds to V in . (c) Through changing therelative phases of three lasers, two different methods of over-lying hexagonal lattice and triangular lattice are formed. Insituation (c1), the lattice depth of point B is higher than A.When the intensity of the triangular reduce to 0, we movethe triangular lattice instantaneously in order to let the lat-tice depth of point A become higher than B, as is the situationin (c2). where δt ij ( t ) is the change of hopping coefficient due toamplitude shaking.Therefore, the Hamiltonian ˆ H of the system can bewritten as: ˆ H = (cid:88) i V i ( t ) c † i c i + (cid:88) i (cid:54) = j t ij c † i c j (5)where t ij ( t ) = t (cid:48) ij + δt ij ( t ) is the hopping coefficientbetween point i and j as a function of time. In sub-sequent calculations, we only keep the nearest and thenext-nearest hopping coefficient.During the shaking process, since each pair of next-nearest sites belong to the same category of points, whichfollow the identical rule of changing in the lattice depth,the difference in well depth between them is always zero.Thus we can consider the next-nearest hopping coefficient t NNN as a constant. And the nearest hopping coefficient t NN equals to one constant t adding one term which isproportional to the lattice depth difference between pointA and B (See Appendix A), so we define t NN as t NN = t + t cos ωt (6)where t is a constant which has the unit of energy. B. Feasible experimental scheme
The above 2D amplitude shaking model can be con-structed with the following experimental scheme. Thisscheme is easy to be expanded to the 3D system byadding a laser beam perpendicular to the lattice plane[44]. As shown in Fig.2(a), we use three elliptical po-larized laser beams with an enclosing angle of 120 ◦ toeach other to form a bipartite optical hexagonal lattice,which has been demonstrated in recent works [45, 46].The total potential of optical lattice is written as V ( (cid:126)r ) = − V out (cid:88) i (cid:48) ,j (cid:48) cos (cid:104) ( (cid:126)k i (cid:48) − (cid:126)k j (cid:48) ) · (cid:126)r − ( θ i (cid:48) − θ j (cid:48) ) (cid:105) (7)+ 12 V in (cid:88) i (cid:48) ,j (cid:48) cos (cid:104) ( (cid:126)k i (cid:48) − (cid:126)k j (cid:48) ) · (cid:126)r (cid:105) where i (cid:48) , j (cid:48) = 1 , , (cid:126)k = ( √ π, − π ) /λ , (cid:126)k =( −√ π, − π ) /λ , (cid:126)k = (0 , π ) /λ are three wavevectors. λ is the wavelength of the laser beam, V out and V in de-note the components perpendicular and parallel to thelattice plane, which can be controlled respectively. Andthe three angles θ , θ , θ represent the relative phases ofthe elliptical polarization of the laser beams.It can be considered as a triangular optical latticeadding to a hexagonal one, where V out corresponds tothe triangular optical lattice and V in corresponds to thehexagonal optical lattice, as shown in Fig.2(b). Throughadjusting the phase θ i , the position of the triangular op-tical lattice can be modulated. At first, the positionof two sets of lattices is as shown in Fig.2(c1), where θ , , = (2 π/ , π/ , V out in the form of the cosine func-tion, the lattice depth at points A and B can changein the form of function | cos | . But in this way, V B isalways higher than V A . So, in order to get the shak-ing curve as shown in Fig.1(b), we need to adjust theposition of triangular optical lattice to make the lattice depth at point A and B reverse, by modulating θ , , to (2 π/ , π/ ,
0) when the intensity of triangular latticereduces to 0. Fig.2(c2) shows the position of triangularoptical lattice at this situation. The time of reversion, intheory, can be reached within several microseconds, andthe shaking period in our protocol is around 500 µs . Aslong as the time of reversion is short enough comparingwith the oscillating period, the continuity of this processremains intact. III. EFFECTIVE HAMILTONIAN ANDENERGY DISPERSION RELATION OF THEFLOQUET SYSTEM
In order to study the nodal-line semimetal, we firstcalculate the band structure of the system. In the abovesection, we derive the Hamiltonian of the system as afunction of time in Equation (5). In this section, weintroduce the method to drive the effective Hamiltonianand get the band structure of the system, which can beconsidered as averaging the ˆ H over time.To start with, we act an unitary transformation on ˆ H .Define unitary operator ˆ U :ˆ U = exp (cid:32) i (cid:126) (cid:88) i (cid:90) t dτ · V i ( τ ) c † i c i (cid:33) (8)Through this unitary transform, we get e -index form ˆ H (cid:48) (See Appendix B):ˆ H (cid:48) = ˆ U (cid:18) ˆ H ( t ) − i (cid:126) ∂∂t (cid:19) ˆ U † − (cid:18) − i (cid:126) ∂∂t (cid:19) (9)= − d (cid:88) i ( a † i a i − b † i b i ) − (cid:88) i (cid:54) = j t ij e i z ij sin ωt c † i c j where z ij = ∆2 ω (cid:126) [ χ ( i ) − χ ( j )] characterizes the responsesof amplitude shaking at different sites.Next, the effective Hamiltonian can be obtained byusing high-frequency expansion. We only keep up to first-order terms, and the effective Hamiltonian takes the form(More details see Appendix C): H eff = H (0) eff + 1 (cid:126) ω H (1) eff = H f + ∞ (cid:88) n =1 [ H n , H − n ] n (cid:126) ω (10)where H f is the zeroth-order term, and H n means the e inωt term Fourier expansion coefficients of Hamiltonianin Eq.(9). Then we get the kernel of effective Hamilto-nian: H eff = H eff, ˆ I + H eff,x ˆ σ x + H eff,y ˆ σ y + H eff,z ˆ σ z (11)where H eff, = − t NNN (cid:88) j =1 cos( (cid:126)k · (cid:126)v j ) H eff,x = − t J ( β ) (cid:88) j =1 cos( (cid:126)k · (cid:126)u j ) H eff,y = − t J ( β ) (cid:88) j =1 sin( (cid:126)k · (cid:126)u j ) H eff,z = 8 t t ∆ J ( β ) (cid:88) j =1 cos( (cid:126)k · (cid:126)v j ) + 32 − δ ˆ σ x , ˆ σ y , ˆ σ z are Pauli matrices, and (cid:126)k is the wave vector ofatomic state function. J n means the n th order Besselfunction of the first kind. Here we define shaking factor β and δ to describe the shaking: β ≡ ∆ (cid:126) ω (12) δ ≡ d/ t t J ( β ) / ∆ (13) β and δ are the main parameters affecting the effec-tive Hamiltonian since they include shaking frequency ω ,shaking amplitude ∆, and the difference of well depth d . Through solving eigenvalues of the kernel of effec-tive Hamiltonian (Eq.(11)), the energy dispersion rela-tion E ( (cid:126)k ) of lowest two bands can be gotten: E ± ( (cid:126)k ) = H eff, ± (cid:113) H eff,x + H eff,y + H eff,z (14)Here we set δ = 3 / H eff,z , which makes this system betterreflect the effects of shaking. Then, the energy disperionrelation in Eq.(14) mainly depends on parameter β =∆ / (cid:126) ω , and hopping coefficients t and t also dependon ∆, (see Appendix A). And further discussion aboutdifferent values of δ will be shown in Section VI. IV. THE TRANSFORMATION FROM DIRACSEMIMETAL TO NODAL-LINE SEMIMETAL
The value of J n ( β ) will change as one of the shakingfactor β changes, as shown in Fig.3(a), where the differentvalue of J and J corresponding to β . With J ( β ) = 0, H eff,z is zero, but H eff,y have a finite quantity. E + ( (cid:126)k )and E − ( (cid:126)k ) equals to each other at six Dirac points, whichperforms as Dirac semimetal, as shown in Fig.3(b1). Itis worth mentioning that although shaking breaks the in-version symmetry, the gap between the upper and lowerband is still closed at Dirac point. Fig.3(b2) is the sec-tional view of as Fig.3(b1). In the figure, the blue linemeans the upper band and the dashed red line representsthe lower band, which touches the former at the edge ofthe Brillouin zone. a FIG. 3.
Band structure of the system. (a) shows thevariation curve of J ( β ) and J ( β ), when β changes from oneof zeros of zeroth-order Bessel function J ( β = 2 . J ( β = 3 . β = 2.43, 2.50, 2.60, 2.90, whichis to show the main character of the transformation processlater in Fig.4. (b) The band structure of Dirac semimetal.(b1) When J ( β ) = 0, the optical lattice performs as Diracsemimetal. (b2) is the sectional view of over-momentum zerocorresponding to (b1), where k y has been set to zero and wefocus on the change of energy with k x . (c) The band structureof nodal-line semimetal. (c1) When J ( β ) = 0, the opticallattice performs as nodal-line semimetal. (c2) is the sectionalview of over-momentum zero corresponding to (c1). In thissituation, conduction band and valence band touch at ring k x + k y = ( π ) /a . In the figure, k = π √ a is the distancefrom center of first Brillouin zone to the Dirac point. In another situation for J ( β ) = 0, H eff,x and H eff,y is zero, but H eff,z is nonzero. E + ( (cid:126)k ) and E − ( (cid:126)k ) equalsto each other at a ring k x + k y = ( π ) /a . This energydispersion relation performs as nodal-line semimetal, asshown in Fig.3 (c1). The upper and lower band touchinside the first Brillouin zone, which differs from the casein Dirac semimetal where the touchpoint is at the verticesof the first Brillouin zone. The most noteworthy featureof nodal-line semimetal is that nodes of two bands arecontinuous and form a so-called nodal-line ring, whereinour model locates at k x + k y = ( π ) /a . Fig.3(c2) issectional view of Fig.3(c1). The two bands touch at thenodal-line ring, and at the edge of the Brillouin zone, thegap between two bands is open.When β changes from J ( β ) = 0 to J ( β ) = 0, thesystem transforms from Dirac semimetal to nodal-linesemimetal. Fig.4 shows the band structure for the othervalue between J ( β ) = 0 and J ( β ) = 0 correspondingto the four vertical dot dashed red lines and two colorfulcurves in Fig.3(a). When the energy spectrum graduallyvaries from nodal-line semimetal to Dirac semimetal, thetouching ring gradually opens as the Dirac points grad-ually close, as shown in Fig.4(a-d). In the experiment, β can be changed continuously by fixing the amplitude ofthe lattice depth shaking and adjusting the rotating fre-quency continuously to observe the transformation fromDirac semimetal to nodal-line semimetal.In Fig.3 and Fig.4, the parameters are chosen from theexperimental parameters of Rb atom: V = 10 . E r , t = 0 . E r , t = 0 . E r , t NNN = 7 . × − E r ,and ∆ = 0 . E r , where E r means atomic recoil energyof electron in our honeycomb lattice.The transformation from Dirac semimetal to nodal-line semimetal can be explained by symmetry as follows.The impact of amplitude shaking is reflected in Hamil-tonian, so we first study the Hamiltonian ˆ H (cid:48) . ThroughJacobi-Anger expansion (ignoring terms unrelated to β ),the Hamiltonian ˆ H (cid:48) can be written as the following form:ˆ H (cid:48) = − (cid:88) (cid:104) i,j (cid:105) and i (cid:54) = j t ij e iz ij sin ωt c † i c j (15)= − (cid:88) (cid:104) ij (cid:105) J ( β ) t NN c † i c j − ( e iωt + e − iωt ) · (cid:88) (cid:104) ij (cid:105) J ( β ) t NN c † i c j − o( t NN c † i c j )Where (cid:80) (cid:104) ij (cid:105) means the sum of nearest-neighbor latticenodes, and o ( t NN c † i c j ) represents the higher order in-finitesimal term considering the pairs of sites with far-ther distance between them than the next-nearest neigh-bor term comparing to t NN c † i c j . So the system mainly iscomposed of two modes, the weights of which are J ( β )and J ( β ).As shown in Fig.5, the expanding term of Hamilto-nian ˆ H (cid:48) can be divided into two parts. The first termcorresponds to Mode I, which is a bipartite hexagonallattice Mode. When J ( β ) = 0, only mode I exists in thesystem, so the system performs like an ordinary Diracsemimetal. In the situation, Dirac points are closed dueto the spatial inversion symmetric potential energy inMode I.The second term of Hamiltonian ˆ H (cid:48) corresponds toMode II, which consists of two symmetrical rotating lat-tices. The Mode II is the same as Mode I, both of whichhave time-reversal symmetry, but different from Mode I,the time inversion symmetry of Mode II is unstable. Inour system, the phase difference is π , so the differenceof well depth can be written as a cosine function, and FIG. 4.
The transformation of band structure as β changes. The band structure for the intermediate states ofDirac and node-line semimetal, at β = 2 . , . , . , . the amplitude of two sub-modes are equivalent to eachother. If we change the phase difference, the form of welldepth difference between point A and B change conse-quently, which causes the two sub-modes of Mode II tobe asymmetric and leads to the broken of time-reversalsymmetry. When J ( β ) = 0, Mode I will vanish whileMode II remains existing in the system, so its valenceand conducting band touch at a ring, performing as anodal-line semimetal. When β is between the zeros of J to J , the system is in the intermediate state of twomodes. V. BERRY CURVATURE AND BERRY PHASE
Considering that this is a two-dimensional system, inthis section, we calculate the Berry curvature and Berryphase rather than winding number to provide a measur-able quantity in the experiment during the transforma-tion from Dirac semimetal to nodal-line semimetal.From Eq.(11), we can define (cid:126)h as: (cid:126)h = H eff,x · ˆ e x + H eff,y · ˆ e y + H eff,z · ˆ e z (16)where ˆ e x , ˆ e y , ˆ e z are the basis vectors at x , y , z directionof momentum space. The role that (cid:126)h plays is similar to amagnetic field coupling with the vector of Pauli matrices(ˆ σ x , ˆ σ y , ˆ σ z ). By solving the eigenequations of effectiveHamiltonian, we get ψ + = 1 (cid:112) h ( h + h ) (cid:18) h + hh − i h (cid:19) (17) + ModeⅠModeⅡ 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑖𝑖 − � <𝑖𝑖𝑖𝑖> 𝒥𝒥 𝛽𝛽 𝑡𝑡 𝑁𝑁𝑁𝑁 𝑐𝑐 𝑖𝑖† 𝑐𝑐 𝑖𝑖 𝑒𝑒 𝑖𝑖𝑖𝑖𝑖𝑖 − ( 𝑒𝑒 𝑖𝑖𝑖𝑖𝑖𝑖 + 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑖𝑖 ) � <𝑖𝑖𝑖𝑖> 𝒥𝒥 𝛽𝛽 𝑡𝑡 𝑁𝑁𝑁𝑁 𝑐𝑐 𝑖𝑖† 𝑐𝑐 𝑖𝑖 FIG. 5.
Two shaking modes of the system.
The resultof Fourier transform of Hamiltonian ˆ H (cid:48) reflects two shakingmodes of the system, where Mode I is a bipartite hexago-nal lattice mode and Mode II is a nontrivial time-reversal-symmetry-unstable mode, which corresponds to nodal-linesemimetal. ψ − = 1 (cid:112) h ( h − h ) (cid:18) h − hh − i h (cid:19) (18)where (cid:126)h = h ˆ e x + h ˆ e y + h ˆ e z , h = (cid:112) h + h + h .Two components A i ( i = 1 ,
2) of Berry connection (cid:126) A of the lowest band can be calculated from ψ − as A i ( (cid:126)k ) = i (cid:104) ψ − | ∂ k i | ψ − (cid:105) = − h ∂ k i h − h ∂ k i h h ( h − h ) (19)Then the Berry curvature of the lowest band can becalculated as (cid:126) Ω( (cid:126)k ) = −∇ × (cid:126) A ( (cid:126)k ) (20)The z component of the Berry curvature isΩ = − ∂ A ∂k + ∂ A ∂k (21)Hence the Berry phase can be calculated from Berrycurvature: γ = (cid:90) S d (cid:126)S · (cid:126) Ω( (cid:126)k ) = (cid:90) S d S Ω ( (cid:126)k ) (22) a bc d 𝛽𝛽 = 2.43 𝛽𝛽 = 2.50 𝛽𝛽 = 2.60 𝛽𝛽 = 2.90 FIG. 6.
The change of Berry curvature during thetransformation process. (a-d) Berry curvature in thereciprocal space in the process of the system transform-ing from nodal-line semimetal to Dirac semimetal. β =2 . , . , . , .
90, which is correspond to the four lines inFig.3(a), respectively. The dashed equilateral triangle denotesthe integral region S of the Berry phase, which goes aroundone vertex of the first Brillouin Zone. Here S denotes the area in the reciprocal space. Fig.6shows the Berry curvature in the transformation fromnodal-line semimetal to Dirac semimetal. In the figure,the values of β are selected according to the intersec-tions of four vertical dot dashed red lines and two col-orful curves in Fig.3 (a). When β is near 2.4048 (oneof the zeros of J ), the system is nearly a nodal-linesemimetal, and Berry curvature of which is nontrivialat the nodal-line ring. The distribution of Berry curva-ture can naturally divide into six parts, where the nu-merical value of adjacent parts is of opposite sign, andthe sum of six parts equals to zero, which is protectedby time-reversal symmetry. With the system transform-ing to Dirac semimetal, while β increases from 2.4048to 3.8317 (the nearest zero of J ), those positive andnegative parts of Berry curvature respectively becomecloser together. Finally, those parts converge on six Diracpoints and keep shrinking to approach the distributionin the case of Dirac semimetal. Because the shrinking ofBerry curvature into one point makes it difficult to dis-tinguish in the figure, here we only show the result up to β = 2 .
90 in Fig.6.The red dashed equilateral triangle which connectsthe center of one site and two adjacent sites in Fig.6shows the path around one Dirac point, and we studythe change of the Berry phase along the path in thetransformation. Generally the path to calculate berry
FIG. 7.
Berry phase around a Dirac point during thetransformation.
Here we set δ = 3 /
2. In the figure, thered circles are simulated data points, and the red line is fittedcurve. The two black dotted lines mark the start and endBerry phase during the transformation. When the system is aDirac semimetal, the Berry phase is equal to π . And when thesystem transforms to nodal-line semimetal, the Berry phasechanges continuously to π/ phase should avoid the passing Berry-curvature singular-ities [47], but in our system, Berry-curvature singularitieswill change from Dirac points to nodal-line ring duringthe transformation. It is impossible to choose a pathalong a Dirac point disjointing with any singularities ina transformation. Hence we choose the equilateral tri-angle path which is along a Dirac point and passes highsymmetry point of Brillouin zone. Fig.7 shows the re-sult of the Berry phase as β increases from 2.408 to 3.70(between zeros of J J π/
3, which corresponds to the situation inFig.6(a). When the system is near Dirac semimetal, ourFloquet system gets the same result as in bipartite hexag-onal lattice, and the Berry phase around Dirac pointcomes close to π . In the transformation procession, asis shown in Fig.6(a)-(d), the Berry phase around Diracpoint increases continuously, with the Berry curvatureshrinking into one point.Above, we focus δ = 1 . β to observethe transformation from Dirac semimetal to Nodal-linesemimetal. Next, we will set β = 2 . δ changes.The band structure of the nodal-line phase with dif-ferent δ is shown in Fig.8(a). In the figure, throughchanging the well difference δ , different band structureof nodal-line phase is gotten. When δ <
0, the upperand lower bands don’t intersect with each other; when δ ∈ (0 , . C symmetry; b1 𝜹𝜹 = −𝟎𝟎 . b2 b3b4 b5 b6 𝜹𝜹 = .
𝟑𝟑 𝜹𝜹 = . = .
𝟕𝟕 𝜹𝜹 = .
𝟓𝟓 𝜹𝜹 = . FIG. 8.
The band structure of nodal-line phase with δ changes. (a1-a6) The front view of the band struc-ture for the nodal-line phase at β = 2 . δ = − . , . , . , . , . , . β = 2 . δ = − . , . , . , . , . , . . when δ ∈ (0 . , . δ ∈ (1 . , . δ increases. When δ = 4 .
5, thereis a single node at the center of the first Brillouin zone,and two bands will be gapped for δ > . δ and symmetry analysis arestudied for fixed β = 2 . VI. CONCLUSIONS
In summary, we propose a feasible scheme to sim-ulate nodal-line semimetal with ultracold atoms in anamplitude-shaken optical lattice. We derive the effec-tive Hamiltonian of the Floquet system, and by calcu-lating the band structure, the transformation from Diracsemimetal to nodal-line semimetal is observed. Whenthe shaking factor β , which is determined by the shak-ing amplitude and frequency, is at zeros of the first-orderBessel function of the first kind, the band structure per-forms as Dirac semimetal, and when shaking factor is atthe zeros of the zeroth-order Bessel function, the bandstructure performs as a nodal-line semimetal. ThroughFourier transform, we divide the shaking into two modes,which explain the above transformation. The change ofBerry curvature and Berry phase during the transforma-tion process shows the topological characteristics of oursystem. However, if we set the phase difference of am-plitude oscillation of sites A and B to not be π exactly,the time-reversal symmetry of the system will be bro-ken, which may lead to the research of new topologicalsemimetals. VII. ACKNOWLEDGEMENT
We thank Xiaopeng Li, Wenjun Zhang, Yuan Zhanfor helpful discussion. This work is supported bythe National Basic Research Program of China (GrantNo. 2016YFA0301501) and the National Natural Sci-ence Foundation of China (Grants No. 61727819, No.11934002, No. 91736208, and No. 11920101004).
Appendix A: hopping coefficient
We estimate the nearest-neighbor hopping coefficientin our Floquet system, taken Rb atom for instance.For a hexagonal optical lattice with different well depthat points A and B, the overlapping integral of Wannierfunction helps us to calculate the nearest-neighbor hop-ping coefficient t NN , which renders result consistent withEquation (6).Fig.9 shows the change of nearest-neighbor hoppingcoefficient as the difference of well depth ( V A − V B ) /V A changes. As a result, the nearest-neighbor hopping coef-ficient t NN equals one constant t adding one term whichis nearly proportional to the lattice depth difference be-tween points A and B in a large range (Correlation coef-ficient is 0 . V A is less than 20%, and the difference of well depth ( V A − V B ) /V A isabout 35%. The linearity will be better if the processhappens in a smaller range. So the static process can beused to estimate our dynamic perturbation process. Inour system, the lattice depth difference between pointsA and B changes over time as a cosine function, whichleads to the conclusion that the nearest-neighbor hop-ping coefficient also changes as a cosine function, so weuse t NN = t + t cos( ωt ) to estimate it in the calculation.Using the fitting result, we can calculate that t = 0 . V A E r (A1) t = 0 . E r − . V A E r (A2)where E r means the atomic recoil energy of electron inour honeycomb lattice. By substituting ∆ = 0 . E r and V A = 10 E r , we get t = 0 . E r , t = 0 . E r (A3)As for the next-nearest-neighbor hopping coefficient,because the lattice depth of the next neighbor point isalways 0, we use a constant to replace the next neighborhopping coefficient t NNN = 0 . E r (A4)In addition, according to Eq.(10), the high-frequencyexpansion requires shaking frequency ω to be far greaterthan H f / (cid:126) . In our system, it requires that difference ofthe lattice depth ∆ is larger than 0.15 V A , which meetsthe linear range. Appendix B: Unitary transformation of Hamiltonian
The unitary operator ˆ U can be obtained by substitut-ing Equation (3) into Equation (8):ˆ U ( t )= exp (cid:34) i (cid:126) (cid:88) i (cid:90) t dτ · χ ( i ) ∆2 cos( ωτ ) c † i c i (cid:35) = exp (cid:34) i (cid:126) (cid:88) i χ ( i ) ∆2 ω sin( ωτ ) c † i c i (cid:35) (B1)Then we getˆ U (cid:18) − i (cid:126) ∂∂t (cid:19) ˆ U † = − i (cid:126) ∂∂t − (cid:88) i ∆2 cos( ωt ) χ ( i ) c † i c i (B2)Define ˆ U (cid:48) k ( t ) asˆ U (cid:48) k ( t ) = exp (cid:20)(cid:18) i (cid:126) ∆2 ω sin ωtχ ( k ) (cid:19) c † k c k (cid:21) (B3)0 FIG. 9.
The theoretical calculation of the nearestneighbor hopping coefficient J as difference of welldepth changes. The nearest neighbor hopping coefficient ina hexagonal optical lattice nearly changes linearly as the welldepth difference between the lattice depth of point A and Bchanges. The correlation coefficient r is 0.9995, which provesthe good linearity. In the figure, the red circle is the calcu-lated data and the line is the fitting line. The well depthat point A is about 10 E r , where E r means the atomic recoilenergy in our honeycomb lattice. According to Baker-Campbell-Hausdorff formula, e ˆ A ˆ Be − ˆ A = ˆ B + [ ˆ A, ˆ B ] + 12! [ ˆ A, [ ˆ A, ˆ B ]] + . . . (B4)We can getˆ U ( − t ij c † i c j ) ˆ U † = (cid:88) k ˆ U (cid:48) k ( − t ij c † i c j ) (cid:88) k ˆ U (cid:48)† k = − t ij e i (cid:126) ∆2 ω sin ωt [ χ ( i ) − χ ( j )] c † i c j (B5)The transformed Hamiltonian was finally obtainedˆ H (cid:48) = ˆ U (cid:18) ˆ H ( t ) − i (cid:126) ∂∂t (cid:19) ˆ U † − (cid:18) − i (cid:126) ∂∂t (cid:19) (B6)= (cid:88) i V i ( t ) c † i c j − (cid:88) i (cid:54) = j t ij e iz ij sin ωt c † i c j − i (cid:126) ∂∂t − (cid:88) i ∆2 cos( ωt ) χ ( i ) c † i c i + i (cid:126) ∂∂t = − d (cid:88) i ( a † i a i − b † i b i ) − (cid:88) i (cid:54) = j t ij e i z ij sin ωt c † i c j where the first term in Equation (2) was neglected, sincethe zero point of energy has been set at V , and z ij ≡ ∆2 ω (cid:126) [ χ ( i ) − χ ( j )]. Appendix C: Fourier expansion and derivation ofeffective Hamiltonian
By a Jacobi-Anger expansion, e i z sin θ = (cid:80) + ∞ n = −∞ J n ( z ) e i nθ , we expand ˆ H (cid:48) in Eq.(9), tak- ing both nearest-neighbor and next-nearest-neighborhopping into consideration, thenˆ H (cid:48) ( t )= − d (cid:88) i ( a † i a i − b † i b i ) − + ∞ (cid:88) n = −∞ e i nωt (cid:88) (cid:104) ij (cid:105) ( t J n ( z ij )+ t J n − ( z ij ) + J n +1 ( z ij )] (cid:19) c † i c j + (cid:88) (cid:104)(cid:104) ij (cid:105)(cid:105) J n (0) t NNN c † i c j (C1)where (cid:80) (cid:104) ij (cid:105) denotes the summation over the nearest-neighbor nodes, and (cid:80) (cid:104)(cid:104) ij (cid:105)(cid:105) denotes the summation overthe next-nearest-neighbor nodes. Then we can obtain theFourier expansion coefficients of Hamiltonian for n > H n = − (cid:88) (cid:104) ij (cid:105) (cid:18) t + t nz ij (cid:19) J n ( z ij ) c † i c j − (cid:88) (cid:104)(cid:104) ij (cid:105)(cid:105) J n (0) t NNN c † i c j (C2)From the Floquet theory, we can derive the formula tocalculate effective Hamiltonian[51, 52] H eff = H eff + 1 (cid:126) ω H eff = H f + ∞ (cid:88) n =1 [ H n , H − n ] n (cid:126) ω (C3)where higher-order terms which consider the pairs of siteswith the farther distance between them than the next-nearest neighbor term have been neglected due to high-frequency approximation. The calculation of substitutingEq.(C2) into Eq.(C3) is shown below.For n = 0 term, H f = − (cid:88) (cid:104) ij (cid:105) t J n ( z ij ) c † i c j − (cid:88) (cid:104)(cid:104) ij (cid:105)(cid:105) t NNN c † i c j − d (cid:88) i ( a † i a i − b † i b i )= H , + H , + H , (C4)where H , denotes the summation over the nearest-neighbor nodes, and H , denotes the summation overthe next-nearest-neighbor nodes. For H , , because J ( − z ij ) = J ( z ij ), we can get H , = − t J n ( β ) (cid:88) i ( a † (cid:126)r i b (cid:126)r i + (cid:126)u + a † (cid:126)r i b (cid:126)r i + (cid:126)u + a † (cid:126)r i b (cid:126)r i + (cid:126)u ) + h.c. (C5)where β ≡ ∆ (cid:126) ω , a i , b j denote the annihilation operators oflattice site A i , B j . The first term describes the hopping1from B j to A i , while the second term is the Hermitianconjugate of the former one, describes the hopping from A i to B j . The (cid:126)r i is the lattice vector for bipartite hexag-onal lattice, and the (cid:126)u i is the nearest-neighbor vector.Considering the creation and annihilation operatorsas periodic functions in real space, we can take theFourier transformation to get the corresponding cre-ation and annihilation operators in momentum space a i = 1 √ N (cid:80) k a k e − i (cid:126)k · (cid:126)r and b i = 1 √ N (cid:80) k b k e − i (cid:126)k · (cid:126)r , where N is the number of lattice sites. By substituting theseequations into the effective Hamiltonian in Eq.(C5), us-ing δ (cid:126)k (cid:126)k (cid:48) = N (cid:80) i e i ( (cid:126)k − (cid:126)k (cid:48) ) · (cid:126)r i , we can obtain that H , = − t J n ( β ) (cid:88) k a † k b k ( e − i (cid:126)k · (cid:126)u + e − i (cid:126)k · (cid:126)u + e − i (cid:126)k · (cid:126)u ) + h.c. = (cid:88) k (cid:0) a † k b † k (cid:1) (cid:18) H ( k ) H (cid:63) ( k ) 0 (cid:19) (cid:18) a k b k (cid:19) (C6)where H ( k ) = − t J ( β ) (cid:80) j =1 e − i (cid:126)k · (cid:126)u j For H , , denote the lattice constant of bipartitehexagonal lattice as (cid:126)v j , as shown in Fig.1(a). Fouriertransform and direct calculation gives that H , = − (cid:88) (cid:104)(cid:104) ij (cid:105)(cid:105) t NNN c † i c j = − t NNN (cid:88) j =1 (cid:88) k ( a † k a k + b † k b k ) cos( (cid:126)k · (cid:126)v j )= (cid:88) k (cid:0) a † k b † k (cid:1) (cid:18) H ( k ) 00 H ( k ) (cid:19) (cid:18) a k b k (cid:19) (C7)where H = − t NNN (cid:80) j =1 cos( (cid:126)k · (cid:126)v j ).As for H , , H = − d (cid:88) i ( a † i a i − b † i b i )= − d (cid:88) k ( a † k a k − b † k b k )= (cid:88) k (cid:0) a † k b † k (cid:1) (cid:18) H −H (cid:19) (cid:18) a k b k (cid:19) (C8)where H = − d/ n > H n = − (cid:88) (cid:104) ij (cid:105) (cid:18) t J n ( z ij ) + t J n − ( z ij ) + J n +1 ( z ij )] (cid:19) c † i c j (C9) According to Eq.(C3), H (1) eff = ∞ (cid:88) n =1 [ H n , H − n ] n (C10)= 4 t t β J ( β ) · (cid:88) i [ a † i a i − b † i b i ] + (cid:88) (cid:104)(cid:104) ij (cid:105)(cid:105) [ a † i a j − b † i b j ] = (cid:88) k (cid:0) a † k b † k (cid:1) (cid:18) H −H (cid:19) (cid:18) a k b k (cid:19) which uses the commutation relationship [ c † k c l , c † p c q ] = δ lp c † k c q − δ kq c † p c l . And H = 8 t t J ( β ) β (cid:88) j =1 cos( (cid:126)k · (cid:126)v j ) + 32 (C11)Combining the results of Eq.(C6), (C7) and (C10), wecan draw the conclusion that H eff = H (0) eff + 1 (cid:126) ω H (1) eff = (cid:88) k (cid:0) a † k b † k (cid:1) (cid:18) H eff, ( (cid:126)k ) H eff, ( (cid:126)k ) H eff, ( (cid:126)k ) H eff, ( (cid:126)k ) (cid:19) (cid:18) a k b k (cid:19) (C12)where H eff, = (cid:20) − t NNN + 8 t t (cid:126) ω J ( β ) β (cid:21) (cid:88) j =1 cos( (cid:126)k · (cid:126)v j )+ 12 t t (cid:126) ω J ( β ) β (C13) H eff, = (cid:20) − t NNN − t t (cid:126) ω J ( β ) β (cid:21) (cid:88) j =1 cos( (cid:126)k · (cid:126)v j ) − t t (cid:126) ω J ( β ) β (C14)and H eff, = − t J ( β ) (cid:88) j =1 e − i (cid:126)k · (cid:126)u j − d H eff, = − t J ( β ) (cid:88) j =1 e i (cid:126)k · (cid:126)u j + d × H eff , which can bespaned by 2-ranked identity matrix and three 2D Paulimatrices according to linear algebraic theory. The finalresult is exactly Eq.(11) in the article.2 FIG. 10.
The change of Berry phase with shaking fac-tor δ . Here we set β = 2 . δ < . γ is alwaysequal to 2 π . When δ > . γ continuously decreases from 2 π to 0. Appendix D: The Berry phase and symmetryanalysis of nodal-line phase
The change of the Berry phase around the closed pathin Fig.6 with the change of δ is calculated and studiedfor fixed β = 2 . δ < .
5, the Berry phase γ isalways equal to 2 π , and a mutation of the Berry phase oc-curs at δ = 0 .
5. The nodal-line in the first Brillouin zoneat δ = 0 . δ > .
5, the area of nodal-line lies all inside thefirst Brillouin zone, so the nodal-line is a complete closedunrounded loop. However, when δ < .
5, the solution( k x , k y ) satisfying E + − E − = 0 will beyond the scope ofthe first Brillouin zone, so the nodal-line becomes incom-plete. Hence δ = 0 . δ is connected to the symme-try of the lattice. When δ = 0 point A is equivalent topoint B, the lattice meets six-fold rotational symmetry.When lattice meets six-fold rotational symmetry, nodallines degenerate into points at the six vertices of the Bril-louin zone. When δ >
0, the six-fold rotational symmetrybreaks and nodal lines appear. If we fixed β at 2.4048,the nodal-line phase is present for a range of parameter δ ,which is δ ∈ (0 , . δ > δ < [1] H. Weyl, Zeitschrift f¨ur Physik , 330 (1929).[2] Q. R. Ahmad, R. C. Allen, T. C. Andersen, andet al (SNO Collaboration), Phys. Rev. Lett. , 071301(2001).[3] L. Lu, L. Fu, J. D. Joannopoulos, and M. Soljacic, NaturePhoton. , 294 (2013).[4] M. M. Vazifeh and M. Franz, Phys. Rev. Lett. ,027201 (2013).[5] S.-M. Huang, S.-Y. Xu, I. Belopolski, C.-C. Lee,G. Chang, B. Wang, N. Alidoust, G. Bian, M. Neupane,C. Zhang, et al., Nat. Commun. , 7373 (2015).[6] C. He, Phys. Lett. A , 440 (2018).[7] J. Noh, S. Huang, D. Leykam, Y. . D. Chong, K. P. Chen,and M. Rechtsman, Nature Phys. , 611 (2017).[8] L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos,and M. Soljaˇci´c, Science , 622 (2015).[9] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane,G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C.Lee, et al., Science , 613 (2015).[10] B. Yan and C. Felser, Annu. Rev. Conden. Ma. P. , 337(2017).[11] Y. Chen, Y.-M. Lu, and H.-Y. Kee, Nat. Commun. ,6593 (2015).[12] A. A. Burkov, M. D. Hook, and L. Balents, Phys. Rev.B , 235126 (2011). [13] R. Yu, H. Weng, Z. Fang, X. Dai, and X. Hu, Phys. Rev.Lett. , 036807 (2015).[14] M. Hirayama, R. Okugawa, T. Miyake, and S. Murakami,Nat. Commun. , 14022 (2017).[15] Y.-H. Chan, C.-K. Chiu, M. Y. Chou, and A. P. Schny-der, Phys. Rev. B , 205132 (2016).[16] L. Li, C. H. Lee, and J. Gong, Phys. Rev. Lett. ,036401 (2018).[17] K. Mullen, B. Uchoa, and D. T. Glatzhofer, Phys. Rev.Lett. , 026403 (2015).[18] Y. Kim, B. J. Wieder, C. L. Kane, and A. M. Rappe,Phys. Rev. Lett. , 036806 (2015).[19] Q. Xu, R. Yu, Z. Fang, X. Dai, and H. Weng, Phys. Rev.B , 045136 (2017).[20] M. Koshino and T. Ando, Phys. Rev. B , 195431(2010).[21] M. Koshino and I. F. Hizbullah, Phys. Rev. B , 045201(2016).[22] G. P. Mikitik and Y. V. Sharlai, Phys. Rev. B , 195123(2016).[23] J. Behrends, J.-W. Rhim, S. Liu, A. G. Grushin, andJ. H. Bardarson, Phys. Rev. B , 245101 (2017).[24] H. Weng, Y. Liang, Q. Xu, R. Yu, Z. Fang, X. Dai, andY. Kawazoe, Phys. Rev. B , 045108 (2015). [25] G. Bian, T.-R. Chang, H. Zheng, S. Velury, S.-Y. Xu,T. Neupert, C.-K. Chiu, S.-M. Huang, D. S. Sanchez,I. Belopolski, et al., Phys. Rev. B , 121113 (2016).[26] J. Hu, Z. Tang, J. Liu, X. Liu, Y. Zhu, D. Graf, K. Myhro,S. Tran, C. N. Lau, J. Wei, et al., Phys. Rev. Lett. ,016602 (2016).[27] J. Dalibard, F. Gerbier, G. Juzeli¯unas, and P. ¨Ohberg,Rev. Mod. Phys. , 1523 (2011).[28] P. Hauke, O. Tieleman, A. Celi, C. ¨Olschl¨ager, J. Si-monet, J. Struck, M. Weinberg, P. Windpassinger,K. Sengstock, M. Lewenstein, et al., Phys. Rev. Lett. , 145301 (2012).[29] X.-J. Liu, K. T. Law, T. K. Ng, and P. A. Lee, Phys.Rev. Lett. , 120402 (2013).[30] X.-J. Liu, K. T. Law, and T. K. Ng, Phys. Rev. Lett. , 086401 (2014).[31] P. Hauke, M. Lewenstein, and A. Eckardt, Phys. Rev.Lett. , 045303 (2014).[32] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat,T. Uehlinger, D. Greif, and T. Esslinger, Nature ,237 (2014).[33] L. Niu, X. Guo, Y. Zhan, X. Chen, W. M. Liu, andX. Zhou, Applied Physics Letters , 144103 (2018).[34] B. Song, L. Zhang, C. He, T. F. J. Poon, E. Hajiyev,S. Zhang, X.-J. Liu, and G.-B. Jo, Sci. Adv. (2018).[35] X. Guo, W. Zhang, Z. Li, H. Shui, X. Chen, and X. Zhou,Opt. Express , 27786 (2019).[36] M. Lohse, C. Schweizer, H. M. Price, O. Zilberberg, andI. Bloch, Nature , 55 (2018).[37] B. K. Stuhl, H.-I. Lu, L. M. Aycock, D. Genkina, andI. B. Spielman, Science , 1514 (2015).[38] M. Leder, C. Grossert, L. Sitta, M. Genske, A. Rosch,and M. Weitz, Nat. Commun. , 13112 (2016).[39] T. Dubˇcek, C. J. Kennedy, L. Lu, W. Ketterle,M. Soljaˇci´c, and H. Buljan, Phys. Rev. Lett. , 225301 (2015).[40] D.-W. Zhang and S. Cao, Laser. Phys. Lett. , 065201(2016).[41] X. Kong, J. He, Y. Liang, and S.-P. Kou, Phys. Rev. A , 033629 (2017).[42] B. Song, C. He, S. Niu, L. Zhang, Z. Ren, X.-J. Liu, andG.-B. Jo, Nature Phys. , 911 (2019).[43] E. Arimondo, D. Ciampini, A. Eckardt, M. Holthaus,and O. Morsch, in Advances in Atomic, Molecular, andOptical Physics , edited by P. Berman, E. Arimondo, andC. Lin (Academic Press, 2012), vol. 61 of
Advances InAtomic, Molecular, and Optical Physics , pp. 515 – 547.[44] S. Jin, X. Guo, P. Peng, X. Chen, X. Li, and X. Zhou,New Journal of Physics , 073015 (2019).[45] N. Fl¨aschner, B. S. Rem, M. Tarnowski, D. Vogel, D.-S. L¨uhmann, K. Sengstock, and C. Weitenberg, Science , 1091 (2016).[46] S. Jin, W. Zhang, X. Guo, X. Chen, X. Zhou, and X. Li,arXiv 1910.11880 (2019).[47] C.-H. Park and N. Marzari, Phys. Rev. B , 205440(2011).[48] M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. T.Barreiro, S. Nascimb`ene, N. R. Cooper, I. Bloch, andN. Goldman, Nature Phys. , 162 (2015).[49] L. Asteria, D. T. Tran, T. Ozawa, M. Tarnowski, B. S.Rem, N. Fl¨aschner, K. Sengstock, N. Goldman, andC. Weitenberg, Nature Phys. , 449 (2019).[50] H. M. Price and N. R. Cooper, Phys. Rev. A , 033620(2012).[51] N. Goldman and J. Dalibard, Phys. Rev. X , 031027(2014).[52] N. Goldman and J. Dalibard, Phys. Rev. X5