Dirac points emerging from flat bands in Lieb-kagom\'e lattices

TL;DR

This paper introduces a model interpolating between Lieb and kagomé lattices, revealing how flat bands generate Dirac points with topological properties, advancing understanding of multiband contact points.

Contribution

It presents a continuous interpolation model between Lieb and kagomé lattices, analyzing the topological and spectral evolution of Dirac points and flat bands.

Findings

Flat band distortion creates Dirac points in the spectrum.

Dirac points evolve and merge, altering their topological winding numbers.

The model offers a systematic way to characterize two-band contact points.

Abstract

The energy spectra for the tight-binding models on the Lieb and kagom\'e lattices both exhibit a flat band. We present a model which continuously interpolates between these two limits. The flat band located in the middle of the three-band spectrum for the Lieb lattice is distorted, generating two pairs of Dirac points. While the upper pair evolves into graphene-like Dirac cones in the kagom\'e limit, the low energy pair evolves until it merges producing the band-bottom flat band. The topological characterization of the Dirac points is achieved by projecting the Hamiltonian on the two relevant bands in order to obtain an effective Dirac Hamiltonian. The low energy pair of Dirac points is particularly interesting in this respect: when they emerge, they have opposite winding numbers, but as they merge, they have the same winding number. This apparent paradox is due to a continuous rotation of their states in pseudo-spin space, characterized by a winding vector. This simple, but quite rich model, suggests a way to a systematic characterization of two-band contact points in multiband systems.

Equations

$$\displaystyle H({\bf k})=\left(\begin{array}[]{ccc}0&1+e^{ik_{x}}&t(1+e^{i(k_{x}+k_{y})})\\ \ldots&0&1+e^{ik_{y}}\\ \ldots&\ldots&0\\ \end{array}\right),$$

$$\mathcal{H}(\mathbf{q})=h_{\mu }(\mathbf{q}){\sigma }_{\mu }+h_{\nu }(\mathbf{q}){\sigma }_{\nu }\equiv \vec{h}(\mathbf{q})\cdot \vec{\sigma }$$

$${\mathcal{H}}_{+-}$$

$${\mathcal{H}}_{++}$$

$$\vec{w}=\frac{1}{2\pi }\int \vec{n}\times d\vec{n}$$

$$\mathcal{H}=v_x{q}_x{\sigma }_{\mu }+v_y{q}_y{\sigma }_{\nu },$$

$$\vec{w}=\text{sgn}(v_x{v}_y){\vec{u}}_{\mu }\times {\vec{u}}_{\nu }.$$

$$\xi {\mathbf{k}}_a=\xi (2{\varphi }_a,-2{\varphi }_a)\text{ and }{E}_a=-2\cos {\varphi }_a,$$

$$\cos {\varphi }_a=t.$$

$${\tilde{H}}_a(k_x,-k_x)=\text{diag}(E_{+},E_{-},E_3)$$

$$\displaystyle\tilde{H}_{M}({\bf q})=\left(\begin{array}[]{c|cc}2t&iq_{\perp}&i\sqrt{2}tq_{\parallel}\\ \hline\cr\ldots&0&-iq_{\parallel}+\frac{q_{\perp}q_{\parallel}}{\sqrt{2}}\\ \ldots&\ldots&-2t\\ \end{array}\right)$$

$${\mathcal{H}}_M(\mathbf{q})$$

$$\frac{q_{\parallel }{q}_{\perp }}{2\sqrt{2}}{\sigma }_x.$$

$$\displaystyle\mathcal{H}_{a}({\bf q})\simeq\bigl{(}-2t-\xi\frac{q_{\perp}}{\sqrt{2}}\bigr{)}\sigma_{0}+q_{\parallel}(\sigma_{y}+\xi t\sigma_{x})-\xi\frac{q_{\perp}}{\sqrt{2}}\sigma_{z}.$$

$$\displaystyle\tilde{H}_{\Gamma}({\bf q})=\left(\begin{array}[]{c|cc}4&-iq_{\perp}&i\sqrt{3}q_{\parallel}\\ \hline\cr\ldots&-2-\frac{2\epsilon}{3}+\frac{q_{\perp}^{2}}{3}&\displaystyle-\frac{q_{\parallel}q_{\perp}}{\sqrt{3}}\\ \ldots&\ldots&\displaystyle-2+2\epsilon+q_{\parallel}^{2}\end{array}\right)$$

$${\mathcal{H}}_{\Gamma }(\mathbf{q})$$

$$\left(\frac{q_{\perp }^2-q_{\parallel }^2}{4\sqrt{6}}-\frac{4ϵ}{3\sqrt{6}}\right)q_{\parallel }{\sigma }_y.$$

$${\mathcal{H}}_a(\mathbf{q})$$

$$\displaystyle\tilde{H}_{a}({\bf q})=\left(\begin{array}[]{c|cc}4t&\mathcal{O}({\bf q})&\mathcal{O}({\bf q})\\ \hline\cr\ldots&E_{a}-\frac{2\sqrt{2}}{3}\sin\phi_{a}q_{\perp}&\sqrt{\frac{2}{3}}\sin\phi_{a}e^{-i\phi_{a}}q_{\parallel}\\ \ldots&\ldots&E_{a}\end{array}\right),$$

$$\displaystyle\mathcal{H}_{a}({\bf q})\simeq\bigl{(}E_{a}-v^{0}_{\perp}q_{\perp}\bigr{)}\sigma_{0}+v_{\parallel}q_{\parallel}\sigma_{\phi_{a}}-v_{\perp}q_{\perp}\sigma_{z}\ ,$$

$$v_{\parallel }=\sqrt{\frac{2}{3}}\sin {\varphi }_a,v_{\perp }=v_{\perp }^0=\frac{v_{\parallel }}{\sqrt{3}}.$$

$${\vec{w}}_{\xi }$$

$$\xi {\mathbf{k}}_d=\xi (2{\varphi }_d,2{\varphi }_d)\text{ and }{E}_d=2\cos {\varphi }_d$$

$$\cos {\varphi }_d=\frac{\sqrt{1+8t^2}-1}{4t}.$$

$${\tilde{H}}_d(k_x,k_x)=\text{diag}(E_1,E_{+},E_{-})$$

$$\displaystyle\tilde{H}_{M}({\bf q})=\left(\begin{array}[]{cc|c}2t&iq_{\perp}+\frac{q_{\parallel}^{2}+q_{\perp}^{2}}{2\sqrt{2}}&i\sqrt{2}tq_{\parallel}\\ \ldots&0&-iq_{\parallel}\\ \hline\cr\ldots&\ldots&-2t\\ \end{array}\right),$$

$${\mathcal{H}}_M(\mathbf{q})$$

$$\frac{q_{\perp }^2}{2\sqrt{2}}{\sigma }_x.$$

$${\mathcal{H}}_d\simeq (2t-\frac{\xi }{\sqrt{2}}{q}_{\parallel }){\sigma }_0-q_{\perp }{\sigma }_y+\frac{\xi }{\sqrt{2}}{q}_{\parallel }{\sigma }_z.$$

$${\mathcal{H}}_d(\mathbf{q})=\left(E_d-\xi v_{\parallel }^0{q}_{\parallel }\right){\sigma }_0+\xi v_{\parallel }{q}_{\parallel }{\sigma }_z-v_{\perp }{q}_{\perp }{\sigma }_y,$$

$$v_{\parallel }^0=\frac{\sqrt{2}\sin (3{\varphi }_d)}{3\cos (2{\varphi }_d)}$$

Full text

I Introduction

Multiband systems with flat band are of current interest thanks to their experimental realizations in a variety of systems Jo12 ; Guzman14 ; Taie15 ; Vicencio15 ; Muk14 ; Baboux16 ; Slot17 . The list includes cold atoms with optical lattices Jo12 ; Taie15 , photonic Guzman14 ; Vicencio15 ; Muk14 and polaritonic Baboux16 systems. Typically, the systems are engineered to faithfully realize exotic tight-binding models with multiple orbitals or atoms per unit cell, as well as complex lattice geometries. In these artificial systems, direct imaging of localized states Vicencio15 ; Slot17 ; Milicevic19 or the study of tunable interaction-induced effects Ozawa17 are few examples which show approaches hitherto not easily realizable in conventional condensed matter systems.

In this work, we study a connection between two two-dimensional tight-binding models - on the Lieb lattice and on the kagomé lattice (also known as the trihexagonal tiling) - both featuring two dispersive bands and a flat band. On the one hand, the Lieb model Dagotto1986 ; Apaja2010 ; Goldman2011 ; Tsai2015 has a band contact between three bands among which a topological flat band Aoki1996 at zero energy (see Fig. 1a, the tight-binding model is defined in Fig. 2a). On the other hand, the kagomé model Ohgushi00 ; Xiao2003 , chosen here with a square symmetry, features a quadratic band contact point Chong2008 ; Sun2009 ; Dora between its lower bands: a dispersive band and a flat band resulting from destructive interferences Aoki1996 (see Fig. 1d). The real-space model is defined in Fig. 2b. Apart from the flat band, the kagomé energy spectrum is very similar to that of graphene with a pair of Dirac points between the upper bands (see Fig. 1d). Both models have found physical realizations in the aforementioned systems Jo12 ; Guzman14 ; Taie15 ; Vicencio15 ; Muk14 ; Baboux16 ; Slot17 . There are also solid state systems expected to behave as a kagomé lattice for itinerant electrons, see e.g. Refs. Mazin2014 ; Ye2018 ; Leykam2018 .

By interpolating between the two lattice models, we define the Lieb-kagomé model (see Fig. 2c). It displays a smooth crossover of the flat band from the middle of the spectrum to the lowest energy. This model was already briefly introduced and its energy spectrum discussed by Asano and Hotta (see Fig. 7 in Ref. Asano2011 ). Of particular interest here is the continuous deformation of the flat band without gap opening when deviating from the two limits. This deformation is accompanied by the creation of pairs of Dirac points – i.e. linear band contact points between two bands.

Very recently, W. Jiang et al. WeiJiang2019 have considered the Lieb-kagomé model in the presence of a time-reversal breaking term (related to intrinsic spin-orbit coupling) that gaps the energy spectrum. The topology of the isolated bands was then characterized by computing the corresponding Chern numbers, which obviously depends on the way the gap is opened. In contrast, in the present article, we do not gap the energy spectrum and instead characterize the band contact points as topological defects by computing the associated charge.

On the one hand, slightly deforming the Lieb lattice, two pairs of Dirac points are generated from the three band contact at the M$=(\pi,\pi)$ point of the first Brillouin zone (BZ) (see Fig. 1a and b). One pair is between the lower two bands and will evolve in the flat band of the kagomé limit. The other pair (between the upper two bands) will separate and turn into the Dirac points of the kagomé lattice. As we will see below, near the M point, each pair has opposite winding numbers ($+-$). On the other hand, deforming the kagomé lattice, a single pair of Dirac points emerges from the quadratic-flat band contact at the $\Gamma$ point between the two lower bands (see Fig. 1f and e). The latter has a $+2$ winding number that produces a pair of Dirac points with identical winding numbers ($++$).

The evolution of the lower two bands then raises the question of how the two distinct fusion scenarios of Dirac points can be smoothly connected, namely, a $(+-)$ pair of Dirac points emerging at M evolving into a $(++)$ pair merging at $\Gamma$. Actually, a similar phenomenon has been studied by us in a two-band model, the staggered Mielke lattice Montambaux2018 . While the integer-valued winding number characterizes the circulation of the Dirac spinor on the great circle of the Bloch sphere, the great circle itself is not fixed generally - it evolves according to the Hamiltonian parameters. We have shown that, in order to describe properly a pair of Dirac points and their merging, the orientation of the great circle around each Dirac point has to be specified by a vector perpendicular to the plane of the great circle, that we named the winding vector. Then during the smooth evolution between the two merging scenarios, the relative direction of winding vectors continuously evolves from a parallel to an antiparallel situation.

In the three-band Lieb-kagomé model, the mechanism of distortion and restoration of the flat band implies the generation of Dirac pairs which follows different merging scenarios. The winding vector description a priori applies only in a two-band system. One of our main task is thus to find a consistent basis to represent the Lieb-kagomé Hamiltonian, local in ${\bf k}$-space, from which effective two-band descriptions for the contact points can be derived.

There is a deep connection with the work by Ahn et al. on the topological characterization of real band structures Ahn2019 . This is best seen by working in a basis in which the Bloch Hamiltonian of the Lieb-kagomé model is real. We detail this comparison in Appendix D.

The paper is organized as follow. In Sect. II, we define the model and give a brief summary of the universal Hamiltonians of Dirac points merging/emergence. In Sect. III and IV, we construct the smooth basis with which effective two-band Hamiltonians can be derived. Then we determine the evolution of the winding vectors. In Sect. V, we give a global picture of the evolving winding vector and follow with conclusions. In the Appendices A and B, we give details on several computations. We also generalize the Lieb-kagomé model in Appendix C. In Appendix D, we discuss the Lieb-kagomé model in an alternative representation in which the Bloch Hamiltonian is real and make contact with Ref. Ahn2019 .

II Lieb-kagomé lattice and contact points

II.1 Tight-binding model

The Lieb-kagomé tight-binding model is schematically shown in Fig. 2c, where the full bonds indicate hopping amplitudes of fixed magnitude (which we take as energy unit) and the dashed bonds correspond to variable hopping of strength $t$. The unit cell contains three orbitals at sites A, B and C, which results in three energy bands. In the present work, we are interested in the evolution of the band structure for the parameter range $0\leq t\leq 1$. For any $t$, this model has time-reversal symmetry, inversion symmetry (with center either on the $B$ sites or in the middle of the unit cell shown in Fig. 2c) and it also possesses mirror symmetries with respect to reflections through the diagonal (d) and antidiagonal (a) lines passing through B-sites of the unit cell. Two well-known limits of the above model are:

(1) the Lieb lattice $t=0$. It is a regular square lattice with B-sites forming the Bravais lattice and additional sites (A- and C-sites) located at the center in between all the nearest-neighbor B-sites (Fig. 2a). In addition to the above symmetries, it is also bipartite (i.e. it has a chiral sublattice symmetry) and it possesses mirror symmetries along the $x$ axis and along the $y$ axis. Since it is bipartite, its spectrum respects the energy inversion $E\rightarrow-E$ symmetry and hence the middle band is necessarily a zero-energy flat band. According to the classification of Aoki et al., this flat band is of topological origin Aoki1996 . For example, its energy is not affected by a perpendicular magnetic field Aoki1996 . Another example of a topological flat band is found in the dice lattice Sutherland1986 .

(2) the kagomé lattice $t=1$. Here it is deformed from the usual kagomé to have a square symmetry but it has the same lattice connectivity. The kagomé lattice belongs to a family of frustrated tight-binding models constructed on line graphs of bipartite parent lattices, of which the Mielke (checkerboard) lattice is another instance Mielke ; Bergman ; Chalker . These models display a flat band at the minimum (or maximum) energy which is in contact with a dispersive band. The electron localization leading to the flat band is due to a destructive interference, which is affected by a perpendicular magnetic field Aoki1996 .

We choose to write the Bloch Hamiltonian in the form $H({\bf k})=e^{-i{\bf k}\cdot\mathbf{R}}He^{i{\bf k}\cdot\mathbf{R}}$ where $\mathbf{R}$ are the positions of the Bravais lattice sites (see Appendix A and D.2). In the momentum space basis $\{|{\bf k},A\rangle,|{\bf k},B\rangle,|{\bf k},C\rangle\}$ the Bloch Hamiltonian can be represented as

[TABLE]

where we have chosen the lattice spacing $a=1$ (see Fig 2a).

In the following, we will use the following convention for Bloch Hamiltonians. Three-band Hamiltonians in the site basis will be denoted $H$, three-band Hamiltonians in another basis $\tilde{H}$, and two-band effective Hamiltonians $\mathcal{H}$.

II.2 Energy spectrum

The energy spectrum for $0\leq t\leq 1$ is shown in Fig. 1, which displays three energy bands and remains gapless. The band touching points and their characteristics depend on $t$.

For the Lieb lattice ($t=0$) there is a triple band crossing at the M point in momentum space ($k_{x}=k_{y}=\pi$ and their equivalents, Fig. 3a). It shows a linear contact point with a flat band in the middle, which is characteristic of a pseudospin-1 Goldman2011 . In deviating from the Lieb limit ($t>0$), the flat band is distorted with the triple crossing point splitting into two pairs of Dirac points and particle-hole symmetry of the energy spectrum is lost. They are formed respectively in between the lower and middle bands, and the middle and the upper bands, in two orthogonal directions in momentum space, see Fig. 4.

Moving away from the Lieb limit $t>0$, the distortion of the middle band gets more pronounced and the lower Dirac pair moves in the antidiagonal direction $(k_{x}=-k_{y})$ towards the $\Gamma$ point of the BZ ($k_{x}=k_{y}=0$) (red arrows in Fig. 3a). These Dirac cones are critically tilted such as to exhibit a zero-velocity line as seen on Figs. 1c) and 5. In approaching the kagomé limit $t\rightarrow 1$, the same Dirac points merge at the $\Gamma$ point forming a single contact point between a quadratic band and a flat band, thus restoring a flat lowest energy band. The upper Dirac pair, on the other hand, moves along the diagonal direction ($k_{x}=k_{y}$, blue arrows in Fig. 3a) towards a fixed position $\xi{\bf K}=\xi(2\pi/3,2\pi/3)$, where $\xi=\pm 1$ is a valley index, while maintaining the Dirac cone structure throughout. The tilt of these Dirac cones diminishes with increasing $t$ and vanishes at $t=1$ [see Fig. 1(c) and (d)].

In summary, the crossover between flat bands happens by producing pairs $({\bf k}_{D},-{\bf k}_{D})$ of Dirac points. The lower pair evolves from the M and the $\Gamma$ point. Its merging/emergence at time-reversal invariant momenta (TRIM) resembles a situation we already encountered in the two-band model on a staggered Mielke lattice Montambaux2018 . There, we provided a complete scenario of a $(+-)$ going into a $(++)$ Dirac pair. In the following subsection we briefly summarize the main results of Ref. Montambaux2018 .

II.3 Characterization of contact points: $2\times 2$ effective Hamiltonians

A contact point between two bands is characterized by an effective $2\times 2$ Hamiltonian. Around a single contact point $\mathbf{k}_{D}$, the local two-band Bloch Hamiltonian takes on a general form

[TABLE]

involving only two among three Pauli matrices $\sigma_{\mu}=\vec{u}_{\mu}(\mathbf{k}_{D})\cdot\vec{\sigma}$ and $\sigma_{\nu}=\vec{u}_{\nu}(\mathbf{k}_{D})\cdot\vec{\sigma}$ with $|\vec{u}_{\mu,\nu}|=1$, $\vec{u}_{\mu}\cdot\vec{u}_{\nu}=0$ (e.g., $\mu,\nu=x,y$) and $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$. The winding number is the number of times the effective magnetic field $\vec{h}(\mathbf{q})$ winds on a great circle of the Bloch sphere in the $(\vec{u}_{\mu},\vec{u}_{\nu})$ plane when ${\bf q}$ encircles the contact point once in the trigonometric direction. For example, for a Dirac point one can have $\vec{h}(\mathbf{q})=(v_{x}q_{x},v_{y}q_{y},0)$ and the winding number is $\text{sign}(v_{x}v_{y})$.

Indeed, the stability of two-band contact points in 2D requires a symmetry protection (such as a chiral symmetry or inversion and time-reversal symmetries). The associated topological charge is the $\mathbb{Z}$ winding number that classifies maps from a loop $S^{1}$ encircling the contact point in reciprocal space to a great circle $S^{1}$ of the Bloch sphere associated to the $2\times 2$ Bloch Hamiltonian $\mathcal{H}(\mathbf{q})$. In the absence of such a symmetry, band contact points are unstable defects in 2D Volovik .

To describe a pair of Dirac points when they come close together in the vicinity of a time-reversal invariant momentum (TRIM), two merging scenarios are possible corresponding to the following two universal Hamiltonians Montambaux2009 ; Gail2012

[TABLE]

where $\sigma_{\mu}$, $\sigma_{\nu}$ are two Pauli matrices and $\Delta$, $c$, $m$ and $m^{\prime}$ are parameters. Actually, because of time-reversal symmetry, there is not much choice in the Pauli matrices: in the $(+-)$ case $\sigma_{\nu}$ should be $\sigma_{y}$ and in the $(++)$ case, the $\sigma_{y}$ matrix should not appear in the Hamiltonian.

Making an expansion around each Dirac point, the first (resp. second) Hamiltonian describes two Dirac points with opposite (resp. equal) winding numbers located at position $(q_{x},q_{y})=(\xi\sqrt{-2m\Delta},0)$ if $\Delta<0$. If $\Delta>0$, Dirac points only exist in the second case and are located at $(q_{x},q_{y})=(0,\xi\sqrt{2m\Delta})$. Here the winding number $\pm 1$ associated with a Dirac point is actually defined with respect to the winding plane $(\sigma_{\mu},\sigma_{\nu})$, which can generally be parameter dependent. This is indeed the case in situations encountered in Ref. Montambaux2018 and in this paper: when varying a parameter of the Hamiltonian (here the parameter $t\in[0,1]$), the orientation of the winding plane (defined by a great circle on the Bloch sphere) varies continuously. We shall see here that a pair of Dirac points emerges following the ${\cal H}_{+-}$ scenario with $(\sigma_{\mu},\sigma_{\nu})=(\sigma_{z},\sigma_{y})$, and merges following the ${\cal H}_{++}$ scenario with $(\sigma_{\mu},\sigma_{\nu})=(\sigma_{z},\sigma_{x})$. In order to describe the continuous evolution of the winding planes attached to each Dirac point, we recently introduced the notion of a winding vector Montambaux2018

[TABLE]

with the normalized pseudospin magnetic field $\vec{n}(\mathbf{k})=\vec{h}/|\vec{h}|$. $\vec{w}$ is the unit vector normal to the oriented winding plane of $\vec{n}$. The winding vector is a parameter-dependent vectorial quantity, supplemental to the more familiar winding number in describing a band crossing. This notion becomes essential when describing a pair of Dirac points in close vicinity since the relative orientation of the winding plane can have physical consequences. For a Hamiltonian of the form

[TABLE]

with $\sigma_{\mu}=\vec{u}_{\mu}.\vec{\sigma}$ and $\sigma_{\nu}=\vec{u}_{\nu}.\vec{\sigma}$, we simply have

[TABLE]

In the above equation, it is implicitly assumed that $\{\vec{u}_{x},\vec{u}_{y},\vec{u}_{z}\}$ is an oriented basis.

Finally we note that a $\sigma_{0}$ term with momentum dependence is important for the flatness of bands in the energy spectrum but not important for wavefunctions and the topological characterization of a band contact point.

II.4 Reduction from three to two-band effective Hamiltonians

In order to characterize the distortion of the flat band and the evolution of Dirac points in the three-band system, we introduce local $2\times 2$ effective Hamiltonians. In order to do this, we proceed in three steps. (1) Because the contact points move along the antidiagonal (lower bands) or the diagonal (upper bands), we will first diagonalize $H(\mathbf{k})$ along these two directions for arbitrary $t$. A great help in this process comes from the mirrors symmetries. In this way we will obtain useful eigenvector bases (see Appendix A). (2) Next, for a specific $\mathbf{k}_{0}$ point along these directions, we will rewrite the Bloch Hamiltonian in the eigenvector basis at $\mathbf{k}_{0}$ and then expand in $\mathbf{q}=\mathbf{k}-\mathbf{k}_{0}$ up to first or second order. In the following, $\mathbf{k}_{0}$ will be chosen either as a TRIM ($\Gamma$ or M points) or one of the Dirac points ${\bf k}_{D}$. (3) The third step consists in projecting from three to two bands to obtain an effective two-band Hamiltonian for Dirac points.

In the following, we study the lower (Sec. III) and upper Dirac points (Sec. IV) in turn.

III Lower Dirac points

The lower Dirac points move along the antidiagonal direction $k_{y}=-k_{x}$ from M to $\Gamma$ point as $t$ is varied from 0 to 1, see Sect. II.2 and Fig. 3. Specifically, their positions and energy are given by

[TABLE]

where $\xi=\pm 1$ is a valley index and the angle $\phi_{a}\geq 0$ is defined by

[TABLE]

III.1 Basis along the antidiagonal

The Bloch Hamiltonian commutes with the mirror symmetry operator $S_{a}({\bf k})$ everywhere along the antidiagonal line ($k_{y}=-k_{x}$) in the Brillouin zone. As detailed in Appendix A, this facilitates the diagonalization of the Hamiltonian using common eigenstates with the symmetry operator, which we denote as $(|u_{+}\rangle,|u_{-}\rangle,|u_{3}\rangle)$. We then form the unitary matrix $L(k_{x})=(|u_{+}\rangle\,\,|u_{-}\rangle\,\,|u_{3}\rangle)$ that allows us to write the Bloch Hamiltonian in its eigenbasis as

[TABLE]

with $E_{\pm}(k_{x})=t\pm\Delta_{a}$, $E_{3}=E_{a}=-2t$, where $\Delta_{a}=\sqrt{t^{2}+8\cos^{2}\frac{k_{x}}{2}}$. The largest eigenvalue is $E_{+}$ and the second is $E_{-}\geq E_{3}$ if $t\geq\cos(k_{x}/2)$ and $E_{3}\geq E_{-}$ otherwise. Along the antidiagonal line, they describe the three energy bands. At the Dirac point $k_{x}=2\phi_{a}$, the lower two bands become degenerate $E_{3}=E_{-}=-2t$ and separated from the third band $E_{+}=4t$. We emphasize that, while there is always a gauge freedom in the choice of the set of eigenvectors, we use here an eigenbasis which connects smoothly the two limiting cases from $t=0$ to $t=1$.

By writing the original Hamiltonian (4) in the new basis $\tilde{H}({\bf k})=L^{\dagger}(k_{x})H({\bf k})L(k_{x})$, we now want to successively describe the vicinity of the M point, the $\Gamma$ point and the Dirac points (and not limited to $k_{y}=-k_{x}$), by appropriate expansions.

III.2 Close to the Lieb limit: emergence of Dirac points at M point

This is the limit $t\to 0$ when the lower Dirac points pair emerge from the triple crossing at $M=(\pi,\pi)$ point of the Lieb lattice along the antidiagonal line. Upon a momentum expansion at M point for small $\mathbf{q}$, we get $\tilde{H}_{M}({\bf q})=L(\pi)^{\dagger}H(M+{\bf q})L(\pi)$ in the form

[TABLE]

with $q_{\parallel}\equiv\frac{q_{x}+q_{y}}{\sqrt{2}}$ (“diagonal” direction) and $q_{\perp}\equiv\frac{q_{y}-q_{x}}{\sqrt{2}}$ (“antidiagonal” direction). We keep the second order in $q_{\parallel}$, $q_{\perp}$ and $t$ in the relevant blocks (indeed, we will see that the Dirac points correspond to $|q_{\perp}|\sim t$). The vertical/horizontal lines in the matrix merely serve as an eye guide separating the subspaces $(|u_{+}\rangle)$ and $(|u_{-}\rangle,|u_{3}\rangle)$ which the matrix operates.

We then eliminate the highest band using second order perturbation theory. Following the Löwdin method (see Appendix B) Lowdin51 with typical energy $E_{0}=E_{a}=-2t$, we obtain an effective $2\times 2$ Hamiltonian acting in $(|u_{-}\rangle,|u_{3}\rangle)$ subspace as

[TABLE]

at first order in $t$. Apart from the identity term $\sigma_{0}$, it has the form of the universal Hamiltonian $\mathcal{H}_{+-}$, see Eq. (6), with $\sigma_{\mu}=\sigma_{z}$ and $\sigma_{\nu}=\sigma_{y}$. It describes the merging of two Dirac points located at $(q_{\parallel},q_{\perp})=(0,\xi 2\sqrt{2}t)$ and of opposite winding numbers in the $\sigma_{y}-\sigma_{z}$ plane. The winding vectors are $\vec{w}_{\xi}=-\xi\vec{u}_{x}$.

At second order in $t$, the effective two-band Hamiltonian (18) gets an additional contribution, which, as we show below, is responsible of a rotation of the winding vectors:

[TABLE]

Upon the substitution $q_{\perp}\to\xi 2\sqrt{2}t+q_{\perp}$, the linearized Hamiltonian near the Dirac points becomes:

[TABLE]

The Dirac cones are critically tilted because of the fine-tuning of the $q_{\perp}\sigma_{0}$ term with respect to the $q_{\perp}\sigma_{z}$ term. Also, there is some mixing between $\sigma_{y}$ and $\sigma_{x}$ indicating that the winding vectors rotate. Indeed $\vec{w}_{\xi}=-\xi\vec{u}_{x}+t\vec{u}_{y}$, which evolve from $-\xi\vec{u}_{x}$ (i.e. anti-parallel) at small $t$ to $\vec{u}_{y}$ (i.e. parallel) when increasing $t$.

We remark that the effective Hamiltonian near the M point has a peculiar structure. It results from an expansion a priori valid for small $\mathbf{q}$. But near the merging when $t\rightarrow 0$, one cannot restrict to the two lowest band, since three bands become degenerate in the Lieb limit. Moreover, in this limit, the curvature of the spectrum goes to infinity, which is natural since the spectrum at the merging is indeed conical. The purpose here is to show that the Dirac points show opposite winding numbers and emerge following the $(+-)$ scenario. It should not be applied strictly at the Lieb limit $t=0$ (see the corresponding discussion in Appendix C).

III.3 Close to the kagomé limit: merging of Dirac points at $\Gamma$ point

Close to the kagomé limit ($t=1$), we consider an expansion of the Hamiltonian at the $\Gamma=(0,0)$ point to give $\tilde{H}_{\Gamma}({\bf q})=L(0)^{\dagger}H(\Gamma+{\bf q})L(0)$ as

[TABLE]

with $t=1-\epsilon$ and $\epsilon\to 0$. Similarly to the previous section we keep up to the second order in the variables $q_{\parallel}$, $q_{\perp}$ and first order in $\epsilon$ (as we will see that the Dirac points correspond to $q_{\perp}^{2}\sim\epsilon$) in the relevant matrix elements. Because we will treat the highest band as a perturbation to the lower sub-space, we only need to consider the off-diagonal blocks at first order and the highest band energy at zeroth order. From Löwdin method (see Appendix B), we arrive at the effective two-band Hamiltonian

[TABLE]

at second order in $q$ and first order in $\epsilon$. Besides the identity term, it has the form of the universal Hamiltonian $\mathcal{H}_{++}$, see Eq. (7) with $\sigma_{\mu}=\sigma_{z}$ and $\sigma_{\nu}=\sigma_{x}$, following the $(++)$ merging scenario. It describes the merging of two Dirac points located at $(q_{\parallel},q_{\perp})=(0,-\xi 4\sqrt{\epsilon})$ if $\epsilon>0$ [beware that the $+{\bf k}_{a}$ Dirac point is at $q_{\perp}=-4\sqrt{\epsilon}$] and of identical winding numbers in the $\sigma_{x}-\sigma_{z}$ plane. Note that for $\epsilon<0$, the merging is along the $q_{\parallel}$ direction as the Dirac points are located at $(q_{\parallel},q_{\perp})=(\xi 4\sqrt{-\epsilon/3},0)$. This agrees with the model of Eq. (7) showing the emergence for $t\leq 1$ as well as $t\geq 1$.

When $\epsilon=0$, this Hamiltonian describes a quadratic band crossing in the kagomé limit. A similar Hamiltonian was obtained by Ref. Xiao2003 in the triangular version of the kagomé lattice. However their Hamiltonian – see Eq. (2.6) in Ref. Xiao2003 with $d_{3}=0$ – breaks time-reversal symmetry which is not possible at the $\Gamma$ point. This is due to taking the Pauli matrices $\sigma_{x}$, $\sigma_{y}$ instead of $\sigma_{z}$, $\sigma_{x}$ in the Hamiltonian $\mathcal{H}_{++}$, see Eq. (7).

Finding a rotation of the winding vector requires performing a computation at third order in $q\sim\sqrt{\epsilon}$ in Eq. (24). The Löwdin method then gives the additional terms as:

[TABLE]

Upon the substitution $q_{\perp}\to-\xi 4\sqrt{\epsilon}+q_{\perp}$, the linearized Hamiltonian near the Dirac points becomes:

[TABLE]

The winding vector is $\vec{w}_{\xi}=\vec{u}_{y}-\xi\frac{2}{3}\sqrt{2\epsilon}\vec{u}_{x}$, when $\epsilon$ is small.

III.4 Motion of the lower Dirac points

We now find the paradoxical situation where a $(+-)$ pair of Dirac points near M (when $t\simeq 0$) continuously evolves into a $(++)$ pair near $\Gamma$ (when $t\simeq 1$). Naively one expects that a topological charge conservation, which holds locally in the momentum space, should also hold globally in the full BZ. Here, the global charge seems to change from [math] (for $t=0$) to $+2$ (for $t=1$).

To solve the apparent paradox, we now follow the evolution of the Dirac points and determine the evolution of the associated winding vector. We make an expansion following the Dirac point at ${\bf k}_{a}=(2\phi_{a},-2\phi_{a})$, where the angle $\phi_{a}$ is defined in Eq. (12), to get $\tilde{H}_{a}({\bf q})=L(2\phi_{a})^{\dagger}H(\mathbf{k}_{a}+\mathbf{q})L(2\phi_{a})$ as

[TABLE]

where $E_{a}=-2\cos\phi_{a}=-2t$. For the other Dirac point at $-{\bf k}_{a}$, one should replace $\phi_{a}$ by $-\phi_{a}$. Because we are now interested in a Dirac Hamiltonian, which is linear in momentum, and in contrast to the two cases previously studied, we only need to keep the momentum variable up to the linear order. Corrections from the third band are necessarily of the second order and can be neglected. The effective Hamiltonian can therefore be read off from the $2\times 2$ block spanned by $(|u_{-}\rangle,|u_{3}\rangle)$ (see Appendix A) as

[TABLE]

with the (positive) velocities

[TABLE]

In the $t\to 0$ limit, $v_{\parallel}\simeq\sqrt{2/3}$ and $v_{\perp}=v_{\perp}^{0}\simeq\sqrt{2}/3$; whereas when $t=1-\epsilon\to 1$, $v_{\parallel}\simeq 2\sqrt{\epsilon/3}$ and $v_{\perp}=v_{\perp}^{0}\simeq 2\sqrt{\epsilon}/3$. These velocities agree with that found in the expansion near the M point (when $t\to 1$) but there is a numerical discrepancy with that found near the $\Gamma$ point (when $t\to 0$).

The Pauli matrix $\sigma_{\phi_{a}}\equiv\cos\phi_{a}\sigma_{x}+\sin\phi_{a}\sigma_{y}$, together with $\sigma_{z}$, defines the pseudospin winding plane. This plane evolves in the pseudospin space from $(\sigma_{y},\sigma_{z})$ to $(\sigma_{x},\sigma_{z})$ with $\phi_{a}$ going from $\pi/2$ ($t=0$) to 0 ($t=1$), see Fig. 6. For the two Dirac points at $\xi{\bf k}_{a}=\xi(2\phi_{a},-2\phi_{a})$, following Eq. (10), the associated winding vectors are:

[TABLE]

In the $t\to 0$ limit, $\vec{w}_{\xi}\simeq-\xi\vec{u}_{x}+t\,\vec{u}_{y}$; whereas when $t=1-\epsilon\to 1$, $\vec{w}_{\xi}\simeq-\xi\sqrt{2\epsilon}\vec{u}_{x}+\vec{u}_{y}$. Apart from a factor of $\frac{2}{3}$ when $t$ is close to $1$, this agrees with the results of the expansions near the M and $\Gamma$ points.

As an aside, we note that the $\sigma_{0}$ term modifies the dispersion of the Dirac cone quite drastically. As $v_{\perp}=v_{\perp}^{0}$ for all $t$, the cones are critically tilted (also known as type III Dirac-Weyl cones, see e.g. Ref. Milicevic19 and references therein): the velocities of excitations along the $\perp$ direction are $v_{min}=v_{\perp}-v_{\perp}^{0}=0$ and $v_{max}=v_{\perp}+v_{\perp}^{0}=2v_{\perp}$. The pair of Dirac points is connected by a flat energy line (see red dots in Figs 4 and 5).

In summary, for the crossings in the lower two bands, Eqs. (18) and (25) describes the emergence/merging of two Dirac cones with different winding scenarios. In the single Dirac cone Hamiltonian (32), we follow the evolving winding vector and show how the apparent winding number gets reversed as $t$ is varied.

IV Upper Dirac points

For the upper Dirac point pair, we perform an analysis similar to the one contained in the last section. The Dirac points move along the diagonal line $k_{y}=k_{x}$, with their positions and energy (Fig. 3) given by

[TABLE]

and specified by the angle $\phi_{d}\geq 0$ such that

[TABLE]

The upper Dirac points that emerge from the M point ($t=0$) move to a fixed position $\xi{\bf K}=\xi(2\pi/3,2\pi/3)$ in the kagomé limit ($t=1$) and do not merge (in contrast to the lower Dirac points pair).

IV.1 Basis along the diagonal

The Bloch Hamiltonian commutes with the mirror symmetry operator $S_{d}({\bf k})$ everywhere along the diagonal line ($k_{y}=k_{x}$) in the Brillouin zone. As detailed in Appendix A, a common basis of eigenvectors of the Hamiltonian and the symmetry operator is found, which we denote as $(|v_{1}\rangle,|v_{+}\rangle,|v_{-}\rangle)$. We emphasize that this basis furnishes a smooth interpolation between the two limits $t=0$ and $t=1$. The unitary matrix is defined by $U(k_{x})\equiv(|v_{1}\rangle\,\,|v_{+}\rangle\,\,|v_{-}\rangle)$ and written in this basis, the Hamiltonian along the diagonal line is given by

[TABLE]

with $E_{1}=-2t\cos k_{x}$ and $E_{\pm}=t\cos k_{x}\pm\Delta_{d}$, where $\Delta_{d}=\sqrt{t^{2}\cos^{2}k_{x}+8\cos^{2}\frac{k_{x}}{2}}$. We see that the two upper bands become degenerate at energy $E_{1}=E_{+}=2\cos\phi_{d}=-2t\cos(2\phi_{d})$ separated from the lowest band $E_{-}=-4\cos\phi_{d}$ at the Dirac point.

IV.2 Close to the Lieb limit: emergence of Dirac points at M point

At $M=(\pi,\pi)$, the unitary matrix $U(\pi)=(|v_{1}\rangle\,\,|v_{2}\rangle\,\,|v_{3}\rangle)$ is equal to $L(\pi)=(|u_{1}\rangle\,\,|u_{2}\rangle\,\,|u_{3}\rangle)$ as $|v_{j}\rangle=|u_{j}\rangle$. Instead of Eq. (17), we focus on the upper two bands to obtain $\tilde{H}_{M}({\bf q})=U(\pi)^{\dagger}H(M+{\bf q})U(\pi)$ as

[TABLE]

at second order in the variables $q_{\parallel}$, $q_{\perp}$ and $t$ (as we will see that the Dirac points correspond to $q_{\parallel}^{2}\sim t^{2}$) in the relevant blocks. Using Löwdin method with typical energy $E_{0}=E_{d}\simeq 2t$, when $t\to 0$, see Appendix B, we arrive at the effective two-band Hamiltonian (first order in $t$)

[TABLE]

in the subspace of $(|v_{1}\rangle,|v_{+}\rangle)$. Apart from the identity $\sigma_{0}$ term, it takes the form of the universal Hamiltonian $\mathcal{H}_{+-}$, see Eq. (6) with $\sigma_{\mu}=\sigma_{z}$ and $\sigma_{\nu}=\sigma_{y}$, describing the emergence of two Dirac points at $(q_{\parallel},q_{\perp})=(-\xi 2\sqrt{2}t,0)$ [beware that the $+{\bf k}_{d}$ Dirac point is at $q_{\parallel}=-2\sqrt{2}t$] of opposite winding vectors $\vec{w}_{\xi}=\xi\vec{u}_{x}$ for $t\simeq 0$. The Hamiltonians (18) and (42) are similar upon exchanging the diagonal $q_{\parallel}$ and the antidiagonal $q_{\perp}$. They both describe a pair of critically tilted Dirac cones. The critical tilt is true for the lower cones but actually not for the upper ones, as seen on Fig. 4(c), and because of neglected higher-order terms [cf. $v_{\parallel}^{0}\neq v_{\parallel}$ in Eq. (46)].

At second order in $t$, the effective two-band Hamiltonian (42) gets an additional contribution:

[TABLE]

This term is similar to Eq. (19) except that $q_{\parallel}q_{\perp}$ is replaced by $q_{\perp}^{2}$. Whereas Eq. (19) leads to a rotation of the winding vector, it is not the case of Eq. (43). Indeed, upon the substitution $q_{\parallel}\to-\xi 2\sqrt{2}t+q_{\parallel}$, the linearized Hamiltonian near the Dirac points reads

[TABLE]

In contrast to the lower Dirac points pair, here there is no rotation of the winding vector, which is $\vec{w}_{\xi}=\xi\vec{u}_{x}$ independently of $t$. Since the Dirac cones do not approach each other again in the kagomé limit $t\rightarrow 1$, this concludes the analysis for double Dirac cones in the upper two bands.

IV.3 Motion of the upper Dirac points

To follow the evolution of a single Dirac cone, we first expand the Hamiltonian $\tilde{H}_{d}({\bf q})=U(\xi 2\phi_{d})^{\dagger}H({\bf k}_{d}+{\bf q})U(\xi 2\phi_{d})$ at the position $\xi{\bf k}_{d}$ keeping linear order in ${\bf q}$. The effective Hamiltonian for the Dirac point is then obtained directly from the block Hamiltonian in the $(|v_{1}\rangle,|v_{+}\rangle)$ subspace (see Appendix A) as

[TABLE]

where $E_{d}=2\cos\phi_{d}=\frac{\sqrt{1+8t^{2}}-1}{2t}$ and we defined the (positive) velocities

[TABLE]

In the $t\to 0$ limit, $v_{\parallel}^{0}\simeq v_{\parallel}\simeq\sqrt{2}/3$ (critical tilt) and $v_{\perp}\simeq\sqrt{2/3}$; whereas when $t\to 1$, $v_{\parallel}^{0}\simeq 0$ (no tilt), $v_{\parallel}\simeq\sqrt{3/2}$ and $v_{\perp}\simeq 1/\sqrt{2}$. Except when $t=0$, $v_{\parallel}^{0}\neq v_{\parallel}$, which means that the tilt is no critical. An effect not captured by Eq. (42).

When $t$ varies from [math] to $1$, the upper pair emerges from the M point and reaches two fixed points $\xi{\bf K}=\xi({2\pi\over 3},{2\pi\over 3})$ [see Fig. (3a)]. Following Eq. (10), the winding vectors associated to the two Dirac points $\xi{\bf k}_{d}$ are $\vec{w}_{\xi}=\xi\vec{u}_{x}$. In this case, there is no rotation of the winding vectors, which remain anti-parallel.

V Discussion and conclusion

We have studied in detail the evolution of Dirac points in a tight-binding model interpolating between two well-known three-bands models, the Lieb and the kagomé lattices. In both limits, the energy spectrum exhibits a flat band, symmetrically positioned in the middle of a conical spectrum for the Lieb lattice, at a quadratic touching with a dispersive band for the kagomé lattice. One of the interests of this study is to interpolate between two different flat bands: in the Lieb case, the flat band is topological (as a result of the lattice being bipartite), whereas in the kagomé case it results from destructive interferences. In the first case, the flat band resists the introduction of a perpendicular magnetic field, whereas it is destroyed in the second case (see the corresponding Hofstadter butterflies in Ref. Aoki1996 ).

Starting from the kagomé lattice, the quadratic contact point is split into two Dirac points characterized by the same winding number, following a well-characterized universal scenario. Upon further variation of the interpolating parameter, these two Dirac points merge, following another well-characterized merging scenario for a pair of Dirac points with opposite winding numbers. To solve this apparent contradiction and to characterize the topological properties of these contact points, we have reduced the full $3\times 3$ Hamiltonian to effective $2\times 2$ Hamiltonians using two different approaches (see below). During the evolution between the two limits, the effective Hamiltonian (or more precisely its effective pseudo-magnetic field $\vec{h}$) rotates in pseudo-spin space, leading to the notion of winding vector. The winding vectors of a pair evolve from a parallel to an anti-parallel alignment.

This work poses questions about a systematic way to obtain effective two-band Hamiltonians describing contact points in multiband models. Here, we have essentially used two different strategies: (i) an expansion of the three-band Bloch Hamiltonian at second (or third) order in the wavevector in the vicinity of a TRIM, followed by a projection (using second-order perturbation theory) in order to capture the pair of contact points located away from the TRIM; (ii) an expansion of the three-band Bloch Hamiltonian at first order in the wavevector directly in the vicinity of the contact point at ${\bf k}_{D}$ followed by a naive projection (without the need of perturbation theory).

Such a reduction from a multiband to an effective $2\times 2$ Hamiltonian for a pair of Dirac points is intrinsically local in ${\bf k}$ space because it is local in energy space: it can not be done in the whole Brillouin zone. It is important to realize that such a procedure carries some arbitrariness as it depends on a choice of eigenbasis. This is especially crucial for procedure (ii), which depends on an eigenbasis at the Dirac point. Indeed any linear combination of degenerate eigenvectors produces a valid alternative eigenbasis. The consequence is that the direction of the winding vector for a single Dirac point has no absolute meaning and can be changed almost at will by choosing a different basis of degenerate eigenvectors.

However, the relative direction between the winding vectors of the two Dirac points within a pair has an absolute meaning. This appears clearly when computing the winding vectors in the same basis for the two Dirac points, as done in procedure (i). In other words, the effective $2\times 2$ Hamiltonian should describe both Dirac points at once. This is typically a “universal” Hamiltonian – similar to $\mathcal{H}_{+-}$ or $\mathcal{H}_{++}$ – built from a TRIM (either $\Gamma$ or M point in the present context). As we have shown, higher-order terms in momentum provide corrections to these universal Hamiltonians that can describe the rotation of the winding vectors.

In the future, it would be interesting to fully characterize a multiband problem by introducing angles that describe band couplings between more than two bands and that are valid in the complete Brillouin zone and not only locally in reciprocal space. Here, we have seen how the Lieb point splits in four Dirac points as the hopping $t$ becomes finite. It would be interesting to systematically study the different splitting scenarios of degeneracy points between three bands.

Acknowledgements

We thank Gilles Abramovici for useful discussions on the Lieb-kagomé model. L.-K. L. is supported by the Thousand Young Talents Program of China. We thank the Institute for Advanced Study of the Tsinghua University in Beijing for hosting us while part of this work was being performed.

Appendix A Basis along the antidiagonal and the diagonal

This appendix provides details of the derivation of the two distinct eigenbases $(|u_{\pm}(k_{x})\rangle,|u_{3}(k_{x})\rangle)$ and $(|v_{1}(k_{x})\rangle,|v_{\pm}(k_{x})\rangle)$ that diagonalize the Bloch Hamiltonian $H({\bf k})$, for any value of the parameter $t$, along respectively the antidiagonal $k_{x}=-k_{y}$ and diagonal $k_{x}=k_{y}$. The key ingredients are the mirror symmetries $S_{a}$ and $S_{d}$ along the diagonal and antidiagonal respectively.

To start with, it is important to remind that the specific form of the Bloch Hamiltonian $H({\bf k})$ written in Eq. (4) is a direct consequence of the implicit choice of the Bloch basis $|{\bf k},\alpha\rangle=\sum_{{\bf R}}e^{i{\bf k}\cdot{\bf R}}|{\bf R},\alpha\rangle$ ($\alpha=A,B,C$) that depends only on the Bravais lattice vectors ${\bf R}$ (choice known as basis I Bena2009 ). The main advantage of this Bloch basis is that the Bloch Hamiltonian matrix $H({\bf k})$ representing the Hamiltonian is periodic $H({\bf k}+{\bf G})=H({\bf k})$ for any reciprocal lattice vectors ${\bf G}$. However an important drawback of this representation is that the matrix $H({\bf k})$ has lost the information of the full position of the different atoms in the unit cell. In fact, the latter information is now encoded in the matrices representing the different space symmetries of the lattice, highlighting their key role.

A.1 Diagonalization of $H({\bf k})$ along the antidiagonal

The Bloch matrix representation of the mirror symmetry $S_{a}$ along the real-space antidiagonal is given by:

[TABLE]

This symmetry translates into the relation

[TABLE]

for the Bloch Hamiltonian in reciprocal space. In particular, $S_{a}({\bf k})$ commutes with $H({\bf k})$ on the antidiagonal $k_{x}=-k_{y}$. We therefore diagonalize $S_{a}(k_{x},-k_{x})$ to get the eigenvectors

[TABLE]

with respective eigenvalues $(+1,+1,-1)$. Written in this basis, the Hamiltonian transforms into a block diagonal matrix:

[TABLE]

in the $(|u_{1}\rangle,|u_{2}\rangle,|u_{3}\rangle)$ basis. The $2\times 2$ submatrix in the subspace $|u_{1}\rangle,|u_{2}\rangle$ can be rewritten

[TABLE]

by introducing the parameters $\Delta_{a}(k_{x})=\sqrt{t^{2}+8\cos^{2}\frac{k_{x}}{2}}$ and the angle $\theta_{a}(k_{x})$ with $\cos\theta_{a}=t/\Delta_{a}$ and $\sin\theta_{a}=2\sqrt{2}\cos{k_{x}\over 2}/\Delta_{a}$. This has the standard form of a two-level problem with its parametric solution.

Therefore the eigenbasis and their eigenvalues which diagonalizes the Hamiltonian $\tilde{H}_{a}(k_{x})$ everywhere along the anti-diagonal are

[TABLE]

and $|u_{3}\rangle$, with their corresponding eigenvalues $E_{\pm}(k_{x})=t\pm\Delta_{a}$, $E_{3}=-2t$. Finally, the chosen eigenbasis for $H(k_{x},-k_{x})$ reads:

[TABLE]

A.2 Diagonalization of $H({\bf k})$ along the diagonal

We now proceed very similarly along the diagonal $k_{x}=k_{y}$. The Bloch matrix representation of the mirror symmetry $S_{d}$ along the real-space diagonal is given by:

[TABLE]

This symmetry translates into the relation

[TABLE]

for the Bloch Hamiltonian in reciprocal space. In particular, $S_{d}({\bf k})$ commutes with $H({\bf k})$ on the diagonal $k_{x}=k_{y}$. We therefore diagonalize $S_{d}(k_{x},k_{x})$ to get the eigenvectors

[TABLE]

with respective eigenvalues $(-1,+1,+1)$, where $s$ is defined as $s=\text{sign }k_{x}$ if $0<|k_{x}|\leq\pi$ and $s=0$ if $k_{x}=0$. The reason for this phase choice is that we want the kets $|v_{j}\rangle$ to be real at $\Gamma$ and M, to be periodic with the BZ and that $|v_{j}(-k_{x})\rangle=(|v_{j}(k_{x})\rangle)^{*}$. Written in this basis, the Hamiltonian thus transforms into a block diagonal matrix:

[TABLE]

in the $(|v_{1}\rangle,|v_{2}\rangle,|v_{3}\rangle)$ basis. Consider the $2\times 2$ effective submatrix $\mathcal{H}_{v_{2},v_{3}}(k_{x})$ in the subspace $(|v_{2}\rangle,|v_{3}\rangle)$; we obtain an effective Hamiltonian of the same structure as in Eq. (59) with the substitution $\Delta_{a}(k_{x})\rightarrow\Delta_{d}(k_{x})=\sqrt{t^{2}\cos^{2}k_{x}+8\cos^{2}\frac{k_{x}}{2}}$ and the angle $\theta_{a}(k_{x})\rightarrow\theta_{d}(k_{x})$ with $\cos\theta_{d}=-t\cos k_{x}/\Delta_{d}\geq 0$ and $\sin\theta_{d}=2\sqrt{2}\cos{k_{x}\over 2}/\Delta_{d}$, and the coefficient for the identity term $t\rightarrow t\cos k_{x}$. Thus, a basis that diagonalizes the Hamiltonian $\tilde{H}_{d}(k_{x})$ everywhere along the diagonal is

[TABLE]

and $|v_{1}\rangle$, with the eigenvalues $E_{\pm}=t\cos k_{x}\pm\Delta_{d}$ and $E_{1}=-2t\cos k_{x}$. Finally, the chosen eigenbasis for $H(k_{x},k_{x})$ reads:

[TABLE]

Appendix B Löwdin partitioning method for the effective $2\times 2$ Hamiltonian

The effective Hamiltonian in a sub-space can be obtained using the method of Löwdin Lowdin51 , also called partitioning. Let the Hilbert space be described as the direct sum of two sub-spaces labelled $\alpha$ and $\beta$. The Hamiltonian then takes the following block form:

[TABLE]

It is then easy to show that the eigen-equation $H|\psi\rangle=E|\psi\rangle$ reads

[TABLE]

in the $\alpha$ sub-space, where $|\psi\rangle=(|\psi_{\alpha}\rangle,|\psi_{\beta}\rangle)^{T}$. The operator in the left-hand side of the above equation acts as an effective Hamiltonian in the $\alpha$ sub-space, except that it depends on the energy $E$, i.e. (79) is a self-consistent equation in the spirit of the well-known Brillouin-Wigner perturbation theory. Replacing $E$ in the left-hand side by a typical relevant energy $E=E_{0}$ gives an effective Hamiltonian

[TABLE]

so that the approximate eigen-equation replacing (79) reads $\mathcal{H}_{\alpha}|\psi_{\alpha}\rangle\simeq E|\psi_{\alpha}\rangle$. Compared to a naive projection $\mathcal{H}_{\alpha}\simeq H_{\alpha\alpha}$, this effective Hamiltonian includes corrections $H_{\alpha\beta}\left(E_{0}-H_{\beta\beta}\right)^{-1}H_{\beta\alpha}$ akin to second-order perturbation theory.

Specifically, for the $3\times 3$ Hamiltonian

[TABLE]

we get the effective $2\times 2$ Hamiltonian after eliminating the third band which reads

[TABLE]

with

[TABLE]

With this method there is some arbitrariness in the choice of $E_{0}$. In the present paper, we are interested in describing Dirac points as seen from the nearest TRIM (e.g., $\Gamma$ or M point in the BZ): this means that $E_{0}$ will be chosen to be the energy of the Dirac points – called $E_{a}$ and $E_{d}$ in the main text – in order to properly describe the vicinity of the Dirac points and not the vicinity of the TRIM.

Appendix C Generalized Lieb-kagomé model

The effective $2\times 2$ Hamiltonian (18) describes the evolution of the two lower Dirac points along the antidiagonal, close to the Lieb limit (red points on Fig. 4). Similarly the Hamiltonian (42) describes the evolution of the two upper Dirac points along the diagonal (blue points on Fig  4). However, when $t\rightarrow 0$, a two-band description is not appropriate, since the four points merge together. The Hamiltonians (18, 42) provide a correct description of the merging of each pair close to the M point, but not too close. They do not correctly describe the $t=0$ limit. This is obvious when considering the divergence of the $q_{\perp}^{2}/8t$ and $q_{\parallel}^{2}/8t$ terms respectively in Eqs. (18) and (42). In order to show that the merging of the two lower Dirac points is actually properly described by the ${\cal H}_{+-}$ Hamiltonian of Eq. (18), we reconsider the same problem by gapping the upper band from the two lower bands. This is done by adding a dimerization along vertical and horizontal bonds (which breaks the mirror symmetry with respect to the diagonal and the inversion center), as shown in Fig. 7(Top). We introduce a parameter $t^{\prime}\neq 1$ on the $A-B$ bond while the coupling is still $1$ on the $B-A$ bond. The Bloch Hamiltonian (basis I) is

[TABLE]

We call it the generalized Lieb-kagomé model with three hopping amplitudes. This generalization is inspired by the dimerization of the Lieb lattice used in Poli2017 to gap the three bands of the Lieb model. It is also close to, but different from, the anisotropic kagomé model discussed in Asano2011 .

Using the same procedure as in section III.2, we construct the $2\times 2$ effective Hamiltonian:

[TABLE]

with

[TABLE]

and the velocity

[TABLE]

The parameters have the following limits:

[TABLE]

The limit $t\rightarrow 0$ is now well behaved when $t^{\prime}\neq 1$, since the parameter $\delta$ stays finite up to the merging transition which is reached when $\delta=3t$, that is $t=(1-t^{\prime})/2$. The full evolution is shown on Fig. 7

Appendix D Real Bloch Hamiltonian and orientability

In this appendix, we study Bloch Hamiltonians that are real and discuss the notion of orientability as defined in section VI.B of Ref. Ahn2019 . Real Bloch Hamiltonians are possible when the system has both space inversion and time-reversal symmetries. For pedagogical purposes, we start with the two-band staggered Mielke model Montambaux2018 before turning to the three-band Lieb-kagomé model.

We first need to introduce two different Bloch Hamiltonians that can be defined for any tight-binding model possessing more than a single site per unit cell. On the one hand, the Bloch Hamiltonian in the so-called basis I Bena2009 is defined as

[TABLE]

where $\mathbf{R}$ is the Bravais lattice position operator. It is periodic with any Bravais lattice vector ${\bf G}$:

[TABLE]

This basis I is also the convention we choose in the core of the article when studying the Lieb-kagomé model. On the other hand, the Bloch Hamiltonian in the so-called basis II Bena2009 is defined as

[TABLE]

where $\mathbf{r}=\mathbf{R}+\boldsymbol{\delta}$ is the complete position operator involving the Bravais lattice position $\mathbf{R}$ and the intra-cell position operator $\boldsymbol{\delta}$. For certain Bravais lattice vectors ${\bf G}$, it is not periodic but obeys:

[TABLE]

D.1 Mielke model with a real Bloch Hamiltonian

Here we discuss the case of a two-band real Bloch Hamiltonian, having in mind the staggered Mielke model Montambaux2018 . This model has inversion and time-reversal symmetries and the inversion center can be chosen on-site. It is defined on the checkerboard lattice. The Bloch Hamiltonian is real and non-periodic in basis II, whereas it is complex and periodic in basis I (it involves the three Pauli matrices). In basis II, under translation by a Bravais lattice vector ${\bf G}$, it behaves as Eq. (101). When ${\bf G}$ is a basis vector of the reciprocal lattice then $e^{i{\bf G}\cdot\boldsymbol{\delta}}=\sigma_{z}$. The fact that $\det\sigma_{z}=-1$ means that the Bloch Hamiltonian is non-orientable (see section VI.B in Ref. Ahn2019 ). In other words, it is not possible to find a representation in which the Bloch Hamiltonian would be both real and periodic with the first BZ.

The real Bloch Hamiltonian can be written as

[TABLE]

which defines an angle $\theta({\bf k})$. It produces a checkerboard pattern in reciprocal space with a unit cell twice larger than the first BZ (see the yellow and blue zones in Fig. 8). This is also the case of the pattern of topological charges (winding numbers) as show by plusses and minuses in the same figure.

Orientability of the Bloch Hamiltonian is related to being able to find a pattern of topological charges which has the periodicity of the first BZ. There is an obstruction to that in the Mielke model (see Fig. 8). Fundamentally, it is this obstruction that allows the evolution of topological charges (i.e. the “++ $\to$ +-” phenomenon).

D.2 Lieb-kagomé model with a real Bloch Hamiltonian

As an alternative to the core of the paper, we present here the evolution of the topological charges of band contact points by writing the Bloch Hamiltonian of the Lieb-kagomé model in basis II instead of basis I Bena2009 , as was done in Eq. (4). In basis II, the Bloch Hamiltonian reads:

[TABLE]

This Bloch Hamiltonian turns out to be real as a consequence of time-reversal symmetry $H_{II}(-\mathbf{k})^{*}=H_{II}(\mathbf{k})$ and inversion symmetry $H_{II}(-\mathbf{k})=H_{II}(\mathbf{k})$ with an on-site inversion center. In basis I, space inversion acts in a different way on the Bloch Hamiltonian and does not force it to be real.

In contrast to $H_{I}(\mathbf{k})$ [see Eq. (4)] that satisfies $H_{I}(k_{x}+2\pi,k_{y})=H_{I}(k_{x},k_{y})=H_{I}(k_{x},k_{y}+2\pi)$, the Bloch Hamiltonian $H_{II}(\mathbf{k})$ does not have the periodicity of the reciprocal lattice but a doubled-periodicity in both $k_{x}$ and $k_{y}$ such that $H_{II}(k_{x}+4\pi,k_{y})=H_{II}(k_{x},k_{y})=H_{II}(k_{x},k_{y}+4\pi)$ and $H_{II}(k_{x}+2\pi,k_{y}+2\pi)=-H_{II}(k_{x},k_{y})$. See Figure 9.

Because $H_{II}(\mathbf{k})$ is real, it means that any local two-band Hamiltonian $\mathcal{H}_{II}(\mathbf{q})$ obtained after projection in the vicinity of a band contact point is also real. Therefore $\mathcal{H}_{II}(\mathbf{q})$ can be decomposed onto $\sigma_{x}$ and $\sigma_{z}$ only, so that the corresponding contact point has a winding vector $\vec{w}=w\vec{u}_{y}$ and the notion of winding number $w$ is sufficient. This already means that, in basis II, there is no such thing as a continuous evolution of the winding vector.

We only consider the contact points between the lower two bands. We find that, close to the kagomé limit and by expanding in the vicinity of the $\Gamma$ point (similar to section III.3), the pair of Dirac point has identical winding numbers (either $++$ or $--$ depending on the basis choice at $\Gamma$). In the Lieb limit and by expanding in the vicinity of the M point (similar to section III.2), we find that the pair of Dirac points has opposite winding numbers (either $+-$ of $-+$ depending on the basis chosen at M). How is this possible? Actually, the Bloch Hamiltonian does not have the periodicity of the reciprocal lattice but a doubled periodicity. It means that the topological charges of the contact points need not have the periodicity of the reciprocal lattice. In other words, the charges attributed to the band contact points within the first BZ are not necessarily the same in “another BZ”. Possible patterns respecting the above rules (same charges near $\Gamma$ opposite charges near M) are shown in red in Fig. 10. The unit cell for the winding numbers is two times larger than the first BZ.

It is not enough to have a Bloch Hamiltonian with enlarged periodicity compared to the BZ in order to have a pattern of topological charges for contact points with an enlarged periodicity. Actually, the periodicity of the Bloch Hamiltonian and of the pattern of topological charges need not be the same (the unit cell in reciprocal space is four times larger than the BZ for the Bloch Hamiltonian, see Fig. 9, and twice larger for the charges, see Fig. 10). There is an extra thing that is needed. Ahn et al. have identified it as non-orientability of the real Bloch Hamiltonian Ahn2019 . This has to do with the property of a real Bloch Hamiltonian under translation by a reciprocal lattice vector ${\bf G}$. If we write $H({\bf k}+{\bf G})=O_{\bf G}H({\bf k})O_{\bf G}^{-1}$, with $O_{\bf G}$ an orthogonal matrix, then if $\det O_{\bf G}=-1$ the model is said to be non-orientable (and orientable otherwise). For the Lieb-kagomé model in basis II and for ${\bf G}=2\pi\mathbf{u}_{x}$, one has

[TABLE]

with the orthogonal matrix $O_{x}=\text{diag}(-1,1,1)$. The fact that $\det O_{x}=-1$ means that it is non-orientable. This is also true of the translation by the reciprocal lattice vector $2\pi\mathbf{u}_{y}$:

[TABLE]

with the orthogonal matrix $O_{y}=\text{diag}(1,1,-1)$ such that $\det O_{y}=-1$. However for the translation by ${\bf G}=(2\pi,2\pi)$, we find that

[TABLE]

with the orthogonal matrix $O_{x+y}=\text{diag}(-1,1,-1)$ such that $\det O_{x+y}=1$.

Such an analysis in terms of a real Bloch Hamiltonian does not apply to the case of the generalized Lieb-kagomé model of Appendix C. In this case, one has to use the winding vector concept.