General corner charge formula in two-dimensional $C_{n}$-symmetric higher-order topological insulators

We consider a crystal of finite size which preserves the $C_{n}$ symmetry. We assume that the system is insulating, in the sense that the Fermi level is within the energy gap for both the bulk and the boundaries. In some cases, degenerate corner states exist at the Fermi energy, which violates the above assumption, when charge neutrality is assumed. However, we can add or remove several electrons so that no states are at the Fermi energy. This deviation from charge neutrality is called filling anomalyPhysRevB.99.245151 .

In order to calculate the filling anomaly, we divide the system conceptually into two $C_{n}$-symmetric regions (I) and (II), as shown in Fig. 1. The regions (I) and (II) can be regarded as surface and bulk regions, respectively, and we explain this division in the following. The width of the surface region (I), $d_{\text{edge}}$, is taken to be large enough so that both the electronic states and the ionic positions in the region (II) are the same as those at the bulk.

To calculate the filling anomaly, we proceed with the following steps. First, we classify the ions into two groups according to which region they belong to, so that no ions are at the junctures of two regions.

Second, we classify the electronic states. We assume that the occupied states of this system can be described by a $C_{n}$-symmetric basis set $\{\phi_{i}\}$ composed of exponentially localized basis orbitalsPhysRevB.48.4442 . Well within the crystal, the electronic states can be well described by bulk Wannier functions $W_{nl}$, and we take these Wannier functions $W_{nl}$ to be among $\{\phi_{i}\}$. The remaining $\phi_{i}$’s are additional localized basis orbitals $\varphi_{i}^{(\text{add})}$ near the boundary. We classify them into regions (I) and (II) so that the set of orbitals in region (I) and that in region (II) are $C_{n}$ symmetric. By construction the numbers of ions and electronic localized basis orbitals in the region (I) are integer multiples of $n$. Here, we impose a condition that all the additional localized basis orbitals $\varphi_{i}^{(\text{add})}$ are classified into the region (I), and we choose the region (I) large enough so that this condition is satisfied. Thus, even if the electronic states near the boundary are largely deformed, our theory is still applicable, by incorporating these deformed wavefunctions into region (I). We note that this construction of localized basis functions is similar to that used in the argument in Ref. [PhysRevB.48.4442, ] to show that the bulk polarization in terms of the Zak phase gives surface charge density. We note that some degrees of freedom remain in this division into two regions, but they do not affect the following argument.

Third, we calculate the total charge of the finite crystal in terms of mod $n|e|$. First, the total charge contribution from the region (I) should be $0$ (mod $n|e|$) because the number of the nuclei and the electronic basis orbitals within region (I) is an integer multiple of $n$. Next, we calculate total charges of electrons and ions in the region (II). First, in the region (II), all the orbitals are Wannier orbitals $W_{nl}$, and therefore positions of the ions and electrons are classified in terms of Wyckoff positions. Because the whole crystal is $C_{n}$ symmetric, the center of the crystal should be located at the Wyckoff position with $C_{n}$ local symmetry. Let $1a$ be the Wyckoff position with $C_{n}$ symmetry located at the center of the crystal. Then, the ions not located at $1a$ contribute $nM|e|$ ($M\in\mathbb{Z}$) to the net charge in the region (II), while those at 1a contribute $nM^{\prime}|e|+n_{1a}^{(\text{ion})}|e|$ ($M^{\prime}\in\mathbb{Z}$), where $n_{1a}^{(\text{ion})}|e|$ is the charge of the ion at the $1a$ position. Therefore, the total ionic charge in the region (II) is equal to $n_{1a}^{(\text{ion})}|e|$ mod $n|e|$ and so is that in the whole crystal. Similarly, the total electron charge is equal to $-n_{1a}^{\text{(e)}}|e|$ (mod $n|e|$) where $n_{1a}^{\text{(e)}}$ is the number of electron Wannier orbitals at the $1a$ position. From the above discussion, the total charge of the system in terms of modulo $n|e|$ is calculated as follows,

$\displaystyle Q_{\text{tot}}\equiv(n_{1a}^{\text{(ion)}}-n_{1a}^{\text{(e)}})|e|\quad(\text{mod}\ n|e|),$

(1)

where $n_{1a}^{\text{(ion)}}|e|$ is an ionic charge at $1a$, and $n_{1a}^{\text{(e)}}$ is the number of electronic Wannier functions at $1a$. This formula represents a filling anomaly in the crystal.

Fourth, we discuss how $Q_{\text{tot}}$ is distributed in the crystal. As a first step, we define the moving average of the total charge density $\overline{\rho}(\bm{x})$ as follows:

$\displaystyle\overline{\rho}(\bm{x})=\frac{1}{|\Omega|}\int_{\Omega}\rho(\bm{x}+\bm{x}^{\prime})\ d\bm{x}^{\prime}.$

(2)

Here, $\rho(\bm{x})$ is the total charge density including ions and electrons, $\Omega$ is a $C_{n}$-symmetric unit cell whose center is at the origin and $|\Omega|$ is the area of $\Omega$. From the definition, $\overline{\rho}(\bm{x})$ have the following properties:

	$\displaystyle\overline{\rho}(\bm{x})$	$\displaystyle=0\quad(\bm{x}\in\text{bulk}),$		(3)
	$\displaystyle\int_{\text{crystal}}d\bm{x}^{\prime}\ \overline{\rho}(\bm{x}^{\prime})\$	$\displaystyle=Q_{\text{tot}}.$		(4)
	$\displaystyle\overline{\rho}(C_{n}\cdot\bm{x})$	$\displaystyle=\overline{\rho}(\bm{x})$		(5)

We note that as long as $\overline{\rho}(\bm{x})$ satisfies Eqs. (3), (4) and (5), we can also use other definitions for moving averages. In a sufficiently large system, we expect that electronic states in the middle part of the edges far away from the corners are almost the same as those on an edge of a semi-infinite system with an open boundary condition in one direction. To calculate the corner charge based on this expectation, we divide the system into three types of regions: the (A) corner, (B) edge, and (C) bulk regions in a $C_{n}$-symmetric way as schematically shown in Fig. 2. The edge region (B) should be chosen to satisfy the following conditions. First, it should be sufficiently far away from the corners, such that the electronic states in the edge region (B) are the same as those on an edge of the corresponding semi-infinite system. Second, it should have a width equal to an integer multiple of the period along the edge, $N_{\text{edge}}a_{\text{edge}}$, and a sufficiently long depth into the bulk, $d_{\text{edge}}$, such that in the remaining bulk region (C), the electronic states are the same as those for a bulk infinite system. Here, we assume that one edge region includes $N_{\text{edge}}$ edge unit cells. Since $\overline{\rho}(\bm{x})$ takes nonzero values only in the (A) corner and (B) edge regions, its integral can be calculated as the sum of the total charges in these regions:

$\displaystyle\int_{\text{crystal}}d\bm{x}^{\prime}\ \overline{\rho}(\bm{x}^{\prime})\$

$\displaystyle=n\Bigl{(}Q_{\text{edge}}+Q_{\text{corner}}\Bigr{)},$

(6)

where $Q_{\text{edge}}$ and $Q_{\text{corner}}$ are the total charges at one of the (B) edge regions and one of the (C) corner regions, respectively. Then $Q_{\text{edge}}$ can be calculated as:

$\displaystyle Q_{\text{edge}}$

$\displaystyle=N_{\text{edge}}\sigma_{\text{edge}}a_{\text{edge}}.$

(7)

Here, $N_{\text{edge}}$ is the number of edge unit cells within one edge region, $\sigma_{\text{edge}}$ is the edge charge density, and $a_{\text{edge}}$ is the period of the edge unit cell.

Finally, we calculate the corner charge. We assume that the edge is charge neutral, $\sigma_{\text{edge}}=0$. Then the extra charge $Q_{\text{tot}}$ is distributed into the $n$ corners. Since the crystal is assumed to be $C_{n}$-symmetric, the distributed charges are equal among the corners, and this gives a corner charge formula:

$\displaystyle Q_{\text{corner}}=Q_{c\_1a}^{(n)}\equiv\frac{(n_{1a}^{\text{(ion)}}-n_{1a}^{\text{(e)}})|e|}{n}\quad(\text{mod}\ |e|).$

(8)

Here, $Q_{c\_1a}^{(n)}$ means the corner charge when the center of the crystal is located at the $1a$ Wyckoff position.

If the edge charge density $\sigma_{\text{edge}}\neq 0$, the extra charge $Q_{\text{tot}}$ is distributed not only to the corners, but also to the edges. Since the separation of charges between the corners and the edges depends on the choice of the corner regions, the corner charge is not well defined.

Our discussion so far clarifies that the formula (8) for the corner charge applies to a general system, as long as the system is $C_{n}$-symmetric and is insulating over the whole system, irrespective of the system details. In particular, our theory applies to (i) systems with boundaries with high Miller indices, (ii) those with surface reconstructions, and (iii) systems with nonzero $\bm{P}$ but with zero edge charge density $\sigma_{\text{edge}}=0$. In the following, we briefly explain these three cases. Details of (i) and (ii) are given in Sec. IV and those of (iii) in Sec. III.

First is the case (i) with boundaries with higher Miller indices. In the previous works, the arguments are limited to boundaries with low Miller indices such as (10) and (01) edges, but our theory applies to general Miller indices. For example, a $C_{4}$-symmetric square system with edges along (34) and ($4\bar{3}$) directions is covered by our theory. Second, we can treat cases with (ii) boundaries with surface reconstructions. It is nontrivial because the surface reconstruction may make the perioidicity at edges different from that in the bulk. In some cases, as we discusss later, the edges without the surface reconstruction are gapless and $\sigma_{\text{edge}}$ is zero, but those with the surface reconstruction may become gapped, which makes our theory applicable.

Third, our theory covers also the cases (iii) with $\bm{P}\neq 0$. Our assumption is $\sigma_{\text{edge}}=0$, which does not necessarily mean $\bm{P}=0$. In contrast, in the previous worksPhysRevB.99.245151 ; PhysRevResearch.1.033074 ; PhysRevResearch.2.043131 ; PhysRevB.102.165120 , the bulk polarization $\bm{P}$ is assumed to be 0, which ensures that the edge charge density $\sigma_{\text{edge}}$ is 0. We find that this assumption of $\bm{P}=0$ in the previous works is too strong. The bulk polarization $\bm{P}$ is related to the edge charge densityPhysRevB.48.4442 ; resta1992theory ; PhysRevB.47.1651 ; RevModPhys.66.899 via

$\displaystyle\sigma_{\text{edge}}\equiv\bm{P}\cdot\bm{n}\ (\text{mod}\ |e|/a_{\text{edge}}),$

(9)

where $\bm{n}$ is the unit normal vector of the edge. It is not an equality, because the bulk polarization $\bm{P}$ is defined only in terms of modulo some unit PhysRevB.48.4442 ; resta1992theory ; PhysRevB.47.1651 ; RevModPhys.66.899 , which we call polarization quantum. Thus, our assumption $\sigma_{\text{edge}}=0$ means $\bm{P}\cdot\bm{n}\equiv 0$ (mod $|e|/a_{\text{edge}}$), which does not necessarily mean $\bm{P}=0$.

In this section, we rewrite our formula for the corner charge in terms of rotation eigenvalues of the Bloch eigenstates at high-symmetry wavevectors following Refs. [PhysRevB.99.245151, ,PhysRevResearch.1.033074, ,PhysRevResearch.2.043131, ], in three different systems with and without time-reversal symmetry, and with and without spin-orbit coupling, separately. We need to rewrite it, because eigenstates in a crystal are Bloch wavefunctions, extended over the crystal, and one cannot directly calculate $n_{1a}^{(e)}$ from the Bloch wavefunctions. Then, we discuss a relation between our formula and the formulas in the previous worksPhysRevB.99.245151 ; PhysRevResearch.1.033074 ; PhysRevResearch.2.043131 . The notable difference from the previous works is that the Wyckoff positions of the ions are not always 1a, and the corner charge generally depends on the filling $\nu$ as in Ref. [PhysRevResearch.2.043131, ]. The detail of the calculation is shown in the Appendix.

By following Refs. [PhysRevB.99.245151, ,PhysRevResearch.1.033074, ,PhysRevResearch.2.043131, ], we introduce topological invariants which distinguish the Wannier functions with different Wyckoff positions. With a sufficient number of topological invariants, there can be a one-to-one correspondence between the values of the topological invariants and Wyckoff positions of the occupied Wannier orbitals. We assume that the occupied bands can be written in terms of Wannier functions. Then, the number of Wannier orbitals at each Wyckoff position is determined from the topological invariant for the occupied bands.

Figure 3 shows the Wyckoff positions under $C_{n}$ rotational symmetry ($n=3,4,6$). In each system with $C_{n}$ symmetry, there are four kinds of Wyckoff positions. We define the number of occupied Wannier functions at Wyckoff position $mX$ as $n_{mX}^{(e)}$. Here, $X=a,b,c,d$ represent the types of the Wyckoff positions, and $m$ represents the multiplicity of a Wyckoff position. The weighted sum of the $n_{mX}^{(e)}$ is equal to the electron filling, i.e. the number of occupied bands $\nu$:

$\displaystyle\sum_{X}mn_{mX}^{(e)}=\nu.$

(10)

For the calculation of the corner charge, the value of $n_{1a}^{(e)}$ is needed. As we show in the appendix, $mn_{mX}^{(e)}$ modulo $n$ ($X\neq a$) is directly calculated from the rotational topological invariants defined from rotation eigenvalues of the Bloch wavefunctions, and $n_{1a}^{(e)}$ is determined from $\nu-\sum_{X\neq a}mn_{mX}^{(e)}$. In the appendix, we show the detailed calculation of $n_{mX}^{(e)}$.

Here we consider systems with surface reconstruction. Conceptual pictures of the surface reconstruction is shown in Fig. 6. We assume that the system has nonzero quantized electric polarization in the bulk. Then the edge have nonzero fractional edge charge $\sigma_{\text{edge}}$. In Fig. 6(a) we show its band structure when the edge states are metallic. However, it can under go Peierls transition or surface reconstruction to multiply edge periodicity and open a gap. In Fig. 6(b), due to the surface reconstruction, the surface periodicity is doubled, and the edge charge becomes an integer in the enlarged surface unit cells. Then, the edge energy spectrum can be gapped and insulating. In this case, one can apply our general formula.

Here, we consider the surface with a high Miller index. First, we consider the $C_{4}$-symmetric system with basis vectors of $\bm{a}_{1}=(a,0)$ and $\bm{a}_{2}=(0,a)$, as shown in Fig. 3(b). In the $C_{4}$-symmetric system, the bulk electric polarization is quantized to be $\bm{P}\equiv(0,0)$ or $(\frac{e}{2a},\frac{e}{2a})$ (mod $\frac{e}{a}(m_{1},m_{2})$), ($m_{1},m_{2}\in\mathbb{Z}$)PhysRevB.99.245151 ; PhysRevB.102.165120 ; PhysRevB.48.4442 . In the following, we consider the case with $\bm{P}\equiv(\frac{e}{2a},\frac{e}{2a})$. Next, for the surface with a Miller index ($\alpha,\beta$) ($\alpha,\beta\in\mathbb{Z}$, and $\alpha$ and $\beta$ are coprime), the surface normal vector $\bm{n}$ is defined as follows:

$\displaystyle\bm{n}=\Bigg{(}\frac{\alpha}{\sqrt{\alpha^{2}+\beta^{2}}},\ \frac{\beta}{\sqrt{\alpha^{2}+\beta^{2}}}\Bigg{)}.$

(37)

Therefore, the edge charge in one edge unit cell is calculated as follows:

$\displaystyle\sigma a_{\text{edge}}\equiv\bm{P}\cdot\bm{n}a_{\text{edge}}\equiv\frac{e}{2}(\alpha+\beta)\ (\text{mod}\ e).$

(38)

Here, we set the size of the edge unit cell to be the minimal one, $a_{\text{edge}}=|\beta\bm{a}_{1}-\alpha\bm{a}_{2}|=\sqrt{\alpha^{2}+\beta^{2}}a$, as shown in Fig. 7. If both $\alpha$ and $\beta$ are odd integers, $\sigma a_{\text{edge}}\equiv 0$ (mod $e$), and otherwise $\sigma a_{\text{edge}}\equiv e/2$ (mod $e$). Obviously, the edge can be charge neutral when $\sigma a_{\text{edge}}\equiv 0$ (mod $e$). Even when $\sigma a_{\text{edge}}\equiv e/2$ (mod $e$), after surface reconstruction with doubling the surface period $a_{\text{edge}}^{\prime}=2a_{\text{edge}}$, we get $\sigma a_{\text{edge}}^{\prime}\equiv 0$ (mod $e$), and the edge can be charge neutral. In particular, when the system has edge states in the gap, and the fractional edge charge is accommodated in the edge states, the edge states are half-filled and it always has Peierls instability, meaning that the edge becomes insulating via spontaneous lattice deformation with doubling the edge unit cell at sufficiently low temperature.

Next, we consider the $C_{3}$-symmetric system with basis vectors of $\bm{a}_{1}=(\frac{a}{2},\frac{\sqrt{3}}{2}a)$ and $\bm{a}_{2}=(-\frac{a}{2},\frac{\sqrt{3}}{2}a)$, as shown in Fig. 3(a). In the $C_{3}$-symmetric system, the bulk electric polarization is $\bm{P}=\bm{0}$, $\frac{e}{3\Omega}(\bm{a}_{1}+\bm{a}_{2})$ or $\frac{-e}{3\Omega}(\bm{a}_{1}+\bm{a}_{2})$ (mod $\frac{e}{\Omega}\bm{R}$)PhysRevB.99.245151 ; PhysRevB.102.165120 ; PhysRevB.48.4442 , where $\Omega=\frac{\sqrt{3}}{2}a^{2}$ is the area of the bulk unit cell and $\bm{R}=m_{1}\bm{a}_{1}+m_{2}\bm{a}_{2}$ ($m_{1},m_{2}\in\mathbb{Z}$) is a lattice vector. In the following, we consider the case with $\bm{P}=\frac{\pm e}{3\Omega}(\bm{a}_{1}+\bm{a}_{2})=(0,\frac{\pm 2e}{3a})$. The surface normal vector $\bm{n}$ for the surface with Miller index ($\alpha,\beta$) ($\alpha,\beta\in\mathbb{Z}$, and $\alpha$ and $\beta$ are coprime) is defined as follows:

$\displaystyle\bm{n}$

$\displaystyle=\frac{1}{\sqrt{\alpha^{2}-\alpha\beta+\beta^{2}}}\Bigg{(}\frac{\sqrt{3}}{2}(\alpha-\beta),\ \frac{1}{2}(\alpha+\beta)\Bigg{)}.$

(39)

Therefore, the edge charge in one edge unit cell is calculated as follows:

$\displaystyle\sigma a_{\text{edge}}\equiv\bm{P}\cdot\bm{n}a_{\text{edge}}\equiv\pm\frac{e}{3}(\alpha+\beta)\ (\text{mod}\ e).$

(40)

Here, we used $a_{\text{edge}}=|\beta\bm{a}_{1}-\alpha\bm{a}_{2}|=\sqrt{\alpha^{2}-\alpha\beta+\beta^{2}}a$. If $\alpha+\beta\in 3\mathbb{Z}$, $\sigma a_{\text{edge}}\equiv 0$ (mod $e$), and otherwise $\sigma a_{\text{edge}}\equiv e/3$ or $(2e)/3$ (mod $e$). Obviously, edge can be charge neutral when $\sigma a_{\text{edge}}\equiv 0$ (mod $e$). Even when $\sigma a_{\text{edge}}\equiv e/3$ or $(2e)/3$ (mod $e$), after surface reconstruction with tripling the surface period $a_{\text{edge}}^{\prime\prime}=3a_{\text{edge}}$, we get $\sigma a_{\text{edge}}^{\prime\prime}\equiv 0$ (mod $e$), and the edge can be charge neutral.

Finally, we consider the $C_{6}$-symmetric system with basis vectors of $\bm{a}_{1}=(a,0)$ and $\bm{a}_{2}=(-\frac{a}{2},\frac{\sqrt{3}}{2}a)$, as shown in Fig. 3(c). In the $C_{6}$-symmetric system, the bulk electric polarization is always $\bm{P}=\bm{0}$ (mod $\frac{e}{\Omega}\bm{R}$)PhysRevB.99.245151 ; PhysRevB.102.165120 ; PhysRevB.48.4442 , where $\Omega=\frac{\sqrt{3}}{2}a^{2}$ is the area of the bulk unit cell and $\bm{R}=(m_{1}a-\frac{m_{2}}{2}a,\frac{\sqrt{3}}{2}m_{2}a)$ ($m_{1},m_{2}\in\mathbb{Z}$) is a lattice vector. Therefore, $\sigma_{\text{edge}}=\bm{P}\cdot\bm{n}$ is 0 for surfaces in any direction (mod $e/a_{\text{edge}}$).

To summarize, in either case, even when $\bm{P}\neq 0$ in the $C_{n}$-symmetric systems, $\sigma a_{\text{edge}}\equiv 0$ (mod $e$) can be achieved with proper surface reconstruction thanks to the fractional quantization of $\bm{P}$.

Here, we will explain how to apply our formula to tight-binding models. In principle, tight-binding models contain only electron degrees of freedom. In order to calculate the corner charge, it is necessary to give appropriate information on the positions of ions.

First, it should be pointed out that the sites of the tight-binding models are different from the positions of ions. The apatite electride is one of the examples showing that these two are independent PhysRevResearch.2.043131 . The electronic states near the Fermi level of the apatite electride are well described by the Wannier functions localized at the one-dimensional hollows in the crystal. Therefore the tight-binding model made from these Wannier orbitals has its sites inside the hollows. Since there are no ions in the hollows, the ion positions and the sites are clearly different in this model.

Next, we apply our formula to the well-known 2D Benalcazar-Bernevig-Hughes (BBH) model Benalcazar61 ; PhysRevB.96.245115 as an example. This model is a spinless tight-binding model on a square lattice with $\pi$-flux, showing a higher-order topological insulating phase. As in Ref. [PhysRevB.96.245115, ], we assume that all the positive ions are located at the center of the $C_{4}$-symmetric unit cell, i.e., $n_{1a}^{(\text{ion})}=\nu$ holds. The $k$-space Hamiltonian is given as follows:

	$\displaystyle\mathcal{H}(\bm{k})$	$\displaystyle=(\gamma+\lambda\cos k_{x})\tau_{x}\sigma_{0}-\lambda\sin k_{x}\tau_{y}\sigma_{z}$
		$\displaystyle\quad-(\gamma+\lambda\cos k_{y})\tau_{y}\sigma_{y}-\lambda\sin k_{y}\tau_{y}\sigma_{x}.$		(41)

Here, $\gamma$ and $\lambda$ are real parameters, and $\tau_{j}$, $\sigma_{j}$ ($j=x,y,z$) are Pauli matrices. For simplicity, we take the lattice constants as $a=1$. The Hamiltonian has $C_{4}$-symmetry $\hat{C_{4}}$:

$\displaystyle\hat{C_{4}}$

$\displaystyle=\begin{pmatrix}0&\sigma_{0}\\ -i\sigma_{y}&0\end{pmatrix}_{\tau}.$

(42)

Here, due to the $\pi$-flux, $\hat{C_{4}}$ satisfies $(\hat{C_{4}})^{4}=-1$, and its eigenvalues are $C_{4}=e^{\frac{(2n+1)\pi i}{4}}$ ($n=0,1,2,3$). In the derivation of the corner charge formula in this paper, we assumed that $(\hat{C_{n}})^{n}=1$ for spinless systems. Therefore, our formula cannot be applied as it is, when the rotation eigenvalues are modified due to the $\pi$-flux. However, we can formally determine the corner charge by applying the formula for the spinful class A system as if the model described a spinful system because $(\hat{C_{4}})^{4}=-1$.

For simplicity, we consider two limits: $\lambda=0$ and $\gamma=0$. When $\lambda=0$, the Hamiltonian becomes $\mathcal{H}(\bm{k})=\gamma(\tau_{x}\sigma_{0}-\tau_{y}\sigma_{y})$, which does not depend on $k$. Then topological invariants are trivial: $[X_{1}^{(2)}]=[M_{1}^{(4)}]=[M_{2}^{(4)}]=0$. When $\gamma=0$, the Hamiltonian becomes

	$\displaystyle\mathcal{H}(\bm{k})$	$\displaystyle=\sqrt{2}\lambda\begin{pmatrix}0&Q^{\dagger}(\bm{k})\\ Q(\bm{k})&0\end{pmatrix},$		(43)
	$\displaystyle Q(\bm{k})$	$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}e^{-ik_{x}}&-e^{ik_{y}}\\ e^{-ik_{y}}&e^{ik_{x}}\end{pmatrix}.$		(44)

Topological invariants can be calculated by taking the simultaneous eigenstates of the Hamiltonian and the rotation operation $\hat{C_{4}}$ or $\hat{C_{2}}=(\hat{C_{4}})^{2}$ at each of the high symmetry points $\Gamma=(0,0)$, $X=(\pi,0)$, and $M=(\pi,\pi)$. The results are summarized in Table. 1, and the corner charge is calculated from the Eq.(12):

$\displaystyle Q_{c\_1a}^{(4)}\equiv\frac{|e|}{4}\Big{(}2-2+0-\frac{1}{2}\cdot 1+\frac{3}{2}\cdot(-1)\Big{)}\equiv\frac{e}{2}\ (\text{mod}\ e).$

(45)

This result agrees with the results of previous studiesBenalcazar61 ; PhysRevB.96.245115 . This value of the corner charge is stable as long as $|\gamma|<|\lambda|$, where the bulk gap is open.