Wannier-Stark flatbands in Bravais lattices

We systematically construct flatbands (FB) for tight-binding models on simple Bravais lattices in space dimension $d \geq 2$ in the presence of a static uniform DC field. Commensurate DC field directions yield irreducible Wannier-Stark (WS) bands in perpendicular dimension $d - 1$ with $d$-dimensional eigenfunctions. The irreducible bands turn into dispersionless flatbands in the absence of nearest neighbor hoppings between lattice sites in any direction perpendicular to the DC field. The number of commensurate directions which yield flatbands is of measure one. We arrive at a complete halt of transport, with the DC field prohibiting transport along the field direction, and the flatbands prohibiting transport in all perpendicular directions as well. The anisotropic flatband eigenstates are localizing at least factorially (faster than exponential).

All these phenomena originate from perturbations lifting the macroscopic degeneracy and breaking the destructive interference, which is at the origin of flatbands. Appearance of destructive interference requires either fine-tuning or symmetries, that enforce the interference [30]. A marked feature of FB models with short-range hopping are strictly compact eigenstates, called compact localized states (CLS) [31]. Their presence greatly simplifies the analysis of the FB models [32,33], * marindam@ibs.re.kr † aalexei@ibs.re.kr and can be used as a foundation for their systematic classification [34][35][36][37].
A completely different scenario unfolds when a Bravais lattice is exposed to a dc field, which generates an infinity of (d − 1)-dimensional bands each supporting eigenstates embedded into d-dimensional spaces. In this work we analyze models with applied dc fields that do not need fine-tuning to achieve band flatness, and FB eigenstates that cannot be arranged into CLS. Applying a static field to a 1D tightbinding chain leads to the appearance of a Wannier-Stark (WS) ladder of equidistant eigenenergies for the spectrum of the chain, with all eigenstates being localized and the dynamics of observables in general displaying time-periodic Bloch oscillations [38]. The quest for nontrivial states in higher dimensional lattices with magnetic and electric fields resulted in the observation of flatbands for a square lattice for certain electric field directions [39]. Similar dispersionless features were later identified for a rectangular lattice [40,41]. Here we present a systematic construction of WS flatbands for the five two-dimensional Bravais lattices: Oblique, rectangular, centered rectangular, triangular, and square lattices. We obtain the dependence of the band structure of the static field direction, analyze the localization properties of the flatband eigenstates, and demonstrate that, unlike conventional FB systems which host CLS, WS flatband states cannot host CLS and are at least factorially localized instead. We then generalize these results to higher lattice dimensions and discuss the impact of longer-range hoppings.
The paper is organized as follows: Sec. II introduces the model for 2D lattices and the definitions that are used later in Sec. III to derive the band structure of the models for different directions of the field. Conditions for the bands flatness are also discussed in the same Sec. III, while the eigenstates are analyzed in Sec. IV. We extend our analysis to higher lattice dimensions in Sec. V. We conclude with a summary and open issues.

II. SETTING THE STAGE
We consider a 2D Bravais lattice (see Fig. 1): whereâ 0,1 are the lattice basis vectors and (n, m) are the translation indices of the lattice point. The vectors need not be orthogonal, for example in the case of the triangular lattice, and can be expressed in terms of Euclidean basis vectors (ê 0 ,ê 1 ) aŝ a 0 = γ 0ê0 ,â 1 = γ 1 (cos θê 0 + sin θê 1 ).
Here 0 < θ ≤ π/2 is a tilting angle of the axis, and γ 0,1 are the lengths of the basis vectors of the lattice unit cell. For square and triangular lattice γ 0 = γ 1 . For all other cases (oblique, rectangular and centered rectangular) they are not equal. For rectangular and square lattices θ = π 2 , for triangular lattice θ = π 3 , and for oblique and centered rectangular lattices θ can take any value except π 2 and π 3 . Next we define a tight-binding Hamiltonian on the lattice in the presence of a static dc field E: which acts on the Hilbert space spanned by the basis vectors {|n, m : (n, m) ∈ Z × Z}. The indices l, j denote the hopping range. Our main goal is to compute and analyze the spectrum of this Hamiltonian. Because of the presence of the dc field, discrete translation invariance is in general broken and we do not expect any eigenenergy band structure. However the character of the spectrum depends on the direction of the dc field, which we set to Here F is the strength of the field and x, y ∈ R. Then, parallel and perpendicular directions to the field are encoded by the respective unit vectors n = (xγ 0 + yγ 1 cos θ)ê 0 + (yγ 1 sin θ)ê 1 We choose the direction of the field such that n ⊥ is parallel to one of the lattice vectors p 1â0 +p 2â1 with p 1 , p 2 ∈ Z. This ensures that the translational invariance of the lattice persists in the perpendicular direction n ⊥ , albeit with a different lattice spacing than the original lattice.
In what follows we coin such field directions as commensurate directions. They ensure the existence of a 1D band structure in the spectrum of the model. Commensurate field directions constrain the possible values of x, y and θ to either of the two possibilities as discussed in Appendix A.
In order to diagonalize the Hamiltonian [Eq. (3)], we exploit the partial translation invariance and introduce a new rotated coordinate system (see Fig. 2): the z coordinate along the dc field, and the w coordinate perpendicular to the dc field that we define as follows z = α γ 2 0 x 2 + γ 2 1 y 2 + 2γ 0 γ 1 xy cos θ r nm · n = nx α + my α , with and β = (γ 0 γ 1 sin θ) −1 . For any commensurate field direction there exists a value of the coefficient α such that both x α , y α become mutually prime integers [as can be directly verified by computing the ratio y α /x α using the conditions (6), see Appendix B]. Then we can parametrize the new coordinates as (z, w) = (z, w 0 (z) + η(xx α + yy α )), z, η ∈ Z. (9) with a function w 0 (z), see Appendix B for details. Finally we note that for any coprime (x, y) (as happens for any dc field commensurate direction on a square lattice) the distances ∆ z and ∆ w in real space covered by changing respectively, z to z + 1 and η to η + 1, amount to ∆ z = γ 2 0 x 2 + γ 2 1 y 2 + 2γ 0 γ 1 xy cos θ,

III. THE 2D SPECTRUM
The convenience of the z, w basis is that now the partial translation invariance of the Hamiltonian H [Eq. (3)] in the direction perpendicular to the field is made explicit: In the above we have also defined The Hamiltonian is invariant under the shifts w → w , but not under the shifts z → z . Therefore we can apply the Bloch's theorem in the direction perpendicular to the field. For that we define a complete orthonormal set of basis vectors where z ∈ Z and k ∈ [0, 2π), the quasi-momentum in the direction perpendicular to the field. Then the action of the Hamiltonian H in Eq. (3) on a basis vector is given by which correspond to an eigenvalue E. The eigenvalue is a function of the quasimomentum k: A. Generating function method Using Eq. (14) we write the eigenproblem and solve it by introducing a generating function We note that the generating function g E (q, k) is nothing but the Fourier transformed ψ E (z, k) from z-space to qspace. The eigenproblem, Eq. (17), transforms into an ordinary differential equation This equation can be solved for any pair of (E, k). The solutions are not periodic functions in q for generic E, k.
The periodicity in the q requirement from Eq. (19) imposes a dispersion relation of WS bands between E and k.
We assume that the hopping networks on each of the five Bravais lattices respect inversion symmetry: This assumption is made for convenience only, the con-clusions presented below hold as well if the symmetry is broken. Then, defining λ E = E/F, s lj = t lj /F we arrive at This differential equation can be easily integrated, and we discuss in the next section the possible solutions.

B. Flatbands and Dispersive Bands
The set of commensurate field directions splits into a subset of flatband directions and the complementary subset of dispersive directions. All flatband directions satisfy the condition This condition is equivalent to requesting the absence of any direct hopping connection between two lattice sites perpendicular to the chosen commensurate field direction in Eq. (11). Indeed, for this case the solution of Eq. (22) reads as Enforcing 2π periodicity of g E in q results in i.e., all bands are flat, equidistant, and labeled by an integer index a. We note that the generating function remains periodic in q and the bands remain flat even in the absence of the inversion symmetry (21). For a dispersive field direction lx α + jy α = 0 for at least one pair (l = l , j = j ), it follows Again enforcing the 2π periodicity of g E in q we get: All bands turn dispersive, equidistant and labeled by an integer index a. For both cases the band gap between two consecutive energy bands is equal to F for any k. The prefactor A(k) is periodic in k: It is instructive to observe that the dc field strength enters completely additive in the dispersion relation (27), leaving the irreducible E(k) dependence invariant. Also, for any short-range network the number of flatband directions is infinite, while the number of dispersive ones is always finite and equal to the number of hopping connections on the network. Adding more hoppings to a given Bravais lattice network will add more dispersive field directions on the expense of the flatband directions. The limiting case of connecting all sites with all (even though the hopping strength may decrease in a suitable way with increasing distance between sites) will eliminate all flatband directions and leave us with dispersive field directions only.

IV. LOCALIZATION OF 2D FLATBAND EIGENSTATES
We next turn to the analysis of the eigenstates of the WS flatbands. Flatbands enjoy macroscopic degeneracy as there is no unique choice of the eigenstate basis. Eigenstates of flatbands in short-range translationally invariant Hamiltonians can be typically arranged into compact localized states [2,31]. We note that in such translationally invariant cases the number of eigenstates of one flatband equals the embedding space dimension of its eigenvectors. For WS flatbands the situation differs, as the embedding space dimension for eigenvectors is infinitely larger than the number of eigenstates of one flatband. It appears impossible to assemble a linear combination of WS flatband eigenstates, which turns compact in real space. Indeed, let us assume that a compact localized state does exist. Then, the generating function g(q, k) can be expanded into a double Fourier series with a finite number of components in both q and k. This contradicts the general solution obtained in Eq. (24). Therefore, compact localized states are ruled out-see Appendix E for details. What is then the best localization which can be achieved with WS flatbands?
Let us attempt to identify the most localized eigenstates. The eigenstates are extracted from the generating function g E (q, k) for a fixed band a. We set a = 0 and E = 0 without loss of generality. Using the property of the Bessel function of first kind we can express the generating function as where R denotes the set of the hoppings l, j up to the inversion/reflection symmetry (21), e.g., for nearest neighbor (n.n.) hopping these sets are {(1, 0), (0, 1), (1, −1)} and {(1, 0), (0, 1)}, respectively, in the case of triangular or centered rectangular lattices and all the other lattices. In this case ν denotes the set of all possible integers ν (l,j) for all the hoppings (l, j) ∈ R. Next we write the flatband basis states in the z, k representation subject to the constraint (l,j)∈R ν (l,j) (lx α + jy α ) = −z .
Because the WS band is flat any linear combination of the basis vectors (31) is an eigenstate of the Hamiltonian H (3). Let us choose c k = [2πA(k)] −1 to remove the normalization factor. It follows subject to the constraints For the five Bravais lattices with n.n. hopping and either four or six n.n. neighbors the above generic expressions can be further simplified. The case of four nearest neighbors-square, oblique, and rectangular latticesthe real space wavefunction is given by products of pairs of Bessel functions: The eigenfunctions on triangular and centered rectangular lattices are given by sums over products of triplets of Bessel functions: The details of the derivations are given in Appendix D.
We are interested in the decay properties of the above wave functions along and perpendicular to the field direction. Recall the asymptotics of Bessel functions J ν (t) ∼ 1 |ν|! t 2 |ν| for large order integer |ν|. Since the above wave functions involve products of Bessel functions, we conclude that the spatial decay will be at least factorial 1/r! in any lattice direction, which is faster than any exponential decay. Let us consider the square lattice with s = t/F and t the nearest hopping strength. From Eq. (10) it follows with ζ ∈ Z. We therefore arrive at the wavefunction asymptotics for |ζ| → ∞ |Φ(z = 0, w = ζ∆ w )| ≈ 1 |ζx|!|ζy|! |s| |x|+|y| |y| |x| |x| |y| |ζ| , (40) Since ∆ w = ∆ z and for nonzero integers |x| = |y| the term |x/y| |x|−|y| > 1, the flatband eigenstates always decay faster along the field direction as compared to the perpendicular one. The only exception is |x| = |y| = 1, for which the decay in both directions is the same. We plot the wave function profile for the square lattice and field direction x = 2, y = 1 in Fig. 3(a) in loglinear scale. The wavefunction decays faster along the field direction than in the perpendicular one. The inset shows that the gradient is monotonically decreasing as a function of position for both directions, i.e., the decay is faster than exponential (for exponential decay it would be a step-function around the origin).
We analyzed numerically a similar case for the triangular lattice, which is shown in Fig. 3(b). The field direction is again x = 2, y = 1. Again we observe that the wavefunction is decaying faster along the field direction than in the perpendicular one, and in both directions faster than exponential (see inset). We therefore expect that other Bravais lattices and commensurate field directions will yield similar localization properties of flatband eigenstates. It appears reasonable that the decay in the perpendicular direction is slower than in the field direction, as the onsite energies of the original lattice vary linearly with distance along the field direction as opposed to the perpendicular one.
It is straightforward to observe that we can define commensurate field directions for higher d-dimensional Bravais lattice in the same way: Choose any two points on the lattice separated by a finite distance, then the connecting line corresponds to an allowed perpendicular direction. In this case the quasimomentum along the perpendicular direction k will be replaced by a d − 1 dimensional vector k in the generating function in Eq. (18), keeping the 2π periodicity of the generating function with respect to q: g E (q + 2π, k) = g E (q, k). The differential equation [Eq. (22)] will remain similar after replacing k by k. In the absence of hopping connections along the perpendicular direction, the solution of the differential equation [Eq. (22)] turns where A( k) and f (q, k) are periodic functions of k and q, and λ E = E/F. Therefore, similar to the 2D case, the periodicity condition g E (q + 2π, k) = g E (q, k) implies λ E takes all possible integer values, and hence the entire spectrum degenerates into an infinite set of WS flatbands.
In the presence of hopping connections along the chosen perpendicular direction, the bands will acquire nonzero dispersion. For all commensurate field directions each band is characterized by (d − 1)-dimensional wavevectors k, while the wavefunctions are embedded into a d-dimensional space. We expect that flatband eigenstates will be localizing faster than exponential in all directions, and slower in the perpendicular directions as compared to the field direction.
As an example, we consider the case of the cubic lattice with nearest-neighbor hopping t that respects the inversion symmetry, Eq. (21). The Hamiltonian then is a 3D version of Eq. (3). The derivation of the generating function, and the spectrum in this case, is a straightforward generalization of that for the square lattice. Therefore, we only provide the final results. For a commensurate field direction E ∝ F (x 1 , x 2 , x 3 ) T with x 1 x 2 x 3 = 0 and gcd(x 1 , x 2 ) = gcd(x 1 , x 3 ) = 1, the generating function reads ( k = (k 1 , k 2 )), where λ E = E/F, s = t/F, and τ 2 x 1 + τ 1 x 2 = 1 for fixed integer values of x 1 , x 2 . Similarly, to the 2D case the 2π periodicity in q of the generating function fixes the eigenenergies and implies that the spectrum consists of 2D flatbands E( k) = aF with a ∈ Z. The inverse Fourier transform of g E yields the eigenfunction (a = 0), that clearly displays factorial decay in the lattice coordinates (n, m, l).

VI. DISCUSSION AND CONCLUSIONS
We considered the impact of a dc field on the spectra of tight-binding models on d-dimensional Bravais lattices. For commensurate field directions (for which the perpendicular direction is parallel to some lattice vector) the spectrum consists of an infinite number of equidistant (d − 1)-dimensional bands. A finite number of commensurate directions yields dispersive bands. The remaining infinite set of commensurate directions leads to dispersionless Wannier-Stark flatbands. The flatband wavefunctions are embedded in a d-dimensional vector space, and cease to form compact localized states, which is the usual scenario for translationally invariant lattices with short-range hoppings. As a result, Wannier-Stark flatband eigenfunctions decay factorially in space, and typically faster along the field direction than perpendicular to it.
Our results are applicable for ultracold atoms in optical lattices, where the electric field is substituted by a tilt of the lattice in the gravitational field [42] or acceleration of the whole lattice [43]. Notably, the same type of perturbations can be arranged in optical waveguide arrays where the electric field is modeled by a curved geometry of the waveguides [44]. In both cases experimental platforms for two-dimensional settings have been developed. Such experiments can test the sensitivity of choosing commensurate field directions, the existence of flatband field directions, Bloch oscillations, and the factorial localization of flatband eigenstates. Another intriguing set of potential applications could be related to electronic transport in strained two-dimensional materials, which is, e.g., a hot topic in graphene research [45].
Some open problems for Wannier-Stark flatbands on Bravais lattices are their fate in the presence of perturbations such as disorder, magnetic fields, and many-body interactions. The method used in this work might not be applicable in these cases, and one would need to resort to the generic analytical and numerical methods. However, a possible way to analyze the impact of the perturbations is to perform perturbation calculations around the limit of strong dc field strength, which allows the approximation of the factorially localized eigenstates by compact ones and use of the methods developed for the analysis of disordered and/or interacting flatband models [8,32,33,46]. An even more interesting problem is the case of non-Bravais lattices where there is more than one site per unit cell. The generating function approach extends naturally to non-Bravais lattices [47] and again leads to a Wannier-Stark ladder of an irreducible band structure, now consisting of several bands. Can some of these bands be tuned to become flat? A positive answer exists for some chiral lattices, which transport the chiral symmetry into the generating function and the irreducible band structure [47]. Similar to translationally invariant chiral flatbands [30], the irreducible Wannier-Stark band structure will contain a chiral flatband. The eigenfunctions will be noncompact as in the Bravais lattice case. At variance to the Bravais case they were observed to localize only exponentially [47], perhaps due to the presence of other dispersive nonflat bands in the irreducible spectrum. Can we finetune non-Bravais lattice hoppings such that the irreducible Wannier-Stark band structure turns one, or several, or even all, bands flatwithout imposing a symmetry like the chiral one? We think these are exciting questions for future research. is parallel to some lattice vector indexed by p 1 , p 2 ∈ Z: p 1â0 + p 2â1 = (p 1 γ 0 + p 2 γ 1 cos θ)ê 0 + p 2 γ 1 sin θê 1 .
After some simple algebra one arrives at the required result: For x α y α = 0 we absorb gcd(x α , y α ) in α to make x α , y α coprime for the convenience of our analysis. In the case x α y α = 0, we can choose either x α = 1 while y α = 0 or y α = 1 while x α = 0.
Therefore z is integer for all (n, m) ∈ Z × Z. On the other hand, w takes discrete values, that are not necessarily integer or can not be made integer by rescaling of w = ny − mx by a constant factor, since x/y is not a rational number in general. Let us denote the set of all possible pairs (z, w) by S and construct its parametrization. We show that z can take any integer value, and we can parameterize w for a fixed value of z. As mentioned before, either x α y α = 0 (for which we can make either x α = 1 or y α = 1) or x α and y α are mutually prime. For fixed z, we pick one lattice point that corresponds to z, w 0 (z), and is indexed by (n 0 , m 0 ) According to Bézout's identity from number theory [48], there is always a choice of τ 1 , τ 2 such that This implies z = λ, and it takes any integer value. The values of τ 1 , τ 2 can be determined for example using the Euclidean division algorithm [49]. Note that despite the fact that the choice of τ 1 , τ 2 is not unique, the condition (B6) remains unchanged under the simultaneous change of τ 1 by τ 1 − px α and τ 2 by τ 2 + py α , where p can be any integer. The corresponding value of the perpendicular coordinate is For a given integer value of z all the other values of w are generated by simultaneous shifts of n from n 0 and m from m 0 by the following way: Therefore the entire set S of valid lattice points is parameterized as [50] (z, w) = (z, w 0 (z) + η(xx α + yy α )), z, η ∈ Z. (B9) We use the parametrization (B9) in the main text. In the following we provide the values of parameters τ 1 , τ 2 for some special cases.
For some simple cases we can even guess restrictions on the set S. For example, when x, y are integers the set S is only a subset of Z × Z and does not contain all of its elements. One can check that by simply looking at the square lattice case, where z = nx + my, w = ny − mx. In that case for the field direction x = 2, y = 1 the point (z, w) = (0, 1) does not exist in the original lattice since it corresponds to fractional indices (n, m) = (1/5, −2/5).
The generating function g E (q, k) is periodic both in k and q, and hence it can be expanded as a Fourier series in the variables k and q: g E (q, k) = p1,p2∈Z g p1,p2 e ikp1+iqp2 .
Compactness of the eigenstate in the field direction implies g p1,p2 = 0 except for a finite number of values of p 2 . But from the solution Eq. (24) it follows that the generating function cannot be expressed as finite polynomial in e iq . Hence a flatband wavefunction cannot be made compact in the direction of the field.
We will now discuss the possibility for the flatband wavefunction to be compact in the perpendicular direction of the field. Since we study lattice eigenvalue problems, C CLS (k) is a 2π-periodic function in k: Therefore, |−p 2 , w 0 (−p 2 ) − (p 1 + p 3 )∆w .
Compactness in the direction perpendicular to the field implies that the product C p3 g p1,p2 = 0 except for a finite number of integer values of the sum (p 1 + p 3 ). Simple inspection of Eq. (24) yields that, for any fixed value of p 3 with a corresponding nonzero C p3 , there always exists an infinite number of p 1 values for which g p1,p2 turns nonzero.
Therefore there exists no function C CLS (k) for which the flatband eigenfunction turns into a CLS. Moreover, we proved that the flatband eigenfunctions are necessarily non-compact in all space directions. The proof can be generalized to higher space dimensions d ≥ 2 by replacing k with a (d − 1) dimensional vector k.