Exceptional points and their coalescence of PT-symmetric interface states in photonic crystals

TL;DR

This paper demonstrates that PT-symmetric interfaces in photonic crystals support unique surface states with real eigenfrequencies, exhibiting complex dispersion features like zig-zag trajectories and exceptional point coalescence, revealing rich non-Hermitian physics.

Contribution

It introduces a non-Hermitian Hamiltonian model for PT-symmetric photonic crystal interfaces, revealing novel dispersion behaviors and exceptional point phenomena in surface electromagnetic states.

Findings

Surface states with real eigenfrequencies exist at PT-symmetric interfaces.

Dispersion features include zig-zag trajectories and closed loops.

Exceptional points can coalesce into higher order EPs or anisotropic EPs.

Abstract

The existence of surface electromagnetic waves in the dielectric-metal interface is due to the sign change of real parts of permittivity across the interface. In this work, we demonstrate that the interface constructed by two semi-infinite photonic crystals with different signs of the imaginary parts of permittivity also supports surface electromagnetic eigenmodes with real eigenfrequencies, protected by $ PT $ symmetry of such loss-gain interface. Using a multiple scattering method and full wave numerical methods, we show that the dispersion of such interface states exhibit unusual features such as zig-zag trajectories or closed-loops. To quantify the dispersion, we establish a non-Hermitian Hamiltonian model that can account for the zig-zag and closed-loop behaviour for arbitrary Bloch momentums. The properties of the interface states near the Brillouin zone center can also be explained within the framework of effective medium theory. It is shown that turning points of the dispersion are exceptional points (EPs), which are characteristic features of non-Hermitian systems. When the permittivity of photonic crystal changes, these EPs can coalesce into higher order EPs or anisotropic EPs. These interface modes hence exhibit and exemplify many complex phenomena related to exceptional point physics.

Equations

$$M(\omega ,k_y){\Psi }^{\mathrm{loc}}={\Psi }^{\mathrm{inc}},$$

$$M(\omega ,k_y){\Psi }^{\mathrm{loc}}=0.$$

$$\det \left[M\left(\omega ,k_y\right)\right]=0.$$

$$\mathrm{Min}\left\{\left|\mathrm{eig}\left[M\left(\omega ,k_y\right)\right]\right|\right\}=0,$$

$$Z=\frac{E_z}{H_y}=-\frac{\omega \mu }{k}=-\frac{k}{\omega \epsilon }$$

$$Z_B=\frac{{\int }_I{E}_z^{\mathrm{PC}}dy}{{\int }_I{H}_y^{\mathrm{PC}}dy},$$

$${\epsilon }_{\mathrm{e}}=-\frac{k}{\omega {\epsilon }_0{Z}_B},{\mu }_{\mathrm{e}}=-\frac{k}{{\mu }_0\omega }{Z}_B.$$

$$\varepsilon(x)=\left\{\begin{array}[]{lll}{\varepsilon_{\mathrm{e}}^{*},}&{x<0}\\ {\varepsilon_{\mathrm{e}},}&{x>0}\end{array},\quad\mu(x)=\left\{\begin{array}[]{ll}{\mu_{\mathrm{e}}^{*},}&{x<0}\\ {\mu_{\mathrm{e}},}&{x>0}\end{array}\right.\right..$$

$$E_{z}(x,y)=C\exp(i\beta y)\left\{\begin{array}[]{ll}{\exp\left(ik_{x}^{*}x\right),}&{x<0}\\ {\exp\left(ik_{x}x\right),}&{x>0}\end{array}\right.$$

$$H_y(x,y)=\frac{-k_x(x)}{\omega \mu (x)}{E}_z(x,y).$$

$$\frac{k_x}{{\mu }_{\mathrm{e}}}=\frac{k_x^{\ast }}{{\mu }_{\mathrm{e}}^{\ast }}.$$

$$\frac{{\mu }_{\mathrm{e}}}{{\mu }_{\mathrm{e}}^{\ast }}=\frac{k_x}{k_x^{\ast }}=\frac{i\sqrt{{\beta }^2-{\epsilon }_{\mathrm{e}}{\mu }_{\mathrm{e}}(\omega /c{)}^2}}{-i\sqrt{{\beta }^2-{\epsilon }_{\mathrm{e}}^{\ast }{\mu }_{\mathrm{e}}^{\ast }(\omega /c{)}^2}},$$

$$\beta =\frac{\omega }{c}\sqrt{\frac{{\mu }_{\mathrm{e}}{\mu }_{\mathrm{e}}^{\ast }\left(-{\epsilon }_{\mathrm{e}}{\mu }_{\mathrm{e}}^{\ast }+{\epsilon }_{\mathrm{e}}^{\ast }{\mu }_{\mathrm{e}}\right)}{{\mu }_{\mathrm{e}}^2-{\left({\mu }_{\mathrm{e}}^{\ast }\right)}^2}}=\frac{\omega \left|{\mu }_{\mathrm{e}}\right|}{c}\sqrt{\frac{1}{2}\left(\frac{{\epsilon }^{\prime }}{{\mu }^{\prime }}-\frac{{\epsilon }^{\prime \prime }}{{\mu }^{\prime \prime }}\right)},$$

$$k_y(\omega )=\frac{\omega \left|{\mu }_{\mathrm{e}}\right|}{c}\sqrt{\frac{1}{2}\left(\frac{p}{q}\frac{\omega -{\omega }_{\mathrm{e}}}{\omega -{\omega }_m}-\rho \right)}$$

$${\omega }_{\beta }=\frac{{\omega }_{\mathrm{e}}\left[1-\eta \left({\omega }_m/{\omega }_{\mathrm{e}}\right)\right]}{1-\eta },$$

$$\left[{\nabla }^2+{\left(\frac{\omega }{c}\right)}^2\epsilon (\mathbf{r})\right]E_z(\mathbf{r})=0,$$

$$\varepsilon^{(0)}(\mathbf{r})=\left\{\begin{array}[]{cccc}{\varepsilon_{r}^{(0)}-i\gamma^{(0)}}&{\left|\mathbf{r}-\mathbf{r}_{ij}\right|<r_{c}}&{\text{ and }}&{x<0}\\ {\varepsilon_{r}^{(0)}+i\gamma^{(0)}}&{\left|\mathbf{r}-\mathbf{r}_{ij}\right|<r_{c}}&{\text{ and }}&{x>0}\\ {1}&{\left|\mathbf{r}-\mathbf{r}_{ij}\right|>r_{c}}\end{array}\right.$$

$$\int {\mathrm{d}}^{2}r{v}_{\mathbf{k}{m}^{\prime }}^{\left(0\right)\ast }\left(\mathbf{r}\right){\epsilon }^{\left(0\right)}\left(\mathbf{r}\right)u_{\mathbf{k}m}^{\left(0\right)}\left(\mathbf{r}\right)={\delta }_{m^{\prime }m},$$

$$\displaystyle\sum\limits_{m}{{\alpha_{n,{\bf{k}}m}}}\left\{\begin{array}[]{l}{\nabla^{2}}\left({u_{{\bf{k}}m}^{\left(0\right)}\left({\bf{r}}\right){e^{i{\bf{k}}r}}}\right)\\ +\varepsilon\left({\bf{r}}\right){\left({\frac{{\omega_{{\bf{k}}n}}}{{\rm{c}}}}\right)^{2}}u_{{\bf{k}}m}^{\left(0\right)}\left({\bf{r}}\right){e^{i{\bf{k}}r}}\end{array}\right\}=0.$$

$${\nabla }^2\left(u_{\mathbf{k}m}^{\left(0\right)}\left(\mathbf{r}\right)e^{i\mathbf{k}r}\right)+{\epsilon }^{\left(0\right)}\left(\mathbf{r}\right){\left(\frac{{\omega }_{\mathbf{k}m}^{\left(0\right)}}{\mathrm{c}}\right)}^2{u}_{\mathbf{k}m}^{\left(0\right)}\left(\mathbf{r}\right)e^{i\mathbf{k}r}=0.$$

$$\sum\limits_{m}{{\alpha_{n,{\bf{k}}m}}}\left[\begin{array}[]{l}\varepsilon\left({\bf{r}}\right){\left({\frac{{\omega_{{\bf{k}}n}}}{{\rm{c}}}}\right)^{2}}\\ -\varepsilon^{\left(0\right)}\left({\bf{r}}\right){\left({\frac{{\omega_{{\bf{k}}m}^{\left(0\right)}}}{{\rm{c}}}}\right)^{2}}\end{array}\right]u_{{\bf{k}}m}^{\left(0\right)}\left({\bf{r}}\right){e^{i{\bf{k}}r}}=0,$$

$$\displaystyle\begin{array}[]{l}\sum\limits_{m=1}^{+\infty}{{\alpha_{n,{\bf{k}}m}}}{\left({\bar{\varepsilon}^{\left(0\right)}\left({\bf{k}}\right)}\right)_{m^{\prime}m}}{(\omega_{{\bf{k}}m}^{\left(0\right)}/c)^{2}}\\ ={({\omega_{{\bf{k}}n}}/c)^{2}}\sum\limits_{m=1}^{+\infty}{{\alpha_{n,{\bf{k}}m}}}{\left({{{\bar{\varepsilon}}}\left({\bf{k}}\right)}\right)_{m^{\prime}m}}\end{array}$$

$${\left({\bar{\epsilon }}^{\left(0\right)}\left(\mathbf{k}\right)\right)}_{m^{\prime }m}=\int {\mathrm{d}}^{2}r{v}_{\mathbf{k}{m}^{\prime }}^{\left(0\right)\ast }\left(\mathbf{r}\right){\epsilon }^{\left(0\right)}\left(\mathbf{r}\right)u_{\mathbf{k}m}^{\left(0\right)}\left(\mathbf{r}\right),$$

$${\left(\bar{\epsilon }\left(\mathbf{k}\right)\right)}_{m^{\prime }m}=\int {\mathrm{d}}^{2}r{v}_{\mathbf{k}{m}^{\prime }}^{\left(0\right)\ast }\left(\mathbf{r}\right)\epsilon \left(\mathbf{r}\right)u_{\mathbf{k}m}^{\left(0\right)}\left(\mathbf{r}\right).$$

$$H_2\cdot p_{\mathbf{k}n}={\left({\omega }_{\mathbf{k}n}/c\right)}^2{H}_1\cdot p_{\mathbf{k}n}$$

$$H\cdot p_n=W_n{p}_n,$$

$${\text{E}}^{\mathrm{inc}}=E_0\exp \left[i\left(k_{p}y+k_{x}x\right)\right]\hat{z},$$

$$H_x^{\mathrm{inc}}$$

$$H_y^{\mathrm{inc}}$$

$${\mathbf{K}}_m^{\pm }=\left(\pm q_m,k_p-g_m\right),$$

$$P\hat{z}$$

Full text

I Introduction

Interface states commonly exist in quantum systems and classical wave systems. A well-known example is the surface plasmon polaritons Barnes et al. (2003); Popov et al. (2008), which are surface waves traveling along a dielectric-metal interface due to a change in sign of the real part of permittivity across the interface. Photonic crystal (PC) systems also have surface modes, and in some cases, the existence of the boundary modes can be explained using topological concepts Wang et al. (2009); Feng et al. (2011); Fang et al. (2012); Rechtsman et al. (2013); Huang et al. (2014); Xiao et al. (2016); Lu et al. (2016). These prior topological-based research on interface states in photonic crystal systems were mainly focused on Hermitian systems. However, recent works show that boundary modes can also be found in non-Hermitian systems Longhi (2014); Savoia et al. (2015); Weimann et al. (2016); Xu et al. (2017); Shen et al. (2018); Yao and Wang (2018); Martinez Alvarez et al. (2018); Ni et al. (2018). For example, surface states can be localized at the gain-loss interface in $\mathcal{PT}$-symmetric systemsLonghi (2014); Savoia et al. (2015).

In this work, we study the formation of interface states in a non-Hermitian PC with a $PT$-symmetric interface, in which one side of the PC has gain and the other side has loss. We find that such a system carries interface states with real eigenvalues. The dispersions of these interface states are rather unusual, as they form ziz-zag trajectories or closed-loops. Exceptional points Bender and Boettcher (1998); Heiss (2012); Hodaei et al. (2014); Zhen et al. (2015); Hodaei et al. (2017); El-Ganainy et al. (2019); Miri and Alù (2019), which are characteristic features in non-Hermitian systems, appear at the turning points of the dispersions of the interface states. EPs are branch point singularities in parameter spaces, at which eigenvalues and eigenvectors coalesce simultaneously. At EPs, the matrix Hamiltonian is defective, and the coalescing eigenvectors are not linearly independentRotter (2009), which is different from degenerate points in Hermitian systems. As the system parameters changes, we found that the EPs can coalesce into higher order EPsGraefe et al. (2008); Demange and Graefe (2012); Ding et al. (2015, 2016); Zhang and Chan (2019) or anisotropic EPs Ding et al. (2018); Zhang and Chan (2018); Xiao et al. (2019). Besides, we find that in the limit of large gain/loss, one band of interface states with real eigenvalues always persists.

We study the interface states using two different computation approaches (numerical package COSMOL and a multiple scattering method) and two different boundary conditions (periodic and open). These computation details are described in Sec. II and Sec. III. In Sec. IV, we attempted to give a simple explanation to the rather exotic dispersions using effective medium theory (EMT), which works well near the Brillouin zone center. In Sec. V, we formulate a non-Hermitian Hamiltonian model for the interface states that works for a general value of the Bloch momentum and the model shows clearly that there are EPs in the band of interface states. In Sec. VI, we show that EPs can coalesce into higher order EPs and anisotropic EPs as system parameters change. In the last section, we give a conclusion.

II PT-symmetric interface states in photonic crystals

We consider a PC comprising of two different two-dimensional (2D) square lattice PCs, as illustrated in Fig. 1(a). In this work, we focus on the transverse magnetic (TM) polarization ( $E_{z}$ polarization). Within the primitive unit cell with a lattice constant $a$, the rod has a radius of $r_{c}$ and a relative permittivity $\varepsilon_{c}=\varepsilon_{r}\pm i\gamma$ embedded in air. The relative permeability $\mu$ is 1.0 everywhere. The positive (negative) sign of $\gamma$ indicates the rod is a lossy (gain) medium. There is a $PT$-symmetric interface at $x=0$ separating cylinders with gain on one side and lossy cylinders on the other side. The cylinders (blue) of the left semi-infinite PC ($x<0$) are active with gain ($\varepsilon_{c}=\varepsilon_{r}-i\gamma$) while the cylinders (orange) at the right semi-infinite PC ($x>0$) are lossy ($\varepsilon_{c}=\varepsilon_{r}+i\gamma$).

At the Hermitian limit $\gamma=0$, the $\mathcal{PT}$-symmetric interface disappears, and this system becomes a simple Hermitian 2D PC with the same material parameters everywhere. The bulk band structure and projected band diagrams along the $k_{y}$ direction of the Hermitian PC with parameters $r_{c}=0.2a,\varepsilon_{c}=12$ are calculated using COMSOL as shown in Figs. 1(b) and 1(c), respectively. It has been shown theoretically and demonstrated experimentally that in 2D Hermitian PCs possessing bands with a Dirac-like cone dispersion at $k=0$, even a small perturbation of the relative permittivity $\varepsilon_{r}$ is sufficient to create interface states Xiao et al. (2014); Huang et al. (2014); Xiao et al. (2016); Yang et al. (2016). The underlying physics for the existence of an interface state is related to the geometric phase of the bulk bands. In this work, we want to study the emergence of interface states upon the introduction of an imaginary part to the relative permittivity. More specifically, we create an interface by introducing $PT$-symmetric non-Hermiticity into this Hermitian 2D PCs as depicted in Fig. 1(a).

In numerical calculations, the number of column layers of the semi-infinite PC is truncated to a finite number $N$. One method of terminating the far ends of the semi-infinite PC is to apply open boundary conditions, such as perfectly matched layers (PML) or scattering boundary along $x$-direction. Thus, the $PT$-symmetric PC is periodic along $y$-direction but is finite along $x$-direction. Another method to study the interface states is to apply periodic boundary condition also along $x$-direction. Figure 2(a) gives a schematic picture of such a lattice, with a unit cell containing $N=3$ cylinders with gain and $N=3$ cylinders with loss. Both boundary conditions give the same result in the limit of large $N$, showing that in this particular case the results are not dependent on boundary conditions even though the system is non-Hermitian Lee (2016). The periodic condition gives us a more heuristic understanding on the formation of the interface bands from the point view of $PT$-symmetry, which has been extensively studied in prior research.

We first calculate the band structure of the $PT$-symmetric PC with periodic boundary conditions. As shown in Fig. 2(a), the unit cell (marked by black dashed rectangle) of the lattice contains $N$ active cylinders (orange) together with $N$ lossy cylinders (blue). The gray shadow region of Figs. 2(b-d) marks the projected bands of Hermitian PC with parameters $r_{c}=0.2a$, $\varepsilon_{c}=12$, which is the same as the projected band shown in Fig. 1(c). We then calculate the projected bands of $PT$-symmetric PC with parameters $\varepsilon_{c}=12\pm i$ using COMSOL. In Fig. 2(b-d), we plot the real parts of eigenfrequencies of PC, and the number $N$ of gain/loss cylinders in the unit cell is labeled in the figure. The imaginary parts of the eigenfrequency of eigenstates marked by blue circles are zero, and those by red are nonzero. For $PT$-symmetric systems, the eigenstates with purely real eigenvalues are said to be in the exact $PT$ symmetry phase, and those with complex-conjugate-pairs eigenvalues are said to be in the broken $PT$ symmetry phase Bender (2007). From Fig. 2(b), we see that when we introduce the $PT$-symmetric non-Hermiticity into Hermitian supercells of the 2D PC, some eigenstates (marked by red circles) acquire imaginary parts and go into the broken $PT$-symmetry phase. When the supercell is small (e.g. $N=3$), most of the eigenstates (marked by blue circles) are purely real, and these eigenstates are still in the exact $PT$ symmetry phase. When we increase the number of cylinders in the unit cell $N$ to 8, an increasing number of the eigenstates go into the broken $PT$ symmetry phase as depicted in Fig. 2(c). When the number of cylinders in the unit cell $N=15$ as shown in Fig. 2(d), nearly all of the eigenstates in the whole Brillouin zone are in the broken $PT$ symmetry phase, except for one band of states with a zigzag dispersion which remains in the exact $PT$ symmetry phase. In Fig. 2(e), we plot one of these states in $PT$ symmetry exact phase (marked by dark star in Fig. 2 (d)), and find that the electric field distribution $|\textbf{E}|$ is localized at the interface $x=\{-15a,0,15a\}$ separating the gain/loss cylinders, i.e., the $PT$-symmetric interface. The interface modes lie inside the continuum of the real-valued bulk modes when $\gamma=0$ (gray shadow region), but when $\gamma\neq 0$, the continuum becomes complex-valued. Therefore, in a certain sense these interface modes are bound states in the continuum with complex-valued energy spectra. This is different from the conventional “bound states in continuum”, where the continuum is real-valuedLonghi (2014); Longhi and Della Valle (2014).

We emphasize that in going from panels (b) to (d) of Fig. 2, the material parameters and in particular the imaginary parts of the permittivity remain unchanged. The only change is the number of cylinders in the unit cell. There are both localized interface states and extended bulk states in this $\mathcal{PT}$-symmetric PC. For a fixed $\gamma=1$, as $N$ increases, the bulk states undergo a phase transition from the exact $PT$ regime to the broken $PT$ regime. For bulk states at different Bloch $k_{y}$, the threshold value of $N$ for the phase transition point is different. When the bulk states are in the broken $PT$ regime, the field distribution is either concentrated in the lossy medium or in the gain medium. If the bulk state eigenfields are mainly concentrated in the left gain domain, the eigenfrequencies have negative imaginary parts. If the bulk states are mainly localized in the right lossy domain, the eigenfrequencies have positive imaginary parts. However, the fields of interface states are always localized at the loss-gain interface, and the field distributions are symmetric about the interface even for large value of $N$. Accordingly, the interface states persist in the exact $PT$ regime and have purely real eigenfrequencies. By studying how the eigenstates change as the number of cylinders $N$ in the unit cell increases, we distinguish the localized interface states and extended bulk states based on whether their eigenvalues are real numbers.

III $\mathcal{PT}$-symmetric interface states using multiple scattering method with open boundary condition

When the column number $N$ is large enough, the results calculated with periodic condition should be the same with open boundary conditions. In this section, we calculate the band structure using a multiple scattering (MS) method with open boundary conditions, and we are expected to observe the interface states plotted in Fig. 2(d) calculated using periodic boundary conditions. We use a Green’s function method to build a scattering problem of cylinders arranged in a square lattice grid, the details of which are given in Appendix A. In the TM polarization (the electric field is along the axes of cylinders), the dominant excitations of the cylinders at low frequencies are the out-of-plane electric monopole and two in-plane magnetic dipoles. The scattering problem for plane waves incident on the PC can be written as below,

[TABLE]

where $\Psi^{\mathrm{inc}}$ represents external incident waves and $\Psi^{\mathrm{loc}}$ represents total local fields. The band structure of a PC can be obtained by solving an eigenvalue problem in the absence of incident field,

[TABLE]

Frequencies $\omega$ for which Eq. (2) has non-trivial solutions for a given Bloch $k_{y}$ are zeros of matrix $M(\omega,k_{y})$, which can be found by locating the zeros of its determinant,

[TABLE]

As the calculation of the determinant becomes numerically unstable if the matrix is large, we adopt the numerically stable approach to solve the equation

[TABLE]

where ${\rm{Min}}\left\{{\left|{{\mathop{\rm eig}\nolimits}\left[{{M}\left({\omega,{k_{y}}}\right)}\right]}\right|}\right\}$ is the minimal-modulus eigenvalue of the matrix $M\left(\omega,k_{y}\right)$. Therefore, equation (4) defines the dispersion curves $\omega(k_{y})$ for photons propagating through the periodic structure. The MS method is very useful as it can handle non-Hermitian systems and is compatible with various types of boundary conditions.

Using the MS method, we calculate the interface states of a $\mathcal{PT}$-symmetric PC with $\varepsilon_{c}=12\pm i$. The column number of the unit cell of the PC is truncated to $N=15$, and scattering boundary conditions applied in the $x$-direction . In Fig. 3, we plot the values of ${\rm{Min}}\left\{{\left|{{\mathop{\rm eig}\nolimits}\left[{{M}\left({\omega,{k_{y}}}\right)}\right]}\right|}\right\}$ in log-10 scale, in the parameter plane $\left(\omega,k_{y}\right)$, where $\omega$ and $k_{y}$ are both real numbers. The value ${\rm{Min}}\left\{{\left|{{\mathop{\rm eig}\nolimits}\left[{{M}\left({\omega,{k_{y}}}\right)}\right]}\right|}\right\}$ of the states with complex eigenfrequencies cannot be zero at real-valued $\left(\omega,k_{y}\right)$ plane. Therefore, the dark orange lines, marking locations of ${\rm{Min}}\left\{{\left|{{\mathop{\rm eig}\nolimits}\left[{{M}\left({\omega,{k_{y}}}\right)}\right]}\right|}\right\}\to 0$ represent interface eigenstates with real eigenfrequencies. In the vicinity, the areas marked by lighter colors (${\rm{Min}}\left\{{\left|{{\mathop{\rm eig}\nolimits}\left[{{M}\left({\omega,{k_{y}}}\right)}\right]}\right|}\right\}>0$) indicate the bulk states with complex eigenfrequencies. Using the MS method, we can locate the interface states with purely real eigenfrequencies.

To verify the interface states obtained by MS method, we also calculated the eigenstates by COMSOL, with PML boundary conditions applied to the $x$-direction, whereas a periodic boundary condition is applied to the $y$-direction for each wave number $k_{y}$. Using COMSOL, we pick the eigenstates with real eigenvalues and plot them in Fig. 3 by cyan triangles. We also plot the interface states calculated by COMSOL with periodic conditions applied to the $x$-directions by blue circles in Fig. 2(d) in the same picture. We can see from Fig. 3 that the results calculated by COMSOL using PML and periodic boundary conditions agree well with the results calculated by the MS method.

IV Effective Medium Theory

We now consider the interface states from the viewpoint of effective medium theory (EMT). We will see whether an effective medium description in the long wavelength limit can provide a simple heuristic picture to understand the formation of interface states near the zone center. The effective parameters of non-Hermitian PCs can be obtained conveniently using a boundary field averaging method Andryieuski et al. (2012). To obtain the effective parameters of the PC, we need to treat the PC as a homogeneous medium. Assuming that a plane wave propagates along the $x$-direction with wave vector $k\hat{x}$ and polarization $E_{z}\hat{z}$ in a homogeneous medium, we define a wave impedance,

[TABLE]

For the inhomogeneous PC system with a micro-structure, we define a corresponding average field ratio as,

[TABLE]

where $E_{z}^{\mathrm{PC}}$ and $H_{y}^{\mathrm{PC}}$ are the eigenfield at the incident boundary $I$ of the unit cell. When using EMT, the frequency $\omega$ in Eq. (5) should take real number, because the effective parameters are functions of real frequency. To obtain eigenfields $E_{z}^{\mathrm{PC}}$ and $H_{y}^{\mathrm{PC}}$, we use COMSOL to solve a complex-valued $\textbf{k}(\omega)$ vs. real-valued $\omega$ dispersion Davanco et al. (2007); Fietz et al. (2011). According to Eq. (5), the effective permittivity and permeability are

[TABLE]

We first calculate the effective parameters $\varepsilon_{\rm e}$ and $\mu_{\rm e}$ for the right lossy semi-infinite PC with $\varepsilon_{c}=\varepsilon_{r}+i\gamma$. Therefore, the effective parameters for the left active semi-infinite PC with $\varepsilon_{c}=\varepsilon_{r}-i\gamma$ are $\varepsilon_{\rm e}^{*}$ and $\mu_{\rm e}^{*}$. From the point of view of EMT, the PT-symmetric PC illustrated in Fig. 1(a) can be treated as a $\mathcal{PT}$-symmetric homogenous slab with permittivity and permeability like below,

[TABLE]

Enforcing the continuity condition at the interface $x=0$ for the electric field, a modal solution bound at the gain-loss interface can be written as Savoia et al. (2015)

[TABLE]

where $C$ denotes a normalization constant, $\beta$ is the propagation constant. We suppose the interface states are propagating along the $y$-direction without attenuation, which means that $\beta$ is a purely real number, and therefore, we obtain that $k_{x}^{2}=\varepsilon_{\mathrm{e}}\mu_{\mathrm{e}}(\omega/c)^{2}-\beta^{2}$ for $x>0$ and $k_{x}^{*2}=\varepsilon_{\mathrm{e}}^{*}\mu_{\mathrm{e}}^{*}(\omega/c)^{2}-\beta^{2}$ for $x<0$. $c$ is the speed of light in vacuum. From Maxwell’s equation, we then calculate the tangential magnetic field,

[TABLE]

By enforcing its continuity at the interface, we obtainLawrence et al. (2010)

[TABLE]

To obtain states bounded at $\mathcal{PT}$-symmetric interface, the complex number $k_{x}$ should satisfy $\operatorname{Im}\left(k_{x}\right)>0$ to ensure that the electric field exponentially decays for $x>0$. Accordingly, we have $\operatorname{Im}\left(k_{x}^{*}\right)<0$ for $x<0$, which also ensures the electric field exponentially decays away from the interface. We then solve the dispersion relationship like below

[TABLE]

and obtain

[TABLE]

where we have set $\mu_{\rm e}=\mu^{\prime}+i\mu^{\prime\prime},\varepsilon_{\rm e}=\varepsilon^{\prime}+i\varepsilon^{\prime\prime}$.

Using Eq. (7), we calculate the effective parameters $\varepsilon_{\rm e}$ and $\mu_{\rm e}$ of the lossy PC with relative permittivity $\varepsilon_{c}=12+1.2i$ and $\varepsilon_{c}=12+1.6i$. The results are plotted in Figs. 4(a) and 4(c), respectively. The horizontal axes represent the real parts or imaginary parts of the effective parameters. Using the dispersion relation described by Eq. (13), we can calculate the band dispersion of the interface states for the $\mathcal{PT}$-symmetric homogenous slab with effective parameters described by Eq. (8). The dispersions are plotted in Figs. 4(b) and 4(d) by cyan lines. The propagation constant $\beta$ in the effective homogeneous medium is the Bloch $k_{y}$ in the PC. For comparison, we further calculate the band structure of interface states near the Brillouin zone center using COMSOL and plot the results in Figs. 4(b) and 4(d) by open blue circles, which agree reasonably well with the cyan lines.

Since EMT is accurate in the long wavelength limit, where $k_{y}$ is a small number, we can give an analytical explanation to the band diagrams of the interface states near the Brillouin zone center. From Figs. 4(a) and 4(c), we see that $\mu^{\prime}(\omega)$ goes from negative to positive as frequency increases, and it becomes zero at a particular frequency, which we call $\omega_{m}$. In Fig. 4, we denote the frequency $\omega=\omega_{m}$ satisfying $\mu^{\prime}(\omega_{m})=0$ by blue dashed lines. It is shown by Eq. (13) that $k_{y}\to\infty$ at frequency $\omega=\omega_{m}$ due to $\mu^{\prime}(\omega_{m})=0$, which is depicted by the cyan lines in Figs. 4(b) and 4(d). From Figs. 4(a) and 4(c), we can see that the imaginary parts (denoted by asterisks) of the effective parameters are nearly constant functions of frequency while the real parts are linear functions of frequency. Therefore, near the Brillouin zone center, we can write a reasonably good approximation that $\varepsilon^{\prime}=p\left(\omega-\omega_{e}\right),\quad\mu^{\prime}=q\left(\omega-\omega_{m}\right)$ and $\varepsilon^{\prime\prime}/\mu^{\prime\prime}=\rho$. $p$, $q$, and $\rho$ are all positive numbers, and $\omega_{e}$ is the frequency satisfying $\varepsilon^{\prime}(\omega_{e})=0$. We label the frequency $\omega=\omega_{e}$ by red dashed lines in Fig. 4. Then Eq. (13) becomes

[TABLE]

By setting $k_{y}(\omega)=0$, we can solve the frequency $\omega_{\beta}$ satisfying $k_{y}(\omega_{\beta})=0$ as

[TABLE]

where $\eta=\rho q/p\ll 1$ in the considered frequency range. In Fig. 4, we label the frequency $\omega_{\beta}$ satisfying $k_{y}\left(\omega_{\beta}\right)=0$ by orange dots

From Eq. 15, we see that if $\omega_{m}>\omega_{e}$, then we have $\omega_{e}>\omega_{\beta}$ , which is the case illustrated in Figs. 4(a) and 4(b). On the other hand, if $\omega_{m}<\omega_{e}$, we have $\omega_{e}<\omega_{\beta}$, which is the case illustrated in Figs. 4(c) and 4(d). The inversion of ordering ( $\omega_{m}>\omega_{\beta}$ in Fig. 4(b) while $\omega_{m}<\omega_{\beta}$ in Fig. 4(d)) leads to a drastic change in the dispersion. The frequency range between $\omega_{m}$ and $\omega_{\beta}$ is a band gap induced by quasi-longitudinal resonance. Hence, $\omega_{e}$, $\omega_{m}$ and $\omega_{\beta}$, which are determined by $\varepsilon_{r}$, govern the band structure’s pattern of the interface states.

V Non-Hermitian Hamiltonian model and exceptional points of the interface states

In the above section, we study the dispersion of interface states as a function of $k_{y}$. In this section, we will formulate a non-Hermitian Hamiltonian model of the interface states as functions of $\varepsilon_{c}=\varepsilon_{r}\pm\gamma$ for a fixed $k_{y}$ Ding et al. (2015); Cui et al. (2019). For a 2D PC of which the cylinders are uniform in the $z$-direction ($\varepsilon(\textbf{r})$ is independent of $z$ coordinate), we consider the TM polarization with electric fields only having $z$ component $E_{z}(\textbf{r})$. The wave vector k is parallel to the 2D $x-y$ plane, and the electromagnetic fields $E_{z}(\textbf{r})$ , $H_{x}(\textbf{r})$ and $H_{y}(\textbf{r})$ are also independent of the $z$ coordinate. The eigenfrequency and eigenstate of this 2D PC with TM polarization can be obtained by solving the following equation Sakoda (2004)

[TABLE]

where $\nabla^{2}=\partial^{2}x+\partial^{2}y$, and r denotes the 2D position vector ($x,y$).

To calculate the eigenfrequencies of the interface states of non-Hermitian PCs with material parameters $\varepsilon(\textbf{r})$, we construct a model Hamiltonian using the Bloch states of the interface bands of a PC with material parameters $\varepsilon^{(0)}(\textbf{r})$ as the bases. Here, $\varepsilon^{(0)}(\textbf{r})$ is the relative permittivity of the original PC, and $\varepsilon(\textbf{r})$ is the modified relative permittivity of the PC under consideration. The relative permittivity of the original $\mathcal{PT}$-symmetric PC in Fig. 1(a) can be described as

[TABLE]

where $\mathbf{r}_{ij}$ is the position vector of the rod’s center. $\left|\mathbf{r}-\mathbf{r}_{ij}\right|<r_{c}$ denotes the rod domain, and $\left|\mathbf{r}-\mathbf{r}_{ij}\right|>r_{c}$ denotes the air domain. Note that $\varepsilon^{(0)}(\mathbf{r})$ has imaginary parts and the system is non-Hermitian.

The Bloch states for the interface states of PC with $\varepsilon^{(0)}(\mathbf{r})$ can be expressed as $E_{z,\mathrm{\textbf{k}n}}^{(0)}(\mathbf{r})=u_{\mathbf{k}n}^{(0)}(\mathbf{r})e^{i\mathbf{k}\cdot\mathbf{r}}$, where $n$ denotes the band index and k is the wave vector in the first Brillouin zone. The periodic function $u_{\mathbf{k}n}^{(0)}(\mathbf{r})$ and the corresponding eigenfrequency $\omega_{\mathrm{\textbf{k}n}}^{(0)}$ can be obtained from COMSOL. The Bloch states of non-Hermitian PC should satisfy a biorthonormal relationship as below

[TABLE]

where $v_{{\bf{k}}m^{\prime}}^{\left(0\right)}\left({\bf{r}}\right)$ is the normalized left eigenfield and can be obtained by the normalized right eigenfield at $-\textbf{k}$, i.e. $v_{{\bf{k}}m}^{\left(0\right)*}\left({\bf{r}}\right)=u_{-{\bf{k}}m}^{\left(0\right)}\left(r\right)$ [See Appendix B].

To construct the model Hamiltonian for a PC with a new permittivity $\varepsilon(\mathbf{r})$, which has the same functional form as Eq. (17). we express Bloch wave functions of the new PC as $E_{z}(\mathbf{r})=E_{z,\mathbf{k}n}(\mathbf{r})=u_{\mathbf{k}n}(\mathbf{r})e^{i\mathbf{k}\cdot\mathbf{r}}$ with $u_{\mathbf{k}n}(\mathbf{r})=\sum_{m=1}^{+\infty}\alpha_{n,\mathbf{k}m}u_{\mathbf{k}m}^{(0)}(\mathbf{r})$. Substituting this expansion into Eq. (16), we arrive at

[TABLE]

For the original PC with $\varepsilon^{\left(0\right)}\left({\bf{r}}\right)$, we know that

[TABLE]

Then Eq. (21) becomes

[TABLE]

where $\omega_{\mathrm{\textbf{k}n}}$ is eigenfrequency for PCs with $\varepsilon(\mathbf{r})$ . Multiplying Eq. (23) by $v_{{\bf{k}}m^{\prime}}^{\left(0\right)*}\left({\bf{r}}\right)$ and integrating within a unit cell, we obtain

[TABLE]

in which

[TABLE]

Therefore, Eq. (26) can be rewritten as a generalized eigenvaule problem for a specific k as

[TABLE]

where $p_{\mathrm{\textbf{k}n}}=\left(\cdots\alpha_{n,\mathbf{k}n},\cdots\right)^{T}$ is the eigenstate. The matrices in Eq. (28) are $\left(H_{1}\right)_{m^{\prime}m}=(\overline{\varepsilon}(\mathbf{k}))_{m^{\prime}m}$ and $\left(H_{2}\right)_{m^{\prime}m}=\delta_{m^{\prime}m}\left(\omega_{\mathrm{\textbf{k}m}}^{(0)}/c\right)^{2}$. For convenience, we omit the subscript k for simplicity and rewrite Eq.(28) as

[TABLE]

where $H=H_{1}^{-1}\cdot H_{2}$, and $W_{n}=\left(\omega_{n}/c\right)^{2}$. From Eqs. (27) and (28), we see that the model Hamiltonian is a function of $\varepsilon(\textbf{r})$ for a fixed $\textbf{k}=k_{y}\hat{y}$. Therefore, using the eigenfunctions of the interface states of a system with permittivity $\varepsilon^{(0)}(\mathbf{r})$ at a specific $k_{y}$ as bases, we build a Hamiltonian model to calculate a PC with new permittivity $\varepsilon(\mathbf{r})$. This approach works well for any value of $k_{y}$.

Using the model Hamiltonian, we can analyze the dispersion of interface states as functions of the material parameters $\gamma$ and $\varepsilon_{r}$. We first consider the transition of the interface band dispersion depicted in Figs. 4(b) and 4(d), in which $\varepsilon_{r}$ changes from 12 to 12.6. We calculate the interface states with parameters $\varepsilon_{r}^{(0)}=12.3$ and $\quad\gamma^{(0)}=1.2$ using COMSOL, and focus on the three interface states with real eigenvalues at a specific $k_{y}a/2\pi=0.05$ (orange dashed lines in Fig. 5(a)). We use the eigenfunctions $u_{\mathbf{k}n}^{(0)}(\mathbf{r})$ of these three interface states as bases to build the Hamiltonian model shown in Eq. (29). The Hamiltonian $H$ is a $3\times 3$ matrix function of $\varepsilon_{r}$ and $\gamma$. We fixed the non-Hermiticity at $\gamma=\gamma^{(0)}=1.2$, and calculate the eigenvalues of the Hamiltonian $H$ as a function of $\varepsilon_{r}$. The real and imaginary parts of the eigenfrequencies are shown respectively in Figs. 5(a) and 5(b) by solid lines. Ordered from lower frequency to higher frequency, the 1${}^{\text{st}}$ (lowest), 2${}^{\text{nd}}$ (middle), and 3${}^{\text{rd}}$ (highest) bands are denoted by green, olive and black lines respectively. We also plot the numerical results calculated directly using COMSOL with PML boundary conditions by open circles. The blue circles denote the states in exact $\mathcal{PT}$ symmetry phase (the eigenvalues are purely real) and the red circles denote the states in broken $\mathcal{PT}$-symmetry phase (the eigenvalues are complex conjugate pairs). The solid lines show excellent agreement with the open circles, indicating the validity of our non-Hermitian $3\times 3$ model Hamiltonian.

The 2${}^{\text{nd}}$ and 3${}^{\text{rd}}$ bands merge together at $\varepsilon_{r}=12.212$ and form an EP marked by $S_{23}^{L}$ (green dot) in Fig. 5(a). The 1${}^{\text{st}}$ band and the 2${}^{\text{nd}}$ band merge at $\varepsilon_{r}=12.423$, and form an EP marked by $S_{12}^{L}$ (yellow square). The subscript index $mn$ of $S_{23}^{L}$ and $S_{12}^{L}$ denotes the index of bands forming the EPs. Using the $3\times 3$ Hamiltonian, we can analyze the EPs in the $\varepsilon_{r}$ parameter space in detail. The eigenstates of the non-Hermitian Hamiltonian matrix become defective at the EP, characterized by a vanishing phase rigidity which is defined as ${r_{m}}=\left\langle{v_{{\bf{k}}m}}\right|\left.{u_{{\bf{k}}m}}\right\rangle/\left\langle{u_{{\bf{k}}m}}\right|\left.{u_{{\bf{k}}m}}\right\rangle$, where $\left|{u_{{\bf{k}}m}}\right\rangle$ and $\left|{v_{{\bf{k}}m}}\right\rangle$ are the right and left eigenstates Bulgakov et al. (2006). We plot the phase rigidity as a function of $\varepsilon_{r}$ in Fig. 5(c). The vanishing phase rigidities, $r_{m}\to 0$, also confirms the existence of the EPs.

Using the Hamiltonian model, we described the EPs in the parameter space of $\varepsilon_{r}$ for a fixed $k_{y}a/2\pi=0.05$ in Fig. 5. Now we will show the turning points denoted by a green dot and a yellow square in Fig. 4(b) and 4(d) are EPs in the parameter space of $k_{y}$. In Fig. 4(b), where $\varepsilon_{r}=12$, the lowest band emerges from the Brillouin zone center with a negative group velocity and the middle and highest bands form an EP $S_{23}^{L}$. In Fig. 5, the EP $S_{23}^{L}$ appears at $\varepsilon_{r}=12.212$ meaning that we can find an EP $S_{23}^{L}$ at $k_{y}a/2\pi=0.05$ for PC with $\varepsilon_{r}=12.212$. Similarly, in Fig. 5 an EP $S_{12}^{L}$ appears at $\varepsilon_{r}=12.423$ meaning that we can find EP $S_{12}^{L}$ at $k_{y}a/2\pi=0.05$ for PC with $\varepsilon_{r}=12.423$. In Fig. 4(d), where $\varepsilon_{r}=12.6$, the lower two bands form an EP $S_{12}^{L}$, and the highest band emerges from the Brillouin zone center with a positive group velocity.

To better describe the movement of EPs in $\varepsilon_{r}$ and $k_{y}$ parameter spaces, in Fig. 6, we plot the interface bands of PC with (a-b) $\varepsilon_{c}=12.34\pm 1.2i$, (c-d) $\varepsilon_{c}=12.3553\pm 1.2i$, and (e-f) $\varepsilon_{c}=12.36\pm 1.2i$, respectively. From Figs. 6(a) and 6(c), we see that as $\varepsilon_{r}$ increases, the middle band get closer to the lowest band and they merge together at $\varepsilon_{r}=12.3553$. This touching point denoted by $\beta_{12}$ is an anisotropic EP, which we will discuss in next section. As $\varepsilon_{r}$ increases further, as shown in Fig. 6(e), the EP $S_{23}^{L}$ disappears and an EP $S_{12}^{L}$ formed by the lower two bands comes out from the touching point $\beta_{12}$. Figure 7 describes this process. We can see that the touching point $\beta_{12}$ first split into two EPs $S_{12}^{N}$ and $S_{12}^{L}$, and then the EPs $S_{12}^{N}$ and $S_{23}^{L}$ coalesce into an order-3 EP $\Pi_{123}$, where three bands coalesce Hodaei et al. (2017). In summary, Figs. 6 and 7 describe the transition from $S_{23}^{L}$ to $S_{12}^{L}$, and show that as $\varepsilon_{r}$ changes, the evolution of the interface bands is related to the development of EPs.

VI Coalescence of exceptional points

In the above section, we build a $3\times 3$ non-Hermitian Hamiltonian model to analyze the dispersion of interface states of $\mathcal{PT}$-symmetric photonic crystals as functions of the material parameters $\varepsilon_{c}$, and describe the movement of EPs in parameter space $\varepsilon_{r}$ and $k_{y}$. But we did not take the change of the non-Hermiticity $\gamma$ into account. In this section, we expand the parameter space and study systematically the evolution of interface states as $\varepsilon_{r}$ and $\gamma$ changes. Interface bands calculated using the MS method with open boundary conditions are plotted in Fig. 8 for different parameters $\varepsilon_{c}=\varepsilon_{r}\pm i\gamma$, with $\varepsilon_{r}$ varying from 11 to 13, and $\gamma$ varying from 1 to 4. The interface bands show ziz-zag or closed-loop dispersion that cannot be seen in Hermitian systems. The turning points of the bands are EPs which are labeled by the symbols $S_{23}^{L}$, $S_{12}^{L}$ and $S_{12}^{R}$ in the figure panels. In this following, we will show that the evolution of the interface bands is closely associated with the coalescence of the EPs as we change the parameters $\varepsilon_{r}$ and $\gamma$. Such coalescence behavior is further highlighted by tracing EPs in the parameter space $(k_{y},\gamma)$, as shown in Fig. 9. The last column of panels in Fig. 8 shows that when gain/loss becomes large, one isolated band of interface states will always persists.

VI.1 Formation of the order-3 EPs in $\varepsilon_{r}=11$

In the first row of panels in Fig. 8, the real part of the relative permittivity of the cylinders is fixed at $\varepsilon_{r}=11$. As shown in Fig. 8(a), where $\gamma=1$, the band of interface states starts from the Brillouin zone center and exhibits a ziz-zag dispersion, turning around twice at EPs denoted by $S_{23}^{L}$ and $S_{12}^{R}$. As we increase the non-Hermitian strength $\gamma$, the EPs $S_{23}^{L}$ and $S_{12}^{R}$ will get closer to each other and eventually disappear as shown in Fig. 8(b-c). This process is also illustrated in Fig. 9 by the blue line, which traces the movement of EPs in the parameter space $(k_{y},\gamma)$. We can see that as $\gamma$ increases, the EPs $S_{23}^{L}$ marked by dots at small $k_{y}$ and EPs $S_{12}^{R}$ marked by open triangles at large $k_{y}$ get closer and eventually merge into an order-3 EP labeled as $\Pi_{123}$ (marked by a pink star) Demange and Graefe (2012). The EPs $S_{23}^{L}$ and $S_{12}^{R}$ are order-2 EPs, where two bands coalesce. At the order-3 EP $\Pi_{123}$, three bands coalesce.

To further illustrate the formation of the order-3 EP $\Pi_{123}$, we plot in Fig. 10 the interface bands calculated using COMSOL with PML boundary conditions applied on the $x$-direction. As shown in Figs. 10(a) and 10(d), when $\varepsilon_{c}=11\pm 3.2i$, there are two typical order-2 EPs $S_{23}^{L}$ (blue dot) and $S_{12}^{R}$ (open triangle) formed by two bands. The dispersion relations for the real and imaginary branches are of the square-root form near the EPs. As we increase the non-Hermiticity to $\gamma=3.3$, as shown in Figs. 10(b) and 10(e), the EPs $S_{23}^{L}$ and $S_{12}^{R}$ merge into an order-3 EP $\Pi_{123}$, where three bands coalesce. The EP $\Pi_{123}$ disappears upon a further increase in $\gamma$. As shown in Figs. 10(c) and 10(f), when $\gamma=3.4$, there is only one interface state with real eigenvalues (blue open circles). Other interface states (red open circles) will acquire larger and larger imaginary parts and merge into the continuum of propagating waves that are not localized at the gain-loss interface.

VI.2 Formation of the anisotropic EPs in $\varepsilon=13$

In the third row of panels in Fig. 8, the real part of the relative permittivity of the cylinders is increased to $\varepsilon=13$. In Fig. 8 (g), we see that the lower two bands form a closed-loop, with the minimum and maximum $k_{y}$ pinned by the EPs $S_{12}^{L}$ and $S_{12}^{R}$ respectively. From Figs. 8(g-i), we observe that when $\gamma$ increases, the two EPs $S_{12}^{L}$ and $S_{12}^{R}$ of the loop tend to approach each other and finally disappear. The disappearance of the loop is described in Fig. 11. EPs $S_{12}^{L}$ and $S_{12}^{R}$ coalesce at a specific $k_{y}$ point and form a new EP denoted by $\alpha_{12}$. We also plot the trace of EPs in the parameter space $(k_{y},\gamma)$ in Fig. 9 by purple line. As $\gamma$ increases, the EPs $S_{12}^{L}$ labeled by solid squares and EPs $S_{12}^{R}$ labeled by open triangles merge into a new EP labeled by $\alpha_{12}$.

It is interesting to note that EP $\alpha_{12}$ is anisotropic Ding et al. (2018). When an anisotropic EP is approached from different directions in the parameter space, it shows different singular behaviors. As shown in Fig. 11(e), the imaginary parts of the eigenfrequencies are linear functions of $k_{y}$ near the EP $\alpha_{12}$. In Figs. 12 (a-c), we plot the eigenfrequencies and phase rigidities of the interface states as a function of $\gamma$ in the vicinity of the EP $\alpha_{12}$. When the EP $\alpha_{12}$ is approached along the $\gamma$ direction, the dispersion shows a square-root behavior, in contrast to the linear behavior along $k_{y}$ direction shown in Fig. 11(e).

VI.3 Order-3 EPs and anisotropic EPs in $\varepsilon_{r}=12$

In Fig. 8, the first row of panels shows the formation of an order-3 EP $\Pi_{123}$, and the third row shows the formation of an anisotropic EP $\alpha_{12}$. We now turn to the second row of panels in Fig. 8, where $\varepsilon_{r}=12$ and find that in addition to the formation of an order-3 EP, there is also formation of anisotropic EPs. We plot the trace of EPs in the parameter space $(k_{y},\gamma)$ in Fig. 9 by green line. When $\gamma$ is small, there are two typical EPs, one is the EP $S_{23}^{L}$ labeled by solid dots at small $k_{y}$ and the other is EP $S_{12}^{R}$ labeled by open triangles at large $k_{y}$. The corresponding interface state bands with two EPs and a zigzag dispersion is shown in Fig. 8(d). As shown in the right inset of Fig. 9, when we keep on increasing $\gamma$, a new EP denoted by $\beta_{12}$ appears, which we will discuss later, and then it splits into two EPs labeled as $S_{12}^{N}$ (marked by open diamond) and $S_{12}^{L}$ (marked by solid square). When $\gamma$ increases further, the EPs $S_{12}^{N}$ and $S_{23}^{L}$ will coalesce into an order-3 EP labeled as $\Pi_{123}$, which is similar to the formation of $\Pi_{123}$ in the blue lines. This EP $\Pi_{123}$ disappears as $\gamma$ increases further, and there are only two EPs $S_{12}^{L}$ and $S_{12}^{R}$ left. These two EPs will coalesce into an anisotropic EP $\alpha_{12}$ as $\gamma$ increases, which is the same as the process we have illustrated in $\varepsilon_{r}=13$ denoted by the purple line.

Comparing these three lines in Fig. 9, we find that the green line ($\varepsilon_{r}=12$) is more than just a composite of the blue line ($\varepsilon_{r}=11$) and the purple line ($\varepsilon_{r}=13$). In the green line, we find the formation of an order-3 EP $\Pi_{123}$ which also appears in the blue line, and the formation of an anisotropic EP $\alpha_{12}$ which also appears in the purple line. But the splitting of an EP $\beta_{12}$ illustrated in the inset of Fig. 9 is unique to the green line.

Note that in Figs. 6 and 7, we have described the transition from EP $S_{23}^{L}$ to EP $S_{12}^{L}$ as $\varepsilon_{r}$ increases. In Fig. 13, we plot the transition from EP $S_{23}^{L}$ to EP $S_{12}^{L}$ as $\gamma$ increases and find the middle band touches the lowest band, forming an EP denoted by $\beta_{12}$, which is shown in Fig. 14. As $\gamma$ increases further, the EP $\beta_{12}$ splits into two EPs labeled as $S_{12}^{N}$ and $S_{12}^{L}$, which is shown in Fig. 13(b). In Figs. 12(d-f), we plot the EP $\beta_{12}$ in the $\gamma$ direction for a fixed $k_{y}a/2\pi=0.0818$. As shown in Figs. 14 and 12(d-f), $\beta_{12}$ is an anisotropic EP, whose dispersion is linear along the $k_{y}$ direction and of square-root form along the $\gamma$ direction. We note that in the green line, two anisotropic EPs (denoted by $\alpha_{12}$ and $\beta_{12}$) appears. Different from EP $\alpha_{12}$, whose imaginary parts of the eigenfrequency show linear dispersion along $k_{y}$ direction, the real parts show linear dispersion around the EP $\beta_{12}$.

VII Conclusion

In this work, we studied the formation of interface modes in $\mathcal{PT}$-symmetric PCs. It is well known that interface modes will exist if the real part of permittivity changes sign across the boundary of a PC, and we see here that interface modes with real eigenfrequencies will also exist if the imaginary parts of permittivity changes sign across a boundary. Best illustrated in Fig. 8, which shows the boundary mode dispersions for various values of the real part and imaginary part of the permittivity, the interface modes show some peculiar dispersion features that are not found in Hermitian systems. These interesting features include ziz-zag dispersions with turning points being exceptional points and there are also closed-loops of boundary modes with the vertices being exceptional points. The peculiar band diagrams are quite similar to the folded bands with infinite group velocity points discussed previouslyDavanco et al. (2007); Chen et al. (2011). These infinite group velocity points in the bands can be treated as exceptional points. We note in passing that a physical system carrying folded band has to be non-Hermitian. Hermitian systems (for example, those with $\varepsilon$ equal to a constant negative value) can exhibit folded bands but those systems are not compatible with causality. The trajectories of the EPs in the parameter space $(k_{y},\gamma)$ for different values of $\varepsilon_{r}$ are summarized in Fig. 9. Close to the Brillouin zone center, when the Bloch momentum $k_{y}\to 0$, the existence and the dispersion of these boundary modes can be obtained semi-analytically using effective medium theory. This is illustrated in Fig. 4 and elaborated in the associated discussion. As the magnitude of Bloch momentum of the boundary modes increases, a more elaborate Hamiltonian model gives a quantitative description of the dispersion, which shows that the dispersion and its dependence on system parameters are closely related to the formation and coalescence of EPs. The details are rather fascinating, as we observe the coalescence of order-2 EPs into higher order EPs and the formation of anisotropic EPs. The results obtained using different computational methods (COSMOL and multiple scattering) and different boundary conditions are consistent with each other. The results obtained with periodic boundary conditions, as shown in Fig. 2, deserves some additional comments. First, it shows that as the size of the supercell increases, the bulk modes acquire an imaginary part in the eigenfrequency (entering the broken phase). In $\mathcal{PT}$-symmetric system, it is usually the increase in non-Hermiticity (“$\gamma$”) that drives the system into the broken phase. The results in Fig. 1 shows that the increase of system size (the number of cylinders “$N$”) alone can take the system into the broken phase while $\gamma$ remains constant. In the large $N$ limit, only the interface modes remain in the exact phase, having zero imaginary part in its eigenfrequency and this band is localized on the loss-gain interface. This mode lies inside the continuum of the real bulk modes when $\gamma=0$, and in a certain sense these boundary modes are bound states in the continuum.

Acknowledgements.

This work is supported by the Research Grants Council, University Grants Committee, Hong Kong, through grants No. AoE/P-02/12, N_HKUST608/17 and C6013-18G and by the National Natural Science Foundation of China (Grant Nos. 11761161002). K.D. acknowledges funding from the Gordon and Betty Moore Foundation. We would like to thank Prof. Zhao Qing Zhang, Prof. Shubo Wang and Dr. Ruoyang Zhang for helpful discussions.

Appendix A Multiple scattering method

In this appendix, we will construct the plane wave scattering problem of a 2D PC. We can determine the band structure of a PC from the solution of an eigenvalue problem Botten et al. (2001). We consider a 2D PC comprising a rectangular lattice (lattice constant is $a$ along $y$-direction and $b$ along $x$-direction) of dielectric cylinders embedded in the air whose relative permittivity is equal to 1. The cylinders have a radius $r_{c}$ , relative permittivity $\varepsilon_{c}$, and relative permeability $\mu_{c}=1$.

We begin with a single infinite grating periodically arranged along $y$-direction, corresponding to one layer of this lattice. A TM polarization (electric field is along the axes of the cylinders) plane wave with wave vector $\textbf{k}_{0}$ is incident onto the grating and perpendicular to the axes of the cylinders ($z$-direction). The wave number is ${k}_{0}^{2}=(\omega/c)^{2}$, where $\omega$ is the angular frequency and $c$ is the speed of light in vacuum. The component parallel to $y$-direction is denoted by $k_{p}$, and therefore, the wave vector along $x$ direction is $k_{x}=\sqrt{k_{0}^{2}-k_{p}^{2}}$. We define the incident electric field as

[TABLE]

where $E_{0}$ is the amplitude of the incident field. Then the incident magnetic field can be written as

[TABLE]

To calculate the band structure, we need firstly establish the scattering matrices of the structure. The scattering field of the cylinders can be expanded as a summation of cylindrical functions, and then we expand these cylindrical functions as summation of plane waves at the interface of the grating. The various scattering plane waves with wave vector in the $x-y$ plane are like below,

[TABLE]

where $g_{m}=2m\pi/a$ is the reciprocal vector parallel to the grating, $q_{m}=\sqrt{{k_{0}}^{2}-|{k_{p}}-{g_{m}}{|^{2}}}$ and the sign ($\pm$) before $q_{m}$ is corresponding to a propagating (if $q_{m}$ is real) or evanescent wave (if $q_{m}$ is imaginary) along the positive or negative $x$-direction.

The full scattering matrices of the grating composed of cylinder are well established in many references, such as Botten et al. (2000). These scattering matrices are usually complicated to calculate. For a 2D PC with TM polarization, we can represent the cylinders using an out-of-plane electric monopole and two in-plane magnetic dipoles Wang et al. (2016)

[TABLE]

$E_{z}^{\text{loc}}$ and ${\textbf{H}}^{\text{loc}}$ are the local fields at the position of the cylinders. The local fields are summations of the external incident fields and the scattering field from other cylinders. $\alpha^{e}$ and $\alpha^{m}$ are the polarizability of the monopole and dipole of a cylinder, respectively, which are represented as Bohren and Huffman (1998)

[TABLE]

where $\beta_{0}$ and $\beta_{1}$ are Mie scattering coefficients and can be represented using Bessel functions like below

[TABLE]

$J_{n}$, $H_{n}^{(1)}$, ${J^{\prime}}_{n}$, ${H^{\prime}}_{n}^{(1)}$ are the first kind of Bessel functions, the first kind of Hankel functions and their derivatives with order $n$. $n_{c}=\sqrt{\varepsilon_{c}\mu_{c}}$ is the refractive index of the cylinder and $r_{c}$ is the radius of the cylinder. At the limit where $k_{0}r_{c}$ is a small number, the scattering properties of the cylinder are well represented by the lowest three modes. The electric fields induced by electric moment $P\hat{z}$ and magnetic moment M are represented using Green’s function like below,

[TABLE]

$\mathop{\bf{G}}\limits^{\leftrightarrow}\left({{\bf{r}},{\bf{r^{\prime}}}}\right)$ is a dyadic Green’s function, representing the response at r induced by a point source at $\textbf{r}^{\prime}$.

Now let us turn to the PC illustrated in Fig. 1(a). We denote the field acting on the polarizable unit locating at $i$th column and $j$th row as ${E}_{ij}^{\text{loc}}\hat{z}$ and $\textbf{H}_{ij}^{\text{loc}}$ , $i=1,2,\cdots,N_{1}+N_{2}$, with $N_{1}$, $N_{2}$ being the number of column layers of the left and right sublattices. According to Eq. (33), the induced monopole/dipole moments are ${{\bf{P}}_{ij}}=\alpha_{i}^{e}{{\bf{E}}_{ij}^{\rm loc}}$, ${{\bf{M}}_{ij}}=\alpha_{i}^{m}{{\bf{H}}_{ij}^{\rm loc}}$. The polarizability $\alpha_{i}^{e}$, $\alpha_{i}^{m}$ representing the cylinders in different columns could be different. The PC is uniform along $z$ direction, and therefore the fields are also uniform along $z$ direction. The electric monopole moment of the cylinder is polarized along $z$ direction, and the scattered electric fields of the cylinders only have $z$ component. The dyadic green function in Eq. (36) just keep $G_{zz}$ component as below

[TABLE]

Then the scattering electric fields induced by a cylinder at $\textbf{r}=\textbf{r}_{ij}$ can be written as

[TABLE]

The magnetic field can be obtained by Maxewell’s equation like below

[TABLE]

We know that the system is periodic along $y$-direction, applying the Bloch condition along $y$-direction into the polarization,

[TABLE]

where $P_{i}$ and $\textbf{M}_{i}$ represent the polarization of cylinder at $i$th column and [math]th row, and $a_{j}=r_{ij}-r_{i0}$ represents the distance between cylinders at $r_{ij}$ and $r_{i0}$. Then the total scattering field of the structure can be written as

[TABLE]

where we have defined

[TABLE]

representing the Green’s function of the $i$th column and $x_{i0}$ is the $x$ coordinate of the cylinder at ${{\bf{r}}_{i0}}=\left({{x_{i0}},0}\right)$. In the last step of Eq. (42), we use the identity for Hankel function to transform the lattice sum in real space into the reciprocal space Gómez-Medina et al. (2006); Botten et al. (2000).

The local field acting on the cylinder at the position $\textbf{r}_{ij}$ is given by the summation of the external incident field and the scattering field induced by other cylinders at $\textbf{r}\neq\textbf{r}_{ij}$, that is,

[TABLE]

where

[TABLE]

The first term in Eq. (43) represents the incident field, the second term represents the scattering fields induced by the column $i$, and the last term represents the scattering fields of other columns at $i^{\prime}\neq i$. $G_{b}$ in Eq. (44) represents the total scattering field at $\textbf{r}_{ij}$ of $i$th column except for the cylinder at $\textbf{r}_{ij}$, which is expressed as Gómez-Medina et al. (2006)

[TABLE]

where $\gamma_{E}$ is the Euler constant.

Substituting ${P_{i}}=\alpha_{i}^{e}E_{i}^{{\rm{loc}}}$ and ${{\bf{M}}_{i}}=\alpha_{i}^{m}{\bf{H}}_{i}^{{\rm{loc}}}$ into Eq. (43), we can build the self-consistent equations like below

[TABLE]

$F_{n}(i,i^{\prime})\quad(n=1,2\cdots 6)$ represents the lattice sums of the cylinder in column $i^{\prime}$ along the $y$-direction and take expressions as below,

[TABLE]

In the above equations, $m$ is the diffraction order, and the sum over $m$ is truncated to certain value $m_{0}$ in numerical calculation. $l-l^{\prime}$ denotes the distance between different column, and sign($l-l^{\prime}$) is the sign function. $a$ and $b$ are the lattice constant along $y$- and $x$-direction, respectively. For a square lattice that we studied in this paper, we have $a=b$.

Equation (46) describes a scattering problem and can be rewritten as below,

[TABLE]

where we denote the content in the first bracket of Eq.(46) as a matrix $M\left(\omega,k_{p}\right)$, ${\Psi^{\rm loc}}$ represents incident wave from outside and ${\Psi^{\text{loc}}}$ represents local field incident on the cylinder.

Appendix B The left eigenstate of a non-Hermitian PC with reciprocal media

In analogy with the inner product of two wave functions in quantum mechanicsSakurai et al. (2014), we can define the inner product of two vector fields $\textbf{F}(\textbf{r})$ and $\textbf{G}(\textbf{r})$ as

[TABLE]

The Hermitian adjoint of an operator $\Theta$ is labeled as $\Theta^{\dagger}$, and is defined as Joannopoulos et al. (2008)

[TABLE]

We say the operator $\Theta$ is Hermitian if $\Theta^{\dagger}=\Theta$.

In PCs, the eigenvalue problem for electric field can be written as Sakoda (2004)

[TABLE]

where $\mathcal{L}\equiv\varepsilon^{-1}\nabla\times\nabla\times$. The Hermitian adjoint of the operator $\mathcal{L}$ is

[TABLE]

For a reciprocal medium $\varepsilon^{T}=\varepsilon$, we have

[TABLE]

In a periodic system, Bloch state ${\bf{E}}_{\bf{k}}^{R}\left({\bf{r}}\right)={e^{i{\bf{k}}\cdot{\bf{r}}}}{\bf{u}}_{\bf{k}}^{R}\left({\bf{r}}\right)$ is a common eigenstate of $\mathcal{L}$ and the lattice translation operator $T_{\mathbf{a}}$ ,

[TABLE]

For the non-Hermitian system, we define the corresponding left eigen state as

[TABLE]

Taking complex conjugate of Eq. (54) and multiply $\varepsilon^{-1}$ on both sides, we obtain

[TABLE]

Therefore, ${\varepsilon^{-1}}{\left({{\bf{E}}_{\bf{k}}^{L}}\right)^{*}}$ is also a right eigenstate of $\mathcal{L}$ due to the reciprocity of the permittivity. Since $\left[{{{\cal L}^{\dagger}},{T_{\bf{a}}}}\right]=0$, the left eigenstate $\mathbf{E}_{\mathbf{k}}^{L}$ should be an eigenstate of $T_{\mathbf{a}}$, i.e. $T_{\mathrm{a}}\mathbf{E}_{\mathrm{k}}^{L}=\lambda(\mathbf{a})\mathbf{E}_{\mathrm{k}}^{L}$. And using the normalization relation of the Bloch states,

[TABLE]

we can obtain $\lambda(\mathbf{a})=e^{i\mathbf{k}\cdot\mathbf{a}}$, which means that $\mathbf{E}_{\mathrm{k}}^{L}$ is a Bloch state at k. We can check that

[TABLE]

which means $\varepsilon^{-1}\left(\mathbf{E}_{\mathrm{k}}^{L}\right)^{*}=\mathbf{E}_{-\mathrm{k}}^{R}$ is the right Bloch state at $-\mathbf{k}$, and therefore,

[TABLE]

The biorthonormal relationship of a non-Hermitian PC becomes

[TABLE]