EP restoration and fast-light edge states in photonic crystal waveguide
with glide and time reversal symmetry

Takahiro Uemura Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan NTT Basic Research Laboratories, NTT Corporation, 3-1 Morinosato-Wakamiya, Atsugi-shi, Kanagawa 243-0198, Japan Taiki Yoda Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan Yuto Moritake Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan Shutaro Otsuka Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan NTT Basic Research Laboratories, NTT Corporation, 3-1 Morinosato-Wakamiya, Atsugi-shi, Kanagawa 243-0198, Japan Kenta Takata NTT Basic Research Laboratories, NTT Corporation, 3-1 Morinosato-Wakamiya, Atsugi-shi, Kanagawa 243-0198, Japan Nanophotonics Center, NTT Corporation, 3-1, Morinosato-Wakamiya, Atsugi-shi, Kanagawa 243-0198, Japan Masaya Notomi Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan NTT Basic Research Laboratories, NTT Corporation, 3-1 Morinosato-Wakamiya, Atsugi-shi, Kanagawa 243-0198, Japan Nanophotonics Center, NTT Corporation, 3-1, Morinosato-Wakamiya, Atsugi-shi, Kanagawa 243-0198, Japan [email protected]

Abstract

Exceptional points (EPs) in the propagation states give rise to the emergence of intriguing properties with the divergence of the group velocity. However, there have been no experimental reports due to the necessity of maintaining high levels of fabrication precision and the requisite high group velocity contrast. In our study, we propose a design of photonic crystal waveguide with glide and time reversal symmetry, and derive an effective Hamiltonian for edge states to realize fast-light edge states. We adopt a systematic method to generate EPs in edge states by introducing non-Hermitian perturbations to Dirac points guaranteed by glide symmetry, which ensures that EP modes are free from out-of-plane radiation losses. Then, our study reveals the conditions for the exact EP restoration and provides an analytical solution to offset the EP smoothing due to symmetry breaking, which drastically reduces the group velocity contrast. A good symmetry property of the photonic crystal waveguide allows us to derive the effective Hamiltonian as a simple form, and the EPs can be restored by adjusting the real part of the permittivity. Furthermore, we design a feasible photonic crystal slab waveguide incorporating graphene as the absorbing material, and numerically demonstrate a group velocity reaching $v_{g}=3.3c$ near the EP, which is up to 25 times that of the original structure. Thanks to the short periodicity of photonic crystals, it’s possible to reach the speed of light in vacuum with group velocity contrasts on the order of one digit. Our study paves an innovative way to manipulate the group velocity of light.

^†^†preprint: APS/123-QED

I Introduction

Non-Hermitian optical systems have gained interest as a new approach for controlling light. PT-symmetric systems, which exhibit symmetry in spatial parity and time inversion, have received considerable attention [1, 2]. Above all, PT-symmetric optical systems are significant for implementing complex potentials through the use of optical gain from lasers [3] or current injection [4], coupled with loss from radiation or material absorption [5, 6, 7]. PT-symmetric systems have been used to realize various optical phenomena, including non-reciprocal transmission [8], single-mode lasing [9, 10] and loss-induced transparency [11].

In the context of periodic potential structures, it has been noted that the slope of the band diverges in the vicinity of the exceptional point (EP) [12]. The control of light propagation stands as a fundamental concern in photonics, and the manipulation of the group velocity of light is a central technology for the advancement of optical devices. Recent advancements in semiconductor microfabrication technology have enabled the realization of on-chip slow light in nanophotonics [13], including resonator-type coupled resonator optical waveguides (CROW) [14, 15, 16] and photonic crystal waveguide [17, 18, 19]. Conversely, light propagating with a group velocity faster than the speed of light in vacuum, known as fast light, has been observed in media with anomalous dispersion near an absorption line [20], atomic systems [21, 22], as well as in artificial media like metamaterials [23, 24]. However, achieving fast light in conventional transparent media for on-chip waveguides remains a challenge. The concept of PT-symmetric optical waveguide system harbors the potential to bring about a breakthrough in achieving fast light in on-chip waveguides.

While such interesting properties in PT-symmetric optical waveguide system have been predicted, the realization of fast light in PT-symmetric systems remain a challenge. The anomalous group velocity near EPs, referred to as fast light, has been discussed in Ref. [25, 26], and Ref. [27] proposes and theoretically investigates fast light in PT-symmetric systems and implements it numerically in CROW structures using photonic crystal cavity arrays. However, the low group velocity and strong light confinement make it difficult to apply CROW to fast-light propagation states. The initial group velocities of the mode in CROW are inherently slow (e.g. $v_{g}=0.02c$ with the periods of cavities 2.1 $\mu$ m [27]). That is, to achieve group velocities exceeding the speed of light in a vacuum requires about an order of $100v_{g}$ , which requires extremely precise control of spatial structures as well as gains and losses. Moreover, the strong light confinement in the cavity makes it increasingly difficult to achieve high fabrication precision. The issues related to group velocity and fabrication precision can be resolved by weakening the light confinement to increase the coupling rate between resonators, even when the periods of cavities is large. However, this approach is contrary to the typical CROW configuration. When we pursue CROWs with wide periods and weak light confinement, the CROW structure can be regarded as a photonic crystal with a large period. In this case, we would be using very high-order bands as a photonic crystal, which is not advantageous compared to conventional PT-symmetric photonic crystal waveguides. In addition to these, it should be noted that the waveguide modes of CROW are positioned above the light line, and their influence depends on the Q factor. That is, it is inherently difficult to eliminate the effects of out-of-plane radiation losses from the waveguide modes of CROW.

Therefore, our research motivation lies in realizing superluminal states with a group velocity contrast of one digit from the initial edge states by exploiting the PT symmetry in photonic crystal waveguide systems, which have much shorter periods of the order of sub-micrometers. The photonic crystal waveguides offer additional advantages over CROW in terms of scalability and device integration, making them favorable for applications in optical devices. However, adopting photonic crystal band systems presents challenges not encountered in CROW, where the smoothing of EPs due to finite Q can be canceled using a tight-binding-based formalism [28]. When creating EPs in photonic crystal systems, unlike in tight-binding systems, it is not known how to restore the EPs when they are smoothed by symmetry breaking. Therefore, the establishment of a method to restore EPs in photonic crystal systems is a significant step towards the realization of superluminal propagation states. Since the use of photonic crystal bulk modes as waveguide modes is not realistic, the focus of our study is on edge states in photonic crystal waveguides. Several methods for realizing photonic crystal waveguides have been proposed, one that utilizes a degenerate point created by imaginary perturbations and unit cell enlargement [29, 30]; another that employs an accidental degenerate point [31, 30]; and a third method that uses symmetry-protected degenerate points [30]. The first and second methods make it challenging to position the EP below the light line. In contrast, the third method allows for the EP to be positioned below the light line, which assures that the propagation modes near the EP are free from out-of-plane radiation loss. Additionally, the third method’s allows for an analytical solution of the effective Hamiltonian thanks to its good symmetry properties, making it easier to develop methods to restore the EP when the symmetry is broken, which has not been noted in Ref. [30] and is a new result we present in this paper.

In this paper, we propose a photonic crystal waveguide with the symmetry under simultaneous glide and time-reversal operations. Our scheme starts with the generation of EPs using a systematic approach with PT-symmetric perturbations at Dirac points, which appear due to the spatial glide symmetry [32], the EPs appear at the boundary of the Brillouin zone, providing two key advantages, both of which address issues inherent in traditional PT-symmetric photonic crystal waveguides [29, 31, 30]. The first advantage is that the EPs can be positioned below the light line, ensuring that the propagation modes near the EP are free from out-of-plane radiation loss. This is in contrast to most structures proposed in previous studies, such as [30], where the smoothing of EP could not be avoided due to the radiation loss, which drastically reduces the group velocity contrast. The second advantage is that the symmetry at the boundary of the Brillouin zone allows the $\bm{k}\cdot\bm{p}$ Hamiltonian around the EPs to be represented in a simplified form, making it easier to model the response of edge states to non-Hermitian perturbations. EPs in our waveguide are sensitive to any deformation of glide and time reversal symmetry, which presents a challenge to achieving high group velocity contrast. Therefore, we derive the precise EP conditions for complex permittivity perturbations, which are essential for realizing superluminal propagation states. In addition to the above advantages, it is remarkable that we adopt the valley photonic crystal waveguide as the basic structure, which exhibits topological properties and has potential for future applications in optical circuits with bends [33]. We also mention that our proposed structure is characterized by its ease of fabrication. Our design can be implemented by selectively loading graphene along the glide plane of a photonic crystal slab waveguide, which is simpler than the embedded heterostructures that require precise gain and loss adjustments. As a whole, our study is a step towards experimental observations of superluminal propagation states in on-chip non-Hermitian nano-optical systems.

II Design of non-Hermitian photonic crystal waveguides with glide and time reversal symmetry

First, we present the design of our photonic crystal heterostructure based on two-dimensional (2D) simulation, where the system extends infinitely along the $z$ -axis. All numerical calculations in this paper were performed using finite element method (FEM), COMSOL Multiphysics.

Refer to caption — Figure 1: (a) Interface geometry of the photonic crystal hetero-structure with the triangular lattice and triangular holes. This structure has glide symmetry along $y=0$ plane as shown by the dashed line. (b) The band dispersion curves in the $k_{x}$ direction of the TE mode in geometry (a). (c) Interface geometry of the glide symmetric photonic crystal hetero-structure with balanced gain $(y>0)$ and loss $(y<0)$ . This structure is invariant under both a glide operation with the $xz$ -plane and a simultaneous time-reversal operation. (d), (e) The real and imaginary part of band dispersion curves along the $k_{x}$ direction near $k_{x}a=\pi$ .

We start from the design of a glide symmetric photonic crystal heterostructure with the triangular lattice and triangular holes as shown in Fig. 1(a). This is also known as the bearded-edge of the photonic crystal [34]. We considered also other types of lattices and hole shapes (See Appendix. A), but we adoped the triangular lattice and triangular holes because we found that Dirac points can be easily located inside the band gap and below the light line. The difference between triangular and circular holes lies in the presence or absence of valley-polarized features near the K(K’) point in the bulk band, and there is no essential difference for constraction of EPs in edge states. We set the lattice constant as $a=470$ nm, the triangular hole side as $s=0.8a$ . The air and slab permittivity are set to $\varepsilon_{\mathrm{air}}=1$ and $\varepsilon_{\mathrm{slab}}=6.76$ , respectively. The value of $\varepsilon_{\mathrm{slab}}$ corresponds to an effective refractive index of typical semiconductor materials such as silicon and compound semiconductor materials. The dashed line in Fig. 1(a) represents the glide symmetric plane, and we define it as $y=0$ and establish our coordinate axis accordingly. The heterostructure is invariant under a glide operation with the $xz$ -plane, which is a combination of a reflection with the $xz$ -plane and a translation along the $x$ -axis. The glide operator is given by

\hat{G}=\left\{m_{y}\middle|\frac{1}{2}\hat{\mathrm{x}}\right\}

(1)

where $m_{y}$ and $\hat{\mathrm{x}}$ represent the mirror reflection with respect to the $xz$ -plane and a unit vector along the $x$ axis, respectively. Let $\varepsilon(\bm{r}),\ \bm{r}={[x,y,z]}^{\top}$ be the permittivity of initial photonic crystal, and $\varepsilon(\bm{r})$ is invariant with respect to the glide operation:

\hat{G}\varepsilon(\bm{r})=\varepsilon\left(x+\frac{a}{2},-y,z\right)=\varepsilon(\bm{r}).

(2)

The band dispersion curves of TE mode along the $k_{x}$ -axis are shown in Fig. 1(b). The region shaded in gray represents the area outside the bulk band and beyond the light line. The existence of Dirac point at $k_{x}a=\pi$ in the band dispersion is guaranteed by the spatial glide symmetry of the photonic crystal [Söllner2015, 36, 18]

The key point of our design is to introduce a gain on the upper side and a loss of the same magnitude on the lower side, bordering the glide plane. Figure 1(c) shows the interface geometry of the photonic crystal heterostructure with balanced gain $(y>0)$ and loss $(y<0)$ . The perturbation of imaginary part of permittivity are set to $\Delta\varepsilon_{\mathrm{i}}(\bm{r})=0.1\times 1_{\bm{r}\in\mathrm{slab}\ \cap\ y>0}$ (gain) and $\Delta\varepsilon_{\mathrm{i}}(\bm{r})=0.1\times 1_{\bm{r}\in\mathrm{slab}\ \cap\ y<0}$ (loss). Let $\hat{T}$ be the time-reversal operator, then the perturbed permittivity $\varepsilon^{\prime}(\bm{r})=\varepsilon(\bm{r})+\Delta\varepsilon(\bm{r})$ has the symmetry with respect to the simultaneous operation of $\hat{G}$ and $\hat{T}$ :

\displaystyle(\hat{G}\hat{T})\varepsilon^{\prime}(\bm{r})=(\hat{T}\hat{G})\varepsilon^{\prime}(\bm{r})=\varepsilon^{\prime}(\bm{r})

(3)

Figures 1(d) and (e) show the real and imaginary part of the band dispersion near $k_{x}a=\pi$ , respectively. The Dirac point splits into two EPs when gain and loss are balanced, thank to the glide and time reversal (GT) symmetry. Note that the glide operation on our system corresponds to the parity inversion operation in the context of PT-symmetric optical systems. The edge states undergo a phase transition between PT-symmetric and PT-broken phases, just as in the general PT-symmetric optical systems. See Appendix. B for the details with the $H_{z}$ distribution of the edge state. Our approach is advantageous because it allows for the systematic generation of EPs in the band dispersion of waveguide edge states. In addition, the initial Dirac point is located in the band gap, which allows us to utilize for waveguide applications.

III Reconstruction of exact EPs based on perturbation theory

In the previous section we discussed the situation with glide and time reversal symmetry. We need to mention that the EP is sensitive to symmetry breaking; any detuning that breaks glide and time-reversal symmetry lifts the EP and smoothens the band dispersion near the EP. For example, as shown in Fig. 2(a), when a detuning of the real part of the permittivity, $\mathrm{Re}(\Delta\varepsilon)=0.02$ , is applied to the region of the dielectric where $y<0$ , the EP lifts and the band dispersion curves near the initial EP are smoothed. In addition, as shown in Fig. 2 (b), introducing only loss in the dielectric region where $y<0$ also breaks the glide and time-reversal symmetry due to the imbalance of gain and loss. The “passive” systems with only loss are advantageous for the simple fabrication processes [37], but it is known that the EPs are smoothed due to the the imaginary coupling of each eigenstate [28]. In any case, the smoothing of EPs drastically reduces the group velocity contrast, which is a critical issue for realizing superluminal propagation states. Therefore, it is essential to develop methods to restore the EP by counteracting the effects of symmetry breaking. While Ref. [28] provides theoretical analysis and methods to offset the smoothing of EPs in CROW systems, the situation in the photonic crystal waveguide systems has not been investigated so far.

Here, we derive an effective Hamiltonian for using the $\bm{k}\cdot\bm{p}$ perturbation method to model the bands of edge states near the original Dirac point, i.e., $k_{x}a=\pi$ . Furthermore, based on the $\bm{k}\cdot\bm{p}$ Hamiltonian, we derive the precise EP conditions for complex permittivity perturbations, and show that EPs can always be restored even in passive photonic crystal waveguides by adjusting of the real part of the permittivity.

III.1 Effective Hamiltonian of edge states near the Dirac point

This part describes the derivation of the effective Hamiltonian which models the band dispersion of edge states in the photonic crystal waveguide near the Dirac point. We consider the two types of perturbation, one is the detuning of the wave vector from the Dirac point $\bm{\delta k}=[k_{x}-k_{0},0,0]^{\top}$ , and the other denotes the permittivity perturbation $\Delta\varepsilon(\bm{r})$ . For the time being, we consider the case the permittivity perturbation is pure imaginary $\Delta\varepsilon(\bm{r})\in i\mathbf{R}$ . We define the two unperturbed bases $\bm{\varphi}_{\alpha k_{0}}^{(0)}(\bm{r})$ and $\bm{\psi}_{\alpha k_{0}}^{(0)}(\bm{r})$ as the normalized electric and magnetic fields at $k_{x}=\pi/a$ with $\Delta\varepsilon(\bm{r})=0$ . Here, $\alpha$ denotes the label of the basis, and $C_{\alpha}$ denotes the expansion coefficients of the basis. The two bases are degenerated at the Dirac point due to the glide symmetry of the unperturbed photonic crystal.

The $\bm{k}\cdot\bm{p}$ Hamiltonian for the edge states near the Dirac point is given by the following $2\times 2$ matrix form (See Appendix. C for details):

		$\displaystyle\left[\begin{matrix}{\left(\frac{\omega_{\alpha}^{(0)}}{c}\right)}^{2}+P^{\prime}_{\alpha\alpha}+Q^{\prime}_{\alpha\alpha}&P^{\prime}_{\alpha\beta}+Q^{\prime}_{\alpha\beta}\\ P^{\prime}_{\beta\alpha}+Q^{\prime}_{\beta\alpha}&{\left(\frac{\omega_{\beta}^{(0)}}{c}\right)}^{2}+P^{\prime}_{\beta\beta}+Q^{\prime}_{\beta\beta}\end{matrix}\right]\left[\begin{matrix}C_{\alpha}\\ C_{\beta}\end{matrix}\right]$
		$\displaystyle\quad\quad={\left(\frac{\omega_{nk}}{c}\right)}^{2}\left[\begin{matrix}1+G_{\alpha\alpha}&G_{\alpha\beta}\\ G_{\beta\alpha}&1+G_{\beta\beta}\end{matrix}\right]\left[\begin{matrix}C_{\alpha}\\ C_{\beta}\end{matrix}\right]$		(4)

Here, $c$ is the speed of light in vacuum, $\omega_{nk}$ represents the eigenfrequency of the $n$ -th mode, $\alpha$ and $\beta$ denote the label of the eigenmodes, $\omega_{\alpha}^{(0)},\omega_{\beta}^{(0)}$ denote the unperturbed frequencies, and $\omega_{\alpha}^{(0)}=\omega_{\beta}^{(0)}$ holds regardless of the choice of the bases due to the degeneracy by the glide symmetry. $P^{\prime}_{\alpha\beta}$ , $G_{\alpha\beta}$ and $Q^{\prime}_{\alpha\beta}$ are defined as follows:

$\displaystyle P_{\alpha\beta}^{\prime}$	$\displaystyle=\int_{\mathrm{u.c.}}\left(\omega_{\beta}^{(0)}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\times\bm{\varphi}_{\beta k_{0}}^{(0)}(\bm{r})\right.$
	$\displaystyle\qquad\qquad\left.-\omega_{\alpha}^{(0)}\bm{\varphi}_{\alpha k_{0}}^{(0)*}(\bm{r})\times\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\right)\cdot\bm{\delta k}\,\mathrm{d}^{3}r$	(5)
$\displaystyle G_{\alpha\beta}$	$\displaystyle:=\int_{\mathrm{u.c.}}\mathrm{d}^{3}r\Delta\varepsilon(\bm{r})\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})$	(6)
$\displaystyle Q^{\prime}_{\alpha\beta}$	$\displaystyle=-\int_{\mathrm{u.c.}}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\left(\bm{\delta k}\times\bm{\delta k}\times\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\right)\mathrm{d}^{3}r$	(7)

The $\bm{k}\cdot\bm{p}$ perturbation generally describes the response of frequencies and eigenstates to the perturbation of the wave vector $\bm{\delta k}$ . However, the response of the band to the perturbation $\Delta\varepsilon(\bm{r})$ is also naturally incorporated in the form of $G_{\alpha\beta}$ . Considering that $\bm{\delta k}$ only has a $k_{x}$ component, $P^{\prime}_{\alpha\beta}$ and $Q^{\prime}_{\alpha\beta}$ can be expressed using scalars $P_{\alpha\beta}$ and $Q_{\alpha\beta}$ , which are independent of the wavenumber $\delta k_{x}$ , as follows:

\displaystyle P_{\alpha\beta}^{\prime}=P_{\alpha\beta}\delta k_{x},\quad Q^{\prime}_{\alpha\beta}=Q_{\alpha\beta}{(\delta k_{x})}^{2}

(8)

The $P_{\alpha\beta}$ and $Q_{\alpha\beta}$ are given by

$\displaystyle P_{\alpha\beta}$	$\displaystyle=\int_{\mathrm{u.c.}}\left(\omega_{\beta}^{(0)}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\times\bm{\varphi}_{\beta k_{0}}^{(0)}(\bm{r})\right.$
	$\displaystyle\qquad\qquad\left.-\omega_{\alpha}^{(0)}\bm{\varphi}_{\alpha k_{0}}^{(0)*}(\bm{r})\times\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\right)_{x}\,\mathrm{d}^{3}r$	(9)
$\displaystyle Q_{\alpha\beta}$	$\displaystyle=\int_{\mathrm{u.c.}}\left({\left(\psi_{\alpha k_{0}}^{(0)*}(\bm{r})\right)}_{y}{\left(\psi_{\beta k_{0}}^{(0)}(\bm{r})\right)}_{y}\right.$
	$\displaystyle\qquad\qquad\qquad\left.+{\left(\psi_{\alpha k_{0}}^{(0)*}(\bm{r})\right)}_{z}{\left(\psi_{\beta k_{0}}^{(0)}(\bm{r})\right)}_{z}\right)\mathrm{d}^{3}r$	(10)

While $G_{\alpha\beta}$ represents the effect of permittivity perturbation on the bands, $P_{\alpha\beta}$ and $Q_{\alpha\beta}$ is determined exclusively by the unperturbed photonic crystal.

Next, we simplify Eq.(III.1) by using the mirror symmetry along the $x$ axis, which is independent of the glide symmetry. We set the unit cell as shown in Fig. 2 (a), and $x$ -axis as the center of the unit cell. Since the initial permittivity and perturbation are symmetric about the $x$ axis, $\varepsilon(\bm{r})$ and $\Delta\varepsilon(\bm{r})$ in the unit cell are invariant under the mirror operation with respect to the $yz$ -plane:

	$\displaystyle\hat{M}_{x}\left(\varepsilon(\bm{r})\right)$	$\displaystyle=\varepsilon(-x,y,z)=\varepsilon(x,y,z)$		(11)
	$\displaystyle\hat{M}_{x}\left(\Delta\varepsilon(\bm{r})\right)$	$\displaystyle=\Delta\varepsilon(-x,y,z)=\Delta\varepsilon(x,y,z)$		(12)

Here, $\hat{M}_{x}$ is the mirror operator with respect to the $yz$ -plane. Note that the assumption that the initial permittivity are symmetric about the $x$ -axis (Eq.(12)) is generally satisfied in typical photonic crystal lattice such as triangular, square, and honeycomb lattices. Then, the two eigenmodes can be classified into the following two categories according to the eigenvalues for operator $\hat{M}_{x}$ :

	$\displaystyle\hat{M}_{x}{\left(\bm{\varphi}_{ik_{0}}^{(0)}(x,y,z)\right)}_{z}$	$\displaystyle={\left(\bm{\varphi}_{ik_{0}}^{(0)}\left(-x,y,z\right)\right)}_{z}$
		$\displaystyle=\pm{\left(\bm{\varphi}_{ik_{0}}^{(0)}(x,y,z)\right)}_{z},\quad i\in\{\mathrm{u},\mathrm{d}\}$		(13)

There is an arbitrary choice of the basis functions of the eigenmodes at $k_{x}=\pi/a$ since the band dispersion curves are degenerate at the Dirac point, and we choose the basis functions $\bm{\varphi}_{\mathrm{u}k_{0}}^{(0)}(\bm{r}),\bm{\varphi}_{\mathrm{d}k_{0}}^{(0)}(\bm{r})$ , which are characterized by the mirror eigenvalues $-1$ and $1$ :

	$\displaystyle\hat{M}_{x}{\left(\bm{\varphi}_{\mathrm{u}k_{0}}^{(0)}(\bm{r})\right)}_{z}$	$\displaystyle=-{\left(\bm{\varphi}_{\mathrm{u}k_{0}}^{(0)}(\bm{r})\right)}_{z}$		(14)
	$\displaystyle\hat{M}_{x}{\left(\bm{\varphi}_{\mathrm{d}k_{0}}^{(0)}(\bm{r})\right)}_{z}$	$\displaystyle={\left(\bm{\varphi}_{\mathrm{d}k_{0}}^{(0)}(\bm{r})\right)}_{z}$		(15)

From the relationship of electric and magnetic fields of the eigenmodes, The symmetry of $H_{z},E_{x},E_{y}$ at the Dirac point is summarized in Table 1.

Table 1: The symmetry along the

y

-axis of the eigenmodes at the Dirac point. The value of 1 or

-1

represent the symmetric and antisymmetric, respectively.

	u	d
$H_{z}$	$-1$	$1$
$E_{x}$	$-1$	$1$
$E_{y}$	$1$	$-1$

According to the even or odd symmetry of the electromagnetic field components as shown in Table 1,

\displaystyle P_{\mathrm{u}\mathrm{u}}=P_{\mathrm{d}\mathrm{d}}=G_{\mathrm{u}\mathrm{d}}=G_{\mathrm{d}\mathrm{u}}=Q_{\mathrm{u}\mathrm{d}}=Q_{\mathrm{d}\mathrm{u}}=0

(16)

holds when the condition (11) and (12) are satisfied. Moreover, we get ${\left(P_{\mathrm{u}\mathrm{d}}\right)}_{x}=-{\left(P_{\mathrm{d}\mathrm{u}}\right)}_{x}=:-iP$ , $P\in\mathbf{R}$ from the time reversal symmetry of the initial structure, and $Q_{\mathrm{u}\mathrm{u}}=Q_{\mathrm{d}\mathrm{d}}=:Q$ from the mirror symmetry. Therefore, the Eq. (III.1) is reduced to the following form:

	$\displaystyle\left[\begin{matrix}F&iP\delta k_{x}\\ -iP\delta k_{x}&F\end{matrix}\right]\left[\begin{matrix}C_{\mathrm{u}}\\ C_{\mathrm{d}}\end{matrix}\right]$
	$\displaystyle\qquad={\left(\frac{\omega_{nk}}{c}\right)}^{2}\left[\begin{matrix}1+G_{\mathrm{u}\mathrm{u}}&0\\ 0&1+G_{\mathrm{d}\mathrm{d}}\end{matrix}\right]\left[\begin{matrix}C_{\mathrm{u}}\\ C_{\mathrm{d}}\end{matrix}\right],$		(17)
	$\displaystyle\qquad F:={\left(\frac{\omega_{D}}{c}\right)}^{2}+Q{(\delta k_{x})}^{2}$		(18)

Here, $\omega_{D}$ is the eigenfrequency of two eigenmodes at the Dirac point, and $P$ corresponds to the slope of the band dispersion curves at $k_{x}a=\pi$ . Furthermore, multiply the inverse of the diagonal matrix from the left to eliminate the matrix from the right side of the equation, and let

	$\displaystyle\Delta G$	$\displaystyle:=\frac{G_{\mathrm{d}\mathrm{d}}-G_{\mathrm{u}\mathrm{u}}}{2}$
		$\displaystyle=\frac{1}{2}\int_{\mathrm{u.c.}}\mathrm{d}^{3}r\Delta\varepsilon(\bm{r})\left({\left\|\bm{\psi}_{\mathrm{d}k_{0}}^{(0)}(\bm{r})\right\|}^{2}-{\left\|\bm{\psi}_{\mathrm{u}k_{0}}^{(0)}(\bm{r})\right\|}^{2}\right)$		(19)

then we get the following form.

		$\displaystyle\left(F\left(1+\frac{G_{\mathrm{u}\mathrm{u}}+G_{\mathrm{d}\mathrm{d}}}{2}\right)I_{2}+H^{\prime}\right)\left[\begin{matrix}C_{\mathrm{u}}\\ C_{\mathrm{d}}\end{matrix}\right]=E_{nk}\left[\begin{matrix}C_{\mathrm{u}}\\ C_{\mathrm{d}}\end{matrix}\right],\$
		$\displaystyle\quad H^{\prime}=\left[\begin{matrix}F\Delta G&i(1+G_{\mathrm{d}\mathrm{d}})P\delta k_{x}\\ -i(1+G_{\mathrm{u}\mathrm{u}})P\delta k_{x}&-F\Delta G\end{matrix}\right]$		(20)

Here, $I_{2}$ is a $2\times 2$ identity matrix, and $\displaystyle E_{nk}:=(1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}}){\left(\omega_{nk}/c\right)}^{2}$ denotes the eigenvalues. $H^{\prime}$ in Eq.(III.1) characterizes the band structure of the edge states. The diagonal term $\pm{\left(\omega_{D}/c\right)}^{2}\Delta G$ in Eq. (III.1) is regarded as the difference between the potentials of $y>0$ and $y<0$ . More specifically, $\Delta G$ is proportional to the difference between two norm of the basis functions in the perturbed region. Furthermore, the off-diagonal terms $i(1+G_{\mathrm{d}\mathrm{d}})P,\ -i(1+G_{\mathrm{u}\mathrm{u}})P$ in $H^{\prime}$ represent the complex couplings of two basis. Note that it is fundamentally important for the simplicity of the Hamiltonian thanks to Eq.(16) that the basis $\mathrm{u}$ and $\mathrm{d}$ are degenerated due to the glide symmetry, and that they also possess (anti-)symmetry with respect to $\hat{M}_{x}$ .

The corresponding eigenvalues $E_{nk}$ are given by the following equation:

	$\displaystyle E_{nk}$	$\displaystyle=F\left(1+\frac{G_{\mathrm{u}\mathrm{u}}+G_{\mathrm{d}\mathrm{d}}}{2}\right)$
		$\displaystyle\ \pm\sqrt{{(F\Delta G)}^{2}+(1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}}){\left(P\,\delta k_{x}\right)}^{2}}$		(21)

The sign of the square root in Eq. (III.1) corresponds to the solution number of the eigenvalues $n=1,2$ . From the form of eigenvalues in Eq. (III.1), The corresponding EP condition is given by

\displaystyle\frac{{P(\delta k_{x})}^{2}}{{\left({\left(\frac{\omega_{D}}{c}\right)}^{2}+Q{(\delta k_{x})}^{2}\right)}^{2}}=\frac{-{(\Delta G)}^{2}}{(1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}})}

(22)

All of the components in the left-hand side of Eq. (22) are real numbers, there exists EPs in the real $\delta k_{x}$ -space when the right-hand side of Eq. (22) is positive. The right-hand side of Eq. (22) is always positive as long as the system has the glide and time reversal symmetry, since $\Delta G$ is pure imaginary and the $(1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}})$ is a real number because $G_{\mathrm{u}\mathrm{u}}^{*}=G_{\mathrm{d}\mathrm{d}}$ holds. Eq. (22) indicates that the EPs under the glide and time reversal symmetry always exist near $k_{x}=\pi/a$ , which means that EP is always located below the light line.

Although the initial form of Hamiltonian shown in Eq.(III.1) has the similar form in Ref. [38], the derived Hamiltonian shown in Eq.(III.1) and the physical interpretation are entirely different. Their model targets the bulk bands in the photonic crystals, and utilizes the band folding with double-period imaginary perturbations. The band folding places the EP above the light line, which makes it difficult to use as a waveguide mode. In contrast, our system is a photonic crystal waveguide with edge states, which allows the EP to utilize the waveguide mode. Furthermore, the design with glide symmetry allows the EP to be located below the light line, which ensures that the EP modes are free from out-of-plane radiation loss. This aspect highlights an advantage of our scheme compared to the waveguide structures proposed in Refs. [29, 31, 30].

The configuration for numerical analysis and the comparison between FEM and $\bm{k}\cdot\bm{p}$ approximation are as follows. Figures 3(a) and (b) show the computational setup and the magnetic field of the two bases, respectively. We can see that ${|\bm{\varphi}_{\mathrm{u}}(\bm{r})|}^{2}$ and ${|\bm{\varphi}_{\mathrm{d}}(\bm{r})|}^{2}$ , which are proportional to ${|\bm{H}|}^{2}$ , are localized at the upper and lower sites, respectively. Figure 3 (c) shows the $H_{z}$ , $E_{x}$ , and $E_{y}$ distribution of the bases at the Dirac point defined by Eq. (14) and (15). It is clear the $H_{z}$ and $E_{x}$ distribution of mode u are antisymmetric along the $y$ -axis, while the $E_{y}$ distribution is symmetric along the $y$ -axis. Moreover, the opposite is true for mode d. These results are consistent with the classification of the eigenmodes based on the mirror symmetry along the $x$ axis in Table 1.

The band dispersion in a heterostructure with the PT-symmetric permittivity perturbation using the FEM method and $\bm{k}\cdot\bm{p}$ Hamiltonian is shown in Fig. 4 (a). The black circle and blue triangle points represent the real and imaginary parts of the eigenvalues obtained by FEM, and the black and blue lines represent the real and imaginary parts by the $\bm{k}\cdot\bm{p}$ Hamiltonian. No fitting is used to obtain the coefficients in the $\bm{k}\cdot\bm{p}$ Hamiltonian; instead, the coefficients are determined using two unperturbed bases $\bm{\varphi}_{\alpha k_{0}}^{(0)}(\bm{r})$ and $\bm{\psi}_{\alpha k_{0}}^{(0)}(\bm{r})$ calculated by FEM and Eq. (6), (III.1) and (III.1). The band dispersion show the good agreement between the $\bm{k}\cdot\bm{p}$ and FEM methods. Moreover, the EPs are clearly observed in the real $\delta k_{x}$ -space, which is consistent with the theoretical prediction.

III.2 Restoration of EP by adjusting the real part of the permittivity

Next, we describe how to restore the EP in the case where the EP is smoothed out by breaking the glide and time reversal symmetry. It is obvious that if there is a mismatch in the real part of the permittivity on both sides of the glide plane, as shown in Fig. 2(a), the EP can be restored by simply adjusting the real part of the permittivity that compensates for the difference. Mathematically, this can be explained as follows: when there is a mismatch in the real part of the permittivity, $\Delta G$ defined in Eq. (III.1) becomes a complex number, and the EP condition shown in Eq. (22) is no longer satisfied. By compensating for the real part of the permittivity, $\Delta G$ becomes a pure imaginary number again, and the EP is restored. In contrast to the situation, the way to restore the EP in the presence of gain-loss imbalance is not trivial. The gain-loss imbalance breaks the EP condition since $G_{\mathrm{u}\mathrm{u}}^{*}\neq G_{\mathrm{d}\mathrm{d}}$ holds and $(1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}})$ becomes a complex number. The situation has the same mathematical origin as the smoothing of EP due to the complex coupling between the two basis in a coupled resonator array [28].

The following discussion will focus on the case where the uniform imaginary part of the permittivity $\Delta\varepsilon_{i}$ is introduced for the dielectric region where $y<0$ ( $V_{G}$ as shown in the blue region in Fig. 1 (b),(c)), which is responsible for the non-Hermitian perturbation. Moreover, the real part of the permittivity perturbation $\Delta\varepsilon_{r}$ is introduced for $V_{G}$ to restore the EP.

Let $I_{\alpha\beta}$ , $J_{\alpha\beta}$ , and $\kappa$ be defined as follows:

	$\displaystyle I_{\alpha\beta}$	$\displaystyle:=\int_{\bm{r}\in V_{G}}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\mathrm{d}^{3}r,$
	$\displaystyle J_{\alpha\beta}$	$\displaystyle:=\Delta\varepsilon_{\mathrm{i}}I_{\alpha\beta},\quad\kappa:=\frac{\Delta\varepsilon_{\mathrm{r}}}{\Delta\varepsilon_{\mathrm{i}}}$		(23)

From the following definition, we get $G_{\alpha\beta}=(\kappa+i)J_{\alpha\beta}$ . The necessary condition for $-{(\Delta G)}^{2}/(1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}})$ to be a positive real number, which is the condition for the existence of EP shown in Eq. (22), is that the vectors of the denominator and numerator are parallel in the complex plane, so the following conditional expression is obtained.

\displaystyle\mathrm{det}\left[\begin{matrix}\mathrm{Re}\left(-{(\Delta G)}^{2}\right)&\mathrm{Re}\left((1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}})\right)\\ \mathrm{Im}\left(-{(\Delta G)}^{2}\right)&\mathrm{Im}\left((1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}})\right)\end{matrix}\right]

\displaystyle=0

(24)

Substitute $\Delta G,\ G_{\mathrm{u}\mathrm{u}}$ and $G_{\mathrm{d}\mathrm{d}}$ into $\kappa$ and $J_{ij}$ , and then we get the relation under EP conditions:

\displaystyle(J_{\mathrm{u}\mathrm{u}}+J_{\mathrm{d}\mathrm{d}})\kappa^{2}+2\kappa+(J_{\mathrm{u}\mathrm{u}}+J_{\mathrm{d}\mathrm{d}})=0

(25)

By solving Eq. (25) for $\kappa$ and choose a physical solution $\kappa\ll 1$ , then we obtain the relationship for $\Delta\varepsilon_{\mathrm{r}}$ and $\Delta\varepsilon_{\mathrm{i}}$ at EP, which is the condition we are looking for.

\displaystyle\Delta\varepsilon_{\mathrm{r}}=-\frac{1-\sqrt{1-(I_{\mathrm{u}\mathrm{u}}+I_{\mathrm{d}\mathrm{d}})^{2}{(\Delta\varepsilon_{\mathrm{i}})}^{2}}}{I_{\mathrm{u}\mathrm{u}}+I_{\mathrm{d}\mathrm{d}}}

(26)

Here, we numerically verify the EP smoothing and restoration by adjusting the real part of the permittivity. Figure 4(b) shows the band dispersion when the imaginary part of the permittivity is $\Delta\varepsilon_{i}(r)=-0.2$ only for $y<0$ . It reveals that the EPs are smoothed out unlike the case in Fig. 4 (a). In this case, $(1+G_{\mathrm{u}\mathrm{u}})(1+G_{\mathrm{d}\mathrm{d}})$ is a complex number and EP does not exist in real $\delta k_{x}$ -space. According to Eq. (26), $\Delta\varepsilon_{\mathrm{r}}=-0.00235$ for $y<0$ is required to cancel the effects of complex coupling in the case of Fig. 4(b) with $\Delta\varepsilon_{\mathrm{i}}=-0.2$ for $y<0$ . Figure 4(c) shows the band dispersion when $\Delta\varepsilon_{\mathrm{r}}=-0.00235$ for $y<0$ is introduced, and the EP smoothing is clearly suppressed compared to Fig. 4(b). Thus, even when the EP is smoothed out by the gain-loss imbalance, appropriate adjustment of the real part of the permittivity offsets the complex coupling and restores the anomalous dispersion. In this way, even for modes in photonic crystal waveguides, where a formalism like tight-binding is not well known, the effective Hamiltonian model can elucidate the conditions for EP restoration in non-trivial gain-loss imbalance scenarios.

IV Superluminal edge states and enlarged propagation length on photonic crystal waveguide

The anomalous dispersion near the EPs leads to the high group velocity of the edge states, which allows the realization of superluminal propagation. Here, we numerically verify the group velocity of the edge states, and show that the singularity in the band dispersion near the EP is restored by the method we proposed in the previous section, resulting in a group velocity exceeding the light speed in vacuum $c$ . In addition, we discuss the fast-light effect on the propagation length of the edge states. The propagation length of the edge mode is limited by the lifetime of the edge states in the passive systems, but is expected to increase as the group velocity increases. In the following sections, we will also discuss the expected increase in propagation length based on the anomalous group velocity near the EPs.

Figure 5 (a) shows the group velocity of the edge states, which is calculated by the slope of the band dispersion in Fig. 4. The group velocity is defined as $v_{g}=\partial\omega_{R}/\partial k_{x}$ , where $\omega_{R}$ is the real part of the eigenfrequency.

I.

without permittivity adjustment
II.

with $\Delta\varepsilon=i0.1$ for $y>0$ and $\Delta\varepsilon=-i0.1$ for $y<0$ (balanced loss-gain)
III.

with $\Delta\varepsilon=-i0.2$ for $y>0$
IV.

with $\Delta\varepsilon=-0.00235-i0.2$ for $y<0$

Here, II., III. and IV. represent the cases in Fig. 4 (a), (b) and (c), respectively. When the permittivity is not adjusted (I.), the maximum group velocity of the edge states is $v_{g}=0.17c$ . Note that the initial group velocity is larger than that of CROW ( $v_{g}=0.02c$ in Ref [28]). Then, the group velocity is increased by the balanced loss-gain perturbation $\Delta\varepsilon=\pm i0.1$ (II.). The maximum group velocity of the edge states reaches the superluminal ( $v_{g}=4.5c$ ), which is determined by the finite mesh of the wavevector $\Delta k_{x}$ in the FEM calculation and should theoretically diverge at the EPs. The graph III. in Fig. 4 (a) shows the group velocity of the edge states with loss-biased perturbation $\Delta\varepsilon=-i0.2$ for $y>0$ . The maximum group velocity near the EPs in FEM calculation is subluminal ( $v_{g}=0.85c$ ), which is 5.0 times larger than that of $\Delta\varepsilon=0$ ( $v_{g}=0.17c$ ), but is limited by the EP smoothing. The graph IV. in Fig. 5 (a) shows the group velocity of the edge states with the permittivity adjustment $\Delta\varepsilon=-0.00235-i0.2$ for $y<0$ . The adjustment of the real part of the permittivity restore the EP, and the maximum group velocity of the edge states is accelerated to $v_{g}=2.8c$ , which is 16.9 times larger than that of $\Delta\varepsilon=0$ . The group velocity calculated from the FEM method is sufficiently close to the singularity predicted by the $\bm{k}\cdot\bm{p}$ method. Consequently, by restoring the EP through adjustment of the real part of the permittivity, the group velocity can be increased from subluminal to superluminal.

Since the photon lifetime determined by the absorption is independent of the group velocity, the propagation length of the edge states is expected to increase as the group velocity increases. The propagation length of the edge state $l$ is defined as the length at which the intensity of the light decreases to $1/e$ :

\displaystyle l=\frac{\partial\omega_{R}}{\partial k_{x}}\frac{1}{2\omega_{I}},

(27)

Here, $\omega_{R}$ and $\omega_{I}$ are the real and imaginary parts of the eigenfrequency. Figure 5 (b) shows the propagation length of the edge states. Here, the case with balanced loss-gain perturbation (Fig. 4 (a)) is not shown since the propagation length is infinite due to the negative imaginary part of the eigenfrequency. The maximum propagation length is 20 $\mathrm{\mu}$ m when the permittivity adjustment $\Delta\varepsilon=-i0.2$ for $y>0$ is introduced. Moreover, the propagation length reaches 62 $\mathrm{\mu}$ m with the real part of permittivity adjustment $\Delta\varepsilon=-0.00235-i0.2$ for $y<0$ , which is due to the suppression of the EP smoothing. The enlarged propagation length despite short lifetime ( $\tau:=1/(2\omega_{I})=0.072$ ps) is achieved by the high group velocity effect near the EPs. Assuming the same group velocity in the unperturbed case and the same lifetime at EP for the perturbed cases (IV. in Fig. 5 (a)), the propagation length is limited to 4.0 $\mathrm{\mu}$ m. In this way, the propagation length of the edge states can be enlarged even in the lossy system, which may allow us to utilize the guided mode even in the passive system.

V implementation into the graphene-loaded photonic crystal slab

Our scheme for generating and restoring EPs in photonic crystal waveguides can be applied to the realistic silicon photonic crystal slabs. In this section, we show how to implement our scheme into the graphene-loaded photonic crystal slab.

Figure 6 (a) shows a schematic of the graphene loaded photonic crystal slab. We assume the lattice constant is $a=470$ nm, the one side of triangular hole is $s=0.8a$ , which is the same as the previous 2D simulations in Fig. 1 (a). The thickness of the slab is set to $h=200$ nm, and the permittivity of the slab is set to $\varepsilon_{\mathrm{slab}}={(3.48)}^{2}$ . Then, the 5-sheet of graphene is placed on the top of the slab and one side of the heterostructures ( $y<0$ ) to generate the EPs. Graphene has remarkable properties for applications in optical devices, such as large absorption in one place, reconfigurable conductivity [39]. It is well known that the optical absorption of graphene depends on its Fermi level, which can be controlled by voltage injection [40, 41, 42, 43, 44]. Our choice of graphene in this study is based solely on its ability to easily induce material absorption in photonic crystals, but the graphene has other advantages, such as the ability to reconfigurably control its permittivity. Now we aim to purely perturb the imaginary part of the permittivity since the large real part of the permittivity breaks the EP condition, and graphene is particularly well suited for perturbing only the imaginary part of the permittivity, as it has a large imaginary part of the permittivity in the near-infrared region [45]. In addition, some previous reports have demonstrated the successful loading and precise patterning of graphene onto photonic crystal structures [46, 47], further supporting its feasibility for our purposes. We have assumed that the conductivity of graphene is $\sigma=(6.0882+i0.1908)\times 10^{-5}$ [S/m] per sheet, which is the value from the theoretical formula derived from the Kubo formula [45] with the relaxation time $\tau=100$ fs, the temperature of graphene $T=300$ K, the angular frequency of light $\omega=2\pi c/\lambda$ , $\lambda=1550$ nm, and the Fermi energy of graphene $E_{F}=0.2$ eV. The real and imaginary parts of the conductivity correspond to the imaginary and real parts of the permittivity, respectively. See the supplementary material E for the effect of the misalignment of the graphene sheet.

Figure 6 (b) shows the band dispersion of the photonic crystal slab without graphene, and we can see that the linear dispersion of edge states appears near $k_{x}=\pi/a$ as in the case in Fig. 1 (b). Here, we focus on two eigenmodes at the Dirac point ( $k_{x}=\pi/a$ ), and the modes localized at the top and bottom are named modes u and d, as same as the case in Fig. 3 (b). The adjustment of the real part of the permittivity for one side of the heterostructure ( $y<0$ ) as shown in the right panel of Fig. 4 (a) can be replaced by the position shift of triangular holes. we choose the triangular hole position in the forth row from the glide symmetry plane to introduce just the right amount of adjustment, as shown by the yellow triangles in figure 6 (a). Figure 6 (c) shows the frequency of the two eigenmodes at $k_{x}=\pi/a$ without graphene as a function of the $y$ -directional triangular hole shift $\Delta y_{\mathrm{hole}}$ . The frequency of two modes are degenerated in $\Delta y_{\mathrm{hole}}=0$ . The positive hole shift $\Delta y_{\mathrm{hole}}$ perturbs only the mode $\mathrm{d}$ , leads to decrease the difference of frequency between two modes and minimizes near $\Delta y_{\mathrm{hole}}=104$ nm. Therefore, the hole shift of the triangular hole can be used to adjust the effective real part of the permittivity. A details is shown in the supplementary material D.

Next, we investigate the effect of the graphene loading on the band dispersion. Figure 6 (d) shows the band dispersion near $k_{x}=\pi/a$ of the graphene-loaded photonic crystal slab without the triangular hole shift ( $\Delta y_{\mathrm{hole}}=0$ ). The EPs is smoothed out by the graphene loading, which is not only due to the imaginary coupling in loss-biased system, but also due to the positive permittivity modulation by the graphene sheet. Figure 6 (e) shows the band dispersion with the triangular hole shift $\Delta y_{\mathrm{hole}}=104$ nm. The proper adjustment of triangular holes cancels the contribution of both the imaginary coupling and the imaginary part of conductivity, and restores EP. Consequently, our scheme for generating and restoring EPs in photonic crystal waveguides can be applied to the silicon photonic crystal slabs by loading graphene and adjusting the position of triangular holes.

Finally, we investigate the group velocity of the edge state and corresponding propagation length in some triangular hole shift $\Delta y_{\mathrm{hole}}$ . Figure 7 (a) shows the group velocity $v_{g}$ of the edge states in the graphene-loaded photonic crystal slab near EPs with some triangular hole shift $\Delta y_{\mathrm{hole}}$ . The group velocity of initial photonic crystal waveguide near the Dirac cone is about $v_{g}=0.13c$ , and increases to $v_{g}=0.38c$ when graphene is loaded. However, without the hole shift, the effect of the EP smoothing strongly limits the peak of group velocity. When a hole shift is added, the EP smoothing gradually mitigates and the group velocity peak is maximized around $\Delta y_{\mathrm{hole}}=104$ nm. A series of changes in the group velocity with $\Delta y_{\mathrm{hole}}$ correspond to the change of frequency difference at $k_{x}a=\pi$ between the two modes provided by the hole shift, as shown in Fig. 6(d). The group velocity reaches $v_{g}=3.3c$ near the restored EPs in the case of $\Delta y_{\mathrm{hole}}=104$ nm, and is accelerated to 25 times of that without graphene. The bandwidth where the group velocity is accelerated over the speed of light in vacuum is up to 0.4 nm as shown in Fig. 7 (a), which corresponds to $1.5$ ps in terms of time width assuming a Gaussian pulse. Importantly, the propagation length is also elongated by the same ratio, directly resulting from the fast light effect. We confirmed that the propagation length $l$ is elongated to 105 $\mathrm{\mu}$ m near the restored EPs as shown in Fig. 7 (b), which is a staggering number for a lifetime of $0.11$ ps. It takes 0.38 ps to propagate the same distance in air, and it is expected that this difference can be detected for a pulse length of 1.5 ps, given the ratio of pulse length to group delay $(0.38-0.11)/1.5\approx 0.18$ : This value is better than the previously reported ratio of pulse width to group delay in superluminal propagation experiments $(23.8\mathrm{\mu s})/(0.05\mathrm{s})\approx 0.048$ [22]. In addition, the group velocity of the edge state in a photonic crystal waveguide without graphene is about 1/7 of that in air, making it even easier to measure the difference in group delay.

VI Conclusion

We have proposed a feasible photonic crystal waveguide for realizing superluminal edge states. A systematic method to produce EPs in non-Hermitian photonic crystal waveguides with glide and time reversal symmetry has been presented, and the exact EP condition about the permittivity adjustment is analytically derived by $\bm{k}\cdot\bm{p}$ Hamiltonian. We have shown that the EPs can be generated by introducing the loss-bias to the system, and the EPs can be restored by appropriately adjusting the permittivity. The anomalous dispersion near the EPs leads to the high group velocity of the edge states, which allows the realization of superluminal propagation. In addition, the high group velocity of the edge states near the EPs enlarges the propagation length of the waveguide, even in the lossy and short lifetime system. We have confirmed that the group velocity contrast of the edge state is increased and corresponding propagation length is enlarged by appropriately adjusting the permittivity to restore EP. Our scheme for generating and restoring EPs in photonic crystal waveguides can be applied to the silicon photonic crystal slabs by loading graphene and adjusting the position of triangular holes. We have confirmed that the group velocity $v_{g}$ is accelerated to $25$ times of that of w/o graphene near the restored EPs, thereby realizing the superluminal edge state ( $v_{g}=3.3c$ ). Importantly, the propagation length is also elongated by the same ratio, directly resulting from the fast light effect.

Our study pave the way to realize the observation of superluminal propagation states in photonic crystal waveguide. The feasibility and integration capabilities, combined with graphene’s reversibility, make it a promising candidate for optical devices. In addition, the fast-light on-tip nanodevices have the potential for the development of control of light-matter interaction as well as the traditional slow-light waveguide, opening up new possibilities for the application of optical computing. Furthermore, in addition to the anomalous dispersion in the edge states, our structure has additional intriguing aspect about the topological properties. The proposed structure is based on a photonic crystal waveguide, known for exhibiting edge states with valley-polarized characteristics and robustness against bending [48, 19]. Our structure involves uniform permittivity perturbations on both sides of a glide plane, which preserves the symmetry of the effective Hamiltonian for the bulk. This results in maintaining a finite Berry curvature near the K and K’ points in the band structure. Consequently, the edge states near the EPs retain valley-polarized topological features. Although a detailed discussion of valley topology is beyond the scope of this work, our proposed structure holds potential for application in valley photonic devices, and provides a pathway for developing innovative photonic devices.

Acknowledgements.

This work was supported by the Japan Society for the Promotion of Science (Grant number JP20H05641, JP21K14551, 24K01377, 24H02232 and 24H00400).

Appendix A The design of heterostructure in other lattice or hole shapes

Our scheme for generating EPs from a symmetry oriented Dirac point can be applied to other lattice structures and hole shapes as long as the spatial photonic crystal geometry has the glide symmetry. Figure 8 (a) shows the schematic and the band dispersion in the honeycomb lattice photonic crystal. The slab permittivity and the amount of the perturbation are the same as in the triangular lattice. It is obvious that the EPs generated by the non-Hermitian perturbation on both sides of the glide plane are the same as in the triangular lattice. On the other hand, the generated EPs are outside the band gap when we adopt the honeycomb lattice, and it is difficult to bring the EPs inside the band gap by some kind of tuning. Figure 8 (b) shows the results when we adopt the triangular lattice with circular hole, and there is no significant difference compared with the use of triangular holes. Thus, our proposed method for generating symmetry-oriented EPs can be applied to photonic crystal with different lattice structures. The structure of the triangular lattice and triangular holes shown in Fig. 1 will be analyzed in the following section because it is suitable for fast-light waveguide applications since the EP is inside the band gap, the band slope (group velocity) near the Dirac point is large, and the band gap is wide.

Appendix B $H_{z}$ distribution of the eigenmodes

In the context of a two-band system with PT symmetry, the phase transition of eigenvalues and eigenstates through the EP is known as the PT phase transition [5]. Since we now focus on the TE mode, the two eigenmodes are characterized by their $H_{z}$ distributions. For three points, namely outside the EP ( $k_{x}=0.98\times 2\pi/a$ ), at the EP ( $k_{x}=0.49297\times 2\pi/a$ ), and at the original Dirac point ( $k_{x}=\pi/a$ ), we present the $H_{z}$ distributions of the two bands at each $k_{x}$ point in Fig. 9. We discuss the presence or absence of PT symmetry for each case. At $k_{x}=0.49\times 2\pi/a$ , the $H_{z}$ distribution is symmetric about $y=0$ , which is consistent with both eigenvalues being real. In addition, the symmetry under simultaneous space and time reversal operations, i.e., $(\hat{G}\hat{T})H_{z}^{(i)}=H_{z}^{(i)}$ and $i=1,2$ denotes the label of eigenmodes, is evident from the figure (PT-symmetric phase). Similar features are observed in the EP ( $k_{x}=0.49297\times 2\pi/a$ ), but the two eigenmodes show complete coalescence, which distinguishes it from the $k_{x}=0.49\times 2\pi/a$ case. At $k_{x}=\pi/a$ , the $H_{z}$ distribution under PT operations is no longer symmetric, but shows an asymmetric distribution with one side favoring gain and the other side favoring loss (PT-broken phase). In summary, the $H_{z}$ distributions of the eigenmodes in the heterostructure shown in Fig. 1 satisfy the properties of a general PT-symmetric system.

Appendix C Derivation of the $\bm{k}\cdot\bm{p}$ Hamiltonian

Here, we derive the $\bm{k}\cdot\bm{p}$ Hamiltonian from the Maxwell’s equations. The details of the derivation are shown in Ref. [38].

The Maxwell’s equations in the photonic crystal are given by [49]

\displaystyle\bm{\nabla}\times\bm{\nabla}\times\bm{E}_{nk}(\bm{r})=\left[\varepsilon(\bm{r})+\Delta\varepsilon(\bm{r})\right]{\left(\frac{\omega_{nk}}{c}\right)}^{2}\bm{E}_{nk}(\bm{r})

(28)

Here, $\varepsilon(\mathbf{r})$ is the permittivity of initial photonic crystal, $\Delta\varepsilon(\mathbf{r}):=\Delta\varepsilon_{\mathbf{r}}+i\Delta\varepsilon_{\mathbf{i}}(\mathbf{r})$ is the perturbation of the permittivity. $\omega_{nk}$ and $\bm{E}_{nk}(\bm{r})$ are the eigenfrequency and electric filed profile of the photonic crystal, respectively. And, $\omega_{n}^{(0)}$ is the unperturbed eigenfrequency. Let $U_{n}^{(0)}$ be the unperturbed eigenenergy, and the bases $\bm{E}_{nk}^{(0)}(\bm{r})$ holds the orthonormality condition:

\displaystyle\frac{\varepsilon_{0}}{4}\int_{\mathrm{u.c.}}\varepsilon(\bm{r})\bm{E}_{nk}^{(0)*}(\bm{r})\cdot\bm{E}_{mk}^{(0)}(\bm{r})\mathrm{d}\bm{r}=U_{n}^{(0)}\delta_{nm}

(29)

Then, we define the normalized eigenmodes $\bm{\psi}_{nk}^{(0)}(\bm{r})$ as

\displaystyle\bm{\psi}_{nk}^{(0)}(\bm{r})=\frac{\sqrt{\varepsilon_{0}}}{2\sqrt{U_{n}^{(0)}}}\bm{E}_{nk}^{(0)}(\bm{r})

(30)

and we get the following equation:

	$\displaystyle\bm{\nabla}\times\bm{\nabla}\times\bm{\psi}_{nk}(\bm{r})=\left[\varepsilon(\bm{r})+\Delta\varepsilon(\bm{r})\right]{\left(\frac{\omega_{nk}}{c}\right)}^{2}\bm{\psi}_{nk}(\bm{r})$		(31)
	$\displaystyle\int_{\mathrm{u.c.}}\varepsilon(\bm{r})\bm{\psi}_{nk}^{(0)*}(\bm{r})\cdot\bm{\psi}_{mk}^{(0)}(\bm{r})\mathrm{d}\bm{r}=\delta_{nm}$		(32)

Like the electric field, the magnetic field is also normalised as

\displaystyle\bm{\varphi}_{nk}^{(0)}(\bm{r})=\frac{\sqrt{\mu_{0}}}{2\sqrt{U_{n}^{(0)}}}\bm{H}_{nk}^{(0)}(\bm{r})

(33)

Next, we expand the wave function $\bm{\psi}_{nk}(\bm{r})$ when $\Delta\varepsilon(\bm{r})$ is finite by the wave function at $\bm{k}={[k_{0},0,0]}^{\top},\ k_{0}=\pi/a$ with $\Delta\varepsilon(\bm{r})=0$ :

\displaystyle\bm{\psi}_{nk}(\bm{r})

\displaystyle=\sum_{j}C_{nj}(k)e^{i\bm{\delta k}\cdot\bm{r}}\bm{\psi}_{jk_{0}}^{(0)}(\bm{r})

(34)

Here, $\bm{\delta k}=\bm{k}-{[k_{0},0,0]}^{\top}$ and $C_{nj}(k)$ is the expansion coefficient. Substitute Eq. (34) into Eq. (31) and we get

	$\displaystyle\sum_{j}C_{nj}(k)e^{i\bm{\delta k}\cdot\bm{r}}A\bm{\psi}_{jk_{0}}^{(0)}(\bm{r})$
	$\displaystyle={\left(\frac{\omega_{nk}}{c}\right)}^{2}\left[\varepsilon(\bm{r})+\Delta\varepsilon(\bm{r})\right]\sum_{j}C_{nj}(k)e^{i\bm{\delta k}\cdot\bm{r}}\bm{\psi}_{jk_{0}}^{(0)}(\bm{r})$
$\displaystyle A$	$\displaystyle:=-\bm{\delta k}\times\bm{\delta k}\times+i\left(\bm{\nabla}\times\bm{\delta k}\times\right.$
	$\displaystyle\left.\qquad\qquad\qquad+\bm{\delta k}\times\bm{\nabla}\times\right)+\bm{\nabla}\times\bm{\nabla}\times$	(35)

Furthermore, by multiplying from the left by $e^{-i\bm{\delta k}\cdot\bm{r}}\bm{\psi}_{ik_{0}}^{(0)*}(\bm{r})$ and integrating over the unit cell, the eigenvalue equation determining the band structure is obtained.

	$\displaystyle\sum_{j}\left[{\left(\frac{\omega_{jk}^{(0)}}{c}\right)}^{2}\delta_{ij}+P_{ij}^{\prime}+Q_{ij}^{\prime}\right]C_{nj}(k)$
	$\displaystyle\qquad\qquad\qquad\qquad={\left(\frac{\omega_{nk}}{c}\right)}^{2}\sum_{j}\left(\delta_{ij}+G_{ij}\right)$		(36)

Here, $\displaystyle\bm{\nabla}\times\bm{\nabla}\times\bm{\psi}_{jk}^{(0)}(\bm{r})=\varepsilon(\bm{r}){\left(\frac{\omega_{nk}^{(0)}}{c}\right)}^{2}\bm{\psi}_{jk}^{(0)}(\bm{r})$ is used in the derivation. $P_{ij}^{\prime}$ , $G_{ij}$ , $Q^{\prime}_{ij}$ are defined as follows:

$\displaystyle P_{ij}^{\prime}$	$\displaystyle=\int_{\mathrm{u.c.}}\left(\omega_{j}^{(0)}\bm{\psi}_{ik_{0}}^{(0)*}(\bm{r})\times\bm{\varphi}_{jk_{0}}^{(0)}(\bm{r})\right.$
	$\displaystyle\qquad\qquad\qquad\left.-\omega_{i}^{(0)}\bm{\varphi}_{ik_{0}}^{(0)*}(\bm{r})\times\bm{\psi}_{jk_{0}}^{(0)}(\bm{r})\right)\cdot\bm{\delta k}\,\mathrm{d}^{3}r$
	$\displaystyle=:\bm{P}_{ij}\cdot\bm{\delta k}$	(37)
$\displaystyle G_{ij}$	$\displaystyle=\int_{\mathrm{u.c.}}\Delta\varepsilon(\bm{r})\bm{\psi}_{ik_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{jk_{0}}^{(0)}(\bm{r})\mathrm{d}^{3}r$	(38)
$\displaystyle Q^{\prime}_{ij}$	$\displaystyle=-\int_{\mathrm{u.c.}}\bm{\psi}_{ik_{0}}^{(0)*}(\bm{r})\cdot\left(\bm{\delta k}\times\bm{\delta k}\times\bm{\psi}_{jk_{0}}^{(0)}(\bm{r})\right)\mathrm{d}^{3}r$
	$\displaystyle={(\delta k_{x})}^{2}\int_{\mathrm{u.c.}}\left({\left(\psi_{ik_{0}}^{(0)*}(\bm{r})\right)}_{y}{\left(\psi_{jk_{0}}^{(0)}(\bm{r})\right)}_{y}\right.$
	$\displaystyle\qquad\qquad\qquad\qquad\left.+{\left(\psi_{ik_{0}}^{(0)*}(\bm{r})\right)}_{z}{\left(\psi_{jk_{0}}^{(0)}(\bm{r})\right)}_{z}\right)\mathrm{d}^{3}r$
	$\displaystyle=:{(\delta k_{x})}^{2}Q_{ij}$	(39)

Appendix D The EP condition with the hole shift in $\bm{k}\cdot\bm{p}$ Hamiltonian

We have shown in Section III.2 of that the EP can be restored by appropriately adjusting the real part of the permittivity, and adopted the shift of the triangular holes in the realistic structure as a method instead of the permittivity adjustment in Section 11. Here, we theoretically investigate the effect of the triangular hole shift and show that the triangular hole shift offset the effect of the permittivity adjustment.

The following describes the setup for treating the triangular hole shift and graphene loading as perturbations in the permittivity. A domain $V_{1},V_{2}$ and surface $S$ are defined as shown in Fig. 10(a). The domain $V_{1}$ and $V_{2}$ are defined as the regions where the permittivity changes from silicon to air and from air to silicon, respectively, when the triangular holes are shifted by $\Delta y_{\mathrm{hole}}$ . Moreover, the surface $S$ is defined as the region where the graphene is loaded. The elements of the matrix $G$ defined by the Eq.(6) are the sum of the contributions from the domains $V_{1},V_{2}$ and the surface $S$ as follows:

\displaystyle G_{\mathrm{\alpha\beta}}

\displaystyle=G_{\mathrm{\alpha\beta}}^{(V_{1})}+G_{\mathrm{\alpha\beta}}^{(V_{2})}+G_{\mathrm{\alpha\beta}}^{(S)},\quad\alpha,\beta=\mathrm{d,u}

(40)

Here, $G_{\mathrm{ud}},G_{\mathrm{du}}$ are zero because all the domains and surfaces have mirror symmetry with respect to the $y$ -axis of the unit cell. The contributions from the domains $V_{1},V_{2}$ are given by

	$\displaystyle G_{\mathrm{\alpha\beta}}^{(V_{1})}$	$\displaystyle=(\varepsilon_{\mathrm{air}}-\varepsilon_{\mathrm{slab}})\int_{V_{1}}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\mathrm{d}^{3}r$		(41)
	$\displaystyle G_{\mathrm{\alpha\beta}}^{(V_{2})}$	$\displaystyle=(\varepsilon_{\mathrm{slab}}-\varepsilon_{\mathrm{air}})\int_{V_{2}}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\mathrm{d}^{3}r$		(42)

The permittivity of the silicon and air are denoted as $\varepsilon_{\mathrm{slab}}$ and $\varepsilon_{\mathrm{air}}$ , respectively. To simplify the notation, we introduce the following quantities:

\displaystyle I_{\mathrm{\alpha\beta}}^{(V)}

\displaystyle=\int_{V}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\mathrm{d}^{3}r

(43)

Then, the contributions from $V_{1}$ and $V_{2}$ are collectively denoted as $Y_{\mathrm{d}},Y_{\mathrm{u}}$ :

	$\displaystyle Y_{\mathrm{d}}:=G_{\mathrm{dd}}^{(V_{1})}+G_{\mathrm{dd}}^{(V_{2})}$	$\displaystyle=(\varepsilon_{\mathrm{slab}}-\varepsilon_{\mathrm{air}})(I_{\mathrm{dd}}^{(V_{2})}-I_{\mathrm{dd}}^{(V_{1})})$		(44)
	$\displaystyle Y_{\mathrm{u}}:=G_{\mathrm{uu}}^{(V_{1})}+G_{\mathrm{uu}}^{(V_{2})}$	$\displaystyle=(\varepsilon_{\mathrm{slab}}-\varepsilon_{\mathrm{air}})(I_{\mathrm{uu}}^{(V_{2})}-I_{\mathrm{uu}}^{(V_{1})})$		(45)

Next, we consider the contribution from the surface $S$ . Since the complex conductivity of the graphene is given by $\sigma=\sigma_{r}+i\sigma_{i}$ , the permittivity perturbation due to the graphene loading is given by

\displaystyle\Delta\varepsilon(\bm{r})=-i\frac{\sigma}{t\omega\varepsilon_{0}},\quad\bm{r}\in S

(46)

Here, $t$ is the thickness of the graphene, $\omega$ is the angular frequency, and $\varepsilon_{0}$ is the vacuum permittivity. For simplicity, we will ignore wavelength dispersion of the frequency. Then, the contribution from the surface $S$ is given by

$\displaystyle G_{\mathrm{\alpha\beta}}^{(S)}$	$\displaystyle=\int_{S}\frac{-i\sigma}{t\omega\varepsilon_{0}}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})t\,\mathrm{d}^{2}r$
	$\displaystyle=-\frac{i\sigma}{\omega\varepsilon_{0}}\int_{S}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\,\mathrm{d}^{2}r$
	$\displaystyle=:-i\sigma L_{\alpha\beta}$	(47)

where

\displaystyle L_{\alpha\beta}=\frac{1}{\omega\varepsilon_{0}}\int_{S}\bm{\psi}_{\alpha k_{0}}^{(0)*}(\bm{r})\cdot\bm{\psi}_{\beta k_{0}}^{(0)}(\bm{r})\,\mathrm{d}^{2}r

(48)

The contributions from the domains $V_{1},V_{2}$ and the surface $S$ are collectively denoted as

	$\displaystyle G_{\mathrm{dd}}$	$\displaystyle=Y_{\mathrm{d}}+(\sigma_{i}-i\sigma_{r})L_{\mathrm{dd}},$		(49)
	$\displaystyle G_{\mathrm{uu}}$	$\displaystyle=Y_{\mathrm{u}}+(\sigma_{i}-i\sigma_{r})L_{\mathrm{uu}}$		(50)

All $V_{1}$ , $V_{2}$ and $S$ depend on the hole shift $\Delta y_{\mathrm{hole}}$ , but the value of $L_{\alpha\beta}$ is almost constant since the contribution of the hole shift in the graphene-loaded region is small.

Next, we calculate the components to evaluate the EP condition in Eq. (24).

$\displaystyle 4\Delta G$	$\displaystyle=G_{\mathrm{dd}}-G_{\mathrm{uu}}$
	$\displaystyle=Y_{\mathrm{d}}-Y_{\mathrm{u}}+(\sigma_{i}-i\sigma_{r})(L_{\mathrm{dd}}-L_{\mathrm{uu}})$
	$\displaystyle=:\Delta y+(\sigma_{i}-i\sigma_{r})\Delta L$	(51)

	$\displaystyle(1+G_{\mathrm{dd}})(1+G_{\mathrm{uu}})$
	$\displaystyle=(1+Y_{\mathrm{d}}+\sigma_{i}L_{\mathrm{dd}}-i\sigma_{r}L_{\mathrm{dd}})(1+Y_{\mathrm{u}}+\sigma_{i}L_{\mathrm{uu}}-i\sigma_{r}L_{\mathrm{uu}})$
	$\displaystyle=(1+Y_{\mathrm{d}}+\sigma_{i}L_{\mathrm{dd}})(1+Y_{\mathrm{u}}+\sigma_{i}L_{\mathrm{uu}})-\sigma_{r}^{2}L_{\mathrm{dd}}L_{\mathrm{uu}}$
	$\displaystyle\quad-i\sigma_{r}\left((1+Y_{\mathrm{d}}+\sigma_{i}L_{\mathrm{dd}})L_{\mathrm{uu}}+(1+Y_{\mathrm{u}}+\sigma_{i}L_{\mathrm{uu}})L_{\mathrm{dd}}\right)$

By substituting the above equations into Eq. (24), we obtain the EP condition as follows:

\displaystyle\sigma_{r}^{2}

\displaystyle=\frac{m^{2}l-2m\Delta L(1+Y_{\mathrm{d}}+\sigma_{i}L_{\mathrm{dd}})(1+Y_{\mathrm{u}}+\sigma_{i}L_{\mathrm{uu}})}{l{(\Delta L)}^{2}-2m\Delta LL_{\mathrm{dd}}L_{\mathrm{uu}}}

(53)

where

	$\displaystyle l$	$\displaystyle=(1+Y_{\mathrm{d}}+\sigma_{i}L_{\mathrm{dd}})L_{\mathrm{uu}}+(1+Y_{\mathrm{u}}+\sigma_{i}L_{\mathrm{uu}})L_{\mathrm{dd}}$		(54)
	$\displaystyle m$	$\displaystyle=Y_{\mathrm{d}}-Y_{\mathrm{u}}+\sigma_{i}\Delta L$		(55)

Eq. (53) assures that the EP can be restored by shifting the triangular holes. From the expression for $Y_{\mathrm{d}}+\sigma_{i}L_{\mathrm{dd}}$ and $Y_{\mathrm{u}}+\sigma_{i}L_{\mathrm{uu}}$ that appears in Eq.(53), it can be seen that the triangular hole shift $\Delta y_{\mathrm{hole}}$ and the imaginary part of the graphene conductivity (i.e., the real part of the permittivity) have the opposite effect on the permittivity perturbation. Generally, $I_{\mathrm{dd}}^{(V_{2})}<I_{\mathrm{dd}}^{(V_{1})}$ , $I_{\mathrm{uu}}^{(V_{2})}<I_{\mathrm{uu}}^{(V_{1})}$ hold since the edge states are localized near the glide plane. Consequently, $Y_{\mathrm{d}},Y_{\mathrm{u}}$ defined by Eqs. (44), (45) are negative. Therefore, the triangular hole shift $\Delta y_{\mathrm{hole}}$ has the same effect as a net negative tuning of the real part of the permittivity. Furthermore, due to the difference in the electric field distribution between modes $\mathrm{u}$ and $\mathrm{d}$ , the permittivity perturbation due to the hole shift affects only the $\mathrm{d}$ mode. Note that since graphene has a positive imaginary part of its electrical conductivity, i.e., a positive real part of its permittivity, it is not possible to restore the EP only using graphene.

Figure 10 (b) shows the relationship between the real part of the conductivity and the triangular hole shift required to satisfy the EP condition when the imaginary part of the conductivity is fixed. When the imaginary part of the conductivity is zero, the triangular hole shift that satisfies the EP condition is approximately proportional to the square of the real part of the conductivity. As previously mentioned, a positive imaginary part of the conductivity has an opposite effect to the triangular hole shift, so the larger the imaginary part of the conductivity, the narrower the range of the real part of the conductivity that can be canceled out by the triangular hole shift. In practice, the imaginary part of the conductivity for graphene is about $\sigma_{i}=2\times 10^{-6}$ [S/m] [45] so within the range of perturbation theory, it is difficult to restore EP using only the hole shift. In actual FEM calculations, however, the EP restoration has been confirmed with a hole shift of about $\Delta y_{\mathrm{hole}}=105$ nm, but this is believed to be due to the large permittivity contrast $\varepsilon_{\mathrm{slab}}-\varepsilon_{\mathrm{air}}\approx 11$ , which causes changes in the electric field distribution, making the contribution of the actual triangular hole shift larger than the linear perturbation effect. The change in the electric field distribution due to the triangular hole shift is confirmed in Fig. 10 (c).

Appendix E Restoration of EPs and group velocity in graphene-Loaded structures with misalignment

The graphene-loaded structure proposed in Section assumes that graphene is placed along the glide-symmetry plane. From the perspective of fabrication and measurement, it is important to consider the effects on EP and group velocity when the graphene alignment is offset. In practice, even if the graphene position is misaligned, our scheme for EP restoration can still be applied without issue. Figure 11 (a) shows the case where the graphene loading position is shifted by $\sqrt{3}a/4$ below the glide plane. For this structure, the band diagram and group velocity for both the case without triangular hole shifts and with a shift of 44 nm are shown in Fig. 11 (b) and (c), respectively. As is evident from the figures, even when the graphene position is misaligned, the EP is restored, and the singularity in group velocity is recovered as shown in Fig. 11 (d). Mathematically, the misalignment of the graphene loading position is equivalent to changing the magnitude of graphene’s complex permittivity. In other words, it ultimately affects the value of the integral $L_{\alpha\beta}$ in Eq. 48. Therefore, by adjusting the value of the triangular hole shift or permittivity shift required to cancel the EP, it is possible to restore the EP in exactly the same manner as in the case without misalignment.

References

Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real spectra in non-hermitian hamiltonians having $\mathcal{P}\mathcal{T}$ symmetry, Phys. Rev. Lett. 80, 5243 (1998).
Bender [2007] C. M. Bender, Making sense of non-hermitian hamiltonians, Reports on Progress in Physics 70, 947 (2007).
Parto et al. [2018] M. Parto, S. Wittek, H. Hodaei, G. Harari, M. A. Bandres, J. Ren, M. C. Rechtsman, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Edge-mode lasing in 1d topological active arrays, Phys. Rev. Lett. 120, 113901 (2018).
Takata et al. [2021] K. Takata, K. Nozaki, E. Kuramochi, S. Matsuo, K. Takeda, T. Fujii, S. Kita, A. Shinya, and M. Notomi, Observing exceptional point degeneracy of radiation with electrically pumped photonic crystal coupled-nanocavity lasers, Optica 8, 184 (2021).
Feng et al. [2017] L. Feng, R. El-Ganainy, and L. Ge, Non-hermitian photonics based on parity–time symmetry, Nature Photonics 11, 752 (2017).
Regensburger et al. [2012] A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel, Parity-time synthetic photonic lattices, Nature 488, 167 (2012).
Özdemir et al. [2019] Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Parity–time symmetry and exceptional points in photonics, Nature Materials 18, 783 (2019).
Peng et al. [2014] B. Peng, Ş. K. Özdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, Parity–time-symmetric whispering-gallery microcavities, Nature Physics 10, 394 (2014).
Feng et al. [2014] L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, and X. Zhang, Single-mode laser by parity-time symmetry breaking, Science 346, 972 (2014).
Hodaei et al. [2014] H. Hodaei, M.-A. Miri, M. Heinrich, D. N. Christodoulides, and M. Khajavikhan, Parity-time-symmetric microring lasers, Science 346, 975 (2014).
Guo et al. [2009] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Observation of $\mathcal{P}\mathcal{T}$ -symmetry breaking in complex optical potentials, Phys. Rev. Lett. 103, 093902 (2009).
Szameit et al. [2011] A. Szameit, M. C. Rechtsman, O. Bahat-Treidel, and M. Segev, $\mathcal{P}\mathcal{T}$ -symmetry in honeycomb photonic lattices, Phys. Rev. A 84, 021806 (2011).
Baba [2008] T. Baba, Slow light in photonic crystals, Nature Photonics 2, 465 (2008).
Yariv et al. [1999] A. Yariv, Y. Xu, R. K. Lee, and A. Scherer, Coupled-resonator optical waveguide:?a proposal and analysis, Opt. Lett. 24, 711 (1999).
Olivier et al. [2001] S. Olivier, C. Smith, M. Rattier, H. Benisty, C. Weisbuch, T. Krauss, R. Houdré, and U. Oesterlé, Miniband transmission in a photonic crystal coupled-resonator optical waveguide, Opt. Lett. 26, 1019 (2001).
Notomi et al. [2008] M. Notomi, E. Kuramochi, and T. Tanabe, Large-scale arrays of ultrahigh-q coupled nanocavities, Nature Photonics 2, 741 (2008).
Vlasov et al. [2005] Y. A. Vlasov, M. O’Boyle, H. F. Hamann, and S. J. McNab, Active control of slow light on a chip with photonic crystal waveguides, Nature 438, 65 (2005).
Yoshimi et al. [2020] H. Yoshimi, T. Yamaguchi, Y. Ota, Y. Arakawa, and S. Iwamoto, Slow light waveguides in topological valley photonic crystals, Opt. Lett. 45, 2648 (2020).
Yoshimi et al. [2021] H. Yoshimi, T. Yamaguchi, R. Katsumi, Y. Ota, Y. Arakawa, and S. Iwamoto, Experimental demonstration of topological slow light waveguides in valley photonic crystals, Opt. Express 29, 13441 (2021).
Brillouin [1960] L. Brillouin, Wave propagation and group velocity, Vol. 8 (Academic Press, New York, London, 1960).
Wang et al. [2000] L. J. Wang, A. Kuzmich, and A. Dogariu, Gain-assisted superluminal light propagation, Nature 406, 277 (2000).
Bigelow et al. [2003] M. S. Bigelow, N. N. Lepeshkin, and R. W. Boyd, Superluminal and slow light propagation in a room-temperature solid, Science 301, 200 (2003).
Ziolkowski [2001] R. W. Ziolkowski, Superluminal transmission of information through an electromagnetic metamaterial, Phys. Rev. E 63, 046604 (2001).
Xu et al. [2012] S. Xu, X. Cheng, S. Xi, R. Zhang, H. O. Moser, Z. Shen, Y. Xu, Z. Huang, X. Zhang, F. Yu, B. Zhang, and H. Chen, Experimental demonstration of a free-space cylindrical cloak without superluminal propagation, Phys. Rev. Lett. 109, 223903 (2012).
Longhi [2013] S. Longhi, Convective and absolute $\mathcal{PT}$ -symmetry breaking in tight-binding lattices, Phys. Rev. A 88, 052102 (2013).
Schomerus and Wiersig [2014] H. Schomerus and J. Wiersig, Non-hermitian-transport effects in coupled-resonator optical waveguides, Phys. Rev. A 90, 053819 (2014).
Takata and Notomi [2017] K. Takata and M. Notomi, $\mathcal{P}\mathcal{T}$ -symmetric coupled-resonator waveguide based on buried heterostructure nanocavities, Phys. Rev. Appl. 7, 054023 (2017).
Takata et al. [2022] K. Takata, N. Roberts, A. Shinya, and M. Notomi, Imaginary couplings in non-hermitian coupled-mode theory: Effects on exceptional points of optical resonators, Phys. Rev. A 105, 013523 (2022).
Cui et al. [2019] X. Cui, K. Ding, J.-W. Dong, and C. T. Chan, Exceptional points and their coalescence of $\mathcal{PT}$ -symmetric interface states in photonic crystals, Phys. Rev. B 100, 115412 (2019).
Mock [2020] A. Mock, Symmetry-engineered waveguide dispersion in $\mathcal{P}\mathcal{T}$ symmetric photonic crystal waveguides, J. Opt. Soc. Am. B 37, 168 (2020).
Fang et al. [2019] Y.-T. Fang, S.-F. Ye, and X.-X. Li, Unique band coalescence and exceptional points from two-dimensional photonic crystal waveguide, Journal of Optics 21, 055103 (2019).
Iglesias Martínez et al. [2022] J. A. Iglesias Martínez, N. Laforge, M. Kadic, and V. Laude, Topological waves guided by a glide-reflection symmetric crystal interface, Phys. Rev. B 106, 064304 (2022).
Plotnik et al. [2014] Y. Plotnik, M. C. Rechtsman, D. Song, M. Heinrich, J. M. Zeuner, S. Nolte, Y. Lumer, N. Malkova, J. Xu, A. Szameit, Z. Chen, and M. Segev, Observation of unconventional edge states in ‘photonic graphene’, Nature Materials 13, 57 (2014).
Yoda and Notomi [2019] T. Yoda and M. Notomi, Air-hole-type valley photonic crystal slab with simple triangular lattice for valley-contrasting physics, in CLEO: Science and Innovations (Optica Publishing Group, 2019) pp. JTh2A–10.
Söllner et al. [2015] I. Söllner, S. Mahmoodian, S. L. Hansen, L. Midolo, A. Javadi, G. Kiršanskė, T. Pregnolato, H. El-Ella, E. H. Lee, J. D. Song, S. Stobbe, and P. Lodahl, Deterministic photon–emitter coupling in chiral photonic circuits, Nature Nanotechnology 10, 775 (2015).
Mahmoodian et al. [2017] S. Mahmoodian, K. Prindal-Nielsen, I. Söllner, S. Stobbe, and P. Lodahl, Engineering chiral light–matter interaction in photonic crystal waveguides with slow light, Opt. Mater. Express 7, 43 (2017).
Liu et al. [2016] Z.-P. Liu, J. Zhang, i. m. c. K. Özdemir, B. Peng, H. Jing, X.-Y. Lü, C.-W. Li, L. Yang, F. Nori, and Y.-x. Liu, Metrology with $\mathcal{PT}$ -symmetric cavities: Enhanced sensitivity near the $\mathcal{PT}$ -phase transition, Phys. Rev. Lett. 117, 110802 (2016).
Cerjan et al. [2016] A. Cerjan, A. Raman, and S. Fan, Exceptional contours and band structure design in parity-time symmetric photonic crystals, Phys. Rev. Lett. 116, 203902 (2016).
Castro Neto et al. [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
Tan et al. [2007] Y.-W. Tan, Y. Zhang, K. Bolotin, Y. Zhao, S. Adam, E. H. Hwang, S. Das Sarma, H. L. Stormer, and P. Kim, Measurement of scattering rate and minimum conductivity in graphene, Phys. Rev. Lett. 99, 246803 (2007).
Tse et al. [2008] W.-K. Tse, E. H. Hwang, and S. Das Sarma, Ballistic hot electron transport in graphene, Applied Physics Letters 93, 023128 (2008).
Breusing et al. [2011] M. Breusing, S. Kuehn, T. Winzer, E. Malić, F. Milde, N. Severin, J. P. Rabe, C. Ropers, A. Knorr, and T. Elsaesser, Ultrafast nonequilibrium carrier dynamics in a single graphene layer, Phys. Rev. B 83, 153410 (2011).
Wang et al. [2010] H. Wang, J. H. Strait, P. A. George, S. Shivaraman, V. B. Shields, M. Chandrashekhar, J. Hwang, F. Rana, M. G. Spencer, C. S. Ruiz-Vargas, and J. Park, Ultrafast relaxation dynamics of hot optical phonons in graphene, Applied Physics Letters 96, 081917 (2010).
Mihnev et al. [2016] M. T. Mihnev, F. Kadi, C. J. Divin, T. Winzer, S. Lee, C.-H. Liu, Z. Zhong, C. Berger, W. A. de Heer, E. Malic, A. Knorr, and T. B. Norris, Microscopic origins of the terahertz carrier relaxation and cooling dynamics in graphene, Nature Communications 7, 11617 (2016).
Hanson [2008] G. W. Hanson, Dyadic Green’s functions and guided surface waves for a surface conductivity model of graphene, Journal of Applied Physics 103, 064302 (2008).
Majumdar et al. [2013] A. Majumdar, J. Kim, J. Vuckovic, and F. Wang, Electrical control of silicon photonic crystal cavity by graphene, Nano Letters 13, 515 (2013).
Satoshi et al. [2023] S. Satoshi, S. Otsuka, Y. Moritake, T. Yoda, T. Uemura, M. Ono, E. Kuramochi, and M. Notomi, Non-hermitian chirality and topological properties of graphene-loaded photonic crystals, in CLEO 2023 (Optica Publishing Group, 2023) p. FM4B.3.
Shalaev et al. [2019] M. I. Shalaev, W. Walasik, A. Tsukernik, Y. Xu, and N. M. Litchinitser, Robust topologically protected transport in photonic crystals at telecommunication wavelengths, Nature Nanotechnology 14, 31 (2019).
Sakoda and Sakoda [2005] K. Sakoda and K. Sakoda, Optical properties of photonic crystals, Vol. 2 (Springer, 2005).

EP restoration and fast-light edge states in photonic crystal waveguide with glide and time reversal symmetry

Abstract

I Introduction

II Design of non-Hermitian photonic crystal waveguides with glide and time reversal symmetry

III Reconstruction of exact EPs based on perturbation theory

III.1 Effective Hamiltonian of edge states near the Dirac point

III.2 Restoration of EP by adjusting the real part of the permittivity

IV Superluminal edge states and enlarged propagation length on photonic crystal waveguide

V implementation into the graphene-loaded photonic crystal slab

VI Conclusion

Acknowledgements.

Appendix A The design of heterostructure in other lattice or hole shapes

Appendix B HzH_{z} distribution of the eigenmodes

Appendix C Derivation of the 𝒌⋅𝒑\bm{k}\cdot\bm{p} Hamiltonian

Appendix D The EP condition with the hole shift in 𝒌⋅𝒑\bm{k}\cdot\bm{p} Hamiltonian

Appendix E Restoration of EPs and group velocity in graphene-Loaded structures with misalignment

References

EP restoration and fast-light edge states in photonic crystal waveguide
with glide and time reversal symmetry

Appendix B $H_{z}$ distribution of the eigenmodes

Appendix C Derivation of the $\bm{k}\cdot\bm{p}$ Hamiltonian

Appendix D The EP condition with the hole shift in $\bm{k}\cdot\bm{p}$ Hamiltonian