# Tight binding models for ultracold atoms in honeycomb optical lattices

## Abstract

We discuss how to construct tight-binding models for ultra cold atoms in honeycomb potentials, by means of the maximally localized Wannier functions (MLWFs) for composite bands introduced by marzari1997. In particular, we work out the model with up to third-nearest neighbors, and provide explicit calculations of the MLWFs and of the tunneling coefficients for the graphene-lyke potential with two degenerate minima per unit cell. Finally, we discuss the degree of accuracy in reproducing the exact Bloch spectrum of different tight-binding approximations, in a range of typical experimental parameters.

###### pacs:

67.85.Hj, 03.75.LmIntroduction. Ultracold atoms in optical lattices are routinely employed as simulator of condensed matter physics, thanks to the possibility of engineering several geometry configurations and of tuning the system parameters with great flexibility and precision bloch2008; lewenstein2007. In particular, honeycomb lattices are attracting an increasing interest as they allow to mimic the physics of graphene, where the presence of topological defects in momentum space, the so-called Dirac points, leads to remarkable relativistic effects wunsch2008; lee2009; soltan-panahi2011; soltan-panahi2011b; gail2012; tarruell2012; lim2012.

Despite the fact that the potentials describing optical lattices are continuous and can be expressed in simple analytic forms as the combination of a number of sinusoidal potentials, from the theoretical point of view it is often convenient to describe the system by means of a tight binding approach on a discrete lattice. In fact, the potential intensity can be tuned to sufficiently high values so to localize the atoms in the lowest vibrational states of the potential wells, justifying a description in terms of tunneling coefficients related to the hopping between neighboring sites (the potential minima), and interaction strengths which characterize the onsite interaction among the atoms bloch2008. This applies also for the case of honeycomb lattices, where a number of tight binding approaches have been considered recently lee2009; gail2012; lim2012; hasegawa2012, in analogy to the case of graphene wallace1947; cloizeaux1963; cloizeaux1964; reich2002.

A crucial ingredient for the connection between the continuos and discrete versions of the system hamiltonian is the existence of a basis of functions localized around the potential minima. This is important not only conceptually - in order to justify the tight binding expansion - but also from the practical point of view, as a precise knowledge of the basis functions is needed to connect the tight binding coefficients with the actual experimental parameters modugno2012. In the case of optical lattices with a cubic-like arrangement (with a single well per unit cell) this basis is provided by the exponentially decaying Wannier functions discussed by Kohn kohn1959; he2001; bloch2008, from which one can derive analytic expressions for the tight-binding coefficients gerbier2005. However, in general this approach may fail when the potential has more than one well per unit cell. For example, for the case of a honeycomb potential with two degenerate minima in the unit cell, the Kohn-Wannier cannot be associated to a single lattice site as they occupy both cells for symmetry reasons cloizeaux1963; cloizeaux1964.

A powerful approach, that is widely used for describing real material structures, is represented by the maximally localized Wannier functions (MLWFs) introduced by Marzari and Vanderbilt marzari1997. The MLWFs are obtained by minimizing the spread of a set of generalized Wannier functions by means of a suitable gauge transformation of the Bloch eigenfunctions for a composite band. Due to their exponential decay brouder2007; panati2011, the MLWFs provide an optimal basis set for tight-binding models and is widely employed in condensed matter physics marzari2012.

In this paper we discuss the tight-binding expansion up to third-nearest neighbors for ultracold atoms in honeycomb lattices, and explicitly calculate the MLWFs and the tunneling coefficients for a potential with two degenerate minima in the unit cell, by using the WANNIER90 package mostofi2008. For the tunneling coefficients we also provide an analytic expression in terms of the lattice intensity, obtained from a fit of the data. Then we discuss the validity of different tight-binding approximations - including only the nearest neighbor tunneling or up to the third-nearest neighbor - in terms of the experimental parameters.

Tight binding approach for the honeycomb potential.

Let us start by considering the two-dimensional graphene-like lattice discussed by lee2009

(1) |

where and ( being the modulus of the laser wavevectors), corresponding to a honeycomb structure, with a diamond-shaped elementary cell with basis and as shown in Fig. 1. The Bravais lattice in real space is generated by two fundamental vectors defined by , so that lee2009.

The starting point for constructing a tight-binding model is the many-body hamiltonian for bosonic or fermionic particles, described by the field operator . In the following we will focus on the single-particle term

(2) |

as the mapping onto the tight-binding model is determined by the spectrum of the single particle hamiltonian jaksch1998. Here represents the potential amplitude in units of the recoil energy .

In general, when the potential wells are deep enough, the hamiltonian (2) can be conveniently mapped onto a tight-binding model defined on the discrete lattice corresponding to the potential minima, by expanding the field operator in terms of a set of functions localized at each minimum, as

(3) |

where labels the cell and is a band index. In Eq. (3) () represent the creation (destruction) operators of a single particle in the cell , and satisfy the usual commutation rules (following from those for the field ).

In order to construct a basis of localized functions at each site of the lattice here we consider the MLWFs for a composite band discussed by marzari1997. These are a set of generalized Wannier functions , defined from a linear combination of Bloch eigenstates , namely marzari1997; modugno2012

(4) |

with indicating the first Brillouin zone, and being a gauge transformation that minimizes the Marzari-Vanderbilt localization functional marzari1997. The real character and the exponential decay of MLWF functions has been demonstrated under very general assumptions brouder2007; panati2011. Thus, these functions represent an optimal basis for tight-binding models. The real character of the calculated Wannier functions has been confirmed in our computations.

In our case, since the honeycomb unit cell contains two potential minima, and , we can construct a basis of MLWFs by considering a composite band consisting of the two lowest Bloch bands, that is . The mentioned Bloch sub-bands have been obtained with a modified version of the QUANTUM-ESPRESSO package espresso intended to solve the single particle Schrödinger equation associated to the external potential (1). We consider a plane wave expansion of the Bloch states, reaching convergence with an energy cutoff corresponding to =10.5 and a k-point mesh of 1414. As a next step, the MLWFs have been computed considering the WANNIER90 program mostofi2008; marzari2012. The typical shape of the calculated MLWFs is shown in Fig. 2.

The strong localization of () around sites and , and their exponential decay are clearly visible in panel (a). Panel (c) shows the distribution of around the original unit cell ; the figure reveals an appreciable overlap of the MLWF with the neighboring and sites, indicated respectively by the yellow and red arrows in Fig. 1. In addition, the MLWFs are characterized by nodes in passing from sites of type A to B, corresponding to a change of sign (see panel (b)).

The approach followed here of including two Bloch bands is the minimal approximation, and corresponds to the generalization of the single band approximation usually employed for cubic-type lattices bloch2008. Within this approximation, the hamiltonian (2) can be written as

(5) |

where the expansion coefficients in the above equation correspond to the onsite energies , and to the tunneling amplitudes between different sites. Both quantities are real thanks to the properties of the MLWFs. Here we will include up to the third-nearest neighbor tunnelings , with (they depend only on the relative distance owing to the uniformity of the lattice). Then, the tunnelings can be divided in three classes, as shown in Fig. 1 for the case .

(i) Terms between and the three nearest neighbors of type , yellow arrows in Fig. 1, e.g.

(6) |

(ii) Terms between sites of the same type ( or ) within neighboring cells, red arrows in Fig. 1, e.g.

(7) |

In general the two tunneling coefficients and are different. As a specific example, here we will explicitly compute them for the case of degenerate minima where . In this case we can set , where the sign is chosen in order to have positive defined (see Fig. 2 and the discussion about the sign of the MLWFs).

(iii) Terms connecting to at opposite corners of the hexagon, blue arrows in Fig. 1, e.g.

(8) |

We remark that the above derivation of the tight-binding model is valid in general for any potential with a honeycomb structure, with two minima per unit cell, and not just for the potential with two degenerate minima in Eq. (1). In the following we will consider explicitly the latter case in order to provide a specific example.

The behavior of the different tunneling coefficients as a function of the lattice intensity is shown in Fig. 3. In order to extract from the numerical values an analytic formula we consider a fit of the type (), in the range , with , , and as fitting parameters. For we find

(9) |

that has to be compared with the semiclassical estimate of lee2009, (see Eq. (38) in lee2009). In the range of considered here, we find that the latter overestimates the actual value in (9) by about for , up to about for . For the other two terms we get

(10) | |||||

(11) |

These three fit are shown as dashed lines in Fig. 3.

Tight-binding spectrum. A convenient way to check the regime of validity of a given tight binding approximation is to compare its prediction for the energy spectrum with the exact Bloch spectrum. The latter can be readily computed by means of a standard Fourier decomposition lee2009. The typical structure of the two lowest bands is shown in Fig. 4. It is characterized by Dirac points at the vertices of the Brillouin zone (a regular hexagon), where the local dispersion is linear and the two bands are degenerate lee2009. For convenience, we fix at the Dirac points.

The tight-binding spectrum can be derived as follows modugno2012. By defining the hamiltonian in Eq. (5) can be written as

(12) |

with . Finally, the matrix turns out of the form lee2009

(13) |

As for the diagonal terms, the leading contribution comes from the onsite energies , that can be conveniently written as , by shifting the total energy by an overall constant. In addition, there is a correction due to , so that eventually we get

(14) |

(15) |

Instead, the off-diagonal terms get the leading contribution from the term proportional to lee2009, and a correction due to , , with

Finally, by diagonalizing the matrix and defining , we get the following expression for the tight-binding spectrum

(16) |

Again, this is valid for a generic honeycomb structure with two minima per unit cell. For the particular case of degenerate minima, as for the potential in (1), this expression further reduces to

(17) |

where the last term has been added in order to make the energy vanishing at the Dirac points, consistently with the definition of the Bloch spectrum. A specific example is shown in Fig. 4(b,c), for . The figure shows that the tight-binding model with just is not sufficient to reproduce the band structure, and that at least the inclusion of is needed. In particular, the latter is necessary to account for the band asymmetry, see Eq. (17). It also remarkable to notice that the values of and can be extracted from a fit of the Bloch spectrum at , by using the expression (17). In fact, we have and . In addition, since is negligible with respect to (see Fig. 3), practically the former expression can be used as an estimate of . Notably, these estimates coincide with the dashed lines shown in Fig. 3, therefore providing an independent check of the tunnelings calculated by means of the MLWFs.

In general, one can evaluate the degree of accuracy in reproducing the exact Bloch spectrum by defining an energy mismatch as follows

(18) |

with being the th bandwidth (). The results are shown in Fig. 5. This figure shows that the tight-binding model with up to third-nearest neighbors accurately reproduces the band structure for , with an error below . In fact, this is a range of lattice intensities where one would expect the MLWFs to localize strongly around each minimum (that is, a proper tight-binding regime). Then, while the inclusion of provides only a minor correction, the model with just the nearest neighbor tunneling lee2009 is clearly less accurate, reaching the level only for . This may be particularly relevant for the range of parameters of current experiments tarruell2012.

Conclusions. In this article we have demonstrated the power of the maximally localized Wannier functions for composite bands in determining the parameters of tight-binding hamiltonians describing ultracold atoms in optical lattices. The application to a honeycomb structure, directly connected to the graphene physics, allows us to accurately parametrize the optimal tight-binding parameters, providing a thorough analysis of the range of validity of different approaches.

Acknowledgments. 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 the Grant No. FIS2010-19609-C02-00 and the Basque Government through the Grant No. IT-472-10.