Nonlinear Bloch-Zener oscillations for Bose-Einstein condensates in a Lieb optical lattice

We investigate Bloch-Zener oscillations and mean-field Bloch bands of a Bose-Einstein condensate (BEC) in a Lieb optical lattice. We find that the atomic interaction will break the point group symmetry of the system, leading to the destruction of the Dirac cone structure, while the flat band is preserved on the highly symmetric lines. Due to the nonlinear effect, a tubular band structure with a flat band will appear in the system. Furthermore, comparing with that the tight-binding (TB) model fails to describe the interacting bosonic systems in the honeycomb lattice, we show that the TB model is applicable to study the nonlinear energy band structures for the Lieb lattice. In addition, we show that the loop structure can be determined by the observation of the chaos of the state in the Bloch-Zener oscillations.


I. INTRODUCTION
Two-dimensional condensed matter such as graphene and similar materials [1][2][3] with Dirac cones at Fermi energy, in which low-energy excitations at high symmetry points can be described by the Dirac equation, are of great theoretical and experimental interests due to their unique electronic properties. Recent advances in exploring spin-orbit coupling and artificial gauge fields for neutral atoms [4][5][6][7][8][9] make possible the creation and manipulation of Dirac cones in optical lattices [10][11][12][13]. Different lattice geometries of one or two dimensions in ultracold atomic gases [10,14,15] have been realized experimentally, such as honeycomb lattice with Dirac points at the corner of Brillouin zone (BZ) [16], the Lieb lattice with a flat band [17], and the kagome lattice [18,19], etc. Furthermore, some well-known models in condensed matter physics such as Harper-Hofstadter model and Haldane model [20][21][22][23][24][25] have been theoretically investigated and experimentally realized in optical lattices with promising controllability. The band structure of the cold atom systems can be reconstructed by combining techniques of Bloch-Zenor oscillations with adiabatic mapping of cold atoms, and the existence of Dirac points can be revealed by the momentum-resolved inter-band transitions [26][27][28][29][30][31][32][33].
On the other hand, with high purity and tunability, the ultracold atomic system is one of the best quantum simulators [4,34,35] to study nonlinear dynamics for a BEC with weak interaction. For instance, the lattice structure can be changed by adjusting the laser intensity, while the interaction can be manipulated by tuning s-wave scattering length [35]. In this paper, a nonlinear Gross-Pitaevskii (GP) theory is derived by applying the mean-field approximation to a weakly coupling system. * Electronic address: lizhiphys@126.com Besides the atomic systems, GP equation can also be used to study other nonlinear systems such as photonic systems with Kerr nonlinearity, exciton-polariton condensates, acoustic cavities, and circuit resonators [36][37][38][39][40][41][42]. Previous studies have shown that the mean-field Bloch bands for a BEC in an optical lattice possess unique loop structures [43][44][45][46]. As predicted in Ref. [47], while Dirac points in a honeycomb lattice are broken by the nonlinearity and an intersecting tubular structure appears around the K point, the TB model, however, fails to describe the nonlinear dynamics in this system [47]. In this article, we study BEC in an optical Lieb lattice. The Lieb lattice features a diabolic single dirac cone with a flat band, distinguishing from most systems in which the Dirac points come up in pairs (the so-called Fermion doubling) [48,[50][51][52][53]. The flat band is attributed to the bipatite nature of the lattice [54,55], leading to exotic physics which is entirely determined by the interaction and topology due to constant dispersion caused by destructive interference [56,57]. Thus the nonlinear Lieb lattice may serve as a natural, and possibly indispensable extension of the honeycomb lattice considered in Ref. [47]. Our results show that the nonlinear Dirac cone for the Lied lattice has a similar tubular structure to that of Ref. [47]. Furthermore, the Dirac points in our model are protected by a plane point group symmetry. Any small interatomic interaction will break the symmetry, then the Dirac points will be lifted. And more importantly, the TB model works well in the vicinity of the M point for our system unexpectedly, which is sharply different with the results in the honeycomb lattice explored in Ref. [47].
The interplay of nonlinear superfluidity and relativistic Dirac dynamics leads to many novel effects, such as adiabatic failure, chaos [45,58] and dynamic instability [59], which can be utilized to reveal information on the nonlinear energy bands. In this article we prove that the Bloch-Zener oscillations for the BEC in an optical lattice can be used as an efficient tool to study the structure of the nonlinear energy bands. The regime of fold bands can be determined by detecting the time evolution of the population of sublattices during Bloch oscillations, while the nonlinearity strength could be measured by the interband transition probability of specific tunneling process.
The paper is organized as follows. In Sec. II, we numerically calculate the energy spectra from both the GP Hamiltonian and the Bloch Hamiltonian of the TB model. And we discuss how the nonlinearity breaks the Dirac cone from the perspectives of symmetry. In Sec. III, we consider atoms performing Bloch-Zener oscillations. The transition probabilities for different tunneling process are given both in the adiabatic limit and sudden limit. Finally, a short conclusion is given in Sec. IV. The mean-filed description of a condensate in superfluid phase leads to the Gross-Pitaevskii (GP) equation,

II. SUPERFLUID IN LIEB LATTICE
where m is the mass of particle, a s is the scattering length and V latt (r) is the optical lattice potential of configurations of Lieb lattice (see Fig. 1). The nonlinear strength is denoted as U ≡ 4π 2 a s /m in following discussions. The lattice V latt (r) could be experimentally realized by superimposing three pairs of laser beams with the following formation [17,48,49], where k L = 2π/λ is the wave number of the lattice and λ the wavelength of the lasers. The potential amplitudes V (x,y) long , V (x,y) short and V (x,y) diag could be tuned by adjusting the laser intensities, φ x,y and ϕ are the phase of laser beams. For simplicity we choose (V The Landau instability of the BEC superfluid is not much affected by the lattice strength [61]. As the lattice strength increases, the system goes into the tight binding regime. Following the methods outlined in Ref. [62][63][64], we map the GP equation (1) onto a discrete nonlinear Scrödinger equation to contruct an effective tight-binding Hamiltonian. We write the bosonic field as a sum over the three sub-lattices is the Wannier function and r i,α (α = 1, 2, 3) are lattice vectors in the three sub-lattices of the i-th cell. Then the tight-binding Hamiltonian for our BEC system readŝ where i, j runs over all nearest-neighboring sites, α/β includes the sublattices in one unit cell, is the hopping constant and g is the nonlinear strength, for simplicity we take t α,β = t and do not consider the next-nearest hopping in the main text. In experiments, one may upload 87 Rb atoms which weigh m = 1.443 × 10 −25 kg into the optical lattice. The lattice can be created by a laser of wavelength λ = 1064 nm. Thus the recoil energy is estimated by E R = 2 k 2 L /2m ∼ 2π × 2 kHz. And the tunneling rate t ∼ 1 kHz for a lattice depth V = 2E R in the absence of nonlinearity [60].
In momentum space, where d(k) = (d x , d y ) denotes the Bloch vector with d x = 2t cos(k x a/2) and d y = 2t cos(k y a/2), and S = (S x , S y ) are spin-one matrices with S x = λ 1 , S y = λ 6 from the Gell-Mann matrices (for explicit form, see Appendix A), is the nonlinear term. The Bloch Hamiltonian H k in absence of nonlinear term gN has energy dispersions, as illustrated in Fig. 1 (c), which are in good agreement with Bloch bands of the continuous system Eq. (1), see Appendix B. Note that the low energy effective Hamiltonian H 0 (q) has aĈ z 4 rotation symmetryĈ z In fact the Lieb lattice has a 4mm plane point group symmetry. The quasi-momentum k is invariant with operations in little group 4mm at M and 2mm at X. And the invariant line L Σ and L ∆ are protected by mirror symmetry m d and m x respectively, thus indegenerate, while the Dirac point at M is protected by the 4mm symmetry featuring a nontrivial two dimensional irreducible representation. However, in presence of a nonzero nonlinear term gN , even weak interaction will violate this symmetry, lifting the degeneracy at M and resulting in two more additional crossings, as shown in Fig. 2.
We numerically calculate the Bloch bands of the GP equation Eq. (1) with the Bloch wave solutions which are of the form ψ k (r) = m,n c mn e i(k+Gmn)·r , here G mn = mb 1 + nb 2 is the reciprocal lattice vector. The results are shown in Fig. 2, which are in good agreement with that of the tight binding model H k in the vicinity of M . As predicted by previous works [47], the Bloch bands of a nonlinear two dimensional system possess tube structures and the Dirac point is extended into a closed self-crossing loop. The bands of our system have similar structures around the M point. As illustrated in Fig. 2 (g), the energy bands consist of three "tubes", which intersect at point M . They have wedged cross-section with size increasing monotonically from L Σ to L ∆ (L ∆ denotes the invariant line cross M and parallel to k x -direction) as shown in Fig. 2 (g), (ii)-(iii), and overlap with each other to produce the structure of Fig. 2 (g), (i/iv). For g = 0 the Dirac point at M accidentally degenerates with the flat band, which produces a triple point. As shown before, the Dirac point at M will be lifted to two crossing points D 1 and D 2 for any g = 0 as the nonlinearity breaks the symmetry. However, the loop structures at the zone of Brillouin zone, as illustrated in Fig. 2 (c) Fig. 2 (a) and (c). Furthermore, the flat bands of the tight binding model for g = 0 and g = 0 on L Σ share the same eigen-states, only with a uniform shift of the eigen-value from E = 0 to E = g/2 for any finite g.

III. LANDAU-ZENER TRANSITION
Now we consider atoms performing Bloch oscillations with a constant applied force F. As these atoms are accelerated in the vicinity of the triple Dirac point, the tunneling probability to other bands is finite. This is a problem first considered by Landau and Zener, usually referred to as Landau-Zener (LZ) transition [65][66][67][68][69].
The effective Hamiltonian can be obtained from Eq. (4) by performing Peierls substitution k = Ft, whereĤ = H k − F ·r. The Peierls substitution can also be understood in a semiclassical way that the quasimomentums are subjected to a constant forcek = F [30]. An initial state |ψ(t 0 ) = |u n (k 0 ) in the n-th band evolves as where k(t) = k 0 + Ft. The force F =k linearly increases momentum k of atoms and leads to a non-adiabatic transition to other bands near the crossing (or anti-crossing) point. From the finial state |ψ(t f ) the transition probability is defined as where m, n = 1, 2, 3 and |u m (k(t)) is m-th eigen-state of For simplicity, we start with considering atoms performing one Bloch oscillation in k x directions. We numerically calculate the interband transition probabilities as a function of force F x (see Fig. 3). The applied force is given by F = F xêx and atoms initially prepared in k = 0. The evolution time is set to be T = 2π/(aF x ). Thus the momentum of atoms increases form k 0 = 0 to k f = 2π/a under the influence of the force, along path of which the system has the spectrum shown in Fig. 2 (f).
For g = 0, if the atoms take shorter time to travel across the anti-crossing than the Zener tunneling time T lz = √ δmax(1, √ δ)/∆, the system will evolve adiabatically, where δ = ∆ 2 /(4F x ) is the adiabaticity parameter and ∆ is the minimal energy splitting [67]. For the parameters adopted in this paper, this corresponds to F x < 1. In presence of nonlinearity, the adiabaticity only holds for the bands whose structure is not modified by the nonlinearity. As seen in Fig. 3 (a) and (c), the interband transition probabilities are vanishing while the force is weak enough. However, as shown in Fig.3 (b), the adiabaticity breaks down when atoms are initially prepared in the second band, since the nonlinearity changes the topological structure of this band.
Furthermore, the tunneling probability shows an irregular oscillation for atoms initially in the band with a loop structure, as observed in similar nonlinear systems studied in previous works [45,71]. This irregularity is thought to be associated with the chaotic behavior of the nonlinear system [45]. As shown in Fig. 3 (d), the evolution of state exhibits chaotic features, which could be directly observed in the laboratory by detecting the particle number density of the sublattices by means of hybrid time-of-flight images. For diagonal sweeping, the oscillation of state emerges shortly after the atoms travel across the band loops. Thus we could determine the regime of the fold in the energy bands. In fact, the time evolution of H k at the loops is unstable, which leads to the chaos in the LZ process. As illustrated in Fig. 4, the norm mode for the evolution of a eigen-state is cyclic. However, the periodicity breakdowns for the states on the fold, which implies the instability of the dynamics. The dynamical instability is attributed to the superfluidity of the BEC, which can be understood by the response of the superflow to a perturbation governed by the Bogoliubov equation [70], or the semiclassic fixed point theory of the Hamiltonian-Jacobi matrix [45]. We adopt the latter method and give a brief introduction in Appendix C. As suggested by a previous work [61], the Landau instability of superfluidity in a two-dimensional lattice is direction independent while the dynamical instability is not. Furthermore, Landau instability occurs beyond a critical wave number k c , and the superfluidity goes into a more stable phase as the nonlinear interaction increases.
Next we consider tunneling process along the diagonal path L Σ , i.e., the applied force is given by F = F xêx + F yêy , where we set F x = F y = F . The spectrum on this sweeping path has the structure shown in Fig.  2 (d). The atoms are still initially prepared in k = 0 and the evolution time still set to be T = 2π/(aF x ). Fig. 5 presents the transition probability as a function of the force F. Different from the linear transition process, the relation of transition probability P jn = P nj (see Appendix D) does not hold for g = 0, since the nonlinear term gN breaks the particle-hole symmetry (or the symmetric structure of the spectrum). And particularly, as shown in Fig. 5 (c) and (f), the results for atoms initially in the middle band remain the same as that of g = 0, as the nonlinearity does not change the structure of the flat band. Moreover, atoms prepared in the up or down band do not tunnel into the middle band when they travel across M, just like their behaviors in the absence of nonlinearity (see Appendix. D). We consider the LZ transition acorss a single M point (gapless for g = 0) in sudden limit, i.e., F x = F y = F 1. As the force F becomes large, the time that atoms take to travel across M point is much shorter than the Zener tunneling time. Thus the atoms stay in the initial state in most of the evolution time. The transition probability that depends on the nonlinear strength g is displayed in Fig. 6. We find that difference between the transition probability for g = 0 and g = 0 narrows as the sweeping velocity increases. And such difference is vanishing in the sudden limit, as shown in Fig. 6 (c), since the transition probability does not relate much to the structure of the spectrum. We could determine the nonlinearity g of the condensate by fitting the experimental data with the interpolation function of P jn numerically.

IV. CONCLUSIONS
In conclusion, we have numerically studied the energy band structure of a BEC in an optical Lieb lattice and have provided a useful TB model. It is shown that under nonlinearity, the Dirac points at M point are broken and extended into a closed curve, while the flat band still exists on the high-symmetry line. We have shown that the Bloch-Zener oscillations could be utilized as a useful tool to explore the nonlinear bands: (i) the size of metastable fold bands can be determined by detecting the time evolution of the population of sublattices during Bloch oscillations, (ii) the nonlinearity strength could be measured through the interband transition probability of specific tunneling process. Since the interatomic interaction can be easily manipulated, the observation of the phenomena studied here is expectable in the table-top ultracold atomic experiments. The Gell-Mann matrices used in this paper is given by At k M = (π, π), an expansion of H k for small momentum shift q gives rise to an effective Hamiltonian which has the universal form where q = (q x , q y ) denotes deviation of the Bloch wave number q away from k M . And by a simple replacement of q j = i∂ j with j = x, y, we obtain the motion equation in real space, where Ψ = (|Ψ 1 | 2 , |Ψ 2 | 2 , |Ψ 3 | 2 ) T and H(i∂ j ) is H q with the replacement q j = i∂ j (j = x, y).
Appendix B: NUMERICAL BAND STRUCTURE CALCULATION FOR g = 0 We numerically study the band structure ab initio with the plane wave truncate approximation. We take Fourier transformation and write down the Hamiltonian in momentum space. The eigen-states have the form ψ k (r) = m,n c mn e i(k+Gmn)·r . And in reciprocal space the potential Eq. (2) is writen as V (r) = mn V mn exp (iG mn · r), with components V mn = − (0.25V Then by diagonalizing the Hamiltonian H k (m1,n1;m2,n2) = δ m1,m2 δ n1,n2 E 0 k+Gm 1 ,n 1 + V m1,m2,n1,n2 , (B2) we obtain the eigen-energies for a specific k, where E 0 q = 2 q 2 /2m denotes the kinetic energy. The results are shown in Fig. 7. Here we choose recoil energy E R = 2 k 2 L /2m as the energy unit. The lowest three bands are quite well approximated by the tight binding model. Bands of Lieb lattice can be gapped for some lattice configurations (V The nonlinear three-level model Eq. 4 can also be described as a classical Hamiltonian system, similar to Ref. [45,71]. [cos(k x ) √ s 1 cos(θ 1 ) + cos(k y ) √ s 2 cos(θ 2 )] , (C1) whereθ 1 = θ 1 − θ 2 ,θ 2 = θ 3 − θ 2 are relative phase, and we treat k x and k y as parameters. s1,θ 1 and s2,θ 2 are two pairs of canonically conjugate variables which are governed by, (C5) By settingṡ 1 =ṡ 2 =θ 1 =θ 2 = 0, the eigenenergy is obtained from = H. And the chaotic nature of the nonlinear system can be revealed by the Hamiltonian-Jacobi matrix, which is obtained by linearizing the above motion equations at the fixed points (corresponding to the eigenstates), We numerically solve the eigenvalues λ of H J along the invariant line L Σ . The results are shown in Fig. 8. In general, the eigenvalues are a complex number and only solutions with pure imaginary values are stable. We can see that near M point the real parts of the eigenvalue become nonzero, which implies the emergency of the unstable fold structure.
And the energy spectrum is given by (k) = 0, ±2t 1 + δ 2 + (1 − δ 2 ) (cos k x a + cos k y a) /2 . (D2) Without loss of generality, we consider the motion along the k x direction. The constant force is given by F x = k x0 + F t, and the evolution time is chosen to be t ∈ [0, 2π/(aF )]. Following the calculations outlined in Ref. [45,65], we can obtain the transition probability for H q , where P jn is the occupation probability on the j band at t = T f for the state initially on the n band at t = T 0 , and ∆ = 2 1 2 (1 − δ 2 ) (cos(k y ) − 1) + δ 2 + 1. And we have P jn = P nj due to the symmetry of the levels.
We calculate the transition probabilities for g = 0 with atoms performing a single Bloch oscillation in k x directions. The results are shown in Fig. 9. Since the parallel momentum k y is conserved, the term 2t cos(k y a/2)S y can be seen as an effective mass term. Thus P jn is modified with different initial momentum k y0 . In particular, for system in gappless phase, i.e., δ = 0, and atoms prepared in initial momentum k = 0, the travel across the triple Dirac point leads to a unit transferred probability P +− = 0 (P +0 = 0) as shown in Fig. 9 (d).