Universal numerical calculation method for the Berry curvature and Chern numbers of typical topological photonic crystals

Abstract

Chern number is one of the most important criteria by which the existence of a topological photonic state among various photonic crystals can be judged; however, few reports have presented a universal numerical calculation method to directly calculate the Chern numbers of different topological photonic crystals and have denoted the influence of different structural parameters. Herein, we demonstrate a direct and universal method based on the finite element method to calculate the Chern number of the typical topological photonic crystals by dividing the Brillouin zone into small zones, establishing new properties to obtain the discrete Chern number, and simultaneously drawing the Berry curvature of the first Brillouin zone. We also explore the manner in which the topological properties are influenced by the different structure types, air duty ratios, and rotating operations of the unit cells; meanwhile, we obtain large Chern numbers from −2 to 4. Furthermore, we can tune the topological phase change via different rotation operations of triangular dielectric pillars. This study provides a highly efficient and simple method for calculating the Chern numbers and plays a major role in the prediction of novel topological photonic states.

Key words： Chern number; topological photonic crystal; finite element method; symmetry

Introduction

In the previous decade, topological photonics has become a rapidly emerging research field in which geometrical and topological concepts are utilized to design and control the behavior of light. Topology was once a geometrical concept; however, it is currently being used in periodic reciprocal space [¹]. Currently, the topological photonic states are generally applied in case of photonic crystals because they have robust features and facilitate the unidirectional transmission of the light waves. The topological photonic states can exist at the interface of two crystals having different topological invariants, allowing the photons of the bandgap to be conducted with negligible dissipation. The topological photonic states can be applied in various photonic device applications owing to their unique properties, including the unidirectionally propagating edge states and robustness against impurities or defects without backscattering [²].

As a topological invariant of the two-dimensional (2D) dispersion band of the photonic crystal, the Chern number is defined based on the Berry curvature within the first Brillouin zone. It is one of the most important criteria used to judge the existence of a topological photonic state. In an electronic system, Thouless et al. showed that the Hall conductance of the 2D systems is proportional to the Chern number, which expresses an integer quantum Hall state on the edge of a topological nontrivial crystal cell, causing wave propagation along the edge [³]. However, until now, few reports have presented a direct method for calculating the Chern number. In majority of the cases, researchers seem to have avoided the calculation of Chern numbers because of their complexity; however, they have often identified the phase transition boundaries of the system via the topological phase diagram to identify the nonzero Chern numbers for evaluating the topological properties [⁴]. At times, the Chern numbers are directly given, including those for some particular topological photonic crystals, such as valley photonic crystals, whose Chern number is 0.5 of half the Brillouin zone for a special band [⁵–⁷]. The lack of a direct numerical calculation method for the Chern number restricts the development of topological photonics and hampers the prediction of novel topological photonic structures.

Herein, we propose a direct and universal numerical method for calculating the Chern number of typical photonic crystals based on the finite element method (FEM). Several difficulties that are associated with the numerical calculation process, including gauge fixing, rhombus Brillouin zone deformation, and discretization, are overcome. Furthermore, the influences on topological properties are studied with respect to different structure types, air duty ratios, and rotating operations of unit cells. Simultaneously, an effective methodology for tuning the topological phase change is developed by constructing different hetero-topological photonic crystals via rotation operations on dielectric triangular columns. This study not only provides a highly efficient and universal method for calculating the Chern number, but also exhibits potential applications in novel topological nanophotonic structures and the design of topological nanophotonic devices.

Calculation

Introduction to basic physical concepts

The topological invariants are defined to exhibit similar intrinsic properties between two isomorphic topological spaces. According to the Gauss-Bonnet theorem [⁴], the integral of the Gaussian curvature K over the whole manifold is a topological invariant, which can be defined in Eq. (1).

(1){\displaystyle \frac{1}{\mathrm{2\pi }}}{\displaystyle {\int }_{\mathrm{surface}}K\mathrm{d}A}=\mathrm{2}(\mathrm{1}-g).

The genus g is a typical topological invariant and denotes the number of holes on the surface. For example, coffee cups can become donuts by continuous deformation using the same genus g. They remain identical until continuous transformation occurs. The topology of a photonic system is defined on the dispersion band in the reciprocal (wave vector) space. The topological invariant of a 2D dispersion band is the Chern number, which characterizes the topological properties of the photonic bands [⁴]. Due to the periodicity, a torus is formed by the curling along the edge of a 2D lattice surface. Analogous to the Gaussian curvature, Berry curvature is defined [⁸], and the integral of the Berry curvature over the entire lattice surface is a constant called the Chern number, which is defined in Eq. (2).

(2)c_n={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }i}}{\displaystyle {\int }_{T^{\mathrm{2}}}{\mathrm{d}}^{\mathrm{2}}\mathit{k}{F}_{\mathrm{12}}(\mathit{k})},

where {{\displaystyle F}}_{\mathrm{12}}(\mathit{k}) is the Berry curvature defined in Eq. (3).

(3){{\displaystyle F}}_{\mathrm{12}}(\mathit{k})={{\displaystyle \partial }}_{\mathrm{1}}{{\displaystyle A}}_{\mathrm{2}}(\mathit{k})-{{\displaystyle \partial }}_{\mathrm{2}}{{\displaystyle A}}_1(\mathit{k}).

{{\displaystyle A}}_{\mu }(\mathit{k}) is the Berry connection [⁸,⁹] and is defined in Eq. (4).

(4){{\displaystyle A}}_{\mu }(\mathit{k})=\langle n(\mathit{k})\left|{{\displaystyle \partial }}_{\mu }\right|n(\mathit{k})\rangle ,

where |n(\mathit{k})\rangle denotes the nth normalized eigenstate of the photonic crystal using the Bloch wave vector k. For electromagnetic systems, the eigenstates are expressed using the six-component vector field u(\mathit{r}), which is defined in Eq. (5).

(5)u(\mathit{r})={({{\displaystyle E}}_x(\mathit{r}),{{\displaystyle E}}_y(\mathit{r}),{{\displaystyle E}}_z(\mathit{r}),{{\displaystyle H}}_x(\mathit{r}),{{\displaystyle H}}_y(\mathit{r}),{{\displaystyle H}}_z(\mathit{r}))}^{\mathrm{T}},

where {{\displaystyle E}}_x,{{\displaystyle E}}_y,{{\displaystyle E}}_z,{{\displaystyle H}}_x,{{\displaystyle H}}_y,{{\displaystyle H}}_z represent the components of the electromagnetic eigenmode normalized by their norms, which are defined by the inner product in Eq. (6) [¹].

(6)\langle u(\mathit{r}),v(\mathit{r})\rangle \equiv {\displaystyle \underset{\mathrm{unit}\hspace{0.17em}\mathrm{cell}}{\iint }\mathrm{d}{r}^{\mathrm{2}}{u}^{†}(\mathit{r})\cdot M(\mathit{r})\cdot v(\mathit{r})},

where M(\mathit{r})=\left[\begin{array}{ll}{{\displaystyle \epsilon }}_r(\mathit{r})& \mathrm{i}\xi (\mathit{r})\\ -\mathrm{i}\xi (\mathit{r})& {{\displaystyle \mu }}_r(\mathit{r})\end{array}\right] is the \mathrm{6}\times \mathrm{6} electromagnetic susceptibility tensor. Here, {{\displaystyle \epsilon }}_r(\mathit{r}) denotes the permittivity, {{\displaystyle \mu }}_r(\mathit{r}) denotes the permeability, \xi (\mathit{r}) describes the electromagnetic coupling [¹⁰], and u^{†}(\mathit{r}) denotes the Hermitian conjugate of u(\mathit{r}). In case of the all-dielectric structure where \xi =\mathrm{0}, the formulas can be simplified. From Eqs. (4) and (6), Berry connection can be expressed as Eq. (7) when \xi =\mathrm{0}.

(7){{\displaystyle A}}_{\mu }(\mathit{k})={\displaystyle \frac{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{E}^{†}(\mathit{r})\cdot {{\displaystyle \epsilon }}_r(r)\cdot {{\displaystyle \partial }}_{\mu }(\mathit{r})E(\mathit{r})+\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{H}^{†}(\mathit{r})\cdot {{\displaystyle \mu }}_r(\mathit{r})\cdot {{\displaystyle \partial }}_{\mu }H(\mathit{r})}{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{E}^{†}(\mathit{r})\cdot {{\displaystyle \epsilon }}_r(\mathit{r})\cdot E(\mathit{r})+\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{H}^{†}(\mathit{r})\cdot {{\displaystyle \mu }}_r(\mathit{r})\cdot H(\mathit{r})}}.

Furthermore, when \xi =\mathrm{0}, the electric and magnetic parts of the Berry connection exhibit the relation presented in Eq. (8) [¹⁰].

(8){\displaystyle \frac{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{E}^{†}(\mathit{r})\cdot {{\displaystyle \epsilon }}_r(\mathit{r})\cdot {{\displaystyle \partial }}_{\mu }E(\mathit{r})}{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{E}^{†}(\mathit{r})\cdot {{\displaystyle \epsilon }}_r(\mathit{r})\cdot E(\mathit{r})}}={\displaystyle \frac{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{H}^{†}(\mathit{r})\cdot {{\displaystyle \mu }}_r(\mathit{r})\cdot {{\displaystyle \partial }}_{\mu }H(\mathit{r})}{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{H}^{†}(\mathit{r})\cdot {{\displaystyle \mu }}_r(\mathit{r})\cdot H(\mathit{r})}}.

Therefore, the Berry curvature can be calculated using its magnetic part or its electric part (Eq. (9)).

(9){{\displaystyle A}}_{\mu }(\mathit{k})={\displaystyle \frac{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{E}^{†}(\mathit{r})\cdot {{\displaystyle \epsilon }}_r(\mathit{r})\cdot {{\displaystyle \partial }}_{\mu }E(\mathit{r})}{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{E}^{†}(\mathit{r})\cdot {{\displaystyle \epsilon }}_r(\mathit{r})\cdot E(\mathit{r})}}={\displaystyle \frac{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{H}^{†}(\mathit{r})\cdot {{\displaystyle \mu }}_r(\mathit{r})\cdot {{\displaystyle \partial }}_{\mu }H(\mathit{r})}{\iint {\mathrm{d}\mathit{r}}^{\mathrm{2}}{H}^{†}(\mathit{r})\cdot {{\displaystyle \mu }}_r(\mathit{r})\cdot H(\mathit{r})}}.

The topological properties of the photonic bandgap can be indicated using the gap Chern number {{\displaystyle C}}_{\mathrm{g}}, which can be defined as the sum of the Chern numbers of all the bands below the bandgap. The gap Chern numbers are useful to predict the existence of robust unidirectional edge states. In electronic systems, according to the bulk-edge correspondence [¹¹], the number of unidirectional edge states that emerge at the interface is equal to the difference between the gap Chern numbers of two photonic crystals when two electronic systems that have overlapped bandgaps with different gap Chern numbers interface with each other. In some photonic systems, such as the magneto–optical systems, the bulk-edge correspondence theory is applicable [¹²–¹⁵].

Calculations for the Chern number

In practical calculation, the electromagnetic field distribution in periodic structures can be conveniently obtained through FEM, where we used the COMSOL Multiphysics software. However, some questions have to be inevitably answered before calculating the Chern number. First, the Chern number expression contains an integral process. The data during the calculation were discrete, which is in conflict with the calculation of the Chern number. The general method is to use summation for replacing the integral. To obtain highly precise results, dense points should be selected even though it will result in a huge computational consumption. Second, the phase of the wave function was uncertain and denoted an arbitrary value after the Maxwell equations were solved using the COMSOL Multiphysics software. Therefore, a method was required to eliminate the effect of phase uncertainty. Third, for the partial derivative of a wave function in the Chern number expression, a function should be constructed to connect the wave function n and the corresponding wave vector \mathit{k} because they were discrete and independent computer data.

Here, we developed a direct and universal numerical method for calculating the discretized Chern number using COMSOL Multiphysics software, which is a commercial FEM software. FEM exhibits a high level of flexibility and reliable error control, enabling numerical calculations for considerably complex structures. The eigenstates were calculated in a discretized Brillouin zone using the numerical method. We discretized the Brillouin zone using straight lines parallel to the reciprocal lattice vectors (Fig. 1). When the Brillouin zone was discretized, the integrals and derivatives could be replaced by finite sums and differences for as long as the mesh of the discretized Brillouin zone was sufficiently fine. Therefore, the lattice Berry curvature could be calculated using Eq. (10).

(10){{\displaystyle F}}_{\mathrm{12}}\left({\mathit{k}}_l\right){{\displaystyle \delta \mathit{k}}}_{\mathrm{1}}{{\displaystyle \delta \mathit{k}}}_{\mathrm{2}}=({{\displaystyle \partial }}_{\mathrm{1}}{{\displaystyle A}}_2({\mathit{k}}_l)-{{\displaystyle \partial }}_{\mathrm{2}}{{\displaystyle A}}_{\mathrm{1}}({\mathit{k}}_l)){{\displaystyle \delta \mathit{k}}}_{\mathrm{1}}{{\displaystyle \delta \mathit{k}}}_{\mathrm{2}}\approx ({{\displaystyle \mathrm{\Delta }}}_{\mathrm{1}}{{\displaystyle A}}_{\mathrm{2}}){{\displaystyle \delta \mathit{k}}}_{\mathrm{2}}-({{\displaystyle \mathrm{\Delta }}}_{\mathrm{2}}{{\displaystyle A}}_{\mathrm{1}}){{\displaystyle \delta \mathit{k}}}_{\mathrm{1}}.

where {{\displaystyle \delta \mathit{k}}}_{\lambda }(\lambda =\mathrm{1},\mathrm{2}) denotes the length of the small displacement in the k space, {{\displaystyle \delta \mathit{k}}}_{\mathrm{1}}{{\displaystyle \delta \mathit{k}}}_{\mathrm{2}} denotes the area of the parallelogram formed by the vectors {{\displaystyle \delta \mathit{k}}}_{\mathrm{1}} and {{\displaystyle \delta \mathit{k}}}_{\mathrm{2}}, and ({{\displaystyle \mathrm{\Delta }}}_{\mu }{{\displaystyle A}}_{\lambda }){{\displaystyle \delta \mathit{k}}}_{\lambda } is defined in Eq. (11).

(11)({{\displaystyle \mathrm{\Delta }}}_{\mu }{{\displaystyle A}}_{\lambda }){{\displaystyle \delta \mathit{k}}}_{\lambda }={{\displaystyle A}}_{\lambda }({\mathit{k}}_l+{{\displaystyle \delta \mathit{k}}}_{\mu }){{\displaystyle \delta \mathit{k}}}_{\lambda }-{{\displaystyle A}}_{\lambda }({\mathit{k}}_l){{\displaystyle \delta \mathit{k}}}_{\lambda },\hspace{1em}(\lambda ,\mu =\mathrm{1},\mathrm{2}),

which is the finite difference of {{\displaystyle A}}_{\lambda } that can be represented by the finite difference of the eigenstate |n(\mathit{k})\rangle, which is defined in Eq. (12).

(12){{\displaystyle A}}_{\lambda }({\mathit{k}}_l){{\displaystyle \delta \mathit{k}}}_{\lambda }\approx \langle n({\mathit{k}}_l)|n({\mathit{k}}_l+{{\displaystyle \delta \mathit{k}}}_{\lambda })-n({\mathit{k}}_l)\rangle =\langle n({\mathit{k}}_l)|n({\mathit{k}}_l+{{\displaystyle \delta \mathit{k}}}_{\lambda })\rangle -\mathrm{1.}

Because the Berry curvature is defined in a discretized Brillouin zone, the Chern number can be calculated using Eq. (13).

(13){{\displaystyle c}}_n={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }i}}\mathrm{\sum }{{\displaystyle F}}_{\mathrm{12}}({\mathit{k}}_l){{\displaystyle \delta \mathit{k}}}_{\mathrm{1}}{{\displaystyle \delta \mathit{k}}}_{\mathrm{2}}.

Subsequently, we present the difficulties and challenges associated with the calculations as well as the solutions to overcome these problems. First, there was a problem called the U(1) gauge ambiguity of the Berry connection. Here, the problem lies in the fact that if the eigenstates are multiplied by a phase factor {\mathrm{e}}^{\mathrm{i}\chi [ik]} (U(1) gauge transformation), they will remain physically similar; however, the Berry connection A will transform into A+\mathrm{i}\nabla \chi (\mathit{k}). Therefore, for a similar physics scenario, there is more than one Berry connection distribution. This problem is referred to as the U(1) gauge ambiguity [¹]. As for the Berry curvature {{\displaystyle F}}_{\mathrm{12}}(\mathit{k}), it remains unchanged under the U(1) gauge transformation. However, it will also be influenced by multiple values of the eigenstates. In the definition of Berry curvature, the partial derivatives of the Bloch vector k require the eigenstates to be smooth functions of the Bloch vector k, but a global smooth function |n(\mathit{k})\rangle of k may not exist. To solve this problem, we obtained a locally smooth patch for each point in the Brillouin zone and calculated the Berry curvature on each of them, making a gauge transformation between two patches at their intersection [¹⁶].

Fig.1 Discretization of the Brillouin zone of the square lattices. Each node {{\displaystyle \mathit{k}}}_l represents a Bloch wave vector; {{\displaystyle \delta \mathit{k}}}_{\mathrm{1}} and {{\displaystyle \delta \mathit{k}}}_{\mathrm{2}} are the finite differences along each reciprocal lattice vector

Full size|PPT slide

In numerical calculations, we cannot ensure that the eigenstate calculated using the eigenvalue solver was a smooth function of the Bloch vector \mathit{k}. There are two ways to solve the problem. The first method is to fix a smooth branch by assigning the phase factor {\mathrm{e}}^{\mathrm{i}\chi [\mathit{k}]} to each eigenstate. Using this method, we could ensure that a component of the eigen electromagnetic field at a point in real space (similar to the x component of the eigen electric field at the origin {{\displaystyle E}}_x(\mathrm{0},\mathrm{0})) has a fixed phase angle (for example, to make them positive real numbers) for the entire Bloch vector \mathit{k} in an open set of momentum space [⁵]. If this process cannot continue, for example, {{\displaystyle E}}_x(\mathrm{0},\mathrm{0}) = 0 for some Bloch wave vector \mathit{k}, then another point in real space should be selected along with a new open set in the momentum space until the Brillouin zone is covered by all the open sets. The second method is to avoid the gauge variants [¹⁶]. In Ref. [¹²], the Berry curvature of each lattice was calculated using the “U(1) link variable” {{\displaystyle U}}_{\mu }{{\displaystyle \mathit{k}}}_i, which is expressed as Eq. (14).

(14){{\displaystyle F}}_{\mathrm{12}}({{\displaystyle \mathit{k}}}_l){{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mathrm{1}}{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mathrm{2}}=\ln \left({\displaystyle \frac{{{\displaystyle U}}_{\mathrm{1}}({{\displaystyle \mathit{k}}}_l){{\displaystyle U}}_{\mathrm{2}}({{\displaystyle \mathit{k}}}_l+{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mathrm{1}})}{{{\displaystyle U}}_{\mathrm{1}}({{\displaystyle \mathit{k}}}_l+{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mathrm{2}}){{\displaystyle U}}_{\mathrm{2}}({{\displaystyle \mathit{k}}}_l)}}\right),

where {{\displaystyle U}}_{\mu }({\mathit{k}}_l) is the U(1) link variable, which is defined in Eq. (15).

(15){{\displaystyle U}}_{\mu }({{\displaystyle \mathit{k}}}_l)\equiv {\displaystyle \frac{\langle n({{\displaystyle \mathit{k}}}_l)\left|n({{\displaystyle \mathit{k}}}_l+{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mu })\right|\rangle }{\left|\langle n({{\displaystyle \mathit{k}}}_l)|n({{\displaystyle \mathit{k}}}_l+{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mu })\rangle \right|}}\mathrm{.}

Further, from Eq. (16), we can observe that {{\displaystyle U}}_{\mu }({{\displaystyle \mathit{k}}}_l) is the normalization of the exponential of {{\displaystyle A}}_{\mu }({{\displaystyle \mathit{k}}}_l){{\displaystyle \delta {\displaystyle \mathit{k}}}}_u when {{\displaystyle \delta {\displaystyle \mathit{k}}}}_u is sufficiently small.

(16){\mathrm{e}}^{{{\displaystyle A}}_{\mathit{\mu }}({{\displaystyle \mathit{k}}}_l){{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mathit{\mu }}}={\mathrm{e}}^{\langle n({{\displaystyle \mathit{k}}}_l)|n({{\displaystyle \mathit{k}}}_l+{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\lambda })\rangle -\mathrm{1}}\approx \mathrm{1}+\langle n({{\displaystyle \mathit{k}}}_l)|n({{\displaystyle \mathit{k}}}_l+{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\lambda })\rangle -\mathrm{1}=\langle n({{\displaystyle \mathit{k}}}_l)|n({{\displaystyle \mathit{k}}}_l+{{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\lambda })\rangle .

Thus, Eq. (14) is equivalent to Eq. (10).

The second problem is associated with the traversal of the hexagonal Brillouin. The rhombus and honeycomb structures exhibit hexagonal first Brillouin zones; therefore, it is difficult to find a parallelogram mesh covering the whole Brillouin zone as in Fig. 1. Fortunately, because of the periodicity of the eigenstates, the hexagonal Brillouin could be deformed into a rhombus region by moving some part of it along the reciprocal vectors (Fig. 2(a)). Then, a uniform parallelogram mesh can be constructed over this rhombus Brillouin zone (Fig. 2(b)).

Another method known as the Wilson loop approach is available to calculate the Chern number. This method is discussed in Refs. [¹⁷,¹⁸] for electronic systems and in Ref. [¹⁰] for photonic systems. In this method, effective one-dimensional (1D) systems are constructed by fixing one component of the Bloch wave vector, and the Wannier centers of each effective 1D system are calculated by the Wilson loops corresponding to the fixed wave vector. The connection between the aforementioned integration method and the Wilson loop approach lies in the fact that the Chern number calculated by the integration of Berry curvature is equal to the winding number of the Wannier centers when the fixed wave vector traverses the 1D Brillouin zone.

We used the aforementioned gauge invariants to calculate the Chern number. The eigenstates were calculated using the commercial FEM software named COMSOL Multiphysics. The data extraction and processing were conducted using MATLAB. In our calculations, the eigenstates were normalized based on the inner product. Subsequently, the electric field data were extracted to calculate the inner product. Finally, the link variables and lattice Berry curvature were calculated using MATLAB. In Section 2.3, we show the numerical results of some typical topological photonic crystals.

Fig.2 Deformation and discretization of the Brillouin zone. (a) Deformation of the Brillouin zone; (b) discretization of the deformed Brillouin zone, where each node {{\displaystyle \mathit{k}}}_l represents the Bloch wave vector, {{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mathrm{1}} and {{\displaystyle \delta {\displaystyle \mathit{k}}}}_{\mathrm{2}} denote the finite differences along each reciprocal lattice vector

Full size|PPT slide

Calculation results for the typical topological photonic crystals

The nontrivial topological states are always attributed to broken symmetries. Two main types of symmetries, i.e., space inversion and time-reversal symmetries, are associated with the Berry curvature distribution. Time-reversal symmetries will lead to the relation {{\displaystyle F}}_{\mathrm{12}}(-{\displaystyle \mathit{k}})=-{{\displaystyle F}}_{\mathrm{12}}({\displaystyle \mathit{k}}), whereas inversion symmetries will lead to the relation {{\displaystyle F}}_{\mathrm{12}}(-{\displaystyle \mathit{k}})={{\displaystyle F}}_{\mathrm{12}}({\displaystyle \mathit{k}}) [¹]. A corollary of these relations is that the Chern number is zero when the time-reversal symmetries are preserved.

Therefore, there are two main methods to develop nontrivial topological structures. One is to break the time-reversal symmetries to obtain nonzero Chern numbers, similar to that observed in the Chern insulators in electronic systems. These structures can be fabricated using special materials such as the gyromagnetic materials [¹²]. The second method is to break the inversion symmetries by changing the shape or the dielectric constant. An important nontrivial topological structure is the valley topological structure, which has a Chern number of zero but nontrivial topological properties of different pseudospin modes. The commonly observed valley topological structures are honeycombs [¹⁹–²³] or rhombus lattices [¹⁶,²⁴,²⁵]. In this section, we will present some numerical results for the broken and preserved time-reversal symmetry structures. In our calculations, \mathrm{32}\times \mathrm{32} Bloch wave vectors are evenly distributed across the Brillouin zone. The relative tolerance of the eigenstate solver is {\mathrm{10}}^{-\mathrm{6}}.

As for the structures with broken time-reversal symmetries, we used common gyromagnetic materials that have an anisotropic permeability

\mathit{\mu }=\left[\begin{array}{lll}\mathrm{1}& \mathrm{i}\kappa & \mathrm{0}\\ -\mathrm{i}\kappa & \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],

where \kappa describes the gyromagnetic properties of the material. Large \kappa indicates a strong gyromagnetic effect in the material. First, we presented a structure of gyromagnetic cylinders in a square lattice (Fig. 3(a)). The radius of the cylinder was 0.13{{\displaystyle a}}_{\mathrm{0}}, where {{\displaystyle a}}_{\mathrm{0}} is the lattice constant. The relative permittivity of the gyromagnetic cylinders is \epsilon =\mathrm{13}, the diagonal elements of permeability are 1, and the off-diagonal element is \mathrm{i}\kappa =\mathrm{0.4}\mathrm{i}. The ambient medium was air, whose permittivity is \mathrm{1}. Figure 3(b) shows the dispersion band of the TM modes of the structure and the Chern number of each band. Figures 3(c)–3(f) show the Berry curvature distributions in the first Brillouin zone of the first four bands of the TM modes. The colormap denotes the value of the Berry curvature with the unit {{{\displaystyle a}}_{\mathrm{0}}}^{\mathrm{2}}/{\mathrm{4}\pi }^{\mathrm{2}}. Here, {{\displaystyle a}}_{\mathrm{0}} is the lattice constant. The red dashed line around the square marks the first Brillouin zone and the irreducible Brillouin zone. This photonic structure has a trivial bandgap between bands 1 and 2 ({{\displaystyle C}}_{\mathrm{g}}=\mathrm{0}), a nontrivial topological bandgap between bands 2 and 3 ({{\displaystyle C}}_{\mathrm{g}}=-\mathrm{1}), and a degenerate point between bands 3 and 4 at the X point. Because the structure preserves inversion symmetry and {{\displaystyle C}}_{\mathrm{4}v} symmetry, the Berry curvature distributions of the photonic bands show both inversion and {{\displaystyle C}}_{\mathrm{4}v} symmetries. Apart from symmetries, the Berry curvature denotes other details. For band 1 (Fig. 3(c)), the peaks can be observed at the X points. However, due to the negative value at the G and M points, the Chern number of band 1 is zero. For band 2 (Fig. 3(d)), negative valleys emerge at the M points, corresponding to the peak of the dispersion band at the M point in Fig. 3(b). For bands 3 and 4 (Figs. 3(e) and 3(f)), sharp peaks and valleys emerge at the X points. This is because bands 3 and 4 degenerate with each other at the X points and result in singular points at the X points in the Berry curvature distributions [¹]. It is interesting that although band 4 has four sharp negative valleys at the X points, it still has a positive Chern number because majority of the Brillouin zone is positive and the integration over the positive part exceeds the impact caused by the four valleys.

Fig.3 A gyromagnetic cylinder model denoting the dispersion relation of the TM mode, band Chern numbers, and Berry curvature. (a) Geometry of a unit cell. The blue part represents the gyromagnetic cylinder with r = 0.13{{\displaystyle a}}_{\mathrm{0}}, e = 13, and κ = 0.4. The gray part symbolizes pure air; (b) first four bands of the TM modes. The Chern numbers are shown at the bottom, whereas the gap Chern numbers are marked at the bandgaps; (c)–(f) Calculated Berry curvature within the first Brillouin zone of the first four bands of the TM modes: (c) for band 1; (d) for band 2; (e) for band 3; and (f) for band 4. The square red dashed line represents the first Brillouin zone, and the triangular red dashed line represents the irreducible Brillouin zone

Full size|PPT slide

Subsequently, we changed the gyromagnetic cylinders from circular columns to square columns using the same material (Fig. 4(a)). The side length of the square is l = 0.26{{\displaystyle a}}_{\mathrm{0}}. The TM modes were investigated. Figure 4(b) shows the dispersion relation, Chern numbers, and gap Chern numbers of the first four TM bands, and the Berry curvature distributions are shown in Figs. 4(c)–4(f). From Fig. 4(b), we can see that the Chern number of each band remains unchanged when compared with that of the structure in Fig. 3 but that the degeneracy at the X point between bands 3 and 4 is lifted.

Fig.4 A gyromagnetic square column model denoting the dispersion relation of the TM mode, band Chern numbers, and Berry curvature. (a) Geometry of a unit cell. The blue part represents the gyromagnetic square column when l = 0.26{{\displaystyle a}}_{\mathrm{0}}, e = 13, and k = 0.4. The gray part symbolizes pure air; (b) first four bands of the TM modes. The Chern numbers are shown at the bottom, whereas the gap Chern numbers are marked at the bandgaps; (c)–(f) Calculated Berry curvature within the first Brillouin zone of the first four bands of the TM modes: (c) for band 1; (d) for band 2; (e) for band 3; and (f) for band 4. The red dashed line represents the first Brillouin zone, and the triangular red dashed line represents the irreducible Brillouin zone

Full size|PPT slide

Therefore, in the Berry curvature distributions, the singular points at the X points of bands 3 and 4 become the broadened Berry curvature extrema. The singular points of band 3 become peaks that are not considerably sharp, and the peaks at the M points become more evident because the Berry curvature at the X points drastically decreases. The singular points of band 4 become negative valleys surrounded by two positive peaks, corresponding to the frequency extrema of the dispersion band (Fig. 4(b), orange dotted line). However, the Berry curvature distributions of bands 1 and 2 are less influenced by the structural change. As for symmetries, the Berry curvature distributions show identical symmetries because the structure does not break either the inversion symmetry or the {{\displaystyle C}}_{\mathrm{4}v} symmetry.

Next, we changed the square lattice into a honeycomb lattice with the same material (Fig. 5). The structure comprised gyromagnetic cylinders with k = 0.4, e = 13, and radius r = 0.13{{\displaystyle a}}_{\mathrm{0}}, where {{\displaystyle a}}_{\mathrm{0}} is the lattice constant of the honeycomb lattice, i.e., \sqrt{\mathrm{3}} times the side length of the hexagonal unit cell. The dispersion relations, Chern numbers, gap Chern numbers (Fig. 5(b)), and Berry curvature distributions of the first four bands of the TM modes (Figs. 5(c)–5(f)) are shown in the figures. The first Brillouin zones were plotted using the red dashed lines to clearly illustrate the special points (G, M, K) in the k space. From Fig. 5(b), we can observe that this structure has a trivial bandgap between bands 1 and 2, a {{\displaystyle C}}_{\mathrm{g}}=-\mathrm{1} bandgap between bands 2 and 3, and a {{\displaystyle C}}_{\mathrm{g}}=\mathrm{1} bandgap between bands 3 and 4. From the Berry curvature distribution of band 1 (Fig. 5(c)), positive peaks and negative valleys can be observed at the M and K points, respectively, causing the Chern number to become zero. From the Berry curvature distribution of band 2 (Fig. 5(d)), six significant valleys emerge at the K points, corresponding to the frequency extremum of its dispersion band (Fig. 5(b), red dash–dot line). From the Berry curvature of band 3 (Fig. 5(e)), the Berry curvature distribution shows positive peaks at the M and K points, corresponding to the extrema of the dispersion band (Fig. 5(b), orange dotted line). Figure 5(f) illustrates the Berry curvature of band 4, which shows gentle positive peaks at the K points and relatively sharp negative valleys at the M points. Although the positive peaks are not as evident as the negative peaks, the Chern number of band 4 is still positive.

Fig.5 Gyromagnetic cylinder model denoting the dispersion relation of the TM mode, band Chern numbers, and Berry curvature. (a) Hexagon cell of a single gyromagnetic cylinder. The blue part represents the gyromagnetic cylinder when r = 0.13{{\displaystyle a}}_{\mathrm{0}}, e = 13, and k = 0.4. The gray part symbolizes pure air; (b) first four bands of the TM modes. The Chern numbers are shown at the bottom, whereas the gap Chern numbers are marked at the bandgaps; (c)–(f) calculated Berry curvature within the first Brillouin zone of the first four bands of the TM modes: (c) for band 1; (d) for band 2; (e) for band 3; and (f) for band 4. The red dashed line represents the first Brillouin zone, and the triangular red dashed line represents the irreducible Brillouin zone

Full size|PPT slide

As for symmetries, this structure preserves both inversion symmetry and {{\displaystyle C}}_{\mathrm{6}v} symmetry. The Berry curvature distributions preserve these symmetries as well.

Fig.6 Dielectric triangular column model denoting the dispersion relation of the TM mode, band Chern numbers, and Berry curvature. (a) Hexagon cell of a single air triangular column. The gray part represents the air triangular column with a side length of l={{\displaystyle \mathrm{0.808}a}}_{\mathrm{0}}. The blue part symbolizes a dielectric slab when e = 8.9 and m = 1; (b) dispersion relation of the first four bands of the TM modes. The Chern number of each band is zero; (c)–(f) calculated Berry curvature within the first Brillouin zone of the first four bands of the TM modes: (c) for band 1; (d) for band 2; (e) for band 3; and (f) for band 4. The red dashed line represents the first Brillouin zone

Full size|PPT slide

As shown above, the structures with broken time-reversal symmetries always have nonzero chern numbers. however, this does not indicate that the berry curvature is meaningless for the all-dielectric structures that preserve the time-reversal symmetries. an important type of all-dielectric photonic structure is the valley photonic structure, which has a pair of frequency extrema for the dispersion bands at the high-symmetry points (K/K′) and exhibits opposite berry curvatures at the K and K′ points [⁵]. The valley photonic crystals support the robust edge states but the direction of each edge state relies on the chiral mode of the source. we will present an all-dielectric model in this section with further discussion being presented in Section 3.3.

The all-dielectric model comprises triangular holes in a dielectric slab with a permittivity of e = 8.9 and a permeability of m = 1 (Fig. 6(a)). The dielectric holes are aligned in a honeycomb lattice. The side length of the triangle is l={{\displaystyle \mathrm{0.808}a}}_{\mathrm{0}}. The TM modes of this structure were investigated. The dispersion relation is shown in Fig. 6(b), and the Berry curvature distributions are shown in Figs. 6(c)–6(f). The Chern numbers of all the bands are zero. In this structure, inversion symmetry is no longer preserved; therefore, the high-symmetry points K and K′ are no longer equivalent, and the irreducible Brillouin zone is no longer the one bounded by polyline G–M–K–G but the one bounded by G–K–K′–G. However, as far as the dispersion relation is concerned, the K and K′ points show the same eigenfrequencies. Therefore, in Fig. 6(b), we choose the polyline G–M–K–G. From Fig. 6(b), we can observe a degeneracy between bands 1 and 2 at the K point, a degeneracy between bands 3 and 4 at the G point, and a degeneracy between bands 4 and 5 at the K point. These degenerate points correlate with the singular points in the Berry curvature distributions.

In Figs. 6(c) and 6(d), there are singular points at the K and K′ points; for each pair of K and K′ points, the Berry curvature has the opposite sign. Furthermore, the integrations of the Berry curvature for bands 1 or 2 over the two halves of the Brillouin zone containing K and K′ are 0.5 and −0.5, respectively. In Fig. 6(e), there are singular points near the G point, corresponding to the degeneracies at the G point in the dispersion band (Fig. 6(b)). In Fig. 6(f), the singular points emerge at the G and K/K′ points because the fourth band degenerates with the fifth band at the K/K′ point and with the third band at the G point. As for symmetries, because of the preservation of the time-reversal symmetry, the Berry curvature shows an antisymmetric distribution {{\displaystyle F}}_{\mathrm{12}}(-\mathit{k})=-{{\displaystyle F}}_{\mathrm{12}}(\mathit{k}).

Based on the aforementioned numerical results, we can conclude that the symmetries of the structure and the degeneracy of the bands affect the Berry curvature distributions. First, when the structures have some types of symmetries, their Berry curvature distributions will exhibit corresponding symmetries. The two most significant symmetries are the time-reversal symmetries and inversion symmetries, which will lead to antisymmetry {{\displaystyle F}}_{\mathrm{12}}(-\mathit{k})=-{{\displaystyle F}}_{\mathrm{12}}(\mathit{k})and inversion symmetry {{\displaystyle F}}_{\mathrm{12}}(-\mathit{k})={{\displaystyle F}}_{\mathrm{12}}(\mathit{k})in Berry curvature distributions, respectively, apart from some critical points. Second, we illustrate that the singular points in the Berry curvature distributions are always correlated with the degeneracies in the dispersion bands. When a small perturbation is added to the structure and the degeneracies are lifted, the singular points in the Berry curvature distributions will transform into less sharp Berry curvature extrema.

Discussion

Different structure types

In the previous section, we presented our calculation results for the all-dielectric photonic crystals as well as the gyromagnetic photonic crystals and verified the validity of our calculation method. Apart from these aforementioned structures, our method is applicable to all types of photonic crystals, including those with anisotropic permittivity and permeability. It is also applicable to photonic crystals with different lattices, including rhombus lattices. In this section, we explore the manner in which the Berry curvature is influenced by the changes in a photonic crystal lattice.

Fig.7 Calculation results of a hexagon cell with three geometric shapes but with the same air duty ratio of 0.192. (a) Hexagon cell with a gyromagnetic cylinder, gyromagnetic square column, and gyromagnetic triangle column at its center; (b) first four band structures of the TM mode and the Chern number of the three structures presented in (a); (c) calculated band structures of the TM mode and Chern number from the fifth to the eighth band of a hexagon cell with a gyromagnetic cylinder, gyromagnetic square column, and gyromagnetic triangle column at its center

Full size|PPT slide

Thus, we transformed the gyromagnetic column in a hexagon lattice to denote cylindrical, square, and triangular shapes, while ensuring an identical filling ratio to the details presented in the upper, middle, and lower parts of Fig. 7(a), respectively. The photonic energy bands of the TM mode for three types of photonic crystals are drawn along with various lines in Figs. 7(b) and 7(c). Figure 7(b) shows the first four bands, whereas Fig. 7(c) shows bands 5–8. This shows that the first three energy bands exhibit very similar line shapes for three different structures and that differences can be observed from the fourth band. If rotation causes the photonic system to become invariant, the rotated mode itself becomes an allowed mode with a similar frequency. It is obvious that the columns in a hexagon lattice exhibiting cylinder, square, and triangle configurations exhibit different rotation symmetries; therefore, their photonic energy bands will be different [²⁶]. In the low-frequency domain, the energy is mainly concentrated in the gyromagnetic column regions instead of the air regions for the TM modes; further, all three columns have the same filling ratio in the same hexagon lattice. Therefore, the energy bands, particularly the first band, denote considerably similar line shapes for the low-frequency domain, where they are almost identical (see Fig. 7(b)). With an increase in frequency, the energy distributions are observed to be increasingly influenced by different structures and the energy bands exhibit increasingly different line shapes, especially from bands 4 to 8, as can be observed from Figs. 7(b) and 7(c). The energy bands are considerably different for different column shapes from a normalized frequency from 0.57 in the fourth band, and the Chern numbers demonstrate large differences from the fifth band. In the high-frequency domain (such as from bands 5 to 8), the interactions can be enhanced and more modes are involved owing to the relatively dense distributions of the energy bands; therefore, different structures will generate different Chern numbers in the high-frequency domain [¹³].

Different air duty ratio

Further, we enlarged the gyromagnetic cylinder in the same hexagon lattice by increasing the radius. The calculated Chern numbers and photonic energy bands are shown in Fig. 8.

Fig.8 Gyromagnetic cylinders with different radii in the center of a hexagon cell, maintaining the hexagon unchanged. (a) Band structures of the TM mode and the corresponding Chern numbers of a gyromagnetic cylinder with a radius of 0.13{{\displaystyle a}}_{\mathrm{0}}; (b) band structures of the TM mode and corresponding Chern numbers of a gyromagnetic cylinder with a radius of 0.23{{\displaystyle a}}_{\mathrm{0}}; (c) band structures of the TM mode and the corresponding Chern numbers of a gyromagnetic cylinder with a radius of 0.33{{\displaystyle a}}_{\mathrm{0}}

Full size|PPT slide

According to the calculated results, drastic changes can be observed with respect to the topological property when the radius is changed. As we can see, with an increase in the radius of the gyromagnetic cylinders, the frequency of each band will decrease in different structures. The normalized frequency of the fourth band is approximately 0.8 when r = 0.13{{\displaystyle a}}_{\mathrm{0}}, as shown in Fig. 8(a), decreases to approximately 0.65 when r = 0.23{{\displaystyle a}}_{\mathrm{0}}, as shown in Fig. 8(b), and finally decreases to become less than 0.5 when r = 0.33{{\displaystyle a}}_{\mathrm{0}}, as shown in Fig. 8(c). The third band is initially between 0.6 and 0.8 for r = 0.13{{\displaystyle a}}_{\mathrm{0}} and becomes less than 0.4 for r = 0.33{{\displaystyle a}}_{\mathrm{0}}. The second band’s behavior is also similar. The ground state beneath the bandgap decreases. The ground state’s decrease is small, which reduces the bandgap. This phenomenon can be explained based on the electromagnetic variational theorem mentioned in Ref. [²⁶], which indicates that low-frequency modes concentrate their energy in high-e regions and that high-frequency modes exhibit a larger energy concentration in low-e regions. Thus, when we increase the radius of a central cylinder that has higher e than that of the surrounding air, more energy from each mode is pushed into the cylinder and the frequency of each corresponding mode naturally decreases. Additionally, along with the evolution of the parameter, the ground state maintains its shape and the other bands above the gap randomly change their appearance. We attribute this to the various energies of different bands. The ground state has the least energy; therefore, it is more likely to maintain a constant energy distribution in case of a structural change, whereas the opposite is true in case of higher bands. Furthermore, the Chern number of the ground state, i.e., the first band, maintains a value of zero for different air duty ratios; however, the states above the bandgap behave considerably differently. The Chern numbers of the higher order bands change in a nonlinear manner. For example, in band 4, the Chern number c= 1 when r = 0.13{{\displaystyle a}}_{\mathrm{0}} (Fig. 8(a)) becomes −2 when r = 0.23{{\displaystyle a}}_{\mathrm{0}} (Fig. 8(b)) and finally zero when r = 0.33{{\displaystyle a}}_{\mathrm{0}} (Fig. 8(c)). Similar to Fig. 7(c), where different energy bands correspond to different Chern numbers, this can be explained by the theorem about large Chern numbers provided in Refs. [¹³,²⁷], which illustrates that the Chern number is directly influenced by the band structure.

Different rotation operations

Here, we explore the manner in which the band structures and Berry curvature change in case of different rotation operations (Fig. 9).

Fig.9 Band structures of the TE mode and Berry curvature of the hexagon cell after the anticlockwise rotation operation. (a) Band structure of the TE mode of the hexagon cell in which the inside triangle is in the state of origin location; (b) band structure of the TE mode of the hexagon cell in which the inside triangle is rotated 15° anticlockwise; (c) band structure of the TE mode of the hexagon cell in which the inside triangle is rotated 30° anticlockwise; (d) calculated Berry curvature of the second band corresponding to the structures without rotation; (e)–(f) calculated Berry curvature of the second band corresponding to the structures that are rotated by 15° and 30°, respectively: (e) for the hexagon cell in which the inside triangle is rotated 15° anticlockwise and (f) for the hexagon cell in which the inside triangle is rotated 30° anticlockwise

Full size|PPT slide

The structure is formed by the triangular columns in a honeycomb lattice. The side length of the triangle is l={{\displaystyle \mathrm{0.808}a}}_{\mathrm{0}}. The band structures and Berry curvature that have been discussed in other sections are in the TM mode. To better understand the distinct change in band structures and the energy-density distribution of the unidirectional propagation edge state after the rotation operations, we calculated the band structures and Berry curvature of a hexagon cell with a silicon triangular column inside the TE mode because majority of the energy is concentrated in the gyromagnetic columns in the TM mode. The energy-density distribution barely changed because the shape of the silicon triangular column did not alter after the rotation operation. There was no obvious difference between the band structures and the energy-density distribution of the unidirectional propagation edge state in the TM mode. However, in the TE mode, the energy was mainly located in the air column. The energy-density distribution drastically changed as the shape of the air column altered considerably after the rotation operation. The band structure and energy-density distribution of the unidirectional propagation edge state were considerably different after the rotation operation in the TE mode.

To investigate the manner in which the band structures changed in the TE modes, we rotated the silicon triangular column in the hexagon cell by 15° anticlockwise each time. The band diagram of a hexagon cell with a silicon triangular column inside without considering any rotation is shown in Fig. 9(a). Two degenerate points are shown. Next, we rotated the silicon triangular column in the hexagon cell by 15°, the band diagram of which is shown in Fig. 9(b). Two degenerate points were obtained as we broke the parity symmetry [²⁸]. When the rotation angle was increased to 30°, the size of the bandgap became larger as shown in Fig. 9(c) because we deepened the destruction of symmetry.

In addition, we also explored the manner in which the Berry curvature changed with varying degree of rotation. The calculated Berry curvature for the second band of the silicon triangular column in the hexagon cell without rotation is shown in Fig. 9(d). There are six peaks in the first Brillouin zone that are not exactly identical because the Berry curvature has a {{\displaystyle C}}_{\mathrm{3}v} symmetry, with one negative peak and positive peak being observed at the K and K′ points, respectively [²⁸]. Interestingly, Fig. 9(e) reveals the Berry curvature pattern with 18 nodes in the first Brillouin zone as a silicon triangular column in the hexagon cell that was rotated 15° anticlockwise. The {{\displaystyle C}}_{\mathrm{3}v} symmetry shown in Fig. 9(e) did not change compared with that in Fig. 9(d). However, the negative peak at the K point in the example without rotation shown in Fig. 9(d) became a positive peak in case of a rotation of 15° shown in Fig. 9(e) by simply increasing the rotation angle from 0° to 15°. The Berry curvature pattern when the rotation angle was increased to 30° is shown in Fig. 9(f) with 18 nodes appearing in the first Brillouin zone. The {{\displaystyle C}}_{\mathrm{3}v} symmetry shown in Fig. 9(f) did not change compared with those in Figs. 9(d) and 9(e), and there was no swap of the K and K′ points (i.e., the negative peak is at the K point) shown in Fig. 9(f) in case of a rotation of 30° when compared with the K and K′ points shown in Fig. 9(d) in case of no rotation. This shows that the rotation angle determines whether the K point is a negative or positive peak. The signs of the peaks in the K points will remain identical when the rotation angle is an integer multiple of 30°.

Berry curvature and boundary states

Berry curvature can provide us important information that cannot be obtained from the dispersion bands. Consider the rotating triangular column model presented in Section 3.3 as an example. Consider three models with anticlockwise rotation angles of −15°, 15°, and 45° (Fig. 10). We calculated their dispersion bands (Figs. 10(a)–10(c)) and the Berry curvature distributions of their second bands (Figs. 10(d)–10(f)). From Figs. 10(a)–10(c), we can see that the dispersion relations of the three models are similar, whereas the Berry curvature distributions of their second bands (Figs. 10(d)–10(f)) are different. The Berry curvature distributions of the models rotated by 15° and 45° are identical, whereas the Berry curvatures of the models rotated by −15° and 15° are opposite to each other. The differences in their Berry curvatures can be attributed to different mirror transformations. Consider the mirror transformation along the x-axis as an example. For the 2D TE modes, the Maxwell equation can be simplified as presented in Eq. (17).

(17){{\displaystyle \nabla }}_{\mathrm{r}}^{\mathrm{2}}{{\displaystyle H}}_z+\epsilon {\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\mathrm{1}}{\epsilon }}\right){\displaystyle \frac{\partial }{\partial x}}{{\displaystyle H}}_z+\epsilon {\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{\mathrm{1}}{\epsilon }}\right){\displaystyle \frac{\partial }{\partial y}}{{\displaystyle H}}_z=-\epsilon {{\displaystyle \left({\displaystyle \frac{\omega }{c}}\right)}}^{\mathrm{2}}{{\displaystyle H}}_z,

where {{\displaystyle H}}_z is the z component of magnetic intensity and \epsilon denotes the permittivity. Substituting the coordinate (x,-y) for (x,y)and using the chain rule for partial derivatives, Eq. (17) becomes Eq. (18).

(18){\displaystyle \frac{\mathrm{1}}{\epsilon (x,-y)}}{{\displaystyle \nabla }}_{\mathrm{r}}^{\mathrm{2}}{{\displaystyle H}}_z(x,-y)+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\mathrm{1}}{\epsilon (x,-y)}}\right){\displaystyle \frac{\partial }{\partial x}}{{\displaystyle H}}_z(x,y)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{\mathrm{1}}{\epsilon (x,-y)}}\right){\displaystyle \frac{\partial }{\partial y}}{{\displaystyle H}}_z(x,-y)=-{{\displaystyle \left({\displaystyle \frac{\omega }{c}}\right)}}^{\mathrm{2}}{{\displaystyle H}}_z(x,-y)\mathrm{.}

This means that {{\displaystyle H}}_z(x-y) satisfies the Maxwell equation if we substitute \epsilon (x,y) with \epsilon (x,-y). Therefore, if the geometry is transformed under a mirror transformation, the eigenstates of the new structure will transform under the same transformation, i.e., they will become a mirror of the original eigenstates about the x-axis. However, under this transformation, the parity changes and the Berry curvature takes the opposite sign. Therefore, the result is that under mirror transformation, the Berry curvature distribution is a mirror of the original distribution that takes the opposite sign. With respect to our example, for constructing a structure with an angle of −15°, the structure with an angle of 15° can be used as the mirror along the x-axis, which will transport the K(K′) points to the K(K′) points and change the sign of the Berry curvature. Due to the {{\displaystyle C}}_{\mathrm{3}} symmetry of the Berry curvature of the {{\displaystyle C}}_{\mathrm{3}} symmetric structures, different K(K′) points exhibit the same Berry curvature. Thus, this mirror transformation will eventually change the sign of the Berry curvature. However, for constructing a structure with an angle of 45°, the structure with an angle of 15° along the y-axis can be mirrored, which will transport the K(K′) points to the K′(K) points and change the sign of the Berry curvature. Due to the antisymmetry of the Berry curvature distribution of the time-reversal symmetric structures, the sign of the Berry curvature at the K(K′) points will revert to the original value.

Fig.10 Band structures of the TE mode and Berry curvature of the hexagon cell after various rotation operations. (a) Band structure of the TE mode of the hexagon cell in which the inside triangle is rotated 15° clockwise; (b) band structure of the TE mode of the hexagon cell in which the inside triangle is rotated 15° anticlockwise; (c) band structure of the TE mode of the hexagon cell in which the inside triangle is rotated 45° anticlockwise; (d)–(f) calculated Berry curvature corresponding to the structures rotated by −15°, 15°, and 45°, respectively: (d) for the hexagon cell in which the inside triangle is rotated 15° clockwise; (e) for the hexagon cell in which the inside triangle is rotated 15° anticlockwise; and (f) for the hexagon cell in which the inside triangle is rotated 45° anticlockwise

Full size|PPT slide

To denote the manner in which the differences between the Berry curvature distributions cause different properties in photonic transport, we constructed two models by joining the photonic crystals rotated by 15° and −15° (Fig. 11(a)) and those rotated by 15° and 45° (Fig. 11(e)). We calculated the dispersion bands of the two models in Figs. 11(b) and 11(f) and marked the bulk states and edge states using blue and red lines, respectively. In Fig. 11(b), there are edge states whose dispersion band lies within the bandgap of the bulk state (red line). From Fig. 11(b), the right side of the dispersion band of the edge states has a positive slope, indicating that these edge states have positive group velocities and propagate toward the positive direction of k. Similarly, the edge states of the left side exhibit negative group velocities and propagate toward the negative direction of k. However, there is no such edge state in Fig. 11(f). Next, we constructed structures in which the two photonic crystals were bounded by a zigzag interface, shown in Figs. 11(c) and 11(d), to show the chiral and robust properties of the edge states. We placed a source in the middle of the interface in each structure. The source comprised six vertical magnetic currents oscillating at a frequency of \mathrm{0.45}c/{{\displaystyle a}}_{\mathrm{0}}, where {{\displaystyle a}}_{\mathrm{0}} is the lattice constant and c is the speed of light. The magnetic currents were aligned to the vertices of a hexagon with a side length of {{\displaystyle \mathrm{0.1}a}}_{\mathrm{0}}/\sqrt{\mathrm{3}}, and each pair of adjacent currents exhibited the same phase difference. For a model comprising photonic crystals with triangles rotated by 15° and -15°, we excited the edge states of different chirality in Figs. 11(c) and 11(d), respectively.

Fig.11 Edge state in the photonic-crystal-based topological insulator. (a) First type of structure whose unit cell comprises four hexagon cells with a black line. The unit cell with the black line is periodic along the x direction. The triangles inside the hexagon cells that were located above the red line were rotated by 15° clockwise. The triangles inside the hexagon cells that were located below the red line were rotated by 15° anticlockwise. (b) Band diagram of the TE mode of the periodic 2D structure located in (a) shows the edge state, which is denoted using a red line. The blue lines represent the bulk state, and the green rectangle represents the bandgap. (c) Energy-density distribution of the edge state propagating unidirectionally toward the left, corresponding to the structure shown in (a). The six yellow points at the top-right corner represent the point source comprising six vertical magnetic currents. The arrow indicates the direction in which the phase increased. The phases of the currents shown in (c) increased counterclockwise. (d) Energy-density distribution of the edge state that propagated unidirectionally toward the right, corresponding to the structure shown in (a). The source is identical with (c) apart from the phases of the currents shown in (d) that increased clockwise. (e) Second type of structure whose unit cell comprised four hexagon cells with a black line. The unit cell with the black line is periodic along the x direction. The triangles inside the hexagon cells that were located above the red line were rotated by 45° anticlockwise. The triangles inside the hexagon cells that were located below the red line were rotated by 15° anticlockwise. (f) Band diagram of the TE mode of the periodic 2D structure located in (b) shows no edge state. The blue lines represent the bulk state, and the green rectangle represents the bandgap. (g) Energy-density distribution for structures shown in (e). The source is the same as (c), where the phase of the currents increased counterclockwise. No edge state existed that propagated unidirectionally. (h) Energy-density distribution for the structure shown in (e). The source is the same as (d), where the phase of currents increased clockwise. No unidirectionally propagating edge state existed

Full size|PPT slide

As shown in Fig. 11(c), the phases of the currents increased anticlockwise and the edge states propagated unidirectionally toward the left, corresponding to the edge states with negative velocities. In Fig. 11(d), the phases of the currents increased clockwise and the unidirectional edge state emerged but propagated in the opposite direction, corresponding to the edge states with positive group velocities. This phenomenon shows that the edge states with different group velocities exhibit different chiral properties. Therefore, the edge states with different chiralities can travel unidirectionally along different directions. Moreover, the edge states were robust against the sharp twisting angles. Figures 11(c) and 11(d) illustrate that little electromagnetic wave is leaked or scattered at the twisting angles. The reason is that the bandgap of the bulk state prevents the leaking of the wave and that the unidirectionality of the edge state restrains the scattering. For a model comprising photonic crystal with triangles rotated by 15 and 45° (Fig. 11(e)), neither the source shown in Fig. 11(c) nor that shown in Fig. 11(d) could excite any of the edge states because from Fig. 11(f), we can observe that this structure had no edge states. The difference between the two modes lies in the Berry curvature of the bulk states. For the structure shown in Fig. 11(a), the bulk states of the two sides exhibited opposite Berry curvatures, as shown in Figs. 10(d) and 10(e). Therefore, similar to the quantum valley Hall effect [²³], edge states emerged at the interface of the two photonic crystals. However, the bulk states of the structure shown in Fig. 11(e) exhibited the same Berry curvature, as shown in Figs. 10(e) and 10(f). Therefore, there was no edge state at the interface. Thus, we have proved our algorithm is suitable for the calculation of the Berry curvature and Chern number, which can enable us to obtain the accurate topological invariant of some typical photonic crystals.

Conclusions

Herein, we proposed a universal method for calculating the Berry curvature and Chern numbers. Our method was based on FEM, exhibiting considerable flexibility and controlled accuracy. Using this method, we calculated the Chern numbers and Berry curvature distributions of different types of trivial and topological photonic crystals. To denote the flexibility of our numerical method, we presented the results of the photonic crystals using gyromagnetic materials and all-dielectric materials. Furthermore, using different structures, we analyzed the manner in which the Berry curvature and Chern numbers are correlated with the geometrical structures and dispersion relations. Next, we discussed the impact of the different shapes of columns, air duty ratios, and rotation angles on the Berry curvature and Chern numbers of the photonic crystals and illustrated the interesting phenomena arising from the novel topological states, including robust unidirectional edge states. In conclusion, our method is very useful to search novel topological structures and find intriguing properties of considerably complicated structures. Further, our method can be easily combined with the intelligent computation and optimization methods, making it helpful for automatically searching topological structures.

Acknowledgements

This paper was supported by the National Natural Science Foundation of China (Grant Nos. 11604378, 91850117, and 11654003), Beijing Institute of Technology Research Fund Program for Young Scholars, and Double First Class University Plan. We would like to thank Prof. Xiangdong Zhang, Dr. Lu He, Dr. Yujing Wang, and Dr. Changyin Ji from the Beijing Institute of Technology for the useful discussion.

References

1 Raghu S, Haldane F D M. Analogs of quantum-Hall-effect edge states in photonic crystals. Physical Review A, 2008, 78(3): 033834

2 Wu Y, Li C, Hu X, Ao Y, Zhao Y, Gong Q. Applications of topological photonics in integrated photonic devices. Advanced Optical Materials, 2017, 5(18): 1700357

3 Den N, Quantized M. Hall conductance in a two dimensional periodic potential. Physica A, 1984, 124(1): 199–210

4 Lu L, Joannopoulos J D, Soljačić M. Topological photonics. Nature Photonics, 2014, 8(11): 821–829

5 Chen X D, Zhao F L, Chen M, Dong J W. Valley-contrasting physics in all-dielectric photonic crystals: Orbital angular momentum and topological propagation. Physical Review B, 2017, 96(2): 020202

6 Gao Z, Yang Z J, Gao F, Xue H R, Yang Y H, Dong J W, Zhang B L. Valley surface-wave photonic crystal and its bulk/edge transport. Physical Review B, 2017, 96(20): 201402

7 Kang Y, Ni X, Cheng X, Khanikaev A B, Genack A Z. Pseudo-spin-valley coupled edge states in a photonic topological insulator. Nature Communications, 2018, 9(1): 3029

PMID

8 Berry M V. Quantal phase factors accompanying adiabatic changes. Proceedings of the Royal Society of London, Series A, 1802, 1984(392): 45–57

9 Tomita A, Chiao R Y. Observation of Berry’s topological phase by use of an optical fiber. Physical Review Letters, 1986, 57(8): 937–940

PMID

10 Wang H X, Guo G Y, Jiang J H. Band topology in classical waves: Wilson-loop approach to topological numbers and fragile topology. New Journal of Physics, 2019, 21(9): 093029

11 Hatsugai Y. Chern number and edge states in the integer quantum Hall effect. Physical Review Letters, 1993, 71(22): 3697–3700

PMID

12 Wang Z, Chong Y D, Joannopoulos J D, Soljacić M. Reflection-free one-way edge modes in a gyromagnetic photonic crystal. Physical Review Letters, 2008, 100(1): 013905

PMID

13 Skirlo S A, Lu L, Soljačić M. Multimode one-way waveguides of large Chern numbers. Physical Review Letters, 2014, 113(11): 113904

PMID

14 Yang B, Zhang H F, Wu T, Dong R, Yan X, Zhang X. Topological states in amorphous magnetic photonic lattices. Physical Review B, 2019, 99(4): 045307

15 Ochiai T, Onoda M. Photonic analog of graphene model and its extension: Dirac cone, symmetry, and edge states. Physical Review B, 2009, 80(15): 155103

16 Fukui T, Hatsugai Y, Suzuki H. Chern numbers in discretized Brillouin Zone: efficient method of computing (spin) hall conductances. Journal of the Physical Society of Japan, 2005, 74(6): 1674–1677

17 Ringel Z, Kraus Y E. Determining topological order from a local ground-state correlation function. Physical Review B, 2011, 83(24): 245115

18 Yu R, Qi X L, Bernevig A, Fang Z, Dai X. Equivalent expression of Z(2) topological invariant for band insulators using the non-Abelian Berry connection. Physical Review B, 2011, 84(7): 075119

19 Ma T, Shvets G. All-Si valley-Hall photonic topological insulator. New Journal of Physics, 2016, 18(2): 025012

20 Ma T, Shvets G. Scattering-free edge states between heterogeneous photonic topological insulators. Physical Review B, 2017, 95(16): 165102

21 Ye L, Yang Y T, Hang Z H, Qiu C Y, Liu Z Y. Observation of valley-selective microwave transport in photonic crystals. Applied Physics Letters, 2017, 111(25): 251107

22 Gao F, Xue H R, Yang Z J, Lai K, Yu Y, Lin X, Chong Y, Shvets G, Zhang B. Topologically protected refraction of robust kink states in valley photonic crystals. Nature Physics, 2018, 14(2): 140–144

23 Xiao D, Yao W, Niu Q. Valley-contrasting physics in graphene: magnetic moment and topological transport. Physical Review Letters, 2007, 99(23): 236809

PMID

24 Shalaev M I, Walasik W, Tsukernik A, Xu Y, Litchinitser N M. Robust topologically protected transport in photonic crystals at telecommunication wavelengths. Nature Nanotechnology, 2019, 14(1): 31–34

PMID

25 He X T, Liang E T, Yuan J J, Qiu H Y, Chen X D, Zhao F L, Dong J W. A silicon-on-insulator slab for topological valley transport. Nature Communications, 2019, 10(1): 872

PMID

26 Joannopoulos J D, Johnson S G, Winn J N, Meade R D. Photonic Crystals Molding the Flow of Light. 2nd ed. America: Princeton University Press, 2008, 1–283

27 Yang B, Wu T, Zhang X. Engineering topological edge states in two dimensional magnetic photonic crystal. Applied Physics Letters, 2017, 110(2): 021109

28 Chan H C, Guo G Y. Tuning topological phase transitions in hexagonal photonic lattices made of triangular rods. Physical Review B, 2018, 97(4): 045422