First-principle calculation of Chern number in gyrotropic photonic crystals

Ran Zhao; Ran Zhao; Guo-Da Xie; Guo-Da Xie; Menglin L. N. Chen; Zhihao Lan; Zhixiang Huang; Zhixiang Huang; Wei E. I. Sha

doi:10.1364/OE.380077

1. Introduction

Topology studies the invariant properties of geometry under continuous deformation [1]. While in mathematics, topological invariants are commonly used to classify topological spaces, in topological physics, topological invariants are explored to distinguish the bulk properties of materials. If physical observables can be expressed as topological invariants, they can only vary discretely and will not be affected by small perturbations of system parameters.

With the discovery of the quantum Hall effects [2,3] and recent advances in the study of topological insulators [4], topological phases of matter have attracted great attention in condensed matter physics. In 2005, Haldane and Raghu transferred the key feature of quantum Hall effect in quantum mechanics to classical electromagnetics [5] and soon after, it was numerically and experimentally verified by using photonic crystals (PCs) [6,7].

Chern number in a photonic system is defined on the dispersion bands in wave-vector space. For a two-dimensional (2D) periodic system, the Chern number is the integration of the Berry curvature over the first Brillouin zone. Once the Chern number is calculated, the topological properties of the system can be identified (trivial or non-trivial). In addition, the Chern number can be used to explain the phenomenon of “topological protection” of the edge state transmission in photonic topological insulators (PTIs). Therefore, the accurate computation of the Chern number is of crucial importance in the PTI design [8–10].

In this paper, we propose a first-principle computation method for Chern number calculation, based on the finite-difference frequency-domain (FDFD) method [11–14]. Compared with the commonly used plane-wave expansion (PWE) method [15–17], the FDFD method is not only accurate and stable but also has lower computational complexity. Firstly, the FDFD method is used to compute the band structure of 2D gyrotropic PCs by solving the generalized eigenvalue equations derived from Maxwell’s equations. Then, the formulae for the numerical calculation of the Chern number are derived in the discretized first Brillouin zone. At last, numerical examples are given to demonstrate the accuracy of the proposed method. Note, in this work, we will focus on 2D PCs with lossless, non-dispersive, local materials, which is a common practice in computing Chern numbers. The effects of material dispersion and loss in a real microwave ferrite has been discussed in [6]. One remark we would like to make is that as the Chern number is a global property of the energy band, it does not depend on the local changes of the band caused by material dispersion as long as the band gap is still open, which also means that the Chern number has a certain degree of built-in robustness against material dispersion. Nevertheless, the extension of the current method to fully include the effect of material dispersion on the energy bands and Chern number calculations is an interesting direction for future research.

2. Generalized eigenvalue problem using FD method

2.1 Maxwell’s equations in the generalized coordinates

Firstly, the Maxwell’s curl equations in free space in the generalized coordinate (Fig. 1) are expressed by

(1)$$\def\adots{\unicode[Times]{x22F0}}\nabla_{q}\times\hat{\mathbf{H}}=\imath k_0\mathbf{\hat{\epsilon}}(\mathbf{r})\hat{\mathbf{E}},\quad \nabla_{q}\times\hat{\mathbf{E}}=-\imath k_0\mathbf{\hat{\mu}}(\mathbf{r})\hat{\mathbf{H}}.$$

Here, $\nabla _{q}$ is the partial differential operator in the generalized coordinate with three unit vectors $\hat {\boldsymbol {u}}_q(x, y, z),(q=1,2,3)$. $\hat {\mathbf {E}}$ and $\hat {\mathbf {H}}$ are the normalized vector fields with $\hat {E}_i=Q_i\sqrt {\frac {\epsilon _0}{\mu _0}}E_i, \quad \hat {H}_i=Q_iH_i(i=1,2,3)$. Actually, in the numerical implementation, $Q_i$s are the sizes of the discrete grids along each direction [19]. $k_0$ is the free-space wavenumber. $\hat {\epsilon }$ and $\hat {\mu }$ are the relative permittivity and permeability in the form of $3\times 3$ tensor which are written as

(2)$$\mathbf{\hat{\epsilon}}_{ij}(\mathbf{r})=[\bar{\bar{\epsilon}}(\mathbf{r})\cdot\mathbf{g}]_{ij} \cdot|\hat{\boldsymbol{u}}_1\cdot \hat{\boldsymbol{u}}_2\times{\hat{\boldsymbol{u}}_3}|\frac{Q_1Q_2Q_3}{Q_iQ_j};\quad \mathbf{\hat{\mu}}_{ij}(\mathbf{r})=[\bar{\bar{\mu}}(\mathbf{r})\cdot\mathbf{g}]_{ij} \cdot|\hat{\boldsymbol{u}}_1\cdot \hat{\boldsymbol{u}}_2\times{\hat{\boldsymbol{u}}_3}|\frac{Q_1Q_2Q_3}{Q_iQ_j};$$

where $\bar {\bar {\epsilon }}(\mathbf {r})$ and $\bar {\bar {\mu }}(\mathbf {r})$ are the original relative permittivity and permeability in the Cartesian basis. $\mathbf {g}$ is the metric tensor matrix which can be expressed as,

(3)$$\mathbf{g}= \left[ \begin{array}{cccc} \hat{\boldsymbol{u}}_1\cdot \hat{\boldsymbol{u}}_1 & \hat{\boldsymbol{u}}_1\cdot \hat{\boldsymbol{u}}_2 & \hat{\boldsymbol{u}}_1\cdot \hat{\boldsymbol{u}}_3\\ \hat{\boldsymbol{u}}_2\cdot \hat{\boldsymbol{u}}_1 & \hat{\boldsymbol{u}}_2\cdot \hat{\boldsymbol{u}}_2 & \hat{\boldsymbol{u}}_2\cdot \hat{\boldsymbol{u}}_3\\ \hat{\boldsymbol{u}}_3\cdot \hat{\boldsymbol{u}}_1 & \hat{\boldsymbol{u}}_3\cdot \hat{\boldsymbol{u}}_2 & \hat{\boldsymbol{u}}_3\cdot \hat{\boldsymbol{u}}_3\\ \end{array} \right]^{-1}.$$

Fig. 1. The 2D Yee’s grid in general coordinates. The dashed arrows are where the periodic boundaries apply.

Download Full Size | PDF

With the metric tensor matrix, the length of a vector in the generalized coordinate is calculated by

(4)$$\left| \mathbf{r} \right|^2= \mathbf{r}^\textrm{T} \cdot\mathbf{g}^{-1}\cdot \mathbf{r}.$$

2.2 Bloch’s theorem

In 2D PCs, due to the translational symmetry, the modes can be written in the Bloch form:

(5)$$\begin{array}{c} \mathbf{E}_{\left(n, k_{z}, \mathbf{k}_{ \|}\right)}(\mathbf{r})=e^{-i \mathbf{k}_{ \|} \cdot \boldsymbol{\rho}} e^{-i k_{z} z} \mathbf{u}_{e \left(n, k_{z}, \mathbf{k}_{ \|}\right)}(\boldsymbol{\rho}),\\ \mathbf{H}_{\left(n, k_{z}, \mathbf{k}_{ \|}\right)}(\mathbf{r})=e^{-i \mathbf{k}_{ \|} \cdot \boldsymbol{\rho}} e^{-i k_{z} z} \mathbf{u}_{h \left(n, k_{z}, \mathbf{k}_{ \|}\right)}(\boldsymbol{\rho}). \end{array}$$

Here, $n$ is the band number, $\boldsymbol {\rho }$ is the mapping of vector $\mathbf {r}$ in the $(u_1, u_2)$ plane. $\mathbf {u}_e( \boldsymbol {\rho })$ and $\mathbf {u}_h(\boldsymbol {\rho })$ are periodic functions satisfying $\mathbf {u}_h(\boldsymbol {\rho })=\mathbf {u}_h(\boldsymbol {\rho }+\mathbf {R})$, $\mathbf {u}_e(\boldsymbol {\rho })=\mathbf {u}_e(\boldsymbol {\rho }+\mathbf {R})$, where $\mathbf {R}$ stands for any lattice vectors.

In 2D PCs, we only need to investigate the in-plane propagating modes, i.e. $k_z=0$. Hence, the mode can be rewritten as

(6)$$\begin{array}{c} \mathbf{E}(\boldsymbol{\rho})=e^{-i \mathbf{k}_{ \|} \cdot \boldsymbol{\rho}} \mathbf{u}_{e}(\boldsymbol{\rho}),\\ \mathbf{H}(\boldsymbol{\rho})=e^{-i \mathbf{k}_{ \|} \cdot \boldsymbol{\rho}} \mathbf{u}_{h}(\boldsymbol{\rho}). \end{array}$$

Then, the Bloch’s theorem can be written as

(7)$$\begin{array}{c} \mathbf{E}(\boldsymbol{\rho}+\mathbf{R})=e^{-i \mathbf{k}_{ \|} \cdot \mathbf{R}} \mathbf{E}(\boldsymbol{\rho}),\\ \mathbf{H}(\boldsymbol{\rho}+\mathbf{R})=e^{-i \mathbf{k}_{ \|} \cdot \mathbf{R} } \mathbf{H}(\boldsymbol{\rho}). \end{array}$$

2.3 Discretization of the eigenvalue problems

From the Maxwell’s curl equations (Eq. (1)), we can derive the governing equations for the modes of the PCs:

(8)$$\nabla\times\hat{\mathbf{\mu}}^{-1}(\textbf{r})\nabla\times \textbf{E}(\textbf{r}) = \omega^2 \hat{\epsilon}(\textbf{r}) \textbf{E}(\textbf{r}),$$

(9)$$\nabla\times\hat{\mathbf{\epsilon}}^{-1}(\textbf{r})\nabla\times \textbf{H}(\textbf{r}) = \omega^2 \hat{\mu}(\textbf{r}) \textbf{H}(\textbf{r}).$$

In 2D PCs, the modes can be separated into two distinct polarizations, transverse-magnetic (TM) modes with $z$-polarized electric fields and transverse-electric (TE) modes with $z$-polarized magnetic fields. The generalized eigenvalue matrix equations for the TM and TE modes of the 2D PCs can be discretized as,

(10)$$\left\{U_{1}\left(\mu_{21}^{-1} V_{2}-\mu_{22}^{-1} V_{1}\right)-U_{2}\left(\mu_{11}^{-1} V_{2}-\mu_{12}^{-1} V_{1}\right)\right\} \hat{E}_{z}=\varepsilon_{33} k_{0}^{2} \hat{E}_{z}$$

(11)$$\left\{V_{1}\left(\varepsilon_{21}^{-1} U_{2}-\varepsilon_{22}^{-1} U_{1}\right)-V_{2}\left(\varepsilon_{11}^{-1} U_{2}-\varepsilon_{12}^{-1} U_{1}\right)\right\} \hat{H}_{z}=\mu_{33} k_{0}^{2} \hat{H}_{z}$$

(12)$${{U}_{1}}=\frac{1}{{{Q}_{1}}}\left[ \begin{matrix} -1 & 1 & {} & {} & {} & {} & {} & {} & {} \\ {} & \ddots & \ddots & {} & {} & {} & {} & {} & {} \\ {{u}_{x}} & {} & -1 & 0 & {} & {} & {} & {} & {} \\ {} & {} & {} & -1 & 1 & {} & {} & {} & {} \\ {} & {} & {} & {} & \ddots & \ddots & {} & {} & {} \\ {} & {} & {} & {{u}_{x}} & {} & -1 & 0 & {} & {} \\ {} & {} & {} & {} & {} & {} & \ddots & 1 & {} \\ {} & {} & {} & {} & {} & {} & {} & -1 & \ddots \\ {} & {} & {} & {} & {} & {} & {{u}_{x}} & {} & -1 \\ \end{matrix} \right]$$

(13)$${{V}_{1}}=\frac{1}{{{Q}_{1}}}\left[ \begin{matrix} 1 & {} & {{v}_{x}} & {} & {} & {} & {} & {} & {} \\ -1 & 1 & {} & {} & {} & {} & {} & {} & {} \\ {} & \ddots & \ddots & {} & {} & {} & {} & {} & {} \\ {} & {} & 0 & 1 & {} & {{v}_{x}} & {} & {} & {} \\ {} & {} & {} & -1 & 1 & {} & {} & {} & {} \\ {} & {} & {} & {} & \ddots & \ddots & {} & {} & {} \\ {} & {} & {} & {} & {} & 0 & 1 & {} & {{v}_{x}} \\ {} & {} & {} & {} & {} & {} & -1 & \ddots & {} \\ {} & {} & {} & {} & {} & {} & {} & \ddots & 1 \\ \end{matrix} \right]$$

(14)$${{U}_{2}}=\frac{1}{{{Q}_{2}}}\left[ \begin{matrix} -1 & {} & {} & 1 & {} & {} & {} & {} & {} \\ {} & -1 & {} & {} & 1 & {} & {} & {} & {} \\ {} & {} & \ddots & {} & {} & \ddots & {} & {} & {} \\ {} & {} & {} & -1 & {} & {} & \ddots & {} & {} \\ {} & {} & {} & {} & \ddots & {} & {} & 1 & {} \\ {} & {} & {} & {} & {} & -1 & {} & {} & 1 \\ {{u}_{y}} & {} & {} & {} & {} & {} & -1 & {} & {} \\ {} & \ddots & {} & {} & {} & {} & {} & \ddots & {} \\ {} & {} & {{u}_{y}} & {} & {} & {} & {} & {} & -1 \\ \end{matrix} \right]$$

(15)$${{V}_{2}}=\frac{1}{{{Q}_{2}}}\left[ \begin{matrix} 1 & {} & {} & {} & {} & {} & {{v}_{y}} & {} & {} \\ {} & \ddots & {} & {} & {} & {} & {} & {{v}_{y}} & {} \\ {} & {} & 1 & {} & {} & {} & {} & {} & {{v}_{y}} \\ -1 & {} & {} & 1 & {} & {} & {} & {} & {} \\ {} & {} & {} & {} & \ddots & {} & {} & {} & {} \\ {} & {} & \ddots & {} & {} & {} & {} & {} & {} \\ {} & {} & {} & \ddots & {} & {} & \ddots & {} & {} \\ {} & {} & {} & {} & {} & {} & {} & {} & {} \\ {} & {} & {} & {} & {} & -1 & {} & {} & 1 \\ \end{matrix} \right]$$

The grid points at the boundaries are treated by the Bloch’s theorem, consequently

(16)$$u_{x}=\exp \left(i \mathbf{k} \cdot a_{1} \hat{\boldsymbol{u}}_{1}\right), v_{x}=-\exp \left(-i \mathbf{k} \cdot a_{1} \hat{\boldsymbol{u}}_{1}\right), \\ u_{y}=\exp \left(i \mathbf{k} \cdot a_{2} \hat{\boldsymbol{u}}_{2}\right), v_{y}=-\exp \left(-i \mathbf{k} \cdot a_{2} \hat{\boldsymbol{u}}_{2}\right),$$

where $a_1$, $a_2$ are the lengths of the unit cell along the directions of $\hat {\boldsymbol {u}}_{1}$ and $\hat {\boldsymbol {u}}_{2}$.

By solving these eigenvalue equations, the band structure together with the eigenstates of 2D PCs can be easily obtained. The matrices

(17)$$\Theta_\mathbf{E}=\left\{U_{1}\left(\mu_{21}^{-1} V_{2}-\mu_{22}^{-1} V_{1}\right)-U_{2}\left(\mu_{11}^{-1} V_{2}-\mu_{12}^{-1} V_{1}\right)\right\}$$

(18)$$\Theta_\mathbf{H}=\left\{V_{1}\left(\varepsilon_{21}^{-1} U_{2}-\varepsilon_{22}^{-1} U_{1}\right)-V_{2}\left(\varepsilon_{11}^{-1} U_{2}-\varepsilon_{12}^{-1} U_{1}\right)\right\}$$

are discretized from the operator $\nabla \times \hat {\mathbf {\mu }}^{-1}(\textbf {r})\nabla \times$ and $\nabla \times \hat {\mathbf {\epsilon }}^{-1}(\textbf {r})\nabla \times$. The condition of the electromagnetic modes $\mathbf {E}_{\omega _1}, \mathbf {E}_{\omega _2}$ (or $\mathbf {H}_{\omega _1}, \mathbf {H}_{\omega _2}$) orthogonality is the Hermitian inner product satisfies $\left \langle \mathbf {E}_{\omega _1} | \mathbf {E}_{\omega _2} \right \rangle =\int \epsilon (\mathbf {r}) \mathbf {E}_{\omega _1} ^{*} \cdot \mathbf {E}_{\omega _2} d^{2} \mathbf {r}=0$ or $\left \langle \mathbf {H}_{\omega _1} | \mathbf {H}_{\omega _2} \right \rangle =\int \mu (\mathbf {r}) \mathbf {H}_{\omega _1} ^{*} \cdot \mathbf {H}_{\omega _2} d^{2} \mathbf {r}=0$. The subscript ${\omega _1}$ and ${\omega _1}$ stand for the different eigen-frequency under a fixed wave-vector $\mathbf {k}$.

3. Chern number calculation

3.1 Single-band Chern number

Chern number is an integer which determines the topological classification of different materials or structures. Taking the TM mode of a 2D PC as an example, the Chern number of the $\textrm {n}^{th}$ band can be computed by integrating the Berry curvature over the first Brillouin zone as following [20],

(19)$$\begin{aligned}[b] C^{(n)} & = \frac{1}{2 \pi} \int_{\textrm{BZ}} \mathbf{F}_{n}(\mathbf{k}) \textrm{d} \mathbf{k}\\ & = \frac{1}{2 \pi} \int_{\textrm{BZ}} \nabla_\mathbf{k}\times\mathbf{A}_{n}(\mathbf{k}) \textrm{d} \mathbf{k}\\ & =\frac{1}{2 \pi i} \oint_{\partial \textrm{BZ}}\left\langle\mathbf{u}_{n,e,\mathbf{k}}\left|\nabla_{\mathbf{k}}\right| \mathbf{u}_{n,e,\mathbf{k}}\right\rangle \textrm{d} \mathbf{k}.\\ \end{aligned}$$

Here, the $\mathbf {F}_{n}(\mathbf {k})$ is the Berry curvature. $\mathbf {A}_{n}(\mathbf {k})=-i\left \langle \mathbf {u}_{n,e, \mathbf {k}}\left |\nabla _{\mathbf {k}}\right | \mathbf {u}_{n,e,\mathbf {k}}\right \rangle$ is the Berry connection. $\mathbf {u}_{n,e,\mathbf {k}}$ denotes the normalized eigenstate satisfying $\left . \left |\mathbf {u}_{n,e,\mathbf {k}}( \boldsymbol {\rho }) \right . \right \rangle =\left . \left |e^{i\mathbf {k} \cdot \boldsymbol {\rho }} \mathbf {E}_{n,e,\mathbf {k}}( \boldsymbol {\rho })\right . \right \rangle$ and $\left . \left |\mathbf {E}_{n,e,\mathbf {k}}( \boldsymbol {\rho })\right . \right \rangle$ is the eigenstate solved from Eq. (10).

For numerical discretization, the Chern number is calculated in the discretized Brillouin zone. In Fig. 2(a), an example for a square-shaped first Brillouin zone is shown. It is discretized into $4\times 4$ plaquettes. Then, the Chern number is calculated using

(20)$$\begin{aligned}[b] C^{(n)} & = \frac{1}{2 \pi} \int_{\textrm{BZ}} \mathbf{F}_{n}(\mathbf{k}) \textrm{d} \mathbf{k}\\ & =\frac{1}{2 \pi \textrm{i}} \oint_{\partial \textrm{BZ}}\left\langle\mathbf{u}_{n,e,\mathbf{k}}\left|\nabla_{\mathbf{k}}\right| \mathbf{u}_{n,e,\mathbf{k}}\right\rangle \textrm{d} \mathbf{k}\\ & = \frac{1}{2 \pi} \sum_{\mathbf{k}\in\textrm{BZ}} F_{\mathbf{k}}^{(n)} \Delta S_{\mathbf{k}} \\ & {= \frac{1}{2 \pi} \sum_{\mathbf{k}\in\textrm{BZ}} \operatorname{Im} \ln\left[{ U_{\mathbf{k}_1 \rightarrow \mathbf{k}_2}^{(n)} U_{\mathbf{k}_2 \rightarrow \mathbf{k}_3 }^{(n)} U_{\mathbf{k}_3\rightarrow \mathbf{k}_4}^{(n)} U_{\mathbf{k}_4 \rightarrow \mathbf{k}_1}^{(n)} }\right]} \end{aligned}$$

where $U_{\mathbf {k}_{\alpha } \rightarrow \mathbf {k}_{\beta }}^{(n)} \equiv \frac {\left \langle \mathbf {u}_{n,e,\mathbf {k}_{\alpha }} | \mathbf {u}_{n,e, \mathbf {k}_{\beta }}\right \rangle }{\left |\left \langle \mathbf {u}_{n,e,\mathbf {k}_{\alpha }} | \mathbf {u}_{n,e, \mathbf {k}_{\beta }}\right \rangle \right |}$, $\alpha ,\beta =1,2,3,4$, and $\left \langle \mathbf {u}_{n,e,\mathbf {k}_{\alpha }} | \mathbf {u}_{n,e, \mathbf {k}_{\beta }}\right \rangle =\int \epsilon (\mathbf {r}) \mathbf {u}_{n,e,\mathbf {k}_{\alpha }}(\mathbf {r}) ^{*} \cdot \mathbf {u}_{n,e, \mathbf {k}_{\beta }}(\mathbf {r}) d^{2} \mathbf {r}.$ $\mathbf {k}_1,\mathbf {k}_2,\mathbf {k}_3,\mathbf {k}_4$ are the vertices of the each plaquette. The integration is replaced by a summation $\sum _{\mathbf {k}\in \textrm {BZ}} F_{\mathbf {k}}^{(n)} \Delta S_{\mathbf {k}}$.

Fig. 2. Computational scheme for the calculation of Chern numbers. (a) A square-shaped first Brillouin zone; (b) Discretization of the first Brillouin zone, and the computational scheme for the Chern number.

Download Full Size | PDF

For the first Brillouin zone in different shape, similar discretization scheme can be used. It should be noted that to simplify the numerical calculation, we may need to use the equivalent Brillouin zone which shares the same size with the first Brillouin zone but in a parallelogram shape with its edges coinciding with the reciprocal lattice vector [18].

3.2 Composite Chern number

When the bands are degenerate, the Chern number can not be assigned to each band, so Eq. (20) will no longer be valid. For this case, these degenerate bands (taking the $\{n, n+1, \ldots , n+N-1\}$ bands as an example) jointly share a composite (first) Chern number, $C^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)}$ which is associated with the multiplets $\mathbf {u}_{n \mathbf {k}} \oplus \mathbf {u}_{n+1, \mathbf {k}} \oplus \cdots \oplus \mathbf {u}_{n+N-1, \mathbf {k}}$. It is calculated using [20,21]

(21)$$\begin{aligned}[b] C^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)} & = \frac{1}{2 \pi} \operatorname{Tr} \oint_{\partial \mathrm{BZ}} \mathbf{A}_{\mathbf{k}}^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)} \mathrm{d} \mathbf{k} \\ & =\frac{1}{2 \pi \mathrm{i}} \oint_{\partial \mathrm{BZ}} \left[ \begin{matrix} {\ddots} & {\vdots} & {\adots } \\ {\cdots} & {\left\langle\mathbf{u}_{n_l,e,\mathbf{k}}\left|\nabla_{\mathbf{k}}\right| \mathbf{u}_{n_m,e,\mathbf{k}}\right\rangle } & {\cdots} \\ {\adots} & {\vdots} & {\ddots}\\ \end{matrix} \right] \mathrm{d} \mathbf{k} \end{aligned}$$

Here, $\mathbf {A}_{\mathbf {k}}^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)}$ is an $N\times N$ matrix with the $(l,m)$ element, $-\textrm {i}{\left \langle \mathbf {u}_{n_l,e,{\mathbf {k}}}\left |\nabla _{\mathbf {k}}\right | \mathbf {u}_{n_m,e,{\mathbf {k}}}\right \rangle }$, and $n_{l}=n+l-1, n_{m}=n+m-1$.

By applying Stokes’ theorem, the composite Chern number can be rewritten in terms of the composite Berry flux:

(22)$$\begin{aligned}[b] C^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)} & =\frac{1}{2 \pi} \int_{\textrm{BZ}} F_{\mathbf{k}}^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)} \textrm{d} S_{\mathbf{k}} \\ & =\frac{1}{2 \pi} \int_{\textrm{BZ}} \operatorname{Tr} \nabla_{\mathbf{k}} \times \mathbf{A}_{\mathbf{k}}^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)} \textrm{d} S_{\mathbf{k}}, \end{aligned}$$

(23)$$\begin{aligned}[b] & F_{\mathbf{k}}^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)} \textrm{d} S_{\mathbf{k}}\\ & =\operatorname{Im} \textrm{Tr} \left\{ \ln \left[ \mathbf{U}_{\mathbf{k}_1 \rightarrow \mathbf{k}_2}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_2 \rightarrow \mathbf{k}_3}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_3 \rightarrow \mathbf{k}_4}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_4 \rightarrow \mathbf{k}_1}^{(n \oplus \cdots \oplus n+N-1)} \right] \right\}\\ & =\operatorname{Im} \ln \left\{ \det \left[ \mathbf{U}_{\mathbf{k}_1 \rightarrow \mathbf{k}_2}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_2 \rightarrow \mathbf{k}_3}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_3 \rightarrow \mathbf{k}_4}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_4 \rightarrow \mathbf{k}_1}^{(n \oplus \cdots \oplus n+N-1)} \right] \right\}. \end{aligned}$$

The element $(l, m)$ of the link matrix $\mathbf {U}_{\mathbf {k}_{\alpha } \rightarrow \mathbf {k}_{\beta }}^{(n \oplus \cdots \oplus n+N-1)}$ is

(24)$$\mathbf{U}_{\mathbf{k}_{\alpha} \rightarrow \mathbf{k}_{\beta}}^{(n \oplus \cdots \oplus n+N-1)}\left(l,m\right)= \frac{\left\langle\mathbf{u}_{n_l,e,\mathbf{k}_{\alpha}} | \mathbf{u}_{n_m,e, \mathbf{k}_{\beta}}\right\rangle}{\left|\left\langle\mathbf{u}_{n_l,e,\mathbf{k}_{\alpha}} | \mathbf{u}_{n_m,e, \mathbf{k}_{\beta}}\right\rangle\right|}.$$

Different from the single-band Chern number, in the composite Chern number, the link matrix $\mathbf {U}_{\mathbf {k}_{\alpha } \rightarrow \mathbf {k}_{\beta }}^{(n \oplus \cdots \oplus n+N-1)}$ is used to replace the link variable $U_{\mathbf {k}_{\alpha } \rightarrow \mathbf {k}_{\beta }}^{(n)}$. In the discretized Brillouin zone, the Chern number is expressed as

(25)$$\begin{aligned}[b] & C^{(n \oplus n+1 \oplus \cdots \oplus n+N-1)} \\ & = \frac{1}{2 \pi}\sum_{\textrm{BZ}} \operatorname{Im} \ln \left\{ \det \left[ \mathbf{U}_{\mathbf{k}_1 \rightarrow \mathbf{k}_2}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_2 \rightarrow \mathbf{k}_3}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_3 \rightarrow \mathbf{k}_4}^{(n \oplus \cdots \oplus n+N-1)} \mathbf{U}_{\mathbf{k}_4 \rightarrow \mathbf{k}_1}^{(n \oplus \cdots \oplus n+N-1)} \right] \right\} \end{aligned}$$

As for the calculation of the Chern number of TE modes, the derived formulae are still valid, except that the normalized eigenstates now become the normalized magnetic field, $\mathbf {u}_{n,h,\mathbf {k}_{\alpha }}(\mathbf {r})$. Finally, we would like to note that the algorithm for the Chern number calculations described above guarantees the calculated Chern number is strictly an integer for arbitrary discretization spacing (for the detailed proof, see [20]). However, this does not mean that any discretization will give the correct Chern numbers, but rather, convergence to the correct Chern numbers still requires a sufficient sampling of the Brillouin zone, which can be achieved using rather coarse discretization as we will show in the following section.

4. Numerical results

In the first numerical example, a 2D PC composed of Yttrium-Iron-Garnet (YIG) rods ($\epsilon =15\epsilon _0$) in a square lattice is analyzed. The radius of the rods is $0.11a$, where $a$ is the lattice constant. An external direct current (DC) magnetic field is applied in the $z$ direction (out-of-plane) to break the time-reversal symmetry. The permeability is a tensor matrix written as

(26)$$\bar{\bar{\mu}}=\left[ \begin{array}{ccc}{\mu} & {i \kappa} & {0} \\ {-i \kappa} & {\mu} & {0} \\ {0} & {0} & {\mu_{0}}\end{array}\right],\quad \mu=1+\frac{\omega_m\omega_0}{\omega^2_0-\omega^2},\quad \kappa=\frac{\omega_m\omega}{\omega^2_0-\omega^2},$$

where $\omega _0=\gamma H_0$ is the precession frequency and $\omega _m=4\pi {\gamma }{M_s}$. When the external magnetic field is 1600 Gauss at 4.28 GHz, $\kappa =12.4 \mu _{0}$, and $\mu =14 \mu _{0}$ [6,24].

Firstly, we calculate the band structure of the PC using the FDFD method. In order to ensure the accuracy, the mesh size is set to be $a/100$. The calculated band structures for the TM modes and TE modes are plotted in Figs. 3(a) and Figs. 3(b), respectively. In Fig. 3(a), $C^{(1)}_{\textrm {TM}}, C^{(2)}_{\textrm {TM}}, C^{(3)}_{\textrm {TM}}$ stand for the Chern numbers of the first, second and third bands. Then, we investigate the computational accuracy of the Chern number by discretizing the Brillouin zone with different numbers of plaquettes. The computational results are shown in Table 1. We can see that the computational results of Chern numbers are converged when the number of plaquettes is larger than $4 \times 4$. By further checking the values of the Chern numbers and composite Chern numbers, we find $C^{(1 \oplus 2)}_{\textrm {TM}}=C^{(1)}_{\textrm {TM}}+C^{(2)}_{\textrm {TM}}$, and $C^{(1 \oplus 2\oplus 3)}_{\textrm {TM}}=C^{(1)}_{\textrm {TM}}+C^{(2)}_{\textrm {TM}}+C^{(3)}_{\textrm {TM}}$, which is as expected.

Fig. 3. The Band structures of a 2D PC composed of YIG rods in a square lattice when a 1600 Gauss $+\textrm {z}$ DC magnetic field is applied. $C^{(1)}_{\textrm {TM}}$, $C^{(2)}_{\textrm {TM}}, C^{(3)}_{\textrm {TM}}$ denote the Chern numbers of the TM modes of the associated bands.

Download Full Size | PDF

Table 1. The TM-mode Chern numbers and composite Chern numbers of different bands calculated with different mesh sizes of the first Brillouin zone.

View Table | View all tables in this article

In the second numerical example, a 2D PC composed of gyrotropic rods in a square lattice is analyzed. The geometric parameters are the same as these of the first example ($r=0.11a$). The relative permittivity and permeability are expressed as following,

(27)$$\bar{\bar{\varepsilon}}=\left[ \begin{array}{ccc}{\varepsilon_{d}} & {i \varepsilon_{f}} & {0} \\ {-i \varepsilon_{f}} & {\varepsilon_{d}} & {0} \\ {0} & {0} & {\varepsilon_\perp}\end{array}\right], \bar{\bar{\mu}}=\left[ \begin{array}{ccc}{\mu_{d}} & { i\mu_{f}} & {0} \\ {-i \mu_{f}} & {\mu_{d}} & {0} \\ {0} & {0} & {\mu_{\perp}}\end{array}\right].$$

Here, the parameters are chosen to satisfy the relationship $\varepsilon _{d}=\mu _{d}, \varepsilon _{f}=-\mu _{f}, \varepsilon _{\perp }=\mu _{\perp }$. Referring to the YIG material under an external magnetic field at 4.28 GHz, the corresponding parameters are set as $\varepsilon _{d}=14, \varepsilon _{f}=-12.4$ and $\mu _{\perp }=15$, $\mu _{d}=14, \mu _{f}=12.4$ and $\varepsilon _{\perp }=15$ [22].

The band structures of TM and TE modes are shown in Fig. 4. Different from the first numerical example, by using the gyrotropic material, there are band gaps for both the TM and TE modes. The $C^{(1)}_{\textrm {TM/TE}}$, $C^{(2)}_{\textrm {TM/TE}}$, $C^{(3)}_{\textrm {TM/TE}}$ stand for the Chern numbers of the first, second and third bands for the TM/TE mode. As shown in the Table 2 and Table 3, $4\times 4$ meshes in the Brillouin zone are also sufficient for the Chern number computation for TM/TE modes. With the increase of the discretization density, the calculated results are stable.

Fig. 4. The Band structures of a 2D PC composed of gyrotopic rods in a square lattice when a 1600 Gauss $+\textrm {z}$ DC magnetic field is applied. $C^{(1)}_{\textrm {TM/TE}}$, $C^{(2)}_{\textrm {TM/TE}}$, $C^{(3)}_{\textrm {TM/TE}}$ denote the Chern numbers of the TM/TE modes of the associated bands.

Download Full Size | PDF

Table 2. The TM-mode Chern numbers and composite Chern numbers of different bands calculated with different mesh sizes of the first Brillouin zone.

View Table | View all tables in this article

Table 3. The TE-mode Chern numbers and composite Chern numbers of different bands calculated with different mesh sizes of the first Brillouin zone.

View Table | View all tables in this article

In the third numerical example, a 2D PC composed of YIG rods in a honeycomb lattice [23] is analyzed. The unit cell of the lattice is shown in the Fig. 5. The corresponding first Brillouin zone (hexagon region) constructed in the reciprocal lattice is also shown in the Fig. 5. As shown in this figure, for the convenience of discretization, the hexagon region needs to be equivalently transformed into a rhombus region, which can be discretized into smaller rhombus cells straightforwardly. The detailed reshaping and discretization procedure for the Brillouin zone can be found in [18]. In Fig. 6(a), the first two TM bands with zero DC magnetic field are given. It can be seen that the two bands are degenerate at the $\mathbf {K}$ point. By applying a DC magnetic field, the time-reversal symmetry is broken and the degeneracy is lifted, as shown in the Fig. 6(b). The calculated Chern numbers and composite Chern numbers corresponding to the two bands are shown in the Table 4. The non-zero Chern numbers can be used to explain the edge-state in [23] at 7.5 GHz. Also, as can be seen from this table, a very coarse mesh ($2\times 2$) can guarantee the accuracy of the results.

Fig. 5. The unit cell and first Brillouin zone of the honeycomb lattice.

Download Full Size | PDF

Fig. 6. The band structure of 2D honeycomb lattice without and with DC magnetic field.

Download Full Size | PDF

Table 4. The honeycomb lattice’s Chern numbers and composite Chern numbers of the first two bands calculated with different mesh sizes of the first Brillouin zone.

View Table | View all tables in this article

5. Conclusion

In conclusion, a first-principle computational method for the Chern number of 2D gyrotropic PCs has been developed [25]. Firstly, the full-wave FDFD method was used to solve the generalized eigenvalue problem based on the Maxwell’s equations. Then, the band structures of the PCs with different material properties were analyzed. The Chern numbers were calculated by integrating the Berry curvature over the first Brillouin zone. All the three numerical examples have demonstrated the accuracy and convergence of this method. Importantly, the method proposed in the current work could be potentially applied to 3D PCs [26] or 3D planar PCs when one cannot separate the TE and TM polarizations, in which case the full vectorial eigenmodes could be used for the Chern number calculations [5].

Funding

Natural Science Foundation of Anhui Province (1808085QF183); Key Natural Science Project of Anhui Provincial Education Department (KJ2018A0015); National Natural Science Foundation of China (61701424, 61722101, 61801002, 61901087, 61975177).

Acknowledgments

After this manuscript was accepted, we became aware of two closely related works [27,28] using commercial software COMSOL Multiphysics and open-source software package MIT Photonic Bands (MPB) for the computation of topological invariants in 2D photonic crystals. Compared with these two works, in which eigensolutions of Maxwell's equations are obtained by either COMSOL or MPB, in our method, the eigensolutions are solved directly via the first-principle FDFD scheme, thus effectively avoiding the tedious data extraction and postprocessing. Moreover, our code is freely available for download [25].

References

1. M. A. Armstrong, Basic Topology (Springer Science & Business Media, 2013).

2. M. Z. Hasan and C. L. Kane, “Colloquium: topological insulators,” Rev. Mod. Phys. 82(4), 3045–3067 (2010). [CrossRef]

3. X. L. Qi and S. C. Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83(4), 1057–1110 (2011). [CrossRef]

4. K. von Klitzing, “The quantized Hall effect,” Rev. Mod. Phys. 58(3), 519–531 (1986). [CrossRef]

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

6. Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, “Reflection-free one-way edge modes in a gyromagnetic photonic crystal,” Phys. Rev. Lett. 100(1), 013905 (2008). [CrossRef]

7. Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, “Observation of unidirectional backscattering-immune topological electromagnetic states,” Nature 461(7265), 772–775 (2009). [CrossRef]

8. L. Lu, J. D. Joannopoulos, and M. Soljačić, “Topological photonics,” Nat. Photonics 8(11), 821–829 (2014). [CrossRef]

9. S. A. Hassani Gangaraj, M. G. Silveirinha, and G. W. Hanson, “Berry phase, Berry connection, and Chern number for a continuum bianisotropic material from a classical electromagnetics perspective,” IEEE J. Multiscale Multiphys. Comput. Tech. 2, 3–17 (2017). [CrossRef]

10. B. Y. Xie, H. F. Wang, X. Y. Zhu, M. H. Lu, Z. D. Wang, and Y. F. Chen, “Photonics meets topology,” Opt. Express 26(19), 24531–24550 (2018). [CrossRef]

11. Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructuredd optical fibers,” Opt. Express 10(17), 853 (2002). [CrossRef]

12. S. Guo, F. Wu, S. Albin, and R. S Rogowski, “Photonic band gap analysis using finite-difference frequency domain method,” Opt. Express 12(8), 1741 (2004). [CrossRef]

13. M. L. N. Chen, L. J. Jiang, and W. E. I. Sha, “Generation of orbital angular momentum by a point defect in photonic crystals,” Phys. Rev. Appl. 10(1), 014034 (2018). [CrossRef]

14. M. Fang, Z. Huang, W. E. I. Sha, and X. Wu, “Maxwell-Hydrodynamic Model for Simulating Nonlinear Terahertz Generation from Plasmonic Metasurfaces,” IEEE J. Multiscale Multiphys. Comput. Tech. 2, 194–201 (2017). [CrossRef]

15. K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic dielectric structures,” Phys. Rev. Lett. 65(25), 3152–3155 (1990). [CrossRef]

16. R. D. Meade, A. M. Rappe, K. D. Brommer, J. D. Joannopoulos, and O. L. Alerhand, “Accurate theoretical analysis of photonic band gap materials,” Phys. Rev. B 48(11), 8434–8437 (1993). [CrossRef]

17. K. M. Leung and Y. F. Liu, “Full vector wave calculation of photonic band structures in FCC dielectric media,” Phys. Rev. Lett. 65(21), 2646–2649 (1990). [CrossRef]

18. D. Jin, T. Christensen, M. Soljacic, N. X. Fang, L. Lu, and X. Zhang, “Infrared topological plasmons in graphene,” Phys. Rev. Lett. 118(24), 245301 (2017). [CrossRef]

19. A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equations,” J. Mod. Opt. 43(4), 773–793 (1996). [CrossRef]

20. T. Fukui, H. Suzuki, and Y. Hatsugai, “Chern numbers in a discretized Brillouin zone: Efficient method to compute (spin) Hall conductances,” J. Phys. Soc. Jpn. 74(6), 1674–1677 (2005). [CrossRef]

21. P. Wang, L. Lu, and K. Bertoldi, “Topological phononic crystals with one-way elastic edge waves,” Phys. Rev. Lett. 115(10), 104302 (2015). [CrossRef]

22. X. C. Sun, C. He, and X. P. Liu, “Photonic topological states in a two-dimensional gyrotropic photonic crystal,” Crystals 9(3), 137 (2019). [CrossRef]

23. X. Ao, Z. Lin, and C. T. Chan, “One-way edge mode in a magneto-optical honeycomb photonic crystal,” Phys. Rev. B 80(3), 033105 (2009). [CrossRef]

24. D. M. Pozar, Microwave Engineering, 2nd ed. (Wiley, 497–517, 1998).

25. http://www.zjuisee.zju.edu.cn/weisha/Publications/Files/Chern_number.zip.

26. S. Oono, T. Kariyado, and Y. Hatsugai, “Section Chern number for a three-dimensional photonic crystal and the bulk-edge correspondence,” Phys. Rev. B 94(12), 125125 (2016). [CrossRef]

27. C. Wang, H. Zhang, H. Yuan, J. Zhong, and C. Lu, “Universal numerical calculation method for the Berry curvature and Chern numbers of typical topological photonic crystals,” Front. Optoelectron., https://doi.org/10.1007/s12200-019-0963-9 (2020). [CrossRef]

28. Y. Li, Z.-A. Xu, and Y. Hatsugai, “Exploring topological superconductivity in topological materials,” Adv. Quantum Technol. 2, 1800112 (2019). [CrossRef]

MeshChern	$C_{TM}^{(2)}$	$C_{TM}^{(3)}$	$C_{TM}^{(1 \oplus 2)}$	$C_{TM}^{(1 \oplus 2 \oplus 3)}$
$2 \times 2$	1	0	1	−1
$4 \times 4$	1	−2	1	−1
$8 \times 8$	1	−2	1	−1
$16 \times 16$	1	−2	1	−1

MeshChern	$C_{TM}^{(2)}$	$C_{TM}^{(3)}$	$C_{TM}^{(1 \oplus 2)}$	$C_{TM}^{(1 \oplus 2 \oplus 3)}$
$2 \times 2$	1	0	1	−1
$4 \times 4$	1	−2	1	−1
$8 \times 8$	1	−2	1	−1
$16 \times 16$	1	−2	1	−1

MeshChern	$C_{TE}^{(2)}$	$C_{TE}^{(3)}$	$C_{TE}^{(1 \oplus 2)}$	$C_{TE}^{(1 \oplus 2 \oplus 3)}$
$2 \times 2$	−1	0	−1	1
$4 \times 4$	−1	2	−1	1
$8 \times 8$	−1	2	−1	1
$16 \times 16$	−1	2	−1	1

MeshChern	$C^{(1)}$	$C^{(2)}$
$2 \times 2$	1	−1
$4 \times 4$	1	−1
$8 \times 8$	1	−1
$16 \times 16$	1	−1

MeshChern	$C_{TM}^{(2)}$	$C_{TM}^{(3)}$	$C_{TM}^{(1 \oplus 2)}$	$C_{TM}^{(1 \oplus 2 \oplus 3)}$
$2 \times 2$	1	0	1	−1
$4 \times 4$	1	−2	1	−1
$8 \times 8$	1	−2	1	−1
$16 \times 16$	1	−2	1	−1

First-principle calculation of Chern number in gyrotropic photonic crystals

Abstract

1. Introduction

2. Generalized eigenvalue problem using FD method

2.1 Maxwell’s equations in the generalized coordinates

2.2 Bloch’s theorem

2.3 Discretization of the eigenvalue problems

3. Chern number calculation

3.1 Single-band Chern number

3.2 Composite Chern number

4. Numerical results

5. Conclusion

Funding

Acknowledgments

References

Cited By

Figures (6)

Tables (4)

Equations (27)

Optics Express