Efficient computation of the $W_3$ topological invariant and application to Floquet-Bloch systems

TL;DR

This paper presents a fast, robust algorithm for computing the $W_3$ topological invariant in unitary maps, applicable to Floquet-Bloch systems, and demonstrates its effectiveness on a driven graphene model to predict edge states.

Contribution

The paper introduces a new efficient algorithm for $W_3$ invariant computation that requires minimal manipulation and applies it to Floquet-Bloch systems, enabling accurate topological analysis.

Findings

Algorithm converges rapidly on coarse grids.

Successfully applied to irradiated graphene model.

Predicts the number of anomalous edge states.

Abstract

We introduce an efficient algorithm for the computation of the $W_3$ invariant of general unitary maps, which converges rapidly even on coarse discretization grids. The algorithm does not require extensive manipulation of the unitary maps, identification of the precise positions of degeneracy points, or fixing the gauge of eigenvectors. After construction of the general algorithm, we explain its application to the $2+1$ dimensional maps that arise in the Floquet-Bloch theory of periodically driven two-dimensional quantum systems. We demonstrate this application by computing the $W_3$ invariant for an irradiated graphene model with a continuously modulated Hamilton operator, where it predicts the number of anomalous edge states in each gap.

Tables

Table 1. Table 1: Vertices 𝐩 1 , … , 𝐩 4 subscript 𝐩 1 … subscript 𝐩 4 \mathbf{p}_{1},\dots,\mathbf{p}_{4} of the faces 𝖥 𝐩 , α subscript 𝖥 𝐩 𝛼 \mathsf{F}_{\mathbf{p},\alpha} used in the algorithm.

	${𝐩}_1$	${𝐩}_2$	${𝐩}_3$	${𝐩}_4$
${𝖥}_{𝐩,1}$:	$𝐩$	$𝐩+{𝜹}_2$	$𝐩+{𝜹}_2+{𝜹}_3$	$𝐩+{𝜹}_3$
${𝖥}_{𝐩,2}$:	$𝐩$	$𝐩+{𝜹}_3$	$𝐩+{𝜹}_1+{𝜹}_3$	$𝐩+{𝜹}_1$
${𝖥}_{𝐩,3}$:	$𝐩$	$𝐩+{𝜹}_1$	$𝐩+{𝜹}_1+{𝜹}_2$	$𝐩+{𝜹}_2$

Equations

$$U(0,{\mu }_2,{\mu }_3)=U(1,{\mu }_2,{\mu }_3),U({\mu }_1,0,{\mu }_3)=U({\mu }_1,1,{\mu }_3),U({\mu }_1,{\mu }_2,0)=U({\mu }_1,{\mu }_2,1).$$

$$W_{3}[U]=\frac{1}{24\pi^{2}}\int\limits_{[0,1)^{3}}\epsilon_{\alpha\beta\gamma}\operatorname{tr}\big{[}(U^{-1}\partial^{\alpha}U)\,(U^{-1}\partial^{\beta}U)\,(U^{-1}\partial^{\gamma}U)\big{]}\,\mathrm{d}^{3}\mu\;.$$

$$W_1^{\alpha }[U]=\frac{1}{2\pi \mathrm{i}}\underset{0}{\overset{1}{\int }}\mathrm{tr}[U^{-1}{\partial }^{\alpha }U]\mathrm{d}{\mu }_{\alpha }$$

$$U(\mu )=S(\mu )D(\mu )S^{†}(\mu ),$$

$$U^{-1}\partial^{\alpha}U=S\big{(}Z^{\alpha}+Y^{\alpha}-X^{\alpha}\big{)}S^{\dagger}\;,$$

$$X^{\alpha }=S^{†}{\partial }^{\alpha }S,Y^{\alpha }=D^{†}{\partial }^{\alpha }D=({\partial }^{\alpha }D)D^{†},Z^{\alpha }=D^{†}{X}^{\alpha }D.$$

$$\operatorname{tr}\big{[}(U^{-1}\partial^{\alpha}U)\,(U^{-1}\partial^{\beta}U)\,(U^{-1}\partial^{\gamma}U)\big{]}=\operatorname{tr}\big{[}(Z^{\alpha}+Y^{\alpha}-X^{\alpha})(Z^{\beta}+Y^{\beta}-X^{\beta})(Z^{\gamma}+Y^{\gamma}-X^{\gamma})\big{]}\;,$$

$$\begin{split}&\epsilon_{\alpha\beta\gamma}\operatorname{tr}\big{[}(U^{-1}\partial^{\alpha}U)\,(U^{-1}\partial^{\beta}U)\,(U^{-1}\partial^{\gamma}U)\big{]}=\\ &\epsilon_{\alpha\beta\gamma}\operatorname{tr}\big{[}-X^{\alpha}X^{\beta}X^{\gamma}+3X^{\alpha}X^{\beta}Y^{\gamma}+3X^{\alpha}X^{\beta}Z^{\gamma}-3X^{\alpha}Y^{\beta}Z^{\gamma}-3X^{\alpha}Z^{\beta}Y^{\gamma}\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad-3X^{\alpha}Z^{\beta}Z^{\gamma}+3Y^{\alpha}Z^{\beta}Z^{\gamma}+Z^{\alpha}Z^{\beta}Z^{\gamma}\big{]}\end{split}$$

$${\partial }^{\alpha }{S}^{†}=-S^{†}({\partial }^{\alpha }S)S^{†},{\partial }^{\alpha }{D}^{†}=-D^{†}({\partial }^{\alpha }D)D^{†}$$

$${\partial }^{\alpha }{X}^{\beta }$$

$${\partial }^{\alpha }{Y}^{\beta }$$

$${\partial }^{\alpha }{Z}^{\beta }$$

$$\displaystyle\partial^{\alpha}\big{(}Z^{\beta}X^{\gamma}\big{)}$$

$$\displaystyle\partial^{\alpha}\big{(}Z^{\beta}Y^{\gamma}\big{)}$$

$$\epsilon_{\alpha\beta\gamma}\operatorname{tr}\big{[}(U^{-1}\partial^{\alpha}U)\,(U^{-1}\partial^{\beta}U)\,(U^{-1}\partial^{\gamma}U)\big{]}=\epsilon_{\alpha\beta\gamma}\Big{(}3\operatorname{tr}\big{[}\partial^{\alpha}\big{(}Z^{\beta}X^{\gamma}\big{)}\big{]}-6\operatorname{tr}\big{[}\partial^{\alpha}\big{(}Z^{\beta}Y^{\gamma}\big{)}\big{]}\Big{)}\;.$$

$$\tilde{S}(\mu )=S(\mu )G(\mu ),$$

$${\tilde{X}}^{\alpha }=G^{†}{X}^{\alpha }G+G^{\alpha },{\tilde{Z}}^{\alpha }=G^{†}{Z}^{\alpha }G+G^{\alpha },$$

$$\operatorname{tr}\big{[}Z^{\beta}X^{\gamma}\big{]}\big{|}_{\textstyle(1,\mu_{2},\mu_{3})}=\operatorname{tr}\big{[}Z^{\beta}X^{\gamma}+Z^{\beta}G^{\gamma}+G^{\beta}X^{\gamma}+G^{\beta}G^{\gamma}\big{]}\big{|}_{\textstyle(0,\mu_{2},\mu_{3})}\;,$$

$$W_{3}[U]=\frac{1}{2\pi\mathrm{i}}\int\limits_{[0,1)^{3}}\operatorname{tr}\big{[}F_{\alpha}Y^{\alpha}\big{]}\,\mathrm{d}^{3}\mu\;,$$

$$F_{\alpha }=\frac{1}{2\pi \mathrm{i}}{ϵ}_{\alpha \beta \gamma }({\partial }^{\beta }{X}^{\gamma }).$$

$$W^{\alpha ,\nu }=\frac{1}{2\pi \mathrm{i}}\underset{0}{\overset{1}{\int }}{Y}^{\alpha ,\nu }\mathrm{d}{\mu }_{\alpha },$$

$$W_1^{\alpha }[U]=\sum _{\nu =1}^n{W}^{\alpha ,\nu }.$$

$$C_{\alpha }^{\nu }=\underset{0}{\overset{1}{\int }}\underset{0}{\overset{1}{\int }}{F}_{\alpha }^{\nu }(\mu )\mathrm{d}{\mu }_{\alpha -1}\mathrm{d}{\mu }_{\alpha +1}$$

$$W_{3}[U]=\frac{1}{2\pi}\left(-\int\limits_{[0,1)^{3}}\operatorname{tr}\big{[}(\partial^{\alpha}F_{\alpha})\Phi\big{]}\,\mathrm{d}^{3}\mu\,+\,\sum_{\alpha=1}^{3}\int\limits_{0}^{1}\!\!\!\int\limits_{0}^{1}\operatorname{tr}\big{[}F_{\alpha}\Phi\big{]}\Big{|}_{\mu_{\alpha}=0}^{\mu_{\alpha}=1}\,\mathrm{d}\mu_{\alpha-1}\,\mathrm{d}\mu_{\alpha+1}\right)\;.$$

$${\hat{F}}_{\mathbf{p},\alpha }^{\nu }=\frac{1}{2\pi \mathrm{i}}\log {\mathcal{U}}^{\nu }({\mathbf{p}}_1,{\mathbf{p}}_2){\mathcal{U}}^{\nu }({\mathbf{p}}_2,{\mathbf{p}}_3){\mathcal{U}}^{\nu }({\mathbf{p}}_3,{\mathbf{p}}_4){\mathcal{U}}^{\nu }({\mathbf{p}}_4,{\mathbf{p}}_1),$$

$${\mathcal{U}}^{\nu }({\mathbf{p}}_i,{\mathbf{p}}_j)=\frac{⟨{\mathbf{s}}^{\nu }({\mathbf{p}}_i),{\mathbf{s}}^{\nu }({\mathbf{p}}_j)⟩}{|⟨{\mathbf{s}}^{\nu }({\mathbf{p}}_i),{\mathbf{s}}^{\nu }({\mathbf{p}}_j)⟩|}$$

$${\hat{C}}_{\mathbf{p}}^{\nu }=\sum _{\alpha =1}^3{\hat{F}}_{\mathbf{p}+{\mathbf{\delta }}_{\alpha },\alpha }^{\nu }-{\hat{F}}_{\mathbf{p},\alpha }^{\nu }.$$

$${\hat{W}}_3=\sum _{i_1,i_2,i_3=1}^N\sum _{\nu =1}^n\left({\hat{C}}_{\mathbf{p}}^{\nu }{M}_{\mathbf{p}}^{\nu }+\sum _{\alpha =1}^3{\hat{F}}_{\mathbf{p},\alpha }^{\nu }{m}_{\mathbf{p},\alpha }^{\nu }\right)$$

$${\varphi }_{\mathbf{p}}^{\nu }=-\mathrm{i}\log d^{\nu }(\mathbf{p})$$

$${\hat{W}}_1^{\alpha }=\sum _{\nu =1}^n\sum _{i_{\alpha }=1}^N{m}_{\mathbf{p},\alpha }^{\nu }$$

$${\hat{C}}_{\alpha }^{\nu }=\sum _{i_{\alpha -1},i_{\alpha +1}=1}^N{\hat{F}}_{\mathbf{p},\alpha }^{\nu }(\text{with }{i}_{\alpha }\text{ fixed})$$

$$|{\Phi }^{\nu }(\mathbf{p}+{\mathbf{\delta }}_{\alpha })-{\Phi }^{\nu }(\mathbf{p})|<\pi \text{ for }1\le i_1,i_2,i_3<N$$

Full text

††journal: J. Phys. A

1 Introduction

Topological invariants have gained considerable interest in solid state physics and related fields through their application to the quantum Hall effect [1, 2], topological insulators [3, 4, 5, 6], and more recently, Floquet-Bloch systems [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. While the relevant topological invariants for the quantum Hall effect are the Chern numbers of the respective bands of the Hamiltonian, the invariants for Floquet-Bloch systems are constructed for unitary maps that are derived from the time propagator of the periodically driven system. These invariants, which can be defined generally for unitary maps in odd dimensions [20, 21, 22, 23, 24, 25], can be understood as the generalization of the winding number of a circle map $\mathbb{S}^{1}\to\mathbb{S}^{1}$. For Floquet-Bloch systems in two spatial and one time dimension, the $W_{3}$ invariant is of primary interest. In particular, it gives the number of (anomalous) edge states that can not be predicted with other invariants such as the simpler $W_{1}$ invariant, which is directly related to the winding number of the eigenvalues of a unitary map, or the Chern numbers of the Floquet bands.

In this paper we develop an algorithm for the efficient computation of the $W_{3}$ invariant of general unitary maps, and then provide the application to Floquet-Bloch systems. Our algorithm is motivated by the algorithm of Fukui, Hatsugai, Suzuki [26], which allows for the direct gauge-invariant computation of Chern numbers. However, while the Chern numbers depend only on the eigenvectors of the respective bands, we must consider both the variation of eigenvalues and eigenvectors of the unitary map. Furthermore, we cannot assume a gap condition: A non-zero $W_{3}$ invariant implies the existence of degeneracy points of the unitary map. Therefore, construction of our algorithm requires a few additional considerations.

The paper is organized as follows. In Sec. 2 we specify the mathematical setting for this paper and fix our notation. We start our construction of the algorithm in Sec. 3 with the derivation of a gauge-invariant expression for the $W_{3}$ invariant that involves the eigenvalue and eigenvector bands of the unitary map, before we explain the algorithmic scheme in Sec. 4. In Sec. 5 we address the application to Floquet-Bloch systems introduced in Ref. [12], which involves only a small modification of the general computational scheme. The actual application to the anomalous edge states in an irradiated graphene model is given in Sec. 6. We conclude in Sec. 7.

2 The $W_{3}$ invariant

The mathematical setting for this paper is as follows. Let $U(\cdot):\mathbb{T}^{3}\to\mathbb{U}(n)$ denote a smooth map from the three-dimensional torus $\mathbb{T}^{3}=\mathbb{S}^{1}\times\mathbb{S}^{1}\times\mathbb{S}^{1}$ to the Lie group $\mathbb{U}(n)$ of unitary $n\times n$ matrices. We can view $U(\cdot)$ as a map defined on the cube $[0,1]^{3}$ that is periodic in every argument, as in

[TABLE]

In this situation, the quantity $W_{3}[U]$ is defined as [11, 12]

[TABLE]

Here, $\epsilon_{\alpha\beta\gamma}$ denotes the antisymmetric Levi-Civita symbol. We write $\partial^{\alpha}A(\mu)=\partial A(\mu)/\partial\mu_{\alpha}$, with $\alpha=1,2,3$, for the derivative of a matrix-valued function $A(\cdot)$. The indices $\alpha,\beta,\gamma$ are counted modulo 3. For example, if one coordinate $\mu_{\alpha}$ is given, the remaining two coordinates are $\mu_{\alpha-1}$ and $\mu_{\alpha+1}$. Summation over repeated indices $\alpha,\beta,\gamma$ is implied, and we often drop the argument $\mu$ to lighten the notation.

$W_{3}[U]$ is an integer that is invariant under continuous deformations of $U(\cdot)$, i.e., it is a topological invariant. Note that $W_{3}[U]=0$ if $U(\cdot)$ is constant along one direction $\mu_{\alpha}$. In addition to $W_{3}[U]$ we have three invariants

[TABLE]

along each direction $\alpha=1,2,3$. These invariants are directly related to the winding numbers of the eigenvalues of $U(\cdot)$, while $W_{3}[U]$ also involves the eigenvectors as we will see below. Note that, for a smooth map, $W^{\alpha}_{1}[U]$ does not depend on the values of the two coordinates $\mu_{\alpha-1}$, $\mu_{\alpha+1}$ kept fixed in the above expression.

3 Gauge-invariant expression for the $W_{3}$ invariant

Let us assume that $U(\cdot)$ is given by its decomposition into eigenvalues and eigenvectors,

[TABLE]

with a unitary matrix map $S(\cdot)$ and a diagonal matrix map $D(\cdot)$. The columns of $S(\mu)$ give the eigenvectors $\mathbf{s}^{\nu}(\mu)$ of $U(\mu)$, while the corresponding eigenvalues $d^{\nu}(\mu)$ appear as the entries of $D(\mu)$.

Since the above decomposition is not unique, the numbering of the eigenvalues and eigenvectors can affect the analytical properties of the maps $S(\cdot)$, $D(\cdot)$. In contrast to the Hermitian case, there is also no natural ordering of the eigenvalues $d^{\nu}(\mu)\in\mathbb{S}^{1}$. In our algorithm, the indices $\nu$ will be assigned according to the distance of eigenvalues and the overlap of eigenvectors on adjacent discretization points (see below).

Away from degeneracy points of $U(\cdot)$, perturbation theory [27] implies that $S(\cdot)$, $D(\cdot)$ can be constructed as smooth maps, essentially by assigning a fixed index $\nu$ to each of the separated ‘bands’ of $U(\cdot)$. In the vicinity of a degeneracy point of $U(\cdot)$, where two eigenvalues coincide, the situation depends on the dimensionality of the domain of $U(\cdot)$. In one dimension, for $\mu\in\mathbb{R}$, it would be possible to continue the maps $S(\cdot)$, $D(\cdot)$ smoothly through the degeneracy point. In dimensions two or higher, as considered here, only continuity of the eigenvalues (but not smoothness) can be achieved in the general case. The eigenvectors may fail to be even continuous at the degeneracy point, as witnessed by the second example in Sec. 4.2. For a generic map $U(\cdot)$ that is derived from a (possibly time-dependent) Hamiltonian, degeneracies occur at isolated points [28]. Therefore, we will assume that the maps $S(\cdot)$, $D(\cdot)$ are smooth, apart from the isolated degeneracy points of $U(\cdot)$.

From the decomposition (3.1) we get

[TABLE]

with

[TABLE]

Note that $D^{\dagger}$ and $\partial^{\alpha}D$ commute, and that $(X^{\alpha})^{\dagger}=-X^{\alpha}$, $(Y^{\alpha})^{\dagger}=-Y^{\alpha}$, $(Z^{\alpha})^{\dagger}=-Z^{\alpha}$. With these definitions, we get

[TABLE]

and so

[TABLE]

for the integrand of Eq. (2.2). Here, we arrange products of operators in lexicographic order, using the cyclicity of $\epsilon_{\alpha\beta\gamma}$ and of the trace $\operatorname{tr}[\cdot]$. Note that terms with two or more $Y$ factors cancel, because these matrices commute among themselves. Note further that the $D$, $D^{\dagger}$ matrices in terms with $Z$ factors partially cancel such that, e.g., $\operatorname{tr}[Y^{\alpha}Z^{\beta}Z^{\gamma}]=\operatorname{tr}[Y^{\alpha}X^{\beta}X^{\gamma}]$. In particular, the first ($X^{\alpha}X^{\beta}X^{\gamma}$) and last ($Z^{\alpha}Z^{\beta}Z^{\gamma}$) term in Eq. (3.5) cancel.

From the derivatives

[TABLE]

of the unitary matrices $S(\mu)$, $D(\mu)$ we get the derivatives

[TABLE]

of the individual factors in the above equations, with $X^{\alpha\beta}=S^{\dagger}\partial^{\alpha\beta}S$, $Y^{\alpha\beta}=D^{\dagger}\partial^{\alpha\beta}D$, $Z^{\alpha\beta}=D^{\dagger}X^{\alpha\beta}D$, and then the relations

[TABLE]

These relations allow us to express Eq. (3.5) in the form

[TABLE]

This expression suggests that one integration in Eq. (2.2), over $\mathrm{d}\mu_{\alpha}$, can be replaced by the difference of the boundary terms for $\mu_{\alpha}=0,1$, and thus may cancel.

However, we cannot expect that the map $S(\cdot)$ in Eq. (3.1) or the factors $X^{\alpha}$, $Z^{\alpha}$ in Eq. (3.2) are compatible with the boundary conditions of the torus (as in Eq. (2.1)). The reason is that the decomposition (3.1) is not unique but depends on the gauge of eigenvectors. Under a gauge transformation $S(\cdot)\mapsto\tilde{S}(\cdot)$, which is defined by a map of unitary diagonal matrices $G(\cdot)$ as

[TABLE]

the above factors change as

[TABLE]

with $G^{\alpha}=G^{\dagger}\partial^{\alpha}G$. The quantity $Y^{\alpha}$ is gauge-invariant.

We must now compare the boundary values of the two terms in Eq. (3.9), which are related through a gauge transformation. From $U(0,\mu_{2},\mu_{3})=U(1,\mu_{2},\mu_{3})$ we conclude that $S(1,\mu_{2},\mu_{3})=S(0,\mu_{2},\mu_{3})G(\mu_{2},\mu_{3})P_{1}$ and $D(1,\mu_{2},\mu_{3})=P_{1}^{\dagger}D(0,\mu_{2},\mu_{3})P_{1}$, with a diagonal gauge matrix $G(\mu_{2},\mu_{3})$ as before, and a permutation matrix $P_{1}\equiv P_{1}(\mu_{2},\mu_{3})$. Because $P$ is a discrete quantity, it does not depend on $\mu_{2}$, $\mu_{3}$. The permutation matrix cancels in all expression involving a trace, but the gauge $G$ does not.

Through the gauge transformation, the boundary values of the first term in Eq. (3.9) are related by

[TABLE]

and similarly for the other directions. The last three terms cancel when combined with $\epsilon_{\alpha\beta\gamma}$ (note that \operatorname{tr}\big{[}Z^{\beta}G^{\gamma}\big{]}=\operatorname{tr}\big{[}X^{\beta}G^{\gamma}\big{]}), and we see that the first term in Eq. (3.9) does not contribute.

In the second term \operatorname{tr}\big{[}\partial^{\alpha}\big{(}Z^{\beta}Y^{\gamma}\big{)}\big{]}=\operatorname{tr}\big{[}\partial^{\alpha}\big{(}X^{\beta}Y^{\gamma}\big{)}\big{]}, the factor $X^{\beta}$ depends on the gauge of $S(\cdot)$ but $Y^{\alpha}$ does not. The boundary values do not cancel and we have to keep this term. Now, we perform the derivative and recognize that in \partial^{\alpha}\big{(}X^{\beta}Y^{\gamma}\big{)}=(\partial^{\alpha}X^{\beta})Y^{\gamma}+X^{\beta}(-Y^{\alpha}Y^{\beta}+Y^{\alpha\beta}) the second term cancels in the final expression because the $Y$ factors commute. In this way, we obtain the manifestly gauge-invariant expression

[TABLE]

with the Berry curvature matrix

[TABLE]

This matrix is diagonal, and gauge-invariant by construction.

For future reference we note that $\partial^{\alpha}F_{\alpha}=0$, unless two eigenvalues coincide at a degeneracy point of $U(\cdot)$. Furthermore, we always have \operatorname{tr}\big{[}F_{\alpha}\big{]}=0.

3.1 Winding and Chern numbers

We can view the decomposition (3.1) as a way to combine the eigenvalues and eigenvectors of $U(\cdot)$ into bands $\mu\mapsto d^{\nu}(\mu)$ and $\mu\mapsto\mathbf{s}^{\nu}(\mu)$ for each $\nu=1,\dots,n$. The corresponding entries of the diagonal matrices $Y^{\alpha}$, $F_{\alpha}$ in Eq. (3.13) contain information about the winding and Chern numbers of these bands.

The diagonal entries $Y^{\alpha,\nu}$ of $Y^{\alpha}$ give the angular velocity $-\mathrm{i}Y^{\alpha,\nu}(\mu)$ of the $\nu$-th eigenvalue, as it moves on the circle $\mathbb{S}^{1}\subset\mathbb{C}$. Therefore, the winding number of band $\nu$ along direction $\alpha$ is

[TABLE]

and the $W_{1}$ invariants of $U$ are given by

[TABLE]

$W^{\alpha,\nu}$ and $W^{\alpha}_{1}[U]$ are both integers. Since the eigenvalues are smooth functions of $\mu$, the value of $W^{\alpha,\nu}$ does not depend on the two coordinates $\mu_{\alpha-1}$, $\mu_{\alpha+1}$ kept fixed in Eq. (3.15).

The diagonal entries $F^{\nu}_{\alpha}(\mu)$ of $F_{\alpha}$ give the Berry curvature of the $\nu$-th eigenvector. Therefore, the integral

[TABLE]

gives the Chern number of band $\nu$ on the surface $\mu_{\alpha}=\mathrm{const.}$, perpendicular to direction $\alpha$. The value of $C^{\nu}_{\alpha}$ can change as a function of $\mu_{\alpha}$ when the integration surface passes through a degeneracy point of $U(\cdot)$, where the eigenvectors fail to be continuous.

3.2 Global expression for the $W_{3}$ invariant

Eq. (3.13) depends only on locally defined quantities. An alternative form of this equation can be obtained by introducing the globally defined eigenvalue bands of $U(\cdot)$. To this end, we note that the integrability condition $\epsilon_{\alpha\beta\gamma}\partial^{\beta}Y^{\gamma}=0$ holds (see Eq. (3.7)) such that the equation $\mathrm{i}\partial^{\alpha}\Phi(\mu)=Y^{\alpha}(\mu)$ has a solution $\Phi(\cdot)$ defined on the entire cube $[0,1]^{3}$. By construction, $\Phi(\mu)$ is a diagonal matrix with $D(\mu)=\exp[\mathrm{i}\Phi(\mu)]$. We can interpret the elements of this matrix as forming the bands $\mu\mapsto\Phi^{\nu}(\mu)$ of $U(\cdot)$. Every band is defined only up to a global shift by an integer multiple of $2\pi$. Especially at the boundaries $\mu_{\alpha}=0,1$ of the cube, where the eigenvalues of $U(\cdot)$ coincide, we have \Phi^{\nu}(\mu)\big{|}_{\mu_{\alpha}=1}=\Phi^{\nu}(\mu)\big{|}_{\mu_{\alpha}=0}+2\pi W^{\alpha,\nu}, with the winding number of the $\nu$-th eigenvalue from Eq. (3.15). In particular, the difference between the $\Phi^{\nu}(\mu)$ is constant on each boundary face.

Partial integration of Eq. (3.13) gives

[TABLE]

The first term contains the divergence of $F^{\alpha}$ that can be non-zero only [29] at degeneracy points of $U(\cdot)$. The second term gives the boundary contributions at $\mu_{\alpha}=0,1$. This equation, which contains the global bands $\Phi(\cdot)$, is the counterpart to Eq. (3.13) that uses only local quantities derived from $U(\cdot)$.

4 Computation of the $W_{3}$ invariant

To compute $W_{3}[U]$ we use a grid of $N\times N\times N$ discretization points $\mathbf{p}=(\delta i_{1},\delta i_{2},\delta i_{3})$, where $\delta=1/N$ for some $N>1$ and $1\leq i_{\alpha}\leq N$. To keep the presentation simple we assume an equidistant discretization grid, but generalization of the following expressions to non-equidistant grids is straightforward. In the following, $\boldsymbol{\delta}_{\alpha}$ denotes the lattice vector in direction $\alpha$, i.e., $\boldsymbol{\delta}_{1}=(\delta,0,0)$ etc., and we count coordinates $\mu_{\alpha}$ modulo $1$, i.e., $\mu+N\boldsymbol{\delta}_{\alpha}=\mu$. At every discretization point, we determine $S(\mathbf{p})$ and $D(\mathbf{p})$ through diagonalization of $U(\mathbf{p})$. We assume that, for given $\nu$, the eigenvalues $d^{\nu}(\mathbf{p})$ and eigenvectors $\mathbf{s}^{\nu}(\mathbf{p})$ on adjacent discretization points $\mathbf{p}$, $\mathbf{p}\pm\boldsymbol{\delta}_{\alpha}$ belong to the same band. To achieve this in practice, we can compare the distance $|d^{\nu}(\mathbf{p})-d^{\nu}(\mathbf{p}\pm\boldsymbol{\delta}_{\alpha})|$ of the eigenvalues or match the eigenvectors according to their overlap $|\langle\mathbf{s}^{\nu}(\mathbf{p}),\mathbf{s}^{\nu}(\mathbf{p}\pm\boldsymbol{\delta}_{\alpha})\rangle|$. Note that we do not need to consistently assign the indices $\nu$ on the entire discretization grid, but only locally on each of the cubes.

Faces

Every discretization point $\mathbf{p}$ is the base point of three rectangular faces $\mathsf{F}_{\mathbf{p},\alpha}$, with vertices as in Tab. 1. To every face $\mathsf{F}_{\mathbf{p},\alpha}$ we assign the real numbers

[TABLE]

which are derived from the eigenvectors of $U(\cdot)$ via the gauge variables

[TABLE]

defined along each edge of the face. Here, $\langle\cdot,\cdot\rangle$ denotes the complex Euclidean scalar product, and we can always avoid $\langle\mathbf{s}^{\nu}(\mathbf{p}_{i}),\mathbf{s}^{\nu}(\mathbf{p}_{j})\rangle=0$ by making the discretization finer. It is $\mathcal{U}^{\nu}(\mathbf{p}_{i},\mathbf{p}_{j})\in\mathbb{S}^{1}$, and $\hat{F}_{\mathbf{p},\alpha}^{\nu}\in\mathbb{R}$. For $\delta\to 0$, we have $\hat{F}_{\mathbf{p},\alpha}^{\nu}=F^{\nu}_{\alpha}(\mathbf{p})\,\delta^{2}+O(\delta^{3})$. This expression for $\hat{F}_{\mathbf{p},\alpha}^{\nu}$ is identical to the construction in Ref. [26].

Cubes

Every discretization point $\mathbf{p}$ is also the base point of a cube $\mathsf{C}_{\mathbf{p}}$, with opposite corner $\mathbf{p}+\boldsymbol{\delta}_{1}+\boldsymbol{\delta}_{2}+\boldsymbol{\delta}_{3}$ (see Fig. 1). To every cube we assign the number

[TABLE]

Every edge of the cube appears twice in this sum, contributing a term $\mathcal{U}^{\nu}(\mathbf{p}_{i},\mathbf{p}_{j})$ and $\mathcal{U}^{\nu}(\mathbf{p}_{j},\mathbf{p}_{i})=\mathcal{U}^{\nu}(\mathbf{p}_{i},\mathbf{p}_{j})^{-1}$. Combining the logarithms in the sum thus results in \exp\big{(}2\pi\mathrm{i}\hat{C}_{\mathbf{p}}^{\nu}\big{)}=1, such that $\hat{C}_{\mathbf{p}}^{\nu}\in\mathbb{Z}$. In the limit $\delta\to 0$, $\hat{C}_{\mathbf{p}}^{\nu}$ equals the volume integral of $\partial^{\alpha}F_{\alpha}$ over the cube $\mathsf{C}_{\mathbf{p}}$, and gives the Chern number of the $\nu$-th band of $U(\cdot)$ on an infinitesimal surface at $\mathbf{p}$. Therefore, $\hat{C}_{\mathbf{p}}^{\nu}\neq 0$ is possible only if $\mathsf{C}_{\mathbf{p}}$ contains a degeneracy point $\mu_{d}$ of $U(\cdot)$, where two eigenvalues $d^{\nu}(\mu_{d})=d^{\nu^{\prime}}(\mu_{d})$ for different $\nu\neq\nu^{\prime}$ coincide. In this situation, we have $\partial^{\alpha}F^{\nu}_{\alpha}(\mu_{d})=-\partial^{\alpha}F^{\nu^{\prime}}_{\alpha}(\mu_{d})$ in agreement with \operatorname{tr}\big{[}\partial^{\alpha}F_{\alpha}\big{]}=0, and $\hat{C}_{\mathbf{p}}^{\nu}=-\hat{C}_{\mathbf{p}}^{\nu^{\prime}}$ in agreement with $\sum_{\nu}\hat{C}_{\mathbf{p}}^{\nu}=0$.

$\hat{W}_{3}$ approximation

The approximation for the $W_{3}$ invariant computed by our algorithm is given by the sum

[TABLE]

over all faces and cubes, where the integers $m_{\mathbf{p},\alpha}^{\nu}$, $M_{\mathbf{p}}^{\nu}$ are obtained from the logarithms

[TABLE]

of the eigenvalues of $U(\cdot)$. On the one hand, we determine $m_{\mathbf{p},\alpha}^{\nu}$ such that $|\phi_{\mathbf{p}}^{\nu}-\phi_{\mathbf{p}-\boldsymbol{\delta}_{\alpha}}^{\nu}+2\pi m_{\mathbf{p},\alpha}^{\nu}|<\pi$. Equality can be excluded by making the discretization finer. On the other hand, if $\hat{C}_{\mathbf{p}}^{\nu}=-\hat{C}_{\mathbf{p}}^{\nu^{\prime}}\neq 0$ for two indices $\nu$, $\nu^{\prime}$, we pick one index, for example $\nu$, and determine $M_{\mathbf{p}}^{\nu}$ such that $|\phi_{\mathbf{p}}^{\nu}-\phi_{\mathbf{p}}^{\nu^{\prime}}+2\pi M_{\mathbf{p}}^{\nu}|<\pi$. Only the term $\hat{C}_{\mathbf{p}}^{\nu}M_{\mathbf{p}}^{\nu}$ is included in the sum, and we set $M_{\mathbf{p}}^{\nu^{\prime}}=0$. If we pick the index $\nu^{\prime}$ instead, we have $M_{\mathbf{p}}^{\nu^{\prime}}=-M_{\mathbf{p}}^{\nu}$ such that $\hat{C}_{\mathbf{p}}^{\nu^{\prime}}M_{\mathbf{p}}^{\nu^{\prime}}=\hat{C}_{\mathbf{p}}^{\nu}M_{\mathbf{p}}^{\nu}$. Therefore, it does not matter which of the two terms is included in the sum (4.4).

Note that the expression (4.4) is manifestly gauge-invariant, due to the construction of the $\hat{F}_{\mathbf{p},\alpha}^{\nu}$. Usually, only few faces and cubes with non-zero $m_{\mathbf{p},\alpha}^{\nu}$, $M_{\mathbf{p}}^{\nu}$ contribute in the sum.

Winding and Chern numbers

In addition to $\hat{W}_{3}$, the algorithm provides us with the approximation

[TABLE]

for the $W_{1}$ invariants of $U(\cdot)$. By construction, $\hat{W}^{\alpha}_{1}$ is an integer, and $\hat{W}_{1}^{\alpha}\to W^{\alpha}_{1}[U]$ for $\delta\to 0$.

Furthermore, we have the approximation

[TABLE]

for the Chern numbers $C^{\nu}_{\alpha}$ as in Eq. (3.17). Repeating the previous argument (for $\hat{C}_{\mathbf{p}}^{\nu}\in\mathbb{Z}$) we see that $\hat{C}^{\nu}_{\alpha}$ is an integer, and that $\hat{C}^{\nu}_{\alpha}\to C^{\nu}_{\alpha}$ for $\delta\to 0$. This expression for $\hat{C}^{\nu}_{\alpha}$ is identical to the construction in Ref. [26].

4.1 Justification of the algorithm

To justify the construction of our algorithm we first note that the value of $\hat{W}_{3}$ does not depend on the position of the branch cut of the complex logarithm in Eq. (4.5). Indeed, if we change one $\phi_{\mathbf{p}}^{\nu}\mapsto\phi_{\mathbf{p}}^{\nu}+2\pi M$ by some $M\in\mathbb{Z}$, the integers in Eq. (4.4) change as $m_{\mathbf{p},\alpha}^{\nu}\mapsto m_{\mathbf{p},\alpha}^{\nu}-M$ and $m_{\mathbf{p}+\boldsymbol{\delta}_{\alpha},\alpha}^{\nu}\mapsto m_{\mathbf{p}+\boldsymbol{\delta}_{\alpha},\alpha}^{\nu}+M$. This gives an additional contribution $M\hat{C}_{\mathbf{p}}^{\nu}$ in the sum (4.4). Now if $\hat{C}_{\mathbf{p}}^{\nu}\neq 0$ we have $M_{\mathbf{p}}^{\nu}\mapsto M_{\mathbf{p}}^{\nu}-M$ (or $M_{\mathbf{p}}^{\nu^{\prime}}\mapsto M_{\mathbf{p}}^{\nu^{\prime}}+M$ for the associated second index $\nu^{\prime}$). Therefore, the sum (4.4) remains unchanged. As one implication, we note that it does not matter whether we determine the integers $m_{\mathbf{p},\alpha}^{\nu}$ using the grid points $\mathbf{p}$, $\mathbf{p}-\boldsymbol{\delta}_{\alpha}$, as we have specified before, or using the points $\mathbf{p}$, $\mathbf{p}+\boldsymbol{\delta}_{\alpha}$.

Furthermore, this freedom in Eq. (4.5) allows us to consider two opposite choices for the $\phi_{\mathbf{p}}^{\nu}$ (even if the evaluation of Eq. (4.4) does not require us to make any such choice). On the one hand, we can choose them in such a way that, whenever a $\hat{C}_{\mathbf{p}}^{\nu}\neq 0$, the two associated $\phi_{\mathbf{p}}^{\nu}$, $\phi_{\mathbf{p}}^{\nu^{\prime}}$ belong to the same branch of the complex logarithm. Then, all $M_{\mathbf{p}}^{\nu}=0$.

On the other hand, we can choose the $\phi_{\mathbf{p}}^{\nu}$ in such a way that they differ by less than $\pi$ on adjacent discretization points $\mathbf{p}$, $\mathbf{p}+\boldsymbol{\delta}_{\alpha}$ in the interior of the discretization grid. For a sufficiently fine111We can understand Eq. (4.8) as an admissibility condition that determines the minimum possible number of discretization points. discretization with

[TABLE]

this choice is equivalent to setting $\phi_{\mathbf{p}}^{\nu}=\Phi^{\nu}(\mathbf{p})$ with the ‘global’ bands of $U(\cdot)$ as in Sec. 3.2. Then, all $m_{\mathbf{p},\alpha}^{\nu}=0$ except for the boundary values at $i_{\alpha}=1$ that account for the total change of $\Phi^{\nu}(\mu)$ from $\mu_{\alpha}=0$ to $\mu_{\alpha}=1$. Explicitly, we can set $m_{\mathbf{p}}^{\nu}=W^{\alpha,\nu}$ at $i_{\alpha}=1$, with the winding number of band $\nu$ as in Eq. (3.15). To prove correctness of our algorithm we now adopt this choice, and then compare Eq. (4.4) with Eq. (3.18). We assume that the discretization is admissible according to Eq. (4.8).

First proposition

The first term in Eq. (4.4) is an integer by construction. With the present choice, the second term can be rewritten as $\sum_{\nu=1}^{n}\sum_{\alpha=1}^{3}\hat{C}^{\nu}_{\alpha}W^{\alpha,\nu}$ with the approximate Chern numbers $\hat{C}^{\nu}_{\alpha}$ from Eq. (4.7), which also gives an integer. This shows our first proposition: $\hat{W}_{3}$ is an integer for arbitrary (admissible) discretizations.

Second proposition

Now consider the limit $\delta\to 0$. We see immediately that the second term in Eq. (4.4) converges to the second term in Eq. (3.18). For the first term in Eq. (4.4), we note that $\hat{C}_{\mathbf{p}}^{\nu}\neq 0$ occurs only if the cube $\mathsf{C}_{\mathbf{p}}$ contains a degeneracy point $\mu_{d}$ of $U(\cdot)$. For two degenerate eigenvalues $d^{\nu}(\mu_{d})=d^{\nu^{\prime}}(\mu_{d})$, with $\partial^{\alpha}F^{\nu}_{\alpha}(\mu_{d})=-\partial^{\alpha}F^{\nu^{\prime}}_{\alpha}(\mu_{d})$ and \big{(}\Phi^{\nu^{\prime}}(\mu_{d})-\Phi^{\nu}(\mu_{d})\big{)}/2\pi\in\mathbb{Z}, the first term in Eq. (3.18) gives the contribution (1/2\pi)\big{(}\partial^{\alpha}F^{\nu}_{\alpha}(\mu_{d})\big{)}\big{(}\Phi^{\nu^{\prime}}(\mu_{d})-\Phi^{\nu}(\mu_{d})\big{)}. On the other hand, the first term in Eq. (4.4) gives the contribution $\hat{C}_{\mathbf{p}}^{\nu}M_{\mathbf{p}}^{\nu}$, and we have $2\pi M_{\mathbf{p}}^{\nu}=\phi_{\mathbf{p}}^{\nu^{\prime}}-\phi_{\mathbf{p}}^{\nu}=\Phi^{\nu^{\prime}}(\mu_{d})-\Phi^{\nu}(\mu_{d})$ for $\mathbf{p}\to\mu_{d}$. Both terms agree in the limit $\delta\to 0$. This shows our second proposition: $\hat{W}_{3}$ converges to $W_{3}[U]$ in the limit $\delta\to 0$.

In combination, the second proposition shows the correctness of the algorithm, while the first proposition justifies our constructions as it promises rapid convergence of the algorithm already on coarse discretization grids.

4.2 Example application: $\mathrm{SU}(2)$ maps

To demonstrate the application of our algorithm, we consider maps $\mathbb{T}^{3}\to\mathbb{U}(2)$ of spin-$\tfrac{1}{2}$ rotations

[TABLE]

Here, $\boldsymbol{\sigma}$ is the vector of Pauli matrices, and $\mathbf{a}(\cdot)$ with $\mathbf{a}(\mu)\in\mathbb{R}^{3}$ defines a vector field on $[0,1]^{3}$. To obtain a map $U(\cdot)$ defined on the torus $\mathbb{T}^{3}$, suitable boundary conditions must be fulfilled by $\mathbf{a}(\cdot)$. In this context, recall that $U(\mu)=(-1)^{n}\,\mathbbm{1}$ is a multiple of the identity if $|\mathbf{a}(\mu)|/2\pi\in\mathbb{N}_{0}$, independent of the direction of $\mathbf{a}(\mu)$. We consider two examples.

First example

For the first example, we choose a map $\mathbf{f}(\cdot)$ from the square $[0,1]^{2}$ to the unit sphere $\mathbb{S}^{2}\subset\mathbb{R}^{3}$ that maps the boundary of the square to the south pole of the sphere, and the interior of the square bijectively onto the punctured sphere without the south pole. Then, we set $\mathbf{a}(\mu)=4\pi w\,\mu_{3}\,\mathbf{f}(\mu_{1},\mu_{2})$. Details of the map $\mathbf{f}(\cdot)$ do not matter for the value of the topological invariant $W_{3}[U]$. For this map, the Berry curvature of either band of $U(\cdot)$ is a function of $\mu_{1},\mu_{2}$, as illustrated in Fig. 2, while the eigenvalues are a function of $\mu_{3}$. In Eq. (3.13) the three-dimensional integral of $F_{\alpha}Y^{\alpha}$ thus separates. Integration of the Berry curvature $F_{\alpha}$ gives the Chern number $\pm 1$ for either band, and we have $W_{3}[U]=2w$.

Second example

For the second example, we choose a bijective map $\mathbf{g}(\cdot)$ from the cube $[0,1]^{3}$ to the unit ball $|\mathbf{r}|\leq 1$ that maps the surface (center) of the cube to the surface (center) of the unit ball. Then, we set $\mathbf{a}(\mu)=2\pi w\,\mathbf{g}(\mu)$. The corresponding map $U(\cdot)$ has a degeneracy point at the center of the cube. Additional degeneracies occur on concentric shells where $w|\mathbf{g}(\mu)|\in\mathbb{N}$. Every additional degeneracy changes the $W_{3}$ invariant by $\pm 1$, and we have $W_{3}[U]=w$.

We now apply our algorithm to both maps. If the values $\phi_{\mathbf{p}}^{\nu}$ are computed with the standard complex logarithm, non-zero integers $m_{\mathbf{p},\alpha}^{\nu}$ contribute in the sum (4.4) at the positions shown in Fig. 2. For the first map, whose eigenvalues of $U(\cdot)$ change as a function of $\mu_{3}$, non-zero $m_{\mathbf{p},\alpha}^{\nu}$ occur in a sheet of the cube where the eigenvalues pass through $e^{\mathrm{i}\pi}$ as they move around the circle $\mathbb{S}^{1}$. For the second map, the eigenvalues depend on the radial distance to the center and a more complicated pattern appears. The cube around the central degeneracy point of $U(\cdot)$, where $\hat{C}_{\mathbf{p}}^{\nu}\neq 0$, is visible in this pattern. As explained before in connection with Eq. (4.4), any other computation of $\phi_{\mathbf{p}}^{\nu}$ would result in a different arrangement of the non-zero $m_{\mathbf{p},\alpha}^{\nu}$ but produce identical results for $\hat{W}_{3}[U]$.

The fast convergence of our algorithm is illustrated in Fig. 3, where it is compared to direct evaluation of Eq. (2.2). The latter, with approximation of the derivatives $U^{-1}\partial^{\alpha}U$ of $U(\cdot)$ by central differences222As in $f^{\prime}(x)\approx(f(x+\delta)-f(x-\delta))/(2\delta)$., converges much more slowly. A simple convergence criterion for our algorithm is to check that the angle between eigenvalues of $U(\cdot)$ on adjacent discretization points $\mathbf{p}$, $\mathbf{p}\pm\boldsymbol{\delta}_{\alpha}$ stays below a certain threshold, say, below $\pi/2=90^{\circ}$. This criterion is modelled after the admissibility condition (4.8). It guarantees that eigenvalues do not move around the entire circle $\mathbb{S}^{1}$ between adjacent discretization points. As the data in Fig. 3 show, correct results are already obtained with even coarser grids. Fig. 3 also shows the required number $N$ to achieve convergence on an $N\times N\times N$ grid. For the example maps, larger $N$ are required with increasing $w$ as the eigenvalues of $U(\cdot)$ change more quickly as a function of $\mu$. In accordance with the above criterion, the required $N$ grows approximately linearly with $w$.

5 The $W_{3}$ invariant for Floquet-Bloch systems

For the application of the $W_{3}$ invariant to Floquet-Bloch systems, as introduced in Ref. [12], we start with the time-dependent Hamiltonian $H(t)$ of a two-dimensional lattice system. With translational symmetry in the $x$, $y$ lattice direction, the Hamiltonian can be parametrized by the $k_{x}$, $k_{y}$ momentum as $H(k_{x},k_{y},t)$. Solution of the Schrödinger equation, either by direct integration [30] or for periodic time-dependence within the Floquet formalism [31], gives the propagator $U(k_{x},k_{y},t)$, with

[TABLE]

The unitary map $U(k_{x},k_{y},t)$ is periodic in the $k_{x}$, $k_{y}$ direction, but not in the $t$ direction. For such $2+1$ dimensional maps, a construction suggested in Ref. [12] provides a link between the $W_{3}$ invariant and the number of edge states in Floquet topological insulators. We start with the definition of the respective invariant $W_{3}[U_{\xi}]$, and then explain its computation with the algorithm from Sec. 4, which requires only a small modification of the previous scheme. The following constructions apply to arbitrary time-dependent Hamiltonians, although applications will often concern Floquet-Bloch systems with a periodic Hamiltonian $H(t+T)=H(t)$. Note that in general the Brillouin zone is not a square, but we can always assume a mapping of the crystal quasimomentum $(k_{x},k_{y})$ in the unit cell of the reciprocal lattice to the coordinates $(\mu_{1},\mu_{2})\in[0,1]^{2}$, and also a mapping from time $t\in[0,T]$ to the coordinate $\mu_{3}\in[0,1]$. Therefore, we will use either set of coordinates wherever appropriate. As before, the topological invariants do not depend on details of the mapping.

5.1 The invariant $W_{3}[U_{\xi}]$

The invariant $W_{3}[U_{\xi}]$ is defined for unitary maps $U(\cdot):\mathsf{[}0,1]^{3}\to\mathbb{U}(n)$ that are periodic in the $\alpha=1,2$ directions, as in the first two equalities in Eq. (2.1), and fulfill $U(\mu_{1},\mu_{2},0)=\mathbbm{1}$. For such a map, choose a $\xi\in\mathbb{S}^{1}$ that differs from all eigenvalues $d^{\nu}(\mu_{1},\mu_{2},1)$ at $\mu_{3}=1$, i.e., that lies in a gap between the bands $(\mu_{1},\mu_{2})\mapsto d^{\nu}(\mu_{1},\mu_{2},1)$ of $U(\cdot,\mu_{3}=1)$. For time-periodic systems, these bands are often called the Floquet bands. Note that the existence of a gap implies that the winding numbers $W^{1,\nu}$, $W^{2,\nu}$ of every such band are zero, which agrees with the assumption $U(\mu_{1},\mu_{2},0)=\mathbbm{1}$.

In this situation, we can consider the map

[TABLE]

where $\log_{\xi}z$ is the complex logarithm with branch cut along the line from zero through $\xi$. As $s=0$ increases to $s=1$, the map $s\mapsto V_{\xi}(s,e^{\mathrm{i}\phi})$ rotates eigenvalues $e^{\mathrm{i}\phi}$ on the two different circle segments from one to $\xi=e^{\mathrm{i}\epsilon}$ clockwise (as $s\mapsto\exp(\mathrm{i}(1-s)\phi)$ for $0\leq\phi<\epsilon\leq 2\pi$) or counterclockwise (as $s\mapsto\exp(\mathrm{i}(1-s)\phi+2\pi\mathrm{i}s)$ for $0\leq\epsilon<\phi\leq 2\pi$) to one (see Fig. 4). The invariant $W_{3}[U_{\xi}]$ is now computed for the map

[TABLE]

which is periodic in all directions $\alpha=1,2,3$.

In Eq. (3.13) the contribution from the $V_{\xi}(\cdot)$ part of Eq. (5.3) can be evaluated explicitly, and we have

[TABLE]

where the value of $W_{3}[U]$ is obtained by evaluation of Eq. (3.13) for the map $U(\cdot)$ (which by itself does not give a topological invariant because $U(\cdot)$ is not periodic in the $\alpha=3$ direction). Note that the $\xi$-dependent correction term in this expression depends only on the Berry curvature and eigenvalues of $U(\cdot)$ at $\mu_{3}=1$.

When the value $\xi$ changes from a gap at $\xi=\xi_{a}$ to a gap at $\xi=\xi_{b}$ it will pass through some bands $\nu_{1},\nu_{2},\dots,\nu_{k}$ of $U(\cdot,\mu_{3}=1)$ (see Fig. 4). For eigenvalues in these bands, the logarithms $\log_{\xi}d^{\nu}(\mu_{1},\mu_{2},1)$ in Eq. (5.4) change by $\pm 2\pi\mathrm{i}$ (for example, $\log_{\xi_{b}}d^{\nu}=\log_{\xi_{a}}d^{\nu}+2\pi\mathrm{i}$ for $\nu=1,2,3$ in Fig. 4). From comparison with Eq. (3.17) we see that the entire expression changes by the Chern numbers $C^{\nu}_{3}$ of these bands. We thus get the relation

[TABLE]

between the $W_{3}[U_{\xi}]$ invariant in two gaps $\xi_{a}$, $\xi_{b}$, and the Chern numbers of the bands lying on the counterclockwise circle segment from $\xi_{a}$ to $\xi_{b}$. If $\xi$ moves once around the entire circle, we recover333This observation also resolves the potential ambiguity that occurs in Eq. (5.2) for $\xi\to 1$. the identity $W_{3}[U_{\xi}]=W_{3}[U_{\xi}]$ since the Chern numbers add up to $\sum_{\nu=1}^{n}C^{\nu}_{\alpha}=0$.

5.2 Computation of $W_{3}[U_{\xi}]$

We can immediately apply the algorithmic scheme of Sec. 4 to the computation of $W_{3}[U_{\xi}]$ if we absorb the $\xi$-dependent term in Eq. (5.4) into a correction to Eq. (4.4). This correction results only from the difference of $\phi_{\mathbf{p}}^{\nu}$ and $-\mathrm{i}\log_{\xi}d^{\nu}(\mathbf{p})$ at $i_{3}=N$.

First, we apply the previous algorithm to the map $U(\cdot)$ and compute $\hat{W}_{3}$ from the sum (4.4) (which by itself does not give an integer because $U(\cdot)$ is not periodic). Because of $U_{\xi}(\cdot)=\mathbbm{1}$ for $\mu_{3}=0,1$, we can set $m^{\nu}_{\mathbf{p},\alpha}=M^{\nu}_{\mathbf{p}}=0$ for $i_{3}=1$, such that the sum starts effectively with $i_{3}=2$, and $M^{\nu}_{\mathbf{p}}=0$ for $i_{3}=N$. Second, to account for the $\xi$-dependent correction, we determine integers $K_{i_{1},i_{2}}^{\nu}$ such that $|-\mathrm{i}\log_{\xi}d^{\nu}(\mathbf{p})-\phi_{\mathbf{p}}^{\nu}+2\pi K_{i_{1},i_{2}}^{\nu}|<\pi$ at $i_{3}=N$. Then, we have

[TABLE]

As before, $\hat{W}_{3}[U_{\xi}]$ is an integer and converges to $W_{3}[U_{\xi}]$ for $\delta\to 0$. In this expression, only the second term depends on $\xi$, and reuses quantities already computed in the algorithm.

Note that we can repeat the argument leading to Eq. (5.5): If $\xi$ moves from gap $\xi=\xi_{a}$ to another gap $\xi=\xi_{b}$ the integers in Eq. (5.6) change by $K_{i_{1},i_{2}}^{\nu}\mapsto K_{i_{1},i_{2}}^{\nu}-1$ for bands $\nu=1,\dots,k$ lying between the two gaps. According to Eq. (3.17), the value of $\hat{W}_{3}[U_{\xi}]$ thus changes by the approximate Chern numbers $-(\hat{C}^{1}_{3}+\cdots+\hat{C}^{k}_{3})$ of the respective bands, which agrees with Eq. (5.5). The Chern numbers $\hat{C}^{\nu}_{3}$ are also provided by the algorithm.

5.3 Tracking the evolution of $W_{3}[U_{\xi}]$ for gapped systems

The invariant $W_{3}[U_{\xi}]$ can be considered as a function of time $t$, and in physics applications we might want to track the change of $W_{3}[U_{\xi}]$ with increasing $t$. Note that $W_{3}[U_{\xi}]$ depends on the entire history of $U(k_{x},k_{y},t^{\prime})$ for $0\leq t^{\prime}\leq t$, not only on the map $U(k_{x},k_{y},t)$ at the final time $t$.

To describe how the change of $W_{3}[U_{\xi}]$ can be tracked we assume that, apart from isolated degeneracy points, the system is fully gapped for all $t>0$. Then, the $n$ bands $\nu=1,\dots,n$ of $U(\cdot,t)$ are separated by $n$ gaps at positions $\xi_{1},\dots,\xi_{n}$, as illustrated in Fig. 5. Because the relation $d^{\nu}\equiv\exp(-\mathrm{i}\epsilon^{\nu}t)$ between eigenvalues $d^{\nu}$ of $U(\cdot,t)$ and (Floquet quasi–) energies $\epsilon^{\nu}$ involves a minus sign, we number bands and gaps clockwise on $\mathbb{S}^{1}$. If we count the index $\nu$ modulo $n$, we have $\xi_{n}=\xi_{0}$.

To each gap, we assign the value $n^{\nu}=W_{3}[U_{\xi_{\nu}}]$. It does not depend on the precise location of $\xi_{\nu}$ within the gap, but only on the gap index $\nu$. According to the bulk-edge correspondence for Floquet-Bloch systems, $n^{\nu}$ gives the number of edge states in the respective gap [12].

According to Eq. (5.5), the different $n^{\nu}$ are related to the Chern numbers of the (Floquet) bands via

[TABLE]

For small times $t$, before the bands can wrap around the circle $\mathbb{S}^{1}$, we have $W_{3}[U_{\xi_{n}}]=0$ and thus

[TABLE]

Especially, we have $C^{1}_{3}+\dots+C^{n}_{3}=0$ in accordance with $n^{n}=n^{0}=0$.

As the systems evolves in time, gaps close and reopen at isolated [28] degeneracy points of $U(\cdot)$. Computing the precise location of these points is impractical. Instead, we want to propagate the system over an entire step $t\mapsto t+\delta t$ (or $\mu_{3}\mapsto\mu_{3}+\delta$ in $\mu$ coordinates), and determine the values $\tilde{n}^{\nu}$ at $t+\delta t$ from the values $n^{\nu}$ at $t$.

Within our algorithm, performing a time step corresponds to extending the discretization grid from $i_{3}=N$ to $i_{3}=N+1$. How the values $n^{\nu}$ change can now be deduced from Eq. (5.6). In the special situation of a fully gapped system considered here, where the same gap $\nu$ closes and reopens at a degeneracy point, all the integers $K_{i_{1},i_{2}}^{\nu}$ stay the same in this equation444To be precise, they stay the same if we set $\phi^{\nu}_{\mathbf{p}}=-\mathrm{i}\log_{\xi}d^{\nu}(\mathbf{p})$. As noted earlier in connection with Eq. (4.4), the choice of the branch cut of the complex logarithm does not affect the algorithm.. The value $n^{\nu}$ changes only if, by extending the grid, we capture a new degeneracy point of $U(\cdot)$ that gives a non-zero $M^{\nu}_{\mathbf{p}}\hat{C}_{\mathbf{p}}^{\nu}=M^{\nu+1}_{\mathbf{p}}\hat{C}_{\mathbf{p}}^{\nu+1}$ contribution at $i_{3}=N$ for the bands $\nu$, $\nu+1$ next to the gap [19]. Recalling the definition of $M^{\nu}_{\mathbf{p}}$ in Eq. (4.4) and of $\log_{\xi}$ in Eq. (5.6), we find that $M^{\nu}_{\mathbf{p}}=-M^{\nu+1}_{\mathbf{p}}=1$ if $\xi=\xi_{\nu}$ lies in the gap between bands $\nu$, $\nu+1$, and $M^{\nu}_{\mathbf{p}}=M^{\nu+1}_{\mathbf{p}}=0$ otherwise.

Therefore, if the algorithm detects a non-zero $\hat{C}_{\mathbf{p}}^{\nu}=-\hat{C}_{\mathbf{p}}^{\nu+1}\neq 0$ at the present $i_{3}=N$, a degeneracy point in gap $\nu$ occurs between $t$ and $t+\delta t$. The values $n^{\nu}$ change as

[TABLE]

At the same time, the Chern numbers of the two bands next to the gap change as $\tilde{C}^{\nu}_{3}=C^{\nu}_{3}+\hat{C}_{\mathbf{p}}^{\nu}$, $\tilde{C}^{\nu+1}_{3}=C^{\nu+1}_{3}-\hat{C}_{\mathbf{p}}^{\nu}$. This replacement leaves the relation (5.7) intact. If more than two non-zero $\hat{C}_{\mathbf{p}}^{\nu}$ occur in the same slice, they have to be summed.

The above scheme to track the change of the values $n^{\nu}$ directly is a certain simplification of the more general algorithm from Sec. 5.2. We should note, however, that the tracking scheme applies primarily to fully gapped systems as in Fig. 5, while the general algorithm is preferable in situations as in Fig. 4 where bands overlap during time evolution. Note that in the practical computation, the cubes with non-zero $\hat{C}_{\mathbf{p}}^{\nu}$ that enter in Eq. (5.9) are simply obtained with Eq. (4.3) evaluated on the discretization grid, as explained previously for our algorithm in Sec. 4. In particular, we do not need to determine the precise positions of the degeneracy points of $U(\cdot)$.

5.4 Specialization to static systems

Let us briefly note how the invariant $W_{3}[U_{\xi}]$ generalizes the bulk-edge correspondence for static systems, which involves only the Chern numbers of the respective bands, and how this fact manifests itself in our algorithm.

For static systems with a constant Hamiltonian $H(k_{x},k_{y},t)\equiv H(k_{x},k_{y})$, the matrix $S(\cdot)$ in the decomposition (3.1) does not depend on $t\equiv\mu_{3}$, while the eigenvalues $d^{\nu}(\mu)=\exp(-\mathrm{i}\mu_{3}E^{\nu}(\mu_{1},\mu_{2}))$ change at a constant rate given by the eigenvalues $E^{\nu}(\mu_{1},\mu_{2})\equiv E^{\nu}(k_{x},k_{y})$ of $H(k_{x},k_{y})$, very much as in the first example map in Sec. 4.2.

For small $t\equiv\mu_{3}$, we have a situation as in Fig. 5 (left panel), and according to the previous considerations the values $n^{\nu}=W_{3}[U_{\xi}]$ in the different gaps differ by the (constant) Chern numbers of the bands of $H(k_{x},k_{y})$. In this way, we recover with Eq. (5.8) the bulk-edge correspondence for static systems, which relates the number of edge states $n^{\nu}$ and the Chern numbers $C^{\nu}_{3}$ of the bands of $H(k_{x},k_{y})$. Our algorithm provides the numbers $n^{\nu}$ immediately, in the first step at $i_{3}=1$: Once as the result of the direct computation of $W_{3}[U_{\xi}]$, and again via summation of the approximate Chern numbers $\hat{C}^{\nu}_{3}$ from Eq. (3.17).

6 Application: Anomalous edge states in irradiated graphene

After introduction of our algorithm we now want to apply it to a representative Floquet-Bloch system. To this end, we consider the Hamilton operator [32, 33, 34, 35] (omitting irrelevant constants)

[TABLE]

of a $2$-band tight-binding model on a graphene lattice, where irradiation with circularly polarized light is included through a time-periodic Peierls [36] phase $A_{ij}(t)=\mathbf{A}(t)\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})$. Here, $\mathbf{R}_{i}-\mathbf{R}_{j}$ is the lattice vector from site $j$ to site $i$, and $\mathbf{A}(t)=A_{0}\left(\sin(\omega t),\cos(\omega t)\right)$ gives the light wave with amplitude $A_{0}$ and frequency $\omega$. Without driving ($A_{0}=0$) the Hamiltonian (6.1) reduces to the standard graphene tight binding model, with linear electron dispersion at the Dirac points. There are no chiral edge states and the system is topologically trivial.

Irradiation with light ($A_{0}\neq 0$) induces a band gap at the Dirac points, and chiral edge states appear in this gap. Through variation of the light amplitude and frequency various topological phase transitions can be observed [35]. In remarkable contrast to the static case, where additional coupling terms such as spin orbit coupling are needed in the Hamiltonian to induce the transition from a topologically trivial to a nontrivial state, simple time-dependent variation of the Peierls phase suffices to establish a Floquet topological insulator.

Here, we consider a parameter regime ($A_{0}=0.7$ and $\omega=3.5$) where anomalous edge states appear in the Floquet band structure (see Fig. 6). The two Floquet bands, between the gaps at $\xi_{1}=1$ and $\xi_{2}=-1$, have Chern numbers $C_{3}^{1}=-3$ and $C_{3}^{2}=3$. For static systems, we would expect from the bulk-edge correspondence (5.8) that three chiral edge states exist in the band gap between the two bands. However, Floquet quasienergies are defined only up to multiples of $2\pi/T$ and no “lowest” or “highest” band exists, such that it is no longer possible to predict the number of edge states from the Chern numbers alone. Indeed, Fig. 6 shows that the correct net-chirality of edge states in the two band gaps is $n^{1}=-1$ and $n^{2}=2$. The relation (5.7) holds instead of the simpler Eq. (5.8). Therefore, computation of the winding number $W_{3}[U_{\xi}]$ is essential for the topological classification of such a Floquet-Bloch system.

We now compute $n^{1}$ and $n^{2}$ directly with the algorithm introduced in Sec. 5. In the following, we map the crystal momentum $k_{x}$, $k_{y}$ to a square with coordinates $\mu_{1},\mu_{2}$, which is equivalent to a rhombic discretization of the hexagonal Brillouin zone of graphene. Fig. 7 shows the cuboids with non-zero $\hat{C}^{\nu}_{\mathbf{p}}$ that occur in the algorithm for a fine (coarse) grid with $16\times 16\times 16$ ($6\times 6\times 6$) discretization points. As discussed previously, such cuboids contain a degeneracy point where two (Floquet) bands touch. A close-up view of such a degeneracy point is shown in Fig. 8, where we zoom into the cuboid marked in Fig. 7 (left panel) with a red circle. At the degeneracy point, local expansion of $U(\cdot)$ would provide the necessary information to predict how the values $W_{3}[U_{\xi}]$ and $n^{\nu}$ change as the gap closes, but isolation of the degeneracy point is impractical and requires high computational effort with many evaluations of $U(\cdot)$. Instead, comparison of the computation on the fine and coarse grid shows that the present algorithm produces the correct values $n^{1}=-1$ and $n^{2}=2$ already on the coarse grid. In the present example, this grid also gives the minimal number of discretization points required to obtain correct results. Coarser grids, for which wrong results would be produced by the algorithm, violate the admissibility conditions formulated in Sec. 4. The simple check from Sec. 4.2 on the maximum angle between eigenvalues on adjacent discretization points (as in Fig. 3) guarantees correctness of the algorithm.

7 Conclusions

In this paper we introduce an algorithm for the efficient computation of the $W_{3}$ invariant, and demonstrate its application to general unitary maps and to the propagators of time-dependent (Floquet-) Bloch systems. In both cases, the $W_{3}$ invariant is computed with moderate effort, requiring evaluation of the unitary map only at a small number of discretization points.

The construction of the algorithm relies on an expression for the $W_{3}$ invariant in terms of the Berry curvature of the eigenvectors and the angular velocity of the eigenvalues of the unitary map. Rapid convergence with the number of discretization points is guaranteed by the fact that the algorithm generates, indepently of the discretization, an integer number. In this respect, our algorithm is related to the algorithm of Fukui, Hatsugai, Suzuki [26], where a similar construction is used for the computation of Chern numbers. We note specifically that our algorithm does not attempt to simply evaluate Eq. (2.2) with some differentiation and integration rule. Instead, it directly extracts topological information about $U(\cdot)$ that is contained in the definition of the $W_{3}$ invariant. In addition to the computation of Chern numbers the algorithm for the $W_{3}$ invariant must also account for the winding of eigenvalues. In this respect, our algorithm extends the algorithm from Ref. [26].

Concerning the physical applications of the algorithm, we note that the $W_{3}$ invariant can be computed for arbitrary time-dependence, not only for the periodic time-dependence of Floquet-Bloch systems. While observation of topological properties for non-periodic time-dependence may be difficult in solid state systems, one should be able to use photonic crystals [37, 38] to implement such scenarios. The present algorithm will then allow us to determine topological properties also in complicated situations where alternative methods of computation fail.

From a mathematical point of view, one might want to note that the $W_{3}$ invariant in Eq. (2.2) is entirely formulated in terms of the unitary Lie group and Lie algebra, while the present algorithm explicitly refers to the matrix realization of their elements. Out of mathematical curiosity it is natural to ask for a “Lie-type” algorithm that does not require spectral decomposition of the unitary map. Such an algorithm can be easily constructed for the $W_{1}$ invariant, but we do not know of a similar algorithm for the $W_{3}$ (or higher dimensional) invariant that preserves the favorable properties of the present algorithm, e.g., by producing integer number results only. Clearly, for the $W_{3}$ invariant the situation complicates because the different derivatives in Eq. (2.2) do not need to commute. In any case, for the physical applications considered here the present formulation in terms of eigenvectors and eigenvalues is both conceptually and computationally adequate. The interpretation that non-trivial topology is generated in time-dependent systems through band “crossings” at degeneracy points gives the correct picture for the physical applications. The implementation of the present algorithm is reasonably easy, and as application to the irradiated graphene system demonstrates the algorithm provides us with a computationally efficient solution to the problem of computing the $W_{3}$ invariant.

Acknowledgments

This work was supported in part by Deutsche Forschungsgemeinschaft through SFB 652 (project no. B5).