Breakdown of the Peierls substitution for the Haldane model with ultracold atoms

We present two independent calculations of the tight-binding parameters for a specific realization of the Haldane model with ultracold atoms. The tunneling coefficients up to next-to-nearest neighbors are computed ab-initio by using the maximally localized Wannier functions, and compared to analytical expressions written in terms of gauge invariant, measurable properties of the spectrum. The two approaches present a remarkable agreement and evidence the breakdown of the Peierls substitution: (i) the phase acquired by the next-to-nearest tunneling amplitude $t_{1}$ presents quantitative and qualitative differences with respect to that obtained by the integral of the vector field A, as considered in the Peierls substitution, even in the regime of low amplitudes of A; (ii) for larger values, also $|t_{1}|$ and the nearest-neighbor tunneling $t_{0}$ have a marked dependence on A. The origin of this behavior and its implications are discussed.

j i Adr} [3]. The latter expression must be evaluated on the straight path connecting sites i and j, as demonstrated under the hypothesis of a same-site, same-orbital interaction with the vector field by Boykin et al. [4].
Despite its popularity, the Peierls substitution is a rather uncontrolled approximation, as already pointed out in Refs. [5,6]. For example, we notice that the integral of the vector field appearing in the Peierls phase factor has been conventionally taken along a straight path (see e.g. [7,8]) long before its formal demonstration [4], just for convenience (in principle, in two and three dimensions there is an ambiguity as the path is not univocally defined [3]). In addition, in the literature the Peierls substitution is often applied as a "magic formula", with little care about its regime of validity.
The Peierls substitution plays a fundamental role in the Haldane model [7], a celebrated two-dimensional periodic tight-binding model, characterized by a quantum Hall effect caused by the breaking of time-reversal symmetry with zero magnetic flux through the unit cell [7]. The model is characterized by exotic quantum phases, with different Chern numbers, depending on the value of the phase ϕ of the next-to-nearest tunneling amplitude t 1 , that is usually computed by the integral of the vector field A cited above. Recently, in the literature there have been proposals for engineering the Haldane model with ultracold atoms in optical lattices by means of artificial gauge fields [9,10], and to study the associated topological quantum states in the presence of sharp boundaries [11]. In fact, these systems represent a very interesting platform for simulating solid state physics [12]. Again, these studies make use of approximate methods to deal with the tunneling amplitudes, by exploiting the Peierls substitution tout court [9] (see also [13][14][15][16][17][18][19][20][21]), or by using approximate atomic orbitals [11].
In this Letter we present two independent calculations of the tight-binding parameters for the Haldane model discussed in Refs. [9,11]. In particular, we show that, within the next-to-nearest neighbors approximation, the tunneling coefficients can be directly written in terms of gauge invariant, measurable properties of the spectrum (namely the gap at the Dirac point and the bandwidths), or computed ab-initio by using the maximally localized Wannier functions (MLWFs) [22,23]. Notably, the two approaches present a remarkable agreement, evidencing the breakdown of the Peierls substitution. As a matter of fact, the phase acquired by the next-to-nearest tunneling amplitude t 1 is quantitatively different from that predicted by the integral of the vector field A, and presents a pronounced dependence on the intensity of the underlying scalar potential. Moreover, both the amplitude of t 1 and of the nearest-neighbor tunneling t 0 turn out to be dependent on the intensity of A.
Let us start from the following single-particle, minimalcoupling Hamiltonian in two-dimensionŝ and under rotations by θ = 2π/3 radians around any vertex of the lattice. The former implies that next-tonearest tunneling amplitudes t1 along the same direction are conjugate pairs (solid and dashed lines); from the latter follows the equivalence of the hopping amplitudes separated by 2π/3 radians. When sites A and B are degenerate, the system is also invariant under rotations by π radians around the center of any elementary cell; this implies that t0 is real.
with r = (x, y), p = −ih∇, and V L being the following honeycomb potential [11,23,24] 3e y ) generate the reciprocal lattice, k L is the laser wavelength and s the amplitude of the potential in units of the recoil energy E R =h 2 k 2 L /2m [25]. Notice that, though this specific realization is characterized by degenerate potential wells, an imbalance can be easily produced by introducing a suitable phase [9,24]. The corresponding Bravais lattice, B = {j 1 a 1 +j 2 a 2 j 1 , j 2 = 0, ±1, ±2 . . . }, with lattice constant a (such that k L = 4π/(3 √ 3a) [24]), is generated by the two basis vectors a 1/2 = (2π/3k L )(e x , ∓ √ 3e y ), obeying a i ·b j = 2πδ ij . As for the vector potential, we consider the same expression discussed in Refs. [9,11] (corresponding to the Coulomb gauge, ∇·A(r) = 0) that has the same symmetry of the underlying honeycomb potential (see Fig. 1). The parameter α represents the amplitude of the vector potential in units ofhk L .
The tight-binding model is constructed from the manybody HamiltonianĤ 0 = drψ † (r)Ĥ 0ψ (r), by expanding the field operator on a basis of localized functions, ψ(r) ≡ jνâ jν w jν (r), with the usual commutation rules [â jν ,â † j ν ] = δ jj δ νν . Then, by restricting the analysis to the two lowest bands,Ĥ 0 can be written as [23] where the matrix elements w jν |Ĥ 0 |w j ν correspond to tunneling amplitudes between different lattice sites (except for the special case j = j, ν = ν , representing the onsite energies). These matrix elements depend only on j − j due to the translational invariance of the lattice. The spectrum ofĤ 0 can be obtained by considering the following transformation from coordinate to momentum with h νν (k) = j e ik·R j w 0ν |Ĥ 0 |w jν , and S B indicating the first Brillouin zone [23]. By truncating the above expression to next-to-nearest neighbors as usual [7,9], we define The first term corresponds to the onsite energies, The second term has only off-diagonal elements, corresponding to the hopping toward the three nearestneighbor sites (see Fig. 1). Thanks to the symmetries of the Hamiltonian (1), the three tunneling amplitudes are equal. By defining t 0 ≡ w 0A |Ĥ 0 |w 0B , we can write and h 21 (k) = z * (k). Finally, by defining and taking again into account the symmetries of the system (see Fig. 1), the last term -corresponding to nextto-nearest tunneling between homologous sites -can be cast in the following form Notice that in general the onsite energies and the tunneling coefficients depend on the amplitudes of both the scalar and vector potentials: . This is a direct consequence of the fact that the optimal choice for the basis of localized functions w jν (r) depends on the properties of overall structure of the Hamiltonian (1). By defining we can write that is equivalent to the expression discussed in Ref. [9]. However, we remark that here we have not made explicit use of the Peierls substitution, and that the dependence of Eq. (10) on the phase ϕ is a consequence of the symmetries of the full potential. Finally, by diagonalizing the matrix h νν (k) and defining f ± (k) ≡ (|t 1A |F A (k) ± |t 1B |F B (k))/2, we get the following expression for the spectrum of the lowest two bands that is a function of |t 0 |, |t 1ν |, and ϕ ν . In the following we will consider for simplicity the degenerate case = 0 (E A = E B ), corresponding to the potential in Eq. (2). In this case, thanks to the symmetries of the system, we have |t 1A | = |t 1B | ≡ |t 1 |, ϕ A = −ϕ B ≡ ϕ and w 0A |Ĥ 0 |w 0B = w 0B |Ĥ 0 |w 0A (when A and B are equivalent the system is invariant under rotation by π radians around the center of any cell, see Fig. 1). The latter implies that t 0 is real. Remarkably, in this case the two tunneling amplitudes t 0 and |t 1 | and the phase ϕ can be expressed in terms of specific properties of the spectrum. Let us start by noticing that f + (0) = 6|t 1 | cos ϕ, f − (0) = 0, |z(0)| = 3t 0 . In addition, we indicate with k D the position of the Dirac points [23], and define ∆ ± ≡ ±( ± (0) − ± (k D )), that correspond to the two bandwidths when the tunneling coefficients satisfy the hierarchy t 1 t 0 . Then, we have with δ D ≡ + (k D ) − − (k D ) being the gap at the Dirac points, due to the presence of the vector potential. Also, at e.g.
Another relation containing |t 1 | and ϕ is Then, by combining Eqs. (15) and (16), we get Eqs. (14), (17) and (18) represent an important contribution of this work: they provide a way to connect the value of the tunneling amplitudes to gauge-invariant, measurable properties of the spectrum. Moreover, they also provide a straightforward method for computing the tunneling amplitudes, as the exact Bloch spectrum can be be readily computed by means of a standard Fourier decomposition [23,24], even in the presence of a vector potential [26].
In addition, we compare these values with those computed ab-initio from their definition in terms of the matrix elements w jν |Ĥ 0 |w j ν . To this end, we make use of the MLWFs for composite bands [22], which are defined through the following unitary mixing of the two lowest Bloch bands U νm (k)ψ mk (r), (19) with R j ∈ B, ψ mk being the eigenfunctions of the Hamiltonian (1) [6], and U νm (k) a 2×2 unitary matrix, periodic in k-space, which minimizes the spread of w jν (r) [22]. In the present case, the MLWFs are obtained by modifying the code discussed in Ref. [23] in order to include a vector potential. The MLWFs turn out to be complex due to the breaking of time-reversal, and this explains the emergence of a phase factor in the tunneling coefficients [26]. The values obtained for t 0 , |t 1 | and ϕ are shown in Figs. 2, 3, along with those extracted from the spectrum. The agreement is remarkable [27]. From these figures we can identify two regimes as a function of the amplitude α of the vector potential: (i) for small enough values, α < ∼ 1, where t 0 and |t 1 | are almost constant and the phase ϕ is linear in α; (ii) for α > ∼ 1 where the dependence on α is less trivial. In particular, in the latter regime, t 0 and |t 1 | present a pronounced dependence on α, in clear contrast with the Peierls substitution (see horizontal lines in Fig. 2) which assumes the phase ϕ to be the only α-dependent quantity. However, this dependence is not surprising, as the presence of the vector potential may significantly affect both the Bloch eigenfunctions ψ mk [6] and the gauge transformation U νm entering Eq. (19) [5], so that the usual implicit assumption that the basis of localized orbitals is not affected by the vector potential (see e.g. [4]) is generally not valid. On the other hand, the calculated phase strongly deviates from the linear behavior expected from the Peierls substitution, namely ϕ = Fig. 3. This figure reveals that the Peierls substitution dramatically fails even in the "linear" regime, as it predicts a slope for the phase far much larger than the actual one. Moreover, it completely neglects its dependence on the amplitude s of the scalar potential (that is appreciable even in the full tight-binding regime, s > 10). This is particularly evident from Fig. 4, where we plot the behavior of the angular coefficient in the linear regime, dϕ/dα| α=0 , as a function of s. This figure provides fur- Note the logarithmic scale on the vertical axis. We remark that the present tight-binding model with up to nearest-neighbor tunnelings is accurate only for s > ∼ 3; for lower values it may be necessary to consider also other next-to-leading tunneling coefficients [23].
ther evidence that the Peierls substitution does not even provide a reasonable estimate for the order of magnitude of ϕ in the linear regime. Essentially, the reason for the breakdown of the Peierls substitution resides in the fact that the hypotheses under which it has been rigorously demonstrated [2,4] cannot be satisfied in the Haldane model. Most importantly, the vector potential can not be considered as slowly varying [8], as it varies on the same length scale as the lattice (see Fig. 1). As a consequence, both the scalar and vector potentials must be treated on equal foot, and all parameters (t 0 , |t 1 | and ϕ) must be considered as dependent on both s and α.
In summary, we have presented two independent calculations of the tight-binding parameters for the Haldane model with ultracold atoms [9], one based on their abinitio definition in terms of the MLWFs, and the other in terms of gauge invariant properties of the spectrum, summarized in Eqs. (14), (17) and (18). The latter provides a straightforward approach whenever the spectrum can be measured or computed with sufficient accuracy. The results obtained with the two methods present a remarkable agreement, and demonstrate the inadequacy of the Peierls substitution, which fails in predicting quantitative and even qualitative properties of the system. The reason for this breakdown is due to the fact that the regime of validity of the Peierls substitution cannot be fulfilled in any realization of Haldane model, regardless of the system, being it cold atoms in optical lattices or electrons in a solid. Our results indicate that a careful revision of the validity of the commonly employed Peierls substitution in tight-binding models is necessary.