# Self-consistent tight-binding description of Dirac points moving and merging in two dimensional optical lattices

## Abstract

We present an accurate ab initio tight-binding model, capable of describing the dynamics of Dirac points in tunable honeycomb optical lattices following a recent experimental realization [L. Tarruell et al., Nature 483, 302 (2012)]. Our scheme is based on first-principle maximally localized Wannier functions for composite bands. The tunneling coefficients are calculated for different lattice configurations, and the spectrum properties are well reproduced with high accuracy. In particular, we show which tight binding description is needed in order to accurately reproduce the position of Dirac points and the dispersion law close to their merging, for different laser intensities.

## I Introduction.

The possibility to simulate graphene-like structures and investigate the physics of Dirac points with ultra cold atoms in optical lattices is attracting an increasing interest in the literature (1); (2); (3); (4); (5); (6); (7); (8); (9); (10). Recently, Tarruel et al. (8) have reported the creation and manipulation of Dirac points in a tunable honeycomb optical lattice, exploring the topological transition occurring at the merging of Dirac points, and comparing their experimental results with ab initio calculations of the Bloch spectrum. The same experiment has also been interpreted by means of an universal tight-binding Hamiltonian defined on a square lattice (9), that describes the merging of Dirac points and the corresponding topological transition between a semi-metallic phase and a band insulator (5). Though this model remarkably captures all the relevant physics, its connection with the optical lattice parameters is indirect and has some limitations, as it relies on a fit of the parameters of the universal Hamiltonian to the two lowest lying energy bands in the vicinity of the Dirac cones (9); (10).

In this paper, we present a comprehensive scheme based on composite maximally localized Wannier functions (MLWFs) (11) for constructing an ab initio tight-binding model corresponding to the tunable honeycomb potential with two minima per unit cell as described in Ref. (8). The MLWFs are obtained by means of a gauge transformation that minimizes their real space spread, and are routinely employed in condensed matter physics (12). As recently demonstrated, these functions represent also an optimal tool for constructing tight-binding models for ultra cold atoms in optical lattices (13); (14); (15), as they allow for an optimized mapping of the system hamiltonian onto the discrete model defined on a lattice. They allow for an ab initio calculation of the tight-binding parameters, with a fine control over next to leading corrections.

The paper is organized as follows. In section II, we review the general approach for mapping a continuous many-body Hamiltonian onto a discrete tight-binding model by means of MLWFs. In section III we apply this scheme in the tunable two-dimensional honeycomb lattice of Ref. (8), discussing the general structure of the associated Bravais lattice in direct and reciprocal space, and presenting the tight-binding expansion up to the third-nearest neighbors. Here we also give explicit numerical results for the calculated MLWFs, the inferred tunneling coefficients and the Bloch band structure. Then, in section IV, we make use of the tight-binding Hamiltonian in reciprocal space - expressed in terms of tunneling coefficients and functions depending on the geometry of the associated lattice - for discussing the behavior of the Dirac points as a function of the lattice parameters. There, we also discuss the dispersion relation close to the merging of Dirac points, refining and improving the analysis of Refs. (5); (9). In addition, we discuss the effect of parity breaking, that destroys the degeneracy of the two potential minima, providing a finite Dirac mass. Some accessory but important details are included in the Appendices, covering the tight-binding expansion, the gauge dependence of the results, and the numerical application.

## Ii Tight binding expansion and MLWFs

Let us consider a many-body system of bosonic or fermionic particles, described by the field operator . In particular, as the physics of Dirac points is determined by the single-particle spectrum, here we will consider the non-interacting many-body Hamiltonian

(1) |

with and the lattice periodic potential , where belongs to the associated Bravais lattice. We remark that the optimal Wannier basis is solely determined by the single particle spectrum, and even though the inclusion of an interaction would be straightforward (16); (13); (15), the derivations following the next lines would not be affected.

The Hamiltonian (1) can be conveniently mapped onto a discrete lattice corresponding to the minima of the potential by expanding the field operator in terms of a set of functions localized at each minimum,

(2) |

where is a band index, and () represent the creation (destruction) operators of a single particle in the -th cell, satisfying the usual commutation (or anti commutation) rules following from those for the field .

The MLWFs, introduced by Marzari and Vanderbilt in (11), are obtained trough a unitary transformation of the Bloch eigenstates

(3) |

with the volume of the first Brillouin zone and a gauge transformation which obeys periodicity conditions in order to preserve the Bloch theorem. This gauge transformation is obtained through the minimization of the Marzari-Vanderbilt localization functional (11). The resulting MLWFs posses the desirable property of being exponentially localized in real space, (18); (19) thus constituting an ideal basis for tight-binding models (14). In this article we consider the MLWFs for composite bands (N1) since we are interested in geometries where the Wigner-Seitz cell has a non trivial basis. This allows each MLWF to be centered on a single potential minimum inside the elementary cell, in contrast to single-band Wannier functions (13); (14); (17). For this work, the MLWFs have been computed by means of the WANNIER90 code (20); (12) and a modified version of the QUANTUM-ESPRESSO package (21) adapted to the case of an optical lattice (14). We also mention that other methods for specific cases have been recently proposed (13); (15).

The Hamiltonian (1) can be written in terms of Wannier states as

(4) |

where the matrix elements depend only on due to the translational invariance of the lattice. They correspond to tunneling amplitudes between different lattice sites, except for the special case , that corresponds to onsite energies. Then, by defining

(5) |

is transformed as

(6) |

with

(7) |

being the Hamiltonian density in quasi-momentum space, whose eigenvalues are in principle equal to the exact energy bands . For practical purposes, however, the expression (7) must be truncated by retaining only a finite number of matrix elements. This corresponds to the tight-binding expansion in -space. The actual number of terms needed to reproduce the energy bands (or any other physical quantity) within a certain degree of accuracy crucially depends on the properties of the basis functions , the MLWFs being the optimal choice due to their minimal spread, see Appendix A.

From now on we will apply the above tight-binding expansion to the tunable two-dimensional honeycomb potential of the experiment (8).

## Iii Tunable honeycomb lattice

The functional form of the potential reproduced experimentally in Ref. (8) is

(8) | |||

where all the parameters can be controlled and tuned in the experiment. In particular, by varying the laser intensities and , several structures can be realized by continuous deformations, ranging from chequerboard to triangular, dimer, honeycomb, and square lattices, including 1D chains.

Let us define the Bravais lattice associated to the potential minima as (see Fig. 1) which is generated by the two basis vectors

(9) |

with being the cartesian unit vectors. Therefore, the basis vectors in reciprocal space are

(10) |

following from . From now on, we can fix , , without loss of generality. This corresponds to measure lengths in units of and energies in units of the recoil energy (14).

As the unit cell in direct space generally contains two basis points and , we consider the mixing of the two lowest bands. It is then customary to write the Hamiltonian density in Eq. (7) as

(11) |

where (see also Eq. (2)) the band index has been traded to as the resulting MLWFs are centered at subwells located in . The two lowest energy bands are then given by the eigenvalues of (11)

(12) |

with .

The matrix elements in (11) can be expanded as

(13) | |||||

(14) |

with

(15) | |||||

(16) |

corresponding to diagonal and off-diagonal matrix elements, respectively. The sign convention is chosen in such a way that all the coefficients appear positive defined (14). Here we truncate the tight-binding expansion by including all possible tunneling between neighboring cells, as indicated in Fig. 1 (22).

Let us start by considering the diagonal terms. By fixing an arbitrary energy offset, we can write

(17) | |||||

(18) |

with

(19) |

and

(20) | |||||

The tunneling coefficients appearing in Eq. (20) are precisely those connecting the minima located at points of the same type or (see Fig. 1), and have been redefined as follows

(21) |

in order to simplify the notations. The form of the function follows from the explicit form of the corresponding lattice vectors .

We notice that when the minima in are degenerate in energy, so that . Also, , so that , and the eigenvalues in Eq. (12) take the following simple form

(22) |

As far as the off-diagonal matrix element , its analytical form is

(23) | |||||

where the tunneling coefficients have been redefined as

(24) |

Notice that the ordering of the tunneling coefficients in Eqs. (24),(21) does not necessarily correspond to the hierarchy of their magnitudes, as this may depend on the regime of the potential parameters (see later on).

In the following, we will use this model to discuss the features of the stretched honeycomb configuration (23), which allows us to analyze the behavior of the Dirac points. We will also discuss the effect of increasing the overall potential intensity (in order to enter a well defined tight-binding regime, but with the same potential structure) and that of breaking the degeneracy between sites of type and .

## Iv MLWFs and tunneling coefficients for the degenerate case

In this section we will discuss the numerical results for the stretched-honeycomb configuration with two degenerate minima per unit cell, obtained with , in (8). This is the most interesting configuration due to the presence of massless Dirac points (see later on). The effect of parity breaking (), that generates a Dirac mass, will be analyzed in the next section. In addition, in Appendix B we will present the results for a wider range of lattice configurations.

Here we start by considering the experimental regime of Tarruel et al. (8), namely , and variable ranging from to . Within this range of parameters, the potential (8) has the stretched-honeycomb structure shown in Fig. 2(a). This configuration determines the shape of the calculated MLWFs, drawn in Fig. 2(b) for the sublattice of type . As shown in this figure, the MLWF is exponentially localized around the site of the central unit cell (note the logarithmic scale), and it presents a non-negligible contribution around the neighboring potential minima, as well. The associated tunneling coefficients are presented in Fig. 3(a). This figure shows the behavior of the diagonal and off-diagonal coefficients, () and () in Eqs. (24) and (21) respectively, as a function of .

With these parameters we then compute the tight-binding energy dispersion in Eq. (22), shown in Fig. 4 for some examples. In particular, we consider two different tight-binding approximations, one including just , and (corresponding to the universal hamiltonian of Ref. (5)), and that including all the coefficients in Fig. 1. In Figs. 4(a-b) we compare the energy dispersion of the different tight-binding approximations with the exact spectrum at . The figures show that the main features, including the band-crossing along the direction (Fig. 4(b)), are well reproduced by both approximations, though the tight binding model with just , and is not capable of approximating the exact bands with sufficient accuracy (this holds in all the range of considered in this paper). In any case, remarkably the model is able to reproduce with sufficient accuracy the position of the Dirac points even in this parameter range, as shown in the next section.

Here we also consider a different set of values for the potential parameters that correspond to a well defined tight-binding regime, while maintaining the stretched-honeycomb structure. In particular, we use the parameter values , and ranging from to , corresponding to twice the values of Tarruel et al. (8). The calculated tunneling coefficients are illustrated in Fig. 3(b), showing the same general structure as the ones in Fig. 3(a), except for minor differences regarding the smallest coefficients. The corresponding energy dispersion is shown in Figs. 4(c-d) for as an example. In this case, already the lowest order approximation with just the coefficients just , and provides a remarkable agreement with the exact data. This has been verified also for the other values of in the range considered.

## V Dirac points

As we have seen in the previous section, the spectrum for a stretched honeycomb configuration with is characterized by points where the two bands are degenerate, with a linear dispersion along at least one direction - the so-called Dirac points. They are defined by and come always in pairs due to time-reversal invariance, with (5). Their existence and position depends on the geometry of the lattice: for a regular honeycomb structure (a graphene-like lattice) they are located at the corners of the Brillouin zone (4); (14), whereas in the present tunable case they can be moved inside the Brillouin zone, as showed in (8). In particular, from the expression in Eq. (23), the position of the Dirac points is obtained by solving the following equation

(25) |

whose imaginary part yields

(26) |

inside the first Brillouin zone. Then, Eq. (25) becomes

(27) |

solved by

(28) |

taking into account the considered hierarchy of the tunneling coefficients indicated in Fig. 3). In the current regimes, this expression can be further approximated as

(29) |

corresponding to the expression of Ref. (5). Both Eq. (28) and its approximate version Eq. (29) provide a valid solution when , which is satisfied also in the range of parameters corresponding to the stretched honeycomb, as shown in Fig. 5.

### v.1 Merging of Dirac points

The merging of Dirac points occurs when the two solutions in Eq. (28) coincide modulo a reciprocal space vector (with , integers, see Eq. (9)), namely at . Therefore, the merging point satisfies (5)

(30) |

In principle, due to the geometry of the lattice, there are four possible inequivalent merging points, at (5), see Fig. 6. However, for the actual values of the tunneling coefficients, only the point and its equivalents are possible. In particular, in our case the Dirac points inside the first Brillouin zone merge (with those of outer cells) at the top and bottom corners and , namely for . For the two examples considered here, see Fig. 5(a) and (b), the merging occurs at (see also (8)) and , respectively.

Following (5); (9), we now expand the hamiltonian density around one of the two merging points, defining . As discussed in Ref. (9), the general form of the off-diagonal component around a merging point is characterized by a linear term in and a quadratic one in , coming respectively from the imaginary and real parts of . Namely, the leading terms of the expansion are

(31) |

Here we also take into account the diagonal term , not included in the approach of (9); (5), as it affects the quadratic behavior introducing an asymmetry between the two bands; neglecting an unimportant constant term we have

(32) |

Therefore, close to the merging point the hamiltonian density can be cast into the form

(33) |

with

(34) | |||||

(35) | |||||

(36) | |||||

(37) |

The corresponding dispersion law is

(38) |

with being vanishing at the merging point, marking the topological transition between semi-metallic and insulating phases driven by a change of sign in the product (5); (9). The expression (38) provides indeed a good approximation of the exact Bloch energies close to the merging point, as shown in Fig. 7 (a similar expansion can be derived around a generic Dirac point). In this picture we consider , that is the top corner of the first Brillouin zone in Fig. 6, and show band-cuts along orthogonal directions at (upper panels) and (lower panels). Panels (a),(c) show two Dirac points belonging to adjacent Brillouin zones, symmetric with respect to (see Eq. (28)). This corresponds to a positive (in the present regime of parameters is always negative). By decreasing the two Dirac points approach each other and eventually merge when , see panels (b,e). At this particular point and the dispersion law is linear along and quadratic along . By further decreasing , a gap opens at the merging point, see panels (c),(f). In this case the mass-like term is characterized by a negative . In all the panels, the bands are compared with the low-energy expansion (34), showing a fine agreement close to the merging point (except panel (d), see caption). In particular, a small asymmetry in the quadratic behavior along is visible and well reproduced close to the merging point, owing to the diagonal term proportional to in Eq. (37) (24).

### v.2 Breaking parity: massive Dirac points

As shown in the experiment (8), a gap can be opened at the Dirac points by breaking the invariance under parity, achieved by tuning the angle away from . In this case, due to the asymmetry of two minima in the unit cell, the diagonal terms and () are no longer degenerate, see Appendix D. This causes the Dirac particles to acquire a mass, as it is evident from Figs. 8,9 where we show the energy bands for two Dirac points, at and at the merging point . This figures show that even small deviations from give rise to a significant mass term (gap) at the Dirac points, and that the the current tight-binding model accurately reproduce the exact energy bands. Note that in this case the full formula (12) has to be used. We also mention that even in this case one can derive an expansion analog to that in Eq. (38).

## Vi Conclusions

Maximally localized Wannier functions for composite bands (11) are a powerful tool for constructing tight binding models for ultra cold atoms in optical lattices. Here, we have considered the tunable honeycomb optical lattice of Ref. (8), and we have shown how to derive the corresponding tight-binding hamiltonian, ab initio. We have calculated the MLWFs and the tunneling coefficients for different lattice configurations, showing that the spectrum properties, including the position of Dirac points and the dispersion law close to their merging, can be reproduced with high accuracy with an expansion up to third-nearest neighbors. We have considered both cases of massless and massive Dirac points, respectively for the case of two degenerate minima per unit cell and for the case of parity breaking. These results provide a direct connection between the experimental results of Ref. (8) and the universal hamiltonian of Refs. (5); (9).

###### Acknowledgements.

This work has been supported by the UPV/EHU under programs UFI 11/55 and IT-366-07, the Spanish Ministry of Science and Innovation through Grants No. FIS2010-19609-C02-00 and No. FIS2012-36673-C03-03, and the Basque Government through Grant No. IT-472-10.## Appendix A Approximate Bloch spectrum

Here we briefly review and comment about the use of the MLWFs for the tight binding expansion of the exact Bloch spectrum, following the discussion in Ref. (11) (see also (12)). Let us start by rewriting the the hamiltonian density (7) as follows (see Eq. (3))

(39) | |||

where are the (exact) Bloch bands, and the are periodic, unitary matrices representing gauge transformations (as a function of quasimomentum) of the Bloch states.

Single bands. For transformations that do not mix the bands, namely when

(40) |

the onsite energies and tunneling coefficients in (39) are independent on the (periodic) phases . Furthermore, if all the terms in the sum are retained, owing to the following formula valid for an infinite lattice

(41) |

one easily recovers the exact diagonal expression . This result is trivial, following directly from the completeness of the Wannier basis. As a consequence, the tunneling coefficients can be expanded in terms of the exact energies (with no reference to the Wannier functions) (27); (28). In addition, the tight binding approximation of the exact Bloch spectrum, namely the truncation of the sum in (39) at a given order, is independent on the choice of the Wannier states. So, in the absence of band mixing (the gauge group being a direct product of groups) the tight binding expansion is gauge independent.

Composite bands. Let us now consider the case of composing bands, via a non-abelian gauge transformation. Again, summing over all lattice sites the hamiltonian density in Eq. (39) takes the form

(42) |

whose eigenvalues coincide with the exact bands owing to the unitarity of the transformation. Moreover, even the trace of in Eq. (39) at a finite order of the tight binding expansion is gauge-independent, as it does not dependent on the transformation matrices. Instead, finite order approximations of individual Bloch bands (or the sum of a subset of them) are gauge-dependent, as they depend on a particular choice of the matrices .

Parallel transport gauge. We recall that the transformation matrices are defined as those that minimize the Wannier spread , and that the latter can be decomposed as (11). The first term is gauge invariant, while – in case of composite bands – the second can be written as the sum of two (non negative) diagonal and off-diagonal components, . Both and can be written in terms of the generalized Berry vector potentials , defined as (25); (26)

(43) |

with the matrix being hermitian. In one-dimension (1D), the gauge in which is minimized corresponds to as a consequence of the vanishing of off-diagonal () Berry connections in Eq. (43), and is called the parallel transport gauge. In this case, the transformation can be obtained directly by requiring the off-diagonal connections to be vanishing (26), so that the optimal tight binding expansion can be seen as a direct property of the Hilbert space. Nevertheless, as the Wannier functions depend on the choice of the gauge, the two points of view are correlated. We remark that this approach is generally limited to 1D cases as in higher dimensions it is not always possible to make vanishing (11), so that the parallel-transport formulation can not be easily generalized. However, though in absence of a formal proof, in general we may assume that the gauge where the spread of Wannier functions is minimal, is the one that provides the best tight binding approximation of individual Bloch bands. In fact, this has been already verified in a number of models (13); (14); (15).

Finally, we remark that the use of composite instead of single band transformations is required in case of a set of almost degenerate bands (well separated from the others), that usually corresponds to more that one minimum per unit cell, as in the present case. A more thorough discussion on this point, for the case of a 1D double well potential, can be found in (13).

## Appendix B MLWFs and tunneling coefficients

In this appendix we analyze the properties of the MLWFs and the associated tunneling coefficients in a range of broader than the one covered in the text. This allows us to analyze the two opposite limits of the potential (8) corresponding to the dimmer and 1D chain structures (8). These are exemplified in Figs. 10(a),(b) ( and , respectively) for the experimental regime of Ref. (8). The dimmer structure is characterized by a relatively low value of the potential in the region between and within the same unit cell. On the opposite, in the 1D chain regime the potential is low along the direction connecting different minima, while it presents a barrier between the and sites of the same unit cell. The stretched-honeycomb regime covered in the text (Fig. 2) represents an intermediate configuration between these two limits.

The structure of the potential in these limits determines the shape of the MLWFs, which we illustrate in Figs. 10(c-d) (results shown for sublattice ). As in the stretched-honeycomb structure, the MLWFs are exponentially localized around the site of the central unit cell and present a non-negligible contribution around the neighboring potential minima. In the case , see Fig. 10(c), we find a large contribution of the MLWF around the site of the central unit cell, consistent with the dimmer structure of the potential. The situation is very different for , in Fig. 10(d), which shows a MLWF highly localized along the axis, resembling the 1D chain structure of the potential.

In order to analyze the degree of localization of the MLWFs, in Fig. 11 we show the spread of the MLWFs, (11), as a function of . The figure shows that, by increasing , rapidly decreases in the regime of low , while it almost saturates in the opposite limit. This indicates that the tight-binding approach is expected to work better in the stretched-honeycomb and the 1D-chain regimes, rather than for the dimmer case.

The behavior of the tunneling coefficients in the whole range from the dimmer to the 1D chain limits are shown in Fig. 12. We first focus on the left hand side of the graphic, . There, we find that the ratio between the two dominant coefficients is . This reflects the dimmer structure of the potential, since connects sites and (see Fig. 1). Noteworthy, is by far the next biggest coefficient, comparable in magnitude to . This reveals that the tunneling between neighboring dimmers in direction is considerable (see in Fig. 1). The rest of the coefficients have a significantly lower value than , and . We note that in Fig. 12 there are two coefficients, and (see Eq. 16), that were not considered in our expansion. In the dimmer regime, these coefficients can be larger than , and , included in our tight-binding model.

As is increased, the various tunneling coefficients evolve in two different ways. Most of them decrease in magnitude, reflecting the stronger localization of the MLWFs as we approach a more tight-binding regime. This could be termed the ‘normal’ behavior, followed for instance by the tunneling coefficients of a perfect honeycomb lattice (14); (23). However, two coefficients, namely and , increase in magnitude as is increased. This ‘inverse’ behavior reflects the evolution of the potential (8) from the dimmer to the 1D chain structure, as these coefficients connect potential minima inside the 1D chains. Owing to this ‘inverse’ behavior, becomes the dominant coefficient for . Similarly, becomes larger than and even for . Thus, it is clear that varying the potential amplitude can modify the role of the different tunneling coefficients.

## Appendix C Accuracy of the tight-binding models

In Figs. 13(a),(b) we compare the exact and tight-binding energy dispersions along () in the dimer and 1D-chain limits (panels (a) and (b), respectively), for the experimental regime (8). Here we have included the results for the two tight-binding approximations considered in the text (with just , , , and with all the coefficients). As in the stretched-honeycomb case (Fig. 4), the tight-binding model reproduces the main features of the exact dispersion, including the approximate position of the Dirac point in the case of Fig. 13(b) (note that there is no such point in the dimmer limit, Fig. 13(a)). In Figs. 13(c),(d) we show the analogous pictures for the tight-binding regime discussed in the text, that is doubling the potential parameters of Ref. (8). In this case, the agreement with the exact energies when all the coefficients are included is remarkable in both limits.

A further way to test the accuracy of the different tight-binding expansions is to analyze the overall mismatch of the tight-binding energies against the exact ones. Here, we evaluate this mismatch using the following expression (13); (14)

(44) |

where are the exact energies, the average bandwidth and the area of the Brillouin zone.

The calculated mismatch is shown in Figs. 14(a),(b) as a function of , for the experimental and tight-binding regimes, respectively. Overall, the mismatch in the tight-binding regime (b) is one order of magnitude smaller than the one of the experimental regime (a). Remarkably, the best approximation in Fig. 14(b) has an error below in all the range of . This further confirms the adequacy of the tight-binding models in terms of the MLWFs.

We identify two different trends in the behavior of the mismatch. Focusing on Fig. 14(a), we find that for , decreases as is increased. This can be expected, since in this region the MLWFs become much more localized as the potential is raised (see Fig. 11), hence a more tight-binding regime is approached. For , in contrast, the mismatch increases with increasing . We recall from Fig. 12 that the tunneling coefficients corresponding to sites inside the 1D-chains grow as is increased. When approaching the 1D-chain limit, some of these coefficients that are not considered in our tight-binding model may become relevant, hence the quality of the approximation may decrease.

## Appendix D Effect of parity breaking

Here we analyze the asymmetric case corresponding to , considering for simplicity just the tight-binding parameter regime. In this configuration, the two potential minima in the unit cell become non-degenerate (8). In Fig. 15(a) we illustrate the structure of the potential for at the merging point, with the deeper minimum at site B. Correspondingly, the associated MLWFs exhibit a higher localization around B, as illustrated in Figs. (b-e). As a consequence of parity breaking, the degeneracy of the diagonal coefficients is also broken (for both the onsite energies - see Eqs. (17)-(19) - and the tunneling coefficients , ), see Fig. 16. This figure shows the splitting of onsite energies and diagonal tunneling coefficients for small deviations from (the off-diagonal tunneling coefficients are weakly affected in these range of values of ). These variations allows to accurately reproduce the exact dispersion law and in particular the opening of a mass gap at the Dirac points, as discussed in the text (see Figs. 8,9).

### References

- S.-L. Zhu, B. Wang, and L.-M. Duan, Phys. Rev. Lett. 98, 260402 (2007).
- C. Wu, D. Bergman, L. Balents, and S. Das Sarma, Phys. Rev. Lett. 99, 070401 (2007); C. Wu and S. Das Sarma, Phys. Rev. B 77, 235107 (2008).
- B. Wunsch, F. Guinea, and F. Sols, New J. Phys. 10, 103027 (2008).
- K. L. Lee, B. Grémaud, R. Han, B.-G. Englert, and C. Miniatura, Phys. Rev. A 80, 043411 (2009).
- G. Montambaux, F. Piéchon, J. N. Fuchs and M.O. Goerbig, Phys. Rev. B 80, 153412 (2009); G. Montambaux, F. Piéchon, J. N. Fuchs, and M. Goerbig, Eur. Phys. J. B 72, 509 (2009); J. N. Fuchs, L. K. Lim, and G. Montambaux, Phys. Rev. A 86, 063613 (2012); R. de Gail, J. N. Fuchs, M. O. Goerbig, F. Piéchon, and G. Montambaux, Physica B 407, 1948 (2012).
- P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick, W. Plenkers, G. Meineke, C. Becker, P. Windpassinger, M. Lewenstein and K. Sengstock, Nat. Phys. 7, 434 (2011); P. Soltan-Panahi, D.-S. Lühmann, J. Struck, P. Windpassinger, and K. Sengstock, Nat. Phys. 8, 71 (2011).
- K. Sun, W. V. Liu, A. Hemmerich, and S. Das Sarma, Nat. Phys. 8, 67 (2012).
- L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu and T. Esslinger, Nature 483, 302 (2012).
- L.-K. Lim, J.-N. Fuchs, and G. Montambaux, Phys. Rev. Lett. 108, 175303 (2012).
- T. Uehlinger, D. Greif, G. Jotzu, L. Tarruell, T. Esslinger, L. Wang, M. Troyer, Eur. Phys. J. Special Topics 217, 121 (2013).
- N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
- N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
- M. Modugno and G. Pettini, New J. Phys. 14, 055004 (2012).
- J. Ibañez-Azpiroz, A. Eiguren, A. Bergara, G. Pettini and M. Modugno, Phys. Rev. A 87, 011602(R) (2013).
- R. Walters, G. Cotugno, T. H. Johnson, S. R. Clark and D. Jaksch, Phys. Rev. A 87, 043613 (2013).
- D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
- B. Vaucher, S. R. Clark, U. Dorner, and D. Jaksch, New J. Phys. 9, 221 (2007).
- C. Brouder, G. Panati, M. Calandra, C. Mourougane, and N. Marzari, Phys. Rev. Lett. 98, 046402 (2007).
- G. Panati and A. Pisante, arXiv:1112.6197 [Commun. Math. Phys. (to be published)].
- A. Mostofi, J. Yates, Y. Lee, I. Souza, D. Vanderbilt, and I. Marzari, Comput. Phys. Commun. 178, 685 (2008).
- P. Giannozzi et al., Journal of Physics: Condensed Matter 21, 395502 (2009).
- The correspondence with the tunneling coefficients of Ref. (9) is the following: , , .
- We notice that in the case of a perfect honeycomb geometry one has , , , whereas can be safely disregarded in the tight-binding regime, see (13).
- This term is not considered in the universal model of Ref. (5); (9).
- M. C. Chang and Q. Niu, Phys. Rev. Lett. 75, 1348 (1995); Phys. Rev. B 53, 7010 (1996).
- G. Pettini and M. Modugno, Phys. Rev. A 83, 013619 (2011).
- L. He and D. Vanderbilt, Phys. Rev. Lett. 86 5341 (2001).
- G. L. Alfimov, P. G. Kevrekidis, V. V. Konotop, and M. Salerno, Phys. Rev. E 66, 046608 (2002).