## Abstract

We developed a selection rule for Dirac-like points in two-dimensional dielectric photonic crystals. The rule is derived from a perturbation theory and states that a non-zero, mode-coupling integral between the degenerate Bloch states guarantees a Dirac-like point, regardless of the type of the degeneracy. In fact, the selection rule can also be determined from the symmetry of the Bloch states even without computing the integral. Thus, the existence of Dirac-like points can be quickly and conclusively predicted for various photonic crystals independent of wave polarization, lattice structure, and composition.

©2013 Optical Society of America

## 1. Introduction

Graphene, a single layer of carbon atoms in a honeycomb lattice, has attracted innumerable attention since a simple and effective means to produce it was found [1]. Many interesting transport behaviors of graphene originate from its unique band structures around the six corners of the hexagonal Brillouin zone, where two bands cross over each other linearly at the Fermi level. Recently, Dirac points were found at the Brillouin zone boundary in classical wave systems like photonic crystals (PCs) and acoustic crystals [2–13]. These Dirac points are associated with intriguing wave transport properties, such as the classical analogs of *Zitterbewegung* and quantum-Hall-effect edge states, as well as the extremal transmission phenomenon. Dirac points have also been achieved at the Brillouin zone center [14–20] in classical wave systems by employing accidental degeneracy. These Dirac points may also lead to interesting wave propagation behaviors, for example, cloaking [14–16] and super-anisotropy [16]. Although both the Dirac points at the zone boundary and the zone center are associated with linear bands, they do not share the same physics. To the best of our understanding, the former is induced by deterministic degeneracy while the latter is a result of accidental degeneracy. Mathematically, the enclosed path integration around the cone center gives a non-zero Berry phase for the former while this is not true for the latter [17]. Thus, to distinguish these two types Dirac points, the latter one is called “Dirac-like” points.

In some PCs, the presence of Dirac-like points at the Brillouin zone center can be predicted by an effective medium approach [14], where the wavelength in the host medium is larger than that in the embedded scatterers, and the Dirac-like points appear when the effective permittivity and permeability simultaneously become zero. Alternatively, Dirac-like points can be predicted by a tight-binding model [18–20], where the electromagnetic transfer integrals between the nearest neighboring sites have been taken into consideration. Very recently, we proposed a perturbation theory based on the $\overrightarrow{k}\cdot \overrightarrow{p}$ method to study the origin of the linear dispersions associated with Dirac/Dirac-like points in two-dimensional (2D) classical wave systems. Different from the effective medium description and the tight-binding model, this theory can be applied to more general cases, as it does not require a long wavelength or wave energy localized in (or near) the scatterers. The validity of the theory has been successfully demonstrated with examples of acoustic crystals, where two material parameters, i.e., the mass density and bulk modulus, can both be altered for the host medium and the scatterers [17].

Though the theory can be leveraged to investigate the properties of Dirac-like points, some fundamental questions, including the following, remain to be answered. How do the electromagnetic (EM) material parameters affect the occurrence of a Dirac-like point? What is the relationship between a Dirac-like point and the polarization of the wave? Does accidental degeneracy always lead to Dirac-like points? The answers to these questions are not only crucial to the engineering of PCs in utilizing Dirac-like points, but also of great scientific interest.

In this work, we adopted the perturbation theory developed earlier and conducted a comprehensive study of various PCs composed of dielectric media under the excitation of waves with different polarizations. The choice of dielectric media would make the fabrication of the PCs more feasible, but it might increase the difficulty of their engineering. In a dielectric medium, only the relative permittivity (${\epsilon}_{r}$) is an alterable parameter and the relative permeability (${\mu}_{r}$) is always set to be unity. Compared with previously studied acoustic systems [15, 17], half of the degrees of freedom are frozen. The other consequence of choosing dielectric media is that they respond distinctly to transverse electric (TE) polarized waves and to transverse magnetic (TM) polarized waves in two dimensions. The differences should reflect the influence of the material parameters and the wave polarization on the characteristics of Dirac-like cones.

After careful study, we found that the basic information on the dispersion relation around the degenerate states, where two or more Bloch states bear the same eigenfrequency, is contained in the corresponding Bloch states as well as their coupling integrals. If the integral vanishes, the dispersion relation is quadratic, whereas a non-zero integral gives rise to linear dispersions and also determines the slope. In fact, the integral itself does not need to be evaluated to determine if its value is zero. The spatial symmetry of the Bloch states is sufficient. We therefore establish a selection rule to define if Dirac-like cone dispersion is allowed at the degenerate states. Since the selection rule requires only the symmetry property of Bloch states, it is easy to implement and the validity does not depend on the composition, lattice structure, or wave polarization.

The remainder of this paper is organized as follows: a brief derivation of our perturbation theory and the selection rule is presented in Section 2. Results of various 2D PC systems together with TE and TM polarizations are shown in Section 3, where examples of both linear and nonlinear dispersion relations associated with accidental degeneracy are given and the underlying physics are explored. Conclusions are drawn in Section 4.

## 2. Theory

The 2D dielectric PCs studied in this work consist of a periodic array of dielectric (or air) cylinders with radii $R$ embedded in air (or dielectric media). The lattice constant is $a$, which is set to be the length unit. The dielectric media considered in this work have permittivity, ${\epsilon}_{d}>{\epsilon}_{0}$, and permeability, ${\mu}_{d}={\mu}_{0}$, respectively, where ${\epsilon}_{0}$ and ${\mu}_{0}$ denote the corresponding parameters of air. For such 2D systems, we start with the TM mode, where the electric field is parallel to the axis of the cylinders, i.e.,$\overrightarrow{E}=(0,0,{E}_{z})$. From Maxwell equations, it is not difficult to find that ${E}_{z}$ satisfies the following scalar wave equation:

Equation (1) in conjunction with the periodic boundary condition forms a standard Sturm-Liouville eigenvalue problem, where ${\omega}^{2}/{c}_{0}^{2}$ is the eigenvalue and the corresponding solution of ${E}_{z}$ is the eigenfunction. According to Bloch’s theorem, the solutions of ${E}_{z}$ can be expressed as Bloch wave functions, ${\Psi}_{n\overrightarrow{k}}$, taking the following form:

In our perturbation theory, the unperturbed eigenfunctions, ${\Psi}_{n{\overrightarrow{k}}_{0}}(\overrightarrow{r})={u}_{n{\overrightarrow{k}}_{0}}(\overrightarrow{r}){e}^{i{\overrightarrow{k}}_{\text{0}}\cdot \overrightarrow{r}}$, and eigenfrequencies, ${\omega}_{n0}$, of the Bloch states at the ${\overrightarrow{k}}_{0}$ point are assumed to be known [22]. These eigenfunctions can form a basis set to expand the Bloch functions at a point near ${\overrightarrow{k}}_{0}$, which is given as:

If we take all Bloch states at ${\overrightarrow{k}}_{0}$ as the basis set, we can always reproduce the band structures around ${\overrightarrow{k}}_{0}$ by solving Eq. (9). In practice, if we want to estimate the dispersion relation near ${\omega}_{n\overrightarrow{k}}^{}$, only those states with eigenfrequencies close to ${\omega}_{n\overrightarrow{k}}^{}$ need to be considered, as the states far away contribute to the higher order of $\Delta k=\left|\overrightarrow{k}-{\overrightarrow{k}}_{0}\right|$ [23]. In dielectric PCs, the lowest bands, which start from zero frequency and the $\Gamma $ point, are always linear. Beyond that, the dispersions around the $\Gamma $ point are in general nonlinear at finite frequencies. To achieve linear dispersions, the degenerate Bloch states are usually used. If the degeneracy is obtained by adjusting the materials and/or geometrical parameters, such as the permittivity and/or radius of the scatterer, it is called “accidental degeneracy”. On the contrary, if the existence of the degeneracy is independent of those parameters, it is called “deterministic degeneracy”. Here, since we are interested in the dispersion relations near the degenerate point, only those degenerate states at that point need to be included to determine the slopes of the dispersion relations. Thus, we need only to deal with the simplified linear system delineated by Eq. (9), which becomes analytically solvable. In this way, the analytic solution to Eq. (9) for small $\Delta k=\left|\overrightarrow{k}-{\overrightarrow{k}}_{0}\right|$ can be expressed as:

Although this formulation is presented here for TM waves, all the derivations and results can be easily carried over to TE waves by making the following substitutions: ${E}_{z}\Rightarrow {H}_{z}$, ${\mu}_{r}\Rightarrow {\epsilon}_{r}$, ${\epsilon}_{r}\Rightarrow {\mu}_{r}$. For this reason, we do not repeat the procedure here. While both cases share the same mathematical structure, their underlying physics is intrinsically distinct and gives rise to significantly different results, which are shown and discussed in detail in the next section.

The method presented here is valid in general, as it does not have lattice structure, frequency, or material parameter restrictions. The accuracy of the method is guaranteed by the fact that all the interactions between any two scatterers have been inherently taken into account because the unperturbed Bloch wave functions are rigorous solutions to EM wave equations in periodic systems. There are many ways to compute the Bloch functions at ${\overrightarrow{k}}_{0}$, and in this work we adopt COMSOL Multiphysics (a commercial package based on the finite element method) to obtain Bloch functions numerically. With the knowledge of the Bloch functions, we can easily evaluate the matrix elements, ${\overrightarrow{p}}_{lj}$ and ${q}_{lj}$, by performing the integrations defined in Eqs. (7) and (8).

## 3. Results and discussion

In this section, we demonstrate, through various examples, the validity of the theory proposed earlier and give a detailed explanation of the selection rule. Basically, we consider two types of PCs: One includes periodic arrays of dielectric cylinders embedded in air and the other is the reciprocal structure, i.e., it includes air boreholes in a dielectric background. We study for both TE and TM waves.

#### 3.1 Linear dispersions induced by accidental degeneracy

Figure 1(a)
shows the band structure of a 2D PC consisting of a triangular lattice of air boreholes with radius $R=0.4429a$ in a dielectric background with ${\epsilon}_{d}=12{\epsilon}_{0}$. It exhibits a triply degenerate point at the Brillouin zone (BZ) center at reduced frequency ${\omega}_{0}a/2\pi {c}_{0}=0.707$. An enlarged view of the bands around the $\Gamma $ point is plotted in Fig. 1(b), which clearly shows that two branches with linear dispersions intersect at the degenerate point at $\overrightarrow{k}=0$, through which point there is another flat band. This point and the associated linear dispersions form a Dirac-like cone, and the point is called a Dirac-like point. This type of dispersion relation is quite special compared with the general cases where quadratic behavior is always found in the dispersion relations near the $\Gamma $ point, such as the ones shown in Figs. 1(c) and 1(d) for smaller (*R = 0.44a*) and larger (*R = 0.445a*) air holes, respectively. It can be seen that there exist separately a single Bloch state and doubly degenerate Bloch states at the $\Gamma $ point. When the radii of the air holes are gradually increased, the single state and the doubly degenerate states would interchange their eigenfrequencies. In the meantime, the corresponding branches also change the signs of their curvatures. The critical point arises when the radii of the air holes are tuned to be *0.4429a*. The single state and the doubly degenerate states merge at the $\Gamma $ point and “accidental” degeneracy occurs. It is the interaction between them that results in the linear dispersion. The Bloch function for the single state is shown in Fig. 1(e), where different colors denote different values of ${H}_{z}$, and where arrows indicate the in-plane electric field, $\stackrel{\rightharpoonup}{E}$, whose magnitude is proportional to the length of the arrows. The doubly degenerate Bloch functions are plotted in Figs. 1(f) and 1(g) in a similar way. Figure 1(e) exhibits a typical pattern for a monopolar mode (corresponding to the *n = 0* angular momentum channel) hybridized with a non-negligible lattice effect. In comparison, the doubly degenerate states shown in Figs. 1(f) and 1(g) correspond to the dipolar mode (corresponding to the *n = 1* angular momentum channel) with lattice structure interference. This implies that the intensive interaction between the monopolar and dipolar modes gives rise to the linear dispersions in the BZ center.

The linear slopes of the band structures near the Dirac-like point can be accurately predicted by our theory. By using the Bloch wave functions shown in Figs. 1(e)-(g) to evaluate the integrals in Eqs. (7) and (8), calculating the matrix elements ${\overrightarrow{p}}_{lj}$ and ${q}_{lj}$, and then solving the determinant in Eq. (9), we finally get the roots for the dimensionless linear slopes, ${\gamma}_{\beta}\equiv \Delta {\omega}_{\beta}/({c}_{0}\Delta k)$. Here, the roots have three values ${\gamma}_{\beta}=0,\text{and}\pm 0.340$(see the Appendix for details). These results are plotted in Fig. 1(b) in solid red curves. Obviously, the first root corresponds to the flat branch and the remaining two represent the two linear bands crossing over each other. It can be proved that the slopes of the linear dispersions are independent of the directions of the wave vector, $\overrightarrow{k}$ (see Appendix), indicating that the linear bands are isotropic around the $\Gamma $ point such that they form two perfect cones and their common vertex is the Dirac-like point. In Figs. 1(a) and 1(b), we have also plotted the band structures obtained by COMSOL in blue dots. Excellent agreement between the finite element calculations and our perturbation theory is evident.

The linear dispersion relation near the degenerate point is a necessary condition for the existence of a Dirac-like point. As stated in the previous section, a non-vanishing ${\overrightarrow{p}}_{lj}$ is sufficient to guarantee the linear behavior of the dispersion relation. When the exact value of the slope is not important, a selection rule can be implemented quickly to predict if there exists a Dirac-like point. This selection rule is established purely by examining the spatial symmetry of the Bloch wave function and does not require evaluating the mode-coupling integrals. According to group theory [24], any Bloch function at the $\Gamma $ point of a triangular (square) lattice corresponds to an irreducible representation of the ${C}_{6v}$ (${C}_{4v}$) point group. In this regard, ${\overrightarrow{p}}_{lj}\ne 0$ is equivalent to the statement that the direct product of the irreducible representations of ${\Psi}_{l{\overrightarrow{k}}_{0}}$, $\overrightarrow{L}$ and ${\Psi}_{j{\overrightarrow{k}}_{0}}$ contains the fully symmetrical representation, ${A}_{1}$. That is to say, if the direct product contains ${A}_{1}$, ${\overrightarrow{p}}_{lj}$ is non-zero, and, therefore, the presence of a Dirac-like point is ensured. This is the “selection rule” for the Dirac-like point.

For the case shown in Fig. 1, the triply degenerate state is formed by accidental degeneracy of a single state and the doubly degenerate states. From the profile of the single Bloch state shown in Fig. 1(e), we know that its symmetry corresponds to the ${A}_{1}$ representation. Similarly, in Figs. 1(f) and 1(g), we find the doubly degenerate states correspond to the ${E}_{1}$ representation. The operator $\overrightarrow{L}$ transforms like a vector [24] and corresponds to the ${E}_{1}$ irreducible representation in the ${C}_{6v}$ point group. Since the direct product, ${A}_{1}\otimes {E}_{1}\otimes {E}_{1}$, contains ${A}_{1}$, according to the selection rule, ${\overrightarrow{p}}_{lj}$ is non-zero, which ensures the presence of a Dirac-like point at the Brillouin zone center. This example illustrates the easy implementation of the selection rule for the existence of Dirac-like points in 2D PCs by simply identifying the spatial symmetry of the Bloch states.

The above example unambiguously shows a Dirac-like cone whose presence can be conclusively predicted by the selection rule and whose linear slope can be evaluated by our perturbation theory. Another prominent example is the TM wave in a PC composed of a square array of dielectric cylinders studied in Ref [14]. A Dirac-like cone similar to the previous case was also found and was interpreted as a result of the double zero effective permittivity and permeability. That Dirac-like cone is also induced by accidental degeneracy of a single and doubly degenerate Bloch states, which transform like ${A}_{1}$ and *E* representations, respectively. Since the operator $\overrightarrow{L}$ corresponds to the $E$ irreducible representation in a ${C}_{4v}$ group, it follows that the direct product, ${A}_{1}\otimes E\otimes E$, contains ${A}_{1}$. Thus, according to the selection rule, there should be a Dirac-like point. Meanwhile, if we use our perturbation theory to calculate the dispersion relation, the results coincide with both the linear bands and the flat band shown in Ref [14]. The validity of our selection rule is again verified.

It is worth mentioning that though both of these cases exhibit similar dispersion behaviors at the $\Gamma $ point, their mechanisms are intrinsically different. For dielectric cylinder arrays, the wavelength inside the scatterer is much shorter than that in the air, and it is clear that the wave field is confined in the scatterer. Thus, the Dirac-like point is mainly determined by the scattering property of a single scatterer and can be well predicted by an effective medium theory. A Dirac-like cone can also be found in a triangular array if the filling ratio of the dielectric cylinder is fixed. However, this does not hold for its reciprocal structure because the wavelength in the air cylinder is larger than that in the dielectric host and the wave field is concentrated in the host, deeming the effective medium description inapplicable. The Dirac-like point is a result of the *interplay* between the scattering property of a single scatterer and the lattice structure. Thus, simply changing the lattice structure from a triangular array to a square array with a fixed filling ratio will not lead to a Dirac-like cone.

Although the mechanisms of the Dirac-like points mentioned earlier are distinct, their presence can be predicted by the selection rule as the selection rule only requires information about the symmetry of the eigenstates and has no restrictions on the wavelength. In principle, our selection rule is applicable to any case at any frequency, independent of lattice structure and material parameters.

#### 3.2 Nonlinear dispersions induced by accidental degeneracy

Accidental degeneracy does not necessarily lead to linear dispersion. According to our theory, when ${\overrightarrow{p}}_{lj}=0$, the dispersion relation near the degenerate point is expected to be quadratic, and, therefore, the degenerate point is not a Dirac-like point.

Figure 2(a) plots the band structure of the same photonic crystal as shown in Fig. 1 but the radius of the air holes is changed to $R=0.3244a$. A triply degenerate point is easily seen at the $\Gamma $ point, and an enlarged view of the band structure near this degenerate point is displayed in Fig. 2(b). Obviously, instead of two linear bands intersecting across a flat band, there are two nonlinear bands tangent to each other and a flat band at the degenerate point. In fact, this triply degenerate point is again the result of accidental degeneracy, which is demonstrated in Figs. 2(c) and 2(d) for larger ($R=0.33a$) and smaller ($R=0.31a$) radii of the air holes, respectively. It can be seen that when the radii deviate from $0.3244a$, the accidental degeneracy is broken and the triply degenerate states split into the single state and the doubly degenerate states. Interestingly, while Fig. 2(c) exhibits similar behavior to that shown in Fig. 1(c), Fig. 2(d) shows a very different characteristic from that of Fig. 1(d). Both the single-state branch and the double-state branch preserve the signs of their curvatures near the $\Gamma $ point after they interchange their eigenfrequencies, which indicates that the interactions between these two groups of states are weak compared with the case shown in Fig. 1. Thus, the quadratic behavior remains even though two bands touch. This conclusion is also supported by our perturbation theory as the calculation gives a zero mode-coupling integral, ${\overrightarrow{p}}_{lj}$, which means that the linear slopes of the bands around the BZ center should be identically zero.

The selection rule is also applicable for these nonlinear dispersions. It is easy to recognize that the single state, shown in Fig. 2(e), and the doubly degenerate states, shown in Figs. 2(f) and 2(g), correspond to ${A}_{1}$ and ${E}_{2}$ representations, respectively. The operator, $\overrightarrow{L}$, corresponds to the ${E}_{1}$ irreducible representation in ${C}_{6v}$. Since the direct product, ${A}_{1}\otimes {E}_{1}\otimes {E}_{2}$, does not contain ${A}_{1}$, the Dirac-like point cannot be found at the BZ center.

Since this case is the same as the one shown in Fig. 1 by simply changing the radius the scatterer, it is meaningful to compare these two cases. Interestingly, it can be seen that the single state for both cases represents the same behavior (i.e., both Fig. 1(e) and Fig. 2(e) exhibit monopolar modes) whereas the doubly degenerate states are absolutely different (dipolar modes shown in Fig. 1 and quadrupolar modes in Fig. 2). When the radius is small, the eigenfrequency of the monopolar mode is “sandwiched” by the eigenfrequencies of quadrupolar (lower) and dipolar (higher) modes. As the radius increases, the eigenfrequencies for these modes at the $\Gamma $ point all move upwards at different speeds, which makes accidental degeneracy possible. When the quadrupolar modes meet the monopolar mode, i.e., the situation shown in Fig. 2, nonlinear dispersions occur because of the accidental degeneracy. On the contrary, when the monopolar mode meets the dipolar modes, i.e., the case shown in Fig. 1, the strong interactions between them lead to linear dispersions around the BZ center. It can be seen that even for the same system, different accidental degeneracies result in different dispersion relations.

Another example that exhibits a similar non-liner behavior in the dispersion relation is a triangular array of dielectric cylinders with radius $R=0.4344a$ under the excitation of a TM wave, whose band structures are shown in Fig. 3 . We found again a triply degenerate point associated with two quadratic bands that are both tangent to a third flat band. According to the field distributions displayed in Figs. 3(c)-3(e), the degenerate Bloch states have the same irreducible representations as those of the previous example. Thus, a non-linear dispersion relation is expected according to the selection rule. On the other hand, it is easy to identify in the field distributions that the single state is a monopolar state and the doubly degenerate states are quadrupolar states. Interestingly, as discussed in the previous subsection, the same system but with a different radius would lead to a Dirac-like point, which is a result of the accidental degeneracy of the monopolar and dipolar modes.

All the examples demonstrated here show that the single and doubly degenerate states preserve their respective irreducible representations before or after the accidental degeneracy is achieved when we tune the radii of the scatterers. Thus, by examining the spatial symmetries of the Bloch states at the BZ center, we can predict which property (linear or nonlinear) of the dispersion relation there is if the accidental degeneracy of those Bloch states can be achieved.

## 4. Summary

In summary, we propose a selection rule for Dirac-like points in photonic crystals. It is based on a perturbation theory, in the spirit of the $\overrightarrow{k}\cdot \overrightarrow{p}$ method in electron systems, and requires only the information about the symmetry property of Bloch states at the degenerate point to determine whether the dispersion relation is linear in the vicinity of the point. We show that the direct product of the irreducible representations for corresponding degenerate Bloch states should contain the fully symmetric representation${A}_{1}$, or equivalently, a non-zero, mode-coupling integral of the degenerate Bloch states is necessary to ensure the existence of linear dispersions. The linear slopes associated with the Dirac-like point can be accurately reproduced by our perturbation formulae, which agree excellently with rigorous band structure calculations. We present examples of various dielectric photonic crystals with TM and TE waves.

## Appendix

In this appendix, we take the TE polarization as an example to show the calculation of the linear slopes around the BZ center. Actually, TM polarization can be dealt with in a similar way. The physical system under consideration is a 2D PC consisting of a triangular lattice of air boreholes in a ${\epsilon}_{d}=12{\epsilon}_{0}$ dielectric background, i.e., the system shown in Fig. 1. The three degenerate Bloch states shown in Figs. 1(e), 1(f) and 1(g) are denoted by ${\phi}_{1}$, ${\phi}_{2}$ and ${\phi}_{3}$, respectively. These Bloch states can be normalized in the following way:

where ${C}_{l}$ are the normalization coefficients that make these ${\tilde{\phi}}_{l}$s orthonormal to each other, i.e.,*a*is the lattice constant and serves as the length unit. According to the definition of ${\overrightarrow{p}}_{lj}$, we have,

After numerically calculating the integrals and retaining terms for the first order in $k$, Eq. (9) becomes

with## Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 11274120), the Fundamental Research Funds for the Central Universities (Grant No. 2012ZZ0077), and the KAUST Baseline Research Fund.

## References and links

**1. **A. H. Castro Neto, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. **81**(1), 109–162 (2009). [CrossRef]

**2. **F. D. M. Haldane and S. Raghu, “Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry,” Phys. Rev. Lett. **100**(1), 013904 (2008). [CrossRef] [PubMed]

**3. **S. Raghu and F. D. M. Haldane, “Analogs of quantum-Hall-effect edge states in photonic crystals,” Phys. Rev. A **78**(3), 033834 (2008). [CrossRef]

**4. **R. A. Sepkhanov, Y. B. Bazaliy, and C. W. J. Beenakker, “Extremal transmission at the Dirac point of a photonic band structure,” Phys. Rev. A **75**(6), 063813 (2007). [CrossRef]

**5. **R. A. Sepkhanov and C. W. J. Beenakker, “Numerical test of the theory of pseudo-diffusive transmission at the Dirac point of a photonic band structure,” Opt. Commun. **281**(20), 5267–5270 (2008). [CrossRef]

**6. **X. Zhang and Z. Liu, “Extremal transmission and beating effect of acoustic waves in two-dimensional sonic crystals,” Phys. Rev. Lett. **101**(26), 264303 (2008). [CrossRef] [PubMed]

**7. **X. Zhang, “Observing Zitterbewegung for photons near the Dirac point of a two-dimensional photonic crystal,” Phys. Rev. Lett. **100**(11), 113903 (2008). [CrossRef] [PubMed]

**8. **Q. Liang, Y. Yan, and J. Dong, “Zitterbewegung in the honeycomb photonic lattice,” Opt. Lett. **36**(13), 2513–2515 (2011). [CrossRef] [PubMed]

**9. **M. Diem, T. Koschny, and C. M. Soukoulis, “Transmission in the vicinity of the Dirac point in hexagonal photonic crystals,” Physica B **405**(14), 2990–2995 (2010). [CrossRef]

**10. **O. Peleg, G. Bartal, B. Freedman, O. Manela, M. Segev, and D. N. Christodoulides, “Conical diffraction and gap solitons in honeycomb photonic lattices,” Phys. Rev. Lett. **98**(10), 103901 (2007). [CrossRef] [PubMed]

**11. **T. Ochiai and M. Onoda, “Photonic analog of graphene model and its extension: Dirac cone, symmetry, and edge states,” Phys. Rev. B **80**(15), 155103 (2009). [CrossRef]

**12. **D. Torrent and J. Sánchez-Dehesa, “Acoustic analogue of graphene: observation of Dirac cones in acoustic surface waves,” Phys. Rev. Lett. **108**(17), 174301 (2012). [CrossRef] [PubMed]

**13. **L.-G. Wang, Z.-G. Wang, and S.-Y. Zhu, “Zitterbewegung of optical pulses near the Dirac point inside a negative-zero-positive index metamaterial,” Europhys. Lett. **86**(4), 47008 (2009). [CrossRef]

**14. **X. Huang, Y. Lai, Z. H. Hang, H. Zheng, and C. T. Chan, “Dirac cones induced by accidental degeneracy in photonic crystals and zero-refractive-index materials,” Nat. Mater. **10**(8), 582–586 (2011). [CrossRef] [PubMed]

**15. **F. M. Liu, X. Q. Huang, and C. T. Chan, “Dirac cones at $\overrightarrow{k}=0$in acoustic crystals and zero refractive index acoustic materials,” Appl. Phys. Lett. **100**(7), 071911 (2012). [CrossRef]

**16. **F. M. Liu, Y. Lai, X. Q. Huang, and C. T. Chan, “Dirac cones at$\overrightarrow{k}=0$in photonic crystals,” Phys. Rev. B **84**(22), 224113 (2011). [CrossRef]

**17. **J. Mei, Y. Wu, C. T. Chan, and Z.-Q. Zhang, “First-principles study of Dirac and Dirac-like cones in phononic and photonic crystals,” Phys. Rev. B **86**(3), 035141 (2012). [CrossRef]

**18. **K. Sakoda and H.-F. Zhou, “Role of structural electromagnetic resonances in a steerable left-handed antenna,” Opt. Express **18**(26), 27371–27386 (2010). [CrossRef] [PubMed]

**19. **K. Sakoda, “Dirac cone in two- and three-dimensional metamaterials,” Opt. Express **20**(4), 3898–3917 (2012). [CrossRef] [PubMed]

**20. **K. Sakoda, “Double Dirac cones in triangular-lattice metamaterials,” Opt. Express **20**(9), 9925–9939 (2012). [CrossRef] [PubMed]

**21. **J. Mathews and R. L. Walker, *Mathematical Methods of Physics*, 2nd ed. (Addison-Wesley, 1970).

**22. **P. M. Hui, W. M. Lee, and N. F. Johnson, “Theory of scalar wave propagation in periodic composites: a$\overrightarrow{k}\cdot \overrightarrow{p}$, ” Solid State Commun. **91**(1), 65–69 (1994). [CrossRef]

**23. **B. A. Foreman, “Theory of the effective Hamiltonian for degenerate bands in an electric field,” J. Phys. Condens. Matter **12**(34), R435–R461 (2000). [CrossRef]

**24. **M. S. Dresselhaus, G. Dresselhaus, and A. Jorio, *Group Theory: Application to the Physics of Condensed Matter* (Springer-Verlag, 2008).