First-principle calculation of Chern number in gyrotropic photonic crystals

As an important figure of merit for characterizing the quantized collective behaviors of the wavefunction, Chern number is the topological invariant of quantum Hall insulators. Chern number also identifies the topological properties of the photonic topological insulators (PTIs), thus it is of crucial importance in PTI design. In this paper, we develop a first principle computatioal method for the Chern number of 2D gyrotropic photonic crystals (PCs), starting from the Maxwell's equations. Firstly, we solve the Hermitian generalized eigenvalue equation reformulated from the Maxwell's equations by using the full-wave finite-difference frequency-domain (FDFD) method. Then the Chern number is obtained by calculating the integral of Berry curvature over the first Brillouin zone. Numerical examples of both transverse-electric (TE) and transverse-magnetic (TM) modes are demonstrated, where convergent Chern numbers can be obtained using rather coarse grids, thus validating the efficiency and accuracy of the proposed method.


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][9][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][12][13][14].Compared with the commonly used plane-wave expansion (PWE) method [15][16][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.

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 Here, ∇ q is the partial differential operator in the generalized coordinate with three unit vectors û u u q (x, y, z), (q = 1, 2, 3).Ê and Ĥ are the normalized vector fields with Êi ).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.ε and μ are the relative permittivity and permeability in the form of 3 × 3 tensor which are written as where ε(r) and μ(r) are the original relative permittivity and permeability in the Cartesian basis.g is the metric tensor matrix which can be expressed as, With the metric tensor matrix, the length of a vector in the generalized coordinate is calculated by Fig. 1.The 2D Yee's grid in general coordinates.The dashed arrows are where the periodic boundaries apply.
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 Then, the Bloch's theorem can be written as

Discretization of the eigenvalue problems
From the Maxwell's curl equations (equation ( 1)), we can derive the governing equations for the modes of the PCs: 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, ¶ The grid points at the boundaries are treated by the Bloch's theorem, consequently where a 1 , a 2 are the lengths of the unit cell along the directions of û u u 1 and û u u 2 .By solving these eigenvalue equations, the band structure together with the eigenstates of 2D PCs can be easily obtained.The matrices are discretized from the operator ∇ × μ−1 (r)∇× and ∇ × ε−1 (r)∇×.The condition of the electromagnetic modes The subscript ω 1 and ω 1 stand for the different eigen-frequency under a fixed wave-vector k.

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 n th band can be computed by integrating the Berry curvature over the first Brillouin zone as following [20], Here, the is the Berry connection.u n,e,k denotes the normalized eigenstate satisfying |u n,e,k (ρ ρ ρ) = e ik•ρ ρ ρ E n,e,k (ρ ρ ρ) and |E n,e,k (ρ ρ ρ) is the eigenstate solved from equation (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 × 4 plaquettes.Then, the Chern number is calculated using where k ∆S k .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].

Composite Chern number
When the bands are degenerate, the Chern number can not be assigned to each band, so equation (20) will no longer be valid.For this case, these degenerate bands (taking the {n, n + 1, . .., n + N − 1} bands as an example) jointly share a composite (first) Chern number, which is associated with the multiplets By applying Stokes' theorem, the composite Chern number can be rewritten in terms of the composite Berry flux: The element (l, m) of the link matrix U Different from the single-band Chern number, in the composite Chern number, the link matrix U . In the discretized Brillouin zone, the Chern number is expressed as 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, u n,h,k α (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.

Numerical results
In the first numerical example, a 2D PC composed of Yttrium-Iron-Garnet (YIG) rods (ε = 15ε 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-ofplane) to break the time-reversal symmetry.The permeability is a tensor matrix written as where ω 0 = γH 0 is the precession frequency and ω m = 4πγM s .When the external magnetic field is 1600 Gauss at 4.28 GHz, κ = 12.4µ 0 , and µ = 14µ 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)

TM , C (3)
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 × 4. By further checking the values of the Chern numbers and composite Chern numbers, we find C (1⊕2) TM , and TM , which is as expected.

Mesh
Chern C (1) TM 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, Here, the parameters are chosen to satisfy the relationship Re- ferring to the YIG material under an external magnetic field at 4.28 GHz, the corresponding parameters are set as ε d = 14, ε f = −12.4and µ ⊥ = 15, µ d = 14, µ f = 12.4 and ε ⊥ = 15 [22].    2 and Table 3, 4 × 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.

Mesh
Chern C (1) TM 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 K point.By applying a DC magnetic field, the time-reversal symmetry is broken and the degeneracy is lifted, as shown in C TM C TM (a) The band structure of the TM modes.
TM/TE , C TM/TE denote the Chern numbers of the TM/TE modes of the associated bands.

Mesh
Chern 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 × 2) can guarantee the accuracy of the results.

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].

Fig. 2 .
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.
The band structure of the TM modes.
The band structure of the TE modes.

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

Fig. 4 .
Fig. 4. The Band structures of a 2D PC composed of gyrotopic rods in a square lattice when a 1600 Gauss +z DC magnetic field is applied.C (1) TM/TE , C

Fig. 5 .
Fig. 5.The unit cell and first Brillouin zone of the honeycomb lattice.
Band structure of the TM modes for a 2D honeycomb lattice.

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

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

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

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

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.