Anomalous Chern-Simons orbital magnetoelectric coupling of three-dimensional Chern insulators: gauge-discontinuity formalism and adiabatic pumping

The topological orbital magnetoelectric coupling for insulators with vanishing Chern number is given by the integral of the Chern-Simons 3-form [1, 6, 8, 7]

$$\theta_{\text{CS}}=-\frac{1}{4\pi}\int d^{3}k\,\epsilon_{abc}\mathrm{Tr}[A_{a}\partial_{b}A_{c}-\frac{2i}{3}A_{a}A_{b}A_{c}]$$

(3)

where $A_{a,mn}=\braket{u_{m\mathbf{k}}}{\partial_{{}_{a}}}{u_{n\mathbf{k}}}$ is the Berry connection matrix element of the occupied Bloch states at wavevector $\mathbf{k}$, and the notation $\partial_{a}\equiv\partial_{k_{a}}$ with $a=x,y,z$. $|u_{n\mathbf{k}}\rangle$ ($|u_{m\mathbf{k}}\rangle$) denotes the periodic part of the Bloch state at $\mathbf{k}$ with band index $n$ ($m$). Mathematically, the integral of the three form’s exterior derivative over a 4-dimensional compact manifold, which can be interpreted as the “four-dimensional Brillouin zone” of a 3D insulator under a cyclic adiabatic parameter $\eta$,

$$C^{(2)}\equiv\frac{1}{32\pi^{2}}\int d^{3}k\,d\eta\,\epsilon^{abcd}\mathrm{Tr}(F_{ab}F_{cd})\;,$$

(4)

is the second Chern number, with $a,b,c,d=x,y,z,\eta$. The generalized Stokes’ theorem would thus show that a non-trivial pumping of $\theta$ through a cyclic adiabatic evolution is dictated by a non-zero second Chern number.

Nevertheless, the Chern-Simons magnetoelectric coupling $\theta$ (Eq. (3)) is well studied only for 3D insulators with vanishing Chern numbers in every 2D cut of the 3D Brillouin zone. Several alternative representations of Chern-Simons coupling $\theta$ have been proposed. For example, it has been shown in [19] that $\theta$ can be expressed in the hybrid Wannier representation,

	$\displaystyle\theta=$	$\displaystyle-\frac{1}{d_{0}}\int d^{2}k\sum_{n}\bar{z}_{n}\,\Omega_{xy,0n,0n}$
		$\displaystyle+\frac{i}{d_{0}}\int d^{2}k\sum_{lmn}(\bar{z}_{lm}-\bar{z}_{0n})A_{x,0n,lm}A_{lm,0n}$		(5)

where $\Omega_{xy,0n,0n}$, $A_{a,0n,lm}$ are the Berry curvature and Berry connection matrix elements defined in the hybrid Wannier basis, $\bar{z}_{n}$ is the $n$th hybrid Wannier center, and $d_{0}$ is the lattice constant in $z$ direction. This draws the connection between Berry curvature (Berry connection), Wannier center and Chern-Simons coupling $\theta$. However, still it only applies to insulators with vanishing Chern numbers.

Since Chern-Simons 3-form involves Berry connection $A_{a}$, it becomes ill-defined at some points in the 3D Brillouin zone when a global smooth and periodic gauge of the occupied Bloch states is missing. This is exactly the case for 3D Chern insulators, where the Chern number of a given 2D $\mathbf{k}$ plane (say, the $(k_{x},k_{y})$ plane) is a nonzero integer $C$. The topological obstruction imposed by nonzero Chern number necessarily introduces singularities or discontinuities of the gauge of occupied Bloch states somewhere in the Brillouin zone.

Therefore, in order to study the Chern-Simons OME coupling of 3D Chern insulators, we are forced to consider the contribution of “gauge discontinuity” somewhere in the 3D Brillouin zone. In Ref. 15, a new formalism to treat the Chern-Simons OME coupling has been proposed. For example, in the case of 3D topological insulators, nontrivial $\mathbb{Z}_{2}$ band topology protected by $\mathcal{T}$ symmetry would cause obstructions in constructing a smooth and periodic gauge respecting $\mathcal{T}$ symmetry, making it hard to numerically achieve a half quantized value $\theta=\pi$. To solve this problem, the authors of Ref. 15 proposed to retain the smoothness of a $\mathcal{T}$-preserving gauge, while relaxing the periodicity condition, which introduces some “gauge discontinuity” at certain 2D boundary of the 3D Brillouin zone. Then, the contributions from the gauge-discontinuity plane is further taken into account in addition to the bulk integral of Chern-Simons 3-form. Since the periodicity condition has been relaxed, the gauge can be made as smooth as possible in the bulk of 3D Brillouin zone, which makes the numerical integral of the bulk Chern-Simons 3-form much easier. Indeed, in the case of 3D topological insulators, it has been shown that the half-quantized $\theta=\pi$ is all contributed by one particular term originating from the gauge-discontinuity plane (known as “vortex loop” term, to be elaborated later) [15], and the bulk integral of Chern-Simons 3-form vanishes due to the choice of a smooth and $\mathcal{T}$-preserving gauge.

In the case of 3D Chern insulators, we are facing a similar situation: one cannot make a smooth and periodic gauge (even for single occupid band) throughout the Brillouin zone due to topological obstructions imposed by nonzero Chern number. Therefore, we would like to make the gauge as smooth as possible in the bulk of 3D Brillouin zone, while relaxing the periodicity condition of the gauge at a certain 2D boundary. Then, we will apply the gauge-discontinuity formalism of Chern-Simons 3-form to the case of 3D Chern insulators, and consider the bulk and gauge-discontinuity contributions separately.

In this subsection, we briefly review the gauge-discontinuity formalism introduced in Ref. 15. We first specify a gauge that is convenient for our analysis. Throughout the paper we use reduced coordinates $k_{x}$, $k_{y}$, $k_{z}\in[0,1)$. Consider a gauge with discontinuity at$k_{y}=1$ but smooth along $k_{z}$ and $k_{x}$.

	$\displaystyle\ket{\psi_{n}(k_{x},k_{y}=1,k_{z})}=$	$\displaystyle\sum_{m=1}^{M}W_{mn}(k_{x},k_{z})$
		$\displaystyle\times\ket{\psi_{m}(k_{x},k_{y}=0,k_{z})}$		(6)

where $W=e^{iB}$ is an $M\times M$ unitary matrix defined in the basis of the $M$ occupied Bloch states $\{\ket{\psi_{n}(k_{x},k_{y},k_{z})}\}$, characterizing the gauge discontinuity across the $k_{y}=0,1$ plane and $B$ is a Hermitian matrix. To handle the discontinuity, the formalism introduces a fictitious $\lambda$-space, in which we linearly interpolate the discontinuity from $k_{y}=1$ to $k_{y}=0$,

$$\ket{\psi_{n}(k_{x},\lambda,k_{z})}=\sum_{m=1}^{M}W_{mn}\ket{\psi_{m}(k_{x},\lambda=0,k_{z})}$$

(7)

where

$$\ket{\psi_{n}(k_{x},\lambda=0,k_{z})}=\ket{\psi_{n}(k_{x},k_{y}=1,k_{z})}$$

(8)

In the $\lambda$-space, the Berry connection and Berry curvature are defined as usual with respect to $\ket{u(k_{x},\lambda,k_{z})}$. With above definition, we have

	$\displaystyle\ket{\psi_{n}(k_{x},\lambda=1,k_{z})}=\ket{\psi_{n}(k_{x},k_{y}=0,k_{z})}$		(9a)
	$\displaystyle\ket{u(k_{x},\lambda=1,k_{z})}=e^{-i2\pi\hat{y}}\ket{u(k_{x},k_{y}=0,k_{z})}$		(9b)

Namely, we have removed the discontinuity by smoothly connecting the $k_{y}=0$ and $k_{y}=1$ planes with a fictitious $\lambda$ space, as schematically shown in Fig. 1(a),

The Chern-Simons OME thus has three contributions

$$\theta=\theta_{\text{BK}}+\theta_{\text{GD}}+\theta_{\text{VL}}$$

(10)

where $\theta_{\text{BK}}$ is the contribution from the integral of Chern-Simons 3-form within the bulk Brillouin zone, but without enforcing periodicity at the $k_{y}=0,1$ plane. $\theta_{\text{GD}}$ and $\theta_{\text{VL}}$ are referred as the “gauge-discontinuity” and “vortex-loop” terms both arising from the contribution of the fictitious $\lambda$ space connecting $k_{y}=0$ and $k_{y}=1$ planes. Following directly from the Eqs. (3), after some integrations by part, the first two terms in Eq. (10) are given by

	$\displaystyle\theta_{\textrm{BK}}=$	$\displaystyle-\frac{1}{4\pi}\int d^{3}k\,\mathrm{Tr}[A_{x}\Omega_{yz}$
		$\displaystyle+A_{y}\Omega_{zx}+A_{z}\Omega_{xy}-2i[A_{x},A_{y}]A_{z}]$		(11a)
	$\displaystyle\theta_{\text{GD}}=$	$\displaystyle-\frac{1}{4\pi}\int dk_{z}\,dk_{x}\,d\lambda\,\mathrm{Tr}[A_{z}\Omega_{x\lambda}$
		$\displaystyle+A_{x}\Omega_{\lambda z}+A_{\lambda}\Omega_{zx}-2i[A_{z},A_{x}]A_{\lambda}]$		(11b)

where the trace is carried out over the occupied bands. We note that the unitary gauge discontinuity matrix $W=e^{iB}$ can be equivalently described by the Hermitian matrix $B$ in its exponent. The eigenvalues of $B$ thus have $2\pi$ ambiguities. Sometimes, it is impossible to insist a branch choice such that all eigenvalues of $B$ would remain smooth in $k_{z}$ and $k_{x}$ plane. Given a fixed branch choice, there can be a sudden jump of $2\pi\nu_{n}$ ($\nu_{n}$ being integer) in the $n$th eigenvalue of $B$ (denoted by $\beta_{n}$) along certain “vortex loop” or “vortex line” in the $k_{z}$-$k_{x}$ plane. Then, we need to account for contributions from such $2\pi$ jumps as well, necessitating the vortex-loop contribution, expressed as [15]

$$\theta_{\text{VL}}=\sum_{n}\frac{\phi_{n}(\mathcal{C})+\xi_{n}({\mathcal{C}})}{2}\nu_{n}\;,$$

(12)

where $n$ labels eigenvalues of $B$.

Eq. (12) is explained as follows. Suppose $V_{mn}$ diagonalize $B$ such that $B=V\,\beta\,V^{\dagger}$. Suppose the $n$th eigenvalue of $B$, denoted as $\beta_{n}$, has a $2\pi\nu_{n}$ jump across a vortex loop $\mathcal{C}$, namely $\beta_{n+}-\beta_{n-}=2\pi\nu_{n}$, where $+$,$-$ labels two sides of the vortex loop and $\beta_{n+}$ and $\beta_{n-}$ are the limit of $\beta_{n}$ from the corresponding sides. For the case of single band 3D Chern insulators, we will shown immediately that one can always make a gauge choice such that at one of the eigenvalues $\beta_{n}$ jumps by $2\pi C$ across an exactly straight line alone $k_{z}$ direction, as schematically shown in Fig. 1(a). In other words, in such a case, $\mathcal{C}$ is a straight line extending from $k_{z}=0$ to $k_{z}=1$ at some given $k_{x}$. Then, $\phi_{n}(\mathcal{C})$ and $\xi_{n}(\mathcal{C})$ are two kinds of Berry phases calculated along $\mathcal{C}$: $\phi_{n}(\mathcal{C})$ is the Berry phase of $n$th eigenvector (of $B$) $\mathbf{v}_{n}\equiv(V_{1n},V_{2n},\cdots,V_{Mn})^{T}$; while $\xi_{n}(\mathcal{C})$ is the Berry phase of a unitarily transformed Bloch state $\sum_{m=1}^{M}V_{mn}\ket{u_{m}(k_{x},\lambda=0,k_{z})}$.

We now start deriving the Chern-Simons OME coupling of a 3D Chern insulator with Chern vector $(0,0,C)$. We note that this is a general situation when discussing bulk properties, since we can always re-define our coordinate system such that the Chern vector $(C_{1},C_{2},C_{3})$ is vanishing in $z$-$x$ and $y$-$z$ plane  [20]. We choose the aforementioned gauge with discontinuity at $k_{y}=0,1$ plane. We further enforce an additional requirement that in $k_{z}$, we adopt a periodic parallel-transport gauge (see Appendix $A$), which corresponds to maximally localized hybrid Wannier representation after Fourier transformation with respect to $k_{z}$  [21].

It is straightforward to attain our aforementioned 3D gauge by first enforcing a periodic parallel-transport gauge in the $k_{x}$, then a smooth gauge obtained from parallel transportation in $k_{y}$ without enforcing periodicity condition, and finally a periodic parallel-transport gauge in $k_{z}$ ¹¹1See Appendix A for detailed explanations of the terminologies such as ‘periodic parallel-transport gauge’ and ‘parallel transportation’ etc.. The procedures for constructing such a 3D gauge is schematically illustrated in Fig. 1(a) and Fig. 1(b). Such a gauge construction would push all the potential discontinuity in $k_{y}$ to the $k_{y}=1$ plane. Note that enforcing a periodic parallel-transport gauge in $k_{z}$ will not spoil the smoothness of gauge in $k_{x}$ and $k_{y}$ because there is no winding of Berry phase in $k_{x}$-$k_{z}$ nor $k_{y}$-$k_{z}$ plane due to the corresponding vanishing Chern numbers, so that the gauge transformation itself is smooth in $\mathbf{k}$.

For simplicity, we will first calculate $\theta$ in the case of a single occupied band, where the essential properties of the topological OME response in 3D Chern insulators are already clearly manifested. In the single-band case, $B$ is reduced to a real number $\beta$, and we note that in principle $\beta$ should depend on both $k_{z}$ and $k_{x}$. However, we will prove in the following that in the single occupied band case, by insisting a periodic parallel-transport gauge to the Bloch states in the $k_{z}$ direction, $\beta$ is only dependent on $k_{x}$. Under such a “periodic parallel transport gauge” in the $k_{z}$ direction, the Fourier transformation of the Bloch states with respect to $k_{z}$ would give rise to the maximally localized Wannier function in $z$ direction.

To see this more explicitly, we first write down the most general gauge-discontinuity condition at the $k_{y}=0,1$ plane

$$\ket{u(k_{x},k_{y}=1,k_{z})}=e^{i\beta(k_{x},k_{z})}e^{-i2\pi\hat{y}}\ket{u_{n}(k_{x},k_{y}=0,k_{z})}\;,$$

(13)

such that $\beta$ depends on both $k_{x}$ and $k_{z}$.

Then, we will show that $\partial_{z}\beta=0$ under parallel-transport gauge in $k_{z}$ direction. Let us define the hybrid Wannier function that is exponentially localized at layer $l$ as

$$\ket{W_{ln}(k_{x},k_{y})}=\int_{0}^{1}dk_{z}\,|u_{n}(k_{x},k_{y},k_{z})\rangle e^{i2\pi k_{z}(\hat{z}-l)}$$

(14)

where $l$ is the layer index (or unit-cell index in $z$ direction). We have chosen the normalization convention that the periodic part of the Bloch function $\ket{u}$ is normalized within the unit cell. For the hybrid Wannier function $\ket{W_{l}(k_{x},k_{y})}$, it is normalized within the unit cell for the $x$, $y$ coordinates, while normalized over the entire space for the $z$ coordinate.

The inverse transformation of Eq. (14) is

$$\ket{u(k_{x},k_{y},k_{z})}=\sum_{l}e^{i2\pi k_{z}(l-\hat{z})}\ket{W_{l}(k_{x},k_{y})}$$

(15)

Note again, that $k_{z}\in[0,1)$ is defined in the reduced coordinate. With periodic parallel-transport gauge in $k_{z}$, $\ket{W_{l}(k_{x},k_{y})}$ is maximally localized in the $z$ direction, which is the eigenstate of $P\,\hat{z}\,P$ operator ($P$ denoting projection operator to occupied subspace), i.e.,

$$\braket{W_{0}(k_{x},k_{y})}{\,\hat{z}\,}{W_{l}(k_{x},k_{y})}=\bar{z}_{0}(k_{x},k_{y})\delta_{0,l}$$

(16)

where $\bar{z}_{0}(k_{x},k_{y})$ is the Wannier center in $z$ direction, which in principle should be dependent on $k_{x}$ and $k_{y}$, and may reflect the nontrivial band topology as in the case of three-dimensional topological insulators. However, for weakly coupled Chern-insulator layers, one would expect to get weak $k_{x}$-$k_{y}$ dispersions of $\bar{z}_{0}(k_{x},k_{y})$. Multiplying by $e^{i2\pi k_{z}^{\prime}l}$ and summing over $l$, we have

$$A_{z}(k_{x},k_{y},k_{z})=2\pi\bar{z}_{0}(k_{x},k_{y})$$

(17)

so that

	$\displaystyle\bar{z}_{0}(k_{x},0)=\frac{i}{2\pi}\braket{u(k_{x},0,k_{z})}{\partial_{z}u(k_{x},0,k_{z})}$		(18a)
	$\displaystyle\bar{z}_{0}(k_{x},1)=\frac{i}{2\pi}\braket{u(k_{x},1,k_{z})}{\partial_{z}u(k_{x},1,k_{z})}$		(18b)

Combining the fact that $\bar{z}_{0}(k_{x},0)=\bar{z}_{0}(k_{x},1)$ and Eq. (15), Eqs. (18b), it follows straightforwardly that

$$\partial_{z}\beta(k_{x},k_{z})=0\;.$$

(19)

This completes our proof that in the single occupied band case under periodic parallel-transport gauge choice in the $k_{z}$ direction, a smooth but non-periodic gauge in the $k_{y}$ direction Eq. (13) is always allowed, with the phase factor $\beta$ (characterizing the gauge-discontinuity at the $k_{y}=0,1$ plane) being independent of $k_{z}$.

Additionally, the nonzero Chern number $C$ within each ($k_{x}$,$k_{y}$) plane and periodicity in $k_{x}$ requires that

$$\beta(1)-\beta(0)=2\pi C$$

(20)

if one makes a proper branch choice such that $\beta$ remains continuous from $k_{x}=0$ to $1$. If a different branch choice is made, there can be a $2\pi C$ jump in $\beta$ at some $k_{x}^{*}$, which creates a “vortex line” or “vortex loop” as explained later. To see this, one only needs to invoke Stokes’ theorem, i.e

$$\int_{0}^{1}dk_{x}{A_{x}(k_{y}=0)}-\int_{0}^{1}dk_{x}{A_{x}(k_{y}=1)}=2\pi C$$

(21)

where we have used the fact that $A_{y}$ is smooth and periodic along $k_{x}$ so that the contributions from the other edges cancel. Then Eq. (20) follows immediately by combining Eq. (15) and Eq. (21).

To put the equation in a more illuminating form, we cast the Berry connection and Berry curvature into the hybrid Wannier representation, the definition of which is given by Eqs. (14)-(15). The Berry connections can be expressed in hybrid-Wannier basis as follows:

$$A_{a}=\sum_{l}e^{i2\pi k_{z}l}A_{a,0l}$$

(22)

where

$$A_{a,0l}=\sum_{l}e^{i2\pi k_{z}l}\,i\braket{W_{0}(k_{x},k_{y})}{\partial_{a}W_{l}(k_{x},k_{y})}\\$$

(23)

where $a=x,y$. Here we use the convention that $\ket{u_{\mathbf{k}}}$ is normalized to 1 over a unit cell and $\braket{u_{\mathbf{k}}}{\partial_{x}u_{\mathbf{k}}}$ is carried out over a unit cell as well. And,

$$A_{z}=2\pi\bar{z}_{0}$$

(24)

where $\bar{z_{l}}(k_{x},k_{y})=\braket{W_{l}(k_{x},k_{y})}{\hat{z}}{W_{l}(k_{x},k_{y})}$ is the hybrid Wannier center in the reduced coordinate. It follows from Eq. (22) and (24) that

	$\displaystyle\Omega_{xy}$	$\displaystyle=\sum_{l}e^{i2\pi k_{z}l}\Omega_{xy,0l}$		(25a)
	$\displaystyle\Omega_{yz}$	$\displaystyle=2\pi\partial_{y}\bar{z}_{0}-i2\pi\sum_{l}le^{i2\pi k_{z}l}A_{y,0l}$		(25b)
	$\displaystyle\Omega_{zx}$	$\displaystyle=i2\pi\sum_{l}le^{i2\pi k_{z}l}A_{x,0l}-2\pi\partial_{x}\bar{z}_{0}$		(25c)
	$\displaystyle A_{x}\Omega_{yz}$	$\displaystyle=2\pi\sum_{l}e^{i2\pi k_{z}l}A_{x,0l}\partial_{y}\bar{z}_{0}$
		$\displaystyle\quad-i2\pi\sum_{l^{\prime},l}e^{i2\pi k_{z}(l^{\prime}+l)}lA_{x,0l^{\prime}}A_{y,0l}$		(25d)
	$\displaystyle A_{y}\Omega_{zx}$	$\displaystyle=i2\pi\sum_{l^{\prime},l}e^{i2\pi k_{z}(l^{\prime}+l)}lA_{y,0l^{\prime}}A_{x,0l}$
		$\displaystyle\quad-2\pi\sum_{l}e^{i2\pi k_{z}l}A_{y,0l}\partial_{x}\bar{z}_{0}$		(25e)
	$\displaystyle A_{z}\Omega_{xy}$	$\displaystyle=2\pi\sum_{l}\bar{z}_{0}e^{i2\pi k_{z}l}\Omega_{xy,0l}$		(25f)

where $\Omega_{xy,0l}=\partial_{x}A_{y,0l}-\partial_{y}A_{x,0l}$. Using Eq. (11a), we put everything together and integrate out the $k_{z}$ dependence, we arrive at

$\displaystyle\theta_{\rm{BK}}=$

$\displaystyle-\frac{1}{2}\int dk_{x}\,dk_{y}\,\Big{(}\bar{z}_{0}\,\Omega_{xy,00}+A_{x,00}\partial_{y}\bar{z}_{0}-A_{y,00}\partial_{x}\bar{z}_{0}\big{)}$

(26)

where we have used

$$\sum_{l}l(A_{x,l0}A_{y,0l}-A_{y,l0}A_{x,0l})=2\sum_{l}l(A_{x,l0}A_{y,0l})=0$$

(27)

In the interlayer decoupled limit, $\bar{z}_{0}$ is independent of the $(k_{x},k_{y})$, so the second and the third terms are exactly zeros. For reasons we will see, we define

$\displaystyle\theta^{\text{anom.}}_{\text{BK}}=$

$\displaystyle-\frac{1}{2}\int dk_{x}\,dk_{y}\,\bar{z}_{0}\,\Omega_{xy,00}$

(28)

For $\theta_{\text{GD}}$, we start with Eq. (11b) and the corresponding gauge choice within such a fictitious space are given in Eqs. (7). As a reminder, we choose a gauge inside the $\lambda$-space that “linearly interpolates” the gauge discontinuity, i.e,

$$\ket{u(k_{x},k_{z},\lambda)}=e^{-i\lambda\beta(k_{x})}\ket{u(k_{x},k_{z},\lambda=0)}$$

(29)

with $\lambda\in{[0,1)}$. We note again that $\ket{u(k_{x},k_{z},\lambda=0}=\ket{u(k_{x},k_{y}=1,k_{z})}$ such that the $\lambda$ linearly interpolates the gauge from $k_{y}=1$ to that of $k_{y}=0$. Define as usual

$$\ket{W_{l}(k_{x},\lambda)}=\int_{0}^{1}dk_{z}\,e^{i2\pi k_{z}(\hat{z}-l)}\ket{u(k_{x},\lambda)}$$

(30)

and

$$\ket{u(k_{x},k_{z},\lambda)}=\sum_{l}e^{i2\pi k_{z}(l-\hat{z})}\ket{W_{l}(k_{x},\lambda)}$$

(31)

We find the Berry connections can be expressed as

	$\displaystyle A_{x}(\lambda)=A^{0}_{x}+\lambda\,\partial_{x}\beta\;$		(32a)
	$\displaystyle A_{z}=2\pi\bar{z}_{0}(k_{x})\;$		(32b)
	$\displaystyle A_{\lambda}=\beta(k_{x})$		(32c)

where $A_{x}^{0}\equiv A_{x}(\lambda=0)$ and $\bar{z}_{0}(k_{x})\equiv\bar{z}_{0}(k_{x},\lambda)=\langle W_{0}(k_{x},\lambda)|\hat{z}|W_{0}(k_{x},\lambda)\rangle$ is independent of $\lambda$. The Berry curvatures read

	$\displaystyle\Omega_{\lambda z}=0\;$		(33a)
	$\displaystyle\Omega_{x\lambda}=0\;$		(33b)
	$\displaystyle\Omega_{zx}=\partial_{z}A^{0}_{x}(k_{x},k_{z})-2\pi\partial_{x}\bar{z}_{0}(k_{x})$		(33c)

Finally, plugging the above equations into Eq. (11b) and taking use of the periodicity in $k_{z}$, we obtain

$$\theta_{\text{GD}}=\frac{1}{2}\int_{0}^{1}dk_{x}\beta\partial_{x}\bar{z}_{0}$$

(34)

For reasons we will soon see, using integration by part, we write it as

	$\displaystyle\theta_{\text{GD}}=$	$\displaystyle\frac{1}{2}\int_{0}^{1}dk_{x}\,\beta\partial_{x}\bar{z}_{0}$
	$\displaystyle=$	$\displaystyle\frac{1}{2}(\beta\bar{z}_{0})\big{|}^{1}_{k_{x}=0}-\frac{1}{2}\int_{0}^{1}\,dk_{x}\,\bar{z}_{0}\partial_{x}\beta$		(35)

where

	$\displaystyle(\beta\bar{z}_{0})\big{|}^{1}_{k_{x}=0}=$	$\displaystyle\beta(1)\bar{z}_{0}(k_{x}=1)-\beta(0)\bar{z}_{0}(k_{x}=0)$
	$\displaystyle=$	$\displaystyle(\beta(0)+2\pi C)\bar{z}_{0}(k_{x}=0)-\beta(0)\bar{z}_{0}(k_{x}=0)$
	$\displaystyle=$	$\displaystyle 2\pi C\bar{z}_{0}(k_{x}=0)$		(36)

For the second equality, we have used the fact that the winding of $\beta$ along $k_{x}$ direction is just $2\pi$ times the Chern number $C$ of the $(k_{x},k_{y})$ plane and that the $(k_{z},k_{x})$ plane has Chern number 0 such that $\bar{z}_{0}(k_{x}=1)=\bar{z}_{0}(k_{x}=0)$. Therefore,

$$\theta_{\text{GD}}=\pi C\bar{z}_{0}(k_{x}=0)-\frac{1}{2}\int_{0}^{1}\,dk_{x}\,\bar{z}_{0}\partial_{x}\beta(k_{x})$$

(37)

For $\theta_{\text{VL}}$, since our “one dimensional matrix” $\beta$ is always diagonalized, we are free to choose the “eigenvector” of $\beta$ (which only has one element), $V_{11}=1$, which certainly has zero Berry phase. This eliminates the first term in Eq. (12), leaving us with

$$\theta_{\text{VL}}=-\frac{1}{2}C\int_{0}^{1}A^{0}_{z}(k_{x}=0)dk_{z}=-\pi C\bar{z}_{0}(k_{x}=0)\;,$$

(38)

where we have used $\bar{z}_{0}=\frac{1}{2\pi}\xi$, with $\xi$ denoting the Berry phase of the physical Bloch state winding around the vortex loop, as explained near Eq. (12). The minus sign comes from the fact that the discontinuity of $\beta$ is manifested as a local $-2\pi C$ shift across the vortex loop (or vortex line) at $k_{x}^{*}$ (see Fig. 1(c)). Combining the two contributions, we arrive at

$$\theta_{\text{GD}}+\theta_{\text{VL}}=-\frac{1}{2}\int_{0}^{1}\,dk_{x}\,\bar{z}_{0}(k_{x}=0,\lambda=0)\partial_{x}\beta\equiv\theta^{\text{anom.}}_{\text{GDVL}}$$

(39)

We see now the boundary term for $\theta_{\text{GD}}$ is canceled by a term in $\theta_{\text{VL}}$. Here we call this term the “anomalous” Chern-Simons term, denoted by $\theta^{\text{anom.}}_{\text{GDVL}}$. The name “anomalous” will be justified later.

With these extensions, Eq. (11b) becomes

	$\displaystyle\theta_{\text{BK}}=$	$\displaystyle-\frac{1}{2}\,\int dk_{x}dk_{y}\Big{(}\,\sum_{n}z_{0n}{\Omega_{xy,0n,0n}}\;$
		$\displaystyle+\sum_{m}(A_{x,0n,0n}\partial_{y}\bar{z}_{0n}-A_{y,0n,0n}\partial_{x}\bar{z}_{0n})\;$
		$\displaystyle+2\sum_{lmn}ilA_{x,lm,0n}A_{y,0n,lm}\;$
		$\displaystyle-2i\sum_{lmn}(\bar{z}_{0m}-\bar{z}_{0n})A_{x,lm,0n}A_{y,0n,lm}\,\Big{)}$		(42)

where

	$$A_{a,0n,lm}=i\langle W_{0,n}(k_{x},k_{y})|\partial_{a}W_{l,m}(k_{x},k_{y})\rangle$$		(43a)
	$$\Omega_{xy,0n,lm}=\partial_{x}A_{y,0n,lm}-\partial_{y}A_{x,0n,lm}$$		(43b)

Here, $\ket{W_{l,m}(k_{x},k_{y})}$ denotes the $m$th hybrid Wannier function localized in layer $l$ and $\bar{z}_{0n}=\langle W_{0,n}(k_{x},k_{y})|\hat{z}|W_{0,n}(k_{x},k_{y})\rangle$. We note the result resembles that presented in Ref. 19, the main difference being the extra $\frac{1}{2}$ on the term related to $\Omega_{xy,00,00}$ and $z_{00}$. We will see that we will be saved by the discontinuity contributing a similar term. Even in the multi-band cases, with non-vanishing interlayer coupling, the bulk term may vanish due to constraint from crystalline symmetries such as inversion symmetry. In the presence of inversion symmetry, it is straightforward to show that

	$\displaystyle A_{a}(\mathbf{k})=-A_{a}(-\mathbf{k})\;$		(44a)
	$\displaystyle\Omega_{ab}(\mathbf{k})=\Omega_{ab}(-\mathbf{k})$		(44b)

Then, clearly Eq. (42) vanishes up to $\pi C$ (see Sec. V).

Following Eq. (11b) and Eq. (12), we are able to find the gauge-discontinuity and vortex-loop contributions can be expressed as

	$\displaystyle\theta_{\textrm{GD}}+\theta_{\textrm{VL}}=$	$\displaystyle-\frac{1}{4\pi}\int dk_{x}d\lambda dk_{z}\,\mathrm{Tr}\,(-iB[\partial_{z}M,\partial_{x}M^{\dagger}]$
		$\displaystyle+iM[\partial_{z}B,\partial_{x}M^{\dagger}]+iM[\partial_{z}M^{\dagger},\partial_{x}B]$
		$\displaystyle+2B[\partial_{z}M^{\dagger},MA_{x}^{0}]+\partial_{z}B[M^{\dagger},MA_{x}^{0}]$
		$\displaystyle+B\partial_{z}A_{x}^{0}-2B[\partial_{x}M^{\dagger},MA_{z}^{0}]$
		$\displaystyle-\partial_{x}B[M^{\dagger},MA_{z}^{0}])$
		$\displaystyle+\int_{0}^{1}dk_{z}\,\sum_{mn}V^{\dagger}_{mn}\partial_{z}V_{nm}C_{n}+\theta^{\text{anom.}}_{\text{GDVL}}$		(45)

where

$$\theta^{\text{anom.}}_{\text{GDVL}}=-\frac{1}{4\pi}\int dk_{z}dk_{z}\,\mathrm{Tr}(\partial_{x}B\,A^{0}_{z})$$

(46)