Dirac-point engineering and topological phase transitions in honeycomb optical lattices

We study the electronic structure and the phase diagram of non-interacting fermions confined to hexagonal optical lattices. In the first part, we compare the properties of Dirac points arising in the eigenspectrum of either honeycomb or triangular lattices. Numerical results are complemented by analytical equations for weak and strong confinements. In the second part we discuss the phase diagram and the evolution of Dirac points in honeycomb lattices applying a tight-binding description with arbitrary nearest-neighbor hoppings. With increasing asymmetry between the hoppings the Dirac points approach each other. At a critical asymmetry the Dirac points merge to open an energy gap, thus changing the topology of the eigenspectrum. We analyze the trajectory of the Dirac points and study the density of states in the different phases. Manifestations of the phase transition in the temperature dependence of the specific heat and in the structure factor are discussed.


I. INTRODUCTION
The isolation of a single layer of graphene 1 , which is a two-dimensional honeycomb lattice of carbon atoms, has attracted considerable attention, both for its basic interest and because it may pave the way to carbon-based electronics. The low energy electronic properties are described by the massless Dirac equation, which causes many anomalies with respect to semiconductor physics [2][3][4] . Close to half-filling, the band structure near the Fermi energy is given by two degenerate Dirac cones, and the centers of these cones (called Dirac points) are located at the two distinct K-points of the first Brillouin zone. A first experimental validation of the relativistic band structure was the measurement of the anomalous sequencing of the Quantum Hall plateaus 5,6 .
In graphene there is little room to tune system parameters such as the interaction strength or the hopping amplitudes, which prevents a systematic study of its phase diagram. In contrast for cold fermionic atoms confined in honeycomb optical lattices 7 the parameters are highly tunable and many of the theoretically predicted phases might be realized there. Recent theoretical works on honeycomb lattices have studied a topological phase transition in the single particle spectrum which is due either to asymmetric hopping energies 8 , or to a Kekulé distortion of the hoppings 9,10 , or to the application of an ac electric field 11 . Further work on honeycomb optical lattices has dealt with the realization of the anomalous Hall and the Spin Hall effects 12 , the Wigner crystallization of π-bands 13,14 , or with the Quantum Hall effect without Landau levels [15][16][17] . Finally, we note recent theoretical work on the appearance and manifestation of Dirac cones in the band structure of triangular photonic crystals [18][19][20] .
There are two main purposes of this work. The first one is to compare the Dirac points present in the lower lying bands of triangular and honeycomb lattices. Therefore we first discuss in section II the potential landscape of these lattices. Then, in section III we discuss the Dirac cones appearing in both lattice types and derive analytical results within the nearly-free-particle and tight-binding limits. The second main goal of this work is to analyze in detail the topological phase transition in honeycomb lattices with asymmetric hopping energies 8 . In section IV, the evolution of the Dirac points and the resulting phase diagram for the honeycomb lattice as a function of the relative tunnel couplings are studied within the tight-binding limit . The density of states for the different phases is derived in section V. Experimental signatures of the phase transition such as the structure factor and the temperature dependence of the specific heat are explained in section VI. Finally, we discuss the topological structure of the phase transition in section VII. Some technical derivations are left to the Appendices.

II. POTENTIAL LANDSCAPES
We consider the optical lattice that was experimentally realized by Grynberg et al. 7 . It is created by three interfering traveling laser beams. The corresponding wave vectors are in the x − y plane and form equal angles between them, as illustrated in Fig. 1a. For simplicity we assume the polarization of the beams to be perpendicular to the plane, although the confining potential is independent of the polarization for large enough detuning 12  The amplitude of the total electric field (polarized in the z-direction) and the light intensity can then be written (1) As illustrated in Fig. 1a b 1 = q 2 − q 3 plus cyclic permutations. We note that the system is essentially insensitive to phase drifts between the lasers, since these just shift the lattice without changing the potential landscape. We can therefore assume E i > 0 without loss of generality. The optical confining potential V (r) is proportional to the intensity of the laser field and the proportionality factor is positive (negative) if the laser frequencies are blue (red) detuned with respect to the transition frequency of the atoms. Up to a constant the confinement is given by with V 1 ∝ E 2 E 3 plus cyclic permutations and V j > 0 (V j < 0) for blue (red) detuning. The potential V (r) has an underlying triangular Bravais lattice, whose lattice spacing a is determined by the wavelength of the traveling lasers a = 4π/3q = 2λ/3. The unit vectors for real and reciprocal lattices are given by: with b i · a j = 2πδ ij . The unit cells of both reciprocal and real space can be chosen of honeycomb form as shown in Figs 1b and c.
In the case of red detuning, the potential minima lie at the centers of the hexagons thus forming a triangular lattice for all values of V j < 0. If the lasers are blue detuned and their intensities are equal, then the minima of the potential V (r) are located at the vertices of the honeycombs while there is a maximum of the potential at the center of each honeycomb, see Fig. 1b. If the intensities V j become different (i.e. the electric field amplitudes E j begin to be different), the two minima per unit cell approach each other until they merge for strong asymmetries, namely, when |E 1 − E 2 | > E 3 or E 3 > E 1 + E 2 . Thus the potential V (r) (3) has two minima per unit cell for |E 1 − E 2 | < E 3 < E 1 + E 2 . It is evident from Eq. (3) that the potential V (r) has inversion symmetry, which is crucial for the stability of the Dirac points 21 . The position of the minima are given by r = ±r min with: We will show in section IV that the energy dispersion E k calculated within the tight-binding approximation has the same structure as the potential V (r). The analog of the motion of the potential minima for the case of the energy spectrum is the evolution of the Dirac points (8), which will finally merge, causing a topological phase transition from a semimetal to an insulator.

III. SYMMETRIC LATTICES
Because of periodic translational invariance, the crystal momentum is conserved and the eigenspectrum of a particle moving in the potential V (r) of Eq. (3) is most conveniently calculated in Fourier space. The eigenfunctions are thus written as ψ q (r) = m,n c q−Gmn exp(i(q − G mn ) · r), where q is restricted to be in the first Brillouin zone, G mn = mb 1 + nb 2 denotes a general reciprocal lattice vector, and the potential is written as V (r) = mn V mn exp(iG mn · r) where V mn = 1 2 [V 3 (δ m,1 δ n,1 + δ m,−1 δ n,−1 ) + V 1 δ n,0 (δ m,1 + δ m,−1 ) + V 2 δ m,0 (δ n,1 + δ n,−1 )] are the Fourier components of the potential. Due to inversion symmetry of the lattice V mn = V −m,−n .
For each q the Fourier components c q−Gmn form a vector v q , where each entry is labeled by a specific pair of integers (v q ) mn = c q−Gmn . The Schrödinger equation then reduces to a set of linear equations for each v q given by: where E 0 q = 2 q 2 /2m denotes the kinetic energy. Figure 2 shows the band structure for symmetric potentials V 1 = V 2 = V 3 = V for various confinement strengths V . Most importantly, we note that, both for triangular and honeycomb lattices, some pairs of bands touch at the K-points and in the vicinity of those touching points the band structure is given by cones, which are called Dirac cones because of the similarity to the relativistic energy dispersion.
There are however important differences between both cases. For the honeycomb lattice the first pair of Dirac points (located at ±K) arises due to the touching of the two lowest bands. Furthermore, the Fermi surface in the relevant energy interval is exclusively determined by these Dirac cones, leading in particular to a vanishing density of states at complete filling of the lowest band.
By contrast, for a triangular lattice the first pair of Dirac points is caused by the touching between the second and third band. Since these bands are also degenerate at the Γ point, the Dirac points are not isolated but resonate with a continuous band that dominates the density of states in the energy regime of the Dirac cones. Furthermore, the filling factor (number of atoms per unit cell) at which the Fermi energy coincides with the Dirac point is not an integer, also in contrast to the honeycomb lattice.
The formation of the Dirac cone can be derived within the nearly-free-particle approximation, where the potential is treated perturbatively as shown in Appendix A. The distinction between honeycomb and triangular lattices is however better appreciated in the tight-binding regime, where the Hamiltonian is written in terms of Wannier states, each localized at a potential minimum as shown in Appendix B.
In the following we concentrate on the honeycomb lattice within the tight-binding approximation.

IV. DIRAC POINTS IN ASYMMETRIC HONEYCOMB LATTICES
We now discuss the eigenspectrum of a honeycomb lattice with asymmetric hoppings, described by the following tight-binding Hamiltonian: where a k , b k destroy Bloch waves on the two different sublattices and t j > 0 denote nearest-neighbor hopping energies. For symmetric hoppings the standard Hamiltonian (B1) is rediscovered. A possible realization of this Hamiltonian is a deep optical potential V (r) (3) with different values of V j . Due to the inversion symmetry of the confinement potential (3) the onsite energies at the two sublattices are equal and thus can be set to zero in Eq. (7). For asymmetric laser intensities the potentials barriers separating a given potential minimum from the three nearest minima are different. In particular, the corresponding distances between neighboring minima are also different, as follows from Eq. (5). Both effects will lead to asymmetric nearest-neighbor hoppings. We note that a pure shift in the position of the potential minima without a change in the magnitude of the tunnel hopping energies only gives rise to a phase factor in ξ k (7) and thus does not affect the energy spectrum E k = ±|ξ k |.
The energy spectrum of the Hamiltonian (7) is symmetric around zero energy and is given by E k = ±|ξ k |. We note the close analogy between the energy amplitude ξ k and the total electric field E(r). This causes a close similarity between E 2 k = |ξ k | 2 and V (r) ∝ E(r) 2 . As discussed above for the symmetric case t 1 = t 2 = t 3 , valence and conduction band (defined by negative and positive energies, respectively) touch at the two K-points, where they form Dirac cones. Introducing an asymmetry between the tunneling amplitudes the two Dirac points move away from the K-points and approach each other. At a critical asymmetry the Dirac points merge at one of the three inequivalent M-points, and by increasing the asymmetry further a band gap opens. Figure 3 illustrates the phase diagram and the displacement of the Dirac points as a function of the tunnel amplitudes. The central pink region in Fig. 3a embraces the parameter regime of the semimetallic phase, where the asymmetry is small enough for the two inequivalent Dirac points to exist. It is defined by the condition |t 1 − t 2 | < t 3 < t 1 + t 2 , which remains invariant under the permutation of the three subindices. The Dirac points are located at k = ±k D , with: If necessary the wavevector k D as defined above is supposed to be folded back to the first Brillouin zone by adding a uniquely defined reciprocal lattice vector. Figure 3b shows the displacement of the Dirac points following three different lines of the phase diagram, each defined by a fixed ratio t 1 /t 2 while spanned by varying t 3 . We focus first on the red line, which corresponds to the partially symmetric case t 1 = t 2 previously studied in Ref. [8]. For t 3 = 0 valence and conduction band touch along the extended dashed line. Increasing t 3 causes the two inequivalent Dirac cones to approach each other along horizontal lines and for symmetric hoppings t 1 = t 2 = t 3 (center of triangle in Fig. 3a) they are located at the K-points. Finally the Dirac cones merge at the symmetry point In the generic case of t 1 = t 2 = t 3 one finds the motion shown by the blue and green lines. We now discuss the motion for increasing t 3 . For t 3 < |t 1 − t 2 | the system is an insulator and there is a single band minimum located at M 1 , if t 1 > t 2 (blue line in region A 1 of Fig. 3a), or M 2 , if t 2 > t 1 (green line in region A 2 of the same figure). While in regions A 1 or A 2 , the minimum stays pinned at the fixed points M 1 or M 2 and motion along the blue or green line only causes the corresponding gap to decrease until it becomes zero when the central pink triangle is touched (phase transition). For |t 1 − t 2 | < t 3 < t 1 + t 2 , in both cases (blue and green) the system enters the central triangle of Fig. 3a defining the semimetallic phase, where the band minimum at M 1 or M 2 separates into two Dirac points. The position of the Dirac points are related to each other by the inversion symmetry around the symmetry points Γ and M 1 , M 2 , M 3 . Finally at t 3 = t 1 + t 2 , where the lines leave the pink triangle in order to enter region A 3 , the Dirac points merge at M 3 and the system undergoes a phase transition from a semimetal to an insulator. Further increase of t 3 only causes the gap at M 3 to increase.

V. DENSITY OF STATES
We now discuss the density of states, which is defined by ρ(E) = A −1 k δ(E − E k ), where A = N A u is the area, N the number of unit cells and A u = a 2 √ 3/2 the size of the unit cell. Characteristic energies are given by the eigenenergies at the symmetry points E Γ is the energy maximum and thus sets the bandwidth. The density has a step at E = ±E Γ .
Without loss of generality we focus on the case t 3 > t 1 , t 2 , so that the energy always has a saddle point at the symmetry points M 2 , M 3 , which results in logarithmic van Hove singularities at E M2 , E M3 . Whether M 3 is a saddle point or a minimum depends on the tunnel amplitudes. For t 3 < t 1 + t 2 the energy has a saddle point at M 3 and there are two Dirac points in the first Brillouin zone. For symmetric hoppings t 3 = t 1 = t 2 the band structure close to zero energy is given by rotationally symmetric cones centered at the Dirac points. As the hopping energies become different from each other, those cones are stretched and the Dirac points move, but still the density of states vanishes linearly close to zero energy (E ≪ E M1 , E M2 ). In particular for t 1 = t 2 > t 3 /2 the low energy spectrum is given by E kD+q = v 2 x q 2 x + v 2 y q 2 y , v x = a 4t 2 1 − t 2 3 /2, v y = at 3 √ 3/2, resulting in a density of states ρ(E) = g v E/2πv x v y , where the degeneracy factor g v = 2 takes into account the two Dirac points. We note that the velocities v x , v y are different if the hoppings are different. While v y increases with increasing t 3 , v x decreases and cancels at the phase transition t 3 = 2t 1 = 2t 2 . Examples of the density of states in the semimetallic phase are given in Fig. 4a.
In the opposite limit t 3 > t 1 + t 2 the energy has its minimum at M 3 and the system is gapped by 2∆ where Correspondingly the density of states is zero for E < ∆ and has a finite step at E = ∆. For t 1 = t 2 < t 3 /2 the low energy spectrum is given by E kD +q = (∆ 2 + v 2 x q 2 x + v 2 y q 2 y ) 1/2 , v x = a ∆t 1 /2, v y = a 3t 1 t 3 /2, leading to a density of states ρ(E) = θ(E − ∆)E/2πv x v y , where θ(x) denotes the Heaviside step function. This situation is depicted in Fig. 4c.
Right at t 3 = t 1 + t 2 the two Dirac cones merge at M 3 with E M3 = 0, so that M 3 is the minimum of the conduction band, but there is still no band gap. In the particular case t 1 = t 2 = t 3 /2, the density of states can be mostly derived analytically: with ǫ = E/2t 3 denoting the energy in units of the bandwidth. This situation is depicted in Fig. 4b. We note that E M3+q ≃ t 3 [3(aq y /2) 2 + (a/2) 4 (q 4 x − 6q 2 x q 2 y − 3q 4 y )/4] 1/2 , which shows that the confinement in the x-direction (along which the Dirac points have merged) is much weaker than in the y-direction. Equation (10) is valid for the whole energy range 0 ≤ ǫ ≤ 1. An expansion around E = 0 results in: Since the square root has a diverging slope at the origin, it can be viewed as a transient behavior between the linear density of states of the semimetallic phase and the step found for the insulating phase. We note that for ǫ = 0.5, which is the energy of the two saddle points E M1 = E M2 , the integrand of Eq. (11) contains the factor (1 − x) −1 which results in a logarithmic divergence in the density of states.

VI. SPECIFIC HEAT AND STRUCTURE FACTOR
Finally we discuss two experimental signatures of the phase transition expected when fermions populate the lattice at half filling, namely, the structure factor accessible by Bragg spectroscopy 8 and the temperature dependence of the specific heat 22 . The specific heat is defined by c V ≡ ∂u/∂T = ∂ ∂T dEρ(E)Ef (E), where f (E) denotes the Fermi distribution, u the energy density and ρ(E) the density of states as calculated in the previous section. We limit ourselves to the case of half-filling, where the chemical potential µ = 0 for all temperatures because of the symmetry of the energy around E = 0 inherent to the employed tight-binding approximation.
The resulting equations are: where we have used ρ(E) = ρ(−E). This results in a different temperature dependence of the specific heat for the different phases: where ζ denotes the Riemann Zeta Function. Interaction induced corrections to the specific heat in the semimetallic phase have recently been studied 23 . As pointed out in Ref. 8 another way to experimentally visualize the phase transition is to measure the structure factor by Bragg spectroscopy. We note that the structure factor is up to a constant given by the imaginary part of the dynamical susceptibility which is well known for graphene [24][25][26] . In fact within the semimetallic phase, the results for the susceptibility of graphene can be easily adapted to include asymmetric hoppings. Noting that E q ≃ v 2 x q 2 x + v 2 y q 2 y 1/2 = κ, with κ x = q x v x , κ y = q y v y , κ = (κ 2 x + κ 2 y ) 1/2 , one can reuse all formulae of the dynamic susceptibility of arbitrary doping by setting v F = 1 and replacing q → κ.
As an illustration, we state the result for half-filling: Here g = 2 is due to the two Dirac points (valley degeneracy). For t 1 = t 2 > t 3 /2 the velocities are given by v x = a 4t 2 1 − t 2 3 1/2 /2, v y = at 3 √ 3/2. We note that in Ref. 8 only the energy spectrum was considered while the overlap factor f λ ′ λ (k, q) resulting from the spinor eigenfunctions was ignored. Therefore the ω 2 term present in the numerator of Eq. (8) of Ref. 8 is absent in Eq. (14) above. We note that f λ ′ λ (k, q) = |u * k+q,λ ′ u kλ | 2 describes how easily the spinor wavefunction u kλ (λ = ∓1 labels valence and conduction band) can be scattered in to u k+q,λ ′ and is responsible for example for the absence of backscattering and the Klein paradox in graphene. 27 The general form of the spinor wavefunction for the tight-binding Hamiltonian (7), is u k λ = 2 −1/2 (−λ ξ k /|ξ k |, 1).
For the gaped case t 1 = t 2 < t 3 /2, an analytical formula for the structure factor at half-filling and for energies close While the step function originates from the single particle eigenspectrum, the whole q-dependence arises from the spinor character of the wavefunctions captured in the factor f λ ′ λ (k, q) described above.

VII. TOPOLOGICAL PROPERTIES
Interestingly, a gap only opens after the two Dirac points have merged. Mathematically it is straightforward to show that an energy minimum of the conduction band that is not located at one of the symmetry points M 1 , M 2 , M 3 , Γ is automatically a Dirac point. To show this, we note that a minimum of the conduction band is a root of Since a 1 , a 2 are linear independent, v 1 , v 2 are also linear independent, unless sin ((a 2 − a 1 ) · k)) = 0, a condition is only fulfilled at the symmetry points M 1 ,M 2 ,M 3 ,Γ. This analysis shows that a minimum which is not located at a symmetry point, is necessarily a Dirac point, since then both the imaginary and the real parts of η k vanish and therefore the energies are E ∓,k = ∓|η k | = 0. Furthermore, these minima always occur in pairs due to time reversal symmetry (inversion symmetry around the Γ-point).
The deeper reason for the stability of the Dirac points is the momentum space topology of the Bloch wavefunctions 21,28 . This is nicely illustrated by noting that the Berry phase, defined as 29 φ B = i S dk u k |∇ k |u k with S denoting a closed path in reciprocal space, takes values ±π, 0 depending on whether S encloses one Dirac point, or the other, or both. Thus each Dirac point is a source of a ±π delta-function flux of Berry curvature, and this flux remains invariant under perturbations that preserve space and time inversion symmetry 21,28 . Only if the two Dirac points merge the Berry curvature vanishes within the whole first Brillouin zone and a gap can open.
We note that the Dirac points are unstable against perturbations that break spatial invariance. Examples are substrate-induced gap opening in graphene due to breaking of sublattice symmetry 30 or honeycomb lattices where nearest-neighbor hoppings are periodically modified to form a Kekulé pattern 9 .

VIII. SUMMARY AND CONCLUSIONS
We have studied the electronic structure and the phase diagram of non-interacting fermions confined to hexagonal optical lattices. In the first part of the paper, we have analyzed the appearance of Dirac points in the band structure of fermionic atoms populating triangular and honeycomb optical lattices for different confinement strengths. Numerical results were complemented by analytical equations both for strong and weak potentials. The Dirac points arising in honeycomb and triangular lattices differ in several important aspects. While in honeycomb lattice the Dirac cone is isolated, it resonates with a continuum band for a triangular lattice. Furthermore, in honeycomb lattices the Fermienergy coincides with the Dirac point for integer filling, while a non-integer filling is needed in the case of a triangular lattice. For the honeycomb lattice the first pair of Dirac points arises from the touching of the two lowest bands, which within a tight-binding description stem from the ground state of the two potential minima per unit cell. By contrast, for a triangular lattice the first pair of Dirac points occurs at the touching of the second and third band, which within a tight-binding description are formed by the doubly degenerate first excited state of the single potential minimum per unit cell.
In the second part of this work, we have focused on the phase diagram and the evolution of the Dirac points in honeycomb lattices adopting a tight-binding description with asymmetric nearest-neighbor hoppings. The semimetallic phase is not only realized for strictly symmetric hoppings but rather for an extended region of the hopping parameter space. This stability of the semimetallic phase is based on topology. With increasing asymmetry the Dirac points approach each other along continuous trajectories and at a critical asymmetry they finally merge and an energy gap opens. We derive analytic formulae for the trajectories of the Dirac points as well as for the density of states. Right at the phase transition the density of states increases like √ E, which interpolates between the linear dependence of the semimetallic phase and the step-like increase of the insulating phase. We also show how the phase transition becomes manifest in a change of the temperature dependence of the specific heat as well as in the properties of the structure factor.
An interesting follow-up of this work concerns the influence of interactions on the phase diagram, since optical lattices permit a systematic control of the effective interaction strength by using Feshbach resonances or modifying the lattice properties 31 .

APPENDIX A: DIRAC CONES IN WEAK SYMMETRIC LATTICES
Here we apply the nearly-free-particle approximation to derive the eigenspectrum determined by Eq. (6) with V 1 = V 2 = V 3 = V close to the K-point. In absence of any potential the dispersion relation at the K-point is three-fold degenerate, where the three states can be labeled by the k-values k ∈ {K, K − G 10 , K − G 11 }. Treating the lattice within lowest order perturbation theory, the eigenvalue problem at k = K + q has then the following form: The solutions of the above eigenproblem for q = 0 are given by E = E 0 K + V with eigenvector 3 −1/2 (1, 1, 1) and E = E 0 K − V /2, which is two fold degenerate with a possible set of orthonormal eigenvectors given by: v 1 = 2 −1/2 (−1, 1, 0) and v 2 = 6 −1/2 (−1, −1, 2). We note that the remaining two-fold degeneracy is protected by topology and thus persists beyond a perturbative description. This degeneracy occurs in the ground state for V > 0 (honeycomb lattice) and in the first excited state for V < 0 (triangular lattice).
The remaining two-fold degenaracy of v 1 , v 2 disappears for finite q. The effective two-band Hamiltonian spanned by v 1 ,v 2 is given by: with v F = 2π /3am. The eigenenergies are given by We note that the effective Hamiltonian can be transformed by a unitary transformation toH eff = E 0 ½ + v F (q x σ x + q y σ y ), which is identical (except for the value of v F ) to the low-energy Hamiltonian obtained in the tight-binding approximation . For the honeycomb lattice there are two degenerate minima per unit cell and for sufficiently deep potential wells [i.e. for V ≪ 2 2ma 2 , with m the atomic mass and V the confinement amplitudes of Eq. (3)] the basis for the two lowest bands can be restricted to the ground state of each minimum. In the symmetric case (equal laser intensities) the potential around each minimum is rotationally symmetric and harmonic and the minima form a honeycomb lattice. Thus the nearest-neighbor hoppings t are all the same.
Here a(r) destroys a particle at the A-site at position r and the d j denotes the vectors connecting an A-site with its neighboring B-sites, see Fig. 1b. Fourier transformation is given by a(r i ) = N −1/2 k exp(ik · r i )a k , where N is the number of unit cells. The eigenenergies of the above Hamiltonian are given by: E k = ±|ξ k | = ±t[3 + 2 cos(k y a) + 4 cos( √ 3ak x /2) cos(ak y /2)] 1/2 . The energy spectrum is symmetric around E = 0, where the Fermi surface only consists of singular (Dirac) points located at the K-points. Close to the K-point one can approximate ξ K+q ≃ −v F (q x − iq y ) with v F = √ 3at/2, so that the Hamiltonian close to the K-point can be written as H eff = v F (q x σ x + q y σ y ) with eigenenergies E q = ± v F q.

Symmetric triangular lattice
The triangular lattice has one minimum per unit cell. For a symmetric lattice the potential around the minimum is rotationally symmetric and harmonic. The ground state of such a minimum is non-degenerate while the first excited state is two fold degenerate and can be represented by p x (x, y) = ψ 1 (x)ψ 0 (y) and p y (x, y) = ψ 0 (x)ψ 1 (y), where x, y are measured with respect of the center of the minimum and ψ 0 (x), ψ 1 (x) denote the ground state and the first excited state of the one-dimensional harmonic oscillator.
The band arising from the ground state is described by the following tight-binding Hamiltonian: [cos(k · a 1 ) + cos(k · a 2 ) + cos(k · (a 1 − a 2 ))] (B2) The hopping between the p x and p y orbitals is no longer isotropic 14 . As visualized in Fig. 5b, we introduce the hopping amplitudes t σ and t π depending on the relative orientation of the orbitals. Non-parallel hoppings (e.g. from p x to p y ) are excluded by the symmetry of the wavefunctions.
The tight-binding Hamiltonian describing the two bands arising from the two-fold degenerate first excited state is