text:
\section{Introduction and Summary of Results}

\subsection{Introduction}\label{intro}

Motivated by the novel and subtle properties of electronic waves in graphene \cites{geim2007rise}, there has been very wide interest in the propagation of waves in two-dimensional structures having the symmetries of a hexagonal tiling of the plane and its applications to electromagnetic and other types of waves; see, for example, \cites{Singha_11,Rechtsman-etal:13,Shvets-PTI:13,MKW:15,yang2015topological,wu2015scheme}. Such physical systems have been dubbed artificial graphene. The present article is motivated by {\it photonic graphene}; the propagation of waves governed by the two-dimensional Maxwell equations in honeycomb media.

Among the remarkable properties of graphene (actual and artificial) being intensively explored in the fundamental and applied scientific communities are the existence of Dirac points \cites{RMP-Graphene:09,novoselov2005two,wallace1947band,FW:12,FLW-CPAM:17,berkolaiko-comech:15} (conical points at the intersections between dispersion surfaces), the implied wave-packet (quasi-particle) dynamics of states which are
 spectrally localized near Dirac points \cites{Segev07prl,ablowitz2009conical,FW:14}, and topologically protected edge states \cites{HR:07,RH:08,RMP-Graphene:09,Delplace-etal:11,poo2011experimental,plotnik2013observation,Shvets:14,chen2014experimental,ma2015guiding,FLW-2DM:15,cheng2016robust,FLW-AnnalsPDE:16}. Conical singularities
 in dispersion relations and their consequences for wave propagation have a long history. For example, they are well-known to occur  for spatially homogeneous anisotropic Maxwell's equations; see
 \cite{Berry-Jeffrey:07} and references cited therein

The goal of this paper is to investigate the phenomena of Dirac points, effective dynamics of wavepackets and the bifurcation theory of topologically protected edge states in the context of a class of elliptic partial differential operators,
which incorporates important physical cases of electromagnetic propagation in honeycomb structures.
 %
Specifically we  consider a class of periodic {\it scalar} divergence-form operators:
\begin{equation}
{\mathcal{L}^A}=-\nabla\cdot A \nabla\ =\ -\sum_{i,j=1}^2 \frac{\partial}{\partial x_i} a_{ij}({\bf x}) \frac{\partial}{\partial x_j} \ ,\ \
A({\bf x})= \left(a_{ij}({\bf x})\right)_{i,j=1,2}\ ,
\label{L_def}
\end{equation}
with structural assumptions on $A({\bf x})$, which we now discuss.

 Introduce the following operations:  $\mathcal{C}[f]({\bf x})=\overline{f({\bf x})}, \ $
$\mathcal{P}[f]({\bf x})=f(-{\bf x})$,  and
$\mathcal{R}[f]({\bf x})\equiv f(R^*{\bf x})$, where $R$ is the $120^\circ$ clockwise rotation matrix of vectors in the plane.
  For the bulk (unperturbed) structure,  specified by a $2\times2$ matrix function $A({\bf x})$, we assume the following properties:
\begin{enumerate}
\item [(i)] $A({\bf x})$ is  periodic with respect to the equilateral triangular lattice, $\Lambda_h$; see Section \ref{triangle}.
\item [(ii)] ${\mathcal{L}^A}$ is self-adjoint and uniformly elliptic on $\mathbb{R}^2$,\\
 and  with respect to some origin of coordinates:
\item [(iii)] ${\mathcal{L}^A}$ is $\mathcal{PC}-$ invariant; $\left[\mathcal{PC},{\mathcal{L}^A}\right]=0$, and
\item [(iv)] ${\mathcal{L}^A}$ is $\mathcal{R}-$ invariant; $\left[\mathcal{R},{\mathcal{L}^A}\right]=0$, where $[\mathcal{L}_1,\mathcal{L}_2]=\mathcal{L}_1\mathcal{L}_2-\mathcal{L}_2\mathcal{L}_1$.
\end{enumerate}
\medskip

\noindent A characterization of matrix constitutive laws, $A({\bf x})$,  which satisfy assumptions (i)-(iv) is presented in Section \ref{honeycomb-structures}.
\medskip

In this article we obtain results on:
\begin{enumerate}
\item Dirac points (conical singularities) in the Floquet-Bloch dispersion surfaces associated with the eigenvalue problem ${\mathcal{L}^A}\psi=E\psi$; see \eqref{Eq_Eigen}.
\item The non-persistence of Dirac points (opening of a local spectral gap) for the perturbed eigenvalue problem for the operator $\mathcal{L}^{A+\delta B}$, where $\delta\ll1$ and  $B({\bf x})$ breaks $\mathcal{P}-$ or $\mathcal{C}-$ invariance, and the persistence of Dirac points for small $\mathcal{P}\mathcal{C}-$ invariant perturbations.
\item The bifurcation of  topologically protected edge states from Dirac points upon introduction of a (non-compact) domain wall modulated periodic perturbation, $\mathcal{L}_{\rm dw}^{\delta}$, of ${\mathcal{L}^A}$; see (d) in Section \ref{main_results}.
   \end{enumerate}

Two examples in 2D electromagnetics which motivate our study of this class of operators, and to which our results apply,  arise for Maxwell's equations considered
 in (1) {\it magneto-optic media} and  (2) in {\it bi-anisotropic meta-materials}.

\noindent {\bf Example 1:}\  In \cites{HR:07,RH:08}, Haldane and Raghu proposed a photonic analogue of the quantum Hall effect. In particular, they demonstrated that uni-directional protected edge states could propagate along a domain-wall for systems governed by Maxwell's equations in, for example, magneto-optic materials which break $\mathcal{C}-$ (complex-conjugation) invariance.
Such photonic edge states were experimentally observed by Wang {\it et. al.} \cites{Soljacic-etal:08}.
Our domain-wall modulated operator $\mathcal{L}_{\rm dw}^\delta$, as well as the class of Schr\"odinger operators considered
in \cites{FLW-PNAS:14,FLW-2DM:15,FLW-AnnalsPDE:16,FLW-MAMS:17}, models the physical setting of \cites{HR:07,RH:08}.
In this setting, the eigenfunction $\psi$
corresponds to $H_3$ in the case of a TE-polarized electromagnetic field:
 ${\bf E}=(E_1,E_2,0)$ and $\vec{H}=(0,0,H_3)$.
\begin{remark}[TM Polarization]
For TM-Polarization: ${\bf E}=(0,0,E_3)$ and ${\bf H}=(H_1,H_2,0)$, and
 $\psi=E_3$ satisfies a scalar eigenvalue problem for a Helmholtz operator. The techniques of this paper
 can be adapted to this case as well, but we focus on the divergence form operator ${\mathcal{L}^A}$.
\end{remark}

 \noindent {\bf Example 2:}\ In Khanikaev {\it et. al.} \cites{Shvets-PTI:13}, bi-anisotropic media which respect time-reversal symmetry but break $\mathcal{C}-$ invariance,  were studied and the counter-propagation of uni-directional states of opposite ``spin'', $\psi_\pm=H_3\pm E_3$, associated with an underlying Kane-Mele model, was computationally demonstrated.
 \medskip

\noindent  A detailed discussion of both examples is presented in Appendix \ref{maxwell}.

\subsection{Main results}\label{main_results}

We give a brief summary of our main results. Our approach makes use of and extends methods developed to treat the case of honeycomb Schr\"odinger operators \cites{FW:12,FW:14,FLW-PNAS:14,FLW-2DM:15,FLW-AnnalsPDE:16,FLW-MAMS:17,FLW-CPAM:17}.
\begin{enumerate}[(a)]
\item {\it Characterization of honeycomb structured media:}\ Theorem \ref{invariance} and Corollary \ref{A-expand} characterize 2D honeycomb structured media, defined by
 the Hermitian operator ${\mathcal{L}^A}=-\nabla\cdot A({\bf x})\nabla$ (see \eqref{L_def}), where $A({\bf x})$ satisfies conditions (i)-(iv) above.
 \item {\it Dirac points:}\ Theorem \ref{prop_conical} gives conditions for the existence of Dirac points at the vertices of the Brillouin zone,
 $\mathcal{B}_h$.
 Theorem \ref{low-dp} states that all sufficiently low-contrast honeycomb structures, for which a distinguished Fourier coefficient of the structure is non-zero, have Dirac points. Depending on the sign of this Fourier coefficient, these Dirac points occur either at intersections of the $1^{st}$ and $2^{nd}$ or $2^{nd}$ and $3^{rd}$ dispersion surfaces.
  Extension of these results to generic honeycomb structured media of arbitrary contrast, using the strategy of
\cites{FW:12} (see also Appendix D of \cites{FLW-MAMS:17}) is discussed in Section \ref{generic_eps}.
 In Figure \ref{E_k_BZ} we plot two dispersion surfaces for a honeycomb structured medium, $A({\bf x})$, which intersect at Dirac points for quasi-momenta located at the size vertices of the Brillouin zone.

 Remark \ref{wavepkts} discusses the effective dynamics of wave-packets for initial data which are spectrally localized near Dirac points. The envelope of such wave-packets evolves slowly and on long time scales according to an effective massless Dirac equation.

 \begin{figure}
\centering
\includegraphics[width=4in]{E_k_BZ_axis}
 \caption{\footnotesize
 Lowest two dispersion surfaces ${\bf k}\equiv(k^{(1)},k^{(2)})\in\mathcal{B}_h\mapsto E({\bf k})$ of the band structure of ${\mathcal{L}^A}$, with $A({\bf x})=a(x)I$ where $a({\bf x})$ is a particular honeycomb lattice function.
 Dirac points occur at the intersection of two dispersion surfaces, at the six vertices of the Brillouin zone, $\mathcal{B}_h$.}
 \label{E_k_BZ}
 \end{figure}

 \item {\it Non-persistence and persistence of Dirac points:}\ Theorem \ref{PT_thm} studies the instability of Dirac points under perturbations of honeycomb media, $-\nabla\cdot A\nabla\ \to\
-\nabla\cdot\ \left(\ A({\bf x}) +\delta B({\bf x})\ \right) \nabla$    which break either
 $\mathcal{P}-$ or $\mathcal{C}-$ invariance; thus, $[\mathcal{PC},{\mathcal{L}^B}]\ne0$.  Here, $\mathcal{P}-$ refers to parity inversion, $f({\bf x})\mapsto\mathcal{P}[f]({\bf x})=f(-{\bf x})$ and $\mathcal{C}-$ refers to complex conjugation, $f({\bf x})\mapsto\mathcal{C}[f]({\bf x})=\overline{f({\bf x})}$.
 If  $\mathcal{PC}-$ invariance is preserved, then the Dirac points persist, although their associated quasi-momenta may no longer be located on the vertices of the Brillouin Zone.

 \item {\it Topologically protected edge states (Theorems \ref{edgemode_thm} and \ref{rig-edge}):}
 %
 %
 We consider a honeycomb structure perturbed by an edge, across which the structure is adiabatically modulated by a domain wall.
 %
 \footnote{We comment briefly on the terms {\it edge} and {\it edge state}.
   An edge is frequently understood to mean an abrupt termination of bulk structure.
The terms ``edge'' for a line-defect across which there is a change in a key characteristic of the structure, and ``edge state'' are also used in the physics literature; see, for example, \cites{HR:07, RH:08,Shvets-PTI:13}.
The edge states we discuss are of the latter type.}
 %

 Starting from a honeycomb operator with Dirac points, ${\mathcal{L}^A}$, we first consider operators
 $\mathcal{L}^\delta_{\pm\infty}=-\nabla\cdot\left[A({\bf x}) \pm\delta\eta_{\infty}B({\bf x})\right]  \nabla$.
 Here, $\delta$ is small non-zero real parameter, $\eta_\infty>0$,  and $B({\bf x})$ breaks $\mathcal{PC}-$ invariance, {\it i.e.} $[\mathcal{PC},B]\ne0$. It follows from (c) that the operators $\mathcal{L}^\delta_{\pm\infty}$ have a local spectral gap, {\it i.e.} a gap for quasi-momenta varying near  Dirac points. Next, transverse to the line $\mathbb{R}{\bm{\mathfrak{v}}}_1$ in the lattice direction ${\bm{\mathfrak{v}}}_1\in\Lambda_h$,  we interpolate between operators $\mathcal{L}^\delta_{\pm\infty}$
with a {\it domain wall function} $\eta(\zeta)$:
 \begin{equation}
  \label{dw_function}
\eta(0)=0,\ \  \eta(\zeta) \rightarrow \pm \eta_{\infty}\ \text{with} \ \ \eta_{\infty}>0 \ \ \text{as} \ \ \zeta\rightarrow\pm\infty,
 \end{equation}
by introducing the operator
\begin{equation}
  \mathcal{L}_{\rm dw}^{(\delta)} \equiv -\nabla\cdot\left[A({\bf x}) +\delta\eta( \delta{\bm{\mathfrak{K}}}_2 \cdot {\bf x})B({\bf x})\right]  \nabla\ .
\label{L_dw} \end{equation}
Here, ${\bm{\mathfrak{K}}}_2\in\Lambda_h^*$ is such that ${\bm{\mathfrak{K}}}_2\cdot{\bm{\mathfrak{v}}}_1=0$.
 %
  The constitutive matrix: $A({\bf x})+\delta\eta(\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})B({\bf x})$ defines a medium which is an interpolation (small and adiabatic), transverse to the  line-defect / ``edge'' $\mathbb{R}{\bm{\mathfrak{v}}}_1$,  between perturbed (``gapped '') honeycomb structures.

  Associated with the  ${\bf x}\mapsto{\bf x}+{\bm{\mathfrak{v}}}_1$ translation invariance of
   $\mathcal{L}_{\rm dw}^{(\delta)}$  is
   a parallel quasi-momentum, denoted ${k_{\parallel}}$. An edge state for parallel quasi-momentum, ${k_{\parallel}}$,  is a non-trivial solution of the eigenvalue problem $\mathcal{L}_{\rm dw}^{(\delta)}\Psi =E\ \Psi$ such that
   $ \Psi({\bf x}+{\bm{\mathfrak{v}}}_1)=e^{i{k_{\parallel}}}\Psi({\bf x})$ and $\Psi({\bf x})\to0$ as $|{\bm{\mathfrak{K}}}_2\cdot{\bf x}|\to\infty$. The edge state eigenvalue problem is naturally posed on $L^2(\Sigma)$, where $\Sigma=\mathbb{R}^2/\mathbb{Z}{\bm{\mathfrak{v}}}_1$ is a cylinder. Edge states are propagating (plane-wave like) parallel to the edge and are localized transverse to the edge.

 Theorem \ref{edgemode_thm} presents the formal multiple scale expansion of such states. These modes bifurcate from the continuous spectrum at the Dirac point energy of the underlying, unperturbed ($\delta=0$) operator as the perturbation parameter $\delta$ is varied away from zero. The bifurcation of these modes is ``topologically protected'' in the sense that these modes persist in the presence of spatially localized (even large) perturbations of the domain-wall function, $\eta(\zeta)$.

 The rigorous formulation of Theorem \ref{edgemode_thm} is given in Theorem \ref{rig-edge}.
 Section \ref{rigor} contains a detailed sketch of the proof of Theorem \ref{rig-edge}, based on the corresponding result for a class of Schr\"odinger operators (see \cites{FLW-AnnalsPDE:16}
 and \cites{FLW-MAMS:17}). A central role is played by the \emph{spectral no-fold condition} which is motivated and discussed. The validity of this condition for low contrast honeycomb structures for zigzag edges follows directly from the analysis for the corresponding Schr\"odinger operator case \cites{FLW-AnnalsPDE:16}.
 Figure \ref{E_delta_zz} illustrates numerical computations of bifurcations of edge states.

 \begin{figure}
 \centering
 \includegraphics[width=\textwidth]{E_delta_zz_L20_extended}
 \caption{\footnotesize
 {\bf Left panels}: Bifurcation curves of topologically protected edge states (red curves) for the zigzag edge, illustrated by the $L^2_{{k_{\parallel}}}(\Sigma)-$ energy spectrum of $\mathcal{L}_{\rm dw}^{(\delta)}$ (see \eqref{L_dw}) vs. the perturbation parameter, $\delta$, for ${k_{\parallel}}={\bf K}\cdot{\bf v}_1=2\pi/3$.
 The lowest two continuous spectral bands are shown in black.
 Here, $A({\bf x})=a(x)I_{_{2\times2}}$ defines a honeycomb structured medium, {\it i.e.} $a(x)$ is real, even and $a(R^*{\bf x})=a({\bf x})$. We take $\eta(\zeta)=\tanh(\zeta)$, a domain-wall function, and $B({\bf x})$ is chosen to break $\mathcal{C}-$  or $\mathcal{P}-$  symmetries, respectively, in the top and bottom panels.
 {\bf Right panels}: Edge states corresponding to red $X$'s in the left panels.
 }
 \label{E_delta_zz}
 \end{figure}

 \item \emph{Group velocity properties of wavepackets concentrated along edges:}
 Theorem \ref{near-kpar1} shows the existence of edge states for all ${k_{\parallel}}$ near
  ${k_{\parallel}}={\bf K}\cdot{\bm{\mathfrak{v}}}_1$ and ${k_{\parallel}}= -{\bf K}\cdot{\bm{\mathfrak{v}}}_1$.
   Taking a weighted and continuous superposition of such edge states in ${k_{\parallel}}$, yields fully localized wave-packets which are concentrated along the edge. The group velocity of these wave-packets is given by $\partial_{k_{\parallel}} E({k_{\parallel}})$ evaluated at ${k_{\parallel}}=\pm{\bf K}\cdot{\bm{\mathfrak{v}}}_1$.

\item \emph{Unidirectional edge propagation in time-reversal invariant and $\mathcal{C}-$ breaking media:} In Section \ref{HR-K} we apply our results to the propagation of edge states in the two examples discussed in Section \ref{intro}: (1) magneto-optic media, and (2) bi-anisotropic media; see also Appendix \ref{maxwell}. We consider the special cases where $\mathcal{PC}-$ invariance is broken by  perturbations which are (a) $\mathcal{P}-$ invariant, but $\mathcal{C}-$ anti-symmetric, and
 (b) $\mathcal{C}-$ invariant, but $\mathcal{P}-$ anti-symmetric. In case (a),  the families of edge states constructed are uni-directional, although Maxwell's equations in this setting is  time-reversal invariant (Section \ref{Maxwell_symmetry}); see the top panel of Figure \ref{E_kpar_zz}. In case (b), the resulting edge states propagate  bi-directionally; see the bottom panel of Figure \ref{E_kpar_zz}.

  We remark that in \cites{HR:07,RH:08}, which considers magneto-optic media, and \cites{Shvets-PTI:13}, which considers bi-anisotropic media,  it was assumed that the unperturbed bulk medium is isotropic ($A({\bf x})=a({\bf x})I_{_{2\times2}}$). In this article we allow for a more general class of  anisotropic media; see our definition of {\it honeycomb structured medium} in Section \ref{honey-media}.

 \begin{figure}
 \centering
 \includegraphics[width=\textwidth]{E_kpar_zz_L20}
 \caption{\footnotesize
 $L^2_{{k_{\parallel}}}(\Sigma)-$ energy spectrum of $\mathcal{L}_{\rm dw}^{(\delta=1)}$ vs. parallel quasi-momentum, ${k_{\parallel}}$. $\mathcal{L}_{\rm dw}^{(\delta=1)}$ is the same as given in each of the respective (top and bottom) left panels in Figure \ref{E_delta_zz}.
 {\bf Top panel:} \emph{$\mathcal{C}-$ symmetry breaking}.
 Topologically protected (domain-wall induced) edge modes are denoted by the red curves. These modes are uni-directional.
 Spurious (gray) ``hard edge'' modes are a result of the finiteness effects associated with the choice of numerical domain; see Appendix \ref{numerical_schemes} for a detailed discussion of the numerics and a discussion of the spurious modes.
 {\bf Bottom panel:} \emph{$\mathcal{P}-$ symmetry breaking}.
 The (red) edge modes are denoted by the red curves are bi-directional.}
 \label{E_kpar_zz}
 \end{figure}

 \end{enumerate}

\subsection{Connections to previous rigorous analytical work}
The existence of Dirac points for Schr\"odinger operators with generic honeycomb lattice potentials was proved in \cites{FW:12,FLW-MAMS:17}. Results for small amplitude potentials were obtained in  \cites{Colin-de-Verdiere:91,Grushin:09}.
A group representation perspective on the existence and persistence of Dirac points is developed in
  \cites{berkolaiko-comech:15}. The case where the potential is a superposition of delta function potentials centered on sites of the honeycomb structure is treated in \cite{Lee:16}. The low-lying dispersion surfaces of honeycomb Schr\"odinger operators in the strong binding regime, where the potential is the superposition of a general class of atomic potential wells, and its relation to the  tight-binding limit, was studied in
  \cites{FLW-CPAM:17}. Tight binding models have been studied extensively in the fundamental and applied physics and mathematics communities;
   see, for example, \cite{Ashcroft-Mermin:76,Dimassi-Sjoestrand:99,ACZ:12,Makwana-Craster:14}.

   A bifurcation theory of topologically protected bound states arising from domain wall perturbations of one-dimensional periodic Schr\"odinger operators with Dirac points (linear band crossings) is developed in
  \cites{FLW-PNAS:14,FLW-MAMS:17}. These results are  applied to a photonic setting in \cites{Thorp-etal:15} and \cites{poo2016observation,poo2017global}. The bifurcation of topologically protected edge states arising from domain wall perturbations of two-dimensional honeycomb Schr\"odinger operators with Dirac points (conical band crossings) was studied in \cites{FLW-AnnalsPDE:16}; see also
    \cites{FLW-2DM:15}.

    For an extensive discussion of Dirac points and edge states for nanotube structures in the context of  quantum graphs, see \cites{Kuchment-Post:07} and \cites{Do-Kuchment:13}.
For results on spectral gaps for elliptic problems modeling high contrast media,  see
 \cites{Figotin-Kuchment:96a,Figotin-Kuchment:96b,Hempel-Lienau:00,Hempel-Post:02,Lipton-Viator_1,Lipton-Viator_2,FLW-CPAM:17}.
  For a general review of mathematical problems arising in photonics see, for example, \cite{kuchment-01}.

\subsection{Outline}  The paper is structured as follows.
 In Section \ref{math-problem}, we review the relevant spectral (Floquet-Bloch) theory and introduce the notions of
 $\mathcal{C}$, $\mathcal{P}$ and $\mathcal{R}$ invariance, associated with complex conjugation, parity and rotational symmetries. Section \ref{HSM} introduces and characterizes honeycomb structured media via Fourier analysis.
 In Section \ref{dirac-pts} we discuss the notion of {\it Dirac point} in the band structure of the operator ${\mathcal{L}^A}$ and prove
 results for the cases of low-contrast and generic honeycomb structures.  Section \ref{dirac_persistence} discusses the stability of Dirac points against small perturbations which preserve $\mathcal{P}\circ\mathcal{C}$ symmetry and their instability when such symmetry is violated.  In  Section \ref{edge_states} we discuss the existence of topologically protected edge states. We first present the detailed formal multi-scale construction of such states and then give a detailed guide to adapting the proof for the case of Schr\"odinger operators \cites{FLW-AnnalsPDE:16} to the case of divergence form operators  for honeycomb structured media. Note that results of this paper apply to large classes of complex-valued and anisotropic matrix-valued constitutive laws; see the examples arising in electromagnetics, discussed in Appendix \ref{maxwell}.
Edge states are constructed for a range of parallel quasi-momenta ${k_{\parallel}}$. In Section  \ref{edge_state_kpar} we discuss
the edge state dispersion curves ${k_{\parallel}}\mapsto E({k_{\parallel}})$, which determine the propagation properties of edge-localized wave-packets; see, in particular, the discussion of group velocity in Subsection \ref{direction}.
Appendix \ref{maxwell} discusses Maxwell's equations for bi-anisotropic media and the application of our general results to edge states and unidirectional edge-propagation in (1) magneto-optic and (2) bi-anisotropic media, both of which
are time-reversal invariant.
Appendix \ref{numerical_schemes} provides a brief discussion of the numerical methods used in our computer simulations. Appendix \ref{fig_potentials} provides the specific potentials used in these simulations.

\subsection{Notation and conventions}
\begin{enumerate}
\item $\Lambda_h\subset\mathbb{R}^2$ denotes the equilateral triangular lattice.\  $\Lambda_h^*\subset (\mathbb{R}^2)^*=\mathbb{R}^2_{\bf k}$
denotes the dual  lattice. ${\bf v}_j, ~j=1,2$ are the basis vectors of  $\Lambda_h$.
${\bf k}_l, ~l=1,2$ are the dual basis vectors of $\Lambda_h^*$, chosen to satisfy ${\bf k}_l\cdot{\bf v}_j=2\pi\delta_{lj}$.
\item ${\bm{\mathfrak{v}}}_1=a_1{\bf v}_1+a_2{\bf v}_2\in\Lambda_h$,\ $a_1, a_2$ co-prime integers. The ${\bm{\mathfrak{v}}}_1$- edge is $\mathbb{R}{\bm{\mathfrak{v}}}_1$.
 ${\bm{\mathfrak{v}}}_j,\ j=1,2$, is an alternate basis for $\Lambda_h$ with corresponding dual basis, ${\bm{\mathfrak{K}}}_\ell, \ell=1,2$, satisfying ${\bm{\mathfrak{K}}}_\ell\cdot{\bm{\mathfrak{v}}}_j=2\pi\delta_{\ell j}$.
  \item For ${\bf m}=(m_1,m_2)\in\mathbb{Z}^2$, ${\bf m}\vec{\bf k}=m_1{\bf k}_1+m_2{\bf k}_2$.
 \item $\mathcal{B}$ denotes the Brillouin Zone, associated with $\Lambda_h$, shown in the right panel of Figure \ref{lattice-and-itsdual}.
\item $z \in \mathbb C \Rightarrow \overline{z}$ denotes the complex conjugate of $z$.
 \item $J=\begin{pmatrix}0&-1\\1&0\end{pmatrix}$ is the counterclockwise $90^\circ$ rotation matrix.
 \item ${\bf k}=(k^{(1)},k^{(2)})^T\in\mathbb R^2$ is represented in $\mathbb{C}$ by $\mathfrak{z}=k^{(1)}+ik^{(2)}$, and ${\bf k}_\perp\equiv J{\bf k}=(-k^{(2)},k^{(1)})^T$ is represented in $\mathbb C$ by $\mathfrak{z}_\perp=i\mathfrak{z}({\bf k})=\mathfrak{z}({\bf k}_\perp)$.
  The corresponding points on the unit circle in $\mathbb{C}$ are denoted $\hat{\mathfrak{z}}({\bf k})$ and $\hat{\mathfrak{z}}({\bf k}_\perp)$.
\item $\nabla=(\partial_{x_1}, \partial_{x_2})^T$.
\item $\inner{f,g}_{D}=\int_{D}\overline{f}g$  is the $L^2(D)$ inner product. If the region of integration is not specified, it is assumed to be a choice of fundamental period cell, $\Omega$.
\item Pauli matrices:
\[ \sigma_1=\begin{pmatrix}0&1\\1&0\end{pmatrix},\qquad \sigma_2=\begin{pmatrix}0&-i\\i&0\end{pmatrix}=iJ,\qquad \sigma_3=\begin{pmatrix}1&0\\0&-1\end{pmatrix}.\]
\item Unless otherwise specified, summation convention over repeated indices is assumed.
\end{enumerate}

\begin{acknowledgement}
The authors wish to thank C.L. Fefferman, L. Lu, M. Rechtsman, D. Ketcheson, V. Quenneville-B{\'e}lair and  N. Yu
 for stimulating discussions.
This research was supported in part by NSF grants: DMS-1412560, DMS-1620418, DGE-1069420 and Simons Foundation Math + X Investigator grant \#376319 (MIW); and the NSF grant DMR-1420073 (JPL-T).
YZ acknowledges the hospitality of the Department of Applied Physics and Applied Mathematics during academic visits to Columbia University,  supported by Tsinghua University Initiative Scientific Research Program \# 20151080424 and NSFC grants \#11471185 and \#11871299.
\end{acknowledgement}

\section{Preliminaries}\label{math-problem}

In this section, we outline the relevant spectral theory
 \cites{Eastham:74, RS4,kuchment2012floquet,kuchment2015overview} and introduce terminology and notation for discussing the symmetry properties of the (unperturbed) bulk operator ${\mathcal{L}^A}$.

\subsection{Fourier analysis}\label{fourier-analysis}

Let $\{{\bf v}_1, {\bf v}_2\}$ be a linearly independent set in $\mathbb{R}^2$. The lattice generated by $\{{\bf v}_1, {\bf v}_2\}$ is the subset of $\mathbb{R}^2$:
\[\Lambda=\{{\bf m}\vec{\bf v}=m_1{\bf v}_1+m_2{\bf v}_2: {\bf m}=(m_1, m_2)\in \mathbb{Z}^2\}=\mathbb{Z}{\bf v}_1\oplus\mathbb{Z}{\bf v}_2 .\]
A choice of fundamental cell is the parallelogram:
\[\Omega=\{\theta_1{\bf v}_1+\theta_2{\bf v}_2: 0\le\theta_j\le1,j=1,2\}.\]
The area of $\Omega$ is denoted $|\Omega|$.

The dual lattice, $\Lambda^*$, is defined to be
\[\Lambda^*=\{{\bf m}\vec{\bf k}=m_1{\bf k}_1+m_2{\bf k}_2: {\bf m}=(m_1, m_2)\in \mathbb{Z}\}=\mathbb{Z}{\bf k}_1\oplus\mathbb{Z}{\bf k}_2,\]
where ${\bf k}_1$ and ${\bf k}_2$ are dual lattice vectors, satisfying the reciprocal relations
\[ {\bf k}_i\cdot{\bf v}_j=2\pi \delta_{ij}.\]
We also introduce the (affine) lattice $\Lambda^*_{\bf k}$:
\begin{equation}
 \label{dual_lattice_k}
 \Lambda^*_{{\bf k}}\equiv {\bf k}+\Lambda^{*}=\left\{{\bf k}+{\bf m}\vec{\bf k}: {\bf m}\in \mathbb Z^2\right\}.
\end{equation}

Let $L^2(\mathbb{R}^2/\Lambda)$ denote the space of $L^2_{loc}$ functions which are $\Lambda-$periodic, i.e.
$f({\bf x}+{\bf v})=f({\bf x})$, for almost all $ {\bf x}\in \mathbb{R}^2$ and all $ {\bf v}\in \Lambda$.
For ${\bf k}\in\mathbb{R}^2$, we denote by $L^2_{{\bf k}}$ the space of $L^2_{loc}$ functions such that
 $e^{-i{\bf k}\cdot{\bf x}}f({\bf x})\in L^2(\mathbb{R}^2/\Lambda)$, {\it i.e.}
 $f({\bf x}+{\bf v})=e^{i{\bf k}\cdot {\bf v}}f({\bf x}), \quad\text{for} \quad {\bf x}\in \mathbb{R}^2,\ {\bf v}\in \Lambda$.
Note that if ${\bf k}=0$, $L^2_{{\bf k}}$ reduces to $L^2(\mathbb{R}^2/\Lambda)$.

Let $f$ and $g$ in $L^2_{{\bf k}}$. Then, $\bar{f} g $ is locally integrable and $\Lambda-$periodic, and their inner product is naturally defined as
\[\inner{f,g}=\intomega{\overline{f({\bf x})}g({\bf x})}.\]
In a standard way, one can also introduce the Sobolev space $H^s_{{\bf k}}$.

Any $f\in L^2_{{\bf k}}$ has the Fourier expansion
\begin{align}
  f({\bf x}) &=\sum_{{\bf m} \in \mathbb{Z}^2}f_{\bf m} e^{i ({\bf k}+{\bf m}\vec{\bf k})\cdot\, {\bf x}},\ \ \textrm{where}\ \
f_{\bf m} = \frac{1}{|\Omega|}\int_\Omega e^{-i({\bf k}+{\bf m}\vec{\bf k})\cdot \,{\bf y}}f({\bf y}) d{\bf y}. \label{fourier}
\end{align}

\subsection{Floquet-Bloch theory of ${\mathcal{L}^A}=-\nabla\cdot A\nabla$}\label{fb_theory}

We summarize definitions and results from the theory of self-adjoint, elliptic and periodic divergence form
operators; see, for example,  \cite{Kuchment_Levendorskii_01,kuchment2012floquet,kuchment2015overview}.

\begin{assumption}\label{A_Assumption}
The $2\times2$ complex-valued matrix function $A({\bf x})$ satisfies
\begin{enumerate}[label=(\subscript{A}{\arabic*})]
\item $A({\bf x})$ is a smooth  and Hermitian, {\it i.e.}, $A({\bf x})^{\dagger}\equiv\overline{A({\bf x})^T}=A({\bf x})$ for all ${\bf x}$.
 \item The mapping ${\bf x}\mapsto A({\bf x})$ is $\Lambda-$ periodic.
\item $A({\bf x})$ is uniformly elliptic, {\it i.e.} there exist constants $c_\pm>0$, such that for all ${\bf x}\in\mathbb{R}^2$ and all $\xi\in\mathbb{C}^2$:
   $c_-|\xi|^2\ \le\ \left\langle \xi, A({\bf x})\xi\right\rangle_{\mathbb{C}^2}\ \le c_+|\xi|^2$.
 \end{enumerate}
\end{assumption}

For each fixed ${\bf k}\in \mathbb{R}^2$, consider the $L^2_{\bf k}-$ Floquet-Bloch eigenvalue problem:
\begin{equation}\label{Eq_Eigen}
\begin{split}
{\mathcal{L}^A}\Phi({\bf x})&= E\Phi({\bf x}), \quad {\bf x} \in \mathbb{R}^2 , \\
\Phi({\bf x}+{\bf v})&=e^{i{\bf k}\cdot {\bf v}}\Phi({\bf x}), \quad  {\bf x}\in\mathbb{R}^2,\quad {\bf v}\in \Lambda ,
\end{split}
\end{equation}
where ${\mathcal{L}^A}=-\nabla\cdot A\nabla$ was defined in \eqref{L_def}.
A solution of \eqref{Eq_Eigen} is called a Floquet-Bloch state with quasi-momentum ${\bf k}$.

Since the ${\bf k}$-pseudo-periodic boundary condition in \eqref{Eq_Eigen} is invariant under translation ${\bf k}\rightarrow {\bf k}+{{\bf k}}^\prime$ for any ${{\bf k}}^\prime\in \Lambda^*$, it suffices to consider ${\bf k}$ varying over a fundamental cell. A common choice is the {\it Brillouin Zone} $\mathcal{B}$, consisting of points ${\bf k}\in \mathbb{R}^2$ which are closer to the origin than to any other  points in $\Lambda^*$; see Figure  \ref{lattice-and-itsdual}.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{lattice-and-itsdual}
\caption{\footnotesize
{\bf Left panel}: ${\bf A}=(0,0)$, ${\bf B}=(\frac{1}{\sqrt3},0)$.
Honeycomb structure, ${\bf H}$,  is the union of two sub-lattices  $\Lambda_{\bf A}={\bf A}+\Lambda_h$ (blue)
and $\Lambda_{\bf B}={\bf B}+\Lambda_h$ (red); several hexagons shown. The lattice vectors  $\{{\bf v}_1,{\bf v}_2\}$ generate $\Lambda_h$.
{\bf Right panel}:
Brillouin zone, $\mathcal{B}_h$, and dual basis $\{{\bf k}_1,{\bf k}_2\}$. ${\bf K}$ and ${\bf K}'$ are labeled.}
\label{lattice-and-itsdual}
\end{figure}

An alternative formulation of the eigenvalue problem \eqref{Eq_Eigen}  is obtained by setting $\Phi({\bf x})=e^{i{\bf k}\cdot {\bf x}}\phi({\bf x})$ for ${\bf k}\in\mathcal{B}$. Let
\begin{equation}
\label{Eq_Eigen_k_op}
{\mathcal{L}^A}({\bf k})= e^{-i{\bf k}\cdot{\bf x}}\ \mathcal{L}^A\ e^{i{\bf k}\cdot{\bf x}}\ =\ -(\nabla+i{\bf k})\cdot A({\bf x})(\nabla+i{\bf k}).
\end{equation}
Then, \eqref{Eq_Eigen} is equivalent to  $\phi$ satisfying the  periodic elliptic boundary value problem:
\begin{equation}\label{Eq_Eigen_k}
\begin{split}
{\mathcal{L}^A}({\bf k})\phi({\bf x}) &= E({\bf k})\phi({\bf x}), \quad {\bf x} \in \mathbb{R}^2 ,\\
\phi({\bf x}+{\bf v}) &=\phi({\bf x}), \quad {\bf x}\in\mathbb{R}^2,\ \ {\bf v}\in \Lambda .
\end{split}
\end{equation}

For each fixed ${\bf k}$, the eigenvalue problem \eqref{Eq_Eigen_k} (equivalently \eqref{Eq_Eigen}) has a discrete spectrum: $ E_1({\bf k})\le E_2({\bf k})\le E_3({\bf k})\le\cdots\le E_b({\bf k})\le \cdots$
with eigenpairs $( \phi_b({\bf x};{\bf k})$, $ E_b({\bf k}))$, $b\ge1$. The set $\left\{\phi_b({\bf x};{\bf k}):b\ge1\right\}$ can be taken to be a complete and orthonormal set in $L^2(\mathbb{R}^2/\Lambda)$.

The mappings ${\bf k} \mapsto E_b({\bf k})$ are called band dispersion functions. They are Lipschitz-continuous, see, for example, \cites{Avron-Simon:78,Kuchment_Levendorskii_01,kuchment2012floquet,kuchment2015overview} and Appendix A of \cites{FW:14}. As ${\bf k}$ varies over $\mathcal{B}$, each function $ E_b({\bf k})$ sweeps out a closed real interval. The union over $b\ge 1$ of these closed intervals is precisely the $L^2(\mathbb{R}^2)$-spectrum of  $\mathcal{L}^A$.
Furthermore,  the set $\left\{\Phi_b({\bf x};{\bf k})\right\}_{b\ge 1, {\bf k} \in \mathcal{B}}$ is complete in $L^2(\mathbb{R}^2)$:
\begin{equation*}
  f({\bf x})=\displaystyle{\sum_{b\ge1}\int_\mathcal{B}\inner{\Phi_b(\cdot;{\bf k}), f(\cdot)}_{L^2(\mathbb R^2)}\Phi_b({\bf x};{\bf k})\,d{\bf k}},
\end{equation*}
where the sum converges in the $L^2$ norm.

\subsection{$\mathcal{C}$, $\mathcal{P}$ and $\mathcal{R}$ invariance}\label{CPR-inv}

In this section we introduce terminology for discussing the symmetry properties of $A({\bf x})$ and ${\mathcal{L}^A}$. Let $g({\bf x})$ denote a function defined on $\mathbb{R}^2$.

{\bf $\mathcal{C}-$ invariance:} The function $g({\bf x})$ is $\mathcal{C}-$ invariant if
\begin{equation}
\left(\mathcal{C}g\right)({\bf x})\ \equiv\ \overline{g({\bf x})}\ =\ g({\bf x})
\label{Cg_eq_g}
\end{equation}
for all ${\bf x}\in\mathbb{R}^2$.

{\bf $\mathcal{P}-$ invariance:} The function $g({\bf x})$ is $\mathcal{P}-$ (or parity-inversion) invariant with respect to ${\bf x}_0\in\mathbb{R}^2$  if
\begin{equation*}
\left(\mathcal{P}g\right)({\bf x})\ \equiv\ g(2{\bf x}_0-{\bf x})\ =\ g({\bf x}).
\end{equation*}
for all ${\bf x}\in\mathbb{R}^2$.

{\bf $\mathcal{R}-$ invariance:} The function $g({\bf x})$ is $\mathcal{R}-$ (or $120^\circ$- rotationally) invariant with respect to ${\bf x}_0\in\mathbb{R}^2$  if
\begin{equation*}
\left(\mathcal{R}g\right)({\bf x})\ \equiv\ g({\bf x}_0+R^*({\bf x}-{\bf x}_0))\ =\ g({\bf x}).
\end{equation*}
for all ${\bf x}\in\mathbb{R}^2$. Here, ${\bf x}\mapsto R{\bf x}$ is the mapping on $\mathbb{R}^2$ which rotates a vector clockwise by
$120^\circ$ ($2\pi/3$) about  ${\bf x}={\bf 0}$:
\begin{equation}\label{Def_R}
R=\begin{pmatrix}-\frac{1}{2} & \frac{\sqrt{3}}{2}\\[1.5 ex ]-\frac{\sqrt{3}}{2} & -\frac{1}{2}\end{pmatrix}.
\end{equation}
For later use, we record the eigenvalues and a choice of normalized eigenvectors of $R$:
\begin{equation}
\label{R_evals_evecs}
R\zeta\ =\ \tau\ \zeta, \quad  R\ \bar{\zeta}\ =\ \bar{\tau}\ \bar{\zeta}, \quad \text{where} \quad \tau=e^{i2\pi/3}\ \textrm{and}\
\zeta=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\i\end{pmatrix} .
\end{equation}

{\it We say the operator  $L$ with domain $D(L)\subset D(\mathcal{S})$, is $\mathcal{S}-$ invariant if
the commutator } $ [\mathcal{S},L]\equiv \mathcal{S}L-L\mathcal{S}$ vanishes on $D(L)$.

\section{Honeycomb structured media}\label{HSM}

\subsection{The equilateral triangular lattice, $\Lambda_h$}\label{triangle}

In this section, we specialize the results of the previous section to the triangular lattice, $\Lambda_h$.
 Introduce the lattice
\[ \Lambda_h=\mathbb{Z}{\bf v}_1 \oplus \mathbb{Z}{\bf v}_2,\]
with lattice basis vectors
 \[
 {\bf v}_1=\begin{pmatrix} \frac{\sqrt{3}}{2} \\  \\ \frac{1}{2}\end{pmatrix},\quad
 {\bf v}_2=\begin{pmatrix}\frac{\sqrt{3}}{2} \\ \\ -\frac{1}{2} \end{pmatrix}.\]
The dual lattice
\[\Lambda_h^* = \mathbb{Z}{\bf k}_1\oplus \mathbb{Z}{\bf k}_2,\] is spanned by the dual basis vectors
\[ {\bf k}_1=\frac{4\pi}{\sqrt{3}}\begin{pmatrix}\frac{1}{2} \\ \\ \frac{\sqrt{3}}{2}\end{pmatrix},\quad
 {\bf k}_2=\frac{4\pi}{\sqrt{3}}\begin{pmatrix}\frac{1}{2}\\ \\ -\frac{\sqrt{3}}{2}\end{pmatrix} ;\]
see Figure  \ref{lattice-and-itsdual}.
Note  ${\bf k}_l\cdot{\bf v}_m=2\pi \delta_{lm}$, $l,m=1,2$.

The following proposition characterizes the high-symmetry points of the Brillouin Zone, $\mathcal{B}$.
\begin{proposition}\label{special_k_points}
The translated dual lattice $\Lambda^*_{h,{\bf k}}= {\bf k} + \Lambda_h^*$ is invariant under $120^\circ$ rotation,  {\it i.e.}
 $R:\Lambda^*_{h,{\bf k}}\to \Lambda^*_{h,{\bf k}}$ is one to one and onto, if and only if ${\bf k}\in\{\bb\Gamma,{\bf K},{\bf K}^\prime\}$, where
\begin{equation}
\label{defn_K_Kp}
\bb\Gamma \equiv 0,\quad {\bf K}\equiv\frac{1}{3}({\bf k}_1-{\bf k}_2),\quad {\bf K}' \equiv-{\bf K} .
\end{equation}
\end{proposition}

The set of six vertices of $\mathcal{B}_h$ is invariant under $R$ and decomposes into the two subsets:
the ${\bf K}\,\, \text{type-points:} \ {\bf K}, \ R{\bf K} ={\bf K}+{\bf k}_2, \ R^2{\bf K} = {\bf K} - {\bf k}_1$ and
the ${\bf K}'\, \,\text{type-points:} \ {\bf K}', \ R{\bf K}' ={\bf K}-{\bf k}_2,\ R^2{\bf K}' = {\bf K}' + {\bf k}_1$.
We note for future reference, the action of $R$ on ${\bf k}_1$ and ${\bf k}_2$:
\begin{equation}
 \label{R_bk_rel}
 R{\bf k}_1={\bf k}_2, \ \ \text{and} \ \ R{\bf k}_2=-{\bf k}_1-{\bf k}_2.
\end{equation}

\subsection{Fourier analysis in $L^2_{{\bf K}_\star}$}

We focus on the vertices of $\mathcal{B}_h$. Consider such a vertex,  ${\bf K}_\star$. If $f$ is ${\bf K}_\star-$ pseudo-periodic, it is easy to check that $\mathcal{R}[f]$ is also ${\bf K}_\star-$ pseudo-periodic and that $\mathcal{R}$ is an isomorphism of $L^2_{{\bf K}_\star}$.
Furthermore,  $\mathcal R$ is unitary with eigenvalues $1$, $\tau$ and $\bar{\tau}$, where $\tau=\exp(2\pi i/3)$. It is therefore natural to split $L^2_{{\bf K}_\star}$ into the direct (orthogonal) sum:
\begin{equation*}
L^2_{{\bf K}_\star} = L^2_{{\bf K}_\star,1}\oplus L^2_{{\bf K}_\star,\tau}\oplus L^2_{{\bf K}_\star,\overline\tau}\ ,
\end{equation*} where  the summand subspaces are given by:
\begin{equation}
\label{L2_K_sigma}
L^2_{{\bf K}_\star,\sigma}~=~\left\{f\in L^2_{{\bf K}_\star}:\mathcal{R}f=\sigma f\right\}, \quad \sigma=1,\tau,\bar\tau.
\end{equation}

\subsection{$\mathcal{PC}-$ invariance and $\mathcal{O}-$ invariance}\label{honeycomb-structures}

Let  $A({\bf x})$ be a smooth $\Lambda_h-$ periodic and Hermitian matrix function, for which ${\mathcal{L}^A}=-\nabla\cdot A\nabla$ is elliptic in the sense of Assumption  \ref{A_Assumption} of Section \ref{fb_theory}. In this section we derive further conditions on $A({\bf x})$ which ensure that ${\mathcal{L}^A}$ is $\mathcal{PC}-$ invariant and $\mathcal R-$ invariant.

We make use of the following general result. Let $O$ denote a real orthogonal matrix and define
\begin{equation}
\left(\mathcal{O}f\right)({\bf x})=f({O}^*{\bf x})\ .
\end{equation}
\begin{theorem}\label{invariance}
Let $A({\bf x})$ be an arbitrary smooth matrix function.
\begin{enumerate}
\item $\left[\mathcal{PC},{\mathcal{L}^A}\right]\ =\ 0$ if and only if  $A({\bf x})=\overline{A(-{\bf x})}$.
\item $\left[\mathcal{O},{\mathcal{L}^A}\right]\ =\ 0$ if and only if  $A(O^*{\bf x})=O^*A({\bf x})O$.
\item Assume $O^TO=I$ and $O\neq \pm I$. If the matrix $A({\bf x})$ is Hermitian and $\mathcal{O}-$ invariant:
\begin{equation}
\mathcal{O}[A]({\bf x})\equiv A(O^*{\bf x})=A({\bf x}),
\label{A-Rinvariant}
\end{equation}
 then,  $[{\mathcal{L}^A},\mathcal{O}]=0$ and $[{\mathcal{L}^A},\mathcal{PC}]=0$ if and only if
 \begin{align}
   A({\bf x})&=a({\bf x})I + b({\bf x}) \sigma_2,\nonumber
   \end{align}
   where $a({\bf x})$ and $b({\bf x})$ are scalar functions satisfying
\begin{align}
   a(O^*{\bf x})&=a({\bf x}),\qquad a(-{\bf x})=a({\bf x}),\nonumber\\
b(O^*{\bf x})&=b({\bf x}),\qquad  b(-{\bf x})=-b({\bf x}).\nonumber
   \end{align}
  %
\end{enumerate}
\end{theorem}

The proof of Theorem \ref{invariance} makes use of the following basic lemma:
\begin{lemma}\label{PCRnabla}
\[
(a)\ \mathcal{P}\nabla=-\nabla\mathcal{P},\qquad
(b)\ \mathcal{C}\nabla=\nabla\mathcal{C},\qquad
(c)\  \mathcal{O}\nabla=O^*\nabla\mathcal{O}.
\]
\end{lemma}
\begin{proof}[Proof of Lemma \ref{PCRnabla}]  We prove part (c). For any differentiable $f$, $\partial_{x_m}\mathcal{O}[f]({\bf x})=\partial_{x_m}f(O^*{\bf x})=\frac{\partial y_l}{\partial x_m}\partial_{y_l}f({\bf y})\Big|_{{\bf y}=O^*{\bf x}}$. Now, $y_l=(O^*{\bf x})_l=O^*_{lr}x_r$ and therefore $\frac{\partial y_l}{\partial x_m}=O_{lr}^*\delta_{rm}=O^*_{lm}$.  Therefore,
$\partial_{x_m}\mathcal{O}[f]({\bf x})= O^*_{lm}\partial_{y_l}f({\bf y})\Big|_{{\bf y}=O^*{\bf x}}=O_{ml}\partial_{y_l}f({\bf y})\Big|_{{\bf y}=O^*{\bf x}} = O_{ml}\mathcal{O}\left[\partial_{y_l}f\right]({\bf x})$. This
implies $O^*_{lm}\partial_{x_m}\mathcal{O}[f]({\bf x})= \mathcal{O}\left[\partial_{y_l}f\right]({\bf x})$ for each $l=1,2$,
which completes the proof.
\end{proof}

We now turn to the proof of Theorem \ref{invariance}.
\begin{proof}[Proof of Theorem \ref{invariance}]
 Part 1 is straightforward. To prove part 2, note by part (c) of Lemma \ref{PCRnabla} that
 $\partial_{x_j}\mathcal{O}[f]({\bf x})=\mathcal{O}\left[O_{jl}\partial_{y_l}f\right]({\bf x})$. Multiplying this relation by $a_{ij}({\bf x})$, using that
  $\mathcal{O}[a_{ij}(O{\bf y})]({\bf x})=a_{ij}({\bf x})$,  and summing over $j$ gives,
 \[ a_{ij}({\bf x})\partial_{x_j}\mathcal{O}[f]({\bf x})=\mathcal{O}\left[a_{ij}(O{\bf y})O_{jl}\partial_{y_l}f\right]({\bf x}).\]
 Furthermore, again using Lemma \ref{PCRnabla}, we have
 \begin{align}
({\mathcal{L}^A}\circ\mathcal{O})[f]({\bf x}) &=\ \partial_{x_i} a_{ij}({\bf x})\partial_{x_j}\mathcal{O}[f]({\bf x})\nonumber\\
                 &=\partial_{x_i} \mathcal{O}\left[a_{ij}(O{\bf y})O_{jl}\partial_{y_l}f\right]({\bf x})\nonumber\\
& = \mathcal{O}\left[O_{im}\partial_{y_m} a_{ij}(O{\bf y})O_{jl}\partial_{y_l}f\right]({\bf x})\nonumber\\
& =  \mathcal{O}\left[\partial_{y_m} O_{mi}^*a_{ij}(O{\bf y})O_{jl}\partial_{y_l}f\right]({\bf x})\nonumber\\
& = (\mathcal{O}\circ\mathcal{L}^{\mathcal{O}_\star A})[f]({\bf x})\ ,
\end{align}
where $[\mathcal{O}_\star A]({\bf y})\equiv O^*A(O{\bf y})O$.
Therefore, $[\mathcal{O},{\mathcal{L}^A}]=0$ if and only if $\mathcal{O}_\star A({\bf y})=A({\bf y})$. Setting ${\bf y}=O^*{\bf x}$, it follows that
 $O^*A({\bf x})O=A(O^*{\bf x})$. This completes the proof of part 2.

Finally, we turn to the proof of part 3. By hypothesis
 $O=\begin{pmatrix}\cos \theta& \sin \theta \\ -\sin \theta &\cos \theta,
 \end{pmatrix}$ where $\theta \in (0,2\pi)$ and $\theta\ne\pi$.  Form the $4-$component vector of entries of $A$: $\vec{a}=(a_{11},a_{12},a_{21},a_{22})^T$. Since $A(O^*{\bf x})=A({\bf x})$, by part 2 we have the relation $A({\bf x})= O^* A({\bf x}) O$, which is equivalent to the system $M\vec{a}=0$,  where $M$ is the $4\times4$ matrix:
 \[
 M=\sin\theta\times\begin{pmatrix}-\sin\theta & -\cos\theta & -\cos\theta & \sin\theta \\
 \cos\theta &-\sin\theta &-\sin\theta & -\cos\theta\\
 \cos\theta&-\sin\theta &-\sin\theta &-\cos\theta \\ \sin\theta&\cos \theta &\cos\theta &-\sin\theta\\
  \end{pmatrix}.
 \]

Note that the second and third rows of $M$ are identical, and  the first and fourth rows are the same up to a factor of minus one. Therefore, since $\theta\neq 0,\ \pi$, the rank of $M$ is 2.  Hence, the dimension of the nullspace of $M$ is 2. In fact, the null space of $M$ is easily seen to be spanned by the vectors: $(1,0,0,1)^T$ and $(0,-1,1,0)^T$. Thus, any element of the nullspace of $M$, may be expressed in the form $a({\bf x})(1,0,0,1)^T+ ib({\bf x})(0,-1,1,0)^T$, for some choice of scalar functions $a$ and $b$. Using that $A^\dagger=A$, we have $A({\bf x})=a({\bf x})I + b({\bf x}) \sigma_2$, where $a({\bf x}), b({\bf x})$ are real-valued. Finally, by part 1 we have $\overline{A(-{\bf x})}=A({\bf x})$, and therefore $a({\bf x})$ is even and $b({\bf x})$ is odd.
 This completes the proof of Theorem \ref{invariance}.
 \end{proof}

\subsection{Characterization of honeycomb structured media}\label{honey-media}

Recall that $R$ denotes the $120^\circ$ rotation matrix \eqref{Def_R} and the mapping $\mathcal{R}:g({\bf x})\mapsto g(R^*{\bf x})$. The bulk structures we consider in this article, modeled by the operator ${\mathcal{L}^A}=-\nabla\cdot A\nabla$ are assumed to satisfy:
\begin{enumerate}
\item ${\mathcal{L}^A}$ self-adjoint and elliptic. Hence, $A^\dagger=A$ and ellipticity condition $(A_3)$ of Assumption \ref{A_Assumption} holds.
\item $A({\bf x}+{\bf v})=A({\bf x})$ for all ${\bf x}\in\mathbb{R}^2$ and ${\bf v}\in\Lambda_h$.
\item  $[\mathcal{PC}, {\mathcal{L}^A}]=0$ and hence $\overline{A(-{\bf x})}=A({\bf x})$, by Theorem \ref{invariance} part 1.
 \item $[\mathcal{R}, {\mathcal{L}^A}]=0$ and hence $A(R^*{\bf x})=R^*A({\bf x})R$,  by Theorem \ref{invariance} part 2. \\
 \end{enumerate}

 We shall also consider the more restrictive case where, in addition to (1)-(4), ${\mathcal{L}^A}=-\nabla\cdot A\nabla$ satisfies:
 \begin{enumerate}[5.]
 \item $A(R^*{\bf x})=A({\bf x})$ for all ${\bf x}\in\mathbb{R}^2$. Hence, by Theorem \ref{invariance} part 3
there exist real-valued functions $a({\bf x})$ (strictly positive) and $b({\bf x})$, such that
 \begin{align}
 A({\bf x})&=a({\bf x})I_{_{2\times2}}+b({\bf x})\sigma_2, \label{ARinv}
 \end{align} where
$a(R^*{\bf x})=a({\bf x})$, $a(-{\bf x})=a({\bf x})$, and
$b(R^*{\bf x})=b({\bf x})$, $b(-{\bf x})=-b({\bf x})$.
\end{enumerate}

We shall loosely refer to media which are modeled by ${\mathcal{L}^A}$, satisfying the above properties (1)-(4) or (1)-(5)
 {\it honeycomb structured media}. When we impose the stronger constraint of $\mathcal{R}-$ invariance of the matrix $A({\bf x})$ itself, {\it i.e.} $A(R^*{\bf x})=A({\bf x})$, (5), we shall point this out explicitly.
 %

 %

We shall next derive the form of the  Fourier series of a general  matrix function giving rise to a honeycomb structured medium. The analysis follows that of the Schr\"odinger operator case, presented in \cites{FW:12}. We begin by noting that if ${\bf k}=m_1{\bf k}_1+m_2{\bf k}_2\in\Lambda_h^*$, then $R{\bf k}=-m_2{\bf k}_1+(m_1-m_2){\bf k}_2$; see \eqref{R_bk_rel}.
Therefore, the mapping $R:\Lambda_h^*\to\Lambda_h^*$ induces a  mapping $\tilde{R}:\ \mathbb Z^2\rightarrow \mathbb Z^2$,  $\tilde{R}: (m_1, m_2) \mapsto (-m_2, m_1-m_2)$.  We note $\tilde{R}^2(m_1, m_2) = (m_2-m_1, m_1)$, $\tilde{R}^3(m_1, m_2) = (m_1, m_2)$, and that $\tilde{R}^{-1}=\tilde{R}^2$.

We say that $\mathbf{m}$ and $\mathbf{n}$ are in the same equivalence class if they lie on the same $\tilde{R}$ orbit.  Introduce the subset $\mathcal{S}\subset \mathbb{Z}^2\setminus\{{\bm 0}\}$,  consisting  of exactly one representative from each equivalence class. With the exception of  $\mathbf{m}=0$, left fixed by $\tilde{R}$, each ${\bf m}\in\mathbb{Z}^2$ lies on a $\tilde{R}-$ orbit of length
exactly three.

If $A({\bf x})$ is a smooth $\Lambda_h-$periodic matrix function, then it can be represented as a convergent Fourier series  $A({\bf x})=\sum_{{\bf m}\in\mathbb{Z}^2}A_{\bf m} e^{i{\bf m}\vec{{\bf k}}\cdot{\bf x}}$, with matrix-valued Fourier coefficients:
\begin{equation}
\label{fourier_coef}
A_{\mathbf{m}}\equiv\frac{1}{|\Omega_h|}\int_{\Omega_h} e^{-i{\bf m}\vec{\bf k}\cdot \,{\bf y}}A({\bf y}) d{\bf y} .
\end{equation}

\begin{proposition}\label{prop_fourier_coef}
\begin{enumerate}
\item Let $A(R^*{\bf x}) = R^*A({\bf x})R$. Then,
\begin{align}
A_{\tilde{R}^{-1}\mathbf{m}}=R A_\mathbf{m}R^*,\quad \mathbf{m}\in \mathbb Z^2,
\label{ARx-fr}
\end{align}
and therefore, $A_{\tilde{R}{\bf m}}=R^*A_{\bf m} R$ and $A_{\tilde{R}^2{\bf m}}=RA_{\bf m} R^*$, $\mathbf{m}\in \mathbb Z^2$.
\item Assume  $A(R^*{\bf x}) = A({\bf x})$, and therefore by part 3 of Theorem \ref{invariance},
 $A(R^*{\bf x}) = a({\bf x})I_{_{2\times2}}+b({\bf x})\sigma_2$, where $a$ is real-valued and even, and $b$ is real-valued and odd.
 Then,
$
A_{\bf m}=a_{\bf m} I_{_{2\times2}}+b_{\bf m}\sigma_2,
$ where the Fourier coefficients of $a$ and $b$ satisfy:
\begin{equation}
a_{\bf m}=a_{\tilde{R}{\bf m}}=a_{\tilde{R}^2{\bf m}}, \quad b_{\bf m}=b_{\tilde{R}{\bf m}}=b_{\tilde{R}^2{\bf m}}\\
\end{equation}
Furthermore, since $a(-{\bf x})=a({\bf x}),\ b(-{\bf x})=-b({\bf x})$, we have $\overline{a_{\bf m}}=a_{\bf m}=a_{-{\bf m}},\ \overline{b_{\bf m}}=-b_{\bf m}=b_{-{\bf m}}$.
\item Let $\mathcal{PC}[A]({\bf x}) = A({\bf x})$. Then,
\begin{align}
 A_\mathbf{m}=\overline{A_\mathbf{m}},\quad  \mathbf{m}\in \mathbb Z^2.
\label{PCA-fr}\end{align}
\item Let $A^\dagger=A$. Then,
\begin{align}
A_\mathbf{m}=(A_{-\mathbf{m}})^\dagger,\ \  \mathbf{m}\in \mathbb Z^2.\label{Asa-fr}
\end{align}
\end{enumerate}
\end{proposition}

\begin{proof}[Proof of Proposition \ref{prop_fourier_coef}]
We first prove \eqref{ARx-fr}.
\begin{equation}
\begin{split}
[\mathcal{R}A]_{\bf m}&=\frac{1}{|\Omega_h|}\int_{\Omega_h} e^{-i{\bf m}\vec{\bf k}\cdot \,{\bf y}}A(R^*{\bf y})d{\bf y}\ =\frac{1}{|\Omega_h|}\int_{\Omega_h} e^{-iR^*({\bf m}\vec{\bf k})\cdot \,\bb z}A(\bb z) d\bb z\\
&=\frac{1}{|\Omega_h|}\int_{\Omega_h} e^{-i(\tilde{R}^{-1}\bb m)\bb k\cdot \,\bb z}A(\bb z) d\bb z\ =\ A_{\tilde{R}^{-1}\mathbf{m}}
\end{split}
\end{equation}
Therefore, $\mathcal{R}[A]({\bf x})\equiv A(R^*{\bf x})=R^*A({\bf x})R$ $\iff$ $A_{\tilde{R}^{-1}\mathbf{m}}=R A_{\mathbf{m}}R^*$.

To prove \eqref{PCA-fr}, note that
\begin{equation}
\begin{split}
(\mathcal{PC})[A]_{{\bf m}}&=\frac{1}{|\Omega_h|}\int_{\Omega_h} e^{-i{\bf m}\cdot\vec{\bf k} \,{\bf y}}\overline{A(-{\bf y})}d{\bf y} =\frac{1}{|\Omega_h|}\int_{\Omega_h} e^{i{\bf m}\vec{\bf k}\cdot \,\bb z}\overline{A(\bb z)} d\bb z=\overline{A_{{\bf m}}}\\
\end{split}
\end{equation}
Therefore,
$(\mathcal{PC})[A]({\bf x})\equiv \overline{A(-{\bf x})}=A({\bf x})$ $\iff$ $A_\mathbf{m}$ is real.

To prove \eqref{Asa-fr} note that $(A_{\bf m})^{\dagger}=(A^\dagger)_{-{\bf m}}$. Hence, if $A$ is Hermitian then $A_{\bf m}=(A_{-{\bf m}})^\dagger$.
\end{proof}

Because $A({\bf x})$ is $\Lambda_h-$ periodic, we can express it in the Fourier series
\begin{equation}
\begin{split}
A({\bf x})&= A_{\bf 0}+\sum_{\mathbf{m}\ne{\bf 0}}A_\mathbf{m}e^{i\mathbf{m}\vec{{\bf k}}\cdot {\bf x}}\\
&= A_{\bf 0}+\sum_{\mathbf{m}\in\tilde{\mathcal{S}}\setminus\{\bf 0\}}
\left(\ A_\mathbf{m}e^{i\mathbf{m}\vec{{\bf k}}\cdot {\bf x}}+A_{\tilde{R}\mathbf{m}}e^{i(\tilde{R}\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}
+A_{\tilde{R}^2\mathbf{m}}e^{i(\tilde{R}^2\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}\ \right).
\end{split}
\end{equation}

Symmetries such as $\mathcal{PC}-$ and $\mathcal{R}-$ invariance imply further constraints on the Fourier series:
\begin{corollary}[Characterization of honeycomb structured media]\label{A-expand}
Let $A({\bf x})$ be a Hermitian and $\Lambda_h-$ periodic matrix function. Then, the following holds:
\begin{enumerate}
\item
Assume  $\mathcal{PC}[A]=A$ (hence $[\mathcal{PC},{\mathcal{L}^A}]=0$) and $A(R^*{\bf x})=R^*A({\bf x})R$ (hence $[\mathcal{R},{\mathcal{L}^A}]=0$).  Then,
\begin{equation}
A({\bf x})=a_{_{\bb 0}}I_{_{2\times2}}+\sum_{\mathbf{m}\in\tilde{\mathcal{S}}\setminus\{\bb 0\}}A_\mathbf{m}\ e^{i\mathbf{m}\vec{{\bf k}}\cdot {\bf x}}+R^* A_\mathbf{m}R\ e^{i(\tilde{R}\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}+R A_\mathbf{m}R^*\ e^{i(\tilde{R}^2\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}.
\label{Afourier}
\end{equation}
where $a_{_{\bb 0}}$ is real and positive; $A_\mathbf{m}$ is real and satisfies $A_{-\mathbf{m}}=A_{\mathbf{m}}^T$.
\item Assume further that $A(R^*{\bf x})=A({\bf x})$. Then,
\begin{align}
A({\bf x})&=a_{_{\bb 0}}I_{_{2\times2}}+
\sum_{\mathbf{m}\in\tilde{\mathcal{S}}\setminus\{\bb 0\}}\left(a_{\bf m} I_{_{2\times2}}+b_{\bf m}\sigma_2\right)
\left(e^{i\mathbf{m}\vec{{\bf k}}\cdot {\bf x}}+e^{i(\tilde{R}\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}+e^{i(\tilde{R}^2\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}\right)\nonumber\\
&= a_{_{\bb 0}}I_{_{2\times2}}+\sum_{\mathbf{m}\in\tilde{\mathcal{S}}\setminus\{\bb 0\}}
\begin{pmatrix}
a_{\bf m} & -\beta_{\bf m}\\ \beta_{\bf m} & a_{\bf m}
\end{pmatrix}\
\left(e^{i\mathbf{m}\vec{{\bf k}}\cdot {\bf x}}+e^{i(\tilde{R}\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}+e^{i(\tilde{R}^2\mathbf{m})\vec{{\bf k}}\cdot {\bf x}}\right),
\label{AAfourier}\end{align}
where $a_{\bf m}$ and $\beta_{\bf m}=ib_{\bf m}$ (${\bf m}\in\mathcal{S}\setminus\{\bb 0\}$) are real.
\item Under the assumptions of (2), $A({\bf x})$ may also be rewritten as
\begin{align}
A({\bf x})&=a_{_{\bb 0}}I_{_{2\times2}}\ +\
I_{_{2\times2}}\cdot \sum_{\mathbf{m}\in\tilde{\mathcal{S}}\setminus\{\bb 0\}}\ a_{\bf m}\
\left(\ \cos({\bf m}\vec{\bf k}\cdot{\bf x})\ +\  \cos((\tilde{R}{\bf m})\vec{\bf k}\cdot{\bf x})\ +\  \cos((\tilde{R}^2{\bf m})\vec{\bf k}\cdot{\bf x})\ \right)\nonumber\\
&\ + \sigma_2\cdot \sum_{\mathbf{m}\in\tilde{\mathcal{S}}\setminus\{\bb 0\}}\ \beta_{\bf m}
\left(\ \sin({\bf m}\vec{\bf k}\cdot{\bf x})\ +\  \sin((\tilde{R}{\bf m})\vec{\bf k}\cdot{\bf x})\ +\  \sin((\tilde{R}^2{\bf m})\vec{\bf k}\cdot{\bf x})\ \right) ,
\end{align}
where $a_{\bf m}$ and $\beta_{\bf m}$ are real.
\end{enumerate}
\end{corollary}

\begin{proof}[Proof of Corollary \ref{A-expand}] From Proposition \ref{prop_fourier_coef}, $A_{\bb 0}$ should be a real matrix satisfying $A_{\bb 0}=R A_{\bb 0} R^*$. Thus $A_{\bb 0}=a_{_{\bb 0}}I_{_{2\times2}}$ following from Theorem \ref{invariance}.
The rest of part 1 of Corollary \ref{A-expand} follows from Proposition \ref{prop_fourier_coef}, using that $R^2=R^*$. Part 2 is a consequence of \eqref{ARinv} and that $R^*\sigma_2R=\sigma_2=R\sigma_2R^*$. Part 3 follows from part 2 using that $A=\frac12(A+A^\dagger)$, since $A$ is Hermitian.
\end{proof}

Define the vector operator
\begin{equation}\label{A_def}
\mathscr{A}=A\ \frac{1}{i}\nabla+\frac{1}{i}\nabla\cdot A,
\end{equation}
with component operators:
\begin{equation}\label{Aj_def}
\mathscr{A}_j=\sum_{l=1}^2a_{jl}\ \frac{1}{i}\partial_{x_l}+\sum_{l=1}^2\ \frac{1}{i}\partial_{x_l}a_{lj} , \quad j=1,2.
\end{equation}

Recall, by part 3 of Lemma \ref{PCRnabla}, the relation $\mathcal{R}\nabla f=R^*\nabla \mathcal{R}f$, for all smooth scalar function $f$. We next show that if $A(R^*{\bf x})=R^*A({\bf x})R$ (which implies that $[\mathcal{R},{\mathcal{L}^A}]=0$, by part 2 of Theorem \ref{invariance})
that this relation holds with $-i\nabla$ replaced by $\mathscr{A}$.

\begin{lemma}\label{RcomA}
Let $A({\bf x})$ be a smooth matrix function and define ${\mathcal{L}^A}=-\nabla\cdot A\nabla$.
 If $\left[\mathcal{R}, {\mathcal{L}^A}\right]=0$, then  $\mathcal{R}\mathscr{A}=R^* \mathscr{A}\mathcal{R}$.
Equivalently, for $j=1,2$:
\[ \mathcal{R}\mathscr{A}_j=\sum_{l=1}^2R^*_{jl}\mathscr{A}_{l}\mathcal{R} .\]
\end{lemma}

\begin{proof}[Proof of Lemma \ref{RcomA}] For all differentiable $f$,
\begin{align}
\left(\mathcal{R}\circ\mathscr{A}_j\right)[f]({\bf x})
&=\ \mathcal{R}\left[a_{jl}({\bf y})\frac{1}{i}\partial_{y_l}f+\frac{1}{i}\partial_{y_l}a_{lj}({\bf y})f\right]({\bf x})\nonumber\\
&=\ a_{jl}(R^*{\bf x})\mathcal{R}\left[\frac{1}{i}\partial_{y_l}f\right]+\mathcal{R}\left[\frac{1}{i}\partial_{y_l}a_{lj}({\bf y})f({\bf y})\right]({\bf x}) .
\label{RcomA-step1}
\end{align}
Recall that by part 2 of Theorem \ref{invariance} and part (c) of Lemma \ref{PCRnabla} we have $A(R^*{\bf x})=R^*A({\bf x})R$ and $\mathcal{R}\nabla f=R^*\nabla\mathcal{R}f$. Therefore, continuing from \eqref{RcomA-step1}, we have
{\footnotesize{
\begin{align}
\left(\mathcal{R}\circ\mathscr{A}_j\right)[f]({\bf x}) &= \
R^*_{js}a_{sm}(R^*{\bf x})R_{ml}\cdot R_{ln}^*\frac{1}{i}\partial_{x_n}\mathcal{R}\left[f\right]+R_{ln}^*\frac{1}{i}\partial_{x_n}a_{lj}(R^*{\bf x})\mathcal{R}\left[f\right]({\bf x})\nonumber\\
&=\ R^*_{js}a_{sm}({\bf x}) \delta_{mn}\  \frac{1}{i}\partial_{x_n}\mathcal{R}\left[f\right]+R_{ln}^*\frac{1}{i}\partial_{x_n} \ R^*_{lm}a_{ms}({\bf x})R_{sj}\ \mathcal{R}\left[f\right]({\bf x})\nonumber\\
&=\ R^*_{js}a_{sm}({\bf x}) \delta_{mn}\  \frac{1}{i}\partial_{x_n}\mathcal{R}\left[f\right]+\delta_{nm}\ \frac{1}{i}\partial_{x_n} \ a_{ms}({\bf x})R_{sj}\ \mathcal{R}\left[f\right]({\bf x})\nonumber\\
&=\ R^*_{js}\left[a_{sm}({\bf x})\  \frac{1}{i}\partial_{x_m}\mathcal{R}\left[f\right]+ \frac{1}{i}\partial_{x_m} \ a_{ms}({\bf x})\ \mathcal{R}\left[f\right] \right]\nonumber\\
&=\ R^*_{js}\left[a_{sm}({\bf x})\  \frac{1}{i}\partial_{x_m}+ \frac{1}{i}\partial_{x_m} \ a_{ms}({\bf x}) \right]\ \mathcal{R}\left[f\right]
\ =\ R^*_{js}\ \mathscr{A}_s\mathcal{R}\left[f\right] .
\label{RcomA-step2}\end{align}
}}
This completes the proof of Lemma \ref{RcomA}.
\end{proof}

\section{Dirac points}\label{dirac-pts}

In this section we review the notion of {\it Dirac points}, conical points at the intersections of two dispersion surfaces.
We then provide sufficient conditions on the multiplicity of the $L^2_{{\bf K}_\star}-$ eigenvalues and certain non-degeneracies for the existence of Dirac points for the operator  ${\mathcal{L}^A}=-\nabla\cdot A\nabla$. We use this result in the Section \ref{low-dp} to prove the existence of Dirac points for ${\mathcal{L}^A}$, under the conditions that $A({\bf x})$ is smooth, Hermitian and $\Lambda_h-$ periodic, and that $[ \mathcal{PC},{\mathcal{L}^A}] = [\mathcal{R},{\mathcal{L}^A}]=0$.

\begin{definition}[Dirac points]\label{dirac_pt_defn}%
The quasi-momentum / energy pair $({\bf K}_\star,E_D)\in\mathcal{B}_h\times\mathbb{R}$ is called a {\it Dirac point} if there exists $b_\star\ge1$  and Floquet-Bloch eigenpairs mappings:
\[ {\bf k}\mapsto (\Phi_{b_\star}({\bf x};{\bf k}),E_{b_\star}({\bf k}))\ \ {\rm and}\ \   {\bf k}\mapsto (\Phi_{b_\star+1}({\bf x};{\bf k}),E_{b_\star+1}({\bf k})),\]
such that:
\begin{enumerate}
\item $E_D=E_{b_\star}({\bf K}_\star)=E_{b_\star+1}({\bf K}_\star)$ is a two-fold degenerate $L^2_{{\bf K}_\star}-$ eigenvalue of $\mathcal{L}$.
\item There exist functions
\[ \Phi_1({\bf x}) \in L^2_{{\bf K}_\star,\tau},\qquad \Phi_2({\bf x}) = \left(\mathcal{P}\mathcal{C}\right)[\Phi_1]({\bf x})=\overline{\Phi_1(-{\bf x})}\in L^2_{{\bf K}_\star,\bar\tau},\ \]
(see \eqref{L2_K_sigma}) and $\inner{\Phi_a, \Phi_b}_{L^2_{{\bf K}_\star}(\Omega_h)} = \delta_{ab}$, $a,b=1,2$.
such that
\[ \textrm{Nullspace}(\mathcal{L}-E_D I) =
{\rm span}\{ \Phi_1({\bf x}) , \Phi_2({\bf x})\}.\]
\item  Let $E_-({\bf k})=E_{b_*}({\bf k})$ and $E_+({\bf k})=E_{b_*+1}({\bf k})$.
There exist constants ${\upsilon_F} >0$ and $\zeta_0>0$, and Lipschitz continuous functions $e_\pm({\bf k})$, defined for $|{\bf k}-{\bf K}_\star|<\zeta_0$, such that $\abs{e_\pm({\bf k})|\le C|{\bf k}-{\bf K}_\star}$, and satisfying
\footnote{In  condensed matter physics  $v_{_F}$ is known as  the ``Fermi velocity''.}
\begin{equation}
\label{cones}
\begin{split}
E_+({\bf k})-E_D\ &=\ + {\upsilon_F}\
\left| {\bf k}-{\bf K}_\star \right|\
\left( 1\ +\ e_+({\bf k}) \right),\\
E_-({\bf k})-E_D\ &=\ - {\upsilon_F}\
\left| {\bf k}-{\bf K}_\star \right|\
\left( 1\ +\ e_-({\bf k}) \right) .
\end{split}
\end{equation}
\end{enumerate}
\end{definition}

\subsection{Multiplicity two $L^2_{{\bf K}_\star}-$ eigenvalues of $\mathcal{L}^A$ and conical singularities}

The following result, the analogue of Theorem 4.1 of \cites{FW:12} for Schr\"odinger operators, provides sufficient conditions for the existence of Dirac points for honeycomb-structured media defined by $\mathcal{L}^A$.

\begin{theorem}\label{prop_conical}
Let $\mathcal{L}^A=-\nabla\cdot A\nabla$, where $A({\bf x})$ defines a honeycomb structured medium.
Fix ${\bf K}_\star$, a vertex of $\mathcal{B}_h$.
Assume  that
\begin{enumerate}[label=(\subscript{A}{\arabic*})]
\item $\mathcal{L}^A$ has an $L^2_{{\bf K}_\star,\tau}-$ eigenvalue, $ E_D$, of multiplicity one, with corresponding eigenvector $\Phi_1({\bf x})$, normalized to have $L^2(\Omega_h)$ norm equal to one.
\item $\mathcal{L}^A$ has an $L^2_{{\bf K}_\star,\bar\tau}-$ eigenvalue, $ E_D$, of multiplicity one, with corresponding eigenvector $\Phi_2({\bf x})=\overline{\Phi_1(-{\bf x})}$.
\item $ E_D$ is not an $L^2_{{\bf K}_\star,1}-$ eigenvalue of $\mathcal{L}^A$.
\item the following non-degeneracy condition holds:
\begin{equation}\label{lamsharp_def}
\upsilon_{_F} =\frac12\left|
\overline{\inner{\Phi_1, \mathscr{A}\Phi_2}}\ \cdot\ \begin{pmatrix}1\\i\end{pmatrix}\right|>0,
\end{equation}
where $\mathscr{A}$ is defined in \eqref{A_def}.
\end{enumerate}

Then, we have the following:
\begin{enumerate}
\item $({\bf K}_\star, E_D)$ is a Dirac point in the sense of Definition \ref{dirac_pt_defn}.\medskip

\item
Assume further that $A({\bf x})$ defines  a real-valued honeycomb structured media.
 %
 %
 %
 Then,
 \begin{itemize}
 \item[(i)]\ $(E,\Phi({\bf x}))$ is an $L^2_{\bf K}$-  eigenpair of $\mathcal{L}^A$ if and only if $\left(E,(\mathcal{C}\Phi)({\bf x})\right)$ is an $L^2_{{\bf K}^\prime}$- eigenpair  of $\mathcal{L}^A$; here ${\bf K}^\prime=-{\bf K}$.
\item[(ii)]\ $({\bf K},E_D)$ is a Dirac point with associated normalized modes $\Phi_1^{\bf K}\in L^2_{{\bf K},\tau}$ and
 $\Phi_2^{\bf K}=\left(\mathcal{C}\circ\mathcal{P}\right) \Phi_1\in L^2_{{\bf K},\overline\tau}$, if and only if
 $({\bf K}^\prime,E_D)$ is a Dirac point with associated  normalized modes $\Phi_1^{{\bf K}^\prime}= \mathcal{C}\Phi_2^{\bf K} \in L^2_{{\bf K}^\prime,\tau}$ and
 $\Phi_2^{{\bf K}^\prime}=\mathcal{C} \Phi_1^{\bf K}\in L^2_{{\bf K}^\prime,\overline\tau}$\ .

\item[(iii)]\ Since $\mathcal{R}$ commutes with $\mathcal{L}^A$, we have that $(E,\Phi({\bf x};R^j{\bf K}))$,\  are  $L^2_{\bf K}-$ eigenpairs and $(E,\Phi({\bf x};R^j{\bf K}^\prime))$ are $L^2_{{\bf K}^\prime}-$ eigenpairs for $j=0,1,2$.
\item[(iv)] \ ${\upsilon_F}^{{\bf K}}={\upsilon_F}^{{\bf K'}}$.
\end{itemize}
\end{enumerate}
\end{theorem}

\begin{remark}\label{rem_upsilon_positive} The choice of normalized eigenfunctions $\{\Phi_1,\Phi_2=\mathcal{PC}\Phi_1\}$ in Theorem \ref{prop_conical} is unique once the phase of $\Phi_1$ is specified.  However, ${\upsilon_F}$ is independent of the choices of phase. We may use these phase degrees of freedom to select $\Phi_1\in L^2_{{\bf K},\tau}$ and $\Phi_2=\mathcal{PC}\Phi_1\in L^2_{{\bf K},\bar\tau}$ such that
\begin{equation}
{\upsilon_F}=\frac12
\overline{\inner{\Phi_1, \mathscr{A}{\Phi}_2}}\ \cdot\ \begin{pmatrix}1\\i\end{pmatrix}> 0.\label{lamsharp_pos}
\end{equation}
Indeed, given  $\tilde{\Phi}_1\in L^2_{{\bf K}_\star,\tau}$
 and $\tilde{\Phi}_2=\mathcal{PC}\tilde{\Phi}_1\in L^2_{{\bf K}_\star,\bar\tau}$ satisfying the properties of Theorem \ref{prop_conical}, define ${\lambda_\sharp}=\frac12
\overline{\inner{\tilde{\Phi}_1, \mathscr{A}{\mathcal{PC}\tilde{\Phi}}_1}}\ \cdot\ \begin{pmatrix}1\\i\end{pmatrix}$ and ${\lambda_\sharp}=|{\lambda_\sharp}|e^{i\arg\lambda_\sharp}$.
Set $\Phi_1=e^{-\frac{1}{2}i\arg\lambda_\sharp}\tilde{\Phi}_1$ and $\Phi_2=\left(\mathcal{C}\mathcal{P}\right)\Phi_1$. Then,
$\frac12\
\overline{\inner{\Phi_1, \mathscr{A}\Phi_2}}\ \cdot\ \begin{pmatrix}1\\i\end{pmatrix}=e^{-i\arg\lambda_\sharp}\ \frac12\ \overline{\inner{\tilde{\Phi}_1, \mathscr{A}\tilde{\Phi}_2}}\cdot \begin{pmatrix}1\\i\end{pmatrix}\
=|{\lambda_\sharp}|>0$.
Henceforth, we shall assume this choice of  $\Phi_1$ and $\Phi_2$.
\end{remark}

\begin{remark}\label{dp-dependence}
In the statement of Theorem \ref{prop_conical}, the mappings ${\bf k}\mapsto E_\pm$, ${\bf k}\mapsto\Phi_\pm({\bf x};{\bf k})$
 as well as ${\upsilon_F}$, depend on ${\bf K}_\star$, but this dependence has been suppressed.
   We shall at times make this dependence explicit by
 writing $E^{{\bf K}_\star}_\pm$, $\Phi^{{\bf K}_\star}_\pm({\bf x};{\bf k})$, and ${\upsilon_F}^{{\bf K}_\star}$.
\end{remark}

\begin{remark}\label{dp-nonreal}
Theorem \ref{prop_conical} states that if $A({\bf x})$ is real-valued honeycomb structured media, then the $E_D^{{\bf K}}=E_D^{{\bf K}'}$ and the associated eigenfunctions are related by symmetry.
For general complex-valued honeycomb structured media, this relation between the eigenspaces
 of Dirac points for ${\bf K}$ and ${\bf K}'$ does not hold.  A numerical illustration is  given in Appendix \ref{nonreal}.
\end{remark}

\begin{remark}[Wave-packets dynamics for data spectrally localized near Dirac points]\label{wavepkts}
The conical behavior of dispersion surfaces in a neighborhood of Dirac points suggests that the dynamics of wave-packets, which are initially spectrally localized near a Dirac point, are governed by a system of Dirac equations.
Specifically, consider the time-dependent 2D wave equation to which Maxwell's equations reduces (see Appendix \ref{maxwell})
\begin{equation}
\partial_t^2\psi+{\mathcal{L}^A}\psi\ =\ \left(\ \partial_t^2-\nabla\cdot A\nabla\ \right) \psi\ =\ 0.
\label{LA-wave}
\end{equation}
Using the approach of \cites{FW:14} we may construct solutions for initial conditions which are spectrally localized in a neighborhood of Dirac point $({\bf K}_\star,E_D)$. To leading order, such solutions have a two-scale structure
\begin{equation}
\psi({\bf x},t)\approx e^{-i\sqrt{E_D}\, t}\ \delta\ \left(\ \alpha_1(\bb \delta{\bf x},\delta t)\Phi_1({\bf x})+\alpha_2(\delta{\bf x},\delta t)\Phi_2({\bf x})\ \right).
\label{2scale-sol}\end{equation}
The expression in \eqref{2scale-sol} is of order $\delta^0$ in $H^s(\mathbb{R}^2),\ s\ge0$ and the corrector can be shown to be $o(1)$ in $H^s(\mathbb{R}^2)$
on large but finite time scales: $t\lesssim \mathcal{O}(\delta^{-2+})$. Here,
 $\alpha({\bf X},T)=\left(\alpha_1({\bf X},T),\alpha_2({\bf X},T)\right)^{\rm tr}$, ${\bf X}=\delta{\bf x}=(X_1,X_2)$ and $T=\delta t$,   satisfies 2-dimensional system of massless Dirac equations
 %
  %
\begin{equation}
i\ 2\sqrt{E_D}\ \partial_T\begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}
 =\ \upsilon_{_F}\begin{pmatrix}0&i\partial_{X_1}-\partial_{X_2}\\ i\partial_{X_1}+\partial_{X_2}&0\end{pmatrix}\ \begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}\ .
\end{equation}

This is in contrast to the dynamics of wave-packets, which are spectral concentrated near a spectral band edge which borders a spectral gap. Consider the case where the energy at the spectral band edge is a simple Floquet-Bloch eigenvalue. In this case, wave packet initial conditions which are spectral localized near the edge energy have envelope dynamics governed by an effective Schr\"odinger equation with non-zero effective mass; see, for example,
\cites{Ashcroft-Mermin:76,Allaire-Piatnitski:05,hoefer-weinstein:11,APR:11,KMOW:18} and \cite{BS:99,BS:01,BS:03,Suslina:04,BS:06,KP:01,KP:07}.
 %
A different asymptotic limit is the  geometrical optics / semi-classical regime, considered for Maxwell's equations, for example, in
 \cite{De_Nittis-Lein:14b,De_Nittis-Lein:17b}.
 \end{remark}
\medskip

\begin{proof}[Proof of Theorem \ref{prop_conical}]
The proof of part 2 of the theorem is straightforward.
By part 2, in order to prove part 1,  it suffices to prove all assertions for the vertex ${\bf K}_\star={\bf K}$.

The proof follows a Lyapunov-Schmidt reduction strategy; see also \cites{FW:12}.
Let ${\bf K}+{\bs \kappa}$, $|{\bs \kappa}|$ small, denote a quasi-momentum in a neighborhood of ${\bf K}$. We seek a non-trivial solution $(\phi,E)$ of the Floquet-Bloch eigenvalue problem:
\begin{equation}\label{k_evp}
\begin{split}
\mathcal{L}^A({\bf K}+{\bs \kappa})\phi({\bf x}) &= E\phi({\bf x}), \quad {\bf x} \in \mathbb{R}^2, \\
\phi({\bf x}+{\bf v}) &=\phi({\bf x}), \quad {\bf x}\in\mathbb{R}^2,\quad {\bf v}\in \Lambda_h,
\end{split}
\end{equation}
where
$\mathcal{L}^A({\bf k})=-(\nabla+i{\bf k})\cdot A({\bf x})(\nabla+i{\bf k})$;
 see \eqref{Eq_Eigen_k}-\eqref{Eq_Eigen_k_op}.

For ${\bs \kappa}$ small, we perturbatively construct solutions of the Floquet-Bloch eigenvalue problem \eqref{k_evp}.
 Expanding $\mathcal{L}^A({\bf K}+{\bs \kappa})$, we have
\begin{equation*}
\mathcal{L}^A({\bf K}+{\bs \kappa})=\mathcal{L}^A({\bf K})+{\bs \kappa}\cdot \mathscr{A}({\bf K})+{\bs \kappa}^T A({\bf x}){\bs \kappa} ,
\end{equation*}
where $\mathscr{A}$ is defined in \eqref{A_def} and
\begin{align}
\mathscr{A}({\bf k})&= e^{-i{\bf k}\cdot{\bf x}} \mathscr{A} e^{i{\bf k}\cdot{\bf x}} = A({\bf x}) \frac{1}{i}(\nabla_{\bf x}+i{\bf k}) + \frac{1}{i}(\nabla_{\bf x}+i{\bf k})\cdot A({\bf x}) . \label{Ak_def}
\end{align}

Let $ E^{(0)} = E_D$ denote the hypothesized double eigenvalue, and
\begin{equation}
 \label{phi1_exp}
 \phi^{(0)}({\bf x})=\alpha \phi_1({\bf x})+\beta\phi_2({\bf x}) ,
\end{equation}
with $\alpha$ and $\beta$ to be determined.
Introduce the orthogonal projections: $P_\parallel$, onto $\textrm{span}\{\phi_1,\phi_2\}$:
\[P_\parallel f\ \equiv\  \ \inner{\phi_1,f}\phi_1({\bf x})\ +\ \inner{\phi_2,f}\phi_2({\bf x}) ,
\] and $P_\perp=I-P_\parallel$.
We seek $ E({\bf K}+{\bs \kappa})$ and $\phi({\bf x};{\bf K}+{\bs \kappa})$ in the form:
\begin{equation}\label{Eq_Expan2}
 E({\bf K}+{\bs \kappa})= E_D+ E^{(1)},\quad \phi({\bf x};{\bf K}+{\bs \kappa})=\phi^{(0)}+\phi^{(1)},\ \ P_\parallel\phi^{(1)}=0.
\end{equation}
Substituting \eqref{Eq_Expan2} into the eigenvalue problem \eqref{k_evp} we obtain
\begin{equation}
\label{eqn_for_phi1}
(\mathcal{L}^A({\bf K})- E_D)\phi^{(1)}=\left(-{\bs \kappa}\cdot \mathscr{A}({\bf K})-{\bs \kappa}^T A {\bs \kappa} + E^{(1)}\right)(\phi^{(0)}+\phi^{(1)}), \quad \phi^{(1)}\in L^2(\mathbb{R}^2/\Lambda_h).
\end{equation}

Equation \eqref{eqn_for_phi1} may equivalently expressed as the system:
\begin{eqnarray}
(\mathcal{L}^A({\bf K})- E_D)\phi^{(1)}=P_\perp(-{\bs \kappa}\cdot \mathscr{A}({\bf K})-{\bs \kappa}^T A {\bs \kappa}+ E^{(1)})(\phi^{(1)}+\phi^{(0)}),\label{Eq_Sol1}\\
0=P_\parallel(-{\bs \kappa}\cdot \mathscr{A}({\bf K})-{\bs \kappa}^T A {\bs \kappa}+ E^{(1)})(\phi^{(1)}+\phi^{(0)})\label{Eq_Sol2}.
\end{eqnarray}

We first solve \eqref{Eq_Sol1}, for small ${\bs \kappa}$ and bounded $E^{(1)}$, to obtain $\phi^{(1)}$, which is linear in $\alpha$, $\beta$, with smooth dependence on ${\bs \kappa}$ and $E^{(1)}$.  We then substitute $\phi^{(1)}[\alpha,\beta,{\bs \kappa},E^{(1)}]$ into \eqref{Eq_Sol2} to obtain a homogeneous linear system for $\alpha, \beta$, depending on $E^{(1)}$ and ${\bs \kappa}$. The solvability condition of this system specifies the dependence of $E^{(1)}$ on ${\bs \kappa}$.

By the elliptic theory, $R_{_{\bf K}}(E_D)=(\mathcal{L}^A({\bf K})- E_D)^{-1}$ is defined as a bounded map from $P_\perp L^2(\mathbb{R}^2/\Lambda_h)$ to $P_\perp H^2(\mathbb{R}^2/\Lambda_h)$.
In addition, the mapping
\begin{align}
f\mapsto \Xi({\bs \kappa},E^{(1)}) f &\equiv R_{_{\bf K}}(E_D)\ P_\perp\left(-{\bs \kappa}\cdot \mathscr{A}({\bf K})-{\bs \kappa}^T A {\bs \kappa}+ E^{(1)}\right) f\nonumber
\end{align}
is a bounded operator on $H^s(\mathbb{R}^2/\Lambda_h)$, for any $s\geq0$.
Furthermore, for $|{\bs \kappa}|+| E^{(1)}|$ sufficiently small, the operator norm of $\Xi=\Xi({\bs \kappa},E^{(1)}) $ is less than one, and therefore $(I-\Xi)^{-1}$ is defined as a bounded operator on $H^s(\mathbb{R}^2/\Lambda_h)$. Therefore,  \eqref{Eq_Sol1} has a unique solution:
\begin{equation}
\phi^{(1)} =\hat{c}[{\bs \kappa}, E^{(1)}]\phi_1\ \alpha\ +\ \hat{c}[{\bs \kappa}, E^{(1)}]\phi_2\ \beta\ ,\ P_\perp\phi^{(1)}=\phi^{(1)},
\label{phi1-def}\end{equation}
where the operator $g\mapsto \hat{c}[{\bs \kappa}, E^{(1)}]\ g$,  defined by:
\begin{equation}
\hat{c}[{\bs \kappa}, E^{(1)}]\ g=\left[I-\Xi({\bs \kappa},E^{(1)})\right]^{-1}\ \Xi({\bs \kappa},E^{(1)})\ g ,
\label{c-def}\end{equation}
maps $H^s(\mathbb{R}^2/\Lambda_h)$ to $H^s(\mathbb{R}^2/\Lambda_h)$,\ $s\ge0$.
For $j=1,2$, $({\bs \kappa}, E^{(1)})\mapsto \hat{c}[{\bs \kappa}, E^{(1)}]\ \phi_j({\bf x})$ are smooth mappings from a neighborhood of $(0,0)$ in $\mathbb{R}^2\times \mathbb C$ into $H^2(\mathbb{R}^2/\Lambda_h)$ satisfying the bound:
\begin{equation}\label{Boundc}
\norm{\hat{c}[ E^{(1)},{\bs \kappa}]\ \phi_j}_{H^2(\mathbb{R}^2/\Lambda_h)}\le C\left(|{\bs \kappa}|+|E^{(1)}|\right), \quad j=1,2.
\end{equation}
and $P_\parallel \hat{c}[{\bs \kappa}, E^{(1)}]\phi_j=0$.

Substituting \eqref{phi1-def} into \eqref{Eq_Sol2}, we obtain a system of two homogeneous linear equations for $\alpha$ and $\beta$:
\begin{equation*}
\mathcal M( E^{(1)}, {\bs \kappa})\begin{pmatrix}\alpha\\ \beta\end{pmatrix}=0 ,
\end{equation*}
where $\mathcal M( E^{(1)}, {\bs \kappa})$ is a $2\times2$ matrix:
\begin{equation}\label{Mek}
\mathcal M( E^{(1)}, {\bs \kappa})= E^{(1)}I_{_{2\times2}}-\mathcal{M}_{\mathscr{A}}({\bs \kappa})-\mathcal{M}_R[ E^{(1)},{\bs \kappa}] .
\end{equation}

Recall that $\phi_j=e^{-i{\bf K}\cdot{\bf x}}\Phi_j$, $j=1,2$ and define
\begin{equation}
\mathcal{\hat{C}}[{\bs \kappa},E^{(1)}] g\ \equiv\ e^{i{\bf K}\cdot{\bf x}}\ \hat{c}[{\bs \kappa},E^{(1)}]\ e^{-i{\bf K}\cdot{\bf x}}\ g .
\end{equation}
Hence, $\inner{\Phi_j,\hat{C}[{\bs \kappa},E^{(1)}] g}=0,\ j=1,2$. Expressed in terms of $\Phi_j,\ j=1,2$, the matrices $\mathcal{M}_{\mathscr{A}}$
and $\mathcal{M}_{R}$ are given by:
\begin{equation}\label{M_A-def}
\mathcal{M}_{\mathscr{A}}({\bs \kappa})
=\begin{pmatrix}\inner{\Phi_1,{\bs \kappa}\cdot \mathscr{A} \Phi_1}&\inner{\Phi_1,{\bs \kappa}\cdot \mathscr{A} \Phi_2}\\
\inner{\Phi_2,{\bs \kappa}\cdot \mathscr{A} \Phi_1}&\inner{\Phi_2,{\bs \kappa}\cdot \mathscr{A} \Phi_2}\end{pmatrix} ,\ \textrm{and}
\end{equation}
\begin{align}\label{M_R-def}
&\mathcal M_R[ E^{(1)},{\bs \kappa}] \\
&=
\begin{pmatrix}
\inner{\Phi_1,\left[\left({\bs \kappa}\cdot\mathscr{A}+{\bs \kappa}^T A{\bs \kappa}\right)\hat{C}+
{\bs \kappa}^T A{\bs \kappa}\right]\Phi_1} &
\inner{\Phi_1,\left[\left({\bs \kappa}\cdot\mathscr{A}+{\bs \kappa}^T A{\bs \kappa}\right)\hat{C}+
{\bs \kappa}^T A{\bs \kappa}\right]\Phi_2} \\
\inner{\Phi_2,\left[\left({\bs \kappa}\cdot\mathscr{A}+{\bs \kappa}^T A{\bs \kappa}\right)\hat{C}+
{\bs \kappa}^T A{\bs \kappa}\right]\Phi_1}
&
\inner{\Phi_2,\left[\left({\bs \kappa}\cdot\mathscr{A}+{\bs \kappa}^T A{\bs \kappa}\right)\hat{C}+
{\bs \kappa}^T A{\bs \kappa}\right]\Phi_2}
 \end{pmatrix},
\nonumber\end{align}
where $\hat{C}=\hat{C}[{\bs \kappa},E^{(1)}]$. Furthermore,
by \eqref{Boundc},
\begin{equation} \mathcal M_{R,lm}[ E^{(1)},{\bs \kappa}] = \mathcal{O}\left(|{\bs \kappa}|\left(|{\bs \kappa}|+E^{(1)}\right)\right),\ l,m=1,2.\label{M_R-est}\end{equation}

We have therefore obtained the following result which characterizes the dispersion surfaces in a neighborhood of Dirac points:
\begin{proposition}\label{det0} For $|{\bs \kappa}|$ sufficiently small,
$ E= E_D+ E^{(1)}$ is a $L^2_{_{{\bf K}+{\bs \kappa}}}-$ eigenvalue problem \eqref{k_evp} if and only if
\begin{equation}\label{Det}
\det\mathcal M( E^{(1)},{\bs \kappa})=0,
\end{equation}
where $M( E^{(1)},{\bs \kappa})$ is given by \eqref{Mek}, \eqref{M_A-def} and \eqref{M_R-def}.
\end{proposition}

We shall next exploit  symmetry to obtain simplified expressions for the entries of $\mathcal{M}_{\mathcal{A}}({\bs \kappa})$.
\begin{proposition}\label{prop_g}
Let ${\bf K}_\star$ denote any vertex of the Brillouin zone, $\mathcal{B}_h$. Let $A({\bf x})$ define a honeycomb structured medium; see Section \ref{honey-media}.
 Let  $\Phi_1$ and $\Phi_2$ be as  in the Theorem \ref{prop_conical}, $\mathscr{A}$ be given by \eqref{A_def} and ${\bs \kappa}=(\kappa^{(1)},\kappa^{(2)})\in\mathbb{R}^2-\{\bs 0\}$. Then,
\begin{enumerate}
\item $\inner{\Phi_1,{\bs \kappa}\cdot \mathscr{A} \Phi_1}
=\inner{\Phi_2,{\bs \kappa}\cdot \mathscr{A} \Phi_2}=0$.
\item $\inner{\Phi_1,{\bs \kappa}\cdot \mathscr{A} \Phi_2}
=\overline{\inner{\Phi_2,{\bs \kappa}\cdot \mathscr{A} \Phi_1}}=\upsilon_{_F}(\kappa^{(1)}+i\kappa^{(2)}),$
where $\upsilon_{_F}$ is given by \eqref{lamsharp_def}.
\end{enumerate}
Thus,
\begin{equation}
\mathcal{M}_{\mathscr{A}}({\bs \kappa})
=\upsilon_{_F}\ \begin{pmatrix}0 & \kappa^{(1)}+i\kappa^{(2)}\\
\kappa^{(1)}-i\kappa^{(2)}&0\end{pmatrix}.
\label{MA-simp}\end{equation}
\end{proposition}

\begin{proof}[Proof of Proposition \ref{prop_g}]
Without loss of generality, let ${\bf K}_\star={\bf K}$.
 Since $\Phi_1\in L^2_{{\bf K},\tau}$ and $\Phi_2\in L^2_{{\bf K},\overline{\tau}}$, we have $\mathcal{R}\Phi_1=\tau\Phi_1$ and $\mathcal{R}\Phi_2=\overline{\tau}\Phi_2$.
Let $\nu_1=\tau$ and $\nu_2=\overline{\tau}$.
 For $l,m=1,2$ consider the column vector $\inner{\Phi_l,\mathscr{A}\Phi_m}$.
Applying  Lemma \ref{RcomA} we have
\begin{equation*}
\inner{\Phi_l,\mathscr{A}\Phi_m}=\inner{\mathcal R \Phi_l, \mathcal R \mathscr{A} \Phi_m}=\inner{\mathcal R \Phi_l, R^*\mathscr{A} \mathcal R \Phi_m}=\overline{\nu_l}\nu_m R^*\inner{\Phi_l,\mathscr{A} \Phi_m},\ l,m=1,2\ .
\end{equation*}
Hence, $R \inner{\Phi_l,\mathscr{A}\Phi_m} = \overline{\nu_l}\nu_m\ \inner{\Phi_l,\mathscr{A}\Phi_m}$.

If $l=m$, then $\overline{\nu_l}\nu_m=|\nu_m|^2=1$. Therefore, $R\inner{\Phi_m,\mathscr{A}\Phi_m}_{L^2_{\bf K}}=\inner{\Phi_m,\mathscr{A}\Phi_m}$. Since $1$ is not an eigenvalue of  $R$, it follows that
\[\inner{\Phi_1,\mathscr{A}\Phi_1}_{L^2_{\bf K}}=\inner{\Phi_2,\mathscr{A}\Phi_2}_{L^2_{\bf K}}=0.\]
 If $l=1$ and $m=2$, then $R\inner{\Phi_1,\mathscr{A}\Phi_2}_{L^2_{\bf K}}=\tau \inner{\Phi_1,\mathscr{A}\Phi_2}_{L^2_{\bf K}}$. By \eqref{R_evals_evecs},
 \begin{equation}
\label{lambda_inner_prod}
\inner{\Phi_1, \mathscr{A}\Phi_2}_{L^2_{\bf K}}=\upsilon_{_F}\begin{pmatrix}1\\i\end{pmatrix}.
\end{equation}
Note, by hypothesis $(A4)$ and Remark \ref{rem_upsilon_positive}, we may take $\upsilon_{_F}>0$.
 Finally, if $l=2$ and $m=1$
\begin{equation*}
  \inner{\Phi_2,\mathscr{A}\Phi_1}=\inner{\mathscr{A}\Phi_2,\Phi_1}=\overline{\inner{\Phi_1,\mathscr{A}\Phi_2}}
  ={\upsilon_{_F}}\begin{pmatrix}1\\-i\end{pmatrix} .
\end{equation*}
This completes the proof of Proposition \ref{prop_g}.
\end{proof}

It follows that the eigenvalue condition  $\det\mathcal{M}( E^{(1)},{\bs \kappa}) =0$ (\eqref{Det}) is of the form:
$\left( E^{(1)}\right)^2={\upsilon_F}^2|{\bs \kappa}|^2+g_{21}(E^{(1)},{\bs \kappa})+g_{12}(E^{(1)},{\bs \kappa})+g_{03}({\bs \kappa}),$
where the functions $g_{rs}$ are smooth and satisfy the bounds
$
|g_{rs}(E^{(1)},{\bs \kappa})|\le C |E^{(1)}|^r\ |{\bs \kappa}|^s,
$
for $|E^{(1)}|\le1,$ $|{\bs \kappa}|\le1$.
 The  proof of Proposition \ref{prop_conical} is  completed along the lines of Proposition 4.2 of \cites{FW:12} and yields
 two locally conical solution branches: $E_+({\bf K}+{\bs \kappa})=E_D+E^{(1)}_{+}({\bs \kappa})=E_D+{\upsilon_{_F}}\ |{\bs \kappa}|(1+e_+({\bf K}+{\bs \kappa}))$ and $E_-({\bf K}+{\bs \kappa})=E_D+E^{(1)}_{-}({\bs \kappa})=E_D-{\upsilon_{_F}}\ |{\bs \kappa}|(1+e_-({\bf K}+{\bs \kappa}))$.
\end{proof}

\subsection{Dirac points in low-contrast media}\label{sec:dp-low}

Any low-contrast honeycomb structured medium in the sense of Section \ref{honey-media}
is of the form  $A^{(\varepsilon)}({\bf x})=a_0I+\varepsilon A^{(1)}({\bf x})$ where $a_0>0$ and where $\varepsilon/a_0 >0$ is small.
Indeed, suppose $A^{(\varepsilon)}({\bf x})=A^{(0)}+\varepsilon A^{(1)}({\bf x})$ with $A^{(0)}$ constant. By Theorem \ref{invariance}, $A^{(0)}=a_0I+b_0\sigma_2$ where $a_0>0$ and $b_0$ is an odd function. But $b_0$ constant and odd implies that $b_0=0$.

In this section we apply Theorem \ref{prop_conical} to the study of Dirac points for low-contrast structures. We consider the Floquet-Bloch eigenvalue problem  for the operator
$\mathcal{L}^{(\varepsilon)}=-\nabla\cdot A^{(\varepsilon)}({\bf x})\nabla$. By linearity, we may make the replacements: $A^{(\varepsilon)}({\bf x}) \to  \frac{1}{a_0}A^{(\varepsilon)}({\bf x})$ and $E\to  E/a_0$. Hence without loss of generality, we assume
\begin{equation}
\label{normalized_honeycomb}
 A^{(\varepsilon)}({\bf x})=I+\varepsilon A^{(1)}({\bf x}) .
\end{equation}
For small $\varepsilon$, we study ${\bf K}-$ pseudo-periodic eigenvalue problem
\begin{equation}
\label{eps_evp}
\mathcal{L}^{(\varepsilon)} \Phi= E\Phi, \quad \Phi\in L^2_{{\bf K}}(\mathbb{R}^2/\Lambda_h),
\end{equation}
where
\begin{align}\label{L_eps_def}
\mathcal{L}^{(\varepsilon)}\equiv\mathcal{L}^{A^{(\varepsilon)}} &=
 -\Delta-\varepsilon\nabla\cdot A^{(1)}\nabla = -\Delta+\varepsilon \mathcal{L}^{(1)}.
\end{align}
We shall solve the eigenvalue problem \eqref{eps_evp} for small $\varepsilon$. We begin by summarizing the  relevant spectral properties of $\mathcal{L}^{(0)}=-\Delta$ in  $L^2_{{\bf K}}(\mathbb{R}^2/\Lambda_h)$; see \cites{FW:12,berkolaiko-comech:15}:

\begin{proposition}
\label{eps_0_prop}
\begin{enumerate}
\item $E^{(0)}=|{\bf K}|^2$ is an $L^2_{\bf K}-$ eigenvalue of $-\Delta$ of multiplicity of three with corresponding three-dimensional eigenspace:
\begin{equation*}
  \emph{span}\left\{e^{i{\bf K}\cdot {\bf x}},e^{iR{\bf K}\cdot{\bf x}},e^{iR^2{\bf K}\cdot{\bf x}}\right\}.
\end{equation*}
\item Considered in each of the three orthogonal $\mathcal{R}$-invariant subspaces, $L^2_{{\bf K},\sigma}$, $\sigma=1,\tau,\bar{\tau}$, $ E^{(0)}=|{\bf K}|^2$ is a simple eigenvalue of
$-\Delta$ with corresponding eigenspace spanned by the $L^2(\mathbb{R}^2/\Lambda_h)$ normalized eigenvector:
\begin{equation}
  \Phi^{(0)}_\sigma \equiv \frac{1}{\sqrt{3|\Omega_h|}}\left[e^{i{\bf K}\cdot{\bf x}}+\bar{\sigma} e^{iR{\bf K}\cdot{\bf x}}+\sigma e^{iR^2{\bf K}\cdot{\bf x}}\right] \in\ L^2_{{\bf K},\sigma} .
\label{Phi0sig}\end{equation}
\end{enumerate}
\end{proposition}

We next turn to the eigenvalue problem \eqref{eps_evp} for $\varepsilon$ small.

\begin{theorem}\label{low-dp}
Let $A^{(\varepsilon)}({\bf x})=I_{_{2\times2}}+\varepsilon A^{(1)}({\bf x})$ define a honeycomb structured medium, with Fourier expansion characterized by \eqref{AAfourier} (Corollary \ref{A-expand}).
Assume that $A^{(\varepsilon)}({\bf x})$ satisfies the non-degeneracy condition:
\begin{equation}
{\bf K}^T\ A_{0,-1}\ R{\bf K}\neq0\label{non-deg}.
\end{equation}

Then, there exists $\varepsilon^0>0$, and mappings
$\varepsilon\mapsto E_D^\varepsilon$ and $\varepsilon\mapsto \widetilde{E}^\varepsilon$,
 $\varepsilon\mapsto \Phi_1^{\varepsilon}\in L^2_{{\bf K},\tau}$,  $\varepsilon\mapsto \Phi_2^{\varepsilon}\in L^2_{{\bf K},\bar\tau}$
  and $\varepsilon\mapsto \widetilde{\Phi}^{(\varepsilon)}\in L^2_{{\bf K},1}$
 such that the following holds for all $\varepsilon\in(-\varepsilon^0,\varepsilon^0)$:
 %
 %
 \begin{enumerate}
\item $ E^{\varepsilon}_D$ is a simple $L^2_{{\bf K},\tau}-$ eigenvalue of $\mathcal{L}^{(\varepsilon)}$ with
\begin{equation}
 E_D^\varepsilon= \ |{\bf K}|^2+\varepsilon\left({\bf K}^T\ A^{(1)}_{0,0}\ {\bf K}-{\bf K}^T\ A^{(1)}_{0,-1}\ R{\bf K}\ \right)+\mathcal{O}(\varepsilon^2) ,\label{dp-eps}
\end{equation}
with eigenspace spanned by states $\Phi^{\varepsilon}_1\in L^2_{{\bf K},\tau}$
\item $ E_D^\varepsilon$ is also a simple $L^2_{{\bf K},\overline\tau}-$ eigenvalue of $\mathcal{L}^{(\varepsilon)}$
with corresponding eigenfunction $\Phi^{\varepsilon}_2=\left(\mathcal{C}\circ\mathcal{P}\right) \Phi_1^\varepsilon\in L^2_{{\bf K},\bar{\tau}}$.
\item $\widetilde{ E}^\varepsilon\ne E_D^\varepsilon$ is a simple $L^2_{{\bf K}, 1}-$ eigenvalue of $\mathcal{L}^{(\varepsilon)}$:
\begin{equation}
\widetilde{ E}^{\varepsilon}=|{\bf K}|^2
+\varepsilon\ \left({\bf K}^T A^{(1)}_{0,0}\ {\bf K}+2{\bf K}^T\ A^{(1)}_{0,-1}\ R{\bf K}\  \right)+\mathcal{O}(\varepsilon^2),
\label{Etilde-eps}\end{equation}
with corresponding eigenspace spanned by $\tilde{\Phi}^{(\varepsilon)}\in L^2_{{\bf K},1}$.
\item
\begin{align}
{\upsilon^{\bf K}_{_F}}(\varepsilon)=\left| \frac12 \overline{\inner{\Phi^{\varepsilon}_1,\mathscr{A^{(\varepsilon)}}\Phi^{\varepsilon}_2}} \cdot \begin{pmatrix}1\\i\end{pmatrix} \right|= \frac{4\pi}{3}+\mathcal{O}(\varepsilon).
\label{vF-eps}
\end{align}
Therefore, $({\bf K},E^\varepsilon_D)$ is a Dirac point in the sense of Definition \ref{dirac_pt_defn}.
\item \subitem If $\varepsilon{\bf K}^T\ A^{(1)}_{0,-1}\ R{\bf K}>0$, the Dirac points occurs at the intersection of the $1^{st}$ and $2^{nd}$ dispersion surfaces at the vertices of $\mathcal{B}_h$. \subitem If $\varepsilon{\bf K}^T\ A^{(1)}_{0,-1}\ R{\bf K}<0$, Dirac points occur at the intersection of the $2^{nd}$ and $3^{rd}$ dispersion surfaces at the vertices of $\mathcal{B}_h$.
\end{enumerate}
\end{theorem}

\begin{proof}
Statements (1)-(4) of the theorem imply that for $\varepsilon\in(-\varepsilon_0,\varepsilon_0)$, assumptions $(A_1)$-$(A_4)$ of Theorem \ref{prop_conical} hold. These statements therefore imply, by Proposition \ref{prop_conical}, the existence of Dirac points $({\bf K}_\star,E^{\varepsilon}_D))$ at the vertices of $\mathcal{B}_h$.

Consider the eigenvalue problem in $L^2_{{\bf K},\sigma}$ for $\sigma=1,\tau, \bar{\tau}$:
\begin{equation}\label{Eq_EigenL}
(-\Delta+\varepsilon \mathcal{L}^{(1)})\Phi_\sigma^{(\varepsilon)}= E_\sigma^{\varepsilon}\Phi_\sigma^{(\varepsilon)}, \quad \Phi_\sigma^{(\varepsilon)} \in L^2_{{\bf K},\sigma},
\end{equation}
where $\mathcal{L}^{(1)}=-\nabla\cdot A^{(1)}\nabla$.
We seek $ E_\sigma^{\varepsilon}$ and $\Phi_\sigma^{(\varepsilon)}$ in form of expansions:
\begin{equation}\label{Eq_Expan1}
   E_\sigma^{\varepsilon}= E^{(0)} + \varepsilon E_\sigma^{(1,\varepsilon)},\quad \Phi_\sigma^\varepsilon=\Phi_\sigma^{(0)} + \varepsilon \Phi_\sigma^{(1)},\quad
   \left\langle \Phi_\sigma^{(0)},\Phi_\sigma^{(1)}\right\rangle=0,
\end{equation}
where $ E^{(0)}=|{\bf K}|^2$ (Proposition \ref{eps_0_prop}) and $\Phi_\sigma^{(0)}\in L^2_{{\bf K},\sigma}$, given by the expression in
\eqref{Phi0sig},
satisfies $-\Delta\Phi_\sigma^{(0)}= E^{(0)}\Phi_\sigma^{(0)}$. Then, $\Phi_\sigma^{(1)}$ satisfies the non-homogeneous equation:
\begin{equation}
\label{F_eqn}
\left(-\Delta- E^{(0)}\right)\Phi_\sigma^{(1)}=\left( -\mathcal{L}^{(1)}+ E^{(1)}\right)\left(\Phi_\sigma^{(0)}+\varepsilon\Phi_\sigma^{(1)}\right) .
\end{equation}

 Introduce the orthogonal projections: $P_\parallel$ onto span$\{\Phi_\sigma^{(0)}\}$, the $L^2_{{\bf K},\sigma}-$ nullspace of $(-\Delta- E^{(0)})$,  and $P_{\perp}=I-P_\parallel$.
Equation \eqref{F_eqn} may then be rewritten as the following equivalent system for $\Phi_\sigma^{(1)}$ and $ E^{(1)}$:
\begin{eqnarray}
\left(-\Delta- E^{(0)}\right)\Phi_\sigma^{(1)}=P_\perp \left( -\mathcal{L}^{(1)}+ E^{(1)}\right)\left(\Phi_\sigma^{(0)}+\varepsilon\Phi_\sigma^{(1)}\right),\label{Eq_Phi1}\\
0=P_\parallel \left( -\mathcal{L}^{(1)}+ E^{(1)}\right)\left(\Phi_\sigma^{(0)}+\varepsilon\Phi_\sigma^{(1)}\right)\label{Solvability}.
\end{eqnarray}
In analogy with the reduction strategy for the proof of Proposition \ref{prop_conical}, we first solve \eqref{Eq_Phi1} to obtain a solution for $\Phi_\sigma^{(1)}$ as a smooth functional of $\varepsilon$ and $E^{(1)}$, and then substitute the solution, $\Phi_\sigma^{(1)}=\Phi_\sigma^{(1)}[\varepsilon,E^{(1)}]$,  into \eqref{Solvability} to obtain a closed equation for $E^{(1)}$ as a function of $\varepsilon$.

Because $\left(-\Delta- E^{(0)}\right)^{-1}$ is a bounded operator from $P_\perp L_{{\bf K}, \sigma}^2$ to $P_\perp H_{{\bf K},\sigma}^2$, equation \eqref{Eq_Phi1} may be rewritten as
\begin{equation}\label{Eq_Phi12}
  \left(I-\varepsilon\mathcal{W}(E^{(1)})\right)\Phi^{(1)}_\sigma
  =  -\left(-\Delta- E^{(0)}\right)^{-1}P_\perp\mathcal{L}^{(1)} \Phi^{(0)}_\sigma,
\end{equation}
where the operator $\mathcal{W}(E^{(1)})$, defined by:
\begin{equation*}
 f\mapsto \mathcal{W}(E^{(1)}) f\equiv  \left(-\Delta- E^{(0)}\right)^{-1} P_\perp\left(-\mathcal{L}^{(1)} + E^{(1)}\right),
\end{equation*}
is a bounded on $H^s_{{\bf K},\sigma}$, for any $s\geq0$.
Furthermore, for $|\varepsilon|+ |E^{(1)}|$ sufficiently small, the operator norm of $\varepsilon\mathcal{W}$ is less than one. Hence,  $(1-\varepsilon\mathcal{W}(E^{(1)}))^{-1}$ exists, is bounded on $H^s_{{\bf K},\sigma}$,  and equation \eqref{Eq_Phi12} is uniquely solvable in $H^2_{{\bf K}}$:
\begin{equation}\label{Eq_Phi123}
\Phi^{(1)}_\sigma[\varepsilon,E^{(1)}]= - \left(I-\varepsilon\mathcal{W}(E^{(1)})\right)^{-1} \left(-\Delta- E^{(0)}\right)^{-1} P_\perp \mathcal{L}^{(1)}\Phi_\sigma^{(0)}\ .
\end{equation}

Substituting \eqref{Eq_Phi123} into the solvability equation \eqref{Solvability}, we obtain a closed equation for $E^{(1)}$:
\begin{equation}
\label{G_sigma_def}
G_\sigma(\varepsilon, E^{(1)}) \equiv \inner{\Phi^{(0)}_\sigma, \left(-\mathcal{L}^{(1)}+ E^{(1)} \right) \left(\Phi^{(0)}_\sigma + \varepsilon\Phi^{(1)}_\sigma[\varepsilon,E^{(1)}] \right)}_{L^2_{{\bf K}}} = 0.
\end{equation}
For $\varepsilon=0$, equation \eqref{G_sigma_def} reduces to
$G_\sigma(0, E^{(1)}) = -\inner{\Phi^{(0)}_\sigma, \mathcal{L}^{(1)} \Phi^{(0)}_\sigma} + E^{(1)} = 0$. Thus,
$G_\sigma(0, E_\sigma^{(1,0)})=0$, where
\begin{equation}
E_{\sigma}^{(1,0)}\equiv \inner{\Phi^{(0)}_\sigma, \mathcal{L}^{(1)} \Phi^{(0)}_\sigma}.
\label{Esigma01}\end{equation}
  Moreover,
$\partial_{E^{(1)}} G_\sigma(0,E_{\sigma}^{(1,0)})=1\ne0$. Hence, by the implicit function theorem, there is a unique mapping $\varepsilon\mapsto E_\sigma^{(1,\varepsilon)}$, defined and analytic in a neighborhood of $\varepsilon=0$, which satisfies $E_\sigma^{(1,\varepsilon=0)}=E^{(1,0)}_{\sigma}$ and  $G_\sigma(\varepsilon, E_\sigma^{(1,\varepsilon)})=0$, for all  $\varepsilon$ in a complex neighborhood of $\varepsilon=0$.

Assertions 1-3 of  Theorem \ref{low-dp} are now a direct consequence of the evaluation of $E^{(1,0)}_{\sigma}$, given in the following:
\begin{proposition}\label{prop_per2}
Recall that $\mathcal{L}^{(\varepsilon)}= -\Delta + \varepsilon\mathcal{L}^{(1)}=-\Delta-\varepsilon\nabla\cdot A^{(1)}\nabla$, and
for ${\bf m}=(m_1,m_2)\in\mathbb{Z}^2$, denote by $A^{(1)}_{\bf m}$ the Fourier coefficients of $A^{(1)}({\bf x})$:
\begin{equation*}
  A^{(1)}_{\bf m} \equiv \frac{1}{|\Omega_h|} \int_{\Omega_h} e^{-i{\bf m}\vec{\bf k}\cdot \,{\bf y}} A^{(1)}({\bf y}) d{\bf y}.
\end{equation*}
We have,
\begin{equation}
E^\varepsilon_\sigma\ =\ |{\bf K}|^2\ +\ \varepsilon\ E_\sigma^{(1,\varepsilon)}\ =\ |{\bf K}|^2\ +\ \varepsilon\ E_\sigma^{(1,0)}+\mathcal{O}(\varepsilon^2),
\label{Eeps-sigma}
\end{equation}
where the following assertions concerning $E^{(1,0)}_{\sigma}$ hold:
\begin{enumerate}
\item If $[\mathcal{PC},\mathcal{L}^{(1)}]=0$ and $[\mathcal{R},\mathcal{L}^{(1)}]=0$ (hence $A^{(1)}(R^*{\bf x})=R^*A^{(1)}({\bf x})R$), then
\begin{equation}
\label{phi_sigma_identity}
E^{(1,0)}_{\sigma}\ =\ \inner{\Phi^{(0)}_\sigma, \mathcal{L}^{(1)} \Phi^{(0)}_\sigma} ={\bf K}^T\ A^{(1)}_{0,0}\ {\bf K}\ +\
 (\sigma+\overline{\sigma})\ {\bf K}^T A^{(1)}_{0,-1}\ R{\bf K}
.
\end{equation}
\item Assume further that $A^{(1)}(R^*{\bf x})=A^{(1)}({\bf x})$. (Therefore, by Theorem \ref{invariance},  $A^{(1)}({\bf x})=a^{(1)}({\bf x})\ I_{_{2\times2}}\ +\ b^{(1)}({\bf x})\sigma_2$ with $a^{(1)}(-{\bf x})=a^{(1)}({\bf x})$ and $b^{(1)}(-{\bf x})=-b^{(1)}({\bf x})$, and then $A^{(1)}_{0,-1}=a_{0,-1}^{(1)}({\bf x})\ I_{_{2\times2}}\ +\ b_{0,-1}^{(1)}({\bf x})\sigma_2$ ) \\
Then,
\begin{equation}
E^{(1,0)}_{\sigma}\ =\ \inner{\Phi^{(0)}_\sigma, \mathcal{L}^{(1)} \Phi^{(0)}_\sigma}=|{\bf K}|^2\left(a^{(1)}_{0,0}-\left(\frac{1}{2}a^{(1)}_{0,-1}-\frac{\sqrt{3}}{2}ib^{(1)}_{0,-1}\right)(\sigma+\overline{\sigma})\right).
\label{ab11} \end{equation}
 by ${\bf K}^T \ R {\bf K} =-\frac{1}{2}|{\bf K}|^2$ and ${\bf K}^T \ \sigma_2 R{\bf K} = \frac{\sqrt{3}}{2} i$.
 %
 \end{enumerate}

\end{proposition}

We prove Proposition \ref{prop_per2} below. We first conclude the proof of Theorem \ref{low-dp} by verifying the assertion in part 4.
  By \eqref{Eq_Expan1} and \eqref{Eq_Phi123},
\begin{equation}
\label{phi_eps}
 \Phi_\sigma^{(\varepsilon)} = \Phi^{(0)}_\sigma + \mathcal{O}(\varepsilon),
\end{equation}
where $\Phi^{(0)}_\sigma$ is given in \eqref{Phi0sig}.
Note that $\Phi_\tau^{(\varepsilon)}=\Phi^{\varepsilon}_1$ and $\Phi_{\overline{\tau}}^{(\varepsilon)}=\Phi^{\varepsilon}_2$.
Substituting \eqref{phi_eps} into the definition of $\upsilon_F^{{\bf K}}$ (\eqref{lamsharp_def}) and recalling that $A({\bf x})=I_{_{2\times2}}+\varepsilon A^{(1)}({\bf x})$ (\eqref{normalized_honeycomb}) yields:
\begin{align*}
\upsilon_{_F}^{{\bf K}}(\varepsilon) &=
 \frac12\left| \overline{\inner{\Phi^{\varepsilon}_1,\mathscr{A}^\varepsilon\Phi^{\varepsilon}_2}_{L^2_{\bf K}}}\cdot\begin{pmatrix}1\\i\end{pmatrix} \right|\\
&= \left|\overline{\inner{\Phi_\tau^{(0)},-i\nabla\Phi_{\overline{\tau}}^{(0)}}}_{L^2_{\bf K}} \cdot \begin{pmatrix}1\\i\end{pmatrix}\right| + \mathcal{O}(\varepsilon) \\
&=\left| \frac{1}{3} \left(1 + \tau R + \bar\tau R^2\right){\bf K} \cdot \begin{pmatrix}1\\i\end{pmatrix}\right| + \mathcal{O}(\varepsilon)
\\
&= \left|\frac{1}{3} \left(\tau{\bf k}_2 - \bar\tau{\bf k}_1 \right) \cdot \begin{pmatrix}1\\i\end{pmatrix} \right|+ \mathcal{O}(\varepsilon) \\
&= \left|\frac{2\pi}{3} \begin{pmatrix}i\\1\end{pmatrix} \cdot \begin{pmatrix}1\\i\end{pmatrix} \right|+ \mathcal{O}(\varepsilon)\ =\ \frac{4\pi}{3} + \mathcal{O}(\varepsilon) .
\end{align*}
It follows that  for $\varepsilon$ sufficiently small, $\upsilon_{_F}^{{\bf K}}(\varepsilon)\ne0$.
This completes the proof of part 4, and therewith Theorem \ref{low-dp}.
\end{proof}

It remains to prove Proposition \ref{prop_per2}.

\begin{proof}[Proof of Proposition \ref{prop_per2}]
Since $[\mathcal{R},\mathcal{L}^{(1)}]=0$,
\begin{equation}
 \alpha \equiv \left\langle e^{i{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{i{\bf K}\cdot {\bf x}}\right\rangle =
 \left\langle e^{iR{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR{\bf K}\cdot {\bf x}}\right\rangle = \left\langle e^{iR^2{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR^2{\bf K}\cdot {\bf x}}\right\rangle .
\end{equation}
Note that
$
\alpha=\left\langle e^{i{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{i{\bf K}\cdot {\bf x}}\right\rangle =\int_{\Omega_h}{\bf K}\cdot A^{(1)}({\bf x}){\bf K} d{\bf x}=|\Omega_h|\ {\bf K}^T\ A^{(1)}_{0,0}\ {\bf K}.
$
Again since  $[\mathcal{R},\mathcal{L}^{(1)}]=0$,
\begin{equation}
\beta\equiv \left\langle e^{i{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR{\bf K}\cdot {\bf x}}\right\rangle = \left\langle
e^{iR{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR^2{\bf K}\cdot {\bf x}}
\right\rangle = \left\langle  e^{iR^2{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{i{\bf K}\cdot {\bf x}}\right\rangle,
\end{equation}
and
\begin{equation}
 \gamma \equiv \left\langle e^{i{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR^2{\bf K}\cdot {\bf x}}\right\rangle = \left\langle e^{iR{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{i{\bf K}\cdot {\bf x}}\right\rangle = \left\langle e^{iR^2{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR{\bf K}\cdot {\bf x}}\right\rangle .
\end{equation}

Since $\mathcal{L}^{(1)}$ is self-adjoint,
 we have
 $
\gamma=\left\langle e^{i{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR{\bf K}\cdot {\bf x}}\right\rangle =\overline{\left\langle e^{iR{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{i{\bf K}\cdot {\bf x}}\right\rangle}=\overline{\beta}$.
Moreover,
\begin{equation}
\beta=\left\langle e^{i{\bf K}\cdot {\bf x}}, \mathcal{L}^{(1)}e^{iR{\bf K}\cdot {\bf x}}\right\rangle=\int_{\Omega_h}{\bf K}\cdot A^{(1)}({\bf x})R{\bf K} e^{i{\bf k}_2\cdot {\bf x}}d{\bf x}=|\Omega_h|{\bf K}^T\ A^{(1)}_{0,-1}\ R{\bf K},
\end{equation}
where we have used $\mathbf{R}{\bf K}-{\bf K}={\bf k}_2$.
From the $\mathcal{PC}$ symmetry, we know that $A^{(1)}_{0,-1}$ is real by Corollary \ref{A-expand}. So $\gamma=\beta$ is real.

The preceding discussion yields
\begin{equation}
\begin{split}
\left\langle \Phi^{(0)}_\sigma, \mathcal{L}^{(1)}\Phi^{(0)}_\sigma \right\rangle
&=\frac{1}{3|\Omega_h|}\begin{pmatrix}1&\sigma&\overline{\sigma}\end{pmatrix}\begin{pmatrix}\alpha&\beta&\beta\\ \beta&\alpha&\beta\\ \beta&\beta&\alpha\end{pmatrix}\begin{pmatrix}1\\ \overline{\sigma}\\ \sigma\end{pmatrix}\\
&=\frac{1}{|\Omega_h|}(\alpha+\beta(\sigma+\overline{\sigma}))\\&={\bf K}\cdot A^{(1)}_{0,0}{\bf K}+{\bf K}\cdot A^{(1)}_{0,-1}R{\bf K}(\sigma+\overline{\sigma}).
\end{split}
\end{equation}
For $\sigma=\tau$ or $\bar\tau$, $E_\sigma^{(1,0)}\equiv \left\langle \Phi^{(0)}_\sigma, \mathcal{L}^{(1)}\Phi^{(0)}_\sigma \right\rangle = {\bf K}^T\ A^{(1)}_{0,0}\ {\bf K}-{\bf K}^T\  A^{(1)}_{0,-1}R\ {\bf K}.$
\\ For $\sigma=1$,
$E_1^{(1,0)}\equiv\left\langle \Phi^{(0)}_1, \mathcal{L}^{(1)}\Phi^{(0)}_1 \right\rangle = {\bf K}^T\ A_{0,0}^{(1)}\ {\bf K}\ +\
 2{\bf K}^T\ A_{0,-1}^{(1)}\ R\ {\bf K}$.

\end{proof}

Note that if ${\bf K}^T\ A_{0,-1}^{(1)}\ R\ {\bf K}\ne0$, then
for all small non-zero $\varepsilon$, in a neighborhood of $E^{(0)}=|{\bf K}|^2$ are two distinct $L^2_{{\bf K}}-$ eigenvalues separated by $\mathcal{O}(\varepsilon)$: a simple $L^2_{{\bf K},1}-$ eigenvalue, $\tilde{E}^\varepsilon_1=E_1^\varepsilon$ and doubly-degenerate
 $L^2_{{\bf K},\tau}\oplus L^2_{{\bf K},\bar\tau}-$ eigenvalue,
  $E^\varepsilon_D=E^\varepsilon_\tau=E^\varepsilon_{\bar\tau}$.

 If $\varepsilon{\bf K}^T\ A_{0,-1}^{(1)}\ R\ {\bf K}>0$, then $E^\varepsilon_D<\tilde{E}^\varepsilon$; Dirac points occur at the intersection of the first and second dispersion surfaces.
 If $\varepsilon{\bf K}^T\ A_{0,-1}^{(1)}\ R\ {\bf K}<0$, then $E^\varepsilon_D>\tilde{E}^\varepsilon$; Dirac points occur at the intersection of the second and third dispersion surfaces. This verifies the assertion in part 5 of Theorem \ref{eps_0_prop} and the proof of Theorem \ref{eps_0_prop} is now complete.

\subsection{Dirac points for arbitrary contrast honeycomb structures}\label{generic_eps}

Theorem \ref{low-dp} studies Dirac point for low contrast honeycomb structures by studying
 $\mathcal{L}^{(\varepsilon)}=-\nabla\cdot\left(\ I\ +\ \varepsilon\ A^{(1)}({\bf x})\ \right)\nabla$ for $\varepsilon\ne0$ and small.
In this section we discuss an extension of these results on Dirac points to arbitrary contrast structures.  \medskip

\noindent {\bf Assumption:} For all $\varepsilon>0$, $A^{(\varepsilon)}({\bf x})\equiv I\ +\ \varepsilon\ A^{(1)}({\bf x}) $ defines a honeycomb structured medium and therefore  $\mathcal{L}^{(\varepsilon)}\equiv -\nabla\cdot\ A^{(\varepsilon)}({\bf x})\ \nabla$ is  self-adjoint, positive definite and uniformly elliptic on $\mathbb{R}^2$.
\medskip

We state the result and then briefly discuss the strategy of proof, implemented in full detail
in the context of Schr\"odinger operators by C. L. Fefferman and one of the authors; see \cites{FW:12} and
Appendix D of \cites{FLW-MAMS:17}.  We do not present the detailed implementation of this strategy in this work.
\medskip

\noindent {\bf Claim:}\ Let $\varepsilon^0$ be as in Theorem \ref{low-dp}. There exists a discrete set
$\tilde{\mathcal{C}}\subset \mathbb{R}\setminus(0,\varepsilon^0)$, such that if $\varepsilon\notin\tilde{\mathcal{C}}$, then the conditions of
 Proposition \ref{prop_conical} hold for some quasi-momentum energy pair $({\bf K}_\star,E_D)$, where ${\bf K}_\star$ is any vertex of the Brillouin zone. It follows that $({\bf K}_\star,E_D)$ is a Dirac point.  \medskip
 \medskip

\noindent {\bf N.B.}\ For small $0<|\varepsilon|<\varepsilon^0$, Theorem
 \ref{low-dp} ensures that these Dirac points occur at the intersections between the first and second, or second and third dispersion surfaces.
  For general $\varepsilon\notin\tilde{\mathcal{C}}$, possibly large, we make no assertions on which dispersion surfaces intersect at Dirac points.
\medskip

The strategy is based on an analytical characterization of the $L^2_{{\bf K}_\star,\sigma}-$ eigenvalues of $\mathcal{L}^{(\varepsilon)}$. By assumptions on $A^{(\varepsilon)}$ (strong ellipticity and symmetries), $T(\varepsilon) \equiv (\mathcal{L}^{(\varepsilon)})^{-1}$
is defined as a bounded operator from $L^2_{{\bf K}_\star,\sigma}$ to $H^2_{{\bf K}_\star,\sigma}$. Hence, the $L^2_{{\bf K}_\star,\sigma}-$ eigenvalue problem \eqref{eps_evp} may formulated equivalently as:
\begin{equation}
(I-E\ T(\varepsilon))\Phi=0, \quad \Phi\in L^2_{{\bf K}_\star, \sigma}.
\label{LipSch}
\end{equation}
We require a global analytical criteria for a complex number, $E$, to be an eigenvalue \eqref{LipSch}.
Although  $T(\varepsilon)$ is compact, it is not trace class, and therefore its determinant is not defined. However, $T(\varepsilon)$ is a Hilbert-Schmidt operator and we proceed by working with its modified determinant $\det_2(I-ET(\varepsilon))$ \cites{Newton:72,Simon:05,GLZ:08}.

\begin{theorem}\label{arbitrary-dp}
Let $\sigma$ take on the values $1, \;\tau$ or $\bar{\tau}$.
\begin{enumerate}
\item $\varepsilon \mapsto T(\varepsilon)$ is an analytic mapping from $\mathbb C$ to the space of Hilbert-Schmidt operators on $L^2_{{\bf K},\sigma}$.
\item For $T(\varepsilon)$, considered as a mapping on $L^2_{{\bf K}_\star,\sigma}$, define
\begin{equation}
\mathcal E_\sigma(E,\varepsilon)={\det}_2(I-ET(\varepsilon)).
\label{det2}\end{equation}
The mapping $(E,\varepsilon)\mapsto \mathcal E_\sigma(E,\varepsilon)$, is analytic.
\item For $\varepsilon$ real, $E$ is an $L^2_{{\bf K}_\star,\sigma}-$ eigenvalue of geometric multiplicity $m$ if and only if $E$ is a root of $\mathcal E_\sigma(E,\varepsilon)=0$ of multiplicity $m$.
\end{enumerate}
\end{theorem}

The complex function theory strategy of \cites{FW:12} (see also Appendix D of \cites{FLW-MAMS:17}) can be used to establish that for all $\varepsilon$ which fall outside of a discrete subset of $\tilde{\mathcal{C}}\subset\mathbb{R}\setminus(0,\varepsilon^0)$ that there exists $E_D^\varepsilon\in\mathbb{R}$ such that
(a)\ $E_D^\varepsilon\in\mathbb{R}$ is a simple zero of  $\mathcal{E}_\tau(E_D^\varepsilon,\varepsilon)$ and $\mathcal{E}_{\overline\tau}(E_D^\varepsilon,\varepsilon)$,\ (b)\ $\mathcal{E}_1(E_D^\varepsilon,\varepsilon)\ne0$, and (c) $\upsilon_{_F}(\varepsilon)\ne0$. Here, $\mathcal{E}_\sigma(E,\varepsilon)$,
 for $\sigma=1,\tau,\overline\tau$, is defined in \eqref{det2}. Therefore, by Proposition \ref{prop_conical}, for all such $\varepsilon$ there exist of Dirac points at the vertices of the Brillouin zone, $\mathcal{B}$.

\section{Dirac points under perturbation - instability and persistence}\label{dirac_persistence}

In Section \ref{dirac-pts} we studied the existence of Dirac points of honeycomb structures ${\mathcal{L}^A}=-\nabla\cdot A\nabla$ in the setting where  ${\mathcal{L}^A}$ commutes with $\mathcal{PC}$ and with $\mathcal{R}$. Theorem \ref{prop_conical} gives criteria for the exists of Dirac points at vertices of the Brillouin zone $\mathcal{B}_h$. Theorem \ref{low-dp} in Section \ref{sec:dp-low}, proved using  Theorem \ref{prop_conical}, studies of Dirac points for low-contrast media. In Section \ref{generic_eps} discuss Dirac points for generic honeycomb structures, without any assumptions on contrast, appealing to the continuation argument developed in \cites{FW:12} for the case of Schr\"odinger operators with honeycomb lattice potential, $-\Delta+V({\bf x})$, where $V({\bf x})$ is $\Lambda_h-$ periodic, and $\mathcal{PC}-$ and $\mathcal{R}-$ invariant.

 In this section we discuss how the locally conical structure of dispersion surfaces near Dirac points deform when $\mathcal{PC}-$ invariance in broken. For results on the case of Schr\"odinger operators with a honeycomb potentials, see \cite{FW:12,berkolaiko-comech:15}.
 %

 Starting with ${\mathcal{L}^A}$, we introduce perturbed operators of the form
\begin{equation}
\label{L_delta}
\mathcal{L}^{(\delta)}=-\nabla\cdot\left(A + \delta B \right)\nabla=-\nabla\cdot A\nabla-\delta \nabla\cdot B\nabla \equiv \mathcal{L}^A + \delta \mathcal{L}^B.
\end{equation}
Here,  $B=B({\bf x})$ is a smooth $2\times2$ Hermitian matrix function, which is $\Lambda_h-$periodic.  The parameter, $\delta$, is the strength of the perturbation, and is taken to be a real number and sufficiently small.

The perturbed operator, $\mathcal{L}^{(\delta)}$, breaks $\mathcal{PC}-$ invariance if $[\mathcal{PC},\mathcal{L}^{(\delta)}]\ne0$, or equivalently if
$[\mathcal{PC},{\mathcal{L}^B}]\ne0$. We shall, in particular, assume
 that ${\mathcal{L}^B}$ is $\mathcal{PC}-$ anti-symmetric:
 \begin{equation}
 \mathcal{PC}{\mathcal{L}^B}\ =\ -{\mathcal{L}^B}\mathcal{PC}. \label{pc-anti}
 \end{equation}

In Appendix \ref{maxwell},
 %
 we show that this class of $\mathcal{PC}-$ invariance breaking perturbation encompasses examples of magneto-optic and bi-anisotropic media introduced discussed in the  introduction.
 %
 %

The perturbed eigenvalue problem in $L^2(\mathbb{R}^2/\Lambda_h)$ is given by
\begin{equation}\label{FB-delta}
\mathcal{L}^{(\delta)}({\bf k})\phi^{\delta}({\bf x};{\bf k})= E^{\delta}({\bf k})\phi^{\delta}({\bf x};{\bf k}), \quad \phi^{\delta}({\bf x}+{\bf v};{\bf k})=\phi^{\delta}({\bf x};{\bf k})\ \  \forall\ {\bf v}\in\Lambda_h,
\end{equation}
where
\begin{equation}
\label{L_delta_k}
\mathcal{L}^{(\delta)}({\bf k})= \mathcal{L}^A({\bf k}) + \delta \mathcal{L}^B({\bf k}) \equiv -(\nabla+i{\bf k})\cdot A(\nabla+i{\bf k})-\delta(\nabla+i{\bf k})\cdot B(\nabla+i{\bf k}).
\end{equation}
To study the deformation of dispersion surfaces near Dirac points under perturbation we study the perturbed Floquet-Bloch eigenvalue problem \eqref{FB-delta} for $\delta$ small and  ${\bs \kappa}\equiv {\bf k}-{\bf K}_\star$ small,
 where ${\bf K}_\star$ is a vertex of $\mathcal{B}_h$.

By Theorem \ref{low-dp}, the operator ${\mathcal{L}^A}({\bf K}_\star)$ has a doubly degenerate $L^2(\mathbb{R}^2/\Lambda_h)$- eigenvalue, which we denote $E_{D}$ (independent of the particular vertex ${\bf K}_\star$) with corresponding two dimensional eigenspace: $\textrm{span}
\{\phi^{{\bf K}_\star}_1({\bf x}),\phi^{{\bf K}_\star}_2({\bf x})\}$.
Expansion of $\mathcal{L}^{(\delta)}({\bf k})$ about ${\bf k}={\bf K}_\star$ gives:
\begin{align*}
\mathcal{L}^{(\delta)}({\bf K}_\star+{\bs \kappa})=\mathcal{L}^A({\bf K}_\star)+\delta \mathcal{L}^B({\bf K}_\star)+{\bs \kappa}\cdot \mathscr{A}({\bf K}_\star)+\mathscr{R}_2({\bs \kappa},\delta).
\end{align*}
Here $\mathscr{A}({\bf K}_\star)$ is defined in \eqref{Ak_def} and $\mathscr{R}_2({\bs \kappa},\delta)$ includes all terms of order $\mathcal{O}(|{\bs \kappa}|^2+\delta|{\bs \kappa}|)$:
\begin{equation*}
 \mathscr{R}_2({\bs \kappa},\delta) \equiv {\bs \kappa}^T A{\bs \kappa}-\delta\left(i{\bs \kappa}\cdot B(\nabla+i{\bf K}_\star)+i(\nabla+i{\bf K}_\star)\cdot(B{\bs \kappa})-{\bs \kappa}^TB{\bs \kappa}\right).
\end{equation*}
Let
\begin{equation*}
\phi^{\delta}({\bf x};{\bf k})=\phi_{{\bf K}_\star}^{(0)}({\bf x})+\phi_{{\bf K}_\star}^{(1)}({\bf x};{\bs \kappa}),  \quad  E^{\delta}({\bf k})= E_D+ E_{{\bf K}_\star}^{(1)}({\bs \kappa}),
\end{equation*}
where $\phi_{{\bf K}_\star}^{(0)}\in \textrm{span}
\{\phi^{{\bf K}_\star}_1,\phi^{{\bf K}_\star}_2\}$
 and $\left\langle\phi^{{\bf K}_\star}_j,\phi_{{\bf K}_\star}^{(1)}(\cdot;{\bs \kappa})\right\rangle=0,\ j=1,2$.
We follow a Lyapunov-Schmidt reduction strategy, analogous to that used in the proof of Theorem \ref{prop_conical}. Since $\phi_{{\bf K}_\star}^{(0)}({\bf x})$ is in the nullspace of ${\mathcal{L}^A}({\bf K}_\star)-E_D$, we have
\begin{equation}
\phi_{{\bf K}_\star}^{(0)}({\bf x})=\alpha_1^{{\bf K}_\star}\phi_1^{{\bf K}_\star}({\bf x})+\alpha_2^{{\bf K}_\star}\phi_2^{{\bf K}_\star},
\end{equation}
where $\alpha_1^{{\bf K}_\star},\ \alpha_2^{{\bf K}_\star}$ are complex constants to be determined. Recall also that
 $\phi_j^{{\bf K}_\star}({\bf x})=e^{-i{\bf K}_\star\cdot{\bf x}}\Phi_j^{{\bf K}_\star}({\bf x})$, where $\Phi_1^{{\bf K}_\star}\in L^2_{{\bf K}_\star,\tau}$ and $\Phi_2^{{\bf K}_\star}\in L^2_{{\bf K}_\star,\bar\tau}$; see Theorem \ref{prop_conical}.

  Calculations which are analogous to those  in the proof of Theorem \ref{prop_conical}, lead to a system of homogeneous linear equations for $\alpha_1^{{\bf K}_\star}$ and $\alpha_2^{{\bf K}_\star}$:

\begin{equation}
\label{PT_K_eqn}
\left( E_{{\bf K}_\star}^{(1)}I-\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bs \kappa})-\delta\mathcal M^{{\bf K}_\star}_{{\mathcal{L}^B}}-\mathcal M^{{\bf K}_\star}_{\mathscr{R}_2}( E^{(1)},{\bs \kappa},\delta) \right)\begin{pmatrix}\alpha_1^{{\bf K}_\star}\\ \alpha_2^{{\bf K}_\star}\end{pmatrix}=0,
\end{equation}
where, by Proposition \ref{prop_g} (all inner products over $L^2_{{\bf K}_\star}$),
\begin{equation}
\label{J_defn}
\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bs \kappa}) =
\upsilon_{_F}\ \begin{pmatrix}0 & \kappa^{(1)}+i\kappa^{(2)}\\
\kappa^{(1)}-i\kappa^{(2)}&0\end{pmatrix},
\end{equation}
\begin{equation}
\label{A_defn}
\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B} =
\begin{pmatrix}
\inner{\Phi^{{\bf K}_\star}_1, \mathcal{L}^B \Phi^{{\bf K}_\star}_1} & \inner{\Phi^{{\bf K}_\star}_1, \mathcal{L}^B \Phi^{{\bf K}_\star}_2}\\
\inner{\Phi^{{\bf K}_\star}_2, \mathcal{L}^B \Phi^{{\bf K}_\star}_1}&\inner{\Phi^{{\bf K}_\star}_2, \mathcal{L}^B \Phi^{{\bf K}_\star}_2}
\end{pmatrix} ;
\end{equation}
and $|\mathcal{M}^{{\bf K}_\star}_{\mathcal{R}}( E^{(1)},{\bs \kappa},\delta)|
\lesssim|{\bs \kappa}|^2+|\delta|\ |{\bs \kappa}|$.

Up to this point we have not used particular properties of the perturbation, ${\mathcal{L}^B}$. In the following subsections, we discuss conditions on $B({\bf x})$ under which Dirac points are unstable (conical behavior perturbs to locally smooth and ``gapped'' dispersion surfaces) and conditions under which the local conical structure of  Dirac points persists. The key, as we  see below, is to determine the character of the perturbation matrix
 $\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}$. We note that $\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}$ is Hermitian since ${\mathcal{L}^B}$ is self-adjoint.

\subsection{Instability of Dirac points for a class of $\mathcal{PC}-$ breaking perturbations}\label{unstable-dirac}

\begin{proposition}\label{MLB}

  If $\mathcal{L}^B$ is $\mathcal{PC}-$ anti-symmetric, {\it i.e.} $\mathcal{PC} \mathcal{L}^B= -\mathcal{L}^B \mathcal{PC}$, then
\begin{equation}
\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}=\begin{pmatrix}{\vartheta_{\sharp}}^{{\bf K}_\star} &0\\0&-{\vartheta_{\sharp}}^{{\bf K}_\star} \end{pmatrix},
\label{A_defn1}\end{equation}
where ${\vartheta_{\sharp}}^{{\bf K}_\star}$ is given by
\begin{equation}
{\vartheta_{\sharp}}^{{\bf K}_\star}\equiv \inner{\Phi^{{\bf K}_\star}_1,\mathcal{L}^B \Phi^{{\bf K}_\star}_1},
\label{thta_shp}
\end{equation}
and is real, by self-adjointness of ${\mathcal{L}^B}$.
\end{proposition}

{\bf Assumption:} Throughout we shall assume that
\begin{equation}
{\vartheta_{\sharp}}^{{\bf K}_\star}\ne0.
\label{thta_ne0}
\end{equation}
\bigskip

\begin{proof}[Proof of Proposition \ref{MLB}]
Suppose $\mathcal{PC}\mathcal{L}^B=-\mathcal{L}^B\mathcal{PC}$. In the calculations below, we frequently use the following equality: for any $f,g\in L^2_{{\bf K}_\star}$,
\begin{equation*}
\inner{\mathcal{PC}f,\mathcal{PC}g}=\inner{\mathcal{C}f,\mathcal{C}g}=\overline{\inner{f,g}}=\inner{g,f}.
\end{equation*}

Consider first  the off-diagonal elements of
 $\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}$, we have
 %
 \begin{align*}
\inner{\Phi^{{\bf K}_\star}_1,\mathcal{L}^B \Phi^{{\bf K}_\star}_2}&=\inner{\mathcal{PC}\Phi^{{\bf K}_\star}_2, \mathcal{L}^B \mathcal{PC}\Phi^{{\bf K}_\star}_1} =-\inner{\mathcal{PC}\Phi^{{\bf K}_\star}_2,\mathcal{PC} \mathcal{L}^B\Phi^{{\bf K}_\star}_1}\nonumber\\
&=
-\inner{\mathcal{L}^B\Phi^{{\bf K}_\star}_1,\Phi^{{\bf K}_\star}_2}=-\inner{\Phi^{{\bf K}_\star}_1, \mathcal{L}^B\Phi^{{\bf K}_\star}_2},
\end{align*}
where we have used the self-adjointness of $\mathcal{L}^B$ at the last equality. Then it follows that $\inner{\Phi^{{\bf K}_\star}_1, \mathcal{L}^B \Phi^{{\bf K}_\star}_2}=0$.

We now turn to  the diagonal entries of $\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}$.  We have
\begin{align}
\label{theta_defn1}
\inner{\Phi^{{\bf K}_\star}_1, \mathcal{L}^B\Phi^{{\bf K}_\star}_1}&=\inner{\mathcal{PC} \Phi^{{\bf K}_\star}_2,\mathcal{L}^B \mathcal{PC} \Phi^{{\bf K}_\star}_2}\nonumber\\
&=-\inner{\mathcal{PC} \Phi^{{\bf K}_\star}_2,  \mathcal{PC} \mathcal{L}^B \Phi^{{\bf K}_\star}_2}=-\inner{\mathcal{L}^B\Phi^{{\bf K}_\star}_2, \Phi^{{\bf K}_\star}_2}={\vartheta_{\sharp}}^{{\bf K}_\star}.
\end{align}
This completes the proof of Proposition \ref{MLB}.
\end{proof}

It follows from Proposition \ref{MLB} that if $\mathcal{PC}{\mathcal{L}^B}=-{\mathcal{L}^B} \mathcal{PC}$, then

\begin{equation}
\label{local_K_matrix}
\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bs \kappa})+\delta\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}=
\begin{pmatrix}\delta{\vartheta_{\sharp}}^{{\bf K}_\star} &\upsilon_{_F}(\kappa^{(1)}+i\kappa^{(2)})\\
\upsilon_F(\kappa^{(1)}-i\kappa^{(2)})&-\delta{\vartheta_{\sharp}}^{{\bf K}_\star} \end{pmatrix}.
\end{equation}
 By \eqref{PT_K_eqn} and \eqref{local_K_matrix}, the energy $E=E_D+E_{{\bf K}_\star}^{(1)}$ is a $L^2_{{\bf K}_\star+{\bs \kappa}}-$ eigenvalue of the perturbed Floquet-Bloch eigenvalue problem for $\mathcal{L}^{(\delta)}$ if and only if
 \begin{equation}
 \label{det_rel}
 \det ( E_{{\bf K}_\star}^{(1)}I-\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bs \kappa})-\delta\mathcal M^{{\bf K}_\star}_{{\mathcal{L}^B}}-\mathcal M^{{\bf K}_\star}_{\mathscr{R}_2}( E^{(1)},{\bs \kappa},\delta))  = 0.
 \end{equation}
 The eigenvalue condition \eqref{det_rel} is of the form
\begin{equation}
\label{local_K_energy}
 ( E_{{\bf K}_\star}^{(1)})^2-\delta^2({\vartheta_{\sharp}}^{{\bf K}_\star})^2 -\upsilon_{_F}^2|{\bs \kappa}|^2
 + g_{120} + g_{111} + g_{012} + g_{021} + g_{003} = 0,
\end{equation}
where $g_{rsl}=g_{rsl}( E_{{\bf K}_\star}^{(1)},\delta,{\bs \kappa} )$ are smooth and satisfy the bound
$g_{rsl}( E_{{\bf K}_\star}^{(1)},\delta,{\bs \kappa} ) \leq C |E_{{\bf K}_\star}^{(1)}|^r|\delta|^s|{\bs \kappa}|^l$.
After some manipulations, we may apply the implicit function theorem to obtain from \eqref{local_K_energy} that
\begin{equation}
\label{PT_antisym_local_energy}
 E_{{\bf K}_\star}^{(1)}=\pm\sqrt{\upsilon_{_F}^2|{\bs \kappa}|^2+\delta^2({\vartheta_{\sharp}}^{{\bf K}_\star})^2}(1+e^{{\bf K}_\star}_\pm({\bs \kappa},\delta)),
\end{equation}
where
$e^{{\bf K}_\star}_\pm({\bs \kappa},\delta)=\mathcal{O}(\delta+{\bs \kappa})$.

Therefore, provided ${\vartheta_{\sharp}}^{{\bf K}_\star}\neq0$, the Dirac point $({\bf K}_\star, E_D)$ does not persist in the presence of $\mathcal{PC}$- anti-symmetric perturbations. That is, if $\mathcal{PC}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{PC}$ and ${\vartheta_{\sharp}}^{{\bf K}_\star}\neq0$, then for all $\delta$ sufficiently small, the dispersion surfaces are locally smooth and a \underline{local} spectral gap up for quasi-momenta near ${\bf K}_\star$.  We summarize the above discussion in the following

\begin{theorem}\label{PT_thm}
Consider the operator
$\mathcal{L}^{(\delta)}={\mathcal{L}^A}+\delta{\mathcal{L}^B}=-\nabla\cdot(A+\delta B)\nabla$ defined in \eqref{L_delta}. For the unperturbed operator we assume that $A$ is smooth, Hermitian, $\Lambda_h-$ periodic, and that $[\mathcal{PC},{\mathcal{L}^A}]=0$ and $[ \mathcal{R},{\mathcal{L}^A}]=0$.
For the perturbed operator we assume that
$B({\bf x})$ is smooth, Hermitian, $\Lambda_h-$ periodic, and that
$\mathcal{PC}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{PC}$ ($\mathcal{PC}-$ anti-symmetry).

 Let ${\bf K}_\star$ be a point of ${\bf K}$ or ${\bf K}'$ type, and let $({\bf K}_\star,E_D)$ denote a Dirac point of ${\mathcal{L}^A}$ (Theorem \ref{low-dp}).
Assume ${\vartheta_{\sharp}}^{{\bf K}_\star} \equiv \inner{\Phi^{{\bf K}_\star}_1,{\mathcal{L}^B} \Phi^{{\bf K}_\star}_1}\neq0$. Then,
there exists a $\delta_0>0$ such that for $0<\delta<\delta_0$,
 the Dirac point $({\bf K}_\star,E_D)$ does not persist, the perturbed dispersion surfaces are locally smooth and a local spectral gap opens up in a neighborhood of $({\bf k},E)=({\bf K}_\star,E_D)$.
\end{theorem}

\begin{example}\label{ex-pcanti} Here we give two typical examples of a family of matrices, $B({\bf x})$, satisfying the hypotheses of Theorem \ref{PT_thm}. Let $\mu_1({\bf x})$ be real-valued, $\mu_1({\bf x})=\mu_1(-{\bf x})$, and $\mu_1(R^*{\bf x})=\mu_1({\bf x})$. Define
\begin{equation}
B({\bf x})\equiv\sigma_2\mu_1({\bf x}).\ \
\label{B-example}\end{equation}
Then, $[\mathcal{R},{\mathcal{L}^B}]=0$ and $\mathcal{PC}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{PC}$, since
$
 \mathcal{C}[B({\bf x})\cdot]=-\sigma_2\mu({\bf x})\mathcal{C}\cdot=-B({\bf x})\mathcal{C}\cdot$,
 while $
 \mathcal{P}[B({\bf x})\cdot]=\sigma_2\mu(-{\bf x})\mathcal{P}\cdot
 = \sigma_2\mu({\bf x})\mathcal{P}\cdot=B({\bf x})\mathcal{P}\cdot$.

 The other example is
 \begin{equation*}
B({\bf x})\equiv\mu_2({\bf x}) I_{2\times 2}.\ \
\label{B-example2}\end{equation*}
 where $\mu_2({\bf x})$ is a real-valued function satisfying $\mu_2(-{\bf x})=-\mu_2({\bf x})$, and $\mu_2(R^*{\bf x})=\mu_2({\bf x})$.
Then, in this case $[\mathcal{R},{\mathcal{L}^B}]=0$ and $\mathcal{PC}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{PC}$, since
$
 \mathcal{C}[B({\bf x})\cdot]=B({\bf x})\mathcal{C}\cdot$,
 while $
 \mathcal{P}[B({\bf x})\cdot]=\mu_2(-{\bf x})\mathcal{P}\cdot
 = -B({\bf x})\mathcal{P}\cdot$.

\end{example}

\begin{remark}\label{thetasharp_bK_bKp}
In our discussion of edge states in Section \ref{edge_states}, we shall consider two important sub-cases, where $\mathcal{L}^B$ is $\mathcal{PC}$ anti-symmetric. Case(a): $[\mathcal{P},{\mathcal{L}^B}]=0$ and $\mathcal{C}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{C}$, as in Example \ref{ex-pcanti}, and Case (b) $\mathcal{P}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{P}$, but $[\mathcal{C},{\mathcal{L}^B}]=0$. We shall see that
for case (a): ${\vartheta_{\sharp}}^{\bf K'}=+{\vartheta_{\sharp}}^{\bf K}$ and for case (b): ${\vartheta_{\sharp}}^{\bf K'}=-{\vartheta_{\sharp}}^{\bf K}$.  Implications  for the directionality of edge states in each of cases are explored in Section \ref{edge_states}.
\end{remark}

\subsection{Remarks on the persistence of Dirac points under $\mathcal{P}\circ\mathcal{C}$ symmetry preserving perturbations}
\label{persistence_dirac}

\begin{proposition}\label{persist-prop}
 Assume $B$ is $\Lambda_h-$ periodic
 and $\mathcal{L}^B$ is $\mathcal{PC}-$-symmetric, {\it i.e.} $[\mathcal{PC}, \mathcal{L}^B]=0$.
  Then,
  \begin{equation}
\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}=
\begin{pmatrix}{\vartheta_{\sharp}}^{{\bf K}_\star} & \varrho_\sharp^{{\bf K}_\star}\\ \overline{\varrho_\sharp^{{\bf K}_\star}} &{\vartheta_{\sharp}}^{{\bf K}_\star} \end{pmatrix},
\end{equation}
 where ${\vartheta_{\sharp}}^{{\bf K}_\star}\in\mathbb{R}$ is given in \eqref{thta_shp} and $\varrho_\sharp^{{\bf K}_\star}\equiv \inner{\Phi^{{\bf K}_\star}_1, {\mathcal{L}^B}\Phi^{{\bf K}_\star}_2}$.
\end{proposition}

\begin{proof}[Proof of of Proposition \ref{persist-prop}]
If $\mathcal{PC}\mathcal{L}^B=\mathcal{L}^B\mathcal{PC}$, then
\begin{align*}
{\vartheta_{\sharp}}^{{\bf K}_\star} &\equiv \inner{\Phi^{{\bf K}_\star}_1,\mathcal{L}^B \Phi^{{\bf K}_\star}_1}=\inner{\mathcal{PC} \Phi^{{\bf K}_\star}_2, \mathcal{L}^B \mathcal{PC} \Phi^{{\bf K}_\star}_2} \\
&=\inner{\mathcal{PC} \Phi^{{\bf K}_\star}_2,\mathcal{PC} \mathcal{L}^B \Phi^{{\bf K}_\star}_2} = {\inner{\Phi^{{\bf K}_\star}_2, \mathcal{L}^B \Phi^{{\bf K}_\star}_2}},
\end{align*}
and
\begin{equation}
\varrho_\sharp^{{\bf K}_\star} \equiv \inner{\Phi^{{\bf K}_\star}_1, {\mathcal{L}^B}\Phi^{{\bf K}_\star}_2}=\inner{{\mathcal{L}^B}\Phi^{{\bf K}_\star}_1, \Phi^{{\bf K}_\star}_2}=\overline{\inner{\Phi^{{\bf K}_\star}_2, {\mathcal{L}^B}\Phi^{{\bf K}_\star}_1}}.
\end{equation}
This completes the proof of Proposition \ref{persist-prop}.
\end{proof}

Therefore, if $\mathcal{PC}{\mathcal{L}^B}={\mathcal{L}^B} \mathcal{PC}$, we have
\begin{equation}
\label{pc_M_matrices}
\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bs \kappa})+\delta\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}=
\begin{pmatrix}\delta{\vartheta_{\sharp}}^{{\bf K}_\star} &\upsilon_F(\kappa^{(1)}+i\kappa^{(2)}) + \delta\varrho_\sharp^{{\bf K}_\star}\\
\upsilon_F(\kappa^{(1)}-i\kappa^{(2)})+ \delta\overline{\varrho_\sharp^{{\bf K}_\star}}&\delta{\vartheta_{\sharp}}^{{\bf K}_\star} \end{pmatrix}.
\end{equation}
By \eqref{PT_K_eqn} and \eqref{pc_M_matrices}, the energy $E=E_D+E_{{\bf K}_\star}^{(1)}$ is a $L^2_{{\bf K}_\star+{\bs \kappa}}-$ eigenvalue of the perturbed Floquet-Bloch eigenvalue problem for $\mathcal{L}^{(\delta)}$ if and only if
 \begin{equation*}
 \det ( E_{{\bf K}_\star}^{(1)}I-\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bs \kappa})-\delta\mathcal M^{{\bf K}_\star}_{{\mathcal{L}^B}}-\mathcal M^{{\bf K}_\star}_{\mathscr{R}_2}( E^{(1)},{\bs \kappa},\delta))  = 0.
 \end{equation*}
 Following the same procedure as in Section \ref{unstable-dirac} but for \eqref{pc_M_matrices}, we find that
 that $E_{{\bf K}_\star}^{(1)}$ is given by
\begin{equation}
\label{PT_sym_local_energy}
 E_{{\bf K}_\star}^{(1)}= \left( \delta{\vartheta_{\sharp}}^{{\bf K}_\star} \pm \abs{\upsilon_F(\kappa^{(1)}+i\kappa^{(2)})+\delta\varrho_\sharp^{{\bf K}_\star}} \right) ( 1 +e^{{\bf K}_\star}_\pm({\bs \kappa},\delta) ) ,
\end{equation}
where $e^{{\bf K}_\star}_\pm({\bs \kappa},\delta)=\mathcal{O}(\delta+|{\bs \kappa}|)$.

The local energy expansion \eqref{PT_sym_local_energy} demonstrates, to finite order in $\delta$, that the Dirac point $({\bf K}_\star,E_D)$ is protected by $\mathcal{PC}-$ symmetry. That is, if $[\mathcal{PC},{\mathcal{L}^A}]=0$, the Dirac point degeneracy shifts to $({\bf K}_\star+\delta\bb b,E_D+\delta{\vartheta_{\sharp}}^{{\bf K}_\star})$, where $\bb b=(b^{(1)},b^{(2)})$ satisfies $\upsilon_F(b^{(1)}+ib^{(2)})+\varrho_\sharp^{{\bf K}_\star}=0$
A rigorous proof can be implemented along the lines of Theorem 9.1  in \cites{FW:12}.

\section{Edge States}\label{edge_states}

Edge states are time-harmonic solutions of the wave equation (obtained via reduction of the 2D Maxwell Equations -- see Remark \ref{wavepkts} and Appendix \ref{maxwell}), which are bounded and oscillatory in the direction parallel to an extended line-defect, and spatially localized transverse to it.

In Sections \ref{rational_edges} we review the mathematical framework of \cites{FLW-AnnalsPDE:16,FLW-2DM:15} for describing ``rational'' edges in a honeycomb structure. We then introduce our model of a honeycomb structure perturbed by an edge in Section \ref{honeycomb_model}. We formally derive edge states in Section \ref{multiple_scales}, and sketch the rigorous justifications, along the lines of \cites{FLW-AnnalsPDE:16}, of our formal results in Section \ref{rigor}. Sections \ref{edge_state_kpar} and \ref{direction} investigate the direction of wavepacket propagation along the edge. {\it For simplicity, we restrict our attention to the case of real-valued  honeycomb structured media, $A({\bf x})$; see Remark \ref{dp-nonreal}.}

\subsection{Rational edges}\label{rational_edges}

Recall, from Section \ref{triangle}, that the lattice and dual lattice are given by $\Lambda_h=\mathbb{Z}{\bf v}_1\oplus\mathbb{Z}{\bf v}_2$ and $\Lambda_h^*=\mathbb{Z}{\bf k}_1\oplus\mathbb{Z}{\bf k}_2$, respectively. An {\it edge} is a line  of the form $\mathbb{R} (a_1{\bf v}_1+a_2{\bf v}_2)$, where $(a_1,b_1)=1$, {\it i.e.} $a_1$ and $b_1$ are relatively prime integers. In Section \ref{honeycomb_model}, we shall introduce a family of operators, which interpolate between two distinct $\Lambda_h-$ periodic structures at ``$+\infty$'' and ``$-\infty$''.

We fix an edge by choosing a vector ${\bm{\mathfrak{v}}}_1=a_1{\bf v}_1+b_1{\bf v}_2$, where  $(a_1, b_1)=1$. Since  $a_1, b_1$ are relatively prime,  there exists a relatively prime pair of integers: $a_2,b_2$  such that $a_1b_2-a_2b_1=1$.  Set ${\bm{\mathfrak{v}}}_2 = a_2 {\bf v}_1 + b_2 {\bf v}_2$.
It follows that $\mathbb{Z}{\bm{\mathfrak{v}}}_1\oplus\mathbb{Z}{\bm{\mathfrak{v}}}_2=\mathbb{Z}{\bf v}_1\oplus\mathbb{Z}{\bf v}_2=\Lambda_h$.
Since $a_1b_2-a_2b_1=1$, we have dual lattice vectors ${\bm{\mathfrak{K}}}_1, {\bm{\mathfrak{K}}}_2\in\Lambda_h^*$, given by
\begin{equation}
\label{ktilde_1_2}
{\bm{\mathfrak{K}}}_1=b_2{\bf k}_1-a_2{\bf k}_2,\ \  {\bm{\mathfrak{K}}}_2=-b_1{\bf k}_1+a_1{\bf k}_2,
\end{equation}
  which satisfy
\begin{equation*} {\bm{\mathfrak{K}}}_\ell \cdot {\bm{\mathfrak{v}}}_{\ell'} = 2\pi \delta_{\ell, \ell'},\ \ 1\leq \ell, \ell' \leq 2.\label{ktilde-orthog}
\end{equation*}
Note that $\mathbb{Z}{\bm{\mathfrak{K}}}_1\oplus\mathbb{Z}{\bm{\mathfrak{K}}}_2=\mathbb{Z}{\bf k}_1\oplus\mathbb{Z}{\bf k}_2=\Lambda^*_h$.

Two typical examples of the edges are:
\begin{enumerate}
\item The zigzag edge: ${\bm{\mathfrak{v}}}_1={\bf v}_1,~{\bm{\mathfrak{v}}}_2={\bf v}_2$ and ${\bm{\mathfrak{K}}}_1={\bf k}_1,~{\bm{\mathfrak{K}}}_2={\bf k}_2$.
\item The armchair edge: ${\bm{\mathfrak{v}}}_1={\bf v}_1+{\bf v}_2,~{\bm{\mathfrak{v}}}_2={\bf v}_2$ and ${\bm{\mathfrak{K}}}_1={\bf k}_1,~{\bm{\mathfrak{K}}}_2= {\bf k}_2-{\bf k}_1$.
\end{enumerate}

\subsection{Model of a honeycomb structure with a rational edge}\label{honeycomb_model}

The edge through a honeycomb material is modeled through a perturbation of the bulk operator ${\mathcal{L}^A}$.
We have the following setup.
\begin{enumerate}
\item \emph{Unperturbed operator, ${\mathcal{L}^A}$}: Let $A({\bf x})$ be smooth, $\Lambda_h-$ periodic and such that $[\mathcal{PC},{\mathcal{L}^A}]=0$ and $[\mathcal{R},{\mathcal{L}^A}]=0$, {\it i.e.}, a honeycomb structured medium, as described in Section \ref{honey-media}.
\item \emph{Dirac points}: Let ${\bf K}_\star$ denote a quasi-momentum at the vertex of the Brillouin zone, $\mathcal{B}_h$ (a point of ${\bf K}$ or ${\bf K}'$ type), and assume that $({\bf K}_\star,E_D)$ is a Dirac point of the operator $\mathcal{L}^A=-\nabla\cdot A\nabla$; see Definition \ref{dirac_pt_defn}, {\it e.g.} the setting assumed in the hypotheses of Theorem \ref{prop_conical}.
\item \emph{Perturbed operator asymptotics far from the edge}: Let $B({\bf x})$ be a $\Lambda_h-$periodic, $2\times 2$ Hermitian matrix such that $\mathcal{L}^B \equiv -\nabla\cdot B\nabla$ is anti-$\mathcal{PC}-$ symmetric: $\mathcal{PC}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{PC}$; {\it i.e.}, part (1) of Theorem \ref{PT_thm} applies to the operator $-\nabla\cdot(A+\delta B)\nabla$. Furthermore, assume that
\begin{equation}
{\vartheta_{\sharp}}^{{\bf K}_\star} \equiv \inner{\Phi^{{\bf K}_\star}_1,{\mathcal{L}^B} \Phi^{{\bf K}_\star}_1}\neq0\ .
\label{thta_shp1}
\end{equation}
\end{enumerate}

Suppose we fix an edge, given by the line $\mathbb{R}{\bm{\mathfrak{v}}}_1$ or, alternatively, ${\bf x}\in\mathbb{R}^2$ such that ${\bm{\mathfrak{K}}}_2\cdot{\bf x}=0$.
Our operator is a smooth and slow interpolation, transverse to the edge,  between the operators
\begin{equation}
\label{asymptotic_ops}
\mathcal{L}_{\pm}^{(\delta)} \equiv -\nabla \cdot \left[ A({\bf x}) \pm \delta \eta_\infty B({\bf x})\right] \nabla,
\end{equation}
associated with periodic structures $A({\bf x})-\delta \eta_\infty B({\bf x})$ for ${\bm{\mathfrak{K}}}_2\cdot{\bf x}\to-\infty$ and $A({\bf x})+\delta \eta_\infty B({\bf x})$
 for ${\bm{\mathfrak{K}}}_2\cdot{\bf x}\to+\infty$. Here, $\eta_\infty$ is a positive constant.
This interpolation is  effected by a \emph{domain wall function}:

\begin{definition} \label{domain_wall_defn}
We call $\eta(\zeta)\in C^{\infty}(\mathbb{R})$ a  domain wall function if $\eta(\zeta)$ tends to $\pm\eta_\infty$ as $\zeta\to\pm\infty$. We take $\eta(0)=0$, and
without loss of generality, we assume $\eta_\infty>0$.
\end{definition}

\noindent Our model of a honeycomb structure with an edge is the domain-wall modulated operator:
\begin{equation}
 \label{dw_ham}
 \mathcal{L}_{\rm dw}^{(\delta)} \equiv -\nabla\cdot\left[A({\bf x}) +\delta\eta( \delta{\bm{\mathfrak{K}}}_2 \cdot {\bf x})B({\bf x})\right]  \nabla.
\end{equation}
The operator $\mathcal{L}_{\rm dw}^{(\delta)}$ breaks translation invariance with respect to arbitrary elements of the lattice, $\Lambda_h$, but  is invariant with respect to translation by ${\bm{\mathfrak{v}}}_1$, parallel to the edge (because ${\bm{\mathfrak{K}}}_2 \cdot{\bm{\mathfrak{v}}}_1=0$ in \eqref{dw_ham}). Associated with this translation invariance is a parallel quasi-momentum, which we denote by ${k_{\parallel}}$.

Edge states are solutions of the eigenvalue problem
\begin{align}
&\mathcal{L}_{_{\rm dw}}^{(\delta)}\Psi({\bf x};{k_{\parallel}}) = E({k_{\parallel}})\Psi({\bf x};{k_{\parallel}}), \label{dw_evp}\\
&\Psi({\bf x}+{\bm{\mathfrak{v}}}_1;{k_{\parallel}})=e^{i{k_{\parallel}}}\Psi({\bf x};{k_{\parallel}}),\qquad \textrm{(propagation parallel to the edge,  $\mathbb{R}{\bm{\mathfrak{v}}}_1$)},\label{pseudo-per}\\
&\Psi({\bf x};{k_{\parallel}}) \to 0\ \ {\rm as}\ \ |{\bf x}\cdot{\bm{\mathfrak{K}}}_2|\to\infty. \qquad  \textrm{(localization tranverse to the edge,  $\mathbb{R}{\bm{\mathfrak{v}}}_1$)}\label{localized} .
\end{align}
We refer to a solution pair $(E({k_{\parallel}}),\Psi({\bf x};{k_{\parallel}}))$ of \eqref{dw_evp}-\eqref{localized} as an edge state or edge mode.
We shall construct edge modes for ${k_{\parallel}}$ near ${\bf K}\cdot{\bm{\mathfrak{v}}}_1$ and near ${\bf K}'\cdot{\bm{\mathfrak{v}}}_1=-{\bf K}\cdot{\bm{\mathfrak{v}}}_1$.

\subsection{Multiple scales construction of the edge state}\label{multiple_scales}

Recall the notations
\begin{equation}
 \label{op_notation}
 \mathcal{L}^A \equiv -\nabla_{\bf x}\cdot A({\bf x})\nabla_{\bf x}, \qquad \mathcal{L}^B = -\nabla_{\bf x}\cdot B({\bf x})\nabla_{\bf x}, \qquad
 \mathscr{A} \equiv \frac{1}{i}A({\bf x}) \nabla_{\bf x} +\frac{1}{i}\nabla_{\bf x}\cdot A({\bf x}) ;
\end{equation}
see \eqref{L_def}, \eqref{L_delta} and \eqref{A_def}, respectively.

We first consider the eigenvalue problem \eqref{dw_evp}-\eqref{localized} for the choice of parallel quasi-momentum ${k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$.
Edge states for $|{k_{\parallel}}-{\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1|$ small can be constructed  perturbatively; see the discussion in Section \ref{edge_state_kpar}.

For $\delta\ll1$, we formally seek solutions of the eigenvalue problem \eqref{dw_evp}-\eqref{localized}, which depend on  fast (${\bf x}$) and slow/transverse ($\zeta=\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x}$) spatial scales:
\begin{align}
E^{\delta}&=E^{(0)}+\delta E^{(1)}+\cdots, \label{formal-E} \\
\Psi^{\delta}&=\Psi^{(0)}({\bf x},\zeta)+\delta \Psi^{(1)}({\bf x},\zeta)+\cdots,\ \ \zeta=\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x} . \label{formal-psi}
\end{align}
The pseudo-periodicity condition \eqref{pseudo-per}, with
 ${k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$, and decaying \eqref{localized} boundary conditions are encoded by requiring, for $j\ge0$:
\begin{align*}
&\Psi^{(j)}({\bf x}+{\bm{\mathfrak{v}}}_1,\cdot)=e^{i{{\bf K}_\star}\cdot{\bm{\mathfrak{v}}}_1}\Psi^{(j)}({\bf x},\cdot) \ \ \ \forall\ {\bm{\mathfrak{v}}}\in\Lambda_h , \quad \text{and} \\
&\zeta\to\Psi^{(j)}({\bf x},\zeta)\in L^2(\mathbb{R}_\zeta).
\end{align*}

We substitute the expansions \eqref{formal-E}-\eqref{formal-psi} into \eqref{dw_evp} and equate terms of equal order in $\delta^j$, $j\ge0$.
At order $\delta^0$ we have that $(E^{(0)},\Psi^{(0)})$ satisfies
\begin{equation}\label{Eq_delta0}
\begin{split}
(\mathcal{L}^A-E^{(0)})\Psi^{(0)}&=0,\\
\Psi^{(0)}({\bf x}+{\bm{\mathfrak{v}}}_1,\zeta)&=e^{i{{\bf K}_\star}\cdot{\bm{\mathfrak{v}}}_1}\Psi^{(0)}({\bf x},\zeta).
\end{split}
\end{equation}

We are interested in constructing solutions which are spectrally localized near a Dirac point $({{\bf K}_\star},E_D)$. Therefore, we solve \eqref{Eq_delta0} by taking
\begin{equation}
\label{delta_zero}
E^{(0)}= E_D, \quad \Psi^{(0)}=\alpha^{{\bf K}_\star}_1(\zeta)\Phi^{{\bf K}_\star}_1({\bf x})+\alpha^{{\bf K}_\star}_2(\zeta)\Phi^{{\bf K}_\star}_2({\bf x}),
\end{equation}
where the amplitudes, $\alpha^{{\bf K}_\star}_1(\zeta)$ and
$\alpha^{{\bf K}_\star}_2(\zeta)$ are to be determined.

Proceeding to order $\delta^1$, we find that $(E^{(1)},\Psi^{(1)})$ satisfies
\begin{equation}\label{Eq_delta1}
\begin{array}{l}(\mathcal{L}^A-E_D)\Psi^{(1)} =G_1^{(1)}({\bf x},\zeta;\Psi^{(0)})+G_2^{(1)}({\bf x},\zeta;\Psi^{(0)})+E^{(1)}\Psi^{(0)},\\[1.2 ex]
\Psi^{(1)}({\bf x}+{\bm{\mathfrak{v}}}_1,\zeta)=e^{i{{\bf K}_\star}\cdot{\bm{\mathfrak{v}}}_1}\Psi^{(1)}({\bf x},\zeta).
\end{array}
\end{equation}
where
\begin{equation}\label{G1}
G_1^{(1)}({\bf x},\zeta;\Psi^{(0)})=i(\partial_\zeta\alpha^{{\bf K}_\star}_1 {{\bm{\mathfrak{K}}}_2}\cdot \mathscr{A} \Phi^{{\bf K}_\star}_1+\partial_\zeta\alpha^{{\bf K}_\star}_2 {\bm{\mathfrak{K}}}_2 \cdot \mathscr{A} \Phi^{{\bf K}_\star}_2),
\end{equation}
and
\begin{equation}\label{G2}
G_2^{(1)}({\bf x},\zeta;\Psi^{(0)})=-\eta(\zeta)(\alpha^{{\bf K}_\star}_1 \mathcal{L}^B \Phi^{{\bf K}_\star}_1+\alpha^{{\bf K}_\star}_2 \mathcal{L}^B \Phi^{{\bf K}_\star}_2) .
\end{equation}

The pseudo-periodic boundary value problem \eqref{Eq_delta1}  is solvable if and only if its right hand side is $L^2_{{\bf K}_\star}-$ orthogonal to the nullspace of ${\mathcal{L}^A}-E_D$, which is spanned by $\Phi^{{\bf K}_\star}_j,\ j=1,2$. This yields the two solvability  relations:
\begin{equation}
\label{Solvability_d}
-E^{(1)}\alpha^{{\bf K}_\star}_j=\inner{\Phi^{{\bf K}_\star}_j,G_1^{(1)}({\bf x},\zeta;\Psi^{(0)})+G_2^{(1)}({\bf x},\zeta;\Psi^{(0)})}, \qquad j=1,2.
\end{equation}

Let ${\bm{\alpha}}^{{\bf K}_\star}(\zeta)=(\alpha^{{\bf K}_\star}_1(\zeta),\alpha^{{\bf K}_\star}_2(\zeta))^T$.
Substituting \eqref{G1} and \eqref{G2} into \eqref{Solvability_d}, gives
\begin{equation}
\left( \mathcal{D}^{{\bf K}_\star} - E^{(1)} \right) {\bm{\alpha}}^{{\bf K}_\star}(\zeta) = 0 , \quad {\bm{\alpha}}^{{\bf K}_\star}\in L^2(\mathbb{R})  , \label{dirac_eqn}
\end{equation}
where $\mathcal{D}^{{\bf K}_\star}$ denotes the 1D Dirac operator:
\begin{align}
 \label{dirac_op}
 \mathcal{D}^{{\bf K}_\star}
 &= -i\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bm{\mathfrak{K}}}_2) \partial_\zeta + \eta(\zeta){\vartheta_{\sharp}}^{{\bf K}_\star} \sigma_3 .
\end{align}
In simplifying \eqref{dirac_op}, we have used the assumption that $\mathcal{L}^B$ is anti-$\mathcal{PC}-$ symmetric and recalled the matrix definitions \eqref{J_defn} and \eqref{A_defn1}:
\begin{align}\label{matr-k2}
\mathcal M^{{\bf K}_\star}_{\mathscr{A}}({\bm{\mathfrak{K}}}_2) &=
\upsilon_{_F} \begin{pmatrix}0 & {\bm{\mathfrak{K}}}_2^{(1)}+i{\bm{\mathfrak{K}}}_2^{(2)}\\
{\bm{\mathfrak{K}}}_2^{(1)}-i{\bm{\mathfrak{K}}}_2^{(2)}&0\end{pmatrix},\quad \text{and} \quad
\mathcal M^{{\bf K}_\star}_{\mathcal{L}^B}
= {\vartheta_{\sharp}}^{{\bf K}_\star} \sigma_3 .
\end{align}

\begin{proposition}
\label{0E_prop}
Let $\eta(\zeta)$ be a domain wall function (Definition \ref{domain_wall_defn}), ${\bf K}_\star$ a point of ${\bf K}$ or ${\bf K}'$ type, and assume that ${\vartheta_{\sharp}}^{{\bf K}_\star}\neq0$. Then:
\begin{enumerate}
\item The Dirac operator $\mathcal{D}^{{\bf K}_\star}$ \eqref{dirac_op} possesses a zero-energy eigenvalue, $E^{(1)}=0$,  with exponentially localized eigenfunction given by:
\begin{align}
\label{0mode_theta}
{\bm{\alpha}}^{{\bf K}_\star}_\star(\zeta) &=
 \left\{ \begin{array}{ll}
\gamma\ e^{-\frac{\abs{{\vartheta_{\sharp}}^{{\bf K}_\star}}}{{\upsilon_{_F}} | {\bm{\mathfrak{K}}}_2|}\ \int_0^\zeta \eta(s) ds}\ \chi_-({\bm{\mathfrak{K}}}_2^\perp)
& \quad \text{if} \quad {\vartheta_{\sharp}}^{{\bf K}_\star}>0; \\
\gamma\ e^{-\frac{\abs{{\vartheta_{\sharp}}^{{\bf K}_\star}}}{{\upsilon_{_F}} | {\bm{\mathfrak{K}}}_2|}\ \int_0^\zeta \eta(s) ds}\ \chi_+({\bm{\mathfrak{K}}}_2^\perp)
& \quad \text{if} \quad {\vartheta_{\sharp}}^{{\bf K}_\star}<0, \\
\end{array} \right. \quad \text{where} \\
& \chi_{\pm}({\bm{\mathfrak{K}}}_2^\perp)\ =\ \frac{1}{\sqrt2}\ \begin{pmatrix}\widehat{\mathfrak{z}}({\bm{\mathfrak{K}}}_2^\perp)\\ \pm1\end{pmatrix}, \nonumber
\end{align}
and $\widehat{\mathfrak{z}}({\bm{\mathfrak{K}}})= ({\mathfrak{K}}^{(1)}+i{\mathfrak{K}}^{(2)})/|{\bm{\mathfrak{K}}}|$.
The normalization constant,
\begin{equation}\label{gamma-def}
\gamma=\left(
\int_{-\infty}^\infty
\exp\left(
-2\ \frac{\abs{{\vartheta_{\sharp}}^{{\bf K}_\star}}}{\upsilon_{F} | {\bm{\mathfrak{K}}}_2|}\ \int_0^\zeta \eta(s) ds\right) d\zeta \right)^{-{\frac{1}{2}}} ,
\end{equation}
 is a real and positive constant and chosen so that $\norm{{\bm{\alpha}}^{{\bf K}_\star}_\star}_{L^2(\mathbb{R})}=1$.
\item The zero-energy mode \eqref{0mode_theta} is ``topologically protected'' in the sense that this mode persists against sufficiently spatially localized (even large) perturbations of  $\eta(\zeta)$,
 which preserve the asymptotic behavior, $\pm\eta_\infty$,  as $\zeta\to\pm\infty$.
\end{enumerate}
\end{proposition}

Proposition \ref{0E_prop} is proved below. This proposition and the formal expansion preceding it imply the following result on the bifurcation of edge states, which is topologically protected against arbitrary localized perturbations of the domain-wall, defined by $\eta$:

\begin{theorem}[Formal expansion of edge states]
\label{edgemode_thm}
Assume:
\begin{enumerate}
\item $A({\bf x})$ is smooth, $\Lambda_h-$ periodic and Hermitian, and such that $[\mathcal{PC},{\mathcal{L}^A}]=0$ and $[\mathcal{R},{\mathcal{L}^A}]=0$.
 \item  $({\bf K}_\star,E_D)$ is a Dirac point of ${\mathcal{L}^A}=-\nabla\cdot A\nabla$ (see Definition \ref{dirac_pt_defn} and Theorems \ref{low-dp} and \ref{arbitrary-dp}), where
  ${\bf K}_\star$ is a vertex of the hexagonal Brillouin zone, $\mathcal{B}_h$.
\item $\mathcal{L}^B \equiv -\nabla\cdot B\nabla$ is anti-$\mathcal{PC}-$ symmetric, and such that ${\vartheta_{\sharp}}^{{\bf K}_\star}\neq0$; see \eqref{thta_shp1}.
\end{enumerate}
Let $\mathcal{L}_{\rm dw}^{(\delta)} \equiv -\nabla\cdot\left[A({\bf x}) +\delta\eta( \delta{\bm{\mathfrak{K}}}_2 \cdot {\bf x})B({\bf x})\right]  \nabla$,
denote the domain-wall operator, associated with the line-defect across the rational edge, $\mathbb{R}{\bm{\mathfrak{v}}}_1$.

Then, the edge-state eigenvalue problem for $\mathcal{L}_{\rm dw}^{(\delta)}$, \eqref{dw_evp}-\eqref{localized}, for fixed parallel quasi-momentum ${k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$, has the formal, \underline{topologically protected} eigenpair solution $(E^\delta_{{\bf K}_\star},\Psi^{\delta}_{{\bf K}_\star})$,
corresponding to a state which propagates in the ${\bm{\mathfrak{v}}}_1-$ direction with parallel quasi-momentum ${k_{\parallel}} = {{\bf K}_\star} \cdot {\bm{\mathfrak{v}}}_1$, and is exponentially decaying in the transverse direction, as ${\bm{\mathfrak{K}}}_2\cdot {\bf x} \to \pm \infty$.

Furthermore, the eigenpair $(E^\delta_{{\bf K}_\star},\Psi^{\delta}_{{\bf K}_\star})$ can be expanded to any finite order in $\delta$ in powers of $\delta$. To leading order $E^\delta_{{\bf K}_\star}=E_D+\mathcal{O}(\delta^2)$ and
\begin{equation}
 \label{ldord}
 \Psi^{\delta}_{{\bf K}_\star}({\bf x}) =
 \left\{ \begin{array}{ll}
 \delta^{1/2}\ \gamma\ \Big[ \chi_-({\bm{\mathfrak{K}}}_2^\perp)\cdot\left(\Phi_1^{{\bf K}_\star}({\bf x}),\Phi_2^{{\bf K}_\star}({\bf x})\right) \Big] \
 &e^{-\frac{\abs{{\vartheta_{\sharp}}^{{\bf K}_\star}}}{\upsilon_F|{\bm{\mathfrak{K}}}_2|}\int_0^{\delta {\bm{\mathfrak{K}}}_2\cdot {\bf x}}\eta(s)d s} + \mathcal{O}(\delta), \\
 & \quad \text{if} \quad {\vartheta_{\sharp}}^{{\bf K}_\star}>0; \\\\
 \delta^{1/2}\ \gamma\ \Big[ \chi_+({\bm{\mathfrak{K}}}_2^\perp)\cdot\left(\Phi_1^{{\bf K}_\star}({\bf x}),\Phi_2^{{\bf K}_\star}({\bf x})\right) \Big] \
 &e^{-\frac{\abs{{\vartheta_{\sharp}}^{{\bf K}_\star}}}{\upsilon_F|{\bm{\mathfrak{K}}}_2|}\int_0^{\delta {\bm{\mathfrak{K}}}_2\cdot {\bf x}}\eta(s)d s} + \mathcal{O}(\delta), \\
 & \quad \text{if} \quad {\vartheta_{\sharp}}^{{\bf K}_\star}<0. \\\\
\end{array} \right.
\end{equation}
The constant $\gamma$ is displayed in \eqref{gamma-def}, and the factor of $\delta^{1/2}$ ensures that the edge state is normalized to unity.
\end{theorem}

\begin{proof}[Proof of Proposition \ref{0E_prop}]
Setting $E^{(1)}=0$, equation \eqref{dirac_eqn} may be rewritten as:
\[\partial_\zeta{\bm{\alpha}}^{{\bf K}_\star}(\zeta)\ =\ -i\ \eta(\zeta)\ \mathcal{M}^{{\bf K}_\star}_{\mathscr{A}}({\bm{\mathfrak{K}}}_2)^{-1}\ \mathcal{M}^{{\bf K}_\star}_{{\mathcal{L}^B}}\ {\bm{\alpha}}^{{\bf K}_\star}(\zeta).\]
Using the expressions displayed in \eqref{matr-k2} we obtain:
\begin{equation}
\label{dirac_0E}
\partial_\zeta {\bm{\alpha}}^{{\bf K}_\star}(\zeta)\ =\
\frac{{\vartheta_{\sharp}}^{{\bf K}_\star}}{\upsilon_F^3|{\bm{\mathfrak{K}}}_2|^2}\ \eta(\zeta)\ \mathcal{M}^{{\bf K}_\star}_{\mathscr{A}}({\bm{\mathfrak{K}}}_2^{\perp})\ {\bm{\alpha}}^{{\bf K}_\star}(\zeta),
\end{equation}
where ${\bm{\mathfrak{K}}}_2=({\bm{\mathfrak{K}}}_2^{(1)},{\bm{\mathfrak{K}}}_2^{(2)})$ and ${\bm{\mathfrak{K}}}_2^\perp=(-{\bm{\mathfrak{K}}}_2^{(2)},{\bm{\mathfrak{K}}}_2^{(1)})$.

The system \eqref{dirac_0E} may be diagonalized using the eigenpairs of $\mathcal{M}^{{\bf K}_\star}_{\mathscr{A}}({\bm{\mathfrak{K}}}_2^{\perp})$.
 Note that $\mathcal{M}^{{\bf K}_\star}_{\mathscr{A}}({\bm{\mathfrak{K}}}_2^{\perp})$ has two eigenpairs:
\begin{align}\label{chip+}
\mu_+ &=+\upsilon_F|\ {\bm{\mathfrak{K}}}_2^{\perp}|,\qquad \chi_+({\bm{\mathfrak{K}}}_2^{\perp})=\frac{1}{\sqrt{2}}
\begin{pmatrix}
\widehat{\mathfrak{z}}({\bm{\mathfrak{K}}}_2^{\perp})\\ 1
\end{pmatrix},\  \text{and}\\
\mu_-&=-\upsilon_F|\ {\bm{\mathfrak{K}}}_2^{\perp}|,\qquad \chi_-({\bm{\mathfrak{K}}}_2^{\perp})=\frac{1}{\sqrt{2}}
\begin{pmatrix}
\widehat{\mathfrak{z}}({\bm{\mathfrak{K}}}_2^{\perp})\\ -1
\end{pmatrix},
\label{chip-}
\end{align}
where $\widehat{\mathfrak{z}}({\bm{\mathfrak{K}}})= ({\mathfrak{K}}^{(1)}+i{\mathfrak{K}}^{(2)})/|{\bm{\mathfrak{K}}}|$.

Therefore, \eqref{dirac_0E} has the general solution
\begin{equation}\label{zeromode_generalsolution}
{\bm{\alpha}}^{{\bf K}_\star}_\star(\zeta)=c_1\ e^{\frac{{\vartheta_{\sharp}}^{{\bf K}_\star}}{\upsilon|{\bm{\mathfrak{K}}}_2|}\ \int^{\zeta}_0\eta(s)ds}\chi_+({\bm{\mathfrak{K}}}_2^{\perp})+c_2\ e^{-\frac{{\vartheta_{\sharp}}^{{\bf K}_\star}}{\upsilon|{\bm{\mathfrak{K}}}_2|}\ \int^{\zeta}_0\eta(s)ds}\chi_-({\bm{\mathfrak{K}}}_2^{\perp})
\end{equation}
If ${\vartheta_{\sharp}}^{{\bf K}_\star}>0$, we obtain an exponentially localized solution of \eqref{dirac_0E} by setting $c_1=0$ in \eqref{zeromode_generalsolution}, and if ${\vartheta_{\sharp}}^{{\bf K}_\star}<0$, we obtain the exponentially localized solution by setting $c_2=0$ in \eqref{zeromode_generalsolution}. Imposing the normalization $\norm{{\bm{\alpha}}^{{\bf K}_\star}_\star}_{L^2(\mathbb{R})}=1$ concludes the proof of Proposition \ref{0E_prop}.
\end{proof}

\subsection{Rigorous formulation of Theorem \ref{edgemode_thm}, the spectral no-fold condition,
and outline of the proof}\label{rigor}

A rigorous reformulation of Theorem \ref{edgemode_thm} and a detailed proof follows the strategy implemented in full detail for the case of Schr\"odinger operators in \cites{FLW-AnnalsPDE:16}. These arguments can be adapted to the current setting of divergence form operators. In this section we provide a detailed sketch of the construction, motivate the key {\it spectral no-fold condition}, and conclude the section with the reformulation in Theorem \ref{rig-edge}.

We begin by formulating the edge state eigenvalue problem  in an appropriate Hilbert space.
For a fixed edge, $\mathbb{R}{\bm{\mathfrak{v}}}_1$, introduce the cylinder $\Sigma\equiv \mathbb{R}^2/ \mathbb{Z}{\bm{\mathfrak{v}}}_1$. A function on $\Sigma$ is invariant under the shift  ${\bf x}\to{\bf x}+{\bm{\mathfrak{v}}}_1$. For $s\ge0$, denote by  $H^s(\Sigma)$,  the  Sobolev space of order $s$  of  functions on $\Sigma$. Note $L^2(\Sigma)=H^0(\Sigma)$.
The pseudo-periodicity and decay conditions \eqref{pseudo-per}-\eqref{localized} are encoded by requiring $ \Psi \in H^s_{k_{\parallel}}(\Sigma)$ for some $s\ge0$, where
\begin{equation}
\label{H_kpar_def}
H^s_{k_{\parallel}}(\Sigma)\equiv \ \left\{f : f({\bf x})e^{-i\frac{{k_{\parallel}}}{2\pi}{\bm{\mathfrak{K}}}_1\cdot{\bf x}}\in H^s(\Sigma) \right\}.
\end{equation}

The eigenvalue problem \eqref{dw_evp}-\eqref{localized} can then be reformulated as:
\begin{equation}
 \label{EVP}
 \mathcal{L}_{\rm dw}^{(\delta)} \Psi = E \Psi , \quad \Psi\in H_{{k_{\parallel}}}^2(\Sigma) .
\end{equation}

Edge states for ${k_{\parallel}}$ in a neighborhood of ${k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$ can be constructed perturbatively.
\medskip

In analogy to Theorem 4.2 of \cites{FLW-AnnalsPDE:16}, any function in $L^2_{{k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1}(\Sigma)$ can be decomposed into a superposition of Floquet-Bloch modes with the pseudo-periodicity condition $f({\bf x}+{\bm{\mathfrak{v}}}_1)=e^{i{\bf K}\cdot{\bm{\mathfrak{v}}}_1} f({\bf x})$. The set of all such modes is given by:
$\{\Phi_b({\bf x};{\bf K}_\star+\lambda{\bm{\mathfrak{K}}}_2)\}_{b\ge1}$ and where $|\lambda|\le1/2$.

\begin{theorem}\label{fourier-edge}
Let $A({\bf x})$ be periodic with respect to the triangular lattice, $\Lambda_h$, and such that ${\mathcal{L}^A}=-\nabla_{\bf x}\cdot A({\bf x})\ \nabla_{\bf x}$ is strictly elliptic on $\mathbb{R}^2$. In particular, we may take $A({\bf x})$ to be a honeycomb structured medium.
Let $f\in L_{\kpar=\bK_\star\cdot\vtilde_1}^2(\Sigma)= L_{\kpar=\bK_\star\cdot\vtilde_1}^2(\mathbb{R}^2/\mathbb{Z}{\bm{\mathfrak{v}}}_1)$. Then,
\begin{enumerate}
\item  $f$ can be represented as a superposition of Floquet-Bloch modes of ${\mathcal{L}^A}$ with quasimomenta in $\mathcal{B}_h$ located on the segment $\lambda\mapsto {\bf K}_\star+\lambda{\bm{\mathfrak{K}}}_2$ for $|\lambda|\le\frac{1}{2}$.
\begin{align}
f({\bf x})
&=\sum_{b\ge1}\int_{-\frac12}^{\frac12} \widetilde{f}_b(\lambda) \Phi_b({\bf x};{\bf K}_\star+\lambda {\bm{\mathfrak{K}}}_2) d\lambda  \nonumber \\
&= e^{i{\bf K}_\star\cdot{\bf x}} \sum_{b\ge1}\int_{-\frac12}^{\frac12} e^{i\lambda{\bm{\mathfrak{K}}}_2\cdot{\bf x}}\widetilde{f}_b(\lambda) p_b({\bf x};{\bf K}+\lambda {\bm{\mathfrak{K}}}_2) d\lambda,\qquad {\rm where} \label{fb-Sigma} \\
\widetilde{f}_b(\lambda) \ &=\ \left\langle \Phi_b(\cdot,{\bf K}_\star+\lambda{\bm{\mathfrak{K}}}_2),f(\cdot)\right\rangle_{L_{\kpar=\bK\cdot\vtilde_1}^2} . \nonumber
\end{align}
Here, the sum representing $e^{-i{\bf K}_\star\cdot{\bf x}}f({\bf x})$, in \eqref{fb-Sigma} converges in the $L^2(\Sigma)$ norm.
\item In the special case where $A({\bf x})=I_{2\times2}$\ (${\mathcal{L}^A}=-\Delta$):
 \begin{align*}
 f({\bf x}) = \sum_{{{\bf m}}\in\mathbb{Z}^2}e^{ i({\bf K}_\star+{\bf m}\vec{\bm{\mathfrak{K}}})\cdot{\bf x} } \int_{-\frac12}^{\frac12}
 \widehat{f}_{\bf m}(\lambda)e^{i\lambda{\bm{\mathfrak{K}}}_2\cdot{\bf x}} d\lambda\ .
  \end{align*}
 Here, for ${\bf m}=(m_1,m_2)\in\mathbb{Z}^2$, we define ${\bf m}\vec{\bm{\mathfrak{K}}}=m_1{\bm{\mathfrak{K}}}_1+m_2{\bm{\mathfrak{K}}}_2$.
 \end{enumerate}
\end{theorem}

We now seek $(\Psi,E)$, a solution of the eigenvalue problem \eqref{EVP}, in the form of a truncated multiscale expansion plus a corrector to be constructed:
\begin{align}
\Psi^\delta &\equiv \Psi^{(0)}({\bf x},\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})+
\delta\Psi^{(1)}({\bf x},\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})+\delta U({\bf x}) , \label{U-def}\\
E^\delta &\equiv E_D + \delta^2\mu \label{mu-def},
\end{align}
where  $ \Psi^{(0)}({\bf x},\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x}) = \alpha^{{\bf K}_\star}_1(\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})\Phi^{{\bf K}_\star}_1({\bf x})+\alpha^{{\bf K}_\star}_2(\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})\Phi^{{\bf K}_\star}_2({\bf x})
$. The pair  ${\bm{\alpha}}^{{\bf K}_\star}(\zeta)=(\alpha^{{\bf K}_\star}_1(\zeta),\alpha^{{\bf K}_\star}_2(\zeta))^T$ is a zero energy eigenstate of a Dirac operator: $\mathcal{D}^{{\bf K}_\star}{\bm{\alpha}}^{{\bf K}_\star}(\zeta) = 0$, ${\bm{\alpha}}^{{\bf K}_\star}\in L^2(\mathbb{R})$,
where $\mathcal{D}^{{\bf K}_\star}$ is displayed in \eqref{dirac_op}; see \eqref{0mode_theta} of Proposition \ref{0E_prop}. $\Psi^{(1)}({\bf x},\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})$ is the solution of \eqref{Eq_delta1} with $E^{(1)}=0$. Both  $\Psi^{(0)}({\bf x},\zeta)$ and $\Psi^{(1)}({\bf x},\zeta)$ are constructed
to be ${\bf K}_\star-$ pseudo-periodic as functions of ${\bf x}$ and decaying as functions of $\zeta$. Finally, we seek the $U\in L^2_{{k_{\parallel}}={\bm{\mathfrak{v}}}_1\cdot{\bf K}}$.

Substitution of \eqref{U-def} and \eqref{mu-def} into \eqref{EVP} yields an equation for the corrector $(\mu,U({\bf x}))$.
 By Theorem \ref{fourier-edge}, we have
 \begin{equation}
 U({\bf x})=\sum_{b\ge1}\int_{-\frac12}^{\frac12} \widetilde{U}_b(\lambda) \Phi_b({\bf x};{\bf K}_\star+\lambda {\bm{\mathfrak{K}}}_2) d\lambda\ .
 \label{U-exp}
 \end{equation}
The expression in \eqref{U-exp} is a superposition of modes corresponding to the {\it band structure slice at quasi-momentum ${\bf K}_\star$, which is dual to the edge $\mathbb{R}{\bm{\mathfrak{v}}}_1$}:
\begin{equation*}
\lambda\mapsto E_b({\bf K}_\star+\lambda{\bm{\mathfrak{K}}}_2),\ \ |\lambda|\le 1/2,\ \ b\ge1.
\end{equation*}
 Figure \ref{E_k_slices} displays the first three dispersion curves of band structure slices for the zigzag and armchair edges, respectively.

The equation for $U({\bf x})$ is equivalent to the coupled system for the amplitudes $\{\widetilde{U}_b(\lambda)\}_{b\ge1},\ |\lambda|\le1/2$, associated with the Floquet-Bloch modes
 $\{\Phi_b({\bf x};{\bf K}_\star+\lambda{\bm{\mathfrak{K}}}_2)\}_{b\ge1},\ |\lambda|\le1/2$:

\begin{equation}
\label{amp_b}
 \begin{split}
& \left( E_b(\lambda)  - E_D \right) \widetilde{U}_b(\lambda) -
\delta \left\langle \Phi_b^{{\bf K}_\star}(\cdot;\lambda) , \nabla_{\bf x}\cdot \left[ \eta(\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})  B({\bf x}) \nabla_{\bf x} U({\bf x}) \right]  \right\rangle_{L^2_{{\kpar=\bK\cdot\vtilde_1}}(\Sigma_{\bf x})} \\
& \qquad = \delta\widetilde{F}_b[\mu,\delta](\lambda)  + \delta^2\ \mu\ \widetilde{U}_b(\lambda) , \qquad b\ge1, \ |\lambda|\le1/2 .
\end{split}
\end{equation}

Here, $\widetilde{F}_b[\mu,\delta](\lambda)
$ is an expression involving  $\Psi^{(j)}({\bf x},\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})$, $j=0,1$, and their derivatives, projected onto
the Floquet-Bloch mode, $\Phi_b^{{\bf K}_\star}({\bf x};{\bf K}_\star+\lambda{\bm{\mathfrak{K}}}_2)$. A solution $\{\widetilde{U}_b(\lambda): b\ge1,\ |\lambda|\le1/2\}$ of the system \eqref{amp_b} is sought which satisfies $\sum_{b\ge1} \int_{-1/2}^{1/2}\ (1+|\lambda|^2)^s\ |\widetilde{U}_b(\lambda)|^2\ d\lambda<\infty$, for some $s\ge0$.

Suppose that the Dirac point $({\bf K}_\star,E_D)$ occurs at the intersection of the $b_*^{th}$ and $(b_*+1)^{st}$ spectral bands of ${\mathcal{L}^A}$.

Following the strategy of \cite{FLW-AnnalsPDE:16}, we decompose the system \eqref{amp_b} into two subsystems consisting of near- and far- quasi-momentum Floquet-Bloch amplitudes:
\begin{itemize}
 \item $\widetilde{U}_{{\rm near}}(\lambda)=\{\widetilde{U}_{{\rm near},-}(\lambda) , \widetilde{U}_{{\rm near},+}(\lambda)\}$, for $|\lambda|\le \delta^\nu$: Floquet-Bloch amplitudes corresponding to quasi-momenta which are near the Dirac point, where
 $\{-,+\} \equiv \{b_\star, b_{\star}+1\}$.
 \item $\{ \widetilde{U}_{{\rm far},b}(\lambda)\}_{_{b\ge1}}$ for $1/2\ge|\lambda|\ge \left(\delta_{b,b_\star}+\delta_{b,b_\star+1}\right) \delta^\nu$: Floquet-Bloch amplitudes corresponding to quasi-momenta which are bounded away from the Dirac point.
\end{itemize}
Here, $ \nu>0$ is appropriately chosen and fixed.

To construct a solution for \eqref{amp_b}, we proceed via a Lyapunov-Schmidt reduction strategy.
The first step is to solve for the far-amplitudes in terms of the near-amplitudes, by constructing the mapping: $\widetilde{U}_{{\rm near}}[\mu,\delta]\mapsto \widetilde{U}_{{\rm far}}[\widetilde{U}_{{\rm near}},\mu,\delta]$.
The second step is to substitute this mapping into the near-amplitude equations to obtain a closed system of two
(nonlocal) equations for $\widetilde{U}_{{\rm near},-}(\lambda)$ and $\widetilde{U}_{{\rm near},+}(\lambda)$, for $|\lambda|\le \delta^\nu$. An appropriate rescaling of this closed system  can be solved for $\delta$ sufficiently small. In this way, the corrector $(U^\delta,\mu^\delta)$ can be constructed.

The construction  of the mapping $\widetilde{U}_{{\rm near}}[\mu,\delta]\mapsto \widetilde{U}_{{\rm far}}[\widetilde{U}_{{\rm near}},\mu,\delta]$ and the necessary estimates on this mapping requires the \emph{spectral no-fold condition},  introduced in \cites{FLW-AnnalsPDE:16}. Specifically, we require that  for $\lambda$ bounded away from $\lambda=0$ or $b\notin\{b_*,b_*+1\}$, the energy difference $E_{b}(\lambda)  - E_D$ is uniformly bounded away from zero.
This condition ensures that, for $\delta\ne0$, the $L^2_{{k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1}-$ spectra of the asymptotic operators $\mathcal{L}_-^{(\delta)}$ and  $\mathcal{L}_+^{(\delta)}$ (see \eqref{asymptotic_ops}) have a full spectral gap containing the energy $E=E_D$. In one-dimensional periodic problems with Dirac points (linear band crossings), the spectral no-fold condition holds automatically; see \cites{FLW-PNAS:14,FLW-MAMS:17}.

The left panel of Figure \ref{E_k_slices} displays the first three bands of the zigzag slice of the band structure for a honeycomb structure, $A({\bf x})$.
By Theorem \ref{PT_thm}, the corresponding slice of the band structure of $\mathcal{L}^{(\delta)}_\pm=-\nabla_{\bf x}\cdot\left(A({\bf x})\pm\delta\eta_\infty B({\bf x})\right)\nabla_{\bf x}$ will have small local ($\lambda$ near 0) gap with energy $E_D$ in its interior.
Since the energies ($E$) attained by the dispersion curves of ${\mathcal{L}^A}$, outside a neighborhood of $\lambda=0$, are bounded away from $E_D$, it follows that for $\delta$ sufficiently small  $\mathcal{L}^{(\delta)}$ has a full $L^2_{{k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1}-$ spectral gap, an open interval about $E_D$. That is, the spectral no-fold condition holds in this case.

The right panel of Figure \ref{E_k_slices} displays the first three bands of the armchair slice for the same choice of $A({\bf x})$.
Again, by Theorem \ref{PT_thm}, the corresponding slices of the band structures of $\mathcal{L}^{(\delta)}_\pm$ have a local spectral gap near ${\bf k}={\bf K}$ ($\lambda=0$). However, the intersections of the horizontal dashed line away from $\lambda=0$ indicate that for $\delta$ small there is no $L^2_{{k_{\parallel}}={\bf K}\cdot{\bm{\mathfrak{v}}}_1}-$ spectral gap; the spectral no-fold condition fails for the armchair edge.

\begin{remark}\label{high_contrast_limit}
A honeycomb Schr\"odinger operator is an operator of the form $H^\lambda=-\Delta +\lambda^2 V({\bf x})$,
where $V$ has the symmetries of a hexagonal tiling of the plane \cites{FW:12}. An important special case is where
$V({\bf x})$ is a superposition of translates of  an `atomic' potential well. For any fixed rational edge,  the spectral no-fold condition as stated in  \cites{FLW-AnnalsPDE:16}, holds for all $\lambda$  sufficiently large.
This follows from the scaled-convergence of the low-lying dispersion surfaces to those of the tight-binding limiting model, proved in \cites{FLW-CPAM:17}.  In the setting of Maxwell's equations, numerical studies
indicate that for high contrast the spectral no-fold condition holds for the zig-zag edge. However, to date we have not found this condition to hold for other edges. It is an interesting question to determine for which edges the no-fold condition holds. This has implications for energy localization and is currently under investigation.  \end{remark}

The above strategy is implemented along the lines of the Schr\"odinger operator setting \cites{FLW-AnnalsPDE:16}.
An essential difference between the Schr\"odinger setting and the present Maxwell setting
 is that the perturbation, $\delta \nabla_{\bf x}\cdot \eta(\delta{\bm{\mathfrak{K}}}_2\cdot{\bf x})  B({\bf x})\nabla_{\bf x} U({\bf x})$, which is projected  in \eqref{amp_b}, is a second order differential operator acting on the unknown $U$; compare with  Equation (7.12) of \cites{FLW-AnnalsPDE:16}.
This difference results in simple changes in the analysis of the closed system for the near-amplitudes, $\widetilde{U}_{{\rm near},-}(\lambda)$ and $U_{{\rm near},+}(\lambda)$ (step two), and to a different definition of ${\vartheta_{\sharp}}^{{\bf K}_\star}$; see \eqref{thta_shp}). For the  Schr\"odinger case,  see Proposition 6.2 of \cites{FLW-AnnalsPDE:16}.
To obtain the  far-amplitudes as a functional of the near-amplitudes (step one), we construct the mapping $\widetilde{U}_{{\rm near}}[\mu,\delta]\mapsto \widetilde{U}_{{\rm far}}[\widetilde{U}_{{\rm near}},\mu,\delta]$ by solving the fixed point equation:
 $\mathcal{Q}_\delta \phi = \phi$, where
 {\footnotesize{
  \begin{align}
&\widetilde{\left[ \mathcal{Q}_\delta \phi \right]}_b(\lambda)\nonumber\\
&\quad \equiv -\delta \frac{ \chi(\abs{\lambda}\geq(\delta_{b,b_{\star}}+\delta_{b,b_{\star}+1})\delta^{\nu})}{E_b(\lambda)-E_{\star}}
\inner{\Phi_{b}({\bf x},\lambda),\nabla_{{\bf x}}\cdot \eta(\delta {\bm{\mathfrak{K}}}_2 {\bf x})B({\bf x}) \nabla_{\bf x} \phi ({\bf x})}_{L_{\kpar=\bK\cdot\vtilde_1}^2} \nonumber \\
&\quad + \delta^2\mu \frac{\widetilde{\phi}_b(\lambda)}{{E_b(\lambda)-E_{\star}}} ;  \label{tQ-def}
\end{align}
}}
For  $s\geq0$, introduce the Sobolev space of order $s$ consisting of  $H^s_{{\kpar=\bK_\star\cdot\vtilde_1}}(\Sigma)$ functions with quasi-momenta $\lesssim\delta^\nu$:
$H^s_{\rm near, \delta^\nu} \equiv \{f \in H^s_{{\kpar=\bK_\star\cdot\vtilde_1}}(\Sigma) : \widetilde{f}_b(\lambda) = \chi(|\lambda| \leq \delta^\nu) \widetilde{f}_b(\lambda)\}$.
One can verify, using that $({\mathcal{L}^A}+I)^{-1}\ \nabla_{{\bf x}}\cdot \eta(\delta {\bm{\mathfrak{K}}}_2 {\bf x})B({\bf x}) \nabla_{\bf x}$  is a bounded operator on $H^2_{{\kpar=\bK_\star\cdot\vtilde_1}}$, that  $\mathcal{Q}_\delta$ maps $H^2_{\rm near, \delta^\nu}$ to $H^2_{\rm near, \delta^\nu}$ with norm bounded by ${\rm constant}\times\mathfrak{e}(\delta)$,%
where
 \begin{align}
&\mathfrak{e}(\delta)\equiv \sup_{b=\pm}\ \ \sup_{|\delta|^\nu\le|\lambda|\le\frac12}\
 \frac{|\delta|}{|E_b(\lambda)-E_{\star}|}\nonumber\\
 &\quad  +\ \sup_{ b\ge1,\ b\ne\pm}\ (\ 1+|b|\ ) \sup_{0\le|\lambda|\le\frac12}
 \frac{|\delta|}{|E_b(\lambda)-E_{\star}|} .
\label{frak-e-def}
\end{align}
The spectral no-fold condition ensures that $\mathfrak{e}(\delta)=o(1)$ as $\delta\rightarrow0$.
Hence, by the contraction mapping principle the fixed point equation $\mathcal{Q}_\delta \phi = \phi$ may be solved, for $\delta>0$ and sufficiently small, on an appropriately sized ball in $H^2_{\rm near, \delta^\nu}$.

With these observations, the above strategy of proof yields the existence of edge states for the Maxwell system:

\begin{theorem}\label{rig-edge}
 Consider the ${k_{\parallel}}={\bf K}\cdot{\bm{\mathfrak{v}}}_1-$ edge state eigenvalue problem
\eqref{EVP}.
There exist positive constants $\delta_0, c_0$ and a branch of solutions of \eqref{EVP},
   \[|\delta|\in(0,\delta_0)\longmapsto (E^\delta,\Psi^\delta)\in (E_D-c_0\ \delta_0\ ,\ E_D+c_0\ \delta_0)\times H^2_{{\kpar=\bK_\star\cdot\vtilde_1}}(\Sigma) , \]
   such that the following holds:
\begin{enumerate}
\item $\Psi^\delta$ is well-approximated by a slow modulation of the
degenerate Floquet-Bloch modes $\Phi_1^{{\bf K}_\star}$ and $\Phi_2^{{\bf K}_\star}$,
which decays to zero transverse to the edge, $\mathbb{Z}{\bm{\mathfrak{v}}}_1$:
\begin{align}
&\left\|\ \Psi^\delta(\cdot)\ -\
\left[\alpha^{{\bf K}_\star}_{\star,1}(\delta{\bm{\mathfrak{K}}}_2 \cdot)\Phi_1^{{\bf K}_\star}(\cdot)+\alpha^{{\bf K}_\star}_{\star,2}(\delta{\bm{\mathfrak{K}}}_2 \cdot )\Phi_2^{{\bf K}_\star}(\cdot)\right]\
\right\|_{H^2_{{\kpar=\bK_\star\cdot\vtilde_1}}}\
\lesssim\ \delta^{\frac12} , \label{TM_Psi-error}\\
& E^\delta = E_D + \mathcal{O}(\delta^2) . \label{TM_E-error}
\end{align}
\item The amplitude vector,  ${\bm{\alpha}}^{{\bf K}_\star}_\star(\zeta)=\left(\alpha^{{\bf K}_\star}_{\star,1}(\zeta),\alpha^{{\bf K}_\star}_{\star,2}(\zeta)\right)$, is an  $L^2(\mathbb{R}_\zeta)- $
normalized, topologically protected zero-energy eigenstate of the Dirac operator, $\mathcal{D}^{{\bf K}_\star}$; see \eqref{dirac_op}.
\end{enumerate}
\end{theorem}

\begin{figure}
\centering
\includegraphics[width=4in]{E_k_slices}
\caption{\footnotesize
Band dispersion slices of ${\mathcal{L}^A}$ along the quasi-momentum segments:  ${\bf K}+\lambda{\bm{\mathfrak{K}}}_2,\ |\lambda|\le1/2$, for ${\bm{\mathfrak{K}}}_2={\bf k}_2$ (zigzag) and ${\bm{\mathfrak{K}}}_2=-{\bf k}_1+{\bf k}_2$ (armchair).
Here, $A({\bf x})=a(x)I$ where $a(x)$ is a particular honeycomb lattice function.
Dirac point energy levels $E=E_D$ are indicated with dotted red lines.
The spectral no-fold condition is satisfied for the zigzag slice (left panel) but not the armchair slice (right panel).}
\label{E_k_slices}
\end{figure}

\section{Edge state dispersion curves}\label{edge_state_kpar}

Under the stated hypotheses, Theorem \ref{edgemode_thm} guarantees the existence of protected edge states for $\delta>0$ sufficiently small and ${k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$ fixed. Building off this result, perturbation theory can be used to construct edge state eigenpairs, $(\Psi(\cdot;{k_{\parallel}}),E({k_{\parallel}}))$, of the eigenvalue problem \eqref{dw_evp}-\eqref{localized} for all ${k_{\parallel}}$ in a neighborhood of ${k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$:

\begin{theorem}
 \label{near-kpar1}
 Assume hypotheses (1)-(3) of Theorem \ref{edgemode_thm}. Let ${\bf K}_\star$ be of the ${\bf K}$ or ${\bf K}'$ type. Fix $\delta$ sufficiently small.
 Then there exists $c_0\ll\delta$ such that for all ${k_{\parallel}}$ satisfying $|{k_{\parallel}} - {\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1| < c_0$, there exists a
 ${k_{\parallel}}\mapsto \left(E_{{\bf K}_\star}({k_{\parallel}}),\Psi_{{\bf K}_\star}({\bf x};{k_{\parallel}})\right)$ edge state solution of the eigenvalue problem \eqref{dw_evp}-\eqref{localized} (equivalently \eqref{EVP}). Here, the energy
 $E^\delta_{{\bf K}_\star}({k_{\parallel}})= E^\delta_{{\bf K}_\star} + \mu^{\delta} ({k_{\parallel}}-{\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1) +\mathcal{O}(|{k_{\parallel}}-{\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1|^2)$, where $E^{\delta}_{{\bf K}_\star}=E_D+\mathcal{O}(\delta^2)$ is the eigenvalue of the edge state given in \eqref{edgemode_thm}, and $\mu^\delta$ is a constant, which is independent of ${k_{\parallel}}$.
\end{theorem}

By taking a continuous superposition of the edge states given by Theorem \ref{near-kpar1}, we obtain dispersive wavepackets which are concentrated along the edge. These wavepackets have parallel-quasi-momentum components concentrated near ${\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$, and their group velocity is determined by $\frac{dE^{\delta}({k_{\parallel}})}{d{k_{\parallel}}}$, for ${k_{\parallel}}$ near ${\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$.

\begin{theorem}\label{gp-vel}
Assume hypotheses (1)-(3) of Theorem \ref{edgemode_thm}. Let ${\bf K}_\star$ be of the ${\bf K}$ or ${\bf K}'$ type.
Then,
\begin{equation}
\label{GV2_thm}
 \frac{d E^{\delta}}{d{k_{\parallel}}}\Big|_{_{{k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1}} =
 \left\{
 \begin{array}{ll}
 - \frac{1}{2\pi}\frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}({\bm{\mathfrak{K}}}_2\wedge{\bm{\mathfrak{K}}}_1) + \mathcal{O}(\delta^{1/2})
 & \quad \text{if}\quad {\vartheta_{\sharp}}^{{\bf K}_\star}>0 ; \\
 \\
 +\frac{1}{2\pi}\frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}({\bm{\mathfrak{K}}}_2\wedge{\bm{\mathfrak{K}}}_1) + \mathcal{O}(\delta^{1/2})
  & \quad \text{if}\quad {\vartheta_{\sharp}}^{{\bf K}_\star}<0 .
  \end{array}\right.
\end{equation}
\end{theorem}

\begin{proof}[Proof of Theorem \ref{gp-vel}]
The edge state eigenvalue problem with ${k_{\parallel}}-$ pseudo-periodicity is given \eqref{dw_evp}-\eqref{localized}.
It is convenient to factor the Bloch phase with parallel quasi-momentum, ${k_{\parallel}}$, in order to replace \eqref{pseudo-per} by an equivalent periodic boundary condition. Thus, we set
$\Psi({\bf x}; {k_{\parallel}}) = e^{i\frac{{k_{\parallel}}}{2\pi} {\bm{\mathfrak{K}}}_1 \cdot {\bf x}} \psi({\bf x}; {k_{\parallel}})$, where $\psi({\bf x}+{\bm{\mathfrak{v}}}_1; {k_{\parallel}})=\psi({\bf x}; {k_{\parallel}})$. Further, we set
$
\frac{{k_{\parallel}}}{2\pi} {\bm{\mathfrak{K}}}_1 \equiv {\bf K}_\star + \frac{\xi}{2\pi} {\bm{\mathfrak{K}}}_1$, with  $ |\xi| \ll 1$ .
Since ${\bm{\mathfrak{K}}}_1\cdot{\bm{\mathfrak{v}}}_1=2\pi$, it follows that
\[{k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1+\xi.\]
Substituting
\[\Psi({\bf x}; {k_{\parallel}}) = e^{i\frac{{k_{\parallel}}}{2\pi} {\bm{\mathfrak{K}}}_1 \cdot {\bf x}} \psi({\bf x}; {k_{\parallel}}) = e^{i ({\bf K}_\star + \frac{\xi}{2\pi} {\bm{\mathfrak{K}}}_1) \cdot {\bf x}} \psi({\bf x}; \xi)\]
into \eqref{dw_evp}-\eqref{localized}, yields an equation for
 $\psi({\bf x};\xi)$:
\begin{align}
&\left[\mathcal{L}_{_{\rm dw}}^{(\delta)}({\bf K}_\star)+\frac{\xi}{2\pi} {\bm{\mathfrak{K}}}_1\cdot\left(
\mathscr{A}({\bf K}_\star)+\delta \mathscr{Q}^\delta({\bf K}_\star) \right) \right. \label{kpar_evp_psi1}  \\
&\qquad \left. + \left(\frac{\xi}{2\pi}\right)^2 {\bm{\mathfrak{K}}}_1 \cdot (A + \delta \eta^\delta B) {\bm{\mathfrak{K}}}_1 \right] \psi^\delta ({\bf x}; \xi)  \nonumber
= E^\delta (\xi) \psi^\delta ({\bf x}; \xi) , \\
&\psi({\bf x}+{\bm{\mathfrak{v}}}_1; \xi) = \psi({\bf x};\xi) , \qquad \psi({\bf x}; \xi) \to 0\ \ {\rm as}\ \ |{\bf x}\cdot{\bm{\mathfrak{K}}}_2|\to\infty . \label{kpar_evp_psi3}
\end{align}

Here, $\mathscr{A}$ and $\mathscr{Q}^\delta$ are first order operators given by:
\begin{align*}
\mathscr{A}&=\frac{1}{i}\nabla\cdot A + A\frac{1}{i}\nabla,\qquad
\mathscr{A}(\bb k)=e^{-i{\bf k}\cdot {\bf x}}\ \mathcal{A}\ e^{-i{\bf k}\cdot {\bf x}} ; \quad \text{see also \eqref{Ak_def}}; \nonumber\\
\mathscr{Q}^\delta&=\frac{1}{i}\nabla\eta^\delta B+\frac{1}{i}\eta^\delta B\nabla,\qquad \mathscr{Q}^\delta(\bb k)=e^{-i{\bf k}\cdot {\bf x}}\ \mathscr{Q}^\delta\ e^{-i{\bf k}\cdot {\bf x}}.
\end{align*}
Recall that $E^\delta(\xi=0)=E_{K_\star}^\delta$ is an edge state eigenvalue with corresponding $L^2(\Sigma)-$ (normalized) eigenstate $\psi^\delta ({\bf x}; 0)=e^{-i{\bf K}_\star\cdot{\bf x}} \Psi^\delta({\bf x};0)=e^{-i{\bf K}_\star\cdot{\bf x}} \Psi^\delta_{{\bf K}_\star}$.
Therefore, differentiating \eqref{kpar_evp_psi1} with respect to $\xi$ and setting $\xi=0$ yields
\begin{align}
\label{ham_xi_diff1}
&\left( \mathcal{L}_{_{\rm dw}}^{(\delta)}({\bf K}_\star) - E_{{\bf K}_\star}^\delta \right) \dot\psi^\delta ({\bf x}; 0)\nonumber\\
&\quad = \dot E^\delta(0) \psi^\delta ({\bf x}; 0)-\frac{1}{2\pi} {\bm{\mathfrak{K}}}_1\cdot
\left( \mathscr{A}({\bf K}_\star)+\delta \mathscr{Q}^\delta({\bf K}_\star) \right) \psi^\delta ({\bf x}; 0) ,
\end{align}
where we use the notation $\dot f(0) \equiv \partial_\xi f(\xi)|_{\xi=0}$.

Taking the $L^2(\Sigma)-$ inner product  of \eqref{ham_xi_diff1}
 with $\psi^\delta ({\bf x}; 0)$, we obtain the following expression for the group velocity at $\xi=0$ (${k_{\parallel}}={\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$):
 \begin{align}
 \frac{d E^\delta}{d\xi}(0)
 & = \frac{1}{2\pi} \inner{\psi^\delta ({\bf x}; 0) , {\bm{\mathfrak{K}}}_1\cdot
 \left( \mathscr{A}({\bf K}_\star)+\delta \mathscr{Q}^\delta({\bf K}_\star) \right) \psi^\delta ({\bf x}; 0)}_{_{L^2(\Sigma)}} \nonumber\\
 & = \frac{1}{2\pi} \inner{ \Psi^\delta_{{\bf K}_\star} , {\bm{\mathfrak{K}}}_1\cdot \mathscr{A}\ \Psi^\delta_{{\bf K}_\star} }_{_{L^2(\Sigma)}} + \mathcal{O}(\delta) .
 \label{E_xi_0}
\end{align}

Since $\delta$ is small, we may further simplify the expression for the group velocity in \eqref{E_xi_0} by substituting the expansion for  $\Psi_{{\bf K}_\star}^{\delta}({\bf x})$, displayed in \eqref{ldord} of Theorem \ref{edgemode_thm}.
N.B. In \eqref{ldord}, the positive sign corresponds to the case  $\vartheta_\sharp<0$ and the negative sign to the case $\vartheta_\sharp>0$.

We substitute \eqref{ldord} into \eqref{E_xi_0}. This inner product is a sum of terms of the form:
\begin{equation}
 \inner{ \Gamma_1({\bf x},\delta{\bf k}_2\cdot{\bf x}), \mathscr{A}({\bf x}) \Gamma_2({\bf x},\delta{\bf k}_2\cdot{\bf x})}_{_{L^2(\Sigma)}},
\end{equation}
where ${\bf x}\mapsto\Gamma_j({\bf x},\zeta)$ is $\Lambda_h-$ periodic and $H^2(\Omega_h)$  with values in $L^2(\mathbb{R}_\zeta)$, and $\mathscr{A}\equiv\mathscr{A}({\bf x})=\frac{1}{i}\nabla_{{\bf x}}\cdot A({\bf x}) + A({\bf x})\frac{1}{i}\nabla_{{\bf x}}$.
Therefore, we may exploit the multiple scale character to evaluate the inner product; see Lemma 7.9 of \cites{FLW-AnnalsPDE:16}. We obtain:

\begin{align}
 \frac{d E^\delta}{d\xi}(0)\ =\
 &\frac{1}{2\pi}\ {\bm{\mathfrak{K}}}_1 \cdot
 \inner{ \chi_\pm({\bm{\mathfrak{K}}}_2^\perp)\cdot\left(\Phi_1^{{\bf K}_\star},\Phi_2^{{\bf K}_\star}\right) ,
 \mathscr{A}\  \chi_\pm({\bm{\mathfrak{K}}}_2^\perp)\cdot\left(\Phi_1^{{\bf K}_\star},\Phi_2^{{\bf K}_\star}\right) }_{L^2(\Omega_{h})}\nonumber\\
 &\qquad + \mathcal{O}(\delta^{1/2}),
 \label{GV1}
 \end{align}

 %
 %
where we have also substituted the value of $\gamma$, displayed in \eqref{gamma-def}.
Again we emphasize that in \eqref{GV1}
the positive sign corresponds to the case
$\vartheta_\sharp<0$ and the negative sign to the case
 $\vartheta_\sharp>0$.
Finally, observe that
\begin{align}
& {\bm{\mathfrak{K}}}_1 \cdot \inner{ \chi_\pm({\bm{\mathfrak{K}}}_2^\perp)\cdot\left(\Phi_1^{{\bf K}_\star},\Phi_2^{{\bf K}_\star}\right) ,
  \mathscr{A}\ \chi_\pm({\bm{\mathfrak{K}}}_2^\perp)\cdot\left(\Phi_1^{{\bf K}_\star},\Phi_2^{{\bf K}_\star}\right) }_{L^2(\Omega_{h})} \nonumber\\
& \quad=\chi_\pm^\dagger({\bm{\mathfrak{K}}}_2^{\perp})\mathcal{M}_\mathscr{A}({\bm{\mathfrak{K}}}_1)\ \chi_\pm({\bm{\mathfrak{K}}}_2^{\perp}) \nonumber \\
&= \pm\frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}
 \Re\left(\overline{\mathfrak{z}({\bm{\mathfrak{K}}}_2^\perp)} \mathfrak{z}({\bm{\mathfrak{K}}}_2)\right) \nonumber \\
 &= \pm \frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}({\bm{\mathfrak{K}}}_1^{(2)}{\bm{\mathfrak{K}}}_2^{(1)}-{\bm{\mathfrak{K}}}_1^{(1)}{\bm{\mathfrak{K}}}_2^{(2)})
 = \left\{ \begin{array}{ll}
 +\frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}({\bm{\mathfrak{K}}}_2\wedge{\bm{\mathfrak{K}}}_1) & \quad \text{if} \quad \vartheta_\sharp<0; \\
 -\frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}({\bm{\mathfrak{K}}}_2\wedge{\bm{\mathfrak{K}}}_1) & \quad \text{if} \quad \vartheta_\sharp>0.
 \end{array} \right.
 \label{phipminner}
\end{align}
Here, ${\bm{\mathfrak{K}}}_2\wedge{\bm{\mathfrak{K}}}_1\equiv ({\bm{\mathfrak{K}}}_1^{(2)}{\bm{\mathfrak{K}}}_2^{(1)}-{\bm{\mathfrak{K}}}_1^{(1)}{\bm{\mathfrak{K}}}_2^{(2)})$.

Substituting \eqref{phipminner} into \eqref{GV1}, we obtain
\begin{equation}
\label{GV2}
 \frac{d E^\delta}{d\xi}(0) =
 \left\{
 \begin{array}{ll}
 - \frac{1}{2\pi}\frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}({\bm{\mathfrak{K}}}_2\wedge{\bm{\mathfrak{K}}}_1) + \mathcal{O}(\delta^{1/2})
 & \quad \text{if} \quad {\vartheta_{\sharp}}^{{\bf K}_\star}>0 ;\\
 \\
 + \frac{1}{2\pi}\frac{\upsilon_F}{|{\bm{\mathfrak{K}}}_2|}({\bm{\mathfrak{K}}}_2\wedge{\bm{\mathfrak{K}}}_1) + \mathcal{O}(\delta^{1/2})
  & \quad \text{if}\quad {\vartheta_{\sharp}}^{{\bf K}_\star}<0 .
  \end{array} \right.
\end{equation}
Equation \eqref{GV2_thm} and Theorem \ref{gp-vel} follow from \eqref{GV2} using the change of variables: $\xi={k_{\parallel}}-{\bf K}_\star\cdot{\bm{\mathfrak{v}}}_1$.
\end{proof}

\subsection{The direction of energy propagation along the edge}\label{direction}

In this section we apply Theorem \ref{gp-vel} to two particular situations where ${\mathcal{L}^B}=-\nabla\cdot B\nabla$ is $\mathcal{PC}$ anti-symmetric: $\mathcal{PC}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{PC}$.
\begin{enumerate}
\item  Preserving $\mathcal{C}-$ invariance and breaking $\mathcal{P}-$ invariance;
$[\mathcal{C},{\mathcal{L}^B}]=0$ and $\mathcal {P}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{P}$.
\item  Preserving $\mathcal{P}-$ invariance and breaking $\mathcal{C}-$ invariance;
$[\mathcal{P},{\mathcal{L}^B}]=0$ and $\mathcal {C}{\mathcal{L}^B}=-{\mathcal{L}^B} \mathcal{C}$.
\end{enumerate}

First note that because $\overline{L^2_{{\bf K},\bar\tau}}=L^2_{{\bf K}',\tau}$, we have that $\Phi_2^{{\bf K}}({\bf x})=\mathcal{C}[\Phi_1^{{\bf K}'}]({\bf x})$. Recall also that $\mathcal{PC}\Phi_2^{{\bf K}}=\Phi_1^{{\bf K}}$.

\emph{Case (1): $[\mathcal{C},{\mathcal{L}^B}]=0$ and  $\mathcal {P}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{P}$; bi-directional propagation along the edge.}
We claim that ${\vartheta_{\sharp}}^{{\bf K}'}=-{\vartheta_{\sharp}}^{{\bf K}}$.  Indeed, we have
\begin{align*}
{\vartheta_{\sharp}}^{{\bf K}'}
& = \inner{\Phi_1^{{\bf K}'},{\mathcal{L}^B}\Phi_1^{{\bf K}'}}_{L^2(\Omega)}
 = \inner{\mathcal{C}\Phi_2^{{\bf K}},{\mathcal{L}^B}\mathcal{C}\Phi_2^{{\bf K}}}_{L^2(\Omega)}
  = \inner{\mathcal{P}\Phi_1^{{\bf K}},{\mathcal{L}^B}\mathcal{P}\Phi_1^{{\bf K}}}_{L^2(\Omega)}\\
& = -\inner{\mathcal{P}\Phi_1^{{\bf K}},\mathcal{P}{\mathcal{L}^B}\Phi_1^{{\bf K}}}_{L^2(\Omega)}
 = -\inner{\Phi_1^{{\bf K}},{\mathcal{L}^B}\Phi_1^{{\bf K}}}_{L^2(\Omega)}
 = -{\vartheta_{\sharp}}^{{\bf K}} .
\end{align*}

In particular, ${\vartheta_{\sharp}}^{{\bf K}'}$ and
${\vartheta_{\sharp}}^{{\bf K}}$ have opposite sign. Moreover, by \eqref{GV2_thm}:

\begin{equation}
\lim_{\delta\to0} \partial_{k_{\parallel}} E^\delta ({\bf K}\cdot{\bm{\mathfrak{v}}}_1) = - \lim_{\delta\to0} \partial_{k_{\parallel}} E^\delta({\bf K}'\cdot{\bm{\mathfrak{v}}}_1) .
\label{E_xi_diff}
\end{equation}

The quantities $\partial_{k_{\parallel}} E^\delta({\bf K}\cdot{\bm{\mathfrak{v}}}_1)$ and $\partial_{k_{\parallel}} E^\delta({\bf K}'\cdot{\bm{\mathfrak{v}}}_1)$ are, respectively, the group velocities of wavepackets, constructed as a superposition of edge states,
 $(E^\delta({k_{\parallel}}),\Psi^\delta({\bf x};{k_{\parallel}}))$ for quasi-momentum ${k_{\parallel}}$ near ${\bf K}\cdot{\bm{\mathfrak{v}}}_1$ and ${\bf K}'\cdot{\bm{\mathfrak{v}}}_1$, respectively. We conclude that such packets with quasi-momentum centered at ${\kpar=\bK\cdot\vtilde_1}$ and ${\kpar=\bKp\cdot\vtilde_1}$, travel in opposite directions and at the same  speed; see Figure \ref{E_kpar_zz} (top panel) and the discussion in the Introduction (Section \ref{intro}).

\emph{Case 2: $[\mathcal{P},{\mathcal{L}^B}]=0$ and  $\mathcal {C}{\mathcal{L}^B}=-{\mathcal{L}^B}\mathcal{C}$; uni-directional propagation along the edge.}
In contrast to Case 1, here we claim that ${\vartheta_{\sharp}}^{{\bf K}'}={\vartheta_{\sharp}}^{{\bf K}}$.  Indeed, we have
\begin{align*}
{\vartheta_{\sharp}}^{{\bf K}'}
& = \inner{\Phi_1^{{\bf K}'},{\mathcal{L}^B}\Phi_1^{{\bf K}'}}_{L^2(\Omega)}
 = \inner{\mathcal{C}\Phi_2^{{\bf K}},{\mathcal{L}^B}\mathcal{C}\Phi_2^{{\bf K}}}_{L^2(\Omega)}
  = \inner{\mathcal{P}\Phi_1^{{\bf K}},{\mathcal{L}^B}\mathcal{P}\Phi_1^{{\bf K}}}_{L^2(\Omega)}\\
& = \inner{\mathcal{P}\Phi_1^{{\bf K}},\mathcal{P}{\mathcal{L}^B}\Phi_1^{{\bf K}}}_{L^2(\Omega)}
 = \inner{\Phi_1^{{\bf K}},{\mathcal{L}^B}\Phi_1^{{\bf K}}}_{L^2(\Omega)}
 = {\vartheta_{\sharp}}^{{\bf K}} .
\end{align*}

Arguing analogously to Case 1 we have  find for Case 2:
\begin{equation}
\lim_{\delta\to0} \partial_{k_{\parallel}} E^\delta ({\bf K}\cdot{\bm{\mathfrak{v}}}_1) = \lim_{\delta\to0} \partial_{k_{\parallel}} E^\delta({\bf K}'\cdot{\bm{\mathfrak{v}}}_1) .
\label{E_xi_diff2}
\end{equation}

Hence, the group velocities of wavepackets, constructed as a superposition of edge states,
 $(E^\delta({k_{\parallel}}),\Psi^\delta({\bf x};{k_{\parallel}}))$ for quasi-momentum ${k_{\parallel}}$ near ${\bf K}\cdot{\bm{\mathfrak{v}}}_1$, respectively, near ${\bf K}'\cdot{\bm{\mathfrak{v}}}_1$ are equal in magnitude and sign.  We conclude, in this case where $\mathcal{C}-$ invariance is broken,  that such packets with quasi-momentum centered at ${\kpar=\bK\cdot\vtilde_1}$ and ${\kpar=\bKp\cdot\vtilde_1}$, travel in the same direction and at the same  speed. See Figure \ref{E_kpar_zz} (bottom panel) and the discussion in the Introduction (Section \ref{intro}). Such unidirectional propagation of topologically protected edge states is well-known to be the hallmark of topological insulators \cites{hasan2010colloquium,ando2013topological,bernevig2013topological,ortmann2015topological}.

redpajama_set_name: RedPajamaArXiv