Edge states in super honeycomb structures with ${\mathcal{P}}{\mathcal{T}}$ symmetric deformations

The existence of interface conducting states is one of the most significant properties of topological insulators [18, 20]. It originates in certain energy gaps caused by symmetry breaking in the bulk. In the past two decades, there have been many efforts made to reveal the underlying mechanism of interface conducting states, mainly on the edge states caused by ${\mathcal{P}}{\mathcal{T}}$ symmetry breaking [14, 13, 11, 19]. However, physicists also find edge states in two-dimensional ${\mathcal{P}}{\mathcal{T}}$ symmetric materials [29]. In this paper, we model such phenomena by analyzing the spectral properties of a domain wall modulated Schrödinger operator.

One of the typical lattices whose deformed bulk has interface conducting states is the honeycomb lattice. The honeycomb lattice structure obsesses single Dirac cones at $K$ and $K^{\prime}$ points and chiral edge states after time-reversal symmetry breaking [1, 14, 13, 9, 11]. The most famous paradigm is graphene, a two-dimensional topological material [18]. Around 2015, a new way to deform the honeycomb lattice has been reported [29, 30, 31]. They considered the lattice in a supercell so that it has folding symmetry - caused by smaller periods than the supercell lattice. We call it a super honeycomb lattice in our previous work [8]. This new highly symmetric structure has fourfold degeneracy and a double Dirac cone at the $\Gamma$ point. Deforming the super honeycomb lattice in two different directions makes the energy bands open a gap with different topologies near the $\Gamma$ point [8, 22]. Physicists found a pair of gapped pseudospin edge states when connecting such two types of deformed materials, which are analogs of the helical edge states [24]. The propagation of electromagnetic waves with frequencies between these two edge states is well confined near the interface [7, 26]. Despite the broad applications of this phenomenon, it is rarely studied rigorously. Motivated by this, we consider edge operators interpolating between two kinds of perturbed operators across a rational edge and try to establish the existence of two gapped edge states by rigorous analysis.

Many models are used to develop analysis on edge states. Ammari and his collaborators studied the related mechanism on the subwavelength scale by considering Helmholtz problems [2, 3]. Bal and his collaborators have analyzed properties and topological descriptions of edge states in Dirac operators [4, 5]. Also, some numerical methods are introduced to associated problems [16, 17].

Among all the models, the one particle non-relativistic Schrödinger equation is one of the most effective models in illustrating edge conducting states [13, 21, 10, 11, 17]. Fefferman and Weinstein have laid solid foundations of rigorous analysis on such equations in both ${\mathcal{P}}$ and ${\mathcal{T}}$ breaking case [14, 12, 13]. In this paper, we consider the following two-dimensional Schrödinger operator:

$${\mathcal{H}}^{\delta}_{edge}=-\Delta+V({\bf x})+\delta\eta(\delta\bm{l}_{2}\cdot{\bf x})W({\bf x}),$$

with $V({\bf x})$ a super honeycomb lattice potential, $\eta(\zeta)$ a domain function, and $W({\bf x})$ a ${\mathcal{P}}{\mathcal{T}}$ symmetry preseving perturbation; see section 2.3. The limiting perturbed bulk operators on two sides:

$${\mathcal{H}}^{\delta}_{\pm}=-\Delta+V({\bf x})\pm\delta W({\bf x})$$

are ${\mathcal{P}}{\mathcal{T}}$ symmetric, but the folding symmetry is missing. Different from the analysis of ${\mathcal{P}}$ or ${\mathcal{T}}$ breaking bulk operators [13, 21, 11], we have to deal with the double Dirac cone [8], or two tangent single Dirac cones on the bands of the unperturbed operator ${\mathcal{H}}_{V}=-\Delta+V({\bf x})$, where a nontrivial second-order degeneracy is hidden; see section 3 . This means that only in higher-order terms can we clearly see the interaction between the four branches when discussing edge states. Such behavior is deeply rooted in the ${\mathcal{P}}{\mathcal{T}}$ symmetry preserving property. In the time-reversal symmetry breaking case, a one-dimensional Dirac operator can be obtained as the effective model by asymptotic expansions. An edge state with energy near the Dirac point can be derived from its exponentially decaying zero-energy states [21, 10, 11]. However, when it comes to the time reversal symmetric case, there are a pair of paralleling Dirac operators and therefore, two indistinguishable zero-energy states. Thus, the energy gap between two edge states is of higher order and deserves precise description; see section 4.

Another exciting perspective of understanding the edge states is the interplay between the symmetries and topology. Similar to how the Chern numbers characterize the topology of quantum materials, some topological indices are introduced to characterize the topology of bulk and edge Hamiltonians and the bulk-edge correspondence [23, 25]. There exist results on some one-dimensional models [27, 28] and higher-dimensional Dirac operators [4, 5, 6]. Different from a quantity of the whole band, in this paper, we provide topological intepretations concentrated on the $\Gamma$ point by calculating the local topological charge and analyzing the parities of the eigenstates of limiting bulk operators ${\mathcal{H}}^{\delta}_{\pm}$ on two sides of the edge, which connects the symmetries and the topology more explicitly; see Section 6.2.

This paper is organized as follows. Section 2 introduces the problem and our model – the domain wall modulated edge operator. Results on super honeycomb lattice potentials and double Dirac cones are briefly reviewed in this section. Section 3 establishes the approximations of the eigenstates near the double Dirac cone. Such near-energy approximations play essential roles in calculating the two gapped edge states asymptotically and rigorously proving their existence. The following three sections focus on the main conclusions obtained from our model. Section 4 calculates two edge states explicitly by multiscale expansions. Section 5 proves the existence of two gapped edge states rigorously. Section 6 provides some physical interpretations. Subsection 6.1 is about the effective local topological charge. Subsection 6.2 is about parties of eigenmodes. Section 7 gives numerical simulation of a typical example.

• •

  The lattice and dual lattice: $\bf{U}={\mathbb{Z}}{\bf u}_{1}\oplus{\mathbb{Z}}{\bf u}_{2}$ denotes the parallelogram lattice in ${\mathbb{R}}^{2}$ expanding by ${\bf u}_{1}$ and ${\bf u}_{2}$. $\Omega=\big{\{}{\bf u}=c_{1}{\bf u}_{1}+c_{2}{\bf u}_{2}:c_{1},c_{2}\in(-1/2,1/2)\big{\}}$ denotes its fundamental cell. ${\bf U}^{*}={\mathbb{Z}}{\bf k}_{1}\oplus{\mathbb{Z}}{\bf k}_{2}$ denotes its dual lattice with ${\bf k}_{l}\cdot{\bf u}_{j}=2\pi\delta_{l,j}$. $\Omega^{*}={\mathbb{R}}^{2}/{\bf U}^{*}$ is the fundamental cell of its dual lattice and is called the Brillouin Zone.

• •

  Symmetry operators:

  1. 1.

     $\frac{2}{3}\pi$ rotation operator ${\mathcal{R}}$: ${\mathcal{R}}[f]({\bf x})=f(R^{*}{\bf x})$ with $R^{*}=\begin{pmatrix}-\frac{1}{2}&-\frac{\sqrt{3}}{2}\\ \frac{\sqrt{3}}{2}&-\frac{1}{2}\end{pmatrix}$;

  2. 2.

     reflection operator ${\mathcal{P}}$: ${\mathcal{P}}[f]({\bf x})=f(-{\bf x})$;

  3. 3.

     time reversal operator ${\mathcal{T}}$: ${\mathcal{T}}[f]({\bf x})=\overline{f({\bf x})}$.

• •

  The folding lattice and operator: when ${\bf u}_{2}=-R^{*}{\bf u}_{1}$, ${\bf v}_{1}$ and ${\bf v}_{2}$ as following are the periods of the folding lattice.

  $${\bf v}_{1}=\frac{1}{3}(2{\bf u}_{1}-{\bf u}_{2}),{\qquad}{\bf v}_{2}=\frac{1}{3}({\bf u}_{1}+{\bf u}_{2}).$$

  The associated folding operators are: ${\mathcal{V}}_{l}[f]({\bf x})=f(x+{\bf v}_{l})$.

• •

  The rational edge: $\bm{w}_{1}=a_{1}{\bf u}_{1}+b_{1}{\bf u}_{2}$; $\bm{w}_{2}=a_{2}{\bf u}_{1}+b_{2}{\bf u}_{2}$. $a_{1},b_{1},a_{2},b_{2}$ are integers and $a_{1}b_{2}-a_{2}b_{1}=1$. $\bm{l}_{1}=b_{2}{\bf k}_{1}-a_{2}{\bf k}_{2}$ and $\bm{l}_{2}=-b_{1}{\bf k}_{1}+a_{1}{\bf k}_{2}$ satisfy $\bm{w}_{m}\cdot\bm{l}_{j}=2\pi\delta_{m,j}$. $\bm{w}_{1}$ is the direction of the rational edge. $\bm{l}_{2}$ is its dual direction. $\tilde{\bm{l}}_{1}=\bm{l}_{1}-\frac{\bm{l}_{1}\cdot\bm{l_{2}}}{\|\bm{l}_{2}\|^{2}}{\bm{l}_{2}}$ is orthogonal to $\bm{l}_{2}$ and satisfies $\tilde{\bm{l}}_{1}\cdot\bm{w}_{1}=2\pi$.

• •

  $\Omega_{e}=\big{\{}{\bf u}=c_{1}\bm{w}_{1}+c_{2}\bm{w}_{2}:c_{1}\in(-1/2,1/2),c_{2}\in{\mathbb{R}}\big{\}}$.

• •

  There are three different inner products:

  1. 1.

     $\langle f({\bf x}),g({\bf x})\rangle_{L^{2}(\Omega)}=\int_{\Omega}\overline{f({\bf x})}g({\bf x})d{\bf x}$;

  2. 2.

     $\langle f({\bf x}),g({\bf x})\rangle_{L^{2}(\Omega_{e})}=\int_{\Omega_{e}}\overline{f({\bf x})}g({\bf x})d{\bf x}$;

  3. 3.

     $\langle f(\zeta),g(\zeta)\rangle_{L^{2}({\mathbb{R}})}=\int_{{\mathbb{R}}}\overline{f(\zeta)}g(\zeta)d\zeta$.

  The first one is mainly used in the bulk problem, and the last two are used in the edge problem.

• •

  There are two particular function spaces often used corresponding to bulk and edge problem respectively:

  $$\chi=\big{\{}f\in L^{2}_{loc}({\mathbb{R}}^{2}):f({\bf x}+{\bf u}_{l})=f({\bf x}),l=1,2\big{\}};$$

  $$\chi_{e}=\big{\{}f\in L^{2}_{loc}({\mathbb{R}}^{2}):f({\bf x}+\bm{w}_{1})=f({\bf x}),\int_{\Omega_{e}}|f({\bf x})|^{2}d{\bf x}<\infty\big{\}}.$$

• •

  Pauli matrices:

  $$\sigma_{1}=\begin{pmatrix}0&\hskip 5.69046pt1\\ 1&\hskip 5.69046pt0\end{pmatrix};\quad\sigma_{2}=\begin{pmatrix}0&\hskip 2.84544pt-i\\ i&\hskip 2.84544pt0\end{pmatrix};{\quad}\sigma_{3}=\begin{pmatrix}1&\hskip 2.84544pt0\\ 0&\hskip 2.84544pt-1\end{pmatrix}.$$

In this subsection, we review the symmetries of super honeycomb lattice potentials. It originates in viewing the honeycomb lattice in a supercell. Thus, the super honeycomb lattice possesses a folding symmetry, which results in fourfold degeneracy at the $\Gamma$ point on energy bands.

Fefferman and Weinstein have summarized the structures of the honeycomb lattice in their paper [14]. It is such a periodic structure in ${\mathbb{R}}^{2}$:

• •

  $\frac{2}{3}\pi$-rotation symmetric;

• •

  ${\mathcal{P}}{\mathcal{T}}$ (parity and time reversal) symmetric.

The corresponding transformations in $L^{2}_{loc}({\mathbb{R}}^{2})$ are:

• •

  $\frac{2}{3}\pi$ rotation operator ${\mathcal{R}}$: ${\mathcal{R}}[f]({\bf x})=f(R^{*}{\bf x})$ with $R^{*}=\begin{pmatrix}-\frac{1}{2}&-\frac{\sqrt{3}}{2}\\ \frac{\sqrt{3}}{2}&-\frac{1}{2}\end{pmatrix}$;

• •

  reflection operator ${\mathcal{P}}$: ${\mathcal{P}}[f]({\bf x})=f(-{\bf x})$;

• •

  time reversal operator ${\mathcal{T}}$: ${\mathcal{T}}[f]({\bf x})=\overline{f({\bf x})}$.

Consider the parallelogram lattice $\bf{U}={\mathbb{Z}}{\bf u}_{1}\oplus{\mathbb{Z}}{\bf u}_{2}$. The corresponding unit cell is

$$\Omega=\big{\{}{\bf u}=c_{1}{\bf u}_{1}+c_{2}{\bf u}_{2},{\quad}c_{1},c_{2}\in(-1/2,1/2)\big{\}}.$$

(2.1)

Denote its dual lattice and dual unit cell $\bf{U}^{*}={\mathbb{Z}}{\bf k}_{1}\oplus{\mathbb{Z}}{\bf k}_{2}$ and $\Omega^{*}={\mathbb{R}}^{2}/\bf{U}^{*}$, where ${\bf k}_{j}$ satisfies ${\bf u}_{l}\cdot{\bf k}_{j}=2\pi\delta_{l,j}$.

The relation and differences between honeycomb lattice and super honeycomb lattice are shown explicitly in Figure 2. Note that the parallelogram lattice discussed above and the hexagonal lattice in the pictures are equivalent. Compared with the honeycomb lattice potentials, the super honeycomb lattice potential obsesses smaller periods and, therefore, a folding symmetry. Namely, the smaller periods are

$${\bf v}_{1}=\frac{1}{3}(2{\bf u}_{1}-{\bf u}_{2}),{\qquad}{\bf v}_{2}=\frac{1}{3}({\bf u}_{1}+{\bf u}_{2}).$$

(2.2)

For simiplicity, denote translation operators: ${\mathcal{V}}_{l}[f]({\bf x})=f(x+{\bf v}_{l})$. Let ${\bf q}_{j}$ be the dual vectors of ${\bf v}_{l}$, i.e., ${\bf q}_{j}\cdot{\bf v}_{l}=2\pi\delta_{j,l}$.

The super honeycomb lattice potential is defined below.

The bulk operators ${\mathcal{H}}_{\pm}^{\delta}=-\Delta+V({\bf x})\pm\delta W({\bf x})$ on two sides of the edge are deformed from the highly symmetric operator ${\mathcal{H}}_{V}=-\Delta+V({\bf x})$, where $V({\bf x})$ is a super honeycomb lattice potential. In this subsection, we give two conclusions about ${\mathcal{H}}_{V}$ and ${\mathcal{H}}^{\delta}=-\Delta+V({\bf x})+\delta W({\bf x})$ [8].

Generally, for a ${\bf u}_{1}$ and ${\bf u}_{2}$ doubly periodic elliptic operator ${\mathcal{H}}$, by Floquet-Bloch theorem, its spectrum can be decomposed into momentum space $\Omega^{*}$:

$$\sigma_{L^{2}({\mathbb{R}}^{2})}({\mathcal{H}})=\bigcup_{{\bf k}\in\Omega^{*}}\sigma_{L^{2}_{{\bf k}}({\mathbb{R}}^{2}/{\bf U})}({\mathcal{H}}|_{L^{2}_{{\bf k}}({\mathbb{R}}^{2}/{\bf U})}).$$

Here the ${\bf k}$-momentum function space is the Hilbert space

$$L^{2}_{{\bf k}}({\mathbb{R}}^{2}/{\bf U})=\big{\{}f\in L^{2}_{loc}({\mathbb{R}}^{2}):f({\bf x}+{\bf u})=e^{i{\bf k}\cdot{\bf u}}f({\bf x}),{\bf x}\in{\mathbb{R}}^{2},{\bf u}\in{\bf U};f({\bf x})\in L^{2}(\Omega)\big{\}}$$

equipped with the inner product:

$$\langle f({\bf x}),g({\bf x})\rangle_{L^{2}(\Omega)}=\int_{\Omega}\overline{f({\bf x})}g({\bf x})d{\bf x}.$$

(2.4)

For our bulk operators, $\sigma_{L^{2}_{{\bf k}}({\mathbb{R}}^{2}/{\bf U})}({\mathcal{H}}|_{L^{2}_{{\bf k}}({\mathbb{R}}^{2}/{\bf U})})$ is a lower bounded discrete set:

$$E_{1}({\bf k})\leq E_{2}({\bf k})\leq E_{3}({\bf k})\leq...$$

(2.5)

${\mathcal{S}}_{n}=\big{\{}\bigl{(}{\bf k},E_{n}({\bf k})\bigr{)}:{\bf k}\in\Omega^{*}\big{\}}$ is the $n^{th}$ energy surface of ${\mathcal{H}}$. Dirac cone describes the phenomenon of conical touch between energy surfaces. It is a kind of degeneracy and singularity on energy surfaces.

The first conclusion is that ${\mathcal{H}}_{V}$ has fourfold degeneracy at the $\Gamma$ point and a double Dirac cone for general super honeycomb lattice potential $V({\bf x})$ due to the symmetry. The key is that the space of ${\bf u}_{1}$ and ${\bf u}_{2}$ periodic functions $\chi=L^{2}_{\mathbf{0}}({\mathbb{R}}^{2}/{\bf U})$ has the decomposition:

$$\chi=\bigoplus_{\xi_{1},\xi_{2}=1,\tau,\overline{\tau}}\chi_{\xi_{1},\xi_{2}},{\quad}\tau=e^{\frac{2}{3}\pi i}.$$

(2.6)

$\chi_{\xi_{1},\xi_{2}}$ is the intersection of the characteristic subspace of ${\mathcal{V}}_{1}$ with eigenvalue $\xi_{1}$ and the characteristic subspace of ${\mathcal{R}}$ with eigenvalue $\xi_{2}$. Because any super honeycomb lattice potential $V({\bf x})$ is in $\chi_{1,1}$, ${\mathcal{H}}_{V}$ has the rotation and folding symmetries. Thus, each $\chi_{\xi_{1},\xi_{2}}$ is an invariant subspace of ${\mathcal{H}}_{V}$. The following theorem shows the existence of fourfold degeneracy and its relationship with the above function space decomposition.

Equation (2.7) indicates the existence of a fourfold degeneracy at the $\Gamma$ point. (2.8) provides insight into the four eigenfunctions corresponding to this degeneracy, revealing that they can be related through symmetries. We will give the result of the double Dirac cone later after we calculate the second-order bifurcation terms; see section 3.1.

The second conclusion is that the energy surfaces of ${\mathcal{H}}^{\delta}$ open a gap of order $O(\delta)$ at the $\Gamma$ point when W(x) is a folding symmetry breaking potential as in the following definition.

The second conclusion is stated precisely as follows.

In this paper, we consider two kinds of bulks ${\mathcal{H}}^{\delta}_{\pm}=-\Delta+V({\bf x})\pm\delta W({\bf x})$ connected along a rational edge ${\mathbb{R}}\bm{w}_{1}$. Here $\bm{w}_{1}=a_{1}{\bf u}_{1}+b_{1}{\bf u}_{2}$, where $a_{1}$ and $b_{1}$ are relatively prime integers. Then there exist relatively prime integers $a_{2},b_{2}$ such that $a_{1}b_{2}-a_{2}b_{1}=1$. Let $\bm{w}_{2}=a_{2}{\bf u}_{1}+b_{2}{\bf u}_{2}$. Take $\bm{l}_{1}=b_{2}{\bf k}_{1}-a_{2}{\bf k}_{2}$ and $\bm{l}_{2}=-b_{1}{\bf k}_{1}+a_{1}{\bf k}_{2}$. Then $\bm{w}_{j}\cdot\bm{l}_{l}=2\pi\delta_{j,l}$. Note that $\bm{w}_{1}$ and $\bm{w}_{2}$ expands a cell equivalent to $\Omega$ and $\bm{l}_{1}$ and $\bm{l}_{2}$ expands a cell equivalent to $\Omega^{*}$. Let $\Omega_{e}$ denote such an area around the edge:

$$\Omega_{e}=\{{\bf x}\in{\mathbb{R}}^{2}\mid{\bf x}=p\bm{w}_{1}+q\bm{w}_{2},p\in(-1/2,1/2),q\in{\mathbb{R}}\}.$$

For the rational edge ${\mathbb{R}}\bm{w}_{1}$, we use the following asymptotic domain wall modulated operator:

$${\mathcal{H}}_{edge}^{\delta}=-\Delta+V({\bf x})+\delta\eta(\delta\bm{l}_{2}\cdot{\bf x})W({\bf x}).$$

(2.10)

The limiting bulk operators ${\mathcal{H}}^{\delta}_{\pm}$ have been introduced in the last subsection. The function $\eta(\delta\bm{l}_{2}\cdot{\bf x})$ is the domain wall modulation function. It is flat along the edge direction $\bm{w}_{1}$ and varies slowly in the direction $\bm{w}_{2}$. $\zeta=\delta\bm{l}_{2}\cdot{\bf x}$ is called the slow variable. The rigorous definition of a domain wall function $\eta(\zeta)$ is below.

Similar to that in bulk case, we have:

$$\sigma_{L^{2}({\mathbb{R}}^{2})}({\mathcal{H}}_{edge}^{\delta})=\bigcup_{k_{\parallel}\in[-\pi,\pi)}\sigma_{L^{2}_{k_{\parallel}}({\mathbb{R}}^{2}/{\mathbb{Z}}\bm{w}_{1})}({\mathcal{H}}_{edge}^{\delta}|_{L^{2}_{k_{\parallel}}({\mathbb{R}}^{2}/{\mathbb{Z}}\bm{w}_{1})}).$$

Here the function space is one-dimensional quasi-periodic:

$$L^{2}_{k_{\parallel}}({\mathbb{R}}^{2}/{\mathbb{Z}}\bm{w}_{1})=\big{\{}f\in L^{2}_{loc}({\mathbb{R}}^{2}):f({\bf x}+\bm{w}_{1})=e^{ik_{\parallel}}f({\bf x}),{\bf x}\in{\mathbb{R}}^{2};f({\bf x})\in L^{2}(\Omega_{e})\big{\}}.$$

The inner product on $L^{2}_{k_{\parallel}}({\mathbb{R}}^{2}/{\mathbb{Z}}\bm{w}_{1})$ is :

$$\langle f,g\rangle_{L^{2}(\Omega_{e})}=\int_{\Omega_{e}}\overline{f({\bf x})}g({\bf x})d{\bf x}.$$

(2.11)

We summarize our asymptotic model in the following definition.

Our aim is to solve the following eigenvalue problem near $k_{\parallel}=0$:

$${\mathcal{H}}_{edge}^{\delta}\psi({\bf x},k_{\parallel})={\mathcal{E}}(k_{\parallel})\psi({\bf x},k_{\parallel}),{\quad}\psi({\bf x},k_{\parallel})\in L^{2}_{k_{\parallel}}({\mathbb{R}}^{2}/{\mathbb{Z}}\bm{w}_{1}).$$

(2.13)

The following sections contribute to proving the existence of two gapped edge states near $k_{\parallel}=0$ and characterizing these two edge states precisely.

This subsection discusses the bifurcation matrices for any ${\bf k}$ sufficiently small to estimate the near-energy eigenstates.

The eigenvalue problem of ${\mathcal{H}}_{V}$ on $L^{2}_{{\bf k}}({\mathbb{R}}^{2}/{\bf U})$:

	$\displaystyle{\mathcal{H}}_{V}\phi({\bf x};{\bf k})$	$\displaystyle=E({\bf k})\phi({\bf x};{\bf k}),$
	$\displaystyle\phi({\bf x};{\bf k})=e$	$\displaystyle{}^{i{\bf k}\cdot{\bf x}}p({\bf x}),\quad p({\bf x})\in\chi;$

is equivalent to

$${\mathcal{H}}_{V}({\bf k})p({\bf x};{\bf k})=E({\bf k})p({\bf x};{\bf k}),\quad p({\bf x})\in\chi,$$

where ${\mathcal{H}}_{V}({\bf k})=-(\nabla+i{\bf k})^{2}+V({\bf x})$.

We only consider the near-energy solution:

$$E({\bf k})=E_{{}_{D}}+\mu,\quad p({\bf x};{\bf k})={\Phi}({\bf x})^{\mathrm{T}}\bf{P}({\bf k})+\psi({\bf x};{\bf k}).$$

where $\bf{P}({\bf k})$ $=(p_{1}({\bf k}),p_{2}({\bf k}),p_{3}({\bf k}),p_{4}({\bf k}))^{\mathrm{T}}$, $\psi({\bf x};{\bf k})\in\operatorname{Ker}({\mathcal{H}}_{V}-E_{{}_{D}})^{\perp}$, the corrector $\mu$ and $\psi({\bf x};{\bf k})$ are small, and

$${\Phi}({\bf x})=(\phi_{1}({\bf x}),\phi_{2}({\bf x}),\phi_{3}({\bf x}),\phi_{4}({\bf x}))^{\mathrm{T}}.$$

(3.1)

Here $\{\phi_{j}({\bf x})\}_{j}$ span $\operatorname{Ker}({\mathcal{H}}_{V}-E_{{}_{D}})$ as in (2.8). Substitute these into the eigenvalue problem and rearrange it into:

$$({\mathcal{H}}_{V}-E_{{}_{D}})\psi({\bf x};{\bf k})=(\mu+2i{\bf k}\cdot\nabla-\|{\bf k}\|^{2})\bigl{(}{\Phi}({\bf x})^{\mathrm{T}}{\bf{P}}({\bf k})+\psi({\bf x};{\bf k})\bigr{)}.$$

(3.2)

Thus, the corrector $\psi({\bf x};{\bf k})$ is the solution to the problem:

	$\displaystyle\bigl{(}I-({\mathcal{H}}_{V}-E_{{}_{D}})^{-1}{\mathcal{Q}}_{\perp}(2i{\bf k}\cdot\nabla$	$\displaystyle+\|{\bf k}\|^{2}+\mu)\bigr{)}\psi({\bf x};{\bf k})$
		$\displaystyle=({\mathcal{H}}_{V}-E_{{}_{D}})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla{\Phi}({\bf x})^{\mathrm{T}}\bf{P}({\bf k}),$

where ${\mathcal{Q}}_{\perp}$ is the projection map to the orthogonal complement space of $\operatorname{Ker}({\mathcal{H}}_{V}-E_{{}_{D}})$ and

$$({\mathcal{H}}_{V}-E{{}_{D}})^{-1}:{\mathcal{Q}}_{\perp}\chi\to{\mathcal{Q}}_{\perp}H^{1}({\mathbb{R}}^{2}/{\bf U}).$$

(3.3)

Here $H^{1}({\mathbb{R}}^{2}/{\bf U})$ is the limitation of $\chi$ in $H^{1}({\mathbb{R}}^{2})$. Therefore,

$$\psi({\bf x};{\bf k})=\bigl{(}({\mathcal{H}}_{V}-E_{{}_{D}})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla{\Phi}({\bf x})^{\mathrm{T}}+O(\|{\bf k}\|^{2})\bigr{)}\bf{P}({\bf k}).$$

Besides, from (3.2), we know that:

$$\langle\phi_{l}({\bf x}),(\mu+2i{\bf k}\cdot\nabla-\|{\bf k}\|^{2})\bigl{(}{\Phi}({\bf x})^{\mathrm{T}}{\bf{P}}({\bf k})+\psi({\bf x};{\bf k})\bigr{)}\rangle=0.$$

Based on all these, solving the original eigenvalue problem can be reduced to solving possible $\bf{P}({\bf k})$ in $(\mu I+B({\bf k})){\bf{P}}({\bf k})=\bf{0}$, with

$$B({\bf k})=B_{1}({\bf k})+B_{2}({\bf k})+O(\|{\bf k}\|)^{3}.$$

Here

	$\displaystyle B_{1}({\bf k})$	$\displaystyle=\biggl{(}\big{\langle}\phi_{l}({\bf x}),2i{\bf k}\cdot\nabla\phi_{j}({\bf x})\big{\rangle}_{L^{2}(\Omega)}\biggr{)}_{l,j};$
	$\displaystyle B_{2}({\bf k})$	$\displaystyle=\biggl{(}\big{\langle}\phi_{l}({\bf x}),\bigl{(}2i{\bf k}\cdot\nabla({\mathcal{H}}_{V}-E_{D})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla-\|{\bf k}\|^{2}\bigr{)}\phi_{j}({\bf x})\big{\rangle}_{L^{2}(\Omega)}\biggr{)}_{l,j}.$

$B({\bf k})$ represents the bifurcation matrix.

1. 1.

   First-order bifurcation matrix $B_{1}({\bf k})$.

   Taking advatange of the symmetry relations between $\phi_{j}({\bf x})$ in (2.8), we can obtain [8]

   $$B_{1}({\bf k})=\begin{pmatrix}0&2i{\bf k}\cdot\bm{v}_{\sharp}&0&0\\ \overline{2i{\bf k}\cdot\bm{v}_{\sharp}}&0&0&0\\ 0&0&0&-2i{\bf k}\cdot\bm{v}_{\sharp}\\ 0&0&-\overline{2i{\bf k}\cdot\bm{v}_{\sharp}}&0\end{pmatrix}.$$

   Here $\bm{v}_{\sharp}$ is such a quantity:

   $$\bm{v}_{\sharp}=\langle\phi_{1}({\bf x}),\nabla\overline{\phi_{1}(-{\bf x})}\rangle_{L^{2}(\Omega)}=\frac{v_{{}_{F}}}{2}\begin{pmatrix}1\\ i\end{pmatrix}e^{i\theta_{\sharp}}.$$

   (3.4)

   Remark 3.1.

   The eigenvalues of $b({\bf k})$ are $\lambda_{\pm}=\pm|2i{\bf k}\cdot\bm{v}_{\sharp}|=\pm v_{{}_{F}}\|{\bf k}\|$. Both $\lambda_{\pm}$ are of multiplicity two, which means that the upper two bands near the fourfold degenerate Dirac point are tangent to each other and so do the lower two bands. To distinguish the tangent bands, we have to calculate higher-order approximations. The Fermi velocity $v_{{}_{F}}$ is the slope of the double Dirac cone. It is invariant under phase transformation of $\phi_{1}({\bf x})$. The corresponding eigenvectors are determined by $\theta_{\sharp}$, which naturally changes with the phase transformation of $\phi_{1}({\bf x})$.

2. 2.

   Second-order bifurcation matrix $B_{2}({\bf k})$.

   The second-order bifurcation matrix is Hermitian. Thus, it is enough to discuss the upper part of

   $$\biggl{(}\big{\langle}\phi_{l}({\bf x}),\bigl{(}2i{\bf k}\cdot\nabla({\mathcal{H}}_{V}-E_{D})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla\bigr{)}\phi_{j}({\bf x})\big{\rangle}_{L^{2}(\Omega)}\biggr{)}_{l,j}.$$

   First, due to

   		$\displaystyle\langle\phi_{l}({\bf x}),2i{\bf k}\cdot\nabla({\mathcal{H}}_{V}-E_{D})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla\phi_{j}({\bf x})\rangle_{L^{2}(\Omega)}$
   	$\displaystyle=$	$\displaystyle\langle{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla\phi_{l}({\bf x}),({\mathcal{H}}_{V}-E_{D})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla\phi_{j}({\bf x})\rangle_{L^{2}(\Omega)},$

   and the fact that $({\mathcal{H}}_{V}-E_{D})^{-1}$ is an elliptic and Hermitian operator on the periodic function space ${\mathcal{Q}}_{\perp}L^{2}(\Omega)$, $\langle\phi_{l}({\bf x}),...\phi_{l}({\bf x})\rangle_{L^{2}(\Omega)}$ are real for all $l$. Besides, using the symmetries between the eigenfunctions in Theorem 2.4, $\langle\phi_{l}({\bf x}),...\phi_{l}({\bf x})\rangle_{L^{2}(\Omega)}$ are the same for all $l$. This means the diagonal elements of the two-order bifurcation matrix are just real multiples of the identity. For other elements, because the translation operator ${\mathcal{V}}_{1}$ is a unitary operator and commutative with both $\nabla$ and $({\mathcal{H}}_{V}-E_{D})^{-1}$ on $L^{2}(\Omega)$, we have

   $\displaystyle\langle\phi_{l}({\bf x}),...\phi_{j}({\bf x})\rangle_{L^{2}(\Omega)}=\langle{\mathcal{V}}_{1}\phi_{l}({\bf x}),...{\mathcal{V}}_{1}\phi_{j}({\bf x})\rangle_{L^{2}(\Omega)}.$

   Because all $\phi_{l}({\bf x})$ are eigenfunctions of ${\mathcal{V}}_{1}$, it turns out that the two-order bifurcation matrix is quasi-diagonal:

   	$\displaystyle 0=$	$\displaystyle\langle\phi_{1}({\bf x}),...\phi_{3}({\bf x})\rangle_{L^{2}(\Omega)}=\langle\phi_{1}({\bf x}),...\phi_{4}({\bf x})\rangle_{L^{2}(\Omega)}$
   	$\displaystyle=$	$\displaystyle\langle\phi_{2}({\bf x}),...\phi_{3}({\bf x})\rangle_{L^{2}(\Omega)}=\langle\phi_{2}({\bf x}),...\phi_{4}({\bf x})\rangle_{L^{2}(\Omega)}.$

   Therefore, the second order bifurcation matrix $B_{2}({\bf k})$ is equal to:

   $$\begin{pmatrix}m({\bf k})-\|{\bf k}\|^{2}&b({\bf k})&0&0\\ \overline{b({\bf k})}&m({\bf k})-\|{\bf k}\|^{2}&0&0\\ 0&0&m({\bf k})-\|{\bf k}\|^{2}&b({\bf k})\\ 0&0&\overline{b({\bf k})}&m({\bf k})-\|{\bf k}\|^{2}\end{pmatrix}.$$

   where

   $$m({\bf k})=\langle\phi_{1}({\bf x}),2i{\bf k}\cdot\nabla({\mathcal{H}}_{V}-E_{D})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla\phi_{1}({\bf x})\rangle_{L^{2}(\Omega)};$$

   (3.5)

   $$b({\bf k})=\langle\phi_{1}({\bf x}),2i{\bf k}\cdot\nabla({\mathcal{H}}_{V}-E_{D})^{-1}{\mathcal{Q}}_{\perp}2i{\bf k}\cdot\nabla\overline{\phi_{1}(-{\bf x})}\rangle_{L^{2}(\Omega)}.$$

   (3.6)

   Remark 3.2.

   $m({\bf k})$ is real. The diagonal part $\bigl{(}m({\bf k})-\|{\bf k}\|^{2}\bigr{)}I$ of this matrix does no contribution to bifurcation. The term resulting in second-order bifurcation of the upper or lower two bands is $b({\bf k})$. The form of $b({\bf k})$ is calculated explicitly in Appendix A.