Analytical study of mode degeneracy in non-Hermitian photonic crystals with TM-like polarization

TL;DR

This paper analytically investigates mode degeneracy in non-Hermitian photonic crystals with $C_{4v}$ symmetry, revealing conditions for Dirac-cone like dispersion, its properties, and topological transitions, using coupled-wave-theory.

Contribution

It provides an analytical framework for understanding accidental triple degeneracy and Dirac-cone like dispersion in non-Hermitian photonic crystals, highlighting their linearity, isotropy, and topological features.

Findings

Derivation of conditions for accidental triple degeneracy.

Identification of Dirac-cone like dispersion with good linearity and isotropy.

Observation of topological band transitions without Berry phase accumulation.

Abstract

We present an analytical study of the mode degeneracy in non-Hermitian photonic crystals (PC) with $C_{4v}$ symmetry, from the perspective of the coupled-wave-theory (CWT). The wave couplings and leakages within the non-Hermitian PCs are depicted, and the condition of accidental triple degeneracy is derived which leads to a Dirac-cone like dispersion. We prove that, similar to the real Dirac-cone, the Dirac-cone like band in non-Hermitian PC possesses good linearity and isotropy in the vicinity of the $\Gamma$ point. Moreover, some topological transit of the band may occur when the parameters adiabatically evolve. However, the Berry phase remains zero at the $\Gamma$ point, indicated the Dirac-like-cone dosen't obey Dirac equation and it is only Dirac-cone.

Tables

Table 1. Table 1: Structural parameters used in Fig. 3

Layer	Thickness	$\epsilon \left(𝐫\right)$
upper-clad (liquid)	$\infty$	${\epsilon }_q$
PC (Si3N4/liquid)	$0.6a$	$4.0804/{\epsilon }_q$
lower-clad (SiO2)	$\infty$	$2.1316$

Table 2. Table 2: Structural parameters used in Fig. 4

Layer	Thickness	$\epsilon \left(𝐫\right)$	$f$
upper-clad (liquid)	$\infty$	$4.00$	-
PC (Si/liquid)	$0.4a$	$12.11/4.00$	$0.28$
lower-clad (Si3N4)	$\infty$	$4.0804$	-

Equations

$$\nabla \times [\frac{1}{\epsilon (\mathbf{r})}\nabla \times \mathbf{H}(\mathbf{r})]=k^2\mathbf{H}(\mathbf{r}).$$

$$\mathcal{H}(n_y{H}_{x,mn}-m_x{H}_{y,mn})$$

$$\mathcal{H}(m_x{H}_{x,mn}+n_y{H}_{y,mn})+{\mathcal{H}}_{\delta }{H}_{z,mn}$$

$${\mathcal{H}}_z{H}_{z,mn}$$

$$\mathcal{H}$$

$${\mathcal{H}}_z$$

$${\mathcal{H}}_{\delta }$$

$$i{\beta }_0(m_x{H}_{x,mn}+n_y{H}_{y,mn})=\frac{\partial }{\partial z}{H}_{z,mn}$$

$$n_y{H}_{x,mn}-m_x{H}_{y,mn}$$

$$\mathcal{H}(m,n){\Theta }_{0,mn}(z)=0$$

$$k^2\mathbf{V}=(\mathbf{K}+\mathbf{C})\mathbf{V}=\mathbf{D}\mathbf{V}$$

$$\mathbf{C}={\mathbf{C}}_d+{\mathbf{C}}_{rad}+{\mathbf{C}}_h+{\mathbf{C}}_{TE\text{-}TM}$$

$$(k^2-k_0^2)\mathbf{V}=\mathbf{C}\mathbf{V}$$

$$\mathbf{C}$$

$${\text{TM}}_1:\mathbf{C}{\varphi }_1$$

$${\text{TM}}_2:\mathbf{C}{\varphi }_2$$

$${\text{TM}}_3:\mathbf{C}{\varphi }_3$$

$${\text{TM}}_4:\mathbf{C}{\varphi }_4$$

$${\varphi }_1$$

$${\varphi }_3$$

$$a-b+2c=a+b$$

$$⟹b=c$$

$$b_{1D}={\int }_{PC}{\kappa }_{2,0}[-\frac{{\partial }^2}{\partial z^2}+{\delta }_A\frac{\partial }{\partial z}-{\beta }_0^2]{\Theta }_0(z){\Theta }_0^{\ast }(z)dz$$

$$c_{1D}={\kappa }_{1,-1}{\beta }_0^2{\int }_{PC}{\Theta }_0(z){\Theta }_0^{\ast }(z)dz$$

$$k_{0,mn}=k_0-\frac{\xi }{\varsigma }\frac{{\beta }_0}{k_0}(m{k}_x+n{k}_y)$$

$$\displaystyle\bm{D}=\left(\begin{array}[]{cccc}a+b+2id&\frac{1}{\sqrt{2}}\mu k_{x}&-\frac{1}{\sqrt{2}}\mu k_{x}&0\\ \frac{1}{\sqrt{2}}\mu k_{x}&a-3b&0&\frac{1}{\sqrt{2}}\mu k_{y}\\ -\frac{1}{\sqrt{2}}\mu k_{x}&0&a+b&\frac{1}{\sqrt{2}}\mu k_{y}\\ 0&\frac{1}{\sqrt{2}}\mu k_{y}&\frac{1}{\sqrt{2}}\mu k_{y}&a+b+2id\end{array}\right)$$

$$\displaystyle\bm{D}_{r}=\left(\begin{array}[]{ccc}a+b+2id&-\frac{1}{\sqrt{2}}\mu k_{x}&0\\ -\frac{1}{\sqrt{2}}\mu k_{x}&a+b&\frac{1}{\sqrt{2}}\mu k_{y}\\ 0&\frac{1}{\sqrt{2}}\mu k_{y}&a+b+2id\end{array}\right)$$

$${\lambda }_{1,3}=(a+b+id)\pm \sqrt{\frac{{\mu }^2}{2}(k_x^2+k_y^2)-d^2}$$

$$|{\mathbf{k}}_{\parallel }|=\frac{\sqrt{2}d}{\mu }$$

$$\gamma =i\oint ⟨\Phi ({\mathbf{k}}_{\parallel })|{\nabla }_{{\mathbf{k}}_{\parallel }}|\Psi ({\mathbf{k}}_{\parallel })⟩d{\mathbf{k}}_{\parallel }+\eta$$

$${\nabla }_{{\mathbf{k}}_p}⟨{\Psi }_1^{\ast }({\mathbf{k}}_{\parallel })|{\Psi }_1({\mathbf{k}}_{\parallel })⟩$$

$$\displaystyle\bm{D}_{r}^{f}=\left(\begin{array}[]{ccc}a+b+2id&-\frac{1}{\sqrt{2}}\mu k_{x}&0\\ -\frac{1}{\sqrt{2}}\mu k_{x}&a-b+2c&\frac{1}{\sqrt{2}}\mu k_{y}\\ 0&\frac{1}{\sqrt{2}}\mu k_{y}&a+b+2id\end{array}\right)$$

$${\lambda }_{1,3}^f=(a+c+id)\pm \sqrt{\frac{{\mu }^2}{2}{k}_x^2+(S+id{)}^2}$$

$$\displaystyle u_{1}=\frac{1}{N_{-}}\left(\begin{array}[]{c}\frac{1}{\sqrt{2}}\mu k_{x}\\ S+id-M\\ 0\end{array}\right)u_{3}=\frac{1}{N_{+}}\left(\begin{array}[]{c}\frac{1}{\sqrt{2}}\mu k_{x}\\ S+id+M\\ 0\end{array}\right)$$

$$\frac{\mu k_x}{\sqrt{2}}=+d+r\cos \phi S=r\sin \phi$$

Full text

††preprint: APS/123-QED

I INTRODUCTION

Photonic crystals (PCs) have been attracted much attention in the past years owing to their great potential in manipulating light at wavelength scale. Various devices and applications, such as optical waveguide Vlasov et al. (2005); Mekis et al. (1996), filter Suh and Fan (2003); Noda et al. (2000); Song et al. (2003), cavity Akahane et al. (2003); Song et al. (2005), and laser Hirose et al. (2014) have been demonstrated. Conventionally, the light is vertically confined within the PC slab because of the total internal reflection at dielectric boundaries, and hence, such closed and lossless systems can be described by Hermitian operators. Recently, some extraordinary phenomena have been observed within the radiation continuum of PC slabs where light is allowed to escape and transport energy away Marinica et al. (2008); Plotnik et al. (2011); Lee et al. (2012); Hsu et al. (2013). Anomalously, narrow resonances can occur in the continuum. Such unique resonances have been interpreted as a photonic analogy of the bound states in the continuum (BICs) in quantum mechanics that were first proposed by von Neumann and Wigner von Neuman and Wigner (1929) and have been intensively studied in decades Raghu and Haldane (2008); Zhen et al. (2014); Haldane and Raghu (2008); Ni et al. (2016). Apparently, such PC systems are non-Hermitian since energy conservation no longer holds.

Owing to the analogy between quantum mechanics and electrodynamics, PCs possess photonic band that is similar to the electronic band of crystalline solids. It was recently found that, some accidental degeneracies can occur at the center of the Brillouin zone of a square Huang et al. (2011); Liu et al. (2011); Sakoda (2012) or triangular Sepkhanov et al. (2007); Zhang (2008); Diem et al. (2010); Peleg et al. (2007); Ochiai and Onoda (2009); Bittner et al. (2010) PC lattice, leading to the observation of Dirac-like cone Bittner et al. (2010) and various counter-intuitive transport properties such as zero-index metamaterials Huang et al. (2011); Liu et al. (2011); Li et al. (2015); Kita et al. (2016) and ring of EPs Zhen et al. (2015).

Previously, the 2D Hermitian PC systems were studied by using effective medium theory Diem et al. (2010); Peleg et al. (2007), $\vec{k}\cdot\vec{p}$ method Mei et al. (2012), and tight-binding approximation Sakoda (2012); Li et al. (2015). In these 2D PC systems, there is no radiation loss since the PCs are assumed to be infinitely thick. The triple degeneracy is found by continuously varying the radii of the dielectric pillar. However, radiation loss exists in the non-Hermitian PC slabs, leading to non-orthogonal eigenfunctions with complex eigenvalues. Therefore, we have to take into account more structural parameters to form the accidental degeneracy. Meanwhile, intriguing phenomena such as EPs have been experimentally demonstrated Zhen et al. (2015) which reveal some essential differences between Hermitian and non-Hermitian systems. As clarified by Mei et al. (2012) for Hermitian PC systems, the Dirac-like cone at Brillouin zone center with $C_{4v}$ symmetry has zero Berry phase, and hence it can not be mapped into the massless Dirac equation like the Dirac cone in Graphene. Consequently, such Dirac-like cone is not expected to give rise to extraordinary properties such as Zitterbewegung Zhang (2008); Liang et al. (2011a) and anti-localization. Yet, for non-Hermitian PC systems with Dirac-like cone, it is still unclear if such a claim holds.

In recent years, we have developed a comprehensive coupled-wave theory (CWT) framework capable of analytically modeling non-Hermitian PC systems Liang et al. (2011b); Peng et al. (2012, 2011); Liang et al. (2012, 2013); Yang et al. (2014). Since the radiation loss has been included in the analysis, CWT can be a promising tool to investigate the mode degeneracy in non-Hermitian systems. In this work, we extend the CWT framework to analytically study the complex band structure near the 2nd-order $\Gamma$ point of PC slabs with TM-like polarization and $C_{4v}$ symmetry. In contrast to our previous works, we include the TE-TM coupling terms which are crucial for the study of mode degeneracy within such non-Hermitian PC systems. Furthermore, we obtain the reduced coupling matrix as well as analytical wave functions, which provide analytical insights into the characteristics of band structure for non-Hermitian systems.

The remainder of this paper is organized as follows. In Section II, we derive the extended CWT formulation and the coupling matrix by including the important TE-TM coupling terms. In Section III, we first analytically derive the explicit condition of forming the triple degeneracy. Then, we obtain a reduced coupling matrix from the CWT equation to analytically study the linearity and isotropy of the band structure. It is found that, a ring of EPs inevitably appears around the 2nd-order $\Gamma$ point if the system is non-Hermitian. The topological properties of the triple degeneracy is discussed in Section IV. In Section V, we conclude with our findings.

II Theory and Formulation

In this section, we present the coupled-wave formulation for a two dimensional (2D) PC slab with periodicity in the $x$ and $y$ directions and multilayered structure in $z$ direction. As shown in Fig. 1(a), the PC layer consists of a square lattice of circular-shaped pillars with two different materials (permittivities: $\varepsilon_{1}$ and $\varepsilon_{q}$). The PC pillar has a finite thickness $2h$ with lattice constant $a$ and filling factor $f$. Under TM-like polarization, the $\mathbf{H}$ component follows:

[TABLE]

According to the Bloch’s theorem, For an off-$\Gamma$ reciprocal wave vector $\bm{k}_{\parallel}=k_{x}\hat{e_{x}}+k_{y}\hat{e_{y}}$, the components $H_{i}(z)$ can be expanded as $\sum_{m,n}H_{i,mn}(z)e(m,n)$ where $e(m,n)=\exp(-im_{x}\beta_{0}x-in_{y}\beta_{0}y),~{}i=x,y,z$, $\beta_{0}={2\pi}/{a}$ and $m_{x}=m+k_{x}/\beta_{0},n_{y}=n+k_{y}/\beta_{0}$. We also expand ${1}/{\varepsilon(\mathbf{r})}$ as $\kappa_{0}+\sum_{m\neq 0,n\neq 0}{\kappa}_{mn}e(m,n)$. For the PC layer, upper-clad and lower-clad, we denote $\kappa_{0}(z)$ as $\kappa_{a}=f/\varepsilon_{1}+(1-f)/\varepsilon_{q}$, $\kappa_{b}=1/\varepsilon_{q}$ and $\kappa_{c}=1/\varepsilon_{2}$, respectively.

For TM-like polarization that we focus on, the magnetic fields are given as $(H_{x},H_{y},H_{z})$. The $H_{z}$ component was neglected in our previous works because of the transverse nature of the TM mode. However, for a realistic and non-Hermitian PC slab, $H_{z}=0$ only rigorously holds at one particular cross section of the slab. The existence of $H_{z}$ component implies that the system contains TE components $\{E_{x},E_{y},H_{z}\}$, which will lead to some extra coupling paths that we denote as TE-TM couplings. Hence, the $H_{z}$ terms are included in the formulation to improve the accuracy (see Appendix A).

The set of coupled-wave equations that depicts the wave interactions within non-Hermitian PC slab is presented as (also see previous works Liang et al. (2011b); Peng et al. (2012); Yang et al. (2014)):

[TABLE]

where

[TABLE]

and $\delta_{K}=(\kappa_{b}-\kappa_{a})\delta(z-h)+(\kappa_{a}-\kappa_{c})\delta(z+h)$. The right hand side of Eqs.(2)-(4) depicts the coupling paths between individual Bloch waves, which can be classified into two groups: the terms $\Pi_{p},\Xi_{p},\Upsilon_{p}$ with subscript $p$ describe the in-plane couplings that are induced by the planar permittivity modulation (the only coupling mechanism for TE-like modes); the terms $\Pi_{s},\Xi_{s}$ with subscript $s$ represent the surface couplings owing to the independent boundary condition in TM-like polarization. By applying the boundary conditions appropriately (as discussed in Appendix A), such coupled-wave equations can be solved.

Submitting Eqs.(4) to (3), we obtain the transverse constraint:

[TABLE]

which is exactly the same as the divergence equation $\nabla\cdot\bm{H}=0$. For vertical homogeneous structure, or the case where the TE-TM coupling is sufficiently weak, Eq.(5) turns out to be a trivial form of $m_{x}H_{x,mn}+n_{y}H_{y,mn}=0$.

We choose a set of Bloch waves $\mathbf{V}$ as the basis to simplify Eqs.(2)-(4). At the 2nd-order $\Gamma$ point, the modal energy is dominated by the Bloch waves with in-plane wave vectors satisfying ${m^{2}+n^{2}}=1$, which are referred to as basic waves. Additionally, we denote the Bloch waves ${m^{2}+n^{2}}>1$ and ${m^{2}+n^{2}}=0$ as high-order waves and radiative wave, respectively; these waves are excited by the basic waves.

For an arbitrary wave vector $\bm{k}_{\parallel}$ near the 2nd-order $\Gamma$ point, the basic waves follow:

[TABLE]

where $\Theta_{0,mn},(m^{2}+n^{2}=1)$ are the vertical (out-of-plane) profile of the basic waves; and we denote $V_{mn}$ as $\{R_{x},S_{x},R_{y},S_{y}\}$ by giving $(m,n)\in\{(1,0),(-1,0),(0,1),(0,-1)\}$. Since the TE-TM coupling terms are relatively weak, we assume that the basic waves are almost in TM polarization, and the trivial transverse constraint $m_{x}H_{x,mn}+n_{y}H_{y,mn}=0$ still holds. We take the basic waves as the unperturbed basis and the couplings as perturbation. Hence, $\Theta_{0,mn}$ should satisfy

[TABLE]

By solving Eqs.(6)-(7) with trivial transverse constraint, the basic waves, including the vertical profile $\Theta_{0,mn}$ and the wavenumber $k_{0,mn}$, can be obtained. Besides, the high-order waves and radiative waves can be solved by using Green function method Liang et al. (2011b); Peng et al. (2012); Liang et al. (2013). Hence, we can rewrite the equations into an eigenvalue problem:

[TABLE]

where $\mathbf{V}=[R_{x}~{}S_{x}~{}R_{y}~{}S_{y}]^{T}$, and matrix $\bm{K}$ is in a diagonal form of $\{k_{0,10}^{2},k_{0,-10}^{2},k_{0,01}^{2},k_{0,0-1}^{2}\}$.

The coupling matrix $\bm{C}$ depicts the coupling strength between the individual basic waves, whereas the other orders of Bloch waves contribute as coupling paths. Here, the matrix $\bm{C}$ can be written as:

[TABLE]

where $\bm{C}_{d}$, $\bm{C}_{rad}$, $\bm{C}_{h}$ and $\bm{C}_{TE\text{-}TM}$ correspond to the direct couplings between basic waves [see Fig. 1(b)], the couplings via radiative waves, couplings via high-order waves, and the TE-TM couplings, respectively. Since the basic waves as well as high-order waves are confined within the PC slab, they do not contribute to energy leakage. As a result, the non-Hermiticity of matrix $\bm{C}$ only originates from $\bm{C}_{rad}$ and $\bm{C}_{TE\text{-}TM}$.

At the 2nd-order $\Gamma$ point ($\bm{k}_{\parallel}=0$), we have $(m_{x}^{2}+n_{y}^{2})=(m^{2}+n^{2})=1$, and hence, the basic waves share identical vertical profile $\Theta_{0}(z)$ and wavenumber $k_{0}$. Matrix $\bm{K}$ turns into a constant $k_{0}^{2}$ and Eq.(8) becomes:

[TABLE]

As proved in Appendix B, due to the $C_{4v}$ symmetry owned by the PC slab, the matrix $\bm{C}$ has a symmetric form as:

[TABLE]

where the matrix elements $a,b,c,d$ rely on a series of structural parameters, such as the filling factor, slab thickness, and the permittivity of the PC layer and upper/lower claddings. The CWT gives the explicit expressions of such matrix elements.

III Mode degeneracy and band structure

III.1 Condition of mode degeneracy at the $\Gamma$ point

From the perspective of CWT, the accidental degeneracy is equivalent to the degeneracy of eigenvectors of matrix $\bm{C}$. The non-Hermiticity of the PC system is described by matrix $\bm{B}$. In order to deal with the radiative waves, we should consider the eigenvalues of matrix $\bm{C}$ rather than its Hermitian part matrix $\bm{A}$.

Solving the eigenvalue problem Eq.(10) leads to four band-edge modes TM1∼4 at the $\Gamma$ point, whose eigenstates are presented in Eq.(III.1). Their field patterns are illustrated in Fig. 2.

[TABLE]

where

[TABLE]

Obviously, $\phi_{1}$ and $\phi_{4}$ are two degenerate eigenstates of matrix $\bm{B}$ with eigenvalue [math], corresponding to real frequencies with infinite Q factor. On the contrary, $\phi_{2}$ and $\phi_{3}$ correspond to a non-zero eigenvalue of $2d$, leading to complex mode frequencies, namely, finite lifetime.

According to Eq.(III.1), the condition of mode degeneracy becomes quite straightforward. Since the two leaky modes TM2 and TM3 are degenerate by nature owing to their in-plane symmetry, accidental triple degeneracy can be realized by tuning the structural parameters to make them degenerate with a single high-Q mode, namely, mode TM1. From the perspective of CWT, the real parts of their eigenvalues should be equal to each other, which gives the condition of mode degeneracy at the 2nd-order $\Gamma$ point:

[TABLE]

The detailed expressions of $b$ and $c$ can be found in Appendix B. For a low-index-contrast structure, $b$ and $c$ can be simplified by only keeping the strongest coupling term $\bm{C}_{1D}$ but neglecting $\bm{C}_{2D}$, $\bm{C}_{rad}$, $\bm{C}_{TE\text{-}TM}$, as:

[TABLE]

where $\kappa_{m,n}=2f(1/\varepsilon_{1}-1/\varepsilon_{q})J_{1}(\zeta)/\zeta$ with $\zeta=2\sqrt{f\pi(m^{2}+n^{2})}$. Apparently, $b$ and $c$ not only depend on the coupling coefficients $\kappa_{mn}\{f,\varepsilon_{1},\varepsilon_{q}\}$, but also rely on the vertical profile $\Theta(z)\{h,f,\varepsilon_{1},\varepsilon_{2},\varepsilon_{q}\}$. For infinitely thick PCs stated in other works, the accidental degeneracy is realized by continuously varying the radius of the hole or pillar. However, the realistic PC slabs actually provide extra freedom to achieve $b=c$.

As presented in Fig. 1(b), the element $b_{1D}$ depicts the feedback between basic waves $R_{x}$ and $S_{x}$ (also $R_{y}$ and $S_{y}$), whereas the element $c_{1D}$ presents the coupling between the basic waves propagating in orthogonal directions, for instance, $R_{x}$ and $R_{y}$. For the TE modes, the polarization of basic waves in $x$ and $y$ directions are normal to each other, which forbids the direct couplings between them. However, for the TM modes, the electric fields of the basic waves are always perpendicular to the transversal PC plane, and hence, can be coupled directly regardless of the propagating directions, leading to a non-zero value of $c$. As a result, the condition $b=c$ can be more easily fulfilled for the TM modes.

III.2 Linearity and isotropy of the Dirac-cone like band near the $\Gamma$ point

From the explicit expression of the triple degeneracy,we can analytically investigate the linearity and isotropy of the band near the 2nd-order $\Gamma$ point. Here we focus on a realistic structure shown in Fig. 1 (a) with structural parameters listed in Table I. The parameter set $(f,\varepsilon_{q})$ is tuned to realize this accidental triple degeneracy.

Fig. 3(a) illustrates the elements $b$ and $c$ in the parameter space $\{f,\varepsilon_{q}\}$. The curvature $(b-c)$ intersects with the zero plane, which gives the line $b-c=0$, on which every point represents a set of $(f,\varepsilon_{q})$ for realizing triple degeneracy. Here we focus on the point $O(f=0.35,\varepsilon_{q}=2.15)$. Fig. 3(b) presents the band structure of point $O$ in the $\bm{k}_{\parallel}$-space, and the details in $X\text{-}\Gamma\text{-}M$ direction are illustrated in Fig. 3(c) and (d). The CWT results agree well with the finite-element method (FEM, COMSOL Multiphysics) simulations. The band structure clearly shows an accidental triple degeneracy at the 2nd-order $\Gamma$ point, which is formed by the degenerate leaky modes TM2, TM3 and high-Q mode TM1.

We investigate the band structure by assuming a small off-$\Gamma$ wave vector $\bm{k}_{\parallel}$. For a low-index-contrast Si3N4 structure, the coupling matrix $\bm{C}$ is almost invariant around the 2nd-order $\Gamma$ point, and hence, matrix $\bm{K}$ dominates the band gap at off-$\Gamma$ points.

It’s easy to prove that (see Appendix C):

[TABLE]

where $\xi$ and $\varsigma$ are the coefficients related to structural parameters. Therefore, the matrix $\bm{K}$ at an off-$\Gamma$ point is linear with respect to $k_{x}$ and $k_{y}$. Taking $\bm{K}$ as perturbation and applying the condition of triple degeneracy $b=c$ in matrix $\bm{C}$, the coupling matrix $\bm{D}$ under basis $\{\phi_{3},\phi_{4},\phi_{1},\phi_{2}\}$ becomes:

[TABLE]

where $\mu=-2\xi\beta_{0}/\varsigma$.

The matrix above clearly depicts the couplingss brought by the off-$\Gamma$ wave vector $\bm{k}_{\parallel}$. As modes TM1,2,3 are triply-degenerate and TM4 is far above them, we reduce the matrix by neglecting the coupling terms with mode TM4 as:

[TABLE]

One eigenvalue of $\bm{D}_{r}$ is $\lambda_{2}=a+b+i2d$, corresponding to the flat band whose real part remains $a+b$ and imaginary part remains $2d$. Despite this trivial solution, the other two eigenvalues are

[TABLE]

Clearly $\lambda_{1,3}$ only depend on $|\bm{k}_{\parallel}|=\sqrt{k_{x}^{2}+k_{y}^{2}}$ regardless of the specific direction, which means the mode frequencies are isotropic in arbitrary direction.

Eq.(52) indicates that, there exists an off-$\Gamma$ ring satisfying $\sqrt{k_{x}^{2}+k_{y}^{2}}=\sqrt{2}d/\mu$ on which $\lambda_{1}$ and $\lambda_{3}$ are fully degenerate. These off-$\Gamma$ degenerate points are EPs that we will elaborate in the following sections. It is noteworthy that, as long as the radiation exists, the ring of EPs would always exist around the 2nd-order $\Gamma$ point, making the band discontinuous and not cone-like. In other words, there’s no rigorous cone dispersion in non-Hermitian photonic crystal slab with $C_{4v}$ symmetry.

However, if $d$ is sufficiently small, the EPs ring would be very close to 2nd-order $\Gamma$ point, and the overall band would still look like a cone. For instance, the band illustrated in Fig. 3(c) clearly exhibits a cone-like dispersion. In this case, back to Eq.(52), $d$ could be neglected and $\lambda_{1,3}$ are both linear with the same slope coefficient of $\mu/\sqrt{2}$, which indicates that the top and bottom halves of the cone share the same linearity.

Some interesting details can be noticed from Fig. 3(c). For example, the linearity of the band in $X\text{-}\Gamma$ direction remarkably degrades when deviating from 2nd-order $\Gamma$ point. Moreover, in $\Gamma\text{-}M$ direction, the leaky mode TM2 no longer remains flat but behaves as ripple-like dispersion which could easily be observed in Fig. 3(b). We believe this band distortion is due to the influence of mode TM4. Such mode has been neglected in previous studies of Hermitian systems, and also been dropped for simplicity in our reduced matrix $\bm{D}_{r}$. The couplings brought by TM4 are proportional to $k_{x}$ or $k_{y}$, and hence, significant band distortion will occur at large off-$\Gamma$ wave vectors. On the other hand, for a high-index-contrast structure, or strong optical confinement, the linear region of the cone-like band will be relatively small since the band could be twisted by the strong coupling strength.

III.3 Ring of exceptional points

As mentioned, the EPs ring where $\lambda_{1}$ and $\lambda_{3}$ are fully degenerate is inevitable in non-Hermitian PCs. As stated above, Eq.(52) indicates the EPs ring appearing at the off-$\Gamma$ wave vectors satisfying:

[TABLE]

Obviously, there is an infinite number of EPs on the ring which separates the whole plane into two sections. Inside this ring, the real frequencies of $\lambda_{1}$ and $\lambda_{3}$ are degenerate as $a+b$, while outside this ring, the imaginary parts are degenerate as a value of $d$. The ring is a collection of singular points and its size depends on the strength of radiation loss $d$ and coefficient $\mu$.

Here we present a realistic structure with relatively strong radiation loss (Q${}_{v}\approx 110$), in order to demonstrate a distinct and clear EPs ring. The structural parameters are listed in Table II.

Fig. 4 (a) and (b) illustrate the complex band structure in $\bm{k}_{\parallel}$-space. The EPs ring can be clear observed from both the real and imaginary parts of the band. The details of the band in $X\text{-}\Gamma\text{-}M$ direction are presented in Fig. 4 (c) and (d). The CWT results are in good agreement with FEM simulations, which confirms the validity of the CWT model.

Eq.(52) is derived from the reduced coupling matrix $\bm{D}_{r}$ neglecting the couplings with mode TM4. As a result, the EPs ring is exactly a circle, in other words, isotropic in the $\bm{k}_{\parallel}$-space. From Fig. 4 (c) and (d), we find the results given by $\bm{D}_{r}$ agree well with the CWT results and FEM results. However, similar to the distortion of band in Fig. 3, the couplings between TM4 and $\text{TM}_{1,2,3}$ also degrade the isotropy of the ring.

IV Topological properties of the band with triple degeneracy

In this section, we investigate the topological properties of band with triple degeneracy from reduced coupling matrix $\bm{D}_{r}$ in several specific parameter spaces.

IV.1 Zero Berry Phase of band structure in $\bm{k}_{\parallel}$-space

As we stated above, in non-Hermitian case, there always exists a EPs ring around the 2nd-order $\Gamma$ point. On the ring, the coupling matrix becomes defective. Due to the position with respect to the EPs ring, there are three typical different loops in $\bm{k}_{\parallel}$-space shown in Fig. 5.

First we consider type I and type II loop. It’s clear that there are no singular points on these loops. We revisit the Berry phase that may be accumulated along these loops, since it is an important indicator to distinguish whether the cone-like band follows the Dirac equation. The Berry phase in $\bm{k}_{\parallel}$-space for a non-Hermitian system is defined as Liang and Huang (2013)

[TABLE]

Where $e^{i\eta(T)}=\left\langle\Phi(\bm{k}_{\parallel}(0))|\Psi(\bm{k}_{\parallel}(T))\right\rangle$ and $T$ refers to the end of loop. The left vector $\Phi$ is the eigenvector in the dual space. From the perspective of CWT, the explicit wave functions of the degenerate modes can be chosen as the eigenvectors of coupling matrix. Here we calculate the $\Psi$ from reduced coupling matrix $\bm{D}_{r}$ to simplify the evaluation. This is valid because we only concern a small region around the 2nd-order $\Gamma$ point where the coupling terms with mode TM4 is negligible.

Here we take $|\Psi_{1}\rangle$ with eigenvalue $\lambda_{1}$ in Eq.(52) as an example, whose left vector $|\Phi_{1}\rangle$ satisfies $\lambda_{1}^{*}|\Phi_{1}\rangle=\bm{D}_{r}^{\dagger}|\Phi_{1}\rangle$. Obviously $\bm{D}_{r}$ is a symmetric non-Hermitian matrix, so we have $\bm{D}_{r}^{\dagger}=\bm{D}_{r}^{*}$ and the left vector can be chosen as $|\Phi_{1}\rangle=|\Psi_{1}^{*}\rangle$. Differentiation of the orthogonalization relation $\left\langle\Psi_{1}^{*}(\bm{k}_{\parallel})|\Psi_{1}(\bm{k}_{\parallel})\right\rangle=1$ yields

[TABLE]

So the Berry connection $\left\langle\Phi_{1}(\bm{k}_{\parallel})|\nabla_{\bm{k}_{\parallel}}|\Psi_{1}(\bm{k}_{\parallel})\right\rangle=0$, and the Berry phase is solely dominated by $\eta$ in Eq.(54). Clearly, there’s no modal interchange along the Type I and Type II loops and no phase change of the instantaneous eigenvector $\Psi_{1}(\bm{k}_{\parallel})$. So phase $\eta$ vanishes and Berry phase remains zero. For eigenvector $|\Psi_{3}\rangle$ with eigenvalue $\lambda_{3}$, the result is the same.

Zero Berry phase indicates that, encircling around the $\Gamma$ point in the $\bm{k}_{\parallel}$-space dose not accumulate any extra geometric phase, no matter the loop is inside the EP ring or outside the EP ring. Therefore, the band structure at the $\Gamma$ point under $C_{4v}$ symmetry can not be described by the Dirac equation. In other words, the cone-like dispersion we observe at 2nd-order $\Gamma$ point with $C_{4v}$ symmetry is actually not a Dirac cone but Dirac-like cone. Compared to Dirac cone in Graphene, it does not give rise to some unique properties like anti-localization when disorder exists. Nevertheless, the Dirac-like cone still possesses good linearity and isotropy near the 2nd-order $\Gamma$ point with small radiation.

IV.2 Modal interchange and non-zero Berry phase in parameter space $\{k_{x},S\}$

Though Berry phase remains zero in $\bm{k}_{\parallel}$-space, non-zero Berry phase may be realized in other parameter spaces. We consider a loop in space $\{k_{\parallel},S\}$ where $S$ refers to any structural parameter and $k_{\parallel}$ refers to a specific direction in $\bm{k}_{\parallel}$-space since the band structure is isotropic in the vicinity of 2nd-order $\Gamma$ point. Here we take $\{k_{y}=0,S=f\}$ as an example. When $f$ varies, the condition of triple degeneracy $b=c$ no longer holds and we denote $S=b-c$, which is equivalent to parameter $f$.

Space $\{k_{x},S\}$ is shown in Fig. 5 on which loop is denoted as Type III. Clearly we find space $\{k_{x},S\}$ contains two EPs at $(\pm\sqrt{2}d/\mu,0)$.

When $f$ varies, the reduced coupling matrix becomes

[TABLE]

the non-trivial eigenvalues of $\bm{D}_{r}^{f}$ with $k_{y}=0$ are

[TABLE]

Fig. 6 gives the Riemann surfaces of $\lambda^{f}_{1,3}$ in parameters space $\{k_{x},S\}$, from which two EPs at $(\pm\sqrt{2}d/\mu,0)$ can also be clearly observed. In fact, the two EPs are the branch points at which the two sheets of Riemann surfaces intersect. On the line $S=0$, the real parts are degenerate for $|k_{x}|<\sqrt{2}d/\mu$, and the imaginary parts are degenerate for $|k_{x}|>\sqrt{2}d/\mu$.

The loop evolving in the parameter space $\{k_{x},S\}$ can be topologically different according to how the EPs are encircled. When no EP is encircled, the eigenvalues simply return to their initial values. When a single EP is encircled, both the real and imaginary parts of the eigenvalues cross to the other sheet after finishing the loop, as a result, the eigenvalues interchange with each other. After two cycles the eigenvalues return to theirselves. Further, when both two EPs are encircled, the real parts of eigenvalues stay on the same sheet while the imaginary parts cross to the other sheet but will then cross back. After finishing the loop, the eigenvalues return to their initial values.

The eigenvectors corresponding to $\lambda_{1,3}^{f}$ of $\bm{D}_{r}^{f}$ are

[TABLE]

where $N_{\pm}=\sqrt{\mu^{2}k_{x}^{2}/2+(S+id\pm M)}=\sqrt{2M^{2}\pm 2M(S+id)}$ and $M=\sqrt{\mu^{2}k_{x}^{2}/2+(S+id)^{2}}$. Since coupling matrix $\bm{D}_{r}^{f}$ is symmetric, as we stated in last section, the left vectors could be chosen as $|v_{1,3}\rangle=|u^{*}_{1,3}\rangle$. In order to evaluate the Berry phase, we analyze the behavior of the eigenvectors $u_{1},u_{3}$ in space $\{k_{x},S\}$ by assuming a small circle of radius $r$ around a single EP in Fig. 6 Keck et al. (2003), chosen here as $(\mu k_{x}/\sqrt{2}=+d,S=0)$:

[TABLE]

With the approximation $r\ll d$, the two eigenvectors $u_{1,3}$ could be represented in the form:

[TABLE]

where

[TABLE]

and further we have

[TABLE]

Eq.(78) indicates the behavior of eigenvectors:

[TABLE]

Obviously, the eigenvectors are interchanged after a single circle encircling one EP and one of them gain extra phase $\pi$. After second cycle both eigenvectors are reconstructed with extra phase $\pi$. Further, we can infer that for a single loop encircling both two EPs, both eigenvectors accumulate a extra phase $\pi$ too. As the coupling matrix $\bm{D}_{r}^{f}$ is symmetric and Berry phase is dominated by phase $\eta$ in Eq.(54), the geometric phase is

[TABLE]

both for two circle encircling a single EP and one circle encircling both EPs.

What we should point out is that the approximation $r\ll d$ we take in considerations above is only used to investigate the behavior of eigenvectors in a simple way. Virtually because of the continuity of the band structure except EPs in Fig.6, the conclusion Eq.(80) are exact regardless of the specific loop.

IV.3 Modal exchange in parameter space $\{f,\varepsilon_{q}\}$

Further, we consider a structural parameter space, i.e space $\{f,\varepsilon_{q}\}$ shown in Fig. 7 (a). In this parameter space exists line $b=c$ corresponding to condition of triple degeneracy. It is noteworthy that, because of the uncertainty in fabricating realistic devices, the working points in this space would always slightly deviate from the exact degenerate condition $b=c$. Due to the uncertainty of the measurement, it is also difficult to distinguish whether the rigorous degenerate condition $b=c$ has been met.

When the parameters slightly deviate from $b=c$, the sign of $b-c$ determines which mode (mode TM1 or TM3) lying on the top half of the cone. More specifically, for the points in parameter space with $b>c$, mode TM1 forms the bottom half; while for the points with $b<c$, mode TM3 turns to be the bottom half.

Such phenomena is illustrated in Fig. 7 (b). When system evolves adiabatically along a loop encircling a rigorous degenerate point $O$, the system would cross the $b-c$ plane twice and return to its initial state. During this loop encircling, we can observe that the two bands forming the up/bottom halves of the cone switch twice, and eventually come back to the initial state. We refer to this evolution as a topological exchange of the modes.

V CONCLUSION

In this work, we present an analytical study of the mode degeneracy with TM-like polarization in a non-Hermitian photonic crystal slab with $C_{4v}$ symmetry. The CWT framework provides an analytical perspective to comprehensively understand the physics of a non-Hermitian photonic crystal system. We introduce the TE-TM couplings to improve the accuracy and extend CWT framework to depict the analytical wave functions and complex band structure through reduced coupling matrix.

From the CWT framework, the elements of the coupling matrix are evaluated , and hence, the condition of realizing accidental triple degeneracy is derived. Comparing to the Hermitian PCs that extend infinitely along $z$ direction, the realistic non-Hermitian PCs provide extra freedoms in tuning the structural and material parameters to fulfill the accidental degeneracy condition. The ring of EPs appears around 2nd-order $\Gamma$ point inevitably with triple degeneracy in non-Hermitian system. The band structure exhibits a cone-like dispersion with weak radiation loss. When the radiation loss become strong, the ring of EPs becomes distinct and easy to observe.

Furthermore, the linearity and isotropy, as well as the topological properties on the triple degeneracy, have been studied with the reduced coupling matrix. The analytical study reveals that, the cone-like band in a non-Hermitian PC with weak radiation loss possesses good linearity and isotropy in the vicinity of the $\Gamma$ point, which is similar to a real Dirac cone. The ring of EPs with strong radiation loss also owns good isotropy.

Since the CWT gives a set of explicit wave functions of the degenerate modes, we can use them to calculate the Berry phase directly. The Berry phase remains zero at the center of the Brillouin zone, even the system is non-Hermitian. It indicates that the band structure induced by the accidental triple degeneracy with weak radiation loss is only Dirac-like cone. Different from Graphene with real Dirac-cone, the PC with Dirac-like cone doesn’t obey Dirac equation, and doesn’t give rise to some unique properties like anti-localization against disorder, too. However, analytical study reveals that a loop in parameter space $\{k_{x},S\}$ continuously and adiabatically encircling EPs gives rise to a modal interchange phenomenon and results in a non-zero Berry phase $\pi$, where S refers to a structural parameter like $f$. Moreover, a modal exchange of the band will occur when the parameters evolve as a loop encircling a rigorous degenerate point in the structural parameter space $\{f,\varepsilon_{q}\}$.

Acknowledgments

This work was partly supported by the National Natural Science Foundation of China under Grant 61320106001, 61575002 and the State Key Laboratory of Advanced Optical Communication Systems and Networks, China. Yong Liang is supported by the ETH Zurich Postdoctoral Fellowship Program (No. FEL-27 14-2) cofounded by the Marie Curie Actions for People COFUND Program.

Appendix A the boundary conditions

The representation in details of in-plane coupling items in Eq.(2)-(4) are

[TABLE]

And the surface coupling items have the representation in details as:

[TABLE]

where $\delta_{h}=\delta(z-h)-\delta(z+h)$, $l_{1}=m_{x}n^{\prime}_{y}-n_{y}m^{\prime}_{x}$ and $l_{2}=m_{x}m^{\prime}_{x}+n_{y}n^{\prime}_{y}$.

The boundary conditions can be derived from the integral on Eq.(2) and (3) over the upper and down boundaries $(h-\delta,h+\delta)$ and $(-h-\delta,-h+\delta)$ with $\delta\to 0$. Since the the Dirac function $\delta_{h}$ terms in $\Pi_{s},\Xi_{s}$ has been treated as surface coupling, the discontinuity only comes from the $\delta_{K}$ terms in operator $\mathcal{H}$ and $\mathcal{H}_{\delta}$. We have:

[TABLE]

[TABLE]

where the terms of $(\kappa_{b}-\kappa_{a})\{\kappa_{b}\beta_{0}^{2}(m_{x}^{2}+n_{y}^{2})/(k^{2}-\kappa_{b}\beta_{0}^{2}(m_{x}^{2}+n_{y}^{2}))\}$ and $(\kappa_{c}-\kappa_{a})\{\kappa_{c}\beta_{0}^{2}(m_{x}^{2}+n_{y}^{2})/(k^{2}-\kappa_{c}\beta_{0}^{2}(m_{x}^{2}+n_{y}^{2}))\}$ come from the $\delta_{K}H_{z,mn}$ terms in $\mathcal{H}_{\delta}$ operator, which is a result of the TE-TM coupling. Because there is no Dirac function in $\mathcal{H}_{z}$ operator, the boundary conditions for Eq.(4) is continuous.

Coupling Eqs.(2)-(4) include TE-TM couplings introduced by component $H_{z,mn}$. To highlight the importance of such TE-TM couplings, we show in Fig. 8 the band structure evaluated by CWT without TE-TM couplings. The structural parameters are the same as those used for Fig. 4, as listed in Table II. By comparing to Fig. 4 (which includes the TE-TM couplings), we find that the CWT results without TE-TM couplings cannot reproduce the FEM results and describe the feature of the EPs. Therefore, including the TE-TM couplings is critical for the accurate modeling of Dirac-like-cone dispersion within non-Hermitian PC systems.

Appendix B the matrix $\bm{C}$ at the $\Gamma$ point

Owing to the $C_{4v}$ symmetry at the $\Gamma$ point, the coupling matrix $\bm{C}$ should be symmetric to some extent. As stated above, the coupling matrix $\bm{C}$ can be divided into four parts $\{\bm{C}_{1D},\bm{C}_{rad},\bm{C}_{2D},\bm{C}_{TE\text{-}TM}\}$. First part $\bm{C}_{1D}$ has a form as:

[TABLE]

where $\chi_{m,n}=\int_{PC}\kappa_{m,n}\{[-\partial^{2}/\partial z^{2}+\delta_{A}\partial/\partial z-\beta_{0}^{2}]\Theta_{0}(z)\}\Theta_{0}^{*}(z)dz$.

The second part $\bm{C}_{rad}$ corresponds to the couplings via the radiative waves, which can be evaluated with Green’s function. Further, we can cacluate $\bm{C}_{rad}$ as

[TABLE]

where $u_{mn,pq}=\iint_{PC}\kappa_{m,n}\kappa_{p,q}\{[-\partial^{2}/\partial z^{2}+\delta_{A}\partial/\partial z]\cdot G_{0,-}(z,z^{\prime})[(\partial^{2}/\partial z^{\prime 2}-\delta^{\prime}_{A}\partial/\partial z^{\prime})\Theta_{0}(z^{\prime})]\}\Theta_{0}^{*}(z)dz^{\prime}dz$ and Green’s function $G_{0,-}(z,z^{\prime})$ is the solution of $\mathcal{H}(m_{x}=0,n_{y}=0)G_{0,-}(z,z^{\prime})=-\delta(z,z^{\prime})$, with applying the boundary conditions Eq.(A) in Appendix A.

The third part $\bm{C}_{TE\text{-}TM}$ depicts the couplings between the basic waves $\{H_{x,10},H_{x,-10},H_{y,01},H_{y,0-1}\}$ and $\{H_{z,10},H_{z,-10},H_{z,01},H_{z,0-1}\}$. Similarly, we have:

[TABLE]

Green’s function $G_{st,+}(z,z^{\prime})$ is the solution of $\mathcal{H}(m_{x}=s+\Delta x,n_{y}=t+\Delta y)G_{st,+}(z,z^{\prime})=-\delta(z,z^{\prime})$ with applying the boundary condition Eq.(A). $G_{st,0}(z,z^{\prime})$ is the solution of the same equation but with continuous boundary condition.

As mentioned, the non-Hermiticity of matrix $\bm{C}$ only comes from the $\bm{C}_{rad}$ and $\bm{C}_{TE\text{-}TM}$ parts. At the symmetric $C_{4v}\text{-}\Gamma$ point, the coupling coefficients $\kappa_{m,n}=2f(1/\varepsilon_{1}-1/\varepsilon_{2})J_{1}(\zeta)/\zeta,\zeta=2\sqrt{f\pi(m^{2}+n^{2})}$ only depends on $(m^{2}+n^{2})$. Therefore, all the non-zero elements in $\bm{C}_{rad}$ and $\bm{C}_{TE\text{-}TM}$ are identical , and hence, the imaginary part of matrix $\bm{C}$ has a symmetric form like:

[TABLE]

The forth part $\bm{C}_{2D}$ corresponds to the in-plane 2D distributed couplings via high-order waves. As the close form of $\bm{C}_{2D}$ is too complex, here we prove that the $\bm{C}_{2D}$ matrix also possesses a symmetric form like $\bm{A}$ in Eq.(19) by analyzing the underlying coupling paths.

We consider the lowest set of the high-order waves that consists of four individual waves $\{(1,1),(1,-1),(-1,1),(1,1)\}$. Via these waves, there are 4 different coupling paths depicted by $\bm{C}^{2}_{2D}$. The superscript $1$ denotes that the matrix depicts the coupling via high-order waves of $m^{2}+n^{2}=2$.

Fig. 9 (a) illustrates one possible coupling path where the basic wave $(0,1)$ couple to the four high-order waves, and then back to itself. This coupling path is depicted by the matrix elements $\bm{C}^{2}_{2D}(3,3)$. Apparently, the four diagonal elements corresponding to this coupling path are identical.

Fig. 9 (b) presents another possible path between basic wave $(0,1)$ and $(0,-1)$. This coupling path is depicted by matrix elements $\bm{C}^{2}_{2D}(4,3)$, $\bm{C}^{2}_{2D}(3,4)$, $\bm{C}^{2}_{2D}(1,2)$ and $\bm{C}^{2}_{2D}(2,1)$.

Moreover, Fig. 9 (c) and (d) give the third and fourth possible paths between the basic waves $(0,1)$ and $(-1,0)$. These two coupling paths are depicted by matrix elements $\bm{C}^{2}_{2D}(2,3)$, $\bm{C}^{2}_{2D}(1,4)$, $\bm{C}^{2}_{2D}(3,2)$, $\bm{C}^{2}_{2D}(4,1)$, and by $\bm{C}^{2}_{2D}(1,3)$, $\bm{C}^{2}_{2D}(2,4)$, $\bm{C}^{2}_{2D}(3,1)$, $\bm{C}^{2}_{2D}(4,2)$, respectively. Apparently, the two paths given by Fig. 9 (c) and (d) have the same coupling strength but in an opposite sign.

Combined all the coupling paths, the matrix $\bm{C}^{2}_{2D}$ have the a symmetric form, as:

[TABLE]

Following the same approach, we can conclude that every $\bm{C}^{M}_{2D},M\geq 1$ has the same symmetric form with $\bm{C}^{1}_{2D}$, as a result, the entire $\bm{C}_{2D}$ owns the symmetric form as shown in Eq.(93). Finally, we prove that coupling matrix $\bm{C}$ has a form of Eq.(19).

Appendix C the matrix $\bm{K}$ at off-$\Gamma$ points

In this section, we prove that in the vicinity of the $\Gamma$ point, the eigenvalue $k_{0,mn},(m^{2}+n^{2}=1)$ of Eq.(7) linearly varies along $k_{x}$ and $k_{y}$. $k_{0,mn}$ is the solution to the dispersion equation Yariv (1973):

[TABLE]

where $t=2h$ is the thickness of PC layer. Coefficient $g,p$ and $q$ are the propagation constants of each layer:

[TABLE]

where $\beta=\sqrt{m_{x}^{2}+n_{y}^{2}}\beta_{0}$. $k_{0}$ is the solution of Eq.(94) when $\beta=\beta_{0}$.

From Eq.(C) we know $g,p$ and $q$ all depend on $k$ and $\beta$. We concern about an arbitrary direction $n_{y}=n+\Delta_{y},m_{x}=m+\Delta_{x}$, where $\Delta_{x,y}=k_{x,y}/\beta_{0}$. Expanding $g,p$ and $q$ in the vicinity of $\beta_{0}$, we can have

[TABLE]

where $g_{0},p_{0}$ and $q_{0}$ denote the values of $g,p$ and $q$ when $\beta=\beta_{0}$. Expand $\tan(gt)$ in the vicinity of $\beta_{0}$:

[TABLE]

Notice that both $\tan(g_{0}t)$ and $\cos^{2}(g_{0}t)$ don’t depend on $\Delta_{x}$ or $\Delta_{y}$ since $g_{0}$ only relies on $\beta_{0}$. Submitting Eqs.(C) and (97) into Eq.(94) and neglecting the high-order terms of $(m\Delta_{x}+n\Delta_{y})$, we obtain

[TABLE]

where

[TABLE]

$g_{0},p_{0},q_{0}$ and $\xi$ in Eq.(98) only rely on $\beta_{0}$ since their dependences to $(m\Delta_{x}+n\Delta_{y})$ have been separated into the linear terms. All of them are functions of unknown quantity $k$ and we can expand them near zero-order solution $k_{0}$ as

[TABLE]

where $\sigma=k-k_{0}$ and $g_{00},p_{00}$ and $q_{00}$ denote the values of $g,p$ and $q$ when $\beta=\beta_{0},k=k_{0}$. Since $\xi$ is a slow varying function with respect to $k$, we approximately take $\xi=\xi_{0}(k_{0})$ which can be evaluated through the structural parameters and $k_{0}$ directly.

Expand $\tan(gt)$ in the vicinity of $k_{0}$ similarly:

[TABLE]

Submitting Eqs.(C) and (100) to Eq.(98), we have

[TABLE]

where

[TABLE]

Apparently, Eq.(101) can be simplified as $\varsigma k_{0}\sigma+\xi\beta_{0}^{2}(m\Delta_{x}+n\Delta_{y})=0$, that is

[TABLE]

where $\eta=-\xi\beta_{0}^{2}/(\varsigma k_{0})$.

$\xi=\xi_{0}(k_{0})$ and $\varsigma$ all are determined by to $\beta_{0}$ and $k_{0}$, which can be evaluated by the zero-order trial solution and structural parameters. As $\bm{k}_{T}=(k_{x}\hat{e_{x}}+k_{y}\hat{e_{y}})=\beta_{0}(\Delta_{x}\hat{e_{x}}+\Delta_{y}\hat{e_{y}})$, Eq.(102) turns to be Eq.(41).

It is noteworthy that, the coefficient $\eta$ depends on a series of structural parameters. It becomes tedious to analytically evaluate $\eta$ according to Eq.(98) and (101). Alternatively, we can obtain $\eta$ by fitting the dispersion curve given by Eq.(94).